title: How to improve Python script processing speed (performance)?

text: I've been working on a script which essentially reads in a large amount of point and polygon data and generates frequency values for a featureclass of 2.5 acre hexes covering California. I have found it necessary to write my own version of the frequency tool, as ESRI's default tool has been crashing ArcMap on an irregular basis, and this script needs to be reliable and predictable. To do this I have been using python dictionaries to store values that are read in from the input datasets. The dictionaries are 2 dimensional (each key contains another dictionary within it) and get quite extensive, as there are 60000+ keys with 1-28 secondary keys. Writing the dictionaries isn't the problem, however, those get written reasonably quickly. My problem is when the dictionaries are done being constructed, running my makeshift frequency tool takes a long time. Here's the code and a short example of the dictionary entries, I'll explain a bit more below: There are two dictionaries I read from, here's an example of 'otherSourceDict'. hexDict is handled a little differently (it uses two letter codes (AA, AB, etc...) instead of taxon names like 'Amphibian'), but basically stores the same information:               {13456: {'Amphibian': 2, 'Bird': 5, 'Mammal': 10, 'Plant': 20}, 43156: {'Fish': 1, 'Plant': 4}}      The primary keys (13456 and 43156) represent hex feature unique IDs. The secondary keys (Amphibian, Mammal...) contain the number of taxon observations that are in that hex. My frequency code:               #Build a list of important correlated values for reference later     taxaList = [['RAR2_AMPH', 'AA', 'Amphibian', 'R2_NRM_A'],                 ['RAR2_BIRD', 'AB', 'Bird', 'R2_NRM_B'],                 ['RAR2_FISH', 'AF', 'Fish', 'R2_NRM_F'],                 ['RAR2_MAMM', 'AM', 'Mammal', 'R2_NRM_M'],                 ['RAR2_REPT', 'AR', 'Reptile', 'R2_NRM_R'],                 ['RAR2_PLNT', 'P', 'Plant', 'R2_NRM_P']]     #Add fields for taxa value by hex     sendmsg("     Adding fields...")     fieldList = gp.listfields(inRare)     for taxa in range(len(taxaList)):         if not taxaList[taxa][0] in fieldList:             gp.addfield_management(inRare, taxaList[taxa][0], "SHORT")     cur = gp.updatecursor(inRare)     row = cur.next()     sendmsg("     Writing dictionaries to hexs")     c=0     while row:         #Count for user visualization of process completion         c+=1         if c%10000 == 0:             sendmsg("          Starting record number " + str(c) + "...")         #For each taxa, calculate the frequency and populate the corresponding field.         for taxa in range(len(taxaList)):             val=0             if str(int(row.HEX25_ID)) in hexDict.keys():                 #Search the CNDDB dictionary and count how many of each taxa per hex                 for elmcode in hexDict[str(int(row.HEX25_ID))].keys():                     #If dealing with taxa other than plants                     if elmcode != "ECO_SECT" and taxaList[taxa][1] != "P":                         if elmcode[:2] == taxaList[taxa][1]:                             val+=1                     #If dealing with plants                     elif elmcode != "ECO_SECT" and taxaList[taxa][1] == "P":                         if elmcode[0] == "P" or elmcode[0] == "N":                             val+=1             if str(int(row.HEX25_ID)) in otherSourceDict.keys():                 #Search the other source dictionary and count how many of each taxa per hex                 for otherTaxa in otherSourceDict[str(int(row.HEX25_ID))].keys():                     if otherTaxa == taxaList[taxa][2]:                         val+=int(otherSourceDict[str(int(row.HEX25_ID))][otherTaxa])             #Populate the total number of each taxa per hex to the rare species hex features             row.setValue(taxaList[taxa][0], val)         cur.updaterow(row)         row = cur.next()          del row, cur      So right now, I have it so that when it goes to populate the hex features, it loops through each row 5 times, one for each field, then runs the frequency calculation for the specified taxon, and populates the field. There are roughly 60000 rows that it does this for. Does anyone have any suggestions or alternative methods for populating the hex features that could speed up the process? The script performs several operations, this by far takes the most time (about 25 minutes) and I would like to get it to be as quick as possible. Thanks for your help and suggestions!


title: How to use field calculator from QGIS Sextante plugin

