title: Query SDE to display records created for 'X' date/time

text: I am looking to query our SDE to display feature classes modified beginning from January 1, 2013 to today. I am using Microsoft SQL Server Management Studio, and as it stands searching through dbo tables but can't find any information that would display feature classes modified. The purpose is to build a business case by identifying a volume of edits/modifications to our data over time. Any help is truly appreciated. ![Below is an image of my work environment](http://i.stack.imgur.com/VUo3k.jpg)


title: How do I connect to geoserver on my system using java code?

text: I am writing the following code... But i am getting an error as... java.io.FileNotFoundException: http:// localhost:8080/geoserver/ at sun.net.www.protocol.http.HttpURLConnection.getInputStream(HttpURLConnection.java:1613) at Try.main(Try.java:25) Am i doing it right? I am new at this geoserver as well as client server coding...               import java.io.*;     import java.net.*;          public class Try{              public static void main(String args[]) throws Exception{                  URL hp= new URL ("http://localhost:8080/geoserver");                       try {                      URL myURL = new URL("http://localhost:8080/geoserver/");                      URLConnection myURLConnection = myURL.openConnection();                      myURLConnection.connect();                      BufferedReader in = new BufferedReader(new InputStreamReader (myURLConnection.getInputStream()));                      String inputLine;                      while((inputLine = in.readLine()) != null )                 {                     System.out.println(inputLine);                      }                      in.close();                  }             catch (Exception e){                 System.out.println("Into CATCH");                 e.printStackTrace();             }                    }     }


title: Why a severe performance difference using deleteSearchedRows

text: The issue I am having involves two versioned tables. A function deletes 48 rows from Table A and and 48 rows from Table B using deleteSearchedRows. However, the delete from Table A takes roughly 4 seconds while the delete from table B takes roughly 67 seconds. The query filter for both deletes is a single number field that has been indexed in the table, objectID indexed as well. Table A has roughly 45,000 records and Table B 60,000 records. Both tables also have cascade deletes set to No, so that the only data being deleted would be that table's. Table structures are very similar. Table A has OID field + 3 Long integer fields, 2 short integer fields and a date field. Table B has OID, 3 long integer, 1 Text, and a date field. Running the delete query from SQL Plus is very quick, so it's probably not a drive efficiency issue. I broke apart the select and delete sections to see if the search was slow. I used the Search(queryFilter, false) and it returns instantly, but looping and deleting every IRow.Delete takes slightly over a second per row for Table B. I am at a loss on what the difference could be to cause such a severe performance difference in two tables.   Any ideas/debug methods are appreciated. Using ArcGIS 10.0 with an Oracle database


title: Run function if item in QComboBox is selected

