title: Where should the data update date be entered into the FGDC metadata form (esri)?

text: if you've edited data recently (a gdb feature class), either spatial or tabular, and you want to capture the "update date" and possibly a note about what was updated, without overwriting the "publication date" or "edition" feilds in the metadata, where should it go? Perhaps i'm interpreting the meaning of the publication/edition dates incorrectly and that's what they're for...? any known/definitive answer to this?


title: I have a dozen years' GIS experience. Is it worth it to pursue a post-grad certificate in Accounting?

text: I've been working in GIS since 1998 - I've been at my present position for nearly eight years now. Right now my wife is suggesting I start a post-graduate certificate in accounting. It would take about two years, and then I could take the CPA exam, etc. Does anyone know of job out there where I could take advantage of both disciplines? Otherwise, I'd simply be starting out as an accountant around age 40 - I'm not sure I see the point of that. Anyone have any suggestions here? Much appreciated.


title: Set Data Driven Scale from Arcpy

text: Is there a way to set the data driven scale from arcpy? Looking at the documentation for the `DataDrivenPages` class, it doesn't seem to be accessible from outside of ArcMap. The script I run within the current map document exports the maps as jpgs based on the `Name` and `Page Number` fields I set, and at the `Data Driven Scale` I specified. When I put the same script into a model (with the necessary adjustments made) and run it from ArcCatalog, I get the same result, except the maps are now at `Best Fit` scale, which is the default for data driven pages.


title: ArcGIS Javascript Legend Is not Showing

text: I am trying to add a legend, and it gets to the console.log, but it fails at on the line after it. It does not show any errors saying it failed, but it wont continue past that line. `map[map_div_id]` holds the map (we have multiple maps on the page). Any Ideas?               //add the legend     dojo.connect(map[map_div_id],'onLayersAddResult',function(results){       var layerInfo = dojo.map(results, function(layer,index){         return {layer:layer.layer,title:index};       });            if(layerInfo.length > 0){         console.log(layerInfo);              var legendDijit = new esri.dijit.Legend({           map:map[map_div_id],           layerInfos:layerInfo         },"legendDiv");              legendDijit.startup();       }     });


title: openlayers - export map to canvas problem

text: I am trying to make a export control to export the map as image file. I found this example: http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/examples/exportMapCanvas.html so i added the control to the map and the function               function exportMap() {         var canvas = OpenLayers.Util.getElement("exportedImage");         exportMapControl.trigger(canvas);     }      I also added the ExportMap.js file:               http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/lib/OpenLayers/Control/ExportMap.js      but, i am getting this error, when i am trying to run exportMap():               Uncaught TypeError: Failed to execute 'drawImage' on 'CanvasRenderingContext2D': No function was found that matched the signature provided. ExportMap.js:128      What that mean? how can i fix it?


title: OpenLayers SelectFeature - Manually set feature as selected/unselected

text: In the OpenLayers SelectFeature, how can I manually set a feature as selected/unselected and have it reflect the appropriate color? For instance, when my application loads, it sends back the existing select state which I want to set on the polygons. Thanks!


title: Spherical Mercator to Cylinder mercator?

text: Hello I want to create an overlay for google maps. As far as I know google maps using Spherical Mercator, so the system where I read the maps use Mercator but is not specified which (spherical or cylinder) type it use. When i create the images, i get some information in common. That means I have an offset. I want to try to solve this using cylinder mercator coordinates, but i don't know how can i transform google coordinates to Cylinder Mercator coordinates. Thank you in advance.


title: How to force sql server to use the spatial index?

text: I'm beginner with SQL Server (2012). I created a spatial index, and I'm trying to force SQL Server to use that index to solve the query. Here's what I have: DECLARE @g geography SET @g = geography::STPointFromText('POINT(-79.25161 43.735965)', 4326) SELECT * FROM CUSTOMER WITH(INDEX(idx_Geography)) where HomeGeog.STDistance(@g) The spatial index (idex_Geography) is created for HomeGeog. I'm getting this error: Msg 102, Level 15, State 1, Line 5 Incorrect syntax near '('. I'm not sure where the error is... Thanks for you help!


title: QGIS & GeoTIFF clipping Error "-srcwin 5897 137 0 0 falls outside raster size of 11083x4158 or is otherwise illegal."

text: I'm trying to clip a GeoTIFF file and unfortunately I'm getting the following               Error "-srcwin 5897 137 0 0 falls outside raster size of 11083x4158 or is otherwise illegal."       From another post here I found out my projection was not set correctly. I did set the projection by clicking CRS status window at the bottom right and selecting Enable "on the fly" and WGS 84 and then apply. I'm still getting the same error. Any idea how to fix that ?


title: Is Geoserver's "styler" plugin dead?