text: I am trying to use the Sextante field calculator to calculate a field for a buffer radius using a categorical text field. The formula works fine when I run it using the standard QGIS field calculator via the attribute table: ![QGIS field calculator](http://i.stack.imgur.com/B1TwZ.jpg) When I use the same formula from the Sextante toolbox however (algorithms for vector layers), it produces only null values. ![enter image description here](http://i.stack.imgur.com/Ndkw2.jpg) ![enter image description here](http://i.stack.imgur.com/ls8II.jpg) This is the formula I'm using:               case when  "PointConfL"  =  '10-50 km' then 50000      when  "PointConfL"  =  '5-10 km' then 10000      when  "PointConfL"  =  '2-5 km' then 5000      when  "PointConfL" =  'centroid of farm' then 10000      else 50000     end      It seems to not matter whether the new field exists before or not: if it doesn't exist, it is created, but in both cases, all values are NULL. Am I using the 'algorithm' incorrectly, or is this perhaps a bug in the plugin? **EDIT** : I found the error log in `~/sextante`, but it doesn't shed much light on the problem:               ALGORITHM|Wed Oct 24 2012 17:08:08|Sextante.runalg("sextante:fieldcalculator","/home/rudi/personal/GIS/projects/transform_count/364-236.shp","buffer_sextante","case when  "PointConfL"  =  '10-50 km' then 50000 when  "PointConfL"  =  '5-10 km' then 10000 when  "PointConfL"  =  '2-5 km' then 5000 when  "PointConfL" =  'centroid of farm' then 10000 else 50000 end")     ERROR|Wed Oct 24 2012 17:08:08|Problem evaluation formula: Wrong field values or formula


title: How to accomplish conditional labeling based on map scale and attribute in QGIS 1.7

text: How can I label with qgis 1.7, based on attribute and scale of map? For example scale>1:1000 then label all those features whose area>100000 at font size 12 and those features whose area<99999 at font size 5.


title: Distance from rivers, taking into account terrain elevation

text: I am interested in building a raster of the distance from rivers in a territory under study. I have digitized the course of the rivers in ArcGis 10.1. I have used the 'Euclidean Distance' tool to create a raster of the distance from rivers. Now, since we are moving in a 3D world, not in a 2D one, I would like to 'weight' the calculated distances in order to take the elevation into account. So, I would like to know how can I create a raster depicting the distance of any point from rivers, also taking into account the elevation variable.


title: How to symbolize vector features using RGB values stored in three fields?

text: I need to symbolize a vector layer in Arcmap based on RGB (or HSB - can do either) values found in fields in the attribute table. I'm working with up to 10 different datasets, each with over 100 features. As a result, I can't manually adjust each feature's fill attribute, nor can I use a template .lyr file and then apply its symbology. Is there a method via Arcpy that I can use to perform this? Looking at the documentation it appears that symbology isn't exposed to Arcpy but I'm not familiar enough to know for sure. I could also possibly use ArcObjects, although I try to avoid that rabbit hole like the plague. Otherwise, are there any third-party scripts or tools that can accomplish this? I could also possibly go with QGIS, but only as a last resort.


title: How to validate whether shapefile exists using ArcPy and Python function?

text: I am very new to Python and am trying to complete my final project with not a lot of success. My code works without functions, the problem is I need it to work with functions. In the main(): function I want the user to enter an .shp file path. Then I have an isValid(inPath): function that I want to call the file path the user entered and validate if the file exists. I have done this in many assignments before with .csv and .gpx files, but I am not having success with .shp files. When I keep everything under the main function, it works no problem. Any help would be greatly appreciated! Here is my code currently:               **# Import Modules:**     import arcpy     from arcpy import env     from arcpy import mapping     import os.path          **# Set workspace**     env.workspace = #this is where I set my workspace     env.overwriteOutput = True          **# Checks if shapefile is valid**          def isValid(inPath):         if arcpy.Exists(inPath):             print 'Exists'         else:             print ' Does NOT exist'          def main():         shpFile = raw_input("Enter a Map Document path: ") # User inputs a shapefile         shpFile = shpFile.replace("\\","/") #replaces backslashes with forward slashes


title: [mapserver]: layer filtering with layerObj.setFilter()

text: I'm trying to filter a layer ( _using mapserver_ ) backed by a Shapefile using **layerObj.setFilter()** but this doesn't seem to have any effect on the returned tile(s). Sample code:               layerObj = mapserver.layerObj()     layerObj.setFilter("id = 1)      What am I doing wrong? Nothing blows up per say and mapserver produces tile(s) but they are the same set of tiles as if I had not applied the filter at all. Even applying the filter in the mapfile using **FILTER "id = 1"** also has no effect. Here is a more detailed ( _in code_ ) example of what I'm doing:               def index(self, *args, **kwargs):         mapObj = mapscript.mapObj(self.mapfile)              layers = [layer.strip().split(":", 1) for layer in kwargs["layers"].split(",")]              # create layers         for i, (layerType, layerName) in enumerate(layers):             layerObj = mapObj.getLayerByName("raster" if layerType == "models" else "vector").clone()             layerObj.name = "{0:s}:{1:s}".format(layerType, layerName)             layerObj.data = str(self.models.join(layerName)) if layerType == "models" else str(self.regions.join(layerName))                  if "filter" in kwargs:                 layerObj.setFilter(kwargs.pop("filter"))                  mapObj.insertLayer(layerObj)              owsReq = mapscript.OWSRequest()              # add query string parameters into map request object         for k, v in kwargs.items():             owsReq.setParameter(k.upper(), v)              # write mapserver results to buffer         mapscript.msIO_installStdoutToBuffer()         try:             mapObj.OWSDispatch(owsReq)         except:             return SException(self.response, str(sys.exc_info()[1]))              # adjust content-type         self.response.headers["Content-Type"] = mapscript.msIO_stripStdoutBufferContentType()         output = mapscript.msIO_getStdoutBufferBytes()         self.response.headers["Content-Length"] = str(len(output))              # output WMS image data         return output


title: Automate map creation for 50k unique routes

text: I have a table with 50k entries. There is a to field and from field with location in either lat/lng or street address. All these location are in the USA I would now like to automate the creation of unique jpg maps for each entry showing the route to get from point a to b. I have a ArcGIS 10 license along with spatial analyst (not sure if these will help). Any suggested workflows to automate the creation of these maps?


title: Querying and summarizing/reporting a multiple tables database

text: I have to create a database with information on Scottish buildings from 17 and 18 centuries. In order for my data to be searchable I need to have a series of tables connected with 1 to many relationships. Each building must be positioned on a map, which will have other layers like geo-referenced historic maps, and other features, like industries that supply materials needed for the building's conservation, etc... The database will be continuously used and updated. I was using Access, but I found that the connection with ArcGIS was not ideal: either I had to export from Access and Import in ArcGIS a single table or query, or I had to create an OLE DB connection, which isn't able to "read" the relations between the tables.   Therefore I, thought doing all the database directly in ArcGIS would be better. I created different tables connected with relationship classes, and one of them has X Y coordinates that give me the position of the buildings. The **problem** is I can't query more than one table at a time, nor create a report from more than one table at a time. what am I doing wrong? I'm new to databases and to both Access and ArcGIS, so if you have any suggestions on how I could make my database work better please tell me. Thank you.


title: How long does it take to load the complete Tiger 2013 dataset into Postgis?

text: I'm loading the complete TIGER 2013 dataset into a Postgis database and finding that the process is taking an extraordinarily long time. I was surprised by how long it took to load the TIGER 2012 dataset locally (over 60 hours), but the TIGER 2013 load is so long as to be almost unusable. Loading a single state's worth of data (whether it be Alaska, Alabama, or Arkansas) is taking 12+ hours from my local machine on a high speed network into an Amazon RDS Postgres 9.3 instance with Postgis 2.1. Note that this doesn't include downloading the files from ftp2.census.gov, since I pre- downloaded all the files. This is all time taking to process the SHP files and insert them into Postgres. I don't seem to be throttled by local resources - my CPU is nearly idle, and my load average is hovering around 2. My network has plenty of spare bandwidth. I can't imagine I'm throttled by network on the RDS instance, although it could be limited by CPU or memory (it's an m1.large RDS). There is literally nothing else running against the database though, so the load should be minimal. This situation is further complicated by the fact that the default TIGER scripts for Postgis don't lend themselves to parallelizability - they drop and create the tiger_staging schema, which makes it impossible to load states in parallel. I'm planning to modify the scripts to allow parallel load, but before I do that I'd like to know if there's a more fundamental problem. Has anyone done a full soup-to-nuts load of TIGER 2013 data? If so, how long did it take? Are there any configuration parameters I should be looking at that might be throttling this process? Any advice or guidance would be appreciated.


title: How to define leaders to a label in QGIS 2

text: In QGIS 1.8 and 1.9 I used a plugin called Easy Labels to develop leader lines that would be saved as a memory layer. I have not been able to replicate this with QGIS 2.0, I cant even find the Easy Label Plugin. Assuming it is no longer valid in QGIS 2. Question: How can I create a leader from a label moved outside a polygon to the centroid of the polygon?


title: Add labels for printed map (GeoTools)

text: I have some problems with labels for printing map with GeoTools. Next code (with some cuts, code on Scala language) can help me for printing map:                   val map = new MapContent()         val parser = new GeometryJSON()         val featureTypeBuilder = new SimpleFeatureTypeBuilder()         featureTypeBuilder.setName("Field")         featureTypeBuilder.setNamespaceURI("http://www.geotools.org/")         featureTypeBuilder.setCRS(DefaultGeographicCRS.WGS84)              featureTypeBuilder.add("geography", classOf[MultiPolygon])         featureTypeBuilder.add("text", classOf[String])              val featureBuilder = new SimpleFeatureBuilder(featureTypeBuilder.buildFeatureType())         val featureCollection = new DefaultFeatureCollection()          ...              val featureLayer = new FeatureLayer(featureCollection, createPolygonStyle())         map.addLayer(featureLayer)         val stream = getImageStream(map, width)         val bytes = new Array[Byte](stream.available())         stream.read(bytes)          ...      Next code set style of polygons, which I want show on printed map:               def createPolygonStyle = {         val styleBuilder = new StyleBuilder()         val styleFactory = CommonFactoryFinder.getStyleFactory         val filterFactory = CommonFactoryFinder.getFilterFactory              // Stroke         val stroke = styleFactory.createStroke(           filterFactory.literal(Color.RED),           filterFactory.literal(3),           filterFactory.literal(1)         )              // Fill         val fill = styleFactory.createFill(           filterFactory.literal(Color.YELLOW),           filterFactory.literal(0.2)         )              // Set anchor label         val anchorPoint = styleBuilder.createAnchorPoint(           filterFactory.literal(0.5), filterFactory.literal(0.5)         )         val pointPlacement = styleBuilder.createPointPlacement(           anchorPoint,           null,           styleBuilder.literalExpression(0)         )              // Set label style         val font = styleFactory.createFont(           filterFactory.literal("Arial"),           filterFactory.literal("normal"),           filterFactory.literal("normal"),           filterFactory.literal("20")         )         val textSymbolizer = styleFactory.createTextSymbolizer(           styleFactory.createFill(             filterFactory.literal(Color.BLACK)           ),           Array(font),           styleBuilder.createHalo(),           styleBuilder.attributeExpression("text"),           pointPlacement,           null         )              textSymbolizer.getOptions.put("maxDisplacement", "150")              val polygonSymbolizer = styleFactory.createPolygonSymbolizer(stroke, fill, null)              val style = styleBuilder.createStyle()         style.featureTypeStyles().add(           styleBuilder.createFeatureTypeStyle(             "Field",             Array[Symbolizer](polygonSymbolizer, textSymbolizer)           )         )              style       }      And then I have next result, but not all polygons have good labels: ![enter image description here](http://i.stack.imgur.com/1ZSTl.png) What can I change in my code for showing all labels for polygons with scaled font size and without outbounds text?


title: Android ESRI map getting Failed to generate sequences for graphic error

text: I am using following code to add graphic layer in my ESRI map.               public void footprintdata(String option)      {         graphicsLayer.removeAll();            if(option.equalsIgnoreCase("ALL"))       {        try         {             Graphic[] resultGraphics = new Graphic[footprintdata.size()];         // Envelope to focus on the map extent on the results         Envelope extent = new Envelope();              for (i = 0; i < footprintdata.size(); i++)          {          double lattitude = 0.0;          double longitude = 0.0;                 if(!footprintdata.get(i).getLatitude().equalsIgnoreCase(""))          {           lattitude = Double.parseDouble(footprintdata.get(i).getLatitude());          }          if(!footprintdata.get(i).getLongitude().equalsIgnoreCase(""))          {           longitude = Double.parseDouble(footprintdata.get(i).getLongitude());          }                    mLocation = new Point(XCoordinateFromLongitude(longitude), YCoordinateFromLatitude(lattitude));               mapPoint = (Point) GeometryEngine.project(mLocation,            SpatialReference.create(4326),mapView.getSpatialReference());                 SimpleMarkerSymbol resultSymbol = new SimpleMarkerSymbol(Color.GREEN, 20, SimpleMarkerSymbol.STYLE.DIAMOND);          // create graphic object for resulting location          //PictureMarkerSymbol imagesymbolicon = new PictureMarkerSymbol(imagemarker(i));               Graphic resultgraphics = new Graphic(mapPoint,resultSymbol);              resultGraphics[i] = resultgraphics;          extent.merge(mapPoint);         }              graphicsLayer.addGraphics(resultGraphics);         mapView.setExtent(extent, 100);         mapView.zoomToResolution(mapPoint, 10);         mapView.invalidate();                  } catch (Exception e) {         // TODO Auto-generated catch block         e.printStackTrace();        }       }       else       {             System.out.println("else of spinner selection");             try         {             int count = 0;         int temp = 1;         for (int i = 0; i < footprintdata.size(); i++) {               if (footprintdata.get(i).getLocationtype()            .equalsIgnoreCase(spinner.getSelectedItem().toString())) {                count++;          }         }              System.out.println("count :"+ count);         Graphic[] resultGraphics = new Graphic[count];         Envelope extent = new Envelope();               for (int i = 0; i < footprintdata.size(); i++) {               //PictureMarkerSymbol imagesymbol = null;          if (footprintdata.get(i).getLocationtype()            .equalsIgnoreCase(spinner.getSelectedItem().toString())) {                System.out.println("result success : "             + footprintdata.get(i).getLocationtype());                double lattitude = Double           .parseDouble(footprintdata.get(i)             .getLatitude());           double longitude = Double           .parseDouble(footprintdata.get(i)              .getLongitude());                mLocation = new Point(XCoordinateFromLongitude(longitude), YCoordinateFromLatitude(lattitude));                mapPoint = (Point) GeometryEngine.project(mLocation,             SpatialReference.create(4326),mapView.getSpatialReference());                SimpleMarkerSymbol resultSymbol = new SimpleMarkerSymbol(Color.GREEN, 20, SimpleMarkerSymbol.STYLE.DIAMOND);                //PictureMarkerSymbol imagesymbolicon = new PictureMarkerSymbol(imagemarker(i));           Graphic resultLocation = new Graphic(mapPoint,resultSymbol);                System.out.println("temp :"+ temp);           resultGraphics[temp] = resultLocation;           extent.merge(mapPoint);           temp++;               }         }              graphicsLayer.addGraphics(resultGraphics);         mapView.setExtent(extent, 100);         mapView.zoomToResolution(mapPoint, 10);         mapView.invalidate();                 } catch (Exception e) {         // TODO Auto-generated catch block         e.printStackTrace();        }       }      }      this code works fine at first time. but after it I am getting `MapCore(24187): Failed to generate sequences for graphic` in My Logcat. please help if some one have any idea about it.


title: Data Interproability Toolset not working ArcGIS 10.2

text: The Data Interoperability Tools .tbx appears in my system toolboxes, however when I double click quick export or import, a red X appears. I have gone into the administrator settings and the tool doesn't appear to be installed. How do I install it? I am using ArcGIS 10.2 for desktop trial. After I attempt to evaluate extensions, I receive a message which states all requested features have already been authorized on this machine.


title: How to convert an IGeometryCollection into a Polyline

text: I am using the CostPathAsPolyline method which returns the Polyline as an IGeometryCollection. I now need to get this converted back into a polyline. But it seems to be more complex than just a cast. Can someone please help me? I cannot find anything online.


title: Patents in GeoSpatial software: impact on practices?

text: I've seen on twitter a tweet saying: > Palantir has done well turning chapters of intro computational geometry, > GIS, and ontology textbooks into patents: Looking into the results, my main concern is about the patent concerning "Simplifying a polygon" (I use polygon simplification with the MapShaper tool or the Douglas Peucker algorithm in PostGIS or GDAL ...) Do you know the impact of this patent in particular? I have to stop using polygon simplification?


title: Would PostGIS offer an advantage over MySQL for a produce farm application?

text: I have a web app that stores the locations of farms in West Michigan. You can search for a product (e.g. "broccoli") and it will show you all the farms that grow that product. Right now I'm using MySQL and using trigonometry to calculate the difference between the user's location and the location of each farm. It's not a bad way to go but it did take some doing. Another thing I want to do soon is to map out growing seasons for different products for different regions. (For example, I want to show that avocados grow at a certain time of year in California but never Ohio.) I realize this is an open-ended and possibly naive question, but might it be worth it for me to make the switch to PostgreSQL/PostGIS to take advantage of its spatial capabilities?


title: What does the "Standard Distance" output from the spatial point pattern analysis in QGIS 2 represent?

text: I have a cluster of points and I am using the Spatial point pattern analysis tool from the processing toolbox in QGIS 2. I am having trouble understanding what the "Standard Distance" output is? It looks to me like it represents two standard deviations from the mean as it covers about 84% of my data, but this is only a guess. I clicked on the help tab of the tool but there were no details there. Any one know exactly what the "Standard Distance" layer (dotted circle) is showing? ![enter image description here](http://i.stack.imgur.com/SyrDF.png) Thanks Ando


title: How to create duplicate features based on related table - many-to-one conversion

text: I need to create duplicate features (in this case parcels), using a related table to populate ID fields. In essence, I need to use a related table of parcel owners to create duplicate parcels, each with its own unique ID from the related owner table. The related table already contains the many-to-one link to parcels, I just want to force a one-to-one relationship between the owners and parcels, by creating a duplicate parcel for each owner record. Thanks in advance for any suggestions! Jodi


title: Source for high-res subscription imagery service

text: I'm developing a diagramming tool using aerial imagery as a reference layer. This tool should allow the user to zoom in quite close to a roadway or intersection. Unfortunately none of the public map services I have tested supports zooming in close enough. Are there subscription map services out there that have imagery with resolution high enough to zoom in quite close to the roadway?


title: Convert polygon feature to line feature

text: I have a polygon shape-file i need to extract the line segments where the polygon features intersect each other and also the segment where they don't. Can we automate this process or can scripting like arcpy be used?


title: creating a buffer in a raster layer

text: I've got a raster layer separated into 3 zones, representative of marine park zones. What i'm trying to do is make it so that where zone 1 is adjacent to zone 3, there is a "buffer" between them, of zone 2. How would I go about this? ![enter image description here](http://i.stack.imgur.com/heMJz.png) i want a small buffer zone of the blue zone, where ever the red and green zones are adjacent to one another. short of cracking open photoshop and drawing it in myself, is there a way of doing this: ![enter image description here](http://i.stack.imgur.com/BD9N2.png) in arcmap?


title: Create moving map application from topographic maps

text: I'm trying to create a moving map application for an embedded ATOM device. Data source is a 400MB tif file (500dpi) plus projection and worldfile. My naive approach was to load that tiff into Global Mapper and export it to "Bing", using 4 zoomlevels. It works fine, but the tilecache uses up way too much diskspace. Can anyone point me to a better way of achieveing something like this? Would it be better to use some kind of runtime decompression?


title: Split polygon into lines

text: I have a polygon layer which looks like this ![enter image description here](http://i.stack.imgur.com/MxgNl.jpg). I used the QGIS `multiparts to singleparts` tool and also `polygon to lines`. In both cases i got this type of result ![enter image description here](http://i.stack.imgur.com/61XJ7.jpg). I want to extract polygon into individual lines like this![enter image description here](http://i.stack.imgur.com/3Hk0c.jpg)


title: What Avenue scripts/extensions do you still keep ArcView 3.x around for?

text: Alternatively, if you've completely moved on, which do you miss the most? Thanks and signed,     Not entirely idly curious. (As a side project, I'm looking into the effort required to (semi-) automate conversion from Avenue to python...) :)


title: STDistance Unit in SQL Server 2008

text: I've got a SPROC in SQL Server 2008 where I'm trying to calculate the distance between a point and a polygon, and I'm getting unexpected results. I am passing a parameter (@Radius) in miles to the SPROC, but the distances don't seem to be in KM. I can't figure out what they are in. My SPROC is below:               DECLARE  @Point geometry;     SET   @Point = geometry::STGeomFromText('POINT(' + CAST(@Longitude AS nvarchar(32)) + ' ' + CAST(@Latitude AS nvarchar(32)) +')', 4326);          SELECT  [FS_Geometries].[Description], ([Geometry].MakeValid()).STDistance(@Point)     FROM  [FS_Geometries]     WHERE  (([Geometry].MakeValid()).STDistance(@Point) <= (@Radius * 1609.344))     AND   [ParentId] = @CustomerId     AND   [ParentType] = 80     ORDER BY (([Geometry].MakeValid()).STDistance(@Point)) ASC      In all my other queries, I need to convert my radius from miles to meters. When I do that here, the results are not correct. I have a result set of about 750 polygons that are separated by about 25 miles of geography, so when I pass "1.0", I'd expect to get less than my 750 polygons. Unfortunately, I get the full 750 in the results. Is the STDistance() method in SQL Server 2008 based on something other than meters? Also, is the STDistance() method in this case measuring from the edge of the polygon or the centroid of the polygon? Thanks! UPDATE --- My updated conversion function...               public static double ConvertMilesToDegrees(double miles, bool isLatitude)     {         return isLatitude             ? miles * (((double)1) / 69)  // 1 Deg. Latitude = 69 Miles...             : miles * (((double)1) / 61); // Avg of 69 and 53...              // NOTES: A degree of longitude is widest at the equator at 69.172 miles (111.321)          // and gradually shrinks to zero at the poles. At 40° north or south the distance          // between a degree of longitude is 53 miles (85 km).     }


title: How To Convert A Set of ARC/Info ASCII .e00, .e01, .e02, etc Files Into .shp Files

text: I have a set files that I downloaded with extensions (.e00, .e01, .e02, etc). I have read that these are files used by ARC Info mainly, but I do not have ARC and nor am I really that interested in downloading it. The funnier thing is that I am able to convert a .e00 file in QGIS to a .shp set, but if the file has been broken into many .e0n files, QGIS does not recognize the others and leaves me with only have my shapefile. I have tried taking the other files and copying and pasting their contents to the end of the .e00 file, but this resulted in a very erroneous shape in QGIS. I have been searching for some while now if there is a way to do this with ogr2ogr, but haven't found anything for a set of files like this. I also downloaded e002shp, but had no idea on how to go about using it and didn't find much online for that one. If anyone can tell me any command-line tools I could use (preferably ogr2ogr) or a way to do this in QGIS, I'd gladly appreciate it :) Actually I'd appreciate any tips, but that is the software which I am more familiar with ;) Thanks, elshae


title: Algorithm to find nearest point

text: I've a list of a few hundred cities with their latitude/longitude. Given another location (also in lat/long) I need to find the nearest city. As I don't use any GIS, by now the obvious algorithm is to make a loop for all the cities, calculating the distance between the points. Making the loop is practicable for me, but there is some easy to implement algorithm to accomplish that more efficiently? Or some light Java library that can help to solve that? **Notes** : I don't need/want a complete GIS solution or heavy/complicated library. I prefer a less good but easiest and lighter solution because that is the only thing that I need to solve.


title: postgis insert query to store geometry and other data

text: I'm using Postgis functions `st_dump` and `st_split` to split a polygon based on line. Now, after splitting; output geometry stores in a database table as well. Here is the query:               insert into aoo(geom)     SELECT ST_AsTEXT ((ST_Dump(ST_Split(circle, line))).geom) As wkt     FROM (SELECT          ST_GeomFromText('LINESTRING(73.33393335342407 33.82005181729661,73.33685159683228 33.80960466650892)',4326) As line,         (select geom As circle from aoo where gid=71) ) As foo;      But, the issue is I want to store other data with these geometries as well. Like               insert into aoo(geom,name)      but I don't what is the way to do it.


title: How to disable a base layer in OpenLayers based on zoom level

text: > **Possible Duplicate:**   >  How can I switch baselayers as a user zooms in to the map in OpenLayers I have a map that has two base layers and several overlays. One of the base layers only has 5 zoom levels. I would like to disable that base layer in the layer switcher when the map is zoomed outside of the range of that base layer. I have found how to do that with overlays in this forum at this link. I changed dataLayers to baseLayers. My problem seems to come when the base layer that I am trying to limit is the base layer selected when I zoom outside the range of that base layer. If I then select the other base layer, it appears as if the zoom control is locked up. Any suggestions are greatly appreciated. Thanks.


title: different maps for different zoom levels in openlayers?

text: I'm trying to create a mixture of maps that will load in depending on zoom level However, it will only show as one map to the user, and depending on the zoom level, it will load in a different map I really like open street map data when zoomed out, and skobbler maps when zoomed in, so would like to create this one combined map I think from what i've read its possible, but all the examples i found refer to detail layers any examples or pointing in the right direction would be appreaciated


title: Problem in calculating length

text: I am Digitizing road line in India using open layer plugin.When I calculate the length of line using field calculator it gives me the result in millimeter . I thought it is the projection problem and defined the projection as wgs 84 UTM 44N . The same problem exixts . but when I measure the line using measure line tool it gives the result in meters . What is the problem? I want the result to be in meters in field calculator itself.


title: Merge shapefiles

text: I am a new user of ArcMap. I have several shape files of various portions of the West coast of India (boundary and contour files). I need to combine all these files into one shape file so that I can extract information from this file and make one bathymetric information file. I have tried using the merge option in ArcMap. But the shape files do not line up like I want them too. How can i solve this issue? Please Help.


title: How to convert ASCII in x,y,z format using gdal_grid?

text: I have a map saved in (AFAIK) Arc/Info ASCII format and I would like to convert it with gdal grid but I am getting errors. My data (`p.csv`) is in the following format:               319053.25,319301.12,544.28     319074.57,319275.26,546.44     319075.00,319250.00,546.34     319075.00,319000.00,468.74     319075.00,318950.00,466.19     319075.00,318925.00,472.06      and I have a `p.vrt` written by myself:               <OGRVRTDataSource>         <OGRVRTLayer name="dem">           <SrcDataSource>p.csv</SrcDataSource>            <GeometryType>wkbPoint</GeometryType>            <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>         </OGRVRTLayer>     </OGRVRTDataSource>      When I launch `gdal_grid` I get the following error :               $ gdal_grid -l dem p.vrt p.tiff     ERROR 1: Failed to find layer 'dem' on datasource 'p.csv'.     No point geometry found on layer dem, skipping.


title: How to convert data type from float to integer in QGIS?

text: Currently I am working on importing data from MSSQL to Qgis 2.4. In imported database table attributes are shown in the format of integer ( in attribute table it showing like 48.23 and 16.20 in column). But when i export it ( save as shapefile ) then it got automatically converted into float ( in attribute table it showing like 48 and 16). I wanted it in integer format only ( 48.23 and 16.20 like that) so my question is does any other way to convert it into integer format in QGIS? Or does any other saving option, any plug in available to solve this issue?


title: best route change the starting time automatically?

text: I want to find out the best route from point A to point Z. But this should always be at a different time. For example:   * How is the best route from point A to point Z at 07:00?   * How is the best route from point A to point Z at 07:05?   * How is the best route from point A to point Z at 07:10? and so on ... Is there a possibility to change the start time every minute automatically? I am working with ArcGis 10.1


title: Download OSM data using "ArcGis Editor for OSM"

text: My actions: 1) Add Basemap -> OpenStreetMap 2) zoom to area (in this area about 300 nodes) 3) OpenStreetMap Toolbox -> PostgreSQL Download, Extract and Symbolize Tools -> Download, Extract and Symbolize OSM Data (PostgreSQL) Then DownloadExtractSymbolizeOSMData222 run There are a lot of messages. But I have 3 warnings:               ...OSMGPAAttributeSelector "...fstutorial.ftuser.TestOSM\fstutorial.ftuser.TestOSM_osm_ply" ALL     The process did not execute because pecondition is false               ...OSMGPAAttributeSelector "...fstutorial.ftuser.TestOSM\fstutorial.ftuser.TestOSM_osm_ln" ALL     The process did not execute because pecondition is false               ...OSMGPAAttributeSelector "...fstutorial.ftuser.TestOSM\fstutorial.ftuser.TestOSM_osm_pt" ALL     The process did not execute because pecondition is false      When this process finished I can not add this "feature dataset" to the Table of Contents. _**Update_** I use ArcMap 10.1, postgres 9.0.5, ArcGIS Editor for OSM 2.1 Final for 10.1 It is the result of downloading OSM. ![not all objects](http://savepic.su/3234373.jpg)


title: How to export postgis points as blocks in dxf with ogr2ogr

text: I'm trying to export postgis points as blocks to DXF format. I have a simple table in postgis with columns:               CREATE TABLE pod_com_azimuth     (       ogc_fid integer,       rotation double precision,       pointtype text,       geom geometry(Point)     )      I'm converting it to DXF with an ogr2ogr utility:               ogr2ogr -f DXF angle_test.dxf      PG:"dbname='coms' host='localhost' port='5432' user='mofo' password='mofo' schemas=public"     -sql "select ogc_fid, rotation as blockangle, pointtype as blockname, 'rblk' as layer, geom      from pod_com_azimuth"      But when I open the resulting file in AutoCAD i get a layer with POINT objects. But I need them to be blocks with a blockname and blockangle values taken from postgresql fields. How do I export them as blocks (inserts) ?


title: How to increase the resolution of a 3D image in QGIS GRASS NVIZ

text: I'm trying to display 3D data in NVIZ within QGIS, but the result appears very low resolution in comparison to the same data in ArcScene. Any ideas as to how to improve the quality would be greatly appreciated. I've read the question on 'How to increase the resolution of 3D image in GRASS NVIZ' - but this is based on GRASS not in QGIS and the solution didn't appear to work for my data. Both raster and elevation model have been resampled to 0.5 m resolution. How the data looks in ArcScene: ![How it should look \(ArcScene\)](http://i.stack.imgur.com/VHc5q.png) How the data looks in NVIZ: ![How NVIZ looks...](http://i.stack.imgur.com/XD61J.png)


title: geotiff conversion to coordinate system with gdal

text: I would like to get the coordinates X Y Z from a geotiff file. I am using the gdal translate for the conversion. Which output format is the more adequate for containing a 3 numpay array X Y Z ?


title: How to show polyline direction in ArcGIS for Desktop (ArcMap)?

text: can someone tell me how I can show the direction of a polyline as it was digitised? I need this to check if my network analysis is obeying one way streets...I know of a tool in Xtools that can do this however my trial period just finshed.


title: How to reuse Thematic layer definitions in Mapinfo?

text: I have a map from table. Let the table be **TABLE1**. Now That I have created Thematic Map from **Column1** of this table. Again I executed the following query               Select * from TABLE1 where TABLE1.column1 ="someVal" into TABLE2     Add Map auto  Layer TABLE2      Now that new table has been created the old thematic layer **does not work**. Again I create new thematic layer from TABLE2 and execute following query               Select * from TABLE1 where TABLE2.column1 = "otherVal" into TABLE2.       Add Map auto  Layer TABLE2      Thus again thematic map I created is lost. I have created buttons to execute these query but each time I click the button the thematic layer is lost. Is there any way that I could preserve the thematic layer and apply it to TABLE2 for **all select query** executed. Actually the motive was to **use thematic map** that was created for **TABLE1**. But as the **points or value change in each selection** and **points of thematic layer remains constant** and lie in same initial coordinate the thematic layer becomes unusable. > What is the right solution for reusing thematic layer in each Query?


title: Get the names of currently visible rasters from a catalog

text: I have a simple win form app with a map control on it. The map has a layer (MADtedLayer) that was initialized from a catalog in a personal database (Access). The files displayed are DTED. I need to be able somehow to detect the file names of the currently visible files (that is, the files that intersect with the current extent of the map) Opening the database with MSAccess, I see that for the catalog, there are 2 tables, "MyCatalog" and "MyCatalog_Shape_Index". In table "MyCatalog_Shape_Index" there are 5 fields, "IndexedObj","MinGX",MinGY","MaxGX","MaxGY". I suppose that maybe I could make use of them but I cannot understand what those fields represent or how their values are calculated. For example, for a file with LowerLeft=19.00'00" - 34.00'00" and UpperRight=20.00'00" - 35.00'00" I get MinGX = 51001618, MinGY=91266936, MaxGX=53688209 and MaxGY=93953528. Any help or clarification would greatly appreciated.


title: Problem opening shape file in qgis

text: I am having a problem opening a shape file using QGIS using 'add vector later'. This is the error message i am getting: /Users/Conrad/AmerindianLands 2/Amerindian_Land_Karasabai.shp is not a valid or recognized data source Can anyone help? I have made sure that it is being extracted from a folder containing the other necessary files with the following extensions (.dbf,.prj,.shp, .sbx,.sbn) Many thanks Conrad


title: How do I use sys.exit(0) in an arcpy script to exit early without having an error message show up?

text: I have an arcpy tool script for ArcGIS 10.0 that has two major functional sections. The user can choose whether or not to run the second section. If the user chooses NOT to run the second section, I simply want to run a cleanup function and exit the script with a sys.exit(0) **without having an error message posted to the tool results window.** There are two major threads here in GIS-SE about exiting arcpy scripts, but the solutions therein don't specifically address the error message. The general structure of the code is as follows:               import sys     ##import arcpy          def CleanUp():         print 'Cleaning up ...\n'          def finish():         CleanUp()         print 'Exiting ...'         sys.exit(0)          do_more = False  #or True ... input from user          #Section 1:  do some stuff     print 'Doing some stuff ...\n'          if not do_more:         finish()          #Section 2:  do more stuff     print 'doing more stuff ...'     CleanUp()      If I run this code in the Python interpreter outside of ArcGIS/arcpy, it works as I'd expect, exiting gracefully with no error message; however, in my arcpy script with this same structure, the script exits, but a SystemExit error message is posted to the tool results window. Is there a way to make the arcpy tool script eat the exception and bury the SystemExit error message?


title: How to generalize coastlines?

text: I have a large set of data points that essentially represent the coastline of the world. I don't know the exact resolution but I'd say it's in the neighbourhood of 5-10m. Imagine a set of points which represent a jagged coastline. At a 1m resolution, we see all the detail and need 50 data points to correctly represent the coast, however, at 10km resolution the coast looks like a straight line and we would only need 2 points. I'm writing a mapping application that will require the data to scale from "full-planet" resolutions to "street" resolutions. I need some help reducing the resolution of my data points, and subsequently the amount of data for wider resolutions. I've been reading about Vector-Tiling and think this will be the best solution (I'm doing something like this already), however, at the low resolutions, my dataset it still huge meaning that rending it takes a [relatively] long time, where as the high resolutions are fast because the 'effective' dataset it small (subset of the whole). I'm trying to determine how to take my large dataset and reduce the scale of it so that my 'full-planet' view's tile data is a manageable size. How do I create the tile sub-datasets from the larger complete dataset? I'm happy to use a tool but I'd rather do it myself in an effort to learn how it's done.


title: Prevent ArcGIS from maintaining extents during resize

text: **Background:** I have a custom ArcGIS application written using ArcGIS Engine 10.0 that at its simplest form is just a map control inside of a resizeable form. **The Problem:** During resizing of the form, the performance is absolutely awful and in some cases the entire application crashes _(depends on quality of machine)_. This is because the `MapControl` constantly refreshes in order to maintain the same extents of the map. **Question:** Is there a way to prevent the `MapControl` from refreshing to maintain its extents?


title: Display DegreesMinutesSeconds on status bar

text: I have used status bar to display x and y coordinates in `esriDirectionUnits.esriDUDegreesMinutesSeconds` form. I have created the code like this statusBarXY.Text = string.Format("{0}, {1} {2}", e.mapX, e.mapY, esriDirectionUnits.esriDUDegreesMinutesSeconds); but it will not display any format like this X=70°0'20.797"E Y= 22°28'.072" N X= 700'20.797"E Y= 2228'.072"N How could I make this? If any one has a solution please send me link or a sample.


title: How to add features to a PostGIS layer using Quantum GIS?

text: I am using Quantum GIS and the PostGIS dbms. I am able to call the PostGIS layer into QGIS but i am unable to add feature to it. I am using the method of adding features to a Shapefile but still it is not getting reflected in the layer, and I don't know how to save it to the PostGIS database if I am able to add the feature to the PostGIS layer. I may not have used some technical term for this problem because i don't know what to call it, I am expressing the problem in my own terms.


title: How to handle plugins of two different QGIS installations on the same computer

text: I have a computer with 2 versions of QGIS. Both versions are running on Win 7 64 bits with 16 Mb of ram. One is Dufour 32 bits and the other one is version 2.1 OSgeo4W 64 bits. I am experimenting using the processing menu on both versions. Nevertheless I noticed that the paths for Grass in my 2.01 32 bits version were the ones of the OSgeo4W install (Processing/Options and configuration/Providers). I changed to the Dufour ones. Then I tried to replicate a procedure in my 64 bit install but the paths were the Dufour ones! My question is, then, how I can handle the settings for each install in order that they are independent from each other and I do not have to change my provider paths each time that I want to use the processing menu.


title: ArcGIS API V 3.5 Navigation toolbar (zoomToPrevExtent) not working

text: I am just getting started with learning Dojo and ArcGIS API so please bear with me. I have an existing prototype of an application that use ArcGIS version 3.5. I am looking at the navigation toolbar sample https://developers.arcgis.com/javascript/jssamples/toolbar_navigation.html and trying to get the Prev Extent and Next Extent functionality work in my application. This sample uses ArcGIS V 3.8. When I copy this sample code and run it locally it works. However since I use V3.5 I changed the below line in this sample and it stopped working.               <script src="http://js.arcgis.com/3.8/"></script>      <script src="//serverapi.arcgisonline.com/jsapi/arcgis/3.5"></script>      Is this feature only supported in V3.8? How do I get it work with V3.5? Any help is immensely appreciated. Thanks very much.


title: arcpy: join feature class (one to one) and get feature id of each joined object

text: I am currently spatial joining two feature classes like this               arcpy.SpatialJoin_analysis(target_features,join_features,output_features,                            join_operation="JOIN_ONE_TO_ONE",                            join_type="KEEP_ALL",                            match_option="CLOSEST")      I need to access the object ID of each join feature that was joined to the target feature. However, in spatial joins, arcpy only adds "JOIN_FID" in one- to-many joins. In a one-to-one join I only get "TARGET_FID" (which seems pretty useless as I knew that already). How can I access (programatically) the feature ID of the joined feature from the output feature class?


title: GRASS GIS: auto-completion fails in commande line

text: I'm running GRASS GIS 6.4.3. on Win 7, SP1, 64 bit. The path to my grass-DB is: D:\grassdb Within the GRASS GIS Commande Line, if I press the tab-key to auto-complete only folders, files and documents from "C:\" are listed and none of the maps from GRASS which should appear: If I type:               r.colors map=D       and press the tab-key after the "D" this appears:               r.colors map="data"     r.colors map="Documents and Settings"      and so on. Instead I'd like to see my layers:               r.colors map=DTM_2012     r.colors map=DTM_2013      and so on. Any hints?


title: Is there any readily available implementation of Mapserver (or Geoserver)- Postgres-Openlayers combination?

text: I am developing an web GIS for map data presentation for public. Is there any readily available implementation of Mapserver (or Geoserver)- Postgres-Openlayers combination or similar? ## More Details: **Overall system requirements:** Uploading, managing and publishing map data **Back-end :**   * Upload geospatial data into the administration interface. create new layers and organize. **Admin should be able to:**   * Organize the mapping layers in a navigation menu   * Edit the name, legend, and   * Choose the order of the mapping layers   * Change the thematic colors of layers   * Customizable icons for map layers in the table of contents ( should reflect in the published map also)   * Should be able to add few attributes from the front end. **Front-end should allow the public to:**   * Add google satellite layer.   * Zoom directly to counties from the navigation pane.   * Turn ON/OFF layer visibility   * Turn ON/OFF layer legend   * Turn ON/OFF layer labels   * Change the layer order(up-down) in the navigation   * Export the map to PDF   * Print layout option(with or without legend, labels etc.) Please suggest.


title: There is web based bicycle routing planner that considers accidents but I cannot find it anymore, can someone help me?

text: I remember it uses Google as underlying map and has also a function to rank preferences of the path using forms in which the user can rate routes characteristics, I also remember that is for Los Angeles (but I am not sure about this). I came across this web map a couple of months ago but I did not save it as favourite. Now I have been searching continually for two days in many ways but I can't find it anymore. Does someone know what I am talking about and please provide me the page URL?


title: Creating Thiessen polygons gives ArcGIS ERROR 000537 error in converting points to Voronoi Polygons?

text: I have a problem with creating Thiessen polygons for a set of points. I am creating them with Create Thiessen Polygons tool in ArcGIS 10.1). I am able to create them in a few steps (for parts of the data) but I would like to create it at once.   The error I get is: > ERROR 000537 error in converting points to Voronoi Polygons It doesn't say much, when click at the error in arcgis help I found this message: > Description A TIN surface is created from input points and is used to > generate the output Thiessen (Voronoi) polygons. > > Solution Check the input point features and make sure that features exist in > this particular feature class. As well, make sure that if the extent > environment is set, the input point features do not fall outside their > coordinate range. But if it works "partially", where the problem may be?


title: Split a feature when intersecting with a feature of another layer using PyQGIS/Python?

text: I have a buffer layer (green polygon) which I want to split to two polygons whenever it crosses a barrier (blue line). I have been trying to use "splitGeometry" method, but I just can't get it to work. My code so far is this :               while ldbuffprovider.nextFeature(feat):       while barprovider.nextFeature(feat2):         if feat.geometry().intersects(feat2.geometry()):             intersection = feat.geometry().intersection(feat2.geometry())             result, newGeometries, topoTestPoints=feat.geometry().splitGeometry(intersection.asPolyline(),True)       Which returns 1 for result( error) and an empty list for newGeometries. Any help is greatly appreciated. ![enter image description here](http://i.stack.imgur.com/FBiNX.png)


title: Converting arcs to straight line segments?

text: I have a large number of geometries, currently stored in Oracle, some of which contain arc segments. I need to convert these geometries into LRS (Linear Referencing System) format so I would like to convert all the arcs to straight line segments. I may have to write something in PL/SQL but, being a lazy sort, I was rather hoping someone had done this before. Doesn't have to be in PL/SQL, I also have ArcGIS for Desktop, MapInfo, GeoMedia and FME if any of those would do it.


title: How to change Overlays checkbox feature in layerlist to radio button

text: How to change the default checkbox button of Overlays in Layerlist to radio buttons. My cliet dont want his users see more than one layer on the map. Is there any way that we can prevernt the posibility to check multiple layers in overlay layers in LayerList. Thanks Karteek


title: Generating vineyard's rows in a polygon

text: Does anyone know if it's possible to generate line (the distance between each one is 2 meters) in a polygon?


Attempting to install rgeos package fails. Detailed error logs as follows:-               > install.packages("rgeos", repos="http://R-Forge.R-project.org")     Installing package(s) into ‘/Users/calvin/R/x86_64-apple-darwin11.2.0-library/2.14’     (as ‘lib’ is unspecified)     Warning: dependency ‘stringr’ is not available     trying URL 'http://R-Forge.R-project.org/src/contrib/rgeos_0.1-15.tar.gz'     Content type 'application/x-gzip' length 166974 bytes (163 Kb)     opened URL     ==================================================     downloaded 163 Kb          ERROR: dependency ‘stringr’ is not available for package ‘rgeos’     * removing ‘/Users/calvin/R/x86_64-apple-darwin11.2.0-library/2.14/rgeos’          The downloaded packages are in         ‘/private/var/folders/4l/rd4wl50x52n504t2vzqskhhc0000gn/T/RtmppBNwSt/downloaded_packages’     Warning message:     In install.packages("rgeos", repos = "http://R-Forge.R-project.org") :       installation of package ‘rgeos’ had non-zero exit status     > install.packages("stringr", repos="http://R-Forge.R-project.org")     Installing package(s) into ‘/Users/calvin/R/x86_64-apple-darwin11.2.0-library/2.14’     (as ‘lib’ is unspecified)     Warning message:     In getDependencies(pkgs, dependencies, available, lib) :       package ‘stringr’ is not available (for R version 2.14.0)      Suggestions?


title: Instantiating spatial polygon without using a shape file in R

text: So, the usual way we read a shapefile in R is via the maptools package, like this:               sfdata <- readShapeSpatial("/path/to/my/shapefile.shp", proj4string=CRS("+proj=longlat"))      However, I have a use case whereby I don't have a shapefile.shp but instead I have a series of polygon coordinates               16.484375 59.736328125,17.4951171875 55.1220703125,24.74609375 55.0341796875,22.5927734375 61.142578125,16.484375 59.736328125      and its corresponding projection               coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0       How do I "instantiate" sfdata (which will be a "polygon object") directly from this data? (without going in a roundabout way of creating a shapefile with these data and then reading from the newly created shapefile)


title: How to convert projection

text: I got sample shape files with the following *prj file content: > > *PROJCS["RGF93_Lambert_93",GEOGCS["GCS_RGF_1993",DATUM["D_RGF_1993",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",700000.0],PARAMETER["False_Northing",6600000.0],PARAMETER["Central_Meridian",3.0],PARAMETER["Standard_Parallel_1",44.0],PARAMETER["Standard_Parallel_2",49.0],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]* trying to convert the projection using the shp2pgsql command _shp2pgsql -s 4326_ kept the original projection and didn't convert it as requested. I also check spatial_ref_sys table srtext field for similar known projection unsuccessfully . How can I convert it to WGS 84 projection?


title: Coloring pixels based on values

text: I have an image where pixel values are integers, they go from 0 to 6 and I also have some 999. I want the value 999 to be transparent, which I could achieve by myself. I would like values 0 to 6 to be different colors. I have no idea how to do this. Cant anyone help me? thank you!


title: How to export shp to kml?

text: I already browsed your FAQs, and could not come up with a good answer, I am trying to export shp files to kml to be placed on Google Earth. On a previous version of QGIS it was as simple as right click and save-as. When doing the same steps on QGIS 1.8 everything works fine until opened in GE. My data is in northern Ontario, Canada, and it ends up being placed in the middle of the Arctic Ocean on GE. Is there anything I am doing wrong, the datas CRS is NAD-83 UTM zone 17N. Could it be the projects CRS. is there a way to change a project CRS. On a previous version of QGIS I had no problems, but seem to be running into difficulties on 1.8. Thanks for the help.


title: How to build a mobile app from Flex web app

text: We have a Flex application based on ArcGIS flex API in the final stage of development and ready to move to testing/production. The application uses some map services and a complex GP tool. My boss wants a quick mobile app out of it. I have never worked with mobile app. I want to learn- what it takes to convert this application to mobile app- software, coding effort. Some useful links related to ESRI or Adobe will be helpful. Thanks


title: ArcGIS Explorer Map Display Issue

text: I bet this is probably an issue with my PC, but I have a really odd error with ArcGIS lately. Refer to the image, it'll explain it best: ![ArcGIS Explorer Issue](http://i.stack.imgur.com/0muiC.png) Notice how their no map displayed? No matter what I do I can't get it to display. If I add data, change the basemap, restart my computer, re-installed Explorer. This started happening about a week ago for no apparent reason. Nothing had crashed or changed recently on my PC. Has anyone ever heard of any bugs that might be the cause of this? My ArcGIS Explorer build is: 1750 (2.0.0.1750) Any help would be greatly appreciated.


title: Calculating slope cell-by-cell in ArcPy

text: I am relatively new to the ArcGIS-Python environment and want to calculate elevation raster slope. But it says "ValueError: setting an array element with a sequence". Below code is what I have. Any help would be highly appreciated.               import arcpy, numpy, math     from numpy import *     from arcpy import env          arcpy.env.workspace = ("D:/input")     inRas = arcpy.RasterToNumPyArray('t1')     [rows,cols]=inRas.shape     slope=zeros((rows,cols), float)          desc = arcpy.Describe('t1')     cellSize = desc.meanCellHeight     extent = desc.Extent     pnt = arcpy.Point(extent.XMin,extent.YMin)          for j in range(0,rows-1):         for i in range(0,cols-1):             slope[i,j]= inRas[math.sqrt(((cols+1)-(cols-1)) ** 2 + ((rows+1)-(rows-1)) ** 2)]             print slope             myArray=arcpy.NumPyArrayToRaster('slope',pnt,cellSize,cellSize)                  myArray.save("D:/output/t1")


title: How to select features in ArcGIS based on their relation to other features within same layer?

text: I want to select all the features in my Waterlines layer that share a line segment with other features in the same layer. I've done a complicated search cursor iteration in python that takes a couple hours, but I'd rather have something quicker. I found this page that says you can do it simply with "select by location" which doesn't seem to work from the selection menu in arcmap, but kind of works in arcpy. http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Select_By_Location_Select_based_on_spatial_relationship_within_the_layer/0017000000tr000000/ When I ran this:               arcpy.SelectLayerByLocation_management("waterlines_layer", "SHARE_A_LINE_SEGMENT_WITH", "", "", "NEW_SELECTION")      ...it selected things that definitely didn't share a line segment with anything else. But it didn't select everything(which would have at least made some sense). Surely this is possible without buying that expensive arcmap extention for QA/QC.


title: Project one layer in qgis using projection CRS of another layer

text: I am using QGIS 1.6 and have imported a set of latitude, longitude points from a csv file as a layer. I also have another layer from a shape file. When I go to layer->properties for the second layer I see that the "Layer Spatial Reference System" is +proj=lcc +lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs I would like to re-project the first layer using the 2nd layer's CRS. How can I do this? When I enable on the fly projection it re-projects the second layer into latitude longitudes. What I want is the opposite: I want the lat longs to be projected using the CRS of the second layer. I know that I can choose a CRS when I am saving the first layer but I do not know what CRS I should select. Thanks


title: using python and arcpy , how can I compute the direction and angle of a line segment?

text: if start_point_x < end_point_x and start_point_y > end_point_y:         slope = ((end_point_y - start_point_y)/(end_point_x - start_point_x))         radian = math.atan(slope)         degrees = math.degrees(radian)         degrees = round(degrees,2)         if degrees < 0:             degrees = 90 + degrees         print "Going SE"     elif start_point_y < end_point_y and start_point_x < end_point_x:         slope = ((end_point_y - start_point_y)/(end_point_x - start_point_x))         radian = math.atan(slope)         degrees = math.degrees(radian)         degrees = round(degrees,2)         if degrees < 0:             degrees = 90 + degrees         print "Going NE"     elif start_point_x > end_point_x and start_point_y < end_point_y:         slope = ((end_point_y - start_point_y)/(end_point_x - start_point_x))         radian = math.atan(slope)         degrees = math.degrees(radian)         degrees = round(degrees,2)         if degrees < 0:             degrees = 90 + degrees         print "Going NW"     elif start_point_x > end_point_x and start_point_y > end_point_y:         slope = ((end_point_y - start_point_y)/(end_point_x - start_point_x))         radian = math.atan(slope)         degrees = math.degrees(radian)         degrees = round(degrees,2)         if degrees < 0:             degrees = 90 + degrees         print "Going SW"     elif start_point_x == end_point_x and start_point_y > end_point_y:         degrees = 270     elif start_point_x < end_point_x and start_point_y == end_point_y:         degrees = 0     elif start_point_x == end_point_x and start_point_y < end_point_y:         degrees = 90     elif start_point_x > end_point_x and start_point_y == end_point_y:      The above code is 5 degrees off on the 20 degrees in the SW and NE directions respectively


title: boundless SDK Limit layer by scale

text: I’m adding some Layers into an existing application using last version of Bounless suite SDK. Those layers represent very "small" features compared to the total extent of the map so I’ll like to force them to be unselectable until the map is past a certain scale, just as you can do using the geoexplorer tool. I´ve tried using minScale and maxscale on the layer configuration and nothing… So far I´ve been unable to do it and can’t find any page telling me what´s the correct attribute to do this. Any help you guys can give me is quite welcome....


title: How to do raster conditional comparisons in .NET?

text: I am comparing a raster with a constant value by raster calculator. Now I want to do this in a c# code. **Question** Is there any process to do this kind of raster conditional comparison in .net?


title: how to read/write XML file using Openlayers?

text: I am new to OpenLayers. I am trying to load data from a datasource and show it on a OSM map, then also add the functions such as tagging and making new features as well as editing their attributes. My questions are: 1- whats the best (efficient in means of speed) solution for data storage? should i store data in DBMS such as MySQL, or how is to store points in XML (GML) file and read and write the file?? How many users can load and save xml files at the same time? 2- If I want to do it with xml, how can I do it? any good posts that explain this? Thanks!


title: Qgis on Ubuntu - change permitions?

text: First time using (QGIS on UBUNTU). i am following some articles in order to publish qgis maps on the web. My problem is (although i am the administrator) i am not able to save or copy my *.qgs in specific directories like /qgis-web- client/project. I don' t even know where are some directories, since i am really confused with the ubuntu file structure philosophy. Is there any way to change permissions? Is that a good idea?


title: updating the extent automatically in the print composer

text: Just started using QGIS (after trying what I need to in GRASS). Using the Print Composer, I've put down a map which is based on the map in my map canvas, which is what I'll want to print. Then put down a few logos as well, which need to be printed on the map. Now, the people who will eventually be using this map will be panning and zooming around merrily on the map canvas, then expect it to print what is on the map canvas, as well as the logos. Using a saved Print Composer file here will be fine. Trouble is, can the Print Composer automatically update it's extents to match the map canvas? If not, I don't think it's a HUGE problem - I'll just tell them to use the "Set to map canvas extent" button, but it would be nice if it could do this automatically.


title: How to showing a Matplotlib figure in a QGIS plugin

text: I created a plugin for QGIS, But I want it to show a pi chart when you press on a button. But When I do this nothing happens, also when I try to show the chart by typing the code into the QGIS python console nothing happens. Not sure if that is the same problem. Do any one have some advice for me? Just a note, I can show the chart in a QT form that run on it's own (not a QGIS plugin) I found something about Customizing matplotlib to make it work, but I have no idea what to do and if this is even the problem. Thanks


title: Configuration for GeoWebCache in Geoserver 2.4.0

text: I have developed a web map application by using Geoserver 2.4.0 and tomcat 7 in Windows XP. I have tried to configure the GeoWebCache for WMS layers to speed up the service but unfortunately it didn't work for me. I did something wrong in configuration steps. Kindly Help me to achieve this. I have attached some screenshots for reference. ![enter image description here](http://i.stack.imgur.com/JyJO3.png) **I have created New Gridset for EPSG:404000** ![enter image description here](http://i.stack.imgur.com/Cwced.png) **I have tested with XXX.XXX.X.XX:8080/geoserver/gwc/demo Link:** ![enter image description here](http://i.stack.imgur.com/fJMjc.png) **My layer setting is:** ![enter image description here](http://i.stack.imgur.com/FOOTN.png) **And My Geoserver Log is:** at org.apache.catalina.valves.AccessLogValve.invoke(AccessLogValve.java:953) at org.apache.catalina.core.StandardEngineValve.invoke(StandardEngineValve.java:118) at org.apache.catalina.connector.CoyoteAdapter.service(CoyoteAdapter.java:408) at org.apache.coyote.http11.AbstractHttp11Processor.process(AbstractHttp11Processor.java:1023) at org.apache.coyote.AbstractProtocol$AbstractConnectionHandler.process(AbstractProtocol.java:589) at org.apache.tomcat.util.net.JIoEndpoint$SocketProcessor.run(JIoEndpoint.java:310) at java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source) at java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source) at java.lang.Thread.run(Unknown Source) 2014-01-18 10:30:31,901 ERROR [geowebcache.GeoWebCacheDispatcher] - Problem communicating with GeoServer http://xxx.xxx.x.xx:8080/geoserver/gwc/service/wms 2014-01-18 10:30:31,901 INFO [geoserver.wms] - Request: getServiceInfo 2014-01-18 10:30:31,901 ERROR [geoserver.ows] - org.geoserver.platform.ServiceException: No such style: ODC_Raster at org.geoserver.wms.map.GetMapKvpRequestReader.parseStyles(GetMapKvpRequestReader.java:1386) **My Default Raster Style:** ![enter image description here](http://i.stack.imgur.com/kuicf.png) **My Custome Raster Style:** ![enter image description here](http://i.stack.imgur.com/ca36V.png) Thanks in advance


title: Can I use the LEFT (string) operator in ExecuteSQL()?

text: In QGIS I can use the `LEFT` operator on a dataset like this:               LEFT( "Name2",3) ='031'       What is the syntax for using this expression with ogr in python? Isn't this the right syntax?               a=water.ExecuteSQL("select left('Name2',3)='031' from 'benlux_w_min' ")     'ERROR 1: SQL Expression Parsing Error: syntax error`      A simple where clause works just fine               In [53]: a=water.ExecuteSQL("select * from 'benlux_w_min' where 'Name2'='0310027' ")          In [54]: a.GetFeatureCount()     Out[54]: 367


title: Mapscript. Converting pixel coordinate to latitude longitude for postgis

text: I have a map created with mapscript. It is a 2D (PNG format) that's being shown on a web browser. I want to get the lat/long coordinates whenever the user click on the image, and then save it in my postgresql database. This is what I have done so far:   * Get the current extent.   * Long: minx + ((maxx-minx)/image_width)).   * Lat: miny + ((maxy-miny)/image_height)). Then using: `SELECT ST_SetSRID(ST_MakePoint({$x}, {$y}), 0)` as point to get the geometry. Unfortunately, the geometry is not on the spot where the mouse is clicked. Any ideas to solve this? thanks,


title: shape file for linked micromap not working

text: I am trying to replicate the example provided in the micromap R package of a linked map with the poverty rate for the 50 US states. I am trying to achieve the same results with a Shapefile from Germany. Here is my code:               library(micromap)     library(ggplot2)          setwd ('C:/Users/Jesus/Dropbox/linked maps/Chapter4')     ger<-readShapePoly("germany3.shp")     edPov<-read.csv('gerpoverty.csv')          statePolys <- create_map_table(ger, 'VARNAME_1')     head(statePolys)     lmplot(stat.data=edPov,             map.data=statePolys,             panel.types=c('labels', 'dot','map'),            panel.data=list('Id1','poverty',NA),            ord.by='poverty',            grouping=5, median.row=F,            map.link=c('Id1','VARNAME_1'))      When I tried to parse it into a flat table for use with ggplot2 the procedure does not work (instead of having 16 observations it has 8000 plus observations). Also when I try to plot the linked map, the following error appears on screen:               Error in `[.data.frame`(DF, , ord.by) : undefined columns selected      Here is the link to the files in a zip file: https://www.dropbox.com/s/c43k755aadvu2z6/germany.rar Any ideas or suggestions? Thanks,


title: Commercial Geocoding API that is as good as Google

text: We are looking for a free or paid version of a Geocoding API that will allow us to use it for commercial purposes that is as good as google. Google's api starts at $10,000 a year which is too much. We expect to do 1000-5000 geocodes per month. I have tried Cloudmade, but have had problems geocoding some addresses that google can do. Does anyone have any recommendations based on experience? Here is how I am using the Cloudmade API: http://beta.geocoding.cloudmade.com/v3/API_KEY/api/geo.location.search.2?format=json&source=OSM&enc=UTF-8&limit=10&locale=en&q=%5Bcountry%3DUnited%20States%5D%5Bstate%3DNY%5D%5Bzip%3DZIP%5D%5Bcity%3DCITY%5DAddressRemoved


title: “ArcGIS encountered a serious application error” after Closing

text: i have developed my own add-in for arcmap 10.1 with VS2010. it works fine in my computer but not for other computers. i mean the program goes fine in everything but after closing Arcmap the error mentioned appears for all other computers except mine. any suggestions ? many thanks


title: Why are some OSM GPX line string endpoints connected?

text: I am working with GPX files from OpenStreetMap to render the public transportation lines of Berlin. I converted them to GeoJSON using the togeojson script. For **some** reason some of the lines are rendered as a "circle". The first and the last point are connected for some reason. The effect can be seen at the GitHub project I linked. What am I doing wrong?


title: How to grow several concave hulls until they completely fill a larger area without overlap?

text: This is my problem: I have a very large point layer with polling station locations across one country; each point belongs to a specific constituency (900k points, 4k constituencies). I wish to use this data to create polygon shapefiles of constituencies, as illustrated in the following image (no overlap between constituencies, and no blank space either). On the left, you see three point clouds (points that represent the green, blue and orange constituencies). You also see the country boundary, three blue points out of place, and one orange point out of place (raw data are not always correct, but the dataset is too large to manually correct such things). On the right hand, you see the expected result. ![from point clouds to adjacent polygons](http://i.stack.imgur.com/Xe3uV.jpg) So far, I have explored two methods. The first is to voronoi all points, subsequently dissolve those with the same constituency variable into one polygon, and cut everything away that is outside country boundaries. Two problems: one, the orange outlier. Second, the voronois do not necessarily extend until the country boundary line, leaving blank space. This is roughly the option chosen in reply to this question: Converting point sets to polygon boundaries, but it still leaves me with the problem of outliers and blank space at some of the edges. The second method was to construct concave hulls (convex ones would overlap). The problem here is again both the orange outlier and even more blank space between the various constituencies. The latter could be solved through some sort of growing algorithm, I suppose - but have no idea of how those work. Overall, my intuition is that the problem could be resolved through yet another strategy, namely some sort of clustering algorithm, which calculates the most likely constituency for every coordinate in the country, on the condition that each constituency needs to be one single area. This would both eliminate the orange outlier and grow each concave hull to maximum extent. I have just no idea how to implement something like this. I would be interested in any automated, open source solution (or strategy) to my problem, be it in QGIS itself, via scripting (python, SQL in PostGIS), custom software. Thank you in advance! **Follow-up questions sparked by first answer:** I split my point data into one file per constituency, and transformed each into a heatmap. I then overlayed these heatmaps in QGIS, with the intention of using Raster Calculator to create a unified heatmap reflecting constituency codes, later to be re-vectorized. So far, all looks promising - until I use Raster Calculator, which only calculates values for overlaps of the original heatmaps - presumably because of NULL values. Any idea how to get around this? The image illustrates the problem. Dark green and red were integrated in raster calculator with "(green@1 >= red@1) * 1 + (red@1 > green@1) * 2" - dark grey is 2 and light grey is 1. The boundary on the left was drawn OK-ish, and more importantly, the single green point to the right was correctly subsumed in the dark grey (=red) constituency. But what about all the areas where either green@1 or red@1 are NULL? ![enter image description here](http://i.stack.imgur.com/IXueN.png) **Follow-up edit:** This can be solved through gdal_calculate.py which does not suffer from the NULL value problem. The final and additional question is whether anybody knows how to automate the heatmap plugin in QGIS (or whether there is a gdal-based equivalent function)? I would hate to do this manually for 4k constituencies... Since this touches on python and qgis qnd the like, the question is separate: Can one batch-automate the creation of QGIS heatmaps? Thank you all for helping me out!


title: How to extract NetCDF raster values to time enabled shapefiles?

text: Does anyone know a good way to extract NetCDF raster values to each data point from a time enabled shapefile at the correct time? I have tried using a python script tool and a loop that steps through time but it seems that the loop only works for the shapefile and not for the NetCDF raster - all shapefile values (from all time step intervals) end up getting assigned raster layer values just from the first time step interval. Both the raster and the shapefile are time enabled with the same time step interval and I have checked with Time Slider that the values of both change through time correctly. Here is what I have tried:               mxd = arcpy.mapping.MapDocument('SST.mxd')     df = arcpy.mapping.ListDataFrames(mxd)[0]     rasterLayerName = "SST"     rasterLayer = arcpy.mapping.ListLayers(mxd, rasterLayerName, df)[0]     plasticLayerName = "Plastic"     plasticLayer = arcpy.mapping.ListLayers(mxd, plasticLayerName, df)[0]     df.time.currentTime = df.time.startTime     while df.time.currentTime <= df.time.endTime:        uniqueName = "Plastic_" + df.time.currentTime.strftime("%Y_%m")        outPointFeatures = "U:/gis/" + uniqueName        arcpy.sa.ExtractValuesToPoints(plasticLayer, rasterLayer,              outPointFeatures, "NONE", "VALUE_ONLY")        df.time.currentTime += df.time.timeStepInterval


title: Uploading maps made using ArcGIS for Desktop to Google Maps?

text: I've made a few maps using ArcGIS for Desktop. I wanted to know if it was possible to upload these maps onto Google Maps. Is there a way to do this?


title: Where can I obtain Water Virtual Flow Segments for Minnesota?

text: This is a cross-border task. I already have this data for Ontario, Canada but I also need to complete the model for a large lake that is essentially split in half between the two countries. Definition of Water Virtual Flow Segments: Depicts the flow of water from one stream to another through a lake or wetland. Here is an example of the data I have (Ontario Ministry of Natural Resources): I am looking to find a comparable Minnesota dataset. ![enter image description here](http://i.stack.imgur.com/Ys4Eo.png) ![enter image description here](http://i.stack.imgur.com/pkfEx.png)


title: Get EditText value with event handler in MapBasic

text: I have a dialog in MapBasic that includes an EditText control, a PopupMenu control and a button. When the user presses the button a sub is called. In this sub I would like to get the values of the EditText and PopupMenu, do some maths with them and then output the results into a second EditText control (which will be disabled so as to be used for output only). My question is how do I get the values from the EditText and Popupmenu controls?


title: Is there a way to syncronize the current view extent in QGIS with Google Earth?

text: I already have the openlayers plugin but it's not working has I need it and I wanted to have the same view on both QGIS and Google Earth. Any ideas? Update: I'm using the openlayers plugin which enables me to add the google iamgery as a layer in the TOC. But it takes a long time to load and messes up the scale i'm working with. So i would like to have qgis with my layers, and simultaneously have google earth displaying the same extent view.


title: ApplyEdits in a loop not working as expected in argis flex API version 3.6

text: I have what is hopefully a simple question about the applyEdits method for editing a featureLayer in the ArcGIS API for flex version 3.6. I have an array of data derived from a csv file that is parsed and ready for submission using the applyEdits method, however the array is sometimes too large (over 40,000 rows to data) which ends up taking an unbelievable amount of time to return during which time the application would timeout automatically. I decided to array slice the array and apply edits per thousand records. Unfortunately, it seems that this has no effect at all in educing the time that it takes to apply those edits. The following is a code snippet of the original               applyEdits(arrayOfManyRecords, null, null, true, asyncResponder(successCallBack, failCallback));      The following is the ineffective solution i came up with               applyEdits(arrayOfManyRecords.slice(currentIndex,limit), null, null, true, asyncResponder(successCallBack, failCallback));      The current index starts at zero and will be incremented by the limit (which is 1000) after each return of the apply edits function in a loop until all records have been exhausted. Unfortunately, this does not work as expected and i'm pretty sure it's not due to any sort of syntactical or logical errors though i am open to any suggestions at this point. The problem is that the apply edits function seems to take as much time on the first iteration of 1000 records as it would if it were apply the edits to the entire 40,000 records. If/once this first operation completes, it then takes just fractions of seconds to apply edits to the remaining sets of 1000 records in the array. I am starting to think that the applyEdits method is optimized in someway to be invoked for each possible applyEdits call within the current activation object of the function from which it is called even though this doesn't make much sense as i certainly would not have written the method to behave this way had i been it's author. I'd really appreciate any help on this issue or even just some discussion that would lead to an eventual conclusion/closure (in the non-technical sense of the word) so that i can understand why this is happening. Aside from ensuring that the application does not timeout, an implementation per 1000 records can provide a hugely significant feed back to a user for an operation that takes a lot of time to complete. This would be done in the form of a notification window with an indication of current progress each time the method returns.


title: How to import/copy geodatabase schema from multiple tables to a single table

text: I have an existing ArcGIS 10 geodatabase which has a complicated schema containing multiple tables and joins, and I wish to simplify this for faster web delivery. For example, Featureclass1 is related to Table1, Table2, TableN, and I wish to copy or import the attributes from each Table onto Featureclass1. This will result in a single Featureclass containing all attributes. **What is the fastest way to do this in a repeatable manner** (ie, without doing it manually in ArcCatalog)? This I'm considering:   * When creating a new featureclass or table in ArcCatalog, there is the option to Import the schema from an existing table/FC. However, this only allows you to import from Table1 (if you then choose Table2, it drops the fields from Table1).   * The Geodatabase Diagrammer is a great tool for reporting on an existing schema, and designing a new one - but I can't see any options to import multiple fields simply (it seems to require manual entry of each new field).   * writing a Python script to do it iteratively. Is this the best option? Note that I'm only concerned with the _schema_ at this stage, not the actual data import. Thanks


title: Using Gui tools and arcpy.mapping

text: For a project I want to create a little program that does simple things like export a map from the python command line, but also users are able to edit said map (without having to open arcmap) e.g change font text, replace a street name with another street name, change map layout. I know you can do something similiar in arcpy.mapping from the command line - but what if I use gui tools to create a small program that allows users to do so by pressing buttons (e.g one button exports map, I could then upload a few example maps that can be worked on) - I'm looking at using tkinter. Does this sound overly complex? I can remove the gui tools part if its too hard - is it possible to test out python programs somewhere? - like can I upload it somewhere and test it out and do I need html knowledge? thanks   * on a second note - does anyone know a way to replace a graphic in all maps - for example, if I want to replace the company logo of all my maps with another logo I've made - I've googled and can't seem to find an answer to it.


title: How to use ColorBrewer in ArcGIS 10?

text: I've been using ColorTool (based on supreme ColorBrewer) for most of my mapping needs in ArcMap. After migrating to ArcGIS 10 I discovered that it will no longer be a case. Are there any other options to include ColorBrewer palettes in ArcGIS 10?


title: How to add address block on Google Mapmaker?

text: Is it possible to mark a block of addresses for example B-block, xyz colony, which is a logical block of addresses on Google Mapmaker? If possible, how do I do it? I want to add a polygon feature and select the block but in the options there is no 'block' choice or anything similar.


title: Define Projection for an autocad transferred data

text: If I had a cadastral map on an Autocad file and i convert it to a dxf extension in order to open it in QGIS. I Added The Parcels borders and text in QGIS. Can you please list the steps which I should take starting with defining projection for the streets and points i have(please know that I need them to be in UTM84-36N) and ending by having 2 UTM shapefiles having a lebanese UTM coordinates (Ex: X=734232 ; Y=3749897) Thank you


title: Flag More Than One Space in SQL

text: I am working on an XP machine with FME'13 and I have a _Address_ field and I have to find the spaces(middle spaces and not the leading and trailing ) in each and every row. **St._Petters_Road** is the valid space(it need not be flagged); but **St.__Petters_Road** is invalid (as it has two consecutive whitespaces). It is assumed that the words are correct. But there should not be two two consecutive spaces. and I need to flag the invalid attributes in a new column as 0 and 1 Can we do that with transformers or in SQL?


title: UpdateCursor on two tables with value from SearchCursor

text: I try to keep myself busy and learn some python. I have two tables test2.shp (for searchcursor) ![1](http://s29.postimg.org/j239hkvtz/test2.png) test.shp (needs to be updated) ![1](http://s18.postimg.org/brc1abasp/test.png) I want to update `test.shp` with a `UpdateCursor`. The column which needs to be updated is `column2` if `column1` of `test.shp` satifies a certain criteria `(column1 = aa)`. So, `if column1 = aa`, insert value in `column2`. I try to take the values for the update from `test2.shp` with a `SearchCursor`. That's my code:               import arcpy          arcpy.env.workspace = "c:\\Users\\python_ex"          scursor = arcpy.SearchCursor("test2.shp", "column1" = 'aa', "", ["column1", "column2"], "")     ucursor = arcpy.UpdateCursor("test.shp", "", "", ["column1", "column2"], "")          for srow in scursor:         if srow == "aa":             a = srow.getValue("column2")         else:             continue             for urow in ucursor:                 if urow == "aa":                     urow.setValue("column2", a)                     cursor.updateRow(urow)                 else:                     continue             del ucursor, urow     del urow, srow      I get the error: `Keyword cant be an expression` where I define the `scursor`. But I guess that is not my only mistake... I am working on Arcgis for Desktop 10.2.2 Any ideas? EDIT 8/7/14 12:20pm Alright, thats what I came up with:               import arcpy          #set to folder where features are located     arcpy.env.workspace = "c:\\Users\\python_ex" #on windows use \\ instead of /               #---------------------------     #define variables for cursor     #---------------------------               """SEARCHCURSOR"""     """SearchCursor (in_table, field_names, {where_clause}, {spatial_reference}, {explode_to_points}, {sql_clause})"""     #field taken as "identifier". i.e. if column1 == a, take the value from the same row of column2     stable = "test2.shp"     sfield = ["column1", "column2"]     #where clause will be written in SQL and not assigned to variable               """UPDATECURSOR"""     """UpdateCursor (in_table, field_names, {where_clause}, {spatial_reference}, {explode_to_points}, {sql_clause})"""     #field taken as "identifier". i.e. if column 1 == a, take the value of the SearchCursor and insert into same row, but column2     utable = "test.shp"     ufield = ["column1", "column2"]     #where clause will be written in SQL and not assigned to variable               #--------------------------     #start the loop     #--------------------------          with arcpy.da.SearchCursor(stable, sfield, """"column1" = 'aa'""") as scursor: #SQL queries in python are build that way http://resources.arcgis.com/en/help/main/10.2/index.html#//002z0000001r000000         for srow in scursor:             print srow             with arcpy.da.UpdateCursor(utable, ufield, """"column1" = 'aa'""") as ucursor: #if used with statements, no need to close all the cursors                 for urow in ucursor:                     ucursor.updateRow(srow)      BUT, that only takes the last the last found value of `aa` in `test2.shp` and updates all the other `aa` in `test.shp` with this value!


title: "Back calculating" to quantify the "atmospheric effect" of water vapor or other gases in a "corrected" hyperspectral image using QUAC?

text: I know that when using FLAASH in ENVI, you can produce a water vapor map (caveats considered) from your atmospheric correction. My question is if it is possible to "back-calculate" a water vapor map if I used the QUAC atmospheric correction? Possible for other gases (e.g. methane or CO2)? I'm not so much interested in the accuracy of the output as much as I'm interested in the methodology. I would like to understand better, if only in theory, how FLAASH calculates its water map and if its possible for me to try my own calculation of atmospheric contributions for each pixel from the QUAC correction, which considers only 'in scene' radiance values.


title: Extract Polygon Data for Each Neighborhood from GIS Shape File

text: I have a dataset that has four types of files with the following file extensions - .DBF,.PRJ,SHP,SHX My goal is to get all the polygon coordinates for each neighborhood and draw the boundary outline on the google map when a user is on a specific Neighborhood page. Currently the datasource for the site is a MySQL Database. I have tried converting the shp file to Geojson but then parsing that file is time taking. I have also tried using tools such as QGIS and shp2sql etc. I would like to know what is the best way to get the main attributes from the dataset along with each neighborhood's Polygon coordinates. please advise


title: Is there a way to take control of the mouse pointer using python in QGIS?

text: i would like to take control of the mouse pointer in QGIS (for example constraining it to a given extent...)   Is there any way to do that using PyQGIS ?   Thanks !


title: Function ST_Transform postgis 2.0

text: I've read many forum posts on the function ST_Transform but unfortunately was unable to resolve my doubts so I posted a new question. I have a spatial table SRID = 101010 and want to transform it to SRID = 101020, both already on the table my spatial_ref_sys Initially I did the following SQL and it worked.               SELECT         ST_Transform (geom, 101020)     FROM         lotacao      Later I tried to update my geom column               UPDATE lotacao     SET geom = ST_Transform (geom, 101020)      And the error occurs:               [Err] ERRO:  Geometry SRID (101020) does not match column SRID (101010)


title: turning multiple vector layers into one raster in arcgis 9.3.1

text: I want to make a raster out of the multitier view I have in ArcMap 9.3.1 . The layers in my ArcMap come from several datasets. I know there is a 'Feature to Raster' tool, but it converts one layer at a time. Is there a way to convert all of the layers to a raster ( asides from export as jpg, and georefrencing the jpg) ? ( I have a faint memory telling me this a simple question with a very simple answer, but I honestly don't remember)


title: AutoCAD to OSM, lossless

text: I'm trying to convert AutoCAD 3D data to OSM lossless, or close to lossless. So far, the only way I can see it's possible is to convert from AutoCAD to SVG to OSM, according to OSM's WiKi. Do someone have experience in using AutoCAD 3D data for OSM? I'm going to use it in OSM3S, own hosted overpass-api. I have the original data in MicroStation, so it's actually going from MicroStation to OSM.


title: Can a layer's style be modified dynamically using an SLD in OpenLayers?

text: I am a complete novice with GeoServer and OpenLayers. I'm trying to ascertain at what point a style is applied to a GeoServer layer which is drawn in OpenLayers, when using an SLD. The Introduction to SLD topic says: > When adding a layer and a style to GeoServer at the same time, the style > should be added first, so that the new layer can be associated with the > style immediately which implies that the style is stored with the layer in GeoServer. Does this mean that I can't change the layer's style dynamically in OpenLayers? The Attribute-Based Polygon example shows how to symbolize a polygon based on its attributes. Could a user change this symbology dynamically, within an OpenLayers map which is running on their browser? The OpenLayers Styling page says: > When a feature is added to a map, its style information can come from one of > three sources: Does this mean that it would be necessary to remove a layer, and add it to the map with a new SLD, in order to change its style? Thanks


title: Adding another style to an OpenStreetMap server using Mapnik and PostGIS etc

text: I have created an OpenStreetMap server using this guide: switch to osm it's very good and works like a charm. I have also been able to add a new theme using TileMill and osm-bright However, I now want to be able to switch between two themes, osm-bright and my own osm-night. Is this an apache question or is there a way to modify mapnik to be able to get a parameter from for example openlayers to switch between these two themes? Gratefull for any advice.


title: Summarizing Table on 2 fields

text: I have a table of Census Blocks linked to ZIP4 values. I'm trying to create another table getting the unique ZIP4 values associated with each Census Block. I've tried using the Frequency, and Summary Statistics tool. These give me what I need, but the ZIP4 values are all added up, I need them separately. Would using Pivot Tables get me what I need? Using ArcMap 10.2.2, all extensions.


title: CartoDB create layer at runtime

text: I'm almost totally new to CartoDB. I've been looking around and read the api reference, but there seems not to be what I'm looking for. What I want to do is to "Create visualizations at runtime". From what I read, it is possible to do it only if you already have a visualization created with the UI and then you can work with sublayers. But what I actually want to do is to create a visualization from scratch, using just the table containing the values I want to represent on the map. Is this possible or I'm totally out of mind?


title: Display Point Labels in QGIS

text: ![Question](http://i.stack.imgur.com/nP5fP.jpg) I have coordinates and points from that coordinates. so what do i need to do to show point numbers on the map? Thanks. It's my second day with QGIS