text: I got another issue writing a plugin to convert three shape files into an .obs input file. I created a UI (Qt) with a QComboBox to select a layer (that works fine). As soon as one of the layers in the QComboBox is selected, a second ComboBox should be filled with the fields of the layer to select from. To do that I'm trying to call a function, that is to be called, as soon, as a layer is selected in the first ComboBox. Unfortuntelly that does not work. Here is what I tried: First I wrote a function to fill the second ComboBox with the field names:                 def build_Layer_Index_Changed():         QMessageBox.information(None, "DEBUG", "Build Layer routine called. :)")         canvas = self.iface.mapCanvas()         dlg2 = ShowDialog2()          build = canvas.layer(int(dlg2.Sel_Build_Layer.currentIndex()))         if int(dlg2.Sel_Build_Layer.currentIndex()) >= 0:           if not build.geomTypeDict == 2:   # Types: 0:"point",1:"line",2:"polygon", Buildings layer has to be of type 2             QMessageBox.critical(self.iface.mainWindow(), ":-/ Selected buillding layer is not of type vector(polygon)! :-/" )            else:             build_fields = self.getFields(build)             dlg2.Sel_buildhight_field.clear()             dlg2.Sel_buildhight_field.addItems(build_fields)      "dlg2.Sel_buildhight_field" is the second QComboBox. I tried to connect it to the first ComboBox's (dlg2.Sel_Build_Layer) Signal by               QtCore.QObject.connect(dlg2.Sel_Build_Layer, QtCore.SIGNAL("currentIndexChanged"), self.build_Layer_Index_Changed())      but I get an error: "TypeError: build_Layer_Index_Changed() takes no arguments (1 given)". If I remove the brackets after self.build_Layer_Index_Changed, I get no eror, but the function is not called (as the MessageBox doesn't appear. I thought I did it exactly as given in the Example at http://techbase.kde.org/Development/Tutorials/Python_introduction_to_signals_and_slots, but obviously I'm doing something wrong…


title: How to speed up a query for nearest linestring to point data

text: I have processed a Strahler stream order layer with a lot of river linestrings (type: Multilinestring) and their corresponding strahler number (see information here). Besides I have a point layer. I want to know which linestring is the nearest to the distinct points from the point layer to extract the strahler number for every point. ![enter image description here](http://i.stack.imgur.com/codah.png) My query here is very slow. For around 500 points and 2000 linestrings it takes over 6 minutes:               SELECT DISTINCT ON         (lname)         lnumber,         lname,         strahler,         min(ST_Length(geom)) AS distance     FROM (         SELECT              ST_MakeLine(ST_ClosestPoint(strahler_streams.geom, landslide.geom), landslide.geom) AS geom,             strahler,             lname,             lnumber     FROM landslide, strahler_streams) AS foo         GROUP BY strahler, lname, lnumber, geom          ORDER BY lname, distance;      I know there is the Indexed Nearest Neighbour Search in PostGIS. But when I try this...               SELECT         landslide.lnumber,         strahler_streams.strahler     FROM landslide a, strahler_streams b         WHERE lnumber=2114             ORDER BY b.geom <#> a.geom LIMIT 1      ... it only works when I select one single point (`WHERE lnumber=2114`) How can I speed up the query for the nearest neighbor search (linestring, point)?


title: From text data create a FeatureClass with table with python

text: I'm searching how the output of one program (text file with some writings, but also with data as numbers, dates, times, coordinates, all separated with commas) can be transformed to a Feature Dataset. **1) .txt -> .csv -> gp.TabletodBase** and then it's nearly done.. But get problems writing csv with Excel dialect - all row is writen in one cell. So, I create the delimiter ';' to get the data into cell by cell.               f= open(rawtable, 'r')     f2 = open(outtable, 'wb')     csv.register_dialect('dotcoma', delimiter=';')     writer=csv.writer(f2, dialect='dotcoma')          for line in itertools.islice(f, 6, 7): #reading labels that are in 6'th row         labels= line.strip('\n').split(",")         writer.writerow(labels)      But well, then Arcgis does not recognize the table.. I have been trying to cast the data, like integers, floats, strings.. but that did not help.. > **2) .txt -> .shp...** \----> THAT WORKS!!! The shapefile does not have a type of time.. So was skiping it. _EDITED_ However tried the Pyshp library. But after creating fields blindly I record the tuple of all data in a row and I'm getting an error: `IndexError: tuple index out of range`.. I simply write fields from the list of data, depending to which group they go:               for fieldname in labelkeys:         if fieldname in t_int:             w.field(fieldname, 'N', 8, 0) #1probl was here as I made '8' vs. 8 ! stpd         ...         else:             w.field (fieldname,'C', 25)      Another problem with coordinates. I was writing floats... and it did not like it.. So I had to write integers. Plus one cannot pass a tuple to w.record(). Thanks to StackExchange I've learned the benefits of asterix :) and one simply has to write w.record(*mylist). **3) .txt -> feature class...** Here I have problems writing the entire row at once..(like the list of values) into a NewRow. Since field names are made automatically I can not for each make like `row.fieldname=somevariable`. I have tried while looping to SetValues `row.setValue(labelkeys[unit], data[unit])` and then at the end of a list `rows.InsertRow(row)` , but it seemed that it sets just one value. Though I tried by hand and got an error in data casting.. The error comes that I read data like string '34', cast it to int(34) and try to write into field with type 'DOUBLE' (since I grouped some data and wanted to make more generic types).               if labelkeys[unit] in t_int:         if data[unit]:             row.setValue(labelkeys[unit], int(data[unit]))      So for now.. nothing is really working.. Maybe there is some other way I have not though about? If not, could anyone help to get the things work?...


title: Combine overlapping rasters with priority to values of specific rasters?

text: I have 5 rasters that cover a large area. I want to combine all 5 into 1 raster. Each raster has data for portions of the total area that none of the other rasters cover. Raster C overlaps A,B and D. While E overlaps D. No other rasters overlap each other. On the portions that are overlapping, only certain Rasters contain the data that is needed to be displayed. The following shows which Raster are needed to be displayed. Raster C > A Raster C > B Raster C > D Raster E > D I've tried `Mosaic to New Raster`, `Mosaic Dataset` and `Mosaic (Data Management)`, all unsuccessful. The result is only the first Raster. The Blank area has a Value of 0, which I believe is overlapping the other rasters value when I combine them. How do I display the Values > 0 when I combine the Rasters giving priority to the values of Raster C when combined with A,B and D and Raster E when combined with D? The edges of Raster D that overlap C and E have missing data as seen below. ![enter image description here](http://i.stack.imgur.com/IlaVG.png) ArcGIS 10.1, advanced arcinfo license. Edit: I've tried all the suggestions below and here are results: ![](http://i.stack.imgur.com/7bYtS.jpg) Above: This is Radouxju's suggestion using raster calculator with IsNULL, the resulting image is Raster C and I tried to overlap only D, D's area is not present. It is an exact image of Raster C. I also tried Raster to New Mosaic with FIRST and LAST, the resulting image was the same as above, but with colour. ![enter image description here](http://i.stack.imgur.com/Ot9GH.jpg) Above: This is Cody's suggestion using Mosaic Dataset and setting priorities, as you can seen the 0 Values overlap the other Rasters depending on their priority. This image is an example/ test to visualize what happens, the priorities here are B > CE > AD.


title: How to compile the source code of module MOD14 of NASA

text: I'm working in a project which uses an Active Fire Detection base on MODIS (MOD14) of NASA. The source code and compiled-file can be downloaded at NASA Now I want to compile the source code, although the runable module is already provided. Have anyone ever working on stuff like this?


title: Inserting strings in language other than English into MariaDB with GeoNetwork?

text: I am new to **GeoNetwork** and I testing how to insert a text string into **MariaDB** using a language other than English. The database is showing the inserted text string as `??????????`. I have set XML encoding to `UTF-8` AND MariaDB to `utf8_general_ci`. Can you help me to fix this problem? ![enter image description here](http://i.stack.imgur.com/wbDEp.png) ![enter image description here](http://i.stack.imgur.com/JGT9t.png)


title: how can i convert an excel table to a shapefile with correct scale in arcmap?

text: i know that there is already one similar question to this one but my problem is different and none of the provided answers did resolve my problem. i've an excel table with some rows and in each row there is the identification and the xy coordinates of an especific point. the xy coordinates were aquired with a gps and are in wgs84. the data layer is in wgs84 too. my problem is that every time and every way i try to import this data to a shapefile, the resulting points are perfect in form but at a ridiculous scale (1:0.89)... exporting the table to dbf and using it gave the same result... how can i make the data import in a way that the scale is correct? thanks in advance =^) \-----EDIT----- if i use the "go to xy" tool to enter the coordinates of two points manually, using the values of x and y in the excel/dbf file, the two points appear in the right places and at the right distance (which implies a correct scale)... \-----EDIT----- excel file: http://dl.dropbox.com/u/3639616/xls%20to%20shp.xls the farthest points are codigo=394 and codigo=297, with about 26.3 km between them.


title: Any difference between EPSG:3003 and EPSG:102091?

text: I am working on QGIS 2.2.0. The PRS I am working on is Monte Mario / Italy zone 1 (EPSG:3003), with the following parameters:               PROJCS["Monte Mario / Italy zone 1",     GEOGCS["Monte Mario",         DATUM["Monte_Mario",             SPHEROID["International 1924",6378388,297,                 AUTHORITY["EPSG","7022"]],             AUTHORITY["EPSG","6265"]],         PRIMEM["Greenwich",0,             AUTHORITY["EPSG","8901"]],         UNIT["degree",0.01745329251994328,             AUTHORITY["EPSG","9122"]],         AUTHORITY["EPSG","4265"]],     UNIT["metre",1,         AUTHORITY["EPSG","9001"]],     PROJECTION["Transverse_Mercator"],     PARAMETER["latitude_of_origin",0],     PARAMETER["central_meridian",9],     PARAMETER["scale_factor",0.9996],     PARAMETER["false_easting",1500000],     PARAMETER["false_northing",0],     AUTHORITY["EPSG","3003"],     AXIS["X",EAST],     AXIS["Y",NORTH]]      A dataset I need to use on the above PRS is in Monte_Mario_Italy_1 (EPSG:102091), with the following parameters:               PROJCS["Monte_Mario_Italy_1",     GEOGCS["GCS_Monte_Mario",         DATUM["Monte_Mario",             SPHEROID["International_1924",6378388,297]],         PRIMEM["Greenwich",0],         UNIT["Degree",0.017453292519943295]],     PROJECTION["Transverse_Mercator"],     PARAMETER["False_Easting",1500000],     PARAMETER["False_Northing",0],     PARAMETER["Central_Meridian",9],     PARAMETER["Scale_Factor",0.9996],     PARAMETER["Latitude_Of_Origin",0],     UNIT["Meter",1],     AUTHORITY["EPSG","102091"]]      Is there any difference between the two PRS? (I see a difference in decimals in parameter UNIT["Degree", ...]). In other words, do I need a re-projection or datum transformation to analyse the data together (e.g. buffering, calculating distances etc) ? Thanks in advance. Alessandro


title: The maximum record length exceeded?

text: I downloaded some census tract data which are stored in .gdb. I am able to open them in ArcGIS as shpfiles. I selected some tracts and tried to export them into another folder. But I kept getting this error: "Error exporting data. The maximum record length has been exceeded". Is this because the attribute table has too many columns? What is the solution?


title: Distance between latlongs varies in different projection systems?

text: I have shapefiles with correct projections. The SRS is EPSG 32735 (WGS 1984 UTM 35S) Now, I have these two coordinates from this projection system: (31.3751,-23.4388) and (31.3779,-23.4355). Obviously the coordinates (format and value) will vary on different projection systems, and there in lies my confusion. I used the following link to calculate the distance between the two points: http://www.movable-type.co.uk/scripts/latlong.html I copy pasted the lat/longs exactly as shown above and the distance it gives me is 0.313 meters. But I'm not sure if this is correct at all, since the calculator didn't ask me for any projection system? So:   1. Does this mean the projection doesn't matter to calculate the distance?   2. Does it differ between different projection systems?   3. Is the tool I'm using incapable? Thanks, Kaushik


title: How to create generic building footprints based on attribute information?

text: I'm trying to create a tool that generically creates building footprints on a given parcel according to certain attribute information of the parcel like building density, max. height, utilisation, etc. Is there a way to generically create polygons in model builder with specified dimensions? Hints and suggestions are much appreciated.


title: Help with Turn-feature-class

text: I aleady have a street network dataset. I have create a turn-feature-class in a feature dataset. I want to add the turn-feature-class to the network dataset. How do I do that in ArcCatalog? Many Thanks


title: raster statistics - regressions between raster data

text: I have the following problem: I have two set of rasters, with data that should correlate...one are measured data, one are modelled data. The correlation is weak, however, when looked at with a normal regression analysis. I suspect that different parts of the base terrain need different modelling input, or that in some parts there is a shift between the correlating cells. So now I try to visualize where these residuals are highest between my two rasters, and where the real outliers are and whether they are in similar terrain. Is there any tool that would allow me to directly compare the rasters in GIS? I run ArcGIS 9.3 with Spatial Analyst and Statistics, and also have QGIS and GRASS on my Linux-based machine. What would be easiest? I am not trying to do a full-on analysis, I am just trying to see where the odd data are and if there is a systematic spacial pattern behind them. Cheers, Peter


title: Debugging SQL in arcpy.da.SearchCursor RuntimeError: An invalid SQL statement was used?

text: I am going through a python script in arcpy 10.1 to Select a field value (i.e. NEAR_DIST) based on three condition as                       1) NEAR_DIST<=radius(i.e.500 or anything from user input)              2) NEAR_DIST<>0               3) NEAR_DIST<> -1      I am using 500 which is to be variable asked from user and to be saved in "radius" variable so, i am only using "radius" variable in the SQL statement in following code-                   import arcpy,os,sys         feature =     "C:\Users\Bakul\Desktop\Od\example_shapefile\Polyline_segments_example.shp"         radius = 500         # create a dictionary to hold the list of near ids for each input id         nearest_dict = dict()         arra = []         whereClause = '''"NEAR_DIST" <= '%{x}%' AND "NEAR_DIST" <>'%{y}%' AND "NEAR_DIST"   <>'%{z}%\''''.format(x=radius,y=-1,z=0)         print radius         with arcpy.da.SearchCursor(feature, "NEAR_DIST", whereClause) as rows: ##NEAR_DIST  is a field, I need sort all NEAR_DIST <=500 (i.e radius) and not equal 0 and -1             for row in rows:## debugger does not pass here,                  nearest_id = row[0]                 arra.append(nearest_id)             del rows             del row         print arra      Debugger Error is                   500         Traceback (most recent call last):           File "C:\Users\Bakul\Desktop\Od\Python\scratch.py", line 89, in <module>             for row in rows:         RuntimeError: An invalid SQL statement was used.      Can you tell where is the problem (i think in SQL whereClause)?


title: Can QGIS open a SDF file created by HEC RAS?

text: I used the "export GIS data" option in HEC RAS and got a SDF file. Is it possible to open it in QGIS?


title: Overview effect on postgis polygon table

text: Is There a Way, with postgis and qgis, to reproduce the overview effect of tiff Raster on a polygon table in order to increase speed? What if we could create simplified versions of this table,gathered under the same view. Until now, i only play with different scale visibility on different simplified version of the table... Best regards


title: Intersecting polygons in Oracle

text: Have a polygon in table 2 with its id, no and geometry. In table 1 I have polygons with the same fields. The polygon in table in intersects with a few polygons from table 1.What I am trying to do is clip of the polygons which overlap polygon from table 2 and insert the same fields and geometry in table 3               insert into table 3       select a.store_id,b.store_id,a.store_number,a.client_id,sdo_geom.sdo_intersection(b.geometry,a.geo metry,0.005)       from table_1 a, table_2 b       where b.store_id=34746       and sdo_anyinteract(b.geometry,a.geometry)='True';        When I do run this I dont get any output. Where would I be going wrong


title: Insert Vertex into Polyline based on Point layer using arcpy

text: This is essentially the same question asked here that was not answered... Using ArcGIS 10.1 Advanced I have a point layer with points that represent the intersection of a driveway with a street centerline. There is not always a vertex on the street centerline where the driveway intersects, so I am looking for a solution (preferably using arcpy) that loops through each point and inserts a vertex on the line that intersects with that point. I do not want to use Split Line at Point method because I don't want to split the line. Looking for any ideas on how to accomplish this...


title: OpenLayers: draw path with only two points

text: I want to draw a path with no more then 2 points. On the second point I want the path to be finished. I'm using the OpenLayers DrawFeature, I set a callbacks for "point" where I count the created points and if the points are two I'm finishing the sketch with `finishSketch()` method. The problem is that when finishSketch() method calls the `finishGeometry()` method which removes the last added point from the path which cause the following error to show `"Uncaught TypeError: Cannot read property 'geometry' of null "` My question is the any legit way for finishing a path programmaticaly or this is a bug in the DrawFeature control.


title: Partially fill a sink - ArcGIS

text: I am trying to fill sinks but only to a certain level. For example - if I want to fill them by 5 meters - this means that a sink that is 15 meters deep - will fill only to 10 meters and a sink that is 3 meters deep would fill completely. Using ArcGIS 10.0 - service pack 4. Any help of this would be appreciated. I thought I created a script to do what I wanted, but my client clarified something different than what I though. Attached is an image of what I have and what I want for a custom fill tool.![enter image description here](http://i.stack.imgur.com/MjEVl.png)


title: Split polylines at polygon boundaries in QGIS

text: I have a set of contour lines and a polygon set containing glaciers. Because I want the glaciers to have blue contours, I need to split the contour lines at the glacier boundaries. I can achieve this using the following method: 1) Intersect the contours with the glacier polygon to create the glacier contours. 2) Create a spatial difference between the glacier polygon and an extent polygon. 3) Intersect the contours with the spatial difference from (2) to create the blue conoturs. 4) Merge the two results together if I want it in a single data set. But the following method would be more convenient, and usable for other types of analyses as well. 1) Combine a polyline dataset (or point dataset for that matter) with a polygon dataset. The lines with be split at the polygon boundaries and each output line (or point) will take the attributes from the polygon it is inside. Such a function would make my problem so much simpler to solve, as I can do it in one single operation. Is there a function for this in QGIS?


title: How to add feature class to MXD with ArcPy (Python)?

text: I am using ArcGIS 10.1 and I am trying to add a feature class that I create to an MXD, so that I can join to it. This all needs to be done using Python. Is this possible?


title: Japan Data: Postcode/Census

text: Hi I'm primarily looking for a GIS data source for postal codes in Japan as well as census data for different scales of geography. Has anyone found a good source of these data sets that is in English? Thank you for your help.


title: Retrieving the proportion of specific pixel class from CHM rasters