text: The "Styler" plugin for Geoserver promised flexible styling of your data, using the GeoExt library: ![GeoExt Styler](http://i.stack.imgur.com/hLOQG.png) The user guide tells you to first install the REST plugin, and then extract the styler zip file into the /data/www Geoserver directory. However, the download page for Geoserver 2.2 does not include the REST plugin. Installing older versions (I'm using Geoserver 2.2) failed for me, with one error saying "This web site needs a different Google Maps API key" (from /geoserver/www/styler/ - I have a Google Maps API). REST seems to have disappeared recently: versions up to 2.0.3 contained the option to download the REST zip package, whereas after that there is no mention of REST, and the installation guide remains. I did find some evidence that REST may have died, however: > This feature is no longer updated as an independent GeoServer download, but > is available as part of the OpenGeo Suite. For more information on its > development see OpenGeo project page If Geoserver's "styler" plugin really is dead, then I should probably change the question to something more general like "how do I make OpenGeo's Styler work on Geoserver?" Just wanted to check with others who've had similar experience. **Context** I've set-up Geoserver on an AWS EC2 instance (testing map can be seen here). What I'm trying to do is to make those boring grey polygons colourful and interesting, based on their attributes. I've tried exporting the sld that discribes the style below into Geoserver, but, at present, this does not work. ![styled-polys](http://i.stack.imgur.com/euvA5.png) Apologies for such a verbose question!


title: upload table to ArcGIS Server

text: We are building a process which will allow users to submit a table of GPS data to update our central repository for all GPS data. Let's say that the target is simply a SQL data table. At first I thought about trying to accomplish this with an SQL INSERT wrapped in a Python script (exposed a GP tool). However, a colleague suggested that we somehow just get the table to the server and let integration services handle the Insert/append to the SQL table. My question is then, how to create a simple file retrieval service? I'm comfortable with the Python/GP service aspect of ArcGIS server. Basically it would be just the user submitting a table and placing it in a known location on the server. EDIT: I should add i'm not looking to upload a file from HTTP location; it must allow access to 'local' files on the user's system.


title: Speed of editing attributes in QGIS from a Python plugin

text: I'm trying to edit the value of an attribute for each feature in a layer using a QGIS Python plugin. I've found that doing this outside of editing mode is much slower than while editing (even including committing the edits). See code below (lines interchangeable at the same point in a loop). The difference in speed for my sample dataset is 2 seconds (edit mode) vs 72 seconds (not edit mode). Modifying an attribute in edit mode:               layer.changeAttributeValue(feature.id(), 17, QtCore.QVariant(value))      Modifying an attribute outside of edit mode:               layer.dataProvider().changeAttributeValues({ feature.id() : { 17 : QtCore.QVariant(value) } })      Is this an expected behavior? I don't need the user to be able to undo the changes, so I don't think I need to use edit mode. **Edit 1:** See full code below with both versions included (but commented out):               def run(self):         try:             # create spatial index of buffered layer             index = QgsSpatialIndex()             self.layer_buffered.select()             for feature in self.layer_buffered:                 index.insertFeature(feature)                  # enable editing             #was_editing = self.layer_target.isEditable()             #if was_editing is False:             #    self.layer_target.startEditing()                  # check intersections             self.layer_target.select()             self.feature_count = self.layer_target.featureCount()             for feature in self.layer_target:                 distance_min = None                 fids = index.intersects(feature.geometry().boundingBox())                 for fid in fids:                     # feature's bounding box and buffer bounding box intersect                     feature_buffered = QgsFeature()                     self.layer_buffered.featureAtId(fid, feature_buffered)                     if feature.geometry().intersects(feature_buffered.geometry()):                         # feature intersects buffer                         attrs = feature_buffered.attributeMap()                         distance = attrs[0].toPyObject()                         if distance_min is None or distance < distance_min:                             distance_min = distance                     if self.abort is True: break                 if self.abort is True: break                      # update feature's distance attribute                 self.layer_target.dataProvider().changeAttributeValues({feature.id(): {self.field_index: QtCore.QVariant(distance_min)}})                 #self.layer_target.changeAttributeValue(feature.id(), self.field_index, QtCore.QVariant(distance_min))                      self.calculate_progress()                  # disable editing             #if was_editing is False:             #    self.layer_target.commitChanges()              except:             import traceback             self.error.emit(traceback.format_exc())         self.progress.emit(100)         self.finished.emit(self.abort)      Both methods produce the same result, but writing via the data provider takes much longer. The function classifies proximity of building features to nearby fields (purple) using pre-created buffers (brown-ish). ![Proximity](http://i.stack.imgur.com/9tBRG.png)


title: Best Way to Handle Many Raster Files in QGIS?

text: I have a project with a large number of raster files, but it is becoming very slow. I wonder which is the best way to handle all this. I'm simply georeferencing these files creating .wld file (geoTiff is better?) and than I open them as raster. Am I too naive?


title: Calculating area

text: I have a polygon layer. I toggled to edit mode. I added a new column, while number (integer) and called it "area". I then use the field calculator - geometry - $area to populate the polygon area for my polgyons. Under Layer Properties - first tab - the unit is millimetre In metadata tab: Extents: In layer spatial reference system units : xMin,yMin 114.53,-38.7888 : xMax,yMax 153.639,-9.22883 Layer Spatial Reference System: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs Can anyone tell me what the area is calculated in? km2, acres, hectares or what?


title: Input data format and read multiple files into Grass 6.4.1

text: I am trying to read in multiple (hundreds, thousand) LiDAR ascii files into Grass. They contain elevation data. I think I have managed one, and I can make a map of elevations, using File, Import Raster Map, Multiple formats using GDAL. The format is as follows:               ncols         1000     nrows         1000     xllcorner     562000.000     yllcorner     7124000.000     cellsize      1     nodata_value  -9999.0     108.93 108.92 108.97 108.99 109 109.01...     ...     ...      1000 lines, 1000 columns, representing 1 sq km with values at every 1 m. Files are called               572000-7144000_dem_1m.asc     572000-7145000_dem_1m.asc     572000-7146000_dem_1m.asc      etc. My question is, how do Y read in multiple files? I have tried concatenating them, no luck. Thank you for your help.. regards, Szilard


title: Finding polygon parents

text: I have a PostGIS table full of polygon geometries. Some of these polygons are contained inside other polygons in the same table. I need to find which polygons have parents and self-relate them to their parent. How do I go about this? Can I use `ST_Contains()`?


title: Batch join non-spatial csv to shape files

text: I have 173 point shape files from soil testing and 173 non spatial csv files with lab ressults (one csv for each shape file) that need joining. The files are located in a library with folders named after each client, and the corresponding csv file is located in a sub folder of the client folder. I can join each csv with the corresponding shape in qgis, but the number of files makes this solution not feasible. What i would like to achieve is a script that loops through each client folder, pics out the csv file and joins the data to the shape file by the sampleID fields. My final goal is to loading all of this data in to a postgresql database table. Can any one point me to a solution for this problem?


title: Error in LineEdit when show parsing value as string

text: I m trying to make QGIS plugin to open xml, parse it and then show value i parse from xml to line edit, i use this code :               from PyQt4 import QtCore, QtGui     from ui_testparse import Ui_testparse     import xml.etree.ElementTree as ETree     # create the dialog for zoom to point               class testparseDialog(QtGui.QDialog):         def __init__(self):             QtGui.QDialog.__init__(self)             # Set up the user interface from Designer.             self.ui = Ui_testparse()             self.ui.setupUi(self)                  opendata = self.ui.btnCari             QtCore.QObject.connect(opendata, QtCore.SIGNAL('clicked()'),self.openxml)                        def openxml(self, event=None):                  #open dialog             openfile = QtGui.QFileDialog.getOpenFileName(self, 'Open File', '*.xml')                  self.ui.lineLokasi.setText(openfile)                  #call XML data             self.isiData(openfile)              def isiData(self, nmsatu):             #open teks with read mode             openteks = open(nmsatu, 'r').read()                  self.ui.textXml.setText(openteks)                  #Parse XML from Above             self.parsenow(openteks)              def parsenow(self, parse):             element = ETree.fromstring(parse)             xml_obj = ETree.ElementTree(element)             for title_obj in xml_obj.findall('.//{http://www.isotc211.org/2005/gmd}dateStamp/'                  '{http://www.isotc211.org/2005/gco}Date'):                     #print element                     list1 = [element]                     Str1 = ''.join(str(e) for e in list1)             self.ui.lineSkala.setText(Str1)      And the xml file i parse have a format like this:               <gmd:dateStamp>           <gco:Date>2013-12-12</gco:Date>        </gmd:dateStamp>        <gmd:metadataStandardName>           <gco:CharacterString>ISO 19115</gco:CharacterString>        </gmd:metadataStandardName>        <gmd:metadataStandardVersion>           <gco:CharacterString>ISO19115:2003</gco:CharacterString>        </gmd:metadataStandardVersion>      I want to show Date Value (2013-12-12) inside `<gco:Date>2013-12-12</gco:Date>` to lineEdit but instead of show date (2013-12-12) it show like some random code like this :               <Element '{http://www.isotc211.org/2005/gmd}MD_Metadata' at 0xd3e3c30>      Can someone help me


title: Calculating std dev for ruggedness using QGIS?

text: I have a DEM raster image and a vector layer with 200+ grids. Using QGIS, I got the terrain ruggedness index for the DEM layer. Now, I want the std dev (i.e. variation in ruggedness) for each of the 200+ grids. Using zonal statistics, I only get the mean. Please help.


title: How to visualize the Xdata attribute in AutoCAD?

text: I have converted my GIS data to CAD using data interoperability. The tags (attruibutes)would appear as Xdata in AutoCAD. How can I Show These attributes there? your assitance is appreciated in advance.


title: Best practice in writing automated tests for QGIS plugins

text: I'm looking for some best practice advice on writing automatic tests for QGIS plugins written in Python. I have written tests for Python scripts in the past using PyUnit (the `unittest` module), but have never done so for an application with a GUI. I've found a page describing how to use PyQt4.QTest to do unit tests on Qt widgets (http://www.voom.net/pyqt-qtest-example), but I'm struggling to see how I can use this with a widget that has been designed to run from within QGIS. The section on "Testing" in the PyQGIS documentation is notably absent. The best I have so far is:   * Keep the actual processing of data in isolated modules or functions, and write unit tests for those;   * Perform basic testing of the UI using QTest;   * Pray that it all holds together when using the plugin from within QGIS. Surely there is a better way?


title: How to skip over errors created by cursors running through empty tables?

text: I am trying to automate the creation of many different feature layers were their content is created by a cursor. Sometimes during the automation, a feature layer may be empty, this is okay, but as the cursor runs through the layer adding in the data, the `del row` line returns the error: `"NameError: name 'row' is not defined"` and stops the script. Question: If a layer I am creating happens to be empty, how can I skip over any errors that might pop up due to the cursor not being able to define itself in an empty table? Already Tried: I tried setting a 'count' variable to how many rows are in each table, and telling python to skip over the table if it has less than 1 row, but the "del row" line still caused a fuss, and I had trouble setting up the try statement. The code below is still returning the error. It is not printing empty cursor yet the `ACount` is indeed 0.               cursor = arcpy.da.UpdateCursor("A_Squares", ["LEGEND"])     ACount = arcpy.GetCount_management("A_Squares")     print ACount     if ACount == 0:         print "empty cursor"     else:         for row in cursor:             row[0] = TableA             cursor.updateRow(row)         del row         del cursor      Any ideas?


title: Find the closest line vertex from the linestring using PostGIS

text: How to find the closest line vertex from a linestring? `ST_ClosestPoint` will give a point which is not from the linestring.


title: Clipping faster (i.e. in bulk/batch) in ArcGIS for Desktop?

text: I have 100 different input vector layers and one clip feature vector layer. I have to clip them 100 times but its a large amount of work. Is there any Python script (ArcPy) in ArcGIS for Desktop to accelerate my work and avoid the one by one clipping?


title: ArcGIS for WPF - load Shapefiles and access attributes (dbf file contents)

text: I'm loading various shapefiles using ArcGIS for WPF and displaying them on a map (the mechanism for doing this is taken from the ArcGIS WPF runtime examples - see http://www.arcgis.com/home/item.html?id=953a530ada30477d95cc69e788268bc9). I have a shapefile that contains Point information, and contained in it's attributes are the text labels that apply to those points. I cannot see how to access that information. I know the DBF file is being used as if I delete that file the shapefile will not render as it has no symbol information. Despite this, I cannot see when debugging the provided example how to access this information. Obviously I could simply open the .dbf file manually but that seems awfully clunky. Does anyone have any help or pointers for me? Many thanks in advance. Lee.


title: How can I calculate the amount of memory required to run osm2pgrouting?

text: I have read that osm2pgrouting loads all the data into memory, which makes it faster. On the other side the size of the .osm file which it can parse is limited. How can I exactly know how much memory would I need to parse for example hamburg.osm into Postgres with osm2pgrouting? (hamburg.osm is 17 mb zipped and more than 100 mb unzipped). Is that possible to calculate the required memory size at least approximately? P.S. I have tried to parse a 300 Mb .osm file on a Linux machine with 512 Mb of RAM. It got to "Classes table created" and after that has not returned any "Finished" message. It also created 3 tables in the database - and populated the Nodes table with about 1.5 million nodes. (I have followed the the recent pgRouting workshop).


title: How can i use snapp point in when i use measurecontrol?

text: I have problem about measurecontrol . I can use measurecontrol and i can show length to use measure control . But my problem is that .i didnt use in snap because measurecontrol is a control not layer .`enter code here`                           measureControl = new OpenLayers.Control.Measure(             OpenLayers.Handler.Path,             {                 immediate: true,                 persist: true             }         );                    measureControl.getCustomLength = function(evt) {             var str = '';             str = 'Toplam: ' + (evt.units === 'km' ? evt.measure.toFixed(3) : evt.measure.toFixed(1)) + ' ' + evt.units + '   ';                // uzunluk =evt.units ==='km';                 uzunluk=evt.measure;             return str;         }             map.addControl(measureControl);             measureControl.events.on({             'measure': function(evt) {             $("#mesafe").html(measureControl.getCustomLength(evt));                 }                 });             measureControl.activate();      This is my measure control . And i want to snap this to point to my another vector.So when i look openlayers snapp example `enter code here`                                        var snap = new OpenLayers.Control.Snapping({                                 layer: measurecontrollayer,                                 targets: [{                                 layer:mylayer,                                 tolerance: 15                                 }]                                                 });                                 snap.activate();       snap is using only vector but my measurecontrol not vector .Plz show me some way i hope i can explain my question ...Also is there another way to use measure in layer linestring layer???? thx everybody


title: How to Add Network Analyst dll in My Web Mapping Application adf?

text: I'm a beginner in ArcGIS & I'm shooting blind really. I want to add the Network Analyst tool box to find the shortest path in my project. I tried this solution: I can't add Network Analyst dll to my web mapping application. I use ArcGIS 10 & Visual Studio 10 and I used web application adf to build my website. I tried the Network Analyst toolbox or networkanalyst.dll from developer kit folder but an error occurs: > The dll not found. Any help will be appreciated.


title: Reclassification of environmental values

text: I've started to use Qgis recently. I have a logical question. I'm trying to reclassify a raster. The values in the raster are depths, obtained by an interpolation. The original depth values (of the vector) are: 1 = <30° 2 = 30°-60° 3 = >60° after the interpolation with IDW I get a set of continuous values and I need to convert them into discrete ones. If I round up/down the continuous values I will obtain an over-estimation of the central values (2) Is there a better (more logical) method to reclassify this values?


title: Invoke PostGIS functions outside of the database? Need to force 3D in GeoDjango

text: I'm looking to ensure that all geometries saved to the database are three- dimensional. It's easy enough to do with a call to the database,               select ST_Force_3D( ST_GeomFromText( blah blah ) )      ...but that means for every object saved in the database, we make two database calls. Yuck! Is there a simple way of invoking PostGIS functions without making a call to the database? We're using PostGIS 1.5 on top of Postgres 9.1 with GeoDjango on top of Python 2.7. I've been working with GEOSGeometry objects, and they will tell you if your geometry is three-dimensional with the .hasz function, but I'm looking to FORCE everything to be three-dimensional without having to manually add a third dimension for every coordinate pair on every geometry we create. Thanks!


title: How to do spatial join which is also based on time using ArcGIS for Desktop?

text: I have a time series plot of house sales for 4 years, but in those 4 years some districts change. I was wondering if I need to split the data into 4 layers for each year in order to join them to their respective district at the time of the sale. For example, house 1 sells 2 times, once in 2007 and again in 2010. During that time the district changed from district 229 to 271. Is there any way to join them so that the right observation takes the right district identifying when doing a spatial join? Districts are polygons and addresses are points. The address points are geolocated based on address listed in the assessor information based on house number, road name, road direction, and zip code. I am using ArcMap 10.1.


title: How to create Near table with time series observations using ArcGIS for Desktop?

text: I have house sales that occur over 4 years. Some of those are called "distressed sales" and what I want to do is calculate all distressed sales within a year prior to the non-distressed sale. I'm not quite sure how to go about doing this. They are 2 different layers, each with time identifiers. I need to find all distressed sales within 500 feet of a non-distressed sale and happened within the last year of the sale. So a non-distress sale occured on June 7th 2010, I need to know all distressed sales within 500 feet since June 7th 2009. Does anyone know how I can go about doing this?


title: Raster Multiplication - Can't Multiply Sequence by Non-Int of Type "Float"

text: I am trying to multiply a list of GRID raster datasets using a different value for each raster. The value is in the raster name. Some examples are [wse_1979_5,wse_1980_0,wse_1980_5], where 1979_5, 1980_0, and 1980_5 will become 1979.5, 1980.0, and 1980.5. Every raster currently is classified where VALUE = 1 for all cells, and multiplying by the above values from the names should make the raster datasets equal to those values. So far I have:               arcpy.env.workspace = *path*     from arcpy.sa import *     rasters = arcpy.ListRasters()     for raster in rasters:         multValue = float(raster[4:-2]+"."+raster[-1])         outRaster = raster * multValue         outRaster.save(raster+"_r")     TypeError: can't multiply sequence by non-int of type 'float'      I am getting "TypeError: can't multiply sequence by non-int of type 'float'". I believe I need to somehow iterate the raster since it is apparently a string(?), but I'm not positive and can't quite figure this one out. It is simple to do in raster calculator, but I have a lot of datasets to multiply. Thanks!


title: How to split/explode ArcGIS annotation of multiple words into single word annotation?

text: We have transfer the labels of a certain layer to annotation, and it seems that the single feature of the annotation is polygon, for example there is a label for a road named "Pennsylvania Avenue Northwest" which is labeled along with the road, while I wonder if we can split this annotation feature to three features like "Pennsylvania","Avenue" and "Northwest" each with its geographic position and label angle. Is this possible? ## **update:** I am sorry but I missed that we have millions of features, so it would be better to implement it by program. :)


title: openlayers change marker icon, size and offset after drawn

text: I'm having much trouble changing my (feature)marker icons 'on the fly'. The problem is that I use much different images and also need to change the offset. So I iterate over every marker and make a new icon object with the latest data:                   if (iconImgUrl.indexOf('img/icons/icon_') != -1) {             var size = config.numberPlateTail;             var offset = new OpenLayers.Pixel(0, -size.h);             var iconImg = new OpenLayers.Icon('/img/icons/marker_' + dataitem.numberplate + '.png', size, offset);         } else {             var size = config.iconSize;             var offset = new OpenLayers.Pixel(-(size.w / 2), -size.h);             var iconImg = new OpenLayers.Icon('/' + iconImgUrl, size, offset);         }      Because I also want to change the offset, I can't use               icon.setUrl();     icon.setSize();      I tried to 'kill' the current icon object and change it for a new one and then try to force it to redraw the icon, but somehow the marker doesn't display after that. In the update loop I have:                       tbPosition.features[i].marker.icon.destroy();             tbPosition.features[i].marker.icon = iconImg;             tbPosition.features[i].marker.draw();             tbPosition.features[i].marker.icon.display(true);      and markers are created like:                var feature = new OpenLayers.Feature(this.markers, new OpenLayers.LonLat(lat, lng));      feature.data.icon = iconImg;      var marker = feature.createMarker();      this.markers.addMarker(marker);      this.features.push(feature);      So, how do I change (or update) the complete icon object?


title: Can ArcHydro Tools be used on a network in an SDE enterprise geodatabase?

text: Most of the documentation for ArcHydro Tools discusses their use on data in ESRI's personal and file geodatabases. Is it possible to use ArcHydro Tools to set network flow direction on a networked data set in an 'enterprise geodatabase' (SDE on Oracle)? Both the client and server are at version 9.3


title: Generate centroids for vector polygons based on raster and land-use data using QGIS

text: Although QGIS has the ability to generate centroids for polygons, I am not sure how to generate centroids based on land-use development using raster data. **Here's an example of what I'm trying to do:**   * I have a vector polygon layer, where each polygon represents a distinct zone for analysis.   * I have a color ortho raster layer that has the same geometric boundaries of my polygon layer.   * I am trying to develop a weighted centroid point layer based on the development density in the raster, yet still related somehow to the polygon it's originally in.   * Currently, I do this visually. It is time consuming to do 200+ zones and not as precise if a mathematical model was performing the analysis.   * I'm not sure if there is a plugin or a method for measuring the amount of roof space in a raster and using that and their density as a basis to develop centroids for the area **My Question:** How Can I generate centroids for vector polygons based on development density from a raster layer


title: Can QGIS display GeoRSS in real time?

text: QGIS has the ability to display GeoRSS by following these steps:   1. Click on "Add Vector Layer..." icon   2. Make sure the Data Source is "File" and enter the GeoRSS URL with "/viscurl/" prefix on the GeoRSS (e.g. /vsicurl/http://earthquake.usgs.gov/earthquakes/catalogs/eqs1day-M0.xml). The problem I'm having is that the canvas is not updated when a new earthquake happens. I have to quit qgis, restart qgis, then load the GeoRSS in order to see any updates. Anyone knows how to show GeoRSS in real time? I'd love to hear of any other FOSS solutions to this problem.


title: Previewing File Geodatabase in QGIS Browser - Is it Possible?

text: Can I Preview/look at file geodatabase attributes in the QGIS browser if I have installed QGIS using OSGEO4W installer and FileGDB API? Currently I can add/edit file geodatabase feature classes to QGIS desktop, but the browser does not allow me to see preview/attributes/metadata. All I get is the Param tab previewing all the indexes, tables, yadda, yadda. Thanks, Mike


title: ModelBuilder routine to remove spatially duplicate features

text: I'm attempting to build a ModelBuilder routine in ArcGIS 10.1 based off of an existing post on Stack Exchange. The response from 'Get Spatial' gave a great screen shot of a model that would run an input through an interative process to remove spatially duplicate features. I've attempted to build this model on my own but have run into two problems.   1. I have been unable to get the 'Output Feature Class' at the end of the routine to recognize the Input Features as a Feedback Variable. I've attempted making this connection through the Class Properties of the Output and the connect tool by picking the output and then the input.  What sort of restrictions on the Feedback Variable could be preventing me from making this connection?   1. I've attempted to run the model several times without the Feedback Variable connected and gotten an error that there is no input or fields that are required. I've entered this information in the parameters.  What reasons would cause this error to occur? Could it be related to the missing feedback Variable? Appreciate the help. I may need to reconsider the use of the Feedback Variable and the design of the model if I can't get around these issues. Thanks, Brian


title: Custom Feature Renderer in ArcGIS Server 10.1?

text: I've written a custom feature renderer. It works great in ArcMap (10.1). Now we want to deploy it to ArcGIS Server (10.1) so we can use it in map services. But it doesn't seem to work. Has anyone else successfully deployed a custom feature renderer to ArcGIS Server? More details: First of all, this is all running on Windows Server 2008 R2, 64 bit. And it's ArcGIS Server for Java (and Desktop) 10.1 with SP1 installed. To build the feature renderer, I followed the steps listed here: http://resources.arcgis.com/en/help/arcobjects- java/concepts/engine/index.html#/Explore/0001000006ws000000/ And like I said, it works great in ArcMap. The docs seem to indicate that you can deploy it to ArcGIS Server. I followed the steps at the following link, namely to build a ".jar" file and copy it to ArcGIS/Server/usr/lib/ext: http://resources.arcgis.com/en/help/arcobjects- java/concepts/engine/index.html#/Deploying_extensions/000100000671000000/ But when I try to deploy a service that uses the custom symbology, I get an error message (from the "Analyze" step in service deployment) that has these details:               Severity: High     Status: Unresolved     Code: 00006     Description: Layer's symbology is not supported.      I tried copying my jar to a few different places, just to see if it'd pick it up, but to no avail. (Using ProcessExplorer, I looked at the command line of the various ArcGIS Server processes and saw what entries it was using for the Java classpath and then put my jar in those spots.) Also, the link above says to deploy the jar at `<ArcGIS Server Home>/user/lib/ext`. This path doesn't exist on a stock ArcGIS Server 10.1 install. But `<ArcGIS Server Home>/usr/lib/ext` (i.e. without the "e" in "user") _does_ exist. Just to be sure, I tried both places (by manually creating the directory with the "e"). So I wonder: what additional steps are required to get ArcGIS Server to recognize our custom feature renderer? Or maybe their documentation is misleading and they don't really support custom feature renderers in ArcGIS Server.


title: trying to find areas within a certain driving distance along a road network WITHOUT a specific destination

text: I was hoping somebody might be able to help me out with this. I'm trying to use Arcmap 10.1 to find all the areas along a road network (a county roads layer in this case) that are a certain driving distance from a single origin point. In essence, this would be akin to a service area analysis, except there is only a single facility used rather than a whole list. Applicability: trying to find where a subject could move to within the county (using census block layers) if they want to live no more than, for instance, an hour's drive from work/school/etc. Any suggestions?


title: PyQgis overwrite shape file and adding layer to Qgis

text: I am creating a shapefile in python by:               driver = ogr.GetDriverByName('ESRI Shapefile')     datasource = driver.CreateDataSource('c:/temp/toke2.shp')     layer = datasource.CreateLayer('layerName',geom_type=ogr.wkbLineString)      If the file toke2.shp doesn't exist, it will be created. If it exist, I can't run this code and overwrite the old file. I need an overwrite function or a delete function. I run Qgis in admin mode. My other problem is when I want to add the shapefile to Qgis by:               vlayer = QgsVectorLayer("c:/temp/toke2.shp", "toke2", "ogr")     QgsMapLayerRegistry.instance().addMapLayer(vlayer)      If I add the file in Qgis by "Layer > add vector layer", it works fine. When I use the code in the python console in Qgis it also works fine. But if I use this code in my plugin the file is added to Qgis, but it is empty. So for some reason it doesn't add the datapoints in the file, it is just adding an empty shapefile. In my plugin I import this:                from osgeo import ogr      from osgeo import osr      import math      from PyQt4.QtCore import *      from PyQt4.QtGui import *      from qgis.core import *      import qgis.utils      import resources_rc      Does anyone have a solution to these problems?


title: Python script to convert Lat Long to ITM (Irish Transverse Mercator)

text: I am currently working on a python script that converts Lat Long (WGS84) to ITM. The formulas used are based on a paper published by OSGB (Appendix B): http://www.leica- geosystems.co.uk/downloads123/gb/gps/gps1200/other/OS%20Transformations%20and%20OSGM02%20user%20guide_en.pdf. I get good results in the Eastings but not in the Northings (from few mm to approx. 1.5m) and I can't understand why. Does anyone know the reasons? Thank for your help. Andrea               # this script convert Latitude and Longitude to ITM          # import modules     import math          # GRS80 Ellipsoid constants     a = 6378137.0     b = 6356752.3141     e2 = (a**2-b**2)/a**2          # ITM projection constants     F = 0.99982     Lat0 = 0.933751149817     Long0 = -0.13962634016     E0 = 600000.0     N0 = 750000.0          # input Lat degrees, minutes and seconds for latitude     dLat = float(raw_input('Enter latitude degrees: '))     mLat = float(raw_input('Enter latitude minutes: '))     sLat = float(raw_input('Enter latitude seconds: '))          # convert Lat to decimal degrees     ddLat = dLat + (mLat/60) + (sLat/3600)     print 'Latitude dd:', ddLat          # convert Lat to radians     rLat = (ddLat*math.pi)/180.0     Lat = rLat     print 'Latitude r:', Lat     print          # input Long degrees, minutes and seconds for longitude     dLong = float(raw_input('Enter longitude degrees: '))  #The degrees have to be entered with a posite sign     mLong = float(raw_input('Enter longitude minutes: '))     sLong = float(raw_input('Enter longitude seconds: '))          # convert Long to decimal degrees     ddLong = dLong + (mLong/60) + (sLong/3600)     print 'Longitude dd:', ddLong          # convert Long to radians     rLong = -(ddLong*math.pi)/180.0 # the negative sign account for west direction     Long = rLong     print 'Longitude r:', Long     print          # calculate constants for transformation     n = (a - b)/(a + b)     v = a*F*(1-e2*(math.sin(Lat)**2))**-0.5     p = a*F*(1-e2)*(1-e2*(math.sin(Lat)**2))**-1.5     n2 = v/p-1     print 'n:', n     print 'v:', v     print 'p:', p     print 'n2:', n2     print          M1 = (1+n+5.0/4.0*n**2+5.0/4.0*n**3)*(Lat-Lat0)     M2 = (3*n+3*n**2+21.0/8.0*n**3)*math.sin(Lat-Lat0)*math.cos(Lat+Lat0)     M3 = (15.0/8.0*n**2+15.0/8.0*n**3)*math.sin(2*(Lat-Lat0))*math.cos(2*(Lat+Lat0))     M4 = 35.0/24.0*n**3*math.sin(3*(Lat-Lat0))*math.cos(3*(Lat+Lat0))     M = b*F*(M1-M2+M3-M4)     print 'M1:', M1     print 'M2:', M2     print 'M3:', M3     print 'M4:', M4     print 'M:', M     print          I = M+N0     print 'I:', I     II = v/2*math.sin(Lat)*math.cos(Lat)     print 'II:', II     III = v/24*math.sin(Lat)*(math.cos(Lat))**3*(5-(math.tan(Lat)**2)+9*n2)     print 'III:', III     IIIA = v/720*math.sin(Lat)*(math.cos(Lat)**5)*(61-58*(math.tan(Lat)**2)+(math.tan(Lat)**4))     print 'IIIA:', IIIA     IV = v*math.cos(Lat)     print 'IV:', IV     V = v/6*(math.cos(Lat)**3)*(v/p-(math.tan(Lat)**2))     print 'V:', V     VI = v/120*(math.cos(Lat)**5)*(5-18*(math.tan(Lat)**2)+(math.tan(Lat)**4)+14*n2-58*(math.tan(Lat)**2)*n2)     print 'VI:', VI     print          # calculate Eastings and Northings     N = I+II*(Long-Long0)**2+III*(Long-Long0)**4+IIIA*(Long-Long0)**6     E = E0+IV*(Long-Long0)+V*(Long-Long0)**3+VI*(Long-Long0)**5     print 'Easting: ', E     print 'Northing:', N


title: Is it possible to extract Coordinates of Polygon Centroid in QGIS?

text: I have created polygon centroids using one of the "Geometry Tools" but I would like to extract the actual coordinates into Excel for further work. Is this possible?


title: QGIS Invalid Data Source in raster calculator

text: I just installed QGIS and am new to it. I want to use the raster calculator, but after assigning a name to the output layer QGIS tells me "Invalid Data Source" and "filename.tif is not a valid or recognized raster data source". What's wrong? Jan


title: Re-project raster image from Mercator to Equirectangular

text: I have a raster image of the world map in Mercator projection, square 1:1 format. I want transform the image into an Equirectangular projection accurately. I assumed from my research that GDAL would be the way to go, specifically gdalwarp. I used gdaltranslate to add gcp points to the image on the four corners, and then used gdalwarp to make the transformation. It produced an image with the correct size – 2:1 ratio – but the projection has not been accurately transformed. It looks exactly the same as if I had taken it into Photoshop and simply stretched it to that shape. The proportional relationship between Mercator and Equirectangular is obviously more complicated than a simple, global stretch like this. So my question is: is there a way to accurately transform between projections for raster images? **Edit / Additional Information** Here is the starting image: http://imgur.com/48wtsgn I then applied to GCP points using the following gdal_translate command:               gdal_translate -gcp 0.0 0.0 -180.0 90.0 -gcp 1000.0 0.0 180.0 90.0 -gcp 1000.0 1000.0 180.0 -90.0 -gcp 0.0 1000.0 -180.0 -90.0 Mercator.tif MercatorGCP.tif      Then I used the following gdalwarp command to make the intended reprojection from Mercator (ESPG:3857) to Equirectangular (ESPG:32662):               gdalwarp -s_srs EPSG:3857 -t_srs EPSG:32662 MercatorGCP.tif Equirectangular.tif      Here is the resulting image: http://imgur.com/5LGk16o Compare it to a map with proper Equirectangular projection, like this: http://imgur.com/3ncIW2R You can see that the proportions have not been properly adjusted, the map not actually reprojected at all, simply squashed uniformly. Am I doing something wrong, or is this simply not possible with gdal? If not, is there another way to do this? Thanks!


title: How can I add 3 axis orientation information to my photos?

text: I have a survey project where I need to record not only the location, but the direction (bearing and angle of elevation/depression) that the camera was pointing at the time the picture was taken. Accuracy is not critical. +/- 10 degrees would do. I would prefer not to have to buy a new camera. I don't need another GPS. I need to use a handheld GPS to get to each photo spot anyway, and integrating the GPS info with the photos is an easy ExifTool script. I don't want external cables. This is an all weather project in bushy terrain; cables catch on things, and require open port covers. Similarly I don't want 'big lumpy things' like FotoMapr. (Which also has GPS that I don't need.) The ideal solution would be a device the size of a pair of stacked nickels that would attach to the camera's hot shoe. The device would have 3 axis so that I could get compass direction and elevation angle. It would record these along with a time stamp at the time the picture was taken. Calibration to turn raw data to true north, true elevation would happen at data merge time. The device would have a mini-USB port for calibration and data retrieval. If it requires power, I'd like it to get it from the hot shoe, but battery is acceptable. The unit should have a cost under $100 I originally asked this over at the Photography exchange. They suggested I try here. Rebuttal to Answer # 1 Comments only allow 500 characters and no paragraphs. It's a bad solution:   1. It does not record dip, but just uses a 3 axis magnetometer to correct for tilted camera angles.   2. It has cables, and requires open ports on the camera.   3. It sticks out. It will catch on brush.   4. It's a GPS. I don't need to pay for another GPS. The 3 axis compass is a single chip under 1 cm square. It retails on a 1 inch printed circuit board for $10 quantity 1. It should be possible to make one that would in essence be a bump on the flash shoe. If it can't get a time sync from the camera, it would also need to have a clock chip. If it has a USB port, then it gets clock reset whenever attached to the computer. The computer could also use the amount of drift to apply a linear correction to the time stamp. So now our chip count is up to 4: Memory, USB controller, 3 axis magetometer, clock. * * *


title: Large loads geojson from mysql to leaflet?

text: I have a mysql database with +25K rows including geometry and other attributes. Using a mysql2geojson php script I then added it to a leaflet map. It is super slow now, what can I do to boost the performance? I'm using www.one.com as host meaning I dont have access to anything else than php/mysql. I can't get postgis/geoserver, etc. I also intend to design a search filter limiting which points to be shown on the map. For this purpose, is it best to do new calls to the database, or should I be working against a "local dump" for the searching? To summarize I wonder how to best work with big point data sets with mysql/php/leaflet including searching. Any tips are welcome!


title: Parameter values to create a vector grid with 1km x 1km cells?

text: I created a shapefile from WKT in QGIS using the values from this link. I want to create a vector grid with 1kmx1km cells but I can't figure out what values need to be set for the X and Y parameters. How do I calculate those?


title: How to close polylines that mark shape of cirque glaciers?

text: I have mapped a lot of cirque glaciers, using polylines to mark the shape of them. Now I need to Close these polylines, but I do not know how to do that. I have tried converting them to polygons (which is what would be the best thing to do for me) but since my polylines are not closed, it does not work. Since I have over 70 cirques mapped, I can not close them all by adding more vertexes to each one. (That would take some time!) I am hoping there is a way to just automatically draw a straight line between each polylines two endpoints... Can anyone help me? I use ArcMap 10.2.1, and am a beginner. ![Picture of the polylines that I need to close](http://i.stack.imgur.com/KX7Ot.png)


title: change spatial resolution

text: I have two types of image, one of them is a TRMM image with 0.25 degree cells and the other one is a MODIS image with 1 km spatial resolution. How can I change my TRMM image resolution to 1 km without GCP points?


title: Problem installing QGIS Valmeira on Windows 7 Professional 64-bit

text: I have been having trouble installing QGIS Valmeira on Windows 7 Professional with SP1. When I installed to the default directory C:\Program Files\QGIS Valmiera\, many of the predefined coordinate projections were not loaded. Several other projections were listed in French. At Andre's suggestion I then reinstalled to a folder without blanks such as C:\QGIS_Valmeira. This worked fine and brought in the previously missing projections. I did not have this problem with earlier versions of QGIS like Lisboa. What is going on here? Is this a bug in version 2? If so, it must be causing problems for many other Windows users of QGIS.


title: Query and Visualize GeoJSON data

text: I am trying to visualize specific GeoJSON data. I am successfully loading the full GeoJSON file (Multipolygon) using Leaflet JS. Is there any query to find specific polygon from that GeoJSON file? For an example I am using USA state GeoJSON data. I am using a dropdown list to call specific states such as California - CA. How can I call that specific GeoJSON data of where state = California?


title: Python - Find the largest extent in a list of extents (tuple format)

text: I am trying to find a way (simplest is the best) to retrieve the largest extent in a list of tuple extents in python. The list is looking like that:               extents = [(-180.0, -90.0, 180.0, 83.624), (-124.731, 24.956, -66.97, 49.372), (-122.42, -     37.818, 151.207, 52.516)]      Could it be a good idea to sum the absolute value of each coordinates and retrieve the highest? Or does a geographic library such OGR could be used? **Edit:** I have created Envelope objects with GDAL library:               envelopes = [Envelope(extent) for extent in extents]      Now I have to compare the area of each? Any suggestion? Can anybody give explanation on this? Thanks a lot.


title: ArcHydro Batch Subwatershed Delineation doesn contain catchment boundaries/falls short

text: ![enter image description here](http://i.stack.imgur.com/M0swe.jpg)When running the Batch Subwatershed Delineation in ArcHydro for 10.2/10.2.1, many of the output subwatershed areas do not completely enclose the catchments/adjointcatchments they are are representing. Instead the subwatershed polygon tends to begin about 1/2 raster cell size above the drainagepoint (which followed all necessary batch procedures before application) and in many cases it appears to run out of area from the nearby catchment boundaries and do not cross over the str grid or drainage point. Any thoughts on why this might be happening and how to resolve it to encompass the entire catchment/adjoint catchment areas is appreciated. -Justin


title: qgis - problem importing raster grid from GMRT

text: I am importing a .asc file into QGIS when doing so a window pops up saying that I need to define the layers coordinate system which is set on default WGS: 84 and the authority ID is EPSG: 4326. The file is downloaded through GeoMapAPP and is a apart of the Marine Science Geoscience Data system (MGDS) see this link http://www.marine-geo.org/portals/gmrt/. When I load the raster file it appears in my layers but there are no visuals at all. Any idea what the problem might be? Any assistance at all would be really helpful. Thanks, Anina


title: merging shapefiles . an error occured during the merge operation

text: i am trying to merge multiple polygons of a feature. i select them via select features by polygon and use the advanced digitalizing toolbar " merge selected features" option. it worked for some but mostly i get: error occurred during merging operation.


How do I check spatial reference on all parameters during Tool Validation? Just looking for recommendations. Here is what I have already. This must work with existing validation.               import arcpy     class ToolValidator(object):       """Class for validating a tool's parameter values and controlling       the behavior of the tool's dialog."""            def __init__(self):         """Setup arcpy and the list of tool parameters."""         self.params = arcpy.GetParameterInfo()         self.fcfield = (None, None)            def initializeParameters(self):         """Refine the properties of a tool's parameters.  This method is         called when the tool is opened."""         return            def updateParameters(self):         """Modify the values and properties of parameters before internal         validation is performed.  This method is called whenever a parmater         has been changed."""         if self.params[0].value and self.params[1].value:           fc, col = str(self.params[0].value), str(self.params[1].value)           if self.fcfield != (fc, col):             self.fcfield = (fc, col)                  # Get the unique values of the field 'col' in the feature class 'fc'             unique_values = [str(val) for val in                                 sorted(                                     set(                                         row.getValue(col)                                             for row in arcpy.SearchCursor(fc, None, None, col)                                     )                                 )                             ]                  # Assign the unique_values list to parameters 2 and 3             self.params[2].filter.list = unique_values             self.params[3].filter.list = unique_values              # Set the default values of parameters 2 and 3 to the first item in the list         if self.params[2].value not in self.params[2].filter.list:           self.params[2].value = self.params[2].filter.list[0]         if self.params[3].value not in self.params[3].filter.list:           if self.params[3].filter.list[0] != self.params[2].value:             self.params[3].value = self.params[3].filter.list[0]           else:             self.params[3].value = self.params[3].filter.list[.1]            def updateMessages(self):         """Modify the messages created by internal validation for each tool         parameter.  This method is called after internal validation."""         if self.params[2].value == self.params[3].value:             self.params[3].setErrorMessage("Values cannot be the same")


title: How to calculate cumulative sum conditionally across raster?

text: I don't know if this is possible without coding, but I need to:   * Create a raster file that, for each value in the input raster, sums all pixel values across the entire extent that are less than or equal to the input pixel. That value then needs to be written to the corresponding cell in the new raster.   * Scale the output between 1 and 100, but that is the easier part of   this. The only spatial analyst tool I can find that looks remotely useful is BlockStatistics (I was thinking I could conditionally run BlockSum) but that tool seems to require that the raster be split into multiple blocks. Is there a way to sort of hack BlockStatistics to only use one block/zone that corresponds to the extent of the entire raster? Or would I need to come up with a python script for this? Unfortunately, I don't think my coding skills are up for that. Thanks in advance for the help!


title: Convert ECF to Suface ECF

text: I'm using the 3D Earth WorldWind Model and I have a path that I want my object to travel. I set the the destination and interpolate to the end point from my current location. The problem is using ECF coords my object makes a straight line to the destination point, going through the earth. Whats a way to Input [X,Y,Z] and output [X,Y,Z] thats on the earths surface with a altitude of 0, through out the route.


title: Integrating ArcGIS for Desktop with Past Perfect for mapping public art objects?

text: I am working on a project where I need to map out several assets that are found around the city. The Arts Department utilizes a software called **Past Perfect** to catalogue and track any art assets. My part of the project involves mapping out where these items are located. Is anyone familiar with this program or know if it can be integrated with ArcGIS 10? I am trying to avoid redundant data, as I know that I could simply create a new database for the items being mapped; however, if there is a way to directly pull date from the database into ArcGIS Desktop, that would save a lot of time and streamline the entire process.


title: Extract single points of interest from OSM when described both as nodes and ways

text: I extracted both nodes and ways from OpenStreetMap in PostgreSQL with Osmosis. I'm looking for all restaurants in a given area. The issue is that some of them are coded in OSM as nodes and others as ways. How is it possible to have only one point when it is represented both as a node and a way? The intuitive answer is to intersect the ways with the nodes. If a node is inside a way, they will both count as one. Otherwise, if there is only a node or only a way without nodes inside, it's fine. The 'geometry' column in the 'ways' table does not contain information (it has something like a geometry, but all values are similar). Or maybe it is possible to use the 'nodes' table? The fundamental question is: How is it possible to have a proper geometry column in PostgreSQL with Osmosis?


title: PostgreSQL Trigger on inherited tables

text: In my PostgreSQL I have a parent table and a number of child tables (inherited from the parent table). How to create a trigger which will be bind to all the current and future child tables derived from the parent table? In a nutshell, what should be at "xxx" in the following code to refer to all the child tables?               create trigger trigger1     after insert or update of a, b, c, d on xxx     for each row     execute procedure trigger_procedure1();


title: ArcMap 10.1 convert layer to CAD - clipping only desired area to send out to DWG

text: Is there a way to covert a layer to DWG such that only the desired area of the layer is exported out to a DWG? My issue is that the layer (parcels) contains all parcels for an entire county, so it is a huge amount of information to convert to DWG. The conversion seems to just freeze up and will not complete. I actually only need a small area (a few parcels) that I'd like to have in DWG format. Is there a way to clip or outline a specific area for export to DWG? Thanks very much for any suggestions. Rob


title: Why does this python OGR code for creating a layer and then adding a field not work?

text: Spot the difference! So this works...               fp = "C:\\LiDAR\\MyLayer.shp"     layer = "MyLayer"     driver = ogr.GetDriverByName("ESRI Shapefile")     data_source = driver.CreateDataSource(fp)     srs = osr.SpatialReference()     srs.ImportFromEPSG(4326)     layer = data_source.CreateLayer(layer, srs, ogr.wkbPoint)          field_name = ogr.FieldDefn("Ben", ogr.OFTString)     field_name.SetWidth(24)     layer.CreateField(field_name)     layer = None      But this does not (Error 6: CreateField() not supported by this layer). Anyone know why and know what I need to do to fix? Thanks.               def createLayer(fp,layer):         driver = ogr.GetDriverByName("ESRI Shapefile")         data_source = driver.CreateDataSource(fp)         srs = osr.SpatialReference()         srs.ImportFromEPSG(4326)         return data_source.CreateLayer(layer, srs, ogr.wkbPoint)               layer = createLayer("C:\\LiDAR\\MyLayer.shp","MyLayer")          field_name = ogr.FieldDefn("Ben", ogr.OFTString)     field_name.SetWidth(24)     layer.CreateField(field_name)     layer = None


title: How to create the attribute table for a raster that was produced by the Slope Tool (ArcGIS10)?

text: I use the Slope Tool (Spatial Analyst Tool – Surface) of ArcGIS10 by ASTER GDEM (that was exported to UTM projection). The program produce new raster (with slope value), but the resulting raster doesn’t have the attribute table (it has it in inactive form). But it’s possible to create a map from this raster by slope value, it’s possible to represent each pixel of the secondary raster in its unique value, but impossible to open the attribute table. The same problem has the Aspect Tool (Spatial Analyst Tool – Surface), maybe some others. Also I can’t convert this raster to a polygon, because Raster to Polygon Tool doesn’t see it as a raster. How I can make the attribute table of the resulting raster active?


title: What is the most efficient method to iterate geoprocessing operations through an attribute table?

text: I want to calculate the clipped raster layer percentage for a shapefile consist of overlapping polygons, since ArcGIS doesn't support tabulate area for overlapping features, so I used a loop to process each of these overlapping feature one by one: 1: Create an independent shapefile for each polygon from the attribute table 2: Clip the raster layer using the single one polygon created from step 1. 3: Calculate the raster grid percentage by counting the numbers of raster grids in each categories. I was able to code this up, and the script run well but slowly for a big dataset. I run this loop for 2000+ big overlapping polygons and it took hours. I tried to optimize this computation by:               1: Use arcpy.da.SearchCursor instead of arcpy.SearchCursor          2: Using in-memory workspace to storage the temperate files.          3: Use Parallel Processing: arcpy.env.parallelProcessingFactor = "100%"      * * * These are all helpful but I still think can improve more. I am thinking about using python parallel processing tool such as map function. Could anyone suggest on how to apply the map function on this loop process? The following is the skeleton of my script:               featureClass = "E:/ArcGIS/access_to_recreation/new.gdb/Portland_ME_service_2hr"           rows = arcpy.da.SearchCursor(featureClass, ['FacilityID'])          # Loop through attribute table          for row in rows:         id = row[0]             # Extract each row as a independent polygon layer throughout the entire table         inputfc = "E:/ArcGIS/access_to_recreation/new.gdb/Portland_ME_service_2hr"            outputfc = "in_memory\\feature"         # edit here!         fieldname = "FacilityID"         fieldvalue = id         # construct the where clause by calling the helper function         where_clause = buildWhereClause(inputfc, fieldname, fieldvalue)         # execute Select analysis, create separate layers          arcpy.Select_analysis(inputfc, outputfc, where_clause)              # # Clip the raster data using each polygon          inRaster = "E:/ArcGIS/PAD/For_Recreational.gdb/PAD_Parks_continental_final"                 inMaskData = "in_memory\\feature"  # edit here!              # Execute ExtractByMask         outExtractByMask = ExtractByMask(inRaster, inMaskData)         outExtractByMask.save("in_memory\\extract_mask.img")          # build raster attribute table         try:             arcpy.BuildRasterAttributeTable_management("in_memory\\extract_mask.img", "Overwrite")                    # Count the 0, 1 pixels in the mask output           raster_rows = arcpy.da.SearchCursor("in_memory\\extract_mask.img", ['OID', 'Value', 'Count'])               try:             for ra_row in raster_rows:                 if ra_row[0] == 0:                     blank = ra_row[2]                 else:                     value = ra_row[2]             percentage = 100 * float(value)/(float(blank)+ float(value))         except:             percentage = 0              Percentage.append(percentage)         ID_list.append(id)               # Delete intermediate files         arcpy.Delete_management("in_memory")


title: How to implement MultiValue Parameter with Multiple Columns?

text: This might be a simple question, but I haven't been able to find a straightforward answer. I'm working on an ArcGIS (9.3) Python script. One of the inputs is a list of features. For each feature, I need to select a corresponding field. Setting the parameter to MultiValue gives a nice table in the UI for all the features being processed. How do I add another column to the MultiValue table to select a field for each feature? Cross-posted: http://forums.arcgis.com/threads/30907-MultiValue-Parameter- with-Multiple-Columns EDIT (I can't answer my own question for 8 hours?): ESRI's Answer from the ArcGIS forum: "Unfortunately, the solution for this was implemented just recently for 10.1 with ValueTable inputs for script tools. You will need to have two multivalue parameters with some custom validation to make this work for 9.3.1."


title: Calculate EVI From Landsat 8 in QGIS

text: I have calculated atmospherically corrected reflectance values for Landsat 8. From that, I am getting correct values of NDVI in QGIS, but for EVI I am not getting the correct values. For some images it goes beyond 1 and for others it shows black output in QGIS. I am using the EVI formula derived for MODIS.


title: Check if shapefiles exists in folder

text: In ArcGIS modelbuilder I want to check if a shapefiles exists, if not I want to create the shapefile in that folder. I guess the calculate value tool within ArcGIS modelbuilder can be used, but I need some python scripting to do the work. Any idea or python script examples? Thanks


title: Can I change properties of several layers simultaneously?

text: I'm working on bulk data. I have a few hundred GeoTIFF images imported into a QGIS project, and I'd like to modify their layer properties - things like the no-data-value, transparency of certain pixel values, colour map assignments... But I can't seem to find a way to apply these changes to several layers at once, short of editing the XML file with a global search-and-replace. Is there a way to do this using the GUI?


title: Use python to add layers to TOC

text: I've been battling with adding layers to the TOC for some time now. I have managed to single layers added but need to loop through a folder and add all layers. Shapefiles I cannot get to add at all. Only layer files. Below is the code of where I'm at:               for layer in shp_List:     mxd = arcpy.mapping.MapDocument("CURRENT")     dataFrame = arcpy.mapping.ListDataFrames(mxd, "*")[0]          outlayer = lyr     layerfile = os.path.join(base_Folder, lyr + ".lyr")          arcpy.MakeFeatureLayer_management(layer, outlayer)     arcpy.SaveToLayerFile_management(outlayer, layerfile, "ABSOLUTE")     addlayer = arcpy.mapping.Layer(layerfile)     arcpy.mapping.AddLayer(dataFrame, addlayer, "BOTTOM")                    #addLayer = arcpy.mapping.Layer(layer)     #arcpy.mapping.AddLayer(dataFrame, addLayer, "BOTTOM")     #arcpy.RefreshTOC()     #arcpy.RefreshActiveView()     del addlayer, mxd      The above code will add each layer to the TOC, then infuriatingly removes it and adds the next one. When the whole script finishes there is nothing left in the TOC. I've tried adding this script to a model and creating a derived output parameter which is then added as a model parameter with "add to display" checked. I've also checked the geoprocessing options to make sure the box to add to display is checked. I'm running Arc 10 (no service packs). Can anyone help?


title: How to get polyline features from a raster file using ArcGIS 10.1？

text: I have a watershed accumulation file look like this：![enter image description here](http://i.stack.imgur.com/OT3OL.png) There are several **white lines** by visual, and I tried the **raster to polyline** tool and got this： ![enter image description here](http://i.stack.imgur.com/63bvM.png) Are there better ways to extract polyline features from a raster file instead of creating an new shapefile and drawing lines by myself？


title: How to add point attributes to containing polygon?

text: I have three layers in a QGIS project: layer 1 is a boundary layer of polygons each with and id I will call "p_id" layer 2 is a point layer each with an id called "p_id" layer 3 is a .csv layer with count data that I am joining to both layer 1 and layer 2 on p_id (both the points and the polygons are the same type of thing) What I would like to do is for every polygon add the count data from each point contained within it to the polygons count data Is there a tool in QGIS that will allow me to do that?


title: Connecting to MSSQL 2008 R2 with QGIS 2.0.1

text: There have been questions related to this, but some of them seem out of date or no one has asked. I'm trying to connect to SQL Server 2008R2 using a MSSQL connection in QGIS. I'm not exactly sure if I'm connecting to the correct site. I'm able to connect through Autocad Civil 2014 so i'm using what I used to connect to that. Any ideas? The last posts I've seen are around 2011. ![enter image description here](http://i.stack.imgur.com/xVsM7.png) ![enter image description here](http://i.stack.imgur.com/iTMDy.png)


title: Create GeoTIFF for a specific scale

text: I am working on creating several GeoTIFFs in QGIS. My workflow was to: 1) Export the PDF map at the desired scale (1:64m) 2) Design the map and export the TIFs from Photoshop (I've used several resolutions to see if this was the issue) 3) Import those into QGIS and used the Georeferencer plugin and exported the geogreferenced TIFFs (Transformation Settings are Transformation Type: Linear, Resampling Method: Nearest Neighbor, Compression: None) When I pull the new TIFFs into QGIS and Zoom to the Best Scale the scale shows up around 1:58m rather than the desired 1:64m. I'm doing this for several scales so I need to figure out the best way to pull those in and get the scale to match up. Any ideas? Thanks!


title: postgis, extrapolate a line

text: I'm trying to extrapolate from a line segment to find a point on the line but a 3rd of the way 'back', i.e. trying to find point `new`, given points `A` and `B` below: ![enter image description here](http://i.stack.imgur.com/FTMi3.png) Given a line, I can interpolate it to get a position at any particular percentage along it:               =# select st_line_interpolate_point(        st_makeline('0101000020E6100000300DC347C49418C03EE8D9ACFAA44A40',                     '0101000020E6100000FB743C66A03218C0CDCCCCCCCC7C4A40'),         0.333);     0101000020E6100000ED45B41D537718C069C6A2E9EC984A40      I tried entering a negative number to find a point along the line in the opposite direction, but that fails as the interpolation argument has to be in the range [0, 1] I thought about first scaling the line, but that doesn't use the centre of the line as the origin, so is useless for my purposes.


title: Creating a map with selectable layers in Cartodb

text: I tried to make a map with selectable layers using the tutorial at http://developers.cartodb.com/tutorials/toggle_map_view.html, but was not successful. I really like this map by the LA Times http://graphics.latimes.com/responsivemap-pollution-burdens/ with selectable layers. I tried working off from this http://bl.ocks.org/andrewxhill/10506396, but when I load my map it appears with just the buttons and a black basemap (nothing displayed). My attempt is here https://dl.dropboxusercontent.com/u/253185126/index.html. Any tips or solutions would be helpful.


title: Polygon to raster conversion failure - workaround for very small polygons?

text: As one of many processes within a workflow/script, I need to convert polygons to rasters. I am using the feature_to_raster tool in ArcGIS 10.1. This is done within a python script, but for the purpose of trouble-shooting I'm also trying these things in ArcMap. For most features the process works fine, but there are a handful of features which result in empty rasters (i.e., only NoData cells). I do need to have at least one data cell created for each polygon; empty rasters are not acceptable. I determined that the problems are very small polygons which don't intersect any cell center of my reference (snap) raster, since the cell size is relatively large in comparison to the polygon. If I don't specify the snap raster, I've found that it does successfully create a raster with a single cell that does have a data value. But I DO need to specify a snap raster so that all my outputs line up, and when I do this, I end up with the empty rasters. I've tried using the alternate tool, polygon_to_raster, and specifying the cell assignment type to either MAXIMUM_AREA or MAXIMUM_COMBINED_AREA, but I still get the empty rasters. Does anyone know of a workaround? Thanks! EDIT 1: I got polygon-to-raster to work when I added an arbitrary priority field.


title: Python Script Producing Empty Cost Path Rasters

text: I'm trying to run a script that will create cost paths for a large number of randomly generated points then assign the cost to those points. But my script is beginning to create cost path rasters that don't contain a path. I'm pretty new at python scripting, and the code does need to be cleaned up a bit and there are a few print statements left over from where I was trying to find bugs, but my main concern is getting the cost paths created, then getting the total cost into the points. Every once in awhile the script will work completetly without error, but it takes 3-4 hours to run on my computer, so running it 5 times to get it to work isn't an option.               # Import system modules for coupling with ArcGIS, create geoprocessing link,     # acquire licenses, and load the necessary toolboxes     # copied from lab 4 script by Dr. Aldrich          import sys, string, os, arcgisscripting, arcpy     arcpy.CheckOutExtension("Spatial")               #set the work location for the script          arcpy.env.workspace = "C:\Users\Ryan\Desktop\GIS\python\project_test\project_test_gdb.gdb"          #script will overwrite existing data          arcpy.env.overwriteOutput = 1          #classes     ##add to this as necessary          class randompoint:         pass                              #variables     ##add to this as necessary     dem1 = "grdn36w089_13"     dest1 = "destination_one"     pointarea = "random_point_polygon"     point_names = []     costpath_names = []               #definitions     ##add to this as necessary                         #creates 100 random points, names each point sequencially, and adds that name to the list "point_names"     def create_points ():         pointcount = 0         while pointcount < 100:             pointcount = pointcount + 1             print "making point " + str(pointcount)             pointname = "random_point" + str(pointcount)             point_names.append(pointname)             arcpy.CreateRandomPoints_management ("", pointname, pointarea, "", 1, "", "")          #creates cost paths for points, names them sequencially, and adds that name to the list "costpath_names"             def create_costpaths ():         for point in point_names:             print "creating costpath for point " + point             cpname = "costpath_for_" + point             costpath = arcpy.sa.CostPath (point, "cost_for_first_destination", "back_link_for_first_destination", "", "")             costpath.save ("costpath_for_" + point)             costpath_names.append (cpname)                  #failed ideas             #costpath = None           #assigns the values from costpaths to their points using lists     #eventually replace the lists with classes that contain the name of the point and costpath     def assignvalues ():              checklists ()              counter = 0              print "assigning values"         while counter <= len(point_names) - 1:                  print counter                  pointname = point_names[counter]                  costpathname = costpath_names[counter]                  print pointname             print costpathname                  arcpy.AddField_management (pointname, "cost1", "DOUBLE", "", "", "", "", "", "", "")                  findcursor = arcpy.SearchCursor (costpathname, "", "", "PATHCOST", "")                  replacecursor = arcpy.UpdateCursor (pointname, "", "", "cost1", "")                       for item in findcursor:                 newvalue = item.getValue ("PATHCOST")                 if newvalue != 0:                     print newvalue                     for row in replacecursor:                         print newvalue                         print row.getValue ("cost1")                         row.setValue ("cost1", newvalue)                         print row.getValue ("cost1")                         replacecursor.updateRow(row)                   counter += 1          #checks the lists to make sure they are the same length     def checklists ():         print "checking lists"         if len(point_names) == len(costpath_names):             print "lists ok"         else:             raw_input ("Error: lists do not match, check script.  Press any key to continue")               #script          #generate slope raster          print "generating slope raster"     slope = arcpy.sa.Slope (dem1, "PERCENT_RISE", "")            slope.save ("slope_for_cost")                 #generate cost distance raster          print "generating cost distance raster"     cost1 = arcpy.sa.CostDistance (dest1, slope, "", "")          cost1.save ("cost_for_first_destination")           #generate back link raster          print "generating backlink raster"     backlink1 = arcpy.sa.CostBackLink (dest1, cost1, "", "")            backlink1.save ("back_link_for_first_destination")                create_points()          create_costpaths ()          assignvalues ()               raw_input('press enter to close')


title: ArcGIS equivalent of "Update Selection Set"

text: In MapInfo's MapBasic window I can perform a mass update of a field in a selection by entering UPDATE SELECTION SET TYPE = "APARTMENT". How do I do that in ArcGIS now that I have a selection?


title: Can Viewshed produce area results per point feature?

text: I want to quantify how much land individual point features can "see" using the Viewshed tool. I can generate a complete/universal raster surface, but not one that itemizes the number of cells visible per each individual point feature. I know that I can do this by selecting individual points and then run Viewshed, but given the number of points that I wish to calculate, this would be extremely time consuming. Does anyone know a way that I can calculate how much land can be seen (cell counts are fine) for individual point features without having to run the calculation individually on every point in a database. I was hoping for an output table with each point feature listed, followed by the number of visible cells (for the point feature). Perhaps a different tool?? Thanks in advance


title: Popup for polygon's points

text: i'm working with open-layers, my question is : how i can display popup for polygon's points ? the user can draw a polygon in a vector layer then i want to add popup if the point clicked ,that contain its coordinates and the user can edit the long/lat to redraw the polygon. is it possible? to return a point location i used this : var topleft = polygons_layer.features[0].geometry.getVertices()[0]; do i need to add another layer to contain these points as a feature , then the popup will be easily created? or its inefficient way ><" thanks.


title: How can I calculate the nearest distance to a landscape element taking the direction into account (i.e. slope direction, prevailing winds)?

text: I am trying to build a layer of distance to relevant landscape elements in a watershed. To be relevant, a landscape element must be downstream with respect to the runoff direction. Clearly, I want to assign to each point belonging to a regularly spaced points layer, the distance to the nearest landscape element, which is in the direction of the slope. To do so, I have :   * A polygon shapefile of « LANDSCAPE ELEMENTS » which are polygons of hedgerows, and forests (that are not distinguished).   * A Digital Elevation Model   * A shapefile of a regular point grid (50 m x 50 m, a total of 500 000 points ) Note : I'm using ArcGIS 10.1 with Spatial Analyst


title: How do you set an AddIn button's Help source and topicID values at runtime?

text: I would like a way to dynamically set the help source/topicID for a button. Example Snippet:                 <Button         id="Acme_FullExtent"         class="FullExtent"         caption="Full Extent"         tip="Zooms the map the the full extent."         category="Acme"         image="Images\FullExtent.png">         <Help source="http://www.Acme.com/help.htm" topicID="SectionB"/>       </Button>      More background: I need to link the button's F1 help to online Help that is published by someone else and which is a configuration nightmare for various levels of testing (e.g.Test, UAT vs Production). I like the idea of loosely coupled online AddIn help content, but the help content still needs to be tightly coupled with the appropriate environment (which is defined via configuration at runtime). I suppose the fallback plan is to compile the help into a chm (like pre-10.1 AddIns) and carry that onboard the AddIn.


title: QGIS 2.2 Print Composer does not display map with some datasets

text: I can't seem to get the Print Composer of QGIS 2.2 to work. When I click on "Add new map", the map square is created, but I can't seem to see my data in it, whatever the scale I use or even clicking on "set to map canvas extent". When I add a Legend, the contents are correct, all the layers I have rendered on the main QGIS window are listed there. Can anyone help me with this, please? Thanks!


title: What is the best way to manage geodata with open source software?

text: We are a small group of GIS users with different expertises embedded in an administrative infrastructure attempting to switch to a small and low cost open source infrastructure (open source GDI). Currently we’re using ArcGis 10.1 and 10.2 for producing thematic maps for schools. The aim is to use QGIS for future work. To become independent from location as well as from the administrative data infrastructure we are looking for solutions to store our data. The data should be managed central and accessible for anyone in our group from home. Until now we thought about an ftp server or a cloud solution. One key point is how to deal with simultaneous data access. For QGIS, are there any database server plugins that don’t require advanced database managing experiences? We are interested in any experiences how to manage this switchover and especially which software is practical.


title: Has Action behaviour changed between 1.7.4 and 1.8.0?

text: I have an Action that loads raster layers (.jp2) into the current QGIS window. This action has been working great in version 1.7.4. I have now upgraded to 1.8.0 for many reasons, but the same Action is not working. Here is the Action I have been using: `qgis.utils.iface.addRasterLayer("E:\\\Projects\\\CityTBay\\\Imagery\\\2008\\\%image")` The Action will load an ortho image from the path location specified. The image name is saved in the filed image. When I use this Action in 1.8.0 the following error occurs: `E:\Projects\CityTBay\Imagery\2008\%image is not a valid or recognized raster data source.` The Action is using the Python option. Thanks...


title: Is there a way to link ArcGIS to MS Access and use that DB?

text: I have an IT background but I am new to GIS. I recently recieved training using ARCMAP. I have about 5000 farmers that I need to put in a database and access the information using GIS. The vision is to have a map and when I click on a particular plot of land I need to see the related farming information (Owner, productivity etc) My question is - could I create a database using ARCGIS or is there a way to link ARCGIS to MSAccess and use that DB. I need someone to point me in the right direction. The resulting Map/Database should be Easy to use and update. Reports amay also neeed to be run on a monthly basis.


title: How can several vector layers have different controls, using Openlayers?

text: Using `openlayers` I have two vector layers. One layer is called `lineLayer` consisting of vector lines, and the other is called `pointLayer` consisting of points. When I hover on the lines I want to be able to highlight the line, and identify the vector object when I click on it. When I hover on the point I want to highlight the point and be able to drag it. My code works separately for lines and points, but not when I put it together. This is how I manage the line layer:               // handle interactions with lines     var highlightCtrl = new OpenLayers.Control.SelectFeature(lineLayer, {         hover: true,         highlightOnly: true,                     renderIntent: "temporary"     });          map.addControl(highlightCtrl);     highlightCtrl.activate();          var selectCtrl = new OpenLayers.Control.SelectFeature(lineLayer, {               clickout: true,          onSelect: function(feature){console.log(feature)}            });          map.addControl(selectCtrl);              selectCtrl.activate();      And this is how I manage the points:               // handle interactions with points      var highlightpoint = new OpenLayers.Control.SelectFeature(pointLayer, {         hover: true,         highlightOnly: true,         renderIntent: "temporary",     });          map.addControl(highlightpoint);     highlightpoint.activate();          var dragStart = new OpenLayers.Control.DragFeature(pointLayer, {         autoActivate: true     });           map.addControl(dragStart);     dragStart.activate();      From reading other examples I have noticed that the order seems to matter when a control is added. However, I've tried all the possible combinations of layer adding, and it does not seem to make a difference. Is there some way to do this that I am overlooking?


title: KMZ vector file over Google Maps

text: I have a KMZ vector file which works fine on my Google Earth application, but i want to overlay this on my phone (android) through Google maps. Is that possible??? thank you very much Mike


title: How to do a spatial search without select() using PyQGIS?

text: Using QGIS 1.9.0-master, I want to do a spatial search for features in a vector layer without selecting them, i.e., without using `QgsVectorLayer.select()`, `QgsVectorLayer.selectedFeatures()` et al. More specifically: I want to get the feature IDs inside a rectangle. That could be accomplished with something similar to `QgsSpatialIndex.intersects()` if I could gain access to the spatial index of the layer's dataprovider. Unfortunately, I didn't find a way to do it. I was only able to find how to create the index using `QgsVectorDataProvider::createSpatialIndex()` Does anybody know it this is possible?


title: How can we stop Manifold from writing TMP files?

text: We work with sizable Manifold Projects that sometimes go over a Gb of vector data. After working in Manifold for a while Gigs of TMP files are getting written to the hard drive. We use manifold 8 64 bit - so we would have assumed everything was done in memory as that's what 64 bit's great for! Is there anyway to stop this? We're running 12Gb x64 Windows 7 machines and there's got to be a performance hit in there... Plus our sysadmins aren't too fond of the Gigs of TMP files we leave lying around.


title: script for automated line drawing procedure for set of start and end points

text: A colleague conducted a field survey and recorded the start and end points of the transects with a GPS. These points have been imported to ArcMap 9.3 as 2 point layers. Both layers have a "transect number" field where each pair shares the same value, in order to match up start and end points. I can manually "connect the dots," but there are about 1000 pairs of points, and I am hoping there's a quicker way. If anyone knows of a simple script to run this procedure, I would greatly appreciate it! Thanks!


title: OpenLayers.Format.GML.parseFeature() with GML v3 data

text: I have a partial _GML_ object with the geometry type `polygon`. The GML version must be 3+ since the coordinates are wrapped in the `gml:posList` tag:               <gml:Polygon srsName="EPSG:25832" xmlns:gml="http://www.opengis.net/gml">            <gml:exterior>                <gml:LinearRing srsName="EPSG:25832">                    <gml:posList>                      N points where N is an even integer.                  </gml:posList>                </gml:LinearRing>            </gml:exterior>      </gml:Polygon>      I am using OpenLayers 2.10 and I am trying to convert the GML object into an OpenLayers feature using the _OpenLayers.Format.GML_ class. When executing the `parseFeature` method, which later calls the `parseGeometry.linestring` I encounter this comment:               /**     * Two coordinate variations to consider:     * 1) <gml:posList dimension="d">x0 y0 z0 x1 y1 z1</gml:posList>     * 2) <gml:coordinates>x0, y0, z0 x1, y1, z1</gml:coordinates>     */      And this line of code which assumes the `dimension` attribute to be present:               var dim = parseInt(nodeList[0].getAttribute("dimension"));       In _OpenLayers.Format.GML.v3_ the same line of code looks like this:               var dim = parseInt(node.getAttribute("dimension")) || 2;      The problematic difference for me is the `|| 2` part which assumes that a missing dimension means 2 dimensional - which is what I want. Then why not use the _v3_ edition, but the _OpenLayers.Format.GML.v3_ does not implement the `parseFeature` method. I am recieving _GML_ features from various services and I do not control them. **What is the best approach to handle this problem?**


title: Identity tool does not work as expected

text: I have two shapefiles. `Shapefile A` has `five attributes` and `shapefile B` has `20 attributes`. I use the **Identity** tool to determine, what of `shapefile A` lies in `shapefile B`. This is how both attribute tables look like: Shapefiel A ![A](http://s21.postimg.org/57la8jfgz/Untitled.png) Shapefile B ![B](http://s29.postimg.org/kkyfiq8iv/Untitled.png) My results are rather confusing ![C](http://s24.postimg.org/5o1zsrd9x/Untitled.png) How comes that it splits `BO`, `FE`, ... up, when that is contained in `Reg1`? Why is it contained two times in there? Working with ArcGis Desktop 10.2.2 here. Here is the snipped of the python code and the screenshot as PolyGeo asked for!               # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script     # The following inputs are layers or table views: "clip4_diss_clipCopy", "dissolve_clipCopy"arcpy.Identity_analysis(in_features="clip4_diss_clipCopy",identity_features="dissolve_clipCopy",out_feature_class="C:/Users//Documents/Scoring/up.shp",join_attributes="NO_FID",cluster_tolerance="#",relationship="NO_RELATIONSHIPS")      The underlaying green shapefile was used as the `Identity Feature`. The cyan lines are the different `Reg X` and the dark outlined features are the `Infeatures` for the tool! ![screen](http://s16.postimg.org/3s6t07prp/Untitled.png)


title: How to check whether an optional FeatureSet parameter was set by the user in a geoprocessing script tool?

text: I have created a custom geoprocessing script tool which contains an optional `FeatureSet` parameter. When unpacking the parameters in the script, the `FeatureSet` is defined as a `RecordSet` object, even though no input was received. It seems to set a default value if it was not entered. I need to check if this parameter was specified by the user or not. For example, when I run this through ArcMap without a `FeatureSet` parameter the following occurs:               exportAOI = arcpy.GetParameter(3)     arcpy.AddMessage('AOI is: ' + str(exportAOI))     # AOI is: <geoprocessing record set object object at 0x14133578>      I was expecting `exportAOI` to be `None`. How can I determine whether the user actually set this type of parameter?


title: plug logistic regression equation into Raster Calculator when one predictor is nominal

text: I am conducting a study that aims at building a Logistic Regression model for land productivity quality (binary dependent variable: good, no good) using different predictors such as slope, elevation, aspect, distance from coast- line, soil types, etc. Once I have got my Logistic Regression equation, to get a raster depicting the land productivity quality probability, I have plug into ArcGIS's (10.1)'s Raster Calculator the logistic regression equation. So far so good. On that respect, I have found a discussion (http://forums.esri.com/Thread.asp?c=93&f=993&t=270065) in which an example has been provided of the formula to put into the Raster Calculator. For example, `1.0 / (1.0 + Exp( 1.0 - 2.5 * [slope] ))` where the third 1.0 is the constant, and -2.5 is the coefficient for the slope, and the [slope] is the slope raster. Now, one of my predictors is nominal, i.e. soil types. It has been derived from a feature class where polygons where used to chart the different soil types. So, my logistic regression equation would run as follows: `1.0 / (1.0 + Exp( 1.0 - 2.5 * [slope] + 1.4 * [sandy] ))` where 'sandy' is one of the categories making up the 'soil types' raster. The problem for me is: how can I plug the 'soil types' raster into the logistic regression formula to be used in the Raster Calculator? I mean, how can ArcGIS multiply a nominal value (i.e., sandy) by 1.4.


title: How to remove the 'Report a map error' message in Google Satellite imagery?

text: Using Google Satellite as a basemap, I see the message 'Report a map error' in the bottom right corner of the screen (when the map is over land). Here is an example of this message on a similar map. I'd prefer not to have this message. Does anyone know whether   1. I can (legally) remove this message, and    2. How (using OpenLayers) I can do so?


title: Geodata Service in ArcGIS Java script api (3.1)

text: Can I use Geodata service in ArcGIS JavaScript API 3.1 ??? Background : I have simple map service and the all attribute data is stored in Sql Server database. I want to perform query operation on map service in JSAPI (and data should come from database.) > > I have only database access and I cannot create the Geo-database. I cant > attach the data to shape file (its too large) I searched on Google and asks few people so they told me that publish the data as "geodata service" but according to API documentation JSAPI (3.1) do not support Geodata service Also in JSAPI examples its not mentioned clearly. So what should I do ? Any help/suggestion will be highly appreciate.. Software's : ArcGIS Server 10.1 ,ArcGIS Java script API 3.1 ,ArcSDE (Not available)


title: merge raster images: error message procedure not found

text: I am using QGIS 1.8.0. Whenever I try to merge raster images I get the following error message:               Traceback (most recent call last):      File "C:\PROGRA~1\Quantum GIS Lisboa\bin\gdal_merge.py", line 36, in      import gdal      File "C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\gdal.py", line 2, in      from osgeo.gdal import deprecation_warn      File "C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\osgeo\__init__.py", line 21, in      _gdal = swig_import_helper()      File "C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\osgeo\__init__.py", line 17, in swig_import_helper      _mod = imp.load_module('_gdal', fp, pathname, description)      ImportError: DLL load failed: Die angegebene Prozedur wurde nicht gefunden.       Does anyone know what I sould do?


title: Qgis "Layers menu from project" plugin doesn't open files

text: Has anyone out there any experience with a plugin to Qgis called _Layers menu from project_? This plugin gives you fast access to layers that you specify in a Qgis-project file. It won't work on my computer. It appears to be set up correctly, but won't open any shapefiles when I click on them. Does this make any sense? :) D.


title: Transforming shp files to projected coordinate system

text: I have some KML files, polylines. I will convert them in grass to shp. How can I transform coordinates to some projected coordinate system? I was never doing it, so please help me. I need to transform them to Albanian 1987 GK Zone 4. I'm using Autocad Map 3D 2010 or ArcGIS 10.


title: Problem Downloading ESRI maps for Sharepoint

text: ![enter image description here](http://i.stack.imgur.com/6Lazq.png) When I try and download ESRI maps for sharepoint this error message comes up. Does anyone know how to resolve it?


title: How to find the Local Government Area name using lat/long?

text: I'd like my iPhone app to respond to the device owner with his/her Australian LGA (local government area) name using their lat/long. I've found this page below with a Shapefile: http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1259.0.30.001July%202010?OpenDocument However, on my Windows server I only have MS-SQL and my programmer's preferred language is c#. Can anyone possibly help with this? My programmer is having some issues with this task and I am pasting his email below as someone may spot what the issue could be: Bad News - SharpMap does not implements features we need (Within, Contains e.t.c.) for polygons our map contains. So, I could not write app logic we planned. Good News - SharpMap loads SHP file more correct then my previous custom application and therefore there is no errors as I see after loading. To keep moving forward I have found the next way: I fill each separate region using different color. In our case we have more than 500 regions. So, I shift B part in RGB color from 1 to 255, then G part, and then R part after all. So, I will have map like this one (in server memory, not displayed to user): http://screencast.com/t/xxMtI181z but more bigger size to avoid errors. Then I will re-calculate user coordinates to map pixel values and determine which pixel color is in that place. Using color value I will get region name from the list of regions (gotten from DBF file). So, looks like we still can solve this task. I will continue today. My next steps are: bigger image, re-calc coordinates, determine color.


title: Is there a maximun number of edges/vertices PGRouting can cope with?

text: We have imported OSM networks to route on. So our IDs sometimes grow abnormally big ;) It seems that PGRouting is failing with IDs like 633719916... We got this error message from postgresql: "terminate called after throwing an instance of 'std::bad_alloc'". Is there another way rather than renumbering our IDs? Is there a "fix" planned from PGRouting? Thanks for your help!