text: CHM = canopy height model I am very new to R (specially for GIS purposes) and I am struggling a bit in getting the following pixel statistics:   1. Proportion of zero value pixels = [nº of zero value pixels/(total nº of pixels - nº of NA pixels)]   2. Proportion of pixels smaller than 3rd quantile = [nº of <3rd quantile pixels /(total nº of pixels - nº of NA pixels)] Here is a reproducible example:               library(raster)       f <- system.file("external/test.grd", package="raster")       r <- raster(f) #r is the object       plot(r)       summary(r)


title: How to loop through a list of shapefiles and rename them?

text: I am trying to upgrade a simple script to list all shapefiles in a folder and if necessary in any folders within and then loop through the list of shapefiles and rename them.               import arcpy          # Set work environment     arcpy.env.workspace = "C:\\Users\\Rachael\\Documents\\Export"          # List shapefiles in folders     shp = arcpy.ListDatasets()          try:       for shp in arcpy.ListDatasets():         if arcpy.Exists("Histori2.shp"):           arcpy.Rename_management("Histori2.shp", "H_BEV_CANS.shp")


title: openlayers multiple layers with different projection

text: My page contains two layers: 1. Open Street Map as a basic layer and WMS layer from mapserver as a overlay layer. My map show this two layers but in a wrong position. I spent a lot of time trying everything but it doesn't work. I changed mapfile and file to EPSG4326 projection but it doesn't work to. This is my mapfile               MAP     IMAGETYPE   "png24"     EXTENT     20.682299 52.148211 21.078846 52.337025     UNITS DD     SHAPEPATH "C:\ms4w\Apache\htdocs\projekt3\dane"     SIZE 800 600     DEBUG ON     FONTSET "C:\ms4w\Apache\fonts\fonts.txt"     PROJECTION     'proj=longlat'     'datum=WGS84'     'no_defs'     END     IMAGECOLOR 255 255 255     IMAGEQUALITY 95     IMAGETYPE png          OUTPUTFORMAT     NAME png     DRIVER 'GD/PNG'     MIMETYPE 'image/png'     IMAGEMODE RGBA     EXTENSION 'png'     END          SYMBOL     NAME "circle"     FILLED true     TYPE ellipse     POINTS 1 1 END     END                    LAYER          NAME kilometraz_popr     DATA kilometraz_popr     STATUS ON     DUMP false     TYPE POINT     TRANSPARENCY 100      EXTENT 458986.052854 392913.684787 757412.357297 535008.411398     PROJECTION     'proj=tmerc'     'lat_0=0'     'lon_0=19'     'k=0.9993'     'x_0=500000'     'y_0=-5300000'     'ellps=GRS80'     'towgs84=0,0,0,0,0,0,0'     'units=m'     'no_defs'     END     CLASS         STYLE          SYMBOL "circle"           SIZE 7.0           OUTLINECOLOR 0 0 0          COLOR 90 23 223         END      END     END          LAYER          NAME etykiety     DATA kilometraz_popr     STATUS ON     TYPE annotation     PROJECTION       "init=epsg:2180"     END     LABELITEM OPIS_KM     TRANSPARENCY alpha     CLASS         LABEL             COLOR    0  0  0             FONT  arial             #BACKGROUNDCOLOR 255 255 255             ENCODING "CP1250"             TYPE  truetype             SIZE  7             POSITION  AUTO             PARTIALS  FALSE             BUFFER  5             OUTLINECOLOR  255  255  255             ANTIALIAS true         END      END     END     END      My index.html file is below:               <html lang='pl'>     <head>     <meta charset='utf-8' />     <title>Wawa</title>     <script type='text/javascript' src='OpenLayers.js'></script>     <script src="proj4js-combined.js"></script>      <script type='text/javascript'>     var map;     function init() {             epsg4326 = new OpenLayers.Projection("EPSG:4326");             epsg900913 = new OpenLayers.Projection("EPSG:900913");             epsg2180 = new OpenLayers.Projection("EPSG:2180");     var option = {         projection: new OpenLayers.Projection("EPSG:900913"),         displayProjection: new OpenLayers.Projection("EPSG:4326")     };            map = new OpenLayers.Map('mapa', option);          olmapnik = new OpenLayers.Layer.OSM("OpenStreetMap Mapnik",    "http://tile.openstreetmap.org/${z}/${x}/${y}.png");     map.addLayer(olmapnik);     map.setBaseLayer(olmapnik);                   var warstwa_wektor = new OpenLayers.Layer.MapServer("OpenLayers Mapserver", "http://localhost/cgi-bin/mapserv.exe", {map: 'C:/ms4w/map/map_test1.map', layers: 'kilometraz_popr', srs:'EPSG:4326', projection:epsg900913,     isBaseLayer: false,         visibility: true, transparent: true, format:'image/png'}, {opacity: 1});     map.addLayer(warstwa_wektor);          map.addControl     (new OpenLayers.Control.Scale());     map.addControl(new OpenLayers.Control.Navigation());     map.addControl(new OpenLayers.Control.MousePosition());     map.addControl(new OpenLayers.Control.PanZoomBar());     map.addControl(new OpenLayers.Control.SelectFeature());     map.addControl(new OpenLayers.Control.LayerSwitcher());     map.addControl(new OpenLayers.Control.Permalink());          extent = new OpenLayers.Bounds(19.50,51.50,21.50,52.7).transform(new  OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));     map.zoomToExtent(extent,5);     }          </script>     </head>          <body onload='init();'>     <div id='mapa' style='width: 800px; height: 400px;'>     </div>     </body>     </html>      Please tell you me my mistake. * * * Yes I know what I should do. I really greatfull for your yelp. My main mistake was projection. I use 3857 and everything works. Thank you very much


title: Why do ECWs go missing when exporting MXD to PDF with ArcPy?

text: I'm trying to export some MXD maps to PDF with arcpy (ArcGIS 10.1). Everything works fine and dandy except for ECW images, which don't appear in the PDF. However, when I execute the same code from the ArcGIS python console, the ECW is exported. So if I execute this from an external script:               import arcpy     mxd = arcpy.mapping.MapDocument("path\\to\\test.mxd")     arcpy.mapping.ExportToPDF(mxd, "path\\to\\test.pdf")      I don't get the ECW layer in the PDF. But if I execute the same lines in the console, it works and the ECW layer is exported properly. Also, if I run               execfile('export_mxd_to_pdf.py')      from the ArcGIS console, I also get the ECW in the PDF. I don't have this problem with other image formats, like JPG. What's going on here? What's different between the environments? How can I debug this?


title: OSRM API - nearest point fed into routing

text: I am using the OSRM API to (1) determine the nearest point on a road segment from a give random point with               http://server:5000/nearest?loc=_lat_,_lon_      The returned coordinates are used then as an input to (2) determine the route between this point and an other give point on a road segment with               http://server:5000/viaroute?loc=_lat1_,_lon1_&loc=_lat2_,_lon2_      Every now and then the point determined in step (1) is located on a non routable road (e.g. private road). I am looking for a way to work around this problem. a) Is there a way to determine in step (1) the next closest point on a different road segment or only query road segments on routable roads. or b) Is there a way that the routing query in (2) automatically determines the next closest routable point? I am also open for different solutions. However I am not interested in setting up my own server with pgrouting or similar services. I like to solve it with an open source routing API (preferable with the OSRM one). The request is done from a python script (requests and json).


title: Google Maps/Earth becoming GIS solutions yet?

text: I've been researching internet lately to make a comparaison between **Google Geo** solutions ( Google Maps, Google Earth) and other **GIS solutions** (Free of charge or commercial, QGIS, ArcGIS, etc.). I can say I have now an idea about how they all compare, still, I wish to confirm my conclusions by whatever your experience have served you, especially of those who dealt with Google Geo in their projects. My question is : is Google Geo (let me underline **_Google Earth_** here) becoming a GIS solution ? How are **_geoprocessing tools_** and **_workflows_** in the latest versions of Google Earth ? What does it offer in term of **_geo-analytical capabilities_** ? Edit : I forgot to mention that I am well aware of this thread, I just need extra infos about GIS-like capabilities of Google Geo solutions.


title: GeoServer + PHP + Single Sign On

text: I'm developing a web app which uses GeoServer side by side with a PHP framework. I've been struggling for a while trying to find out an acceptable way of implementing a Single Sign On. I would like to centralize authentication in the web app so the user isn't asked for a second login when accessing the GeoServer admin interface (as I prefer instructing the client on how to use the admin interface rather than building a new one to achieve a subset of the functionality). In my mind I imagine to proxy Tomcat through Apache/PHP. I guess my plan implies creating a wrapper with my PHP framework around GeoServer/Tomcat. What I would like to know is:   * isn't there anything ready for wrapping GeoServer with PHP?   * is the approach I described a good one?   * would it be better to implement something like O-Auth?   * should I eventually leave some routes (say WMS ones, which mainly request images) bypass the proxy (say via an Apache Reverse Proxy)?   * in the case of completely wrapping WMS what would be a good approach to minimize performance impact (speed, cpu and ram overhead) using a php proxy? Any suggestion would be greatly appreciated! This question is strictly linked to my previous question.


title: QGIS Scripting problem updating contrast enhancements for raster layers using bandStatistics

text: Found this function to update statistics for a raster layer :               rlayer.setContrastEnhancement(QgsContrastEnhancement.StretchToMinimumMaximum,                               QgsRaster.ContrastEnhancementMinMax)      Which works great for minmax and cumulative but stddev is hardcoded to 1 standard deviation. For cumulative it reads the lower and upper values from the settings but no setting for stddev. So I have looked at the C++ code in the api and tried to convert that into a QGIS script. The only problem is that the band statistics and cumulative cut functions are returning "nan". Take a look at my code segment below and see if you can help:                      myEnhancement = QgsContrastEnhancement(myType)            myEnhancement.setContrastEnhancementAlgorithm(ContrastEnhancement,True)                 if Limits == 0:               myRasterBandStats = rlayer.renderer().bandStatistics( myBand,                                          QgsRasterBandStats.Min | QgsRasterBandStats.Max)               myMin = myRasterBandStats.minimumValue               myMax = myRasterBandStats.maximumValue                 elif Limits == 1:                 myRasterBandStats = rlayer.renderer().bandStatistics( myBand,                                        QgsRasterBandStats.Mean | QgsRasterBandStats.StdDev)                 myMin = myRasterBandStats.mean - ( StdDev * myRasterBandStats.stdDev )                 myMax = myRasterBandStats.mean + ( StdDev * myRasterBandStats.stdDev )                 elif Limits == 2:                   myMin,myMax = rlayer.renderer().cumulativeCut( myBand, CumulativeLower,                                                               CumulativeUpper)                 raise GeoAlgorithmExecutionException( str(myBand) + " " + str(myMin) + " " + str(myMax))            myEnhancement.setMinimumValue( myMin )            myEnhancement.setMaximumValue( myMax )                 if renderType == "singlebandgray":               rlayer.renderer().setContrastEnhancement(myEnhancement)            elif renderType == "multibandcolor":               if Band == 0:                  rlayer.renderer().setRedContrastEnhancement(myEnhancement)               elif Band == 1:                  rlayer.renderer().setGreenContrastEnhancement(myEnhancement)               elif Band == 2:                  rlayer.renderer().setBlueContrastEnhancement(myEnhancement)            rlayer.triggerRepaint()      * * * Heres the full script using the first function which works great but does not enable you to have more than 1 standard deviation               ##[Example scripts]=group     ##Limits=selection "MinMax";"StdDev";"Cumulative"     ##Stretch=selection "NoStretch";"StretchToMinMax";"StretchAndClipToMinMax";"ClipToMinMax"          import os     from qgis.core import *     from qgis.gui import *     import qgis.utils          layers = qgis.utils.iface.legendInterface().layers()          for layer in layers:         layerType = layer.type()         if layerType == QgsMapLayer.RasterLayer:            if Stretch==0:               ContrastEnhancement = QgsContrastEnhancement.NoEnhancement            elif Stretch==1:                 ContrastEnhancement = QgsContrastEnhancement.StretchToMinimumMaximum            elif Stretch==2:                 ContrastEnhancement = QgsContrastEnhancement.StretchAndClipToMinimumMaximum            elif Stretch==3:                 ContrastEnhancement = QgsContrastEnhancement.ClipToMinimumMaximum                 if Limits==0:               layer.setContrastEnhancement(ContrastEnhancement,QgsRaster.ContrastEnhancementMinMax)            elif Limits==1:               layer.setContrastEnhancement(ContrastEnhancement,QgsRaster.ContrastEnhancementStdDev)            elif Limits==2:                layer.setContrastEnhancement(ContrastEnhancement,QgsRaster.ContrastEnhancementCumulativeCut)            layer.triggerRepaint()


title: gdal read whole geotif band into vector

text: i have another c++ GDAL question. I try to read alle elevation data out of an GeoTiff file into a c++ vector so that i can get elevation at x,y like this: vector[x][y]. I also like to use this data to process a new image out of this data. Getting the data with gedal is quite simple but it seems like it is wrong ordered. If i bake the mapped infos straight to an image it is rotated and mirrored. Is there a way to sort the data in a different way in GDAL? How can I get my data right? I am bit stuck. :(               poBand->RasterIO(   GF_Read, 0, 0, nXSize, nYSize,                              pafScanline2, nXSize, nYSize, GDT_Float32,                              0, 0 );         int totalCount = 0;         for (int i = 0; i < poBand->GetXSize(); i++){                  vector< float > yVal;             for (int j = 0; j < poBand->GetYSize(); j++){                 yVal.push_back(pafScanline2[totalCount]);                 totalCount++;             }             altData.push_back(yVal);         }      Here is the gdalinfo output:               Driver: GTiff/GeoTIFF     Files: test5.tif            test5.tif.aux.xml     Size is 6001, 6001     Coordinate System is:     GEOGCS["WGS 84",         DATUM["WGS_1984",             SPHEROID["WGS 84",6378137,298.257223563,                 AUTHORITY["EPSG","7030"]],             AUTHORITY["EPSG","6326"]],         PRIMEM["Greenwich",0],         UNIT["degree",0.0174532925199433],         AUTHORITY["EPSG","4326"]]     Origin = (74.999583696539048,40.000416763514124)     Pixel Size = (0.000833333333333,-0.000833333333333)     Metadata:       AREA_OR_POINT=Area     Image Structure Metadata:       INTERLEAVE=BAND     Corner Coordinates:     Upper Left  (  74.9995837,  40.0004168) ( 74d59'58.50"E, 40d 0' 1.50"N)     Lower Left  (  74.9995837,  34.9995834) ( 74d59'58.50"E, 34d59'58.50"N)     Upper Right (  80.0004170,  40.0004168) ( 80d 0' 1.50"E, 40d 0' 1.50"N)     Lower Right (  80.0004170,  34.9995834) ( 80d 0' 1.50"E, 34d59'58.50"N)     Center      (  77.5000004,  37.5000001) ( 77d30' 0.00"E, 37d30' 0.00"N)     Band 1 Block=6001x1 Type=Int16, ColorInterp=Gray       Min=1055.000 Max=8569.000        Minimum=1055.000, Maximum=8569.000, Mean=2993.143, StdDev=1711.792       NoData Value=-32768       Metadata:         STATISTICS_MAXIMUM=8569         STATISTICS_MEAN=2993.1430406214         STATISTICS_MINIMUM=1055         STATISTICS_STDDEV=1711.7917180302


title: Create a point C from point A in the direction of point B

text: I am very new to QGIS and GIS in general. I hope someone could help me figure this out. I have pairs of coordinates which signify the start and end of a transect line. I am trying to adjust the end point such that it is x meters from the starting point. What I have:               Coordinates of Points A and B.      What I want to accomplish:               Create a point C such that it is x meters from point A in the direction of B.     or     Move point B such that it is x meters from point A in the same direction.      Of course I could do this manually since I have just more than a dozen points, but I am trying to be a bit more accurate... Thank you very much for your help!


title: Proxy to work with OpenLayers and Geoserver Locally with IIS express

text: I need to set a proxy to work with OpenLayers and Geoserver, which are both using a different port on my pc. I'm the site with openlayers is hosted **IIS express** and I'm not using Apache (so the proxy.cgi solution is not for me). I need it to perform a wmsgetfeatureinfo and, at this point, I really don't know what to put in `OpenLayers.ProxyHost =`. Can you help?


title: Place arcgis layer on top of google map

text: This is my old code that loaded arcgis layers on top of the arcgis map:                   function init(mapServer) {                 map = new esri.Map("map", {                     basemap: "topo",                     //center: [25.45, 45.42],                     zoom: 3                 });                 layer = new esri.layers.ArcGISDynamicMapServiceLayer(mapServer);                 layer.on("load", displayLayerExtent);                    map.addLayer(layer);             }         function displayLayerExtent(evt) {                 map.setExtent(evt.layer.initialExtent, true);             }      where mapserver is the url to a mapserver, I used the test map servers offered by arcgis. But the client doesn't want the map behind the layers to be the arcgis map, but the google map. How can this be solved?


title: How to Compress/Stack Landscape features to examine spatial distribution and fragmentation?

text: I was wondering if anyone could help me figure out how to compress landscape information from multiple rasters (or polygons) into a single output. Basically, I want to explore how landscape characteristics affect the spatial distribution of a population. As of now I envision a raster with unique values for each combination of slope, aspect, forest type, elevation, and population number (all variables separately re-classified then compressed into one output). Does anyone have any ideas? I don't think this is the same 'stacking' as making a multi-band mosaic because all the information covers the same area...this is not jigsaw-puzzling separate locations together, they all have the same surface area and shape. I could be wrong though. Regarding the statistics, I was also wondering if after compressing the information into one raster/polygon, I can look at the population as a whole rather than randomly sampling from it across the landscape because each rastercell has a value. This is also time sensitive (due next week!) so any guidance would really help!


title: Mercator projection image data

text: I'm trying to create a map which I can place points on within my application I found some mecator jpg images but in order to accurately place points I need to convert the lat/long to pixels so that the data is accurately displayed on screen. I've been following this tutorial http://jkfill.com/2011/03/18/placing-geo-coordinates-on-an-image/ When trying to calculate the globe width and centre point of the map I simply took the pixel width and the center point of the pixels. i.e if the image is 100 x 100 then the globe width is 100 and the center point (50,50) although I now think this is incorrect? How can I calculate these points? Alternatively can anyone recommend map images with this data available? Thanks


title: ArcGIS Viewer for Flex providing login

text: I have ArcGIS server and ArcGIS Viewer for Flex installed up and running. I use ArcMap 10.2 to publish my shapefiles. There is a lot of users who is using my maps. How can I provide login for them in the web page for map? Should I use it as a widget or there is another way?


title: Raster on-the-fly projection in ArcGIS has different resolutions

text: I have a map document authored in ArcMap. The spatial reference for each of the raster layers is Alaska State Plane. If I set the Data Frame projection to Web Mercator (Auxiliary Sphere), then my raster layers have a resolution of about 2 meters, when the actual data has a resolution of about 2 feet. I'm trying to publish the map document as a Map Service (ArcGIS Server 10.0) but the quality of the imagery is not acceptable for the use case. Is there any way to make the imagery draw at 2-ft. resolution, or do I need to reproject the data (about 500 GB of raster imagery)?


title: Annotation in ArcView 3.2

text: I've been asked to take a parcel layer (shapefile) that contains house numbers and create an annotation layer so the house numbers are displayed in a separate annotation layer rather than a label to the parcel layer. They'd like the ability to move and rotate the house numbers for readability. Is there a way to create an annotation layer for ArcView 3.2 that will do this? I've experimented with exporting a couple different ways and it loses it's usability once it comes outside the geodatabase in which I originally created it. Will ArcView 3.2 even read a geodatabase?


title: SQL Server data modelling + spatial data = what book?

text: What recommended book can you provide about data modelling for geographical spatial data (geometry) in SQL server?


title: Open street map overlay white background

text: I have been working through trying to get the NOAA ridge radar image to act as an overlay. Looking at other entries on the site I have worked with the image and have it projecting correctly but the png has a semi opaque white background I can't seems to get rid of. I have tried setting transparent = true it appears that the warp is adding the white background               gdalwarp -t_srs '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs' radartemp1.png radartemp2.png      i tried also adding -dstalpha but then the layer comes out black... anyone know how to solve this problem


title: Overlaying a mbtiles tile layer produced with tilemill over another layer maintaining transperency

text: So using leaflet, I'm overlaying my base tile layer as google satellite imagery. I want to overlay another layer (Mbtiles layer) produced with tilemill over this satellite layer. I was successful in doing this. I set the background of the mbtiles file as transparent on tilemill. Im using tilestache to serve the tiles and leaflet for a JS API. Now, I want this mbtiles layer to only show the coloured layers (roads, polygons etc.) and the background to be transparent (to show the google satellite layer wherever the mbtiles layer is transparent). But this is not happening. The background of the MBtiles layer is white and there is no transparency. Where am i going wrong?


title: How to convert Geometry type column to degree Minutes Second in SQL Server 2008 R2

text: I have shape column of type geometry and I want to convert X and Y coordinates into Degree Minutes Second format using SQL query and i am using sql server 2008 R2. Currently I got Lat and long using this query               Select W.SHAPE.STY as latitude, W.SHAPE.STX as longitude from tbdata


title: How to merge overlapped shapefiles without one layer's overlapping area?

text: I created a forest type map from 50 scene satellite images as shapefile under the following conditions.   * Some shapefiles are overlapped, because scene overlaps each other.   * There is an order of shapefiles in the table of contents.   * Merged file has to be follow the order. It means below shape file is not necessary in the overlapping area. It is possible merging files, if I clip one by one, even though it takes a time. Please let me know, if someone knows a more efficient way to resolve this issue.


title: Trouble formatting SQL expression for where clause in "Search layer by attribute" arcpy function

text: I am trying to pass a variable into a where clause in the "Search Layer by Attribute" arcpy function. I cannot seem to figure out the correct format for the SQL expression so that the function reads the value of the given variable, and doesn't just read the variable name as a string. The variable I am trying to pass is called "currentRaster" and is a list of dates in the 'YYYY-MM-DD 00:00:00' format. The expression I am having trouble with is: whereClause = """"DATE" = date 'currentRaster'""" If I run the python script and tell it to "print whereClause", the result is: whereClause = """"DATE" = date 'currentRaster'""" instead of, for example: whereClause = """"DATE" = date '1998-01-31 00:00:00'""" If I substitute one of the dates for the variable name, e.g. 1998-01-31 00:00:00, and run the script, the expression works. Does anyone know the proper way to format this SQL expression so that it reads the variable as the date it contains? The rest of the script it pasted below. Thank you! Liza               #Interpolate points from temporal ArcGIS rasters          #import modules      import arcpy     from arcpy import env     import datetime     from datetime import timedelta          #environment settings     env.workspace = "Z:\\sharks\\Models\\Summer2013\\Data_From_Jason\\compiled Chl"     arcpy.env.overwriteOutput = True          #enable MGET modules and enable verbose logging     import GeoEco.Logging     from GeoEco.Logging import Logger     Logger.Initialize(activateArcGISLogging=True)     import GeoEco.SpatialAnalysis.Interpolation     from GeoEco.SpatialAnalysis.Interpolation import Interpolator          #set up variables     rasters = arcpy.ListRasters()     points = "Z:\\sharks\\Models\\Summer2013\\summer2013.gdb\\bll_wgs84"     fields = ["Chla_Daily"]     method = "Automatic"          #define function that extracts the gregorian date from the raster filename (which is in julian days)     def get_date(fname, nbeg, nend):         yr = int(fname[nbeg:nbeg+4])         j = int(fname[nbeg+4:nend])         beg = datetime.datetime(yr,1, 1)+ timedelta(days=(j-1))         return beg          #create feature layer for select layer by attributes analysis     arcpy.MakeFeatureLayer_management(points,"points_lyr")          #set up loop to step through all rasters     for i in rasters:         currentRaster = get_date(fname=i, nbeg=16, nend=23)         #set up cursor to step through all records in attribute table         with arcpy.da.SearchCursor("points_lyr", ("Date",)) as cursor:             for row in sorted(cursor):                 value = row[0]                 print value                 print currentRaster                 #if this record's date matches the date of the current raster, then select that record and perform interpolation                 if value == currentRaster:                     whereClause = """"DATE" = date 'currentRaster'"""                     print whereClause                     arcpy.SelectLayerByAttribute_management("points_lyr", "NEW_SELECTION",whereClause)                        Interpolator.InterpolateArcGISRasterValuesAtPoints([i], "points_lyr", fields, method,None, None,False,None,None,None,None)                 else:                     print "no"          arcpy.Delete_management("points_lyr")     del row, cursor


title: Interpolation of multibeam bathymetry

text: I'm trying to make a 3D map with multi-beam bathymetry (65536 features, xyz), but it seems that the interpolation methods are too rigid. I need to smooth the data so it looks more natural and realistic. Which parameters do I have to vary so I can get this?


title: How to update Bing base map data?

text: Can someone tell me how, or if it's possible, to send street centerline updates to Bing for inclusion in their base maps.


title: ITopologicalOperator simplify takes more time

text: The objective is to cut the polygon where the lines are intersect. The polygon has 183182 vertices and taking more time while using simply method so process time is increased very huge. Please suggest the best process to cut the polygon.               ITopologicalOperator2 pTopologicalOperator = PolygonFeature.Shape as ITopologicalOperator2;          pTopologicalOperator.IsKnownSimple_2 = false;     pTopologicalOperator.Simplify();          IFeatureEdit featureEdit = PolygonFeature as IFeatureEdit;     ISet newFeaturesSet = featureEdit.Split(NewPolyline);


title: Searching 'select by color' for raster or map data

text: I'm planing to create an online map component with which the user can create/draw polygons based on underlying raster data. Is there some library / tool / procedure existing which enables the creation of polygons based on the colors of the raster data (or map data) - something like the 'select by color' in photoshop? Example: the user clicks in the center of a lake and the tool creates a polygon which covers the complete lake at the coastline.


title: ArcGIS (ArcObjects) Windows 7 - 64 bit -support

text: I have been developing with ArcGIS Cross Platform C++ SDK (with Qt 4.8.2) on 32 bit Windows XP SP3 (development machine). I used to use MapTool Control within ArcObjects on the development machine. Moreover I compile my application with Visual Studio 2010 in both release and debug mode. I installed Engine Runtime on Windows 7 (both 32 bit and 64 bit) on the test machines, however my application simply crashes on both machines whenever ` new QAxCtl(AoPROGID_MapControl, &mapctl, Map Control"); ` is called in my application. However, my application works on another machine running Windows XP SP3 which makes me believe that the cause of the problem is Windows 7 (or something I forgot about Windows 7). I am not familiar with ArcGIS and even I am not sure if I installed everything necessary for the test machines (I just installed just engine runtime and use a valid license). I wonder if using QAxCtl for controls is possible on Windows 7?


title: How do I merge two networks if they don't share nodes?

text: I have a street network and a subway network and I want to create a multi- modal network to analyze for Service Areas. The problem is that the two networks do not share nodes, and if a point is on one network, it cannot jump to the other. I've been able to manually add lines between subway stations and the nearest intersection, and also manually specify a time to enter/exit the station, but is there a more automated way to do this? The task is rather mundane and takes a lot of time, especially when performing analyses on many networks.


title: how can quantum be installed without administrative rights?

text: I have the impression from searching that this may have been possible with earlier versions but I don't know if it is possible with version 1.8. Thanks


title: SDA4PP not found? How to do Kriging in QGIS?

text: I need to make kriging analysis in QGIS (v1.7.4). I found out, that it should be a part of SDA4PP plugin, but I just couldn't find it in any repository (even third party). It looks, like it just disappeared from the web. Is it possible that I did something wrong? Thanks


title: mapping number of days a site is flooded from daily inundated areas maps

text: I have a set of daily inundated area polygons and would like combine them into one map that can give me a map showing the number of days a site remained inundated. I am unable to figure out any ideas ... any platform is ok with me I have ArcInfo, GRASS GIS and R and can use any of the three pretty well. Today is just not my day I guess. Thanking you in advance.


title: Appropriate distance metric for spatial clustering of geographic coordinates

text: I have a set of locations in geographic coordinates, and I would like to group the points using hierarchical clustering followed by tree-cutting at various "heights" in order to calculate group-wise means of variables at recorded at each location. Hierarchical clustering of the distance matrix of geographic coordinates, I presume, may be a misleading way to form groups because latitude and longitude are not equally spaced. I can then imagine two ways forward:   1. Using the great circle distance for the distance metric.   2. Converting the geographic coordinates to an equally-scaled projection and then finding the Euclidean distance. Apart from option two being more complicated to perform, are these approaches equivalent? And what exactly is the meaning of the tree cutting height in these cases?


title: Run time error for GP objects

text: I have a code with toolbars which I tried to run it in another PC but it failed. I get a run time error "Method Addtoolbox of object I GDispatch failed". The line that the code stacks is the last one shown below. Could you help me please? I work with ArcGIS 9.3 in VBA. Thanks           Dim GP As Object      Set GP = CreateObject("esriGeoprocessing.GpDispatch.1")          '''''  'Set the toolbox     GP.AddToolbox "C:/Program Files /ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx"


title: Overlaying user raster datasets on nasa world wind

text: Currently I am developing a desktop application with Nasa world wind JDK for 3d visualization of my data sets, but when i overlay my raster data sets those are vertically exaggerating with the default elevation layer.Can anyone help me how i can visualize my DEM data(SRTM) on nasa world wind in 3d with its own vertical exaggeration. Thanks in advance


title: Is there an open source international GeoCoding database if you only need city level resolution?

text: I'm writing an application that needs to plot a user's location on the map, and the user could be anywhere in the world, but I only need city/town level resolution, not street address. Is there a single open source database that contains this information?


I'm currently using QGIS 1.8.0 I would like to "merge" a raster data set on precipitation with a shapefile of administrative districts. Basically my the aim is to figure out the precipitation in different administrative areas. The data I have on precipitation is formatted in a lat/lon 0.5 grid. Each administrative district however intersects several of these 0.5x0.5 squares. How can I connect each administrative district to multiple corresponding 0.5x0.5 squares and then weigh the value of the squares that are only partially contained in the administrative district/ split up the squares that are divided between several districts? Thanks!


title: How to get hold of the Activated Event for the ArcMap window?

text: I am using ArcObjects 9.3 .NET. In my ArcMap dockable window, I have a single element host hosting a WPF user control (That's the only means of using WPF for ArcMap extensions as far as I know). While most things work fine, I am experiencing a known issue with WPF/Winforms interoperability that I can't quite solve. In particular, when I have focus on my dockable window (a textbox in my WPF control in the element host to be exact) and then alt-tab from ArcMap and go back into ArcMap again, the cursor continues to blink in the textbox within my dockable window, but no text appears when I type. In fact, the focus seemed as if it went to the main ArcMap window, and anything that I type appears to trigger the ArcMap menu shortcuts (F -> File, E-> Edit) etc. Further research led me to this thread, with a similar problem: http://social.msdn.microsoft.com/Forums/en- US/wpf/thread/1b50fec6-6596-4c0a-9191-32cd059f18f7/ The workaround proposed was to listen to the Activated Event of the parent form. Upon the Activated event, do a .Focus() on the WPF user control. However as the parent form in this case is probably the ArcMap window itself, I am having trouble accessing it from the DockableWindow. I can get the dockable window's parent control via .Parent or .ContainerControl, but these controls don't expose the .Activated event as they are technically not Form objects. Any ideas how I can listen to the activated event (aka when ArcMap gets refocused from an alt-tab event)? Thanks in advance!


title: Making polygon boundaries of different feature classes coincident?

text: I have two polygon feature classes:   1. Cadastral layer and    2. another polygon layer "zones".  In both feature classes the polygon boundaries are coincident inside their feature but not between the features. What I want to do is leave the cadastral layer boundaries fixed and snap the 'zones' layer boundaries to them as close as possible within a distance of 2 metres to remove all the small silvers. I have tried snapping zones to cadastral with a distance of 2 metres and that worked for some boundaries but not all slivers 2m of less. Why do some boundary lines snap to each other when inside the 2m distance and some do not? I have also tried integrate where you can rank one class over another and generally it works well for distances 1 metre or less but over 1m can even move the rank 1 class boundaries when snapping the boundaries together which is not suitable. What tool or combination of tools can I use to snap the boundaries of the different features together when they are within 2 metres of each other and boundaries outside of 2 metres stay as they are? The only strong rule I want to use is that the cadastral layer boundaries cannot move at all.


title: Problem loading spatialite extension with sqlite3 Python connector

text: I am trying to load the Spaltialite extension using SQLite3 python module. I am getting the following error when I run my script. > EOFError: [Errno 10054] An existing connection was forcibly closed by the > remote host               import sqlite3          conn = sqlite3.connect(r'C:\folder1\TestDB.sqlite')     conn.enable_load_extension(True)          conn.execute("select load_extension('libspatialite-4.dll')")               for row in conn.execute('SELECT AsText(Geometry) FROM alaska'):         print row      I am able to connect to this db and run a non-spatial query without any issues.


title: How to create a date field in QGIS shapefile

text: When I create a new shapefile in QGIS, I'm presented with the "New Vector Layer" window, where I specify the Shapefile type (Point/Line/Polygon) and am then able to create attribute fields. I'm only able to select the field type as Text/Integer/Decimal - how do I go about creating a date field type?


title: How much other information can be learned from this geocoordinate post?

text: At 5:30 a.m. CDT there was a post on reddit allegedly from Omegle with geocoordinates of a body of a girl who has been missing 5 months. Here is the imgur link http://i.imgur.com/ppqmbeD.png Here is the link to the reddit post. http://www.reddit.com/r/creepy/comments/2cod1q/an_anonymous_person_just_confessed_the_murder_of/ Can this community assist in identifying the person who made this post? Can you assist in verifying whether there is a body there? If you cannot, can you pass this inquiry along to other tech people who will? For Stack Exchange purposes, how much information can be learned from this small amount of data?


title: No result from pgRoutingLayer tsp(euclid) algorithm

text: I use QGIS plugin pgRoutingLayer for displaying the shortest route. Dijkstra algorithm works fine. I tried tsp (euclid) algorithm but it drew nothing (even no error message). I thought it may be something wrong with my data. So I uploaded pgRouting sample data into my database (http://docs.pgrouting.org/dev/doc/src/developer/sampledata.html) and tried tsp on them. But result was the same (no route and no error message). So I debugged the pgRoutingLayer plugin, caught the sql query. It looks like this and it results in 0 records               SELECT seq, id1 AS internal, id2 AS node, cost FROM pgr_tsp('                    SELECT DISTINCT id, x, y FROM                        (SELECT DISTINCT source AS id, x1::float8 AS x, y1::float8 AS y FROM edge_table                        UNION                        SELECT DISTINCT target AS id, x2::float8 AS x, y2::float8 AS y FROM edge_table)                        AS node WHERE node.id IN (7,8,11)',                    2, 12)      I may be missing some basics, I am a routing newbie. Can anyone give me a hand please? Thanks.


title: The difference when querying raster and vector

text: Hi everyone I have two queries have the same meaning but different results. Why is there a difference?. Queries are true? query 1: ![enter image description here](http://i.stack.imgur.com/ZZX6E.png) query 2: ![enter image description here](http://i.stack.imgur.com/Nffnt.png)


title: Buffering a PointGeometry in arcpy only returns a line?

text: I am trying to buffer a `PointGeometry` object in arcpy to a `Polygon` geometry object with ArcGis 10.2. The `PointGeometry` object is in lat/lon WGS84 coordinates and as far as I understand, I can only use `PointGeometry.buffer` by specifying the buffer in the same coordinates. To specify the buffer in metres, I would need to use: `arcpy.Buffer_analysis`. So I setup an empty geometry to do the buffer in like so:               spatialRef = u"GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],\                PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];\                -450359962737.05 -450359962737.05 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision"          pointGeom = arcpy.PointGeometry(arcpy.Point(-27.607361, 152.768311), spatialRef)     temPol = arcpy.Geometry()          pol = arcpy.Buffer_analysis(pointGeom, temPol, "50 meters")[0]      The resulting `Polygon` object however, is a long thin line, rather than the expected circular polygon. Any ideas what I am doing wrong?


title: Instructor led Web mapping development courses?

text: Hey guys I have been messing around with programing for a while now and have gotten proficient with Python and have some interest in HTML 5 and Javascript. I have recently received a request for an IRD (essentially a proposal to do research and development) from my employer for developing a trial web mapping application to be used as a tool for marketing web mapping services to clients and have functional use within the company. The API I have targeted is the ArcGIS Javascript API which relies heavily on the Dojo toolkit and have a targeted budget of $20,000 which includes time and materials. Since I am very interested in this subject area I have no problem spending most of the budget on training and developing in my free time; although I will really be the only person working on this. My questions for you guys are: What exactly will I have to know to be able to produce something that is functional and presentable? Also, if you can recomend any instructor led programs out there that I can put in my budget? I know there are a lot of free courses offered out there but since I'm kinda putting my ass on the line and have to have something to show at the end of this I would rather go with some training that is instructor led. EDIT: Basically my idea is to create a pretty standard web map with some of our larger clients data displayed in layers that can be turned on and off. One of the key concepts is the ability to print the current view to a standardized template. The demand for this functionality comes from our project managers who often times will need a hard copy map as an exhibit in a meeting but wait until the last minute to tell us in the GIS dept. With just 3 full time analysts this can be a major issue. EDIT2: I also plan to use this project as a jumping off point to get into much more advanced development, so the more I can learn the better.


title: Where do I save the Options file for ArcGIS 10.1 to release hanging licenses

text: I have created an options file for ArcGIS 10.1 so that I will be able to release a hanging ArcGIS license after a set amount of inactivity. I found that for ArcGIS 10.0 the correct directory to store the Options file on is C:\program files\ArcGIS\License10.0\bin. Unfortunately for ArcGIS 10.1 once I open up the ArcGIS folder in C:\Program Files\ I can not figure out where to save the Options file. Any help would be much appreciated.


title: C# geometry hex decoder

text: When i am storing in the database things like: Point(90 90). I get back 010100000000000000008056400000000000805640 which is a hex encoding of the geometric type. Does anyone know a library for c# that does the decoding of this hex back into a geometric type or into a string? I tried using Npgsql with no luck and gdal. EDIT: As mentinoned in the comments, i have considered doing the mapping directly from the database, but as its a very complicated db structure, i used NHibernate to map all other fields successfuly. The only problem is the geometric type. So i thought of creating a custom mapping in NHibernate. The get method works but comes back as a hex. I need to convert that to a string of some sort to work with. POstregSQLUserType(custom mapping type)               public object NullSafeGet(IDataReader rs, string[] names, object owner)             {                      object prop1 = NHibernateUtil.String.NullSafeGet(rs, names[0]);                      if (prop1 == null)                          return null;                                return prop1.ToString(); --->>> THIS HERE IS HEX. NEED TO CONVERT TO POINT.                  }             public void NullSafeSet(IDbCommand cmd, object value, int index)             {                      if (value == null)                          ((IDataParameter)cmd.Parameters[index]).Value = DBNull.Value;                      else                           ((IDataParameter)cmd.Parameters[index]).Value = "ST_GeomFromText(" + value.ToString() + ')';  --->>> This is converted in the database directly and converted there.                  }      Thanks M


title: How to Find the Average Raster Value of an area defined by a shapefile using R?

text: I have a set of raster images which represent a specific month over the years, and I want to make a timeline of the average values of an area using a shapefile. How do I extract the values from the rasters and import them in R in a manner that I can use them?


title: Copy single feature into another feature class

text: I do a lot of editing, in this case areas. Sometimes I add an area to one shapefile, only to realise a bit later that it should have been in another. Is there an easy way of copying only that polygon into the proper shapefile? Something along the lines of rightclick the row in attribute table, copy and then paste in the next attribute table. I do know that I can select the polygon and use Append, but this seem unnecessary complicated?


title: egeoxml vs geoxml

text: What is the difference between egeoxml and geoxml? What i know is both are client-side parsers and both allow me to edit the kml files on the run. I guess egeoxml is an initiative by google itself and would require the KML file to be at a publicly accesible location which is not the case in geoxml. Also the documentation available for egeoxml(http://econym.org.uk/gmap/egeoxml.htm) is more as compared to geoxml. But which one of them is better? I mean Performance-wise ? And an extension to this question is that if i have a website which is public, then there is nothing called locally accesible KML, because even geoxml requires the KML files to be present on the same server where the maps are.http://code.google.com/p/geoxml3/ (Last point in General Notes.)


title: Service URL for Web Feature Service (WFS) using GML Simple Features profile?

text: Does anyone know of an "almost always available" Service URL for a Web Feature Service (WFS) using GML Simple Features profile? I'm keen to test WFS client capabilities without having to configure a WFS server. # UPDATE The Answer to this Question provides advice on how to determine if a given WFS appears to be using SFGML.


title: Generate random points in selected areas for

text: I think my question is answered in different posts but I can't really make sense of them. **Background:** We are designing a survey of a sample population of farmers to identify the seeds they use. I have a vector layer indicating the land use (irrigated or rainfed). This layer is composed by 16000 polygons spread across the country. **What I need:** I want to randomly generate points within this polygons where to survey. More specifically, I want to generate 66 points in every province (provinces are indicated by another vector layer). In other words I want to generate 66 points under two conditions: 1) within province boundaries, 2) within the land use polygons. I already tried QGIS plugin "Random points" but it seems to crash the program... How should I proceed? Thank you.


title: How to create a link on a marked item?

text: I'm working with ArcGIS 10.x and QGIS 1.74 on a Windows XP SP3 System. In ArcGIS I get XML-Data which can be converted to a file-gdb using "IP NAS- Manager". Included in the XML-files there a owner informations which are stored in a kind of database, labelled filename.data.db, filename.index.db, filename.log.db and so on. Using ArcGis 10.x and the extensions "IP ALKIS- Karte" and "IP-ALKIS Buch" it is possible, to show owner information in a separeate window. To do this, you must start "IP ALKIS-Buch" before you start ArcGis 10.x and log on to a Tomcat-Server to verify, that you are a registered user. In ArcGIS 10.x now you select the desired extension, point and click to an item and will get the associated information of owner name, creation-date, length, width etc. of the selected item. All works fine, no problem. Now I'm trying to get the same information using QGIS, but I don't know the kind of database, where the owner information is stored. Fortunately, there is another way to get the information. The administration "Kreis" offers "Map Solution" to get the owner information. They also sent me a link to their web-service, where I can log in and search for information. This link can be modified by the item-number I want to search for. And now comes the problem: How can I realise, that   1. there is a way to use the link especially to connect to the desired website?   2. the selected item can be handed over to the link? The result should be similar to the result of clicking on an item with the "information" button selected and the associated table-information will be displayed. I hope, I could get you a glimpse of what I'm intending to do, otherwise feel free to ask. Greetings


title: Is there a map-level event to signal the start/end of drawing?

text: In OpenLayers, is there an event which fires when the map starts drawing **any layer** , and another for when the map has finished drawing **all layers**? I wish to display a "waiting" animation while any layer is still drawing. The help file for LoadingPanel mentions something which sounds promising, but the sample doesn't work and it mentions tiles (I'd like it to work on vector layers too). I'm looking for the equivalent of dojo.connect(map,"onUpdateEnd"... in the ArcGIS JS API: > Fires after layers that are updating their content have completed. This > event is often used in combination with onUpdateStart to display a "Map is > busy" or "Loading? " message as visual feedback to the end-user. (As a workaround, I could count the layers in the map, and monitor the status of each layer using `layer.events.register("loadend"...` until they have all finished. Are there other options?) Thanks


title: Thematic maps on Non Spatial Data in MapGuide

text: I need to create the thematic map on a polygon layer. But the data is not available in the layer; I am using Postgres database for that. It means I need to create the themes on layer using non spatial data. Please give me any hint to this, it is very urgent.


title: A strange OSM tiles in Mapproxy

text: Could someone please help me with my issue? I try to use Mapproxy with OpenStreetMap server of tiles in order to get output WMS for my qgis client. Something like this: ...               sources:        arcgis_tms:         type: tile         url:  http://c.tile.openstreetmap.org/%(tms_path)s.png     grids:       arcgis_online:         base: GLOBAL_MERCATOR         num_levels: 18         origin: nw      But I see something strange with my tiles(North and West are confused):http://img-fotki.yandex.ru/get/5804/5372191.0/0_76c21_4563963e_orig I change it to http://c.tile.openstreetmap.org/%s/%s/%s.png and get many mistakes like: > TypeError: not enough arguments for format string, finally I add > supported_format: image/png but the problem remained. Can someone explain how it works with tiles ( %(z)s/%(y)s/%(x)s.png ) ?


title: Looking for a colormap for SRTM like OCM Landscape

text: In Quantum GIS there's an OpenLayers plugin. I would very much like to use the same colormap/scheme/theme used in the OCM Landscape layer and apply it to SRTM image I have. I checked the the open cycle map website but didn't find any info regarding the color map they used so I wonder if anyone knows how to generate a similar colormap or QML file? TIA.


title: Grid and Graticule in QGIS Print Composer

text: Is there any other option to display DMS grid and graticule in QGIS Print Composer, since there is no plugins available yet?


title: "shutil.rmtree" doesn't get rid of my intermediate folders

text: I have a script that works fine, but now I'd like to remove all the intermediate folders I created therein in order to save memory and to be clean. I've been trying using the "shutil.rmtree ()" function, as I've always done before, but it seems not working this time. Why? Here is my script (the function is at the very end):               import arcpy, os, sys, shutil          try:         newpath_01= r'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_01'         if not os.path.exists (newpath_01):         os.makedirs (newpath_01)          arcpy.env.workspace = 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\MODIS_Clip\\Sample_test'         arcpy.env.overwriteOutput= True         from arcpy.sa import *              listMODIS = arcpy.ListRasters ()         arcpy.CheckOutExtension ('Spatial')              for MODIS in listMODIS:             outSnowmask = ExtractByAttributes(MODIS, "VALUE = 200")             outSnowmask.save (newpath_01+ '//'+ 'snow_'+ MODIS)               newpath_02= r'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_02'         if not os.path.exists (newpath_02):             os.makedirs (newpath_02)              arcpy.ClearEnvironment ('workspace')         arcpy.env.workspace = newpath_01              listMask = arcpy.ListRasters ()         aspect_reclass= 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\DEM\\per_MODIS\\aspect\\asp_rcls_cknp'              for mask in listMask:             arcpy.CheckOutExtension ('Spatial')             outExtractByMask = ExtractByMask (aspect_reclass, mask)             outExtractByMask.save (newpath_02+ '\\'+ mask)                   newpath_03= r'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_03'         if not os.path.exists (newpath_03):             os.makedirs (newpath_03)              arcpy.ClearEnvironment ('workspace')         arcpy.env.workspace = newpath_02         arcpy.env.overwriteOutput= True              listExtraction= arcpy.ListRasters ()         elev_range= 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\DEM\\per_MODIS\\recls_CKNPdem.lyr'              for extraction in listExtraction:             arcpy.CheckOutExtension ('Spatial')             name_MODIS= str (extraction [18:25])             ZonalHistogram (extraction, 'ASPECT', elev_range, newpath_03+ '\\'+ name_MODIS+ '.dbf')              arcpy.ClearEnvironment ('workspace')         arcpy.env.workspace = newpath_03              listTable = arcpy.ListTables ()         print listTable         out_dir= 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\test_CKNP_aspect\\'              for table in listTable:            arcpy.AddField_management (table, 'FILENAME', 'TEXT')            arcpy.CalculateField_management (table, 'FILENAME', '"'+ table+ '"')              arcpy.Merge_management (listTable, out_dir+ 'Snow_cknp_sample.dbf')              shutil.rmtree (newpath_01, newpath_02, newpath_03)          except:         print arcpy.GetMessages ()      Any help will be very appreciated.


title: Where can I find a good Sextante Modeller tutorial?

text: I'd like to enrich my knowledge of the Sextante Modeller. Do you have any good suggestion of a book/workshop etc. to dig into the modeller? Best Regards and thank you


title: open layers: how to acces an event attribute value from WMSGetFeatureInfo

text: My open layers script tries to get some features values from a WMS (which is a shapefile), to open a page in function of the attribute value returned. My problem is the following: \- the code returns a GET which contains the html code for displaying a table with each parameter names/values, and the value do really correspond to the clicked polygon. I can print the same values from event.text \- event.features is empty! Obviously that's not a proxy problem (since I get the values back). How can I fix the code to access the attributes values? Here is the code                  var extent = new OpenLayers.Bounds(-180,-90,180,90);         var options = {restrictedExtent:extent};              var map = new OpenLayers.Map("map-id", options);         var imagery = new OpenLayers.Layer.WMS(             "Global Imagery",             "http://maps.opengeo.org/geowebcache/service/wms",             {layers: "bluemarble", isBaseLayer:false}         );              var lmes=new OpenLayers.Layer.WMS(             "LMEs",             "http://aaa.bbb.ccc.dd/geoserver/ocean/wms",             {layers:'ocean:LME66', transparent:true, styles:'LMES_colors'},             {isBaseLayer:false, opacity:1, singleTile:true, visibility:true}         );         map.addLayers([lmes, imagery]);              map.zoomToMaxExtent();              function toggleRestrictedExtent() {           if(map.restrictedExtent == null) {             map.setOptions({restrictedExtent: extent});           } else {             map.setOptions({restrictedExtent: null});           }         }         var info = new OpenLayers.Control.WMSGetFeatureInfo({            url:'http://aaa.bbb.ccc.ddd/geoserver/ocean/wms',            title:'identify feature by clicking',            output:'features',                 eventListeners: {             getfeatureinfo: function(event) {               console.log(event);               console.log(event.features);               console.log(event.text);             }            }         });                   map.addControl(info);         info.activate();


title: FME - set symbology in TAB files in a directory

text: I want to set the symbology for all files in a directory. I found a way to do it with FME Workbench but it only works for one file. When I try to run it using batch deploy, it seems to be working for each file but I don't see any output. I need to do this for about 600 files in 12 directories. It's fine if it just replaces the existing files with the new symbology. The fmw is in https://drive.google.com/file/d/0BymNZZfQQ3a0VklIMnA3a215TG8/edit?usp=sharing ![enter image description here](http://i.stack.imgur.com/BZkKH.jpg) ![enter image description here](http://i.stack.imgur.com/3VmMT.jpg) ![enter image description here](http://i.stack.imgur.com/XMWfJ.jpg)                   Command-line to run this workspace:                   fme.exe tab_symbology_conversion.fmw               --SourceDataset_MITAB P:\2013\567_Convert_XYZ_to_TAB\Working\Speewah\*.tab               --DestDataset_MITAB P:\2013\567_Convert_XYZ_to_TAB\Working\symbolized\speewah               --SourceDataset_MITAB_2 P:\2013\567_Convert_XYZ_to_TAB\Working\Speewah\345000_8131000_1k_1m_DEM_ESRI.tab               Starting translation...     FME 2013 SP1 (20130325 - Build 13450 - WIN64)     FME_HOME is 'C:\Program Files\FME\'     FME Professional Edition (node locked-crc)     Serial Number: 0     Temporary License: 13 days left.     Machine host name is: ATGIS25     START - ProcessID: 9120, peak process memory usage: 38420 kB, current process memory usage: 38420 kB     FME Configuration: Command line arguments are `C:\Program Files\FME\fme.exe' `P:/2013/567_Convert_XYZ_to_TAB/Working\wb-xlate-1386293055262_10812' `LOG_STANDARDOUT' `YES' `LogCount' `d:\Temp\wb-xlate-counts-1386293055440_10812'     Shared directories for formats are : C:\Users\georgec\Documents\FME\Formats;C:\Program Files\FME\datasources     Shared directories for transformers are : C:\Users\georgec\Documents\FME\Transformers;C:\Program Files\FME\transformers     Shared directories for coordinate systems are : C:\Users\georgec\Documents\FME\CoordinateSystems     Shared directories for coordinate system exceptions are : C:\Users\georgec\Documents\FME\CoordinateSystemExceptions     Shared directories for coordinate system grid overrides are : C:\Users\georgec\Documents\FME\CoordinateSystemGridOverrides     Shared directories for CS-MAP transformation exceptions are : C:\Users\georgec\Documents\FME\CsmapTransformationExceptions     Shared directories for transformer categories are : C:\Users\georgec\Documents\FME\TransformerCategories     FME Configuration: Reader Keyword is `MULTI_READER'     FME Configuration: Writer Keyword is `MULTI_WRITER'     FME Configuration: Writer Group Definition Keyword is `MULTI_WRITER_DEF'     FME Configuration: Reader type is `MULTI_READER'     FME Configuration: Writer type is `MULTI_WRITER'     FME Configuration: No destination coordinate system set     FME Configuration: Current working directory is `P:\2013\567_Convert_XYZ_to_TAB\Working'     FME Configuration: Temporary directory is `d:\Temp'     FME Configuration: FME_HOME is `C:\Program Files\FME\'     FME Configuration: FME_BASE is 'no'     FME Configuration: FME_MF_DIR is 'P:\2013\567_Convert_XYZ_to_TAB\Working/'     FME Configuration: FME_MF_NAME is 'tab_symbology_conversion.fmw'     FME Configuration: FME_PRODUCT_NAME is '2013 SP1'     System Status: 768766MB of disk space available in the FME temporary directory (d:\Temp\)     System Status: 8388419MB of virtual memory available     Operating System: Microsoft Windows 7 64-bit Service Pack 1 (Build 7601)     FME Platform: WIN64     Locale: en_US     Code Page: 1252  (ANSI - Latin I)     FME Configuration: Process limit is 16340 MB (out of 16340 MB physical memory and 8388607 MB address space)     FME Configuration: Start freeing memory when process usage exceeds 9804 MB of memory or 8388007 MB of address space     FME Configuration: Stop freeing memory when process usage is below 7353 MB of memory and 6291005 MB of address space     FME Configuration: Autodetermining optimal maximum number of objects in memory     MULTI_READER(MULTI_READER): Will fail with first member reader failure     MULTI_READER(MULTI_READER): Adding MITAB Reader with keyword MITAB_1     MULTI_READER(MULTI_READER): Adding MITAB Reader with keyword MITAB_2     Using Multi Reader $Revision: 93943 $ ( $Date: 2012-10-30 14:39:34 -0700 (Tue, 30 Oct 2012) $ ) with keyword `MULTI_READER' to read multiple datasets     FME Configuration: Reading 64 dataset(s) with MITAB reader identified by MITAB_1     MULTI_READER(MITAB_1): Will continue past member reader failures     MULTI_READER(MITAB_1): Adding MITAB Reader with keyword MITAB_1     Using Multi Reader $Revision: 93943 $ ( $Date: 2012-10-30 14:39:34 -0700 (Tue, 30 Oct 2012) $ ) with keyword `MITAB_1' to read multiple datasets     Trying to find a DYNAMIC plugin for reader named `MITAB'     Loaded module 'MITAB' from file 'C:\Program Files\FME\plugins/mitab_fme.dll'     FME API version of module 'MITAB' matches current internal version (3.7 20121025)     Using MapInfo (MITAB) Native Reader $Revision: 94856 $ ( $Date: 2012-11-14 11:13:57 -0800 (Wed, 14 Nov 2012) $ ) to open dataset `P:\2013\567_Convert_XYZ_to_TAB\Working\Speewah'     Using MultiWriter $Revision: 59162 $ ( $Date: 2009-04-07 10:04:59 -0700 (Tue, 07 Apr 2009) $ ) with keyword `MULTI_WRITER' to output data (ID_ATTRIBUTE is `multi_writer_id')     Writer output will be ordered by value of multi_writer_id     Loaded module 'LogCount_func' from file 'C:\Program Files\FME\plugins/LogCount_func.dll'     FME API version of module 'LogCount_func' matches current internal version (3.7 20121025)     Set encoding 'fme-system' from input feature     Opened native MapInfo file `P:\2013\567_Convert_XYZ_to_TAB\Working\Speewah\345000_8131000_1k_1m_DEM_ESRI.tab'     The OGC definition of the FME coordinate system 'MGA94-55' is 'Earth Projection 8, 116, "m", 147, 0, 0.9996, 500000, 10000000'     MITAB_1 Reader: Using Enhanced Geometry     FME Configuration: Source coordinate system for reader MITAB_1[MITAB] set to `MGA94-55' as read from input data     Coordinate System `MGA94-55' parameters: CS_NAME=`MGA94-55' DESC_NM=`Map Grid of Australia Zone 55 (GDA 94)' DT_NAME=`GDA94' GROUP=`AUSNZ' MAP_SCL=`1' MAX_LNG=`150.741246785833' MIN_LAT=`-90' MIN_LNG=`143.258753214167' PARM1=`147' PROJ=`TM' QUAD=`1' SCL_RED=`0.9996' SOURCE=`National Mapping COuncil of Australia, Special Publication 10' UNIT=`METER' X_OFF=`500000' Y_OFF=`10000000' ZERO_X=`0.001' ZERO_Y=`0.001'      Reading source feature # 2500     Reading source feature # 5000     Reading source feature # 7500     <snip>     Reading source feature # 995000     Reading source feature # 997500     Reading source feature # 1000000     Closing native MapInfo reader     Trying to find a DYNAMIC plugin for reader named `MITAB'     FME API version of module 'MITAB' matches current internal version (3.7 20121025)     Using MapInfo (MITAB) Native Reader $Revision: 94856 $ ( $Date: 2012-11-14 11:13:57 -0800 (Wed, 14 Nov 2012) $ ) to open dataset `P:\2013\567_Convert_XYZ_to_TAB\Working\Speewah'     Set encoding 'fme-system' from input feature     Opened native MapInfo file `P:\2013\567_Convert_XYZ_to_TAB\Working\Speewah\345000_8132000_1k_1m_DEM_ESRI.tab'     MITAB_1 Reader: Using Enhanced Geometry     FME Configuration: Source coordinate system for reader MITAB_1[MITAB] set to `MGA94-55' as read from input data     Coordinate System `MGA94-55' parameters: CS_NAME=`MGA94-55' DESC_NM=`Map Grid of Australia Zone 55 (GDA 94)' DT_NAME=`GDA94' GROUP=`AUSNZ' MAP_SCL=`1' MAX_LNG=`150.741246785833' MIN_LAT=`-90' MIN_LNG=`143.258753214167' PARM1=`147' PROJ=`TM' QUAD=`1' SCL_RED=`0.9996' SOURCE=`National Mapping COuncil of Australia, Special Publication 10' UNIT=`METER' X_OFF=`500000' Y_OFF=`10000000' ZERO_X=`0.001' ZERO_Y=`0.001'      Reading source feature # 1002500     Reading source feature # 1005000     Reading source feature # 1007500


title: How to correctly use LIKE in the code block of a Reclass expression in ModelBuilder Calculate Field

text: I have a Reclass expression that works when i am search for single instance of a data field. I also wish the Reclass to work with a SQL type LIKE selection with Wildcards. e.g. I need to reclassify any example where the field include the text 'Inland Water' but i do not know what text may precede or follow this text term. I have tried the following but it does not work. I get a syntax error.               Expression     Reclass(!DescGroup!)          Code Block     def Reclass (DescGroup):       if (DescGroup == 'Sea') :         return 'K'       elif (DescGroup LIKE '%Inland water%') :         return 'G'      Any thoughts, tips or help would be very useful. Thanks. ![error message](http://i.stack.imgur.com/zRQKd.png)


title: Is ERDAS 2011 or ArcGIS 10.1 best to use when mosaicking 100's of ECW image files?

text: I have to perform hundreds of imagery files and they are in ECW files. They are 1 ft resolution and are in quarter quad and came from BLM but they are product by NGA. I wonder which software is the best method or best practice to perform them in ERDAS 2011 or ArcGIS 10.1 I have both software on my machine. I ran both software to perform mosaic; however, I found that ArcGIS 10.1 runs better than ERDAS 2011. I'd like to hear your suggest about which one should I use. I have to perform about 12,000 individual quarter-quad ECW files that will take me a couple of weeks to complete them in one big file. The image will be cover from bit stretch from Eastern Arizona all the way to south of El Paso, Texas. On the second hand, would I run model builder run this way be a much faster than the standard mosaic ?


title: Field does not get updated with Geoprocessing Service

text: Can anybody tell why the field gets updated when I run this script in ArcMap but not when I run a geoprocessing service? Thanks               #Import modules     import os, sys, arcpy, traceback, arcgisscripting          gp = arcgisscripting.create(10.1)          Stewardship = gp.GetParameterAsText(0)               Officevar = "\"College Station\""                    # This is to send the data input to the database     with arcpy.da.UpdateCursor(Stewardship, ("City")) as rows:             # row comes back as a tuple in the order specified here, so Office is row[0], Forester is row[1]             for row in rows:                 row[0] = Officevar                      rows.updateRow(row)


title: Leaflet Custom Control Fails

text: Why does the myMenu control not respond to a click?                   var map = L.map('map').setView([51.505, -0.09], 13);         L.tileLayer('https://{s}.tiles.mapbox.com/v3/{id}/{z}/{x}/{y}.png', {             maxZoom: 18, id: 'examples.map-i86knfo3' }).addTo(map);              var myMenu = L.control({position: 'topleft'});         myMenu.onAdd = function (map) {             this._div = L.DomUtil.create('div', 'myMenu');             this.update();             return this._div;             };              myMenu.update = function () {             this._div.innerHTML = "<h2>≡</h2>" ;             };              L.DomEvent.on(myMenu, "click", function (e) {             alert(50);             L.DomEvent.stop(e);             } )              myMenu.addTo(map);


title: Does ArcSDE require Oracle Spatial rather than Oracle Locator licensing for any functionality when working with SDO_GEOMETRY?

text: I'm looking for specific information on whether the ESRI suite requires Oracle Spatial for any of the functionality, or if Locator is always sufficient for working with SDO_GEOMETRY objects? Do you know of any ESRI documentation that discusses the minimum requirements in Oracle for working with SDO_GEOMETRY objects in the ESRI suite?


title: Calculate a polygon around a point - and conver to WKT

text: The use case I have is this: I want to take points (in my case a `data.frame` of lat/long coordinates) and calculate a polygon around each point of a given radius (in km, miles, whatever). To be clear, not a polygon around all points at once, but a polygon around each point separately. This seems like it would be straightforward, but how does one do this in R? The next thing I want to do is take those polygons and get well-known text (WKT) polygons because the web API I am querying against wants WKT. I think I can get a WKT polygon using `rgeos` like               library(rgeos)     g3 <- readWKT("POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))")     writeWKT(g3)          [1] "POLYGON ((1.0000000000000000 1.0000000000000000, 1.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 5.0000000000000000 1.0000000000000000, 1.0000000000000000 1.0000000000000000), (2.0000000000000000 2.0000000000000000, 3.0000000000000000 2.0000000000000000, 3.0000000000000000 3.0000000000000000, 2.0000000000000000 3.0000000000000000, 2.0000000000000000 2.0000000000000000))"


title: How do I use scipy.ndimage.gaussian_filter1d on an arbitrary axis?

text: I have a GDAL raster that looks like this: ![Unfiltered raster of priority values](http://i.stack.imgur.com/fsT8G.png) And I would really like to blur this raster along an _arbitrary_ transect. I see that scipy.ndimage has a gaussian filter that allows me to blur the raster along a given axis, as long as the axis is valid given the raster's numpy matrix. Doing so with the default axis of `-1` (and `sigma=2`) looks like this: ![Same raster with 1-dimensional gaussian filtering](http://i.stack.imgur.com/k3OBG.png) This raster is **pretty close** to what I'm looking for, I'd just really like to have this filtering not strictly happen along the horizontal axis, if that's possible. So my questions are thus:   * Is it possible to filter the matrix along an arbitrary axis, such as an axis that is not necessarily along the vertical or horizontal?   * Is it possible to do so just with Numpy/Scipy and GDAL?


title: How to find the maximum-area-rectangle inside a convex polygon?

text: In this post we are looking for **algorithms** / **ideas** on how to find the **_maximum-area-rectangle_** _inside a **convex polygon_**. In the following figure, numbers are the areas of the fitted rectangles. As shown a desired rectangle can vary in each dimension and can be in any angle. **Edit:** > We don't have any clear idea how to deal with the problem mentioned as so we > asked here. Nevertheless, we **guess** the maximum-area-rectangle may be one > of those that has one edge **_aligned_** on (not necessarily the same length > edge, of course) an edge of the polygon. ![enter image description here](http://i.stack.imgur.com/VMcT4.png)