title: Reading zipped ESRI BIL files with Python

text: I asked the question below on StackOverflow, but haven't received any feedback there yet: I have precipitation data from the PRISM Climate Group which are now offered in .bil format (ESRI BIL, I think) and I'd like to be able to read these datasets with Python. I've installed the spectral package, but the `open_image()` method returns an error:               def ReadBilFile(bil):         import spectral as sp         b = sp.open_image(bil)     ReadBilFile(r'G:\truncated\ppt\1950\PRISM_ppt_stable_4kmM2_1950_bil.bil')      `IOError: Unable to determine file type or type not supported.` The documentation for spectral clearly says that it supports BIL files, can anyone shed any light on what's happening here? I am also open to using GDAL, which supposedly supports the similar/equivalent ESRI EHdr format, but I can't find any good code snipets to get started.


title: Active Directory Authentication & MS SQL Server in ArcReader

text: I'm having difficulty with a single user (That I know of) being able to access geodatabase files in our ArcGIS Server using ArcReader. Our setup is as follows: Enterprise Geodatabase within MS SQL Server 2012, Active Directory integration within SQL Server (GIS Edit group & GIS View group), mixed environment of ArcGIS users and viewers. Every other GIS user with either ArcReader or ArcView can access the geodatabase data, but this new user I have added to the ArcView group cannot see the layers. The PMF map shows all of them as unavailable. Her ArcReader has the server connection in the proper folder to allow for OS authentication, and I even tried putting her in the edit group but she still cannot access the layers. We're running 10.2 server and she has 10.2.1 ArcReader.


title: QGIS Small multiples / animate across columns

text: I want to turn each column of a shapefile into choropleth. I have over 100 columns, so rather than manually classifying each column and exporting an image, I was wondering if there is a way to automatically produce small multiples (one map for each column) or animate through each column (one frame per column). MMQGIS does not do this, to my understanding. The columns are not time based, so TimeManager wouldn't work either. This question has been asked before (here and here) without answers. It's possible that I'm approaching the problem the wrong way too. How can I easily create many choropleths from the same shapefile with varying values?


title: pyshp problem with naturalearthdata-shapefile

text: trying to follow this example: https://code.google.com/p/pyshp/wiki/ShapeRecords entering the following code:               import shapefile     r = shapefile.Reader("ne_110m_admin_0_countries")     sr = r.shapeRecords()      gives this error-message under python 3,3 64 bit on windows 7:               Traceback (most recent call last):       File "D:\stuff\beispiel1.py", line 3, in <module>         sr = r.shapeRecords()       File "C:\Python33\lib\site-packages\shapefile.py", line 553, in shapeRecords         for rec in zip(self.shapes(), self.records())]       File "C:\Python33\lib\site-packages\shapefile.py", line 525, in records         r = self.__record()       File "C:\Python33\lib\site-packages\shapefile.py", line 501, in __record         value = u(value)       File "C:\Python33\lib\site-packages\shapefile.py", line 60, in u         return v.decode('utf-8')     UnicodeDecodeError: 'utf-8' codec can't decode byte 0xf4 in position 1: invalid continuation byte      Any ideas, what I am doing wrong (just used another shapefile from www.naturalearthdata.com) ??? **UPDATE (May, 8th. 2014)** Today I received an email from the developer of pyshp, Mr. Joel Lawhead. The problem I mentioned was a bug, which only occurs, when you use python vers 3.x . python vers. 2.x was not affected. He kindly send me a patched version, which works fine for me. He also told me, that this **new version will be released soon**.


title: Is there an escape option for single quotes & apostrophes in QGIS?

text: I'm trying to use the Data Defined symbology in QGIS ver 2.2. Some of my data has single quotes/apostrophes. Since the string must use single quotes, is there an escape option I can use for the field calculator or for the data defined settings field? So, I'm trying to use something like this:               CASE     WHEN "NAME" = 'WILSON'S' THEN '0,75,224'     END      How do you build an expression around a string that may have a single quote or apostrophe? Thanks. QGIS Valmiera 2.2 Windows 7 - 32 bit


title: Display one of many rasters in Web Map

text: I would like to interactively add a georeferenced tif to my webmap. Use case is: I have a raster catalog of georeferenced tax maps. I would like to be able to select those tax maps (there are multiple) that intersect that tax lot, then choose the tax map of my choice of those selected to display in my web application (JavaScript API). Alternatively (and more simply), I would like to be able to have one of many rasters display based on the name field of the raster catalog. Any suggestions would be useful. My raster catalog is in a published service. I am open to using mosaic datasets as well (w/ Image Service).


title: Point in Polygon operation for 100G of data

text: I have lat/lon point locations for over a billion points, all based in the US. For each point I would like to generate a dataset that links them to their county and to their MSA according to US Census Shapefiles. My question is how to do this efficiently? I have access to a research computing server, postgis, shapely and statistical software.


title: Unable to load shapefile as feature class into geodatabase file, with only 3 of the fields

text: I'm trying to automate importing census boundary shapefiles to feature class in a Geodatabase File, but somehow could not get it to work with this code:               ### ======================== load boundary map ===============================          outLocation = "C:/python/RA506Final/RA506FinalPro.gdb"     outFeatureClass = "RamBoundaryMap"     inputTable = "tl_2010_27123_tract10.shp"          if arcpy.Exists(inputTable):         # List of field headers to convert         inFields = ['FID', 'Shape', 'GEOID10']         fm= arcpy.FieldMappings()         fieldMap= arcpy.FieldMap()         for field in inFields:             fieldMap.addInputField(inputTable, field)             fm.addFieldMap(fieldMap)         #creates a fieldmap with name of input field             vars()[field] = arcpy.FieldMap()             fm.addFieldMap(vars()[field])             #fm.AddFieldMap(vars()[field])              #convert it now         arcpy.FeatureClassToFeatureClass_conversion(inputTable, outLocation, outFeatureClass,field_mapping=fm)              # Confirm Ramsey/Hennepin Boundary map shapefiles imported into Geodatabase file         x = arcpy.ListFeatureClasses("*", "polygon")         for i in x:             print i     else:         print inputTable + " does not exist"          print arcpy.GetMessage      Here is the error message: > FieldMappings: Error in adding a field map to field mappings <function > GetMessage at 0x0C4E8330>


title: How do I find vector line bearing in QGIS or GRASS?

text: How can one extract the bearings from a set of lines in a vector file? I'm trying to find the bearings of several straight vector line segments (simplified contour lines, actually) and can't figure out how to do it. The QGIS measure tool can tell me individual line bearing if I use a N/S grid overlay to provide a 0º reference, but my mouse precision is probably questionable, and doing that for every line in our sample is going to take longer than I can devote to this project. Even if there is just a way to automatically pair line end-point nodes, get the XY data, and then let excel calculate the bearings for me, I would love to know. I really like using QGIS, but am also OK with GRASS (though shaky on text interface), SAGA, and ArcGIS (but only if I really have no other choice), so if any of the programs have a clever way to do this, let me know. -Nick


title: Where can I get good GIS data with height profile of Malaysia?

text: I want to create a 3-dimensional map of Malaysia with grass, but first I need the required data. Does somebody know a good source? EDIT: s a final result, I want to make a 3-dimensional map, that you can use in a flying plane for example. NASA SRTM has a precision from about 90m (I read that on wiki). I would be very happy with a precision from about 10m. I have a Budget from about 500$ (it's a hobby project of mine), but I would prefer free data. It doens't have to be the latest or most precise data.


title: How to show only some labels sharpmap

text: I am new to Sharpmap and GIS. Going through the sample code and building simple programs for better understanding. I have a water bodies layer. I want to display only a particular river. I used CustomTheme. How to show label for only this river. To display a river (named Ganga), I have written following code.               SharpMap.Data.Providers.ShapeFile shapeFileProvider = new SharpMap.Data.Providers.ShapeFile(@"C:\water.shp", false);     SharpMap.Layers.VectorLayer waterLayer = new SharpMap.Layers.VectorLayer("Water", shapeFileProvider);          SharpMap.Rendering.Thematics.CustomTheme waterTheme = new  SharpMap.Rendering.Thematics.CustomTheme(GetWaterStyle);     waterLayer.Theme = waterTheme;     mapBox1.Map.Layers.Add(waterLayer);      The GetWaterStyle method:               private SharpMap.Styles.VectorStyle GetWaterStyle(SharpMap.Data.FeatureDataRow row)     {         SharpMap.Styles.VectorStyle style = new SharpMap.Styles.VectorStyle();         if (row["NAME"].ToString().Equals("Ganga"))         {             style.Fill = Brushes.Aqua;             return style;         }         else             return null;     }      Now I want to show a label "Ganga" (the river name) as well. How to approach? Please help. Thanks and regards, Abhi


title: When digitizing Lines, vertices are automatically generalized when finished feature

text: I am creating lines(feature class) and when the feature is finished Arcmap 10 will automatically generalize/smoothens the line. Is there any reason why this is happening? is this a bug of ArcMap 10? Any answer will be much appreciated. Thank you very much.


title: How to keep polygon from becoming polyline when converting polygon feature class to .dwg

text: Is it possible to convert a polygon feature class from a geodatabase to a CAD .dwg while keep the features as a polygon in AutoCAD? Currently the polygons appear only as polylines. I am currently using ArcGIS 10.2 and AutoCAD Map 3D 2010. The Conversion Tool - Export to CAD appears to produce the necessary point, polyline and polygon data in the ArcMap Table of Contents after the tool has run, but only the polyline appears when I try to open the .dwg file into AutoCAD Map 3D. I have attempted to convert the polygon to a variety of CAD types including DWG_R2010 and DWG_2007. The DXF and DGN files do not appear to be recognized by AutoCAD when converted. It has been several years since I have touched AutoCAD, and any hints (is it even possible?) would be greatly appreciated. Perhaps I am just not going about it in the correct manner?


title: Split by attribute automatically on several layers in QGIS

text: I need to do a 2 stage split by attribute, using QGIS 2.0.1-Dufour. After the first split, I got 30 layers. Is there any way to do the second split automatically? I don't really like the idea to make the split manually :) It would be great to select all the opened layers, and make the split by the - common - attribute with one click. The result would go in a separate directory for each layer. Thank you very much for your help!


title: Refresh Database automatically

text: I'm a new user of QGIS. I have a problem trying to update layers from a database. I have a MySQL database in a server that I add as a layer (or more than one, depends on the project) with no problem. That database changes each five minutes, and I want to update the associate/s layer/s. Now, I push the refresh buttom, but I guess that this could be done automatically, without my interference, as a programmed routine or plugin. Any have tried it? Do you know how to achieve it? Thanks in advance! Laura


title: How to detect burned areas from Landsat or MODIS datasets using ILWIS?

text: I'm trying to find a simple tutorial on how to detect and map burned areas from/with ILWIS . I use Landsat 5/7 or MODIS for detection. Any idea?


title: Measure Displacement in ERDAS Imagine

text: How do I measure the displacement between two images in ERDAS Imagine. What tool do I use? I have looked in the ERDAS field guide.


title: Large shapefile to raster

text: I have a large shapefile (1 gb) and need to rasterise it. I have tried the following already. 1.) Import it into GRASS using `v.in.ogr` this failed with the error message: ERROR: G_realloc: unable to allocate 498240036 bytes at break_polygons.c:188 2.) My second idea was to use PostGIS. Import the shapefile, resample it at the x,y location of the grid and then export these points and create a grid from xyz. I successfully imported the shapefiles (polygons and points) but intersecting 1 million polygons with 300k points seems to be very slow. I used the following PostGIS satement, maybe there is room for improvement.               select polygons.land_id,grid.geom from grid,polygons where grid.geom && polygons.geom and within(grid.geom,polygons.geom)      3.) I did also try to use `simplify()` in PostGIS. But I lost to many small polygons (i.e. some areas that were covered only with small polygons became null). Any ideas would be greatly appreciated.


title: Locating Wlan accesspoints-how?

text: Having a Netbook I would like to map routers geographically-how? Starting from scratch. Thanks a lot for your suggestions!


title: Converting GPX track_points layer to CSV with lat/lng with ogr2ogr

text: What's the correct command to convert a GPX track_points layer to CSV **with** lat/lng coordinates using ogr2ogr?


title: UnicodeDecodeError: While installing gdal in a new Python Virtual Environment

text: I was trying to create a Python Virtual Environment for some initial GIS development. I installed `pip` and `virtualenv` and then tried the following commands from the Linux-Terminal to accomplish the task               $ cd Test/virtualenvs_for_test     $ virtualenv --no-site-packages postgis-cb-env     $ source postgis-cb-env/bin/activate     $ pip install gdal      * * * But found the following error message: > error: **command 'gcc' failed with exit status 1** > > Cleaning up... Command /home/gislinux/Test/virtualenvs/postgis-cb- > env/bin/python -c "import setuptools, tokenize; **file** > ='/home/gislinux/Test/virtualenvs/postgis-cb- > env/build/gdal/setup.py';exec(compile(getattr(tokenize, 'open', open)( > **file** ).read().replace('\r\n', '\n'), **file** , 'exec'))" install > --record /tmp/pip-h8eZSu-record/install-record.txt --single-version- > externally-managed --compile --install-headers > /home/gislinux/Test/virtualenvs/postgis-cb-env/include/site/python2.7 failed > with error code 1 in /home/gislinux/Test/virtualenvs/postgis-cb- > env/build/gdal Traceback (most recent call last): File > "/home/gislinux/Test/virtualenvs/postgis-cb-env/bin/pip", line 11, in > sys.exit(main()) File "/home/gislinux/Test/virtualenvs/postgis-cb- > env/local/lib/python2.7/site-packages/pip/ **init**.py", line 185, in main > return command.main(cmd_args) File "/home/gislinux/Test/virtualenvs/postgis- > cb-env/local/lib/python2.7/site-packages/pip/basecommand.py", line 161, in > main text = '\n'.join(complete_log) **UnicodeDecodeError: 'ascii' codec > can't decode byte 0xe2 in position 42: ordinal not in range(128)** I have no sure idea what this error is exactly. Although initial web-hunting suggested a need to change encoding from `ascii` to `utf8`, but how?


title: Is there a way to save images with none stretching in ArcGIS

text: I have some satellite images, when I imported these in ArcGIS and went to properties and then symbology tab, I saw stretch type was set to standard deviation. I want to save these images with no stretch as I import these images into ArcMap, to display the images without stretching , how can I do so?


title: Geoserver layer group problem!

text: I seem to be having a problem previewing a layer group in geoserver. I have 3 shapefiles that i converted from a dxf and i can preview each successfully however when i try to preview the layer group nothing appears apart from the open layers border and nav controls. Any ideas? Thanks in advance


title: How to Read a Shapefile using C#?

text: I am new to ArcGIS and trying a few things..   1. I'm trying to read shape file(which is supply as an input) using c# and   2. I want to generate a .pdf file which contains underlying metadata information (similar to the Identify tool) ? Thanks in advance.


title: Why don't some layers line up afer converting Mapinfo TAB to Shapefile?

text: I have 6 MapInfo layers / TAB files which I need to convert to ESRI ShapeFiles, preferably using QGIS. When a colleague opens these in MapInfo 10, they all look fine and line up correctly. When I open them in QGIS, however, some layers are at completely different X / Y coordinates, and are tiny when compared to other layers. When my colleague uses the universal translater in MapInfo 10 it creates all the shape files (it actually creates more, there are polyline, rectangle, region and ellipse files), but again when I open up the resulting shape files in QGIS (I click on "user defined projection" during opening) it has the same result. In MapInfo under Map -> Options -> Projection I can see the projection as "Longitude / Latitude". Any ideas how I can export all my MapInfo layers to ESRI shape files, and retain their shape / size and location? I can get access to MapInfo 10 if changes in the source TAB files are required. Thanks in advance, really stuck on this one! Dominik EDIT: I know the process to convert MapInfo to ShapeFiles, and have tried this both in QGIS and MapInfo 10. The problem remains, no matter which application I use to make the conversion, there are problematic layers which are not in the correct location, and are tiny compared to the source layers. If someone has some insight as to how to see the projection for each TAB file, and how to change it that would be the starting point, I guess. This is NOT a question on how to convert from one format to the other!


title: Method to find "if map data is available or not"

text: I am using Arcgis javascript api 3.2. I use the mapservice,                https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer.       I realized that map data is not available when the scale(slider) is zoomed to the fullest at certain places on the map. Is there a method to find out if map data is available or not? Could you please share the steps in doing so.


title: heat map / density map from dynamic points table in mapserver?

text: how can i create a heat map (some would call a density map) from a points layer in mapserver? (preferably able to choose interpolation method, NN, IDW, etc and color map with a transparent color) the points are stored in postgis table and the table is dynamic, meaning i need to create the heat map on the fly. thanks allot!


title: How to draw lines originating from a point for a given angle and length in attribute table of point layer

text: I have a set of points from which i want to draw lines , the point layer has information about the angle and length of the line to be drawn. This point layer is generated from a polygon layer as a result of Directional distribution ellipse tool of ArcGIS. Now i want to draw lines originating from the point layer which will divide the ellipse into 4 equal parts. any ideas on which tool to use or how do i do it?


title: problems setting up postgis and importing data

text: Using Postgres 9.2.3.1 on a new Mac. I installed Postgres from the dmg file that I downloaded from the Postgres site. I then installed PostGIS using stackbuilder. This all worked flawlessly so far as I can tell. I have tweaked the pg_hba.conf file to trust connections from local and localhost and created a database "db". Then I tried::               Mac-mini:pg postgres$ bin/createlang plpgsql db     createlang: language "plpgsql" is already installed in database "db"      so it would appear that postgres knows all about postgis. When I tried:               Mac-mini:pg postgres$ bin/shp2pgsql -c ~rful011/tmp/weed_boundaries weed_plot_boundaries | psql db     Shapefile type: Arc     Postgis type: MULTILINESTRING[2]     SET     SET     BEGIN     NOTICE:  CREATE TABLE will create implicit sequence "weed_plot_boundaries_gid_seq" for serial column "weed_plot_boundaries.gid"     CREATE TABLE     NOTICE:  ALTER TABLE / ADD PRIMARY KEY will create implicit index "weed_plot_boundaries_pkey" for table "weed_plot_boundaries"     ALTER TABLE     ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist     LINE 1: SELECT AddGeometryColumn('','weed_plot_boundaries','geom','0...                    ^      I must have missed something but what? Thanks, Russell


title: Is it possible pull multiple documents out of a geodatabase?

text: I am working on a sliverlight application in which a municipality will be able to view multiple assets for a water treatment plant within a GIS. Within this application, the user will be able to pull up linked documents such as operations and maintenance docs, record drawings, photos, vendor info, etc. Is there a way to export multiple documents at once (say into a zip file onto a local drive)? All documents are stored in a geodatabase, and I have researched the "Attachments" geoprocessing tools, but that is more for adding/removing multiple records, not extracting multiple records into a zip file.


title: Where can I get a trial version of ArcGIS Publisher?

text: I would like to save my MXD files as PMFs. I don't have a license of ArcGIS, only the 60-day trial. Is there a similar trial for the Publisher? I understand that's the only way to save PMF files (through the Publisher).


title: How to represent floors in a building

text: Is there any way to represent floors in a building using GIS? (Also, it's stupid that I have to use existing tags cos none suit my question so I had to choose one at random.)


title: How to delete Field From Shapefile by Arcpy

text: I am using the following code to delete a field from a shapefile but I am getting an error!               def main():         # Import arcpy module        import arcpy          # Local variables:        point_shp = "G:\\Temp\\All_Provinces.shp"          # Process: Add Field        #arcpy.AddField_management(point_shp, "NEWFIELD", "TEXT", "", "", "25", "", "NON_NULLABLE", "NON_REQUIRED", "")             arcpy.DeleteField_management(point_shp,NEWFIELD)     if __name__ == '__main__':         main()      can you please let me know what I am doing wrong?


title: How do i create a flow line map for my data?

text: I am trying to illustrate a passenger travel flow map from a road terminal to several destinations on a road layer. I have QGis and ArcGIS, but I don't know much on how to start the work. I already have the road terminal and all the various destinations datasets. I also have data on the number of passengers moving to all the various destinations. Can someone assist me with some guide?


title: Need of an imagery for Mexico

text: I am doing Rapid Watershed Assessment for my agency and I am responsible for the Southern New Mexico and northwest Texas. I have 2 NAIP added to my ArcMap already and using it with the "Add Image Server Connection" button. The NAIP I have here is New Mexico and Texas. Now I would like to know if there is one for Mexico that I'd like to add to my ArcMap. Is there one that can I get one using the add image service connection button? Thanks


title: Google's Direction API - Indoor Directions

text: I am using the Google's Android MAPs API V2 for indoor mapping in our app. It displays the indoor map correctly. Now, I am calling the Directions API to get the co-ordinates between two points to plot them on the indoor map for the way finding. However the Directions API always returns outside walking co- ordinates for two points selected within the indoor map. isIndoorEnabled is enabled too. Is Indoor directions API only available for Google's Maps app ?


title: why is map tiling computationally intensive?

text: I recently starting using `gdal2tiles.py` (and its parallel implementation) to generate tiled images from raster datasets for kml interpretation and I'm shocked by how computationally intensive the process is. For example, rendering a 3.8MB PNG image of the US (35000 x 15000 px - or 2500 dpi @ 14" x 6") with tile levels 0-7 takes over an hour on a modern quad-core system. Even then, the resolution isn't all that great! What is happening during the tiling process that is so computationally intensive?


title: How to add a legend for raster (GRID) layers in ESRI JS app

text: I have a grid in an MXD that I am publishing using ArcGIS Server 9.3.1. In the TOC I have made a legend with labels against each of my 5 classifications. I was expecting this to appear in my JS web app - but they never render. Does anyone know how to make a raster legend to work with the ESRI JS API. The code I am using is direct from the example pages - and it shows a legend for the basemap so I know the code is working to a certain extent:               //add the legend             dojo.connect(map,'onLayersAddResult',function(results){               var layerInfo = dojo.map(results, function(layer,index){                 return {layer:layer.layer,title:layer.layer.name};               });               if(layerInfo.length > 0){                 var legendDijit = new esri.dijit.Legend({                   map:map,                   layerInfos:layerInfo                 },"legendDiv");                 legendDijit.startup();               }             });      UPDATE - Begining to think that this is only available at v10.x.....


title: Load Leaflet map in canvas element

text: I am working with leaflet tiles loading tiles to a map div http://jsfiddle.net/tBCd5/ I am now playing with some html5 canvas distortions like below http://jsfiddle.net/UF4PK/ It would be cool to load base tiles to a canvas and then be able to use the above distortions on the base map. Is it possible to load my map tiles into a canvas element rather than a div? I have seen examples of creating canvas overlays but not an example where base tiles are loaded directly to the canvas element.


title: How print overviewMap using MapFishPrint?

text: How create overviewMap MapFishPrint? In my spefc json, i have tried insert "overview":true for example:               app/pdf/print.pdf?spec=     {"title": "Printing Demo","layout":"A4 portrait","srs":"EPSG:4326","units":"degrees","dpi":254,"outputFilename": "ggg", "pages": [         {             "center": [-48.634135544908 ,-22.501847553849 ],             "scale": 7000000,             "rotation": 0,             "comment":"",             "mapTitle":"",         }], "layers":[         {             "baseURL":"http://xxx.xxx.xx.xx.xx/geoserver/datageo/wms",             "opacity":1,             "singleTile":true,             "type":"WMS",             "layers":["AERODROMO_ANAC_2013","RodoanelTrechoConstruido"],             "format":"image/png",             "styles":[],             "overview":true         }        }       In may config have:               - !columns       # columns can have an absolute position. In that case, they need a width.       absoluteX: 28       absoluteY: 583       width: 578       config:         borderWidth: 0.3         cells:           - padding: 2       items:         - !map           width: 276           height: 462      But the overviewMap does not appear...


title: Clip_management to raster extent outputs larger raster than raster extent

text: Using Python and ArcMap 10.0 SP5, I'm trying to reproject and clip a large input raster to exactly the same grid as a smaller target raster (same extent, pixel sizes and projection). The problem is that the script produces an output raster that has 1 more row than the target raster (e.g. 1324 rows instead of 1323 -- columns match perfectly, pixels appear to be on the same grid). There is an extra row of pixels that overhangs the target image - everything else appears to be okay. The following is my code for projecting and clipping to a target raster. What's needed to ensure that the output extent is exactly the same as in target raster? I'm new to ArcMap so it could be something obvious.               import arcpy     import tempfile     import os     import shutil          sourceImage = "C:/data/WorldClim_alt/WorldClim_Alt.tif"     inputImage = "E:/subset/nCom_2000/nCom2000101to107.tif"     outputFilename = "E:/subset/test/nComTest_clipped.tif"          # Snap to source image and source image extent     arcpy.env.extent = sourceImage     arcpy.env.snapRaster = sourceImage          # Get cell size     cellx = arcpy.GetRasterProperties_management(sourceImage, "CELLSIZEX")          # Get temporary directory and filename     tempDir = tempfile.mkdtemp()     tempFilename = 'E:/subset/test/' + basename(outputFilename) +'_proj.tif'          # Reproject     arcpy.ProjectRaster_management(inputImage,tempFilename,sourceImage, "#", cellx)          #     # Clip to source file     #          # Get envelope for clipping     top = arcpy.GetRasterProperties_management(sourceImage, "TOP")     left = arcpy.GetRasterProperties_management(sourceImage, "LEFT")     bottom = arcpy.GetRasterProperties_management(sourceImage, "BOTTOM")     right = arcpy.GetRasterProperties_management(sourceImage, "RIGHT")     envelope = str(left) + ' ' + str(bottom) + ' '+ str(right) + ' ' + str(top)          # Clip      arcpy.Clip_management(tempFilename, envelope, outputFilename, sourceImage,"#","NONE")      **Edit in Response to Paul's comment on Clip help:** > From the help on `Clip`, it mentions > >> If the clip extent specified is not aligned with the input raster dataset, this tool makes sure that the proper alignment is used. This may cause the output to have a slightly different extent than specified in the tool. I think that the input to clip (the output from `ProjectRaster`) may be aligned with the `in_template_dataset`. The `snapRaster` environment variable was set to the `in_template_dataset` before `ProjectRaster` was called. Also, the output cell size of `ProjectRaster` was set to the `in_template_dataset` cell size. Visual inspection of the grids seems to indicate that they are aligned (except for an extra row and column in the `ProjectRaster` output). The only effect of `Clip` seems to be the removal of a single column.


title: SLD filter function usage for placing symbol at a particular point on the geometry (the_geom)

text: I am applying a symbol to a point on ellipse. For this i used the startPoint function, but my ellipse's coordinates start from the right side on the curve so it wrongly places the symbol at the right side. I found numPoints and pointN functions which I thought of using for this by getting no. of points on ellipse from numPoints and using pointN to place the symbol at 3/4* numPoints. This is just logic which I am not sure if could be applied. Please help if someone knows how to write this in SLD , please suggest or some other way of achieving this through other functions. Thanks


title: How can I determine which polygons contain other polygons?

text: I am using only free, open source tools. I have two shape files: a shape file of voting precincts and a shape file of census tracts. I want to know: which voting precincts fit inside which census tracts? What is the mapping vocabulary describing what I am trying to answer? I am trying to know which polygons from the precinct shape file are bounded by which polygons from the census shape file? Is that right? How do I answer that question programatically? What tools would I use? I have a vague sense that I could do this in postgis -- but don't really know. Also, is it possible in qgis?


title: openlayers code PanZoomBar is not a constructor and GetFeatureInfo error

text:
I’m having trouble with my openlayers code displaying postgis data in geoserver Why would this line (also PanZoomBar) map.addControl(new OpenLayers.Control.Scale($('scale'))); Produce this error in FF error consol Scale / PanZoomBar is not a constructor Also the GetFeatureInfo code is only returning one feature is this because the way I’ve constructed the layers Please note this is just a template once I’ve worked out all the glitches I’ll be creating a perl script to build this page on the fly each time it loads will be unique based upon how many layers/ VIEWPARAMS the user wishes to view at once. Thanx Holy               <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">                         <link rel="stylesheet" type="text/css" href="http://localhost:8080/geoserver/openlayers/theme/default/style.css"/>
OpenLayers PostGis Layers

location

             Click on the map to get feature info


title: Openlayers dynamic filtering/styling of geojson layer

text: I'm working with Openlayers and overlaying a geojson file representing species type over OSM. I have an HTML form that allows users to select a species type and I then need to filter the geojson layer based on the selection. I then want the map to redraw, displaying only the filtered features. I've tried CQL and messed around with styling but I'm not getting very far. Does anyone have any experience of doing this? Thanks in advance. James


title: Why does IActiveView.PartialRefresh not work?

text: I have a problem with the method IActiveView.PartialRefresh which is not working as it should be (or I'm using it not correctly). I'm writing an ArcGIS 10 Add-In using C# and want to clear a selection. Here comes the code (taken from http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#//0049000000q5000000):               ESRI.ArcGIS.Carto.IFeatureSelection featureSelection = featureLayer as ESRI.ArcGIS.Carto.IFeatureSelection;                                                activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, null, null);     featureSelection.Add(feature);     activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, null, null); // <-- this works          // perform some operation     ...          // remove the selection     featureSelection.Clear();     activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, null, null); // <-- does not work     // activeView.Refresh(); // <-- works      In the middle (where the ... are) I show a dialog and while this dialog is open the feature should be highlighted, so I select it. After closing the dialog it should no longer be selected, so I clear the selection (actually this bit is working, as I see when looking in the Selection Menu of ArcMap). But if I only do a PartialRefresh after the dialog has been closed, the selection is still visible. It only disappears if I make a full Refresh which doesn't look nice and is much slower. Actually I also have the same problem with PartialRefresh (apart from any selection issues) when I only want to redraw a single layer               ArcMap.Document.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGeography, featureLayer, null); <- doesn't work!!!      , but I'm not sure if these two problems are interlinked to each other. I hope someone can help me with this problem...


title: Does anyone use TRS (Township, Range, Section) data?

text: I recently had a request for TRS/PLSS data and I didn't know anything about it. I have researched it and it looks like it is purely mathematical (for the most part) based on geo-coordinates but I'm wondering if anyone would actually benefit from this information. For example, in our address verification service, we could provide the TRS data as well as the lat/lon data for each address. Im wondering what industries actually use the TRS/PLSS data scheme.


title: Is snapping working when changing from NAD83(decimal degree) to NAD83 Zone UTM 11N(meters)?

text: I've used the following ESRI snippet to Convert the display extents in Pixels (at the current map scale) and then return out the map units.               public System.Double ConvertPixelsToMapUnits(ESRI.ArcGIS.Carto.IActiveView activeView, System.Int32 pixelUnits)     {            if(activeView == null)       {         return -1;       }       //Get the ScreenDisplay       ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = activeView.ScreenDisplay;            //Get the DisplayTransformation        ESRI.ArcGIS.Display.IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;            //Get the device frame which will give us the number of pixels in the X direction       ESRI.ArcGIS.Display.tagRECT deviceRECT = displayTransformation.get_DeviceFrame();       System.Int32 pixelExtent = (deviceRECT.right - deviceRECT.left);            //Get the map extent of the currently visible area       ESRI.ArcGIS.Geometry.IEnvelope envelope = displayTransformation.VisibleBounds;       System.Double realWorldDisplayExtent = envelope.Width;            //Calculate the size of one pixel       if(pixelExtent == 0)       {         return -1;       }       System.Double sizeOfOnePixel = (realWorldDisplayExtent / pixelExtent);            //Multiply this by the input argument to get the result       return (pixelUnits * sizeOfOnePixel);     }      and then this to check if the snapping is within tolerance ITolerance which is in Pixels(Example 10 pixels) :               Double dblToleranceX;     dblToleranceX = snaptolerance(ITolerance);// call ESRI Snippet      dblProximity = dblToleranceX;     ouble dblToleranceY = dblToleranceX;     pInputpoint = gpMxDocument.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(intX, intY);     pRectangle.PutCoords(pInputpoint.X - dblToleranceX, pInputpoint.Y - dblToleranceY, pInputpoint.X + dblToleranceX, pInputpoint.Y + dblToleranceY);     pSpatialfilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;     IFeatureClass pFeatureClass = gpFeatureWorkSpace.OpenFeatureClass(strFeaturecls);     pSpatialfilter.Geometry = pRectangle;     pSpatialfilter.GeometryField = pFeatureClass.ShapeFieldName;     pFeatureCursor = pFeatureClass.Search(pSpatialfilter, false);     pFeature = pFeatureCursor.NextFeature();          if (pFeatureClass.FeatureCount(pSpatialfilter) < 1)      {         MessageBox.Show("Snapping is not within tolerance");                         return;                     }      The code work perfectly on NAD83(decimal degree) but does not work when Changing the data frame's coordinate system to NAD 83 UTM Zone 11N (Meters). Appreciate your help if you can figure out a way around this issue in order to make the code work for these kind of situation .


title: google geocoding service combined with openlayers - not working

text: I'm trying to combine google's geocoding service with openlayers. So when a user types an address and hits enter, map zooms there and renders a simple point. I edited my code according to the previous link. Does not work. I'm using a vector layer as basic (comes from Geoserver/PostGIS), a google's layer as second basic. I also added a simple vector layer to render the adrress's point. Do I have to bind google's layer and geocoding service somehow? Here is my code. Can you spot what's wrong? Thank you...               var geocoder;     function codeAddress()      {alert ('inside function');//works     var address = document.getElementById("address").value;              //tell to user, he did not type anything     if (address==""){document.getElementById("geocodMsgG").innerHTML = "You gave no input";}     else{          alert(address);//works          geocoder.geocode( { 'address': address}, function(results, status)        {alert(status);//not work             if (status == google.maps.GeocoderStatus.OK) {               alert('inside 2 function');//not work          alert(results[0].geometry.location);//not work          //create a geometry point          pointgc = new OpenLayers.Geometry.Point(results[0].geometry.location);               //create a feature point          feature_point_gc = new OpenLayers.Feature.Vector(pointgc);                 //lfl contains bounds that i set earlier in my code            if (lfl.containsLonLat(results[0].geometry.location)){            //render the point on the layer named geoco            geoco.addFeatures([feature_point_gc]);            geoco.refresh({force:true});            //zoom there            map.setCenter(results[0].geometry.location);              }           //if point is not inside the bounds (lfl) then dont add - print msg to user           else{document.getElementById("geocodMsgG").innerHTML = "out of bounds";}       }             //refers to the second if(geocoding's status)       else {document.getElementById("geocodMsgG").innerHTML = "error";}              });//geocode      }//major else (if address empty or not)     }//function


title: Create line from two points, two different feature class

text: I'm attempting to automate this in ModelBuilder in ArcMap. I have one feature class with one point (Start) and one feature class with multiple points (Destinations). I want to iterate through the Destinations so I create multiple lines from the Start point to each Destination point. I was looking at possibly extracting the XY coordinates and then using "XY To Line", but that wants the xy to be in one table (or feature), not two. Is there a simpler way of connecting two points?


title: How to get selected Features/Geometry from Windows Clipboard in ArcMap/C#?

text: I have a custom tool in ArcMap 10. The Workflow is the following:   * Select Features in a Layer with the "Select Features by Rectangle" tool   * Copy the Features to the Clipboard using STRG+c   * Select the custom tool   * The tool checks the Clipboard. In case of finding a copied Feature in the Clipboard the Workflow changes  I found an example on http://kiwigis.blogspot.com/2011/07/how-to-create-drop- target-for-layers-in.html where similar stuff is implemented using drag and drop. Following and adapting the example and trying to copy layers in the TOC in ArcMap to the clipboard using STRG+c is working for me:               //Get the Data from the Clipboard     IDataObject clipBoardDataObject = Clipboard.GetDataObject();          //The format of the object data is "ESRI Layers" (this can be checked using clipBoardDataObject.GetFormats();)      MemoryStream geomStream = clipBoardDataObject.GetData("ESRI Layers") as MemoryStream;                     byte[] bytes = geomStream.ToArray();          IMemoryBlobStreamVariant memoryBlobStreamVariant = new MemoryBlobStreamClass();     memoryBlobStreamVariant.ImportFromVariant(bytes);     IMemoryBlobStream2 memoryBlobStream = memoryBlobStreamVariant as IMemoryBlobStream2;     IStream stream = memoryBlobStream as IStream;     IObjectStream objectStream = new ObjectStreamClass();     objectStream.Stream = stream;          byte pv;     uint cb = sizeof(int);     uint pcbRead;          objectStream.RemoteRead(out pv, cb, out pcbRead);     int count = Convert.ToInt32(pv);          //Guid you find in the ESRI ILayer Interface     Guid iLayerGuid = new Guid("34C20002-4D3C-11D0-92D8-00805F7C28B0");          for (int i = 0; i < count; i++)     {       object o = objectStream.LoadObject(ref iLayerGuid, null);       ILayer layer = o as ILayer;           }      I can cast the object to ILayer, this is working! Trying to achieve the same with copied Features is not working for me, at least i dont really know what to do.               //Get the Data from the Clipboard     IDataObject clipBoardDataObject = Clipboard.GetDataObject();          MemoryStream geomStream = clipBoardDataObject.GetData("ESRI Geometry List") as MemoryStream;                     byte[] bytes = geomStream.ToArray();          IMemoryBlobStreamVariant memoryBlobStreamVariant = new MemoryBlobStreamClass();     memoryBlobStreamVariant.ImportFromVariant(bytes);     IMemoryBlobStream2 memoryBlobStream = memoryBlobStreamVariant as IMemoryBlobStream2;     IStream stream = memoryBlobStream as IStream;     IObjectStream objectStream = new ObjectStreamClass();     objectStream.Stream = stream;          byte pv;     uint cb = sizeof(int);     uint pcbRead;          objectStream.RemoteRead(out pv, cb, out pcbRead);     int count = Convert.ToInt32(pv);          //What to do next? What object do i need to unpack to?      You see here that the object format is "ESRI Geometry List" (of which i cant find any infos online). So how can i "unpack" this object? IFeatureSelection, IFeature, IGeometry wont work. Any Ideas? I would be grateful for any leads. Thx!


title: Query against FeatureClass for a layer

text: This is some code that I am using to try to get a working IFeatureLayer from a subset of data from a database table. The fourth line containing "fcWkspc.OpenFeatureQuery" throws and error complaining about an "The operation was attempted on an empty geometry.". Should also mention that this code was lifted directly from ESRI documentation example code at the following link: http://resources.esri.com/help/9.3/arcgisdesktop/arcobjects/esrigeodatabase/IQueryDef_Example.htm Can anyone tell me what I might be doing wrong? Thanks,               // m_wkSpace2 is of type IWorkspace2     IFeatureWorkspace fcWkspc = (IFeatureWorkspace)m_wkSpace2;          IQueryDef qDef = fcWkspc.CreateQueryDef();          Def.Tables = "dbo.TABLE_NAME";     qDef.WhereClause = "Name1 = '" + var1 + "' AND Name2 = '" + var2 + "'";          IFeatureDataset featureDataset = fcWkspc.OpenFeatureQuery("MyJoin", qDef);     IFeatureClassContainer featureClassContainer = (IFeatureClassContainer)featureDataset;     IFeatureClass featureClass = featureClassContainer.get_ClassByName("MyJoin");          IFeatureLayer layerForMap = (IFeatureLayer)featureClass;


title: Issues when using pgrouting's shortest path

text: I am having some issues when using pgrouting's shortest path. I am following this tutorial available at http://workshop.pgrouting.org/chapters/shortest_path.html#dijkstra. This query returns a list of edges. Then I get the start and end points of each edge. Then I just join these points to show the shortest path. Just to let you know that I don't necessarily get the ending point of the first edge in the list to be the starting point of the second edge in the list returned. So I join, the ending and starting points of the two edges as well. However, sometimes the points are repeated and I have lines crossing over each other and the path looks a mess. What could be the reason behind this?


title: Cannot read CSV-file from ArcMap

text: when trying to import a CSV from my hardtrive to my map-document using "add Data"-dialog I get the error **"Selected element could not be added to map. Generic error"** (or something similar, I am using german version, so actual message may vary). When doing the same via TableToTable_conversion-tool in GP I get the error **000732: Dataset 'myFile.csv' not existing or not supported**. I also tried to create a file-GDB and import the table there, result is the same. Here is my CSV-structure:               SID,UUID     1,DENW22AL60001dVy     1,DENW22AL30000FLl      I also tried with semicolons rather then single comma as delimiter, but can´t get it working. Also notice that I only want a standalone-table rather then a featureclass. Therefor no coords are given within the file. EDIT: As of Error: 000732: Dataset does not exist or is not supported this error may occur if the file-name itself is incorrect. So I did the proposed corrections for filenames (use forward-slashes, no spaces) but with no avail. Also since I chosed the file by the file-dialog rather then by typing its name I did not even expect this would do the trick (just for completeness). I´m using ArcGIS 10.0 for Windows 7 (64Bit)


title: If then statement in model builder for several features after a clip

text: ![enter image description here](http://i.stack.imgur.com/uZbUI.png)I'm merging several databases and clipping out certain areas of interest. Sometimes there are no points that are clipped because there is nothing in the clip. I have each of the features output to a geodatabase. Is there a way for it not to export the feature class to the geodatabase if the file is empty?


title: Are there limitations with gdal virtual raster option?

text: I am new at Postgis and I'm groping... I pretend to rasterize some weighty vector layers in Postgis, do some map algebra and reclassification with them and visualize intermediate as well as final rasters in Qgis. I don’t want to have physical raster files (like geotiff or jpg) outside of DB, except at the end of the work flow. Has anyone experienced a similar case? I am thinking about using gdal virtual raster. Can I use Postgis raster spatial functions with this option, or just gdal functions? Has anyone experienced problems in Qgis viewing with gdal virtual raster from postgis db? Thank you.


title: Display county for random points

text: I have a county level shape file for South Carolina. I generated random points using Vector->Research Tools-> Random Points. The Attribute tables for the random points layer has only the ID included. I would like to export the latitude, longitude and corresponding county for all the generated points (there is a COUNTYFP10 variable in the attributes table of the original state layer). I managed to export the latitude and longitude using "Save as" , choosing "Comma Separated Value". Then, in the "Layers" part of the "OGR creation options" I typed "GEOMETRY=AS_XY". I am not sure how to add the corresponding county to the saved CSV file. Any help would be greatly appreciated. Thank you, Georgiana


title: Why might layer appear not to turn off in ArcMap?

text: I was called over to try and fix an issue with a coworkers ArcGIS software (ArcView 9.3.1 SP1, Windows 7) and couldn't fix it. When a certain layer is turned off in the table of contents it never turns off in the map. Refreshing the map doesn't work, jumping between data and layout view doesn't work either. Saving the map with the layer off, closing and then re-opening it doesn't make a difference. I removed and re-added the data, copied the shapefile to a different location, removed joins, etc. still the same problem. Has anyone had this issue and know of a fix? I don't have any control over the version and service pack they're using (they work for the state and I work at the county level, so we're on different systems). Other people in the office using the same data don't have the problem. I appreciate the help in advance.


title: Way to add 1km x 1km UTM grid onto online map service

text: Is there any exist services or way to add 1km x 1km UTM grid onto online map service such as Google Maps or OpenStreetMap?


title: How to create an Azimuthal Equidistant map for long-distance radio?

text: I'm trying to produce world or continent scale maps that use the Azimuthal Equidistant projection. I'm a radio amateur, and to an antenna, the whole world looks like spread out on a flat plane, with the antenna at the centre and the whole world at an (r, Θ) from that point. There is software to do this - AZ_PROJ \- but it uses its own database, and is entirely written in PostScript, so it's not exactly hackable with the usual geo tools. OGR has the 'aeqd' projection, which can be invoked like this:               ogr2ogr -t_srs "+proj=aeqd  +lat_0=43.7308 +lon_0=-79.2647" out.shp in.shp      for a centre at 43.7308 N, 79.2647 W. Unfortunately, this wraps everything beyond the current hemisphere "behind" the map. It also projects the coordinates into cartesian metres, where I'd prefer polar. Any suggestions for alternatives, please?


title: Different results from Spatial Analyst GP Tool and Python

text: I am writing a script that has to subtract one raster from another. The rasters I have represent 1) Parking Lots (1 = Parking Lot, 0 = not Parking Lot) and 2) a buffer around the Parking Lot (1 = Parking Lot or 20ft Buffer around Parking Lot, 0 = not Parking Lot or Buffer). What I need is a single raster representing just the buffer area without the actual Parking Lot, basically a ring around each Parking Lot (1 = ring around Parking Lot, 0 = not ring around Parking Lot). When I perform this task manually in ArcMap (v10.1) using the Minus tool in the Spatial Analyst Extension, it works appropriately. When I try the same thing in Python (testing from the command window in ArcMap), the subtraction does not occur and the buffer with Parking Lot raster is returned unchanged. The same result occurs whether I use `arcpy.sa.Minus(buffer, parkinglots)` or `buffer - parkinglots`. This particular workflow is part of a larger script so doing it manually in ArcMap is just not an option. Any suggestons? **UPDATE:** I tested this a slightly different way. I create a new model that performed the `Minus` operation the same way I did it manually in ArcMap. It worked. According to the documentation, there are 2 ways to accomplish a subtraction. You can either use the `-` operator or the equivalent `Minus` GP tool within the `Math` toolset of the `Spatial Analyst Tools` toolbox (`arcpy.sa.Minus`). So, I tried something else, I exported that model to script. What I found was fairly interesting. In the script that was exported from my model, the Python code was `arcpy.gp.Minus_sa`. I haven't found any documentation about this, but it worked.


title: Splitting polylines according to coincident (or touching) points with other polylines?

text: The question is unfortunately not answered yet. I'm looking for a method or python script to split polylines according to **coincident points** with other polylines, but **avoid splitting the ones which just intersect** (without a sharing point like the planarize or intersection tool in ArcGIS would produce). The "touches (second_geometry)" function might be a way but I have at the moment no clue how to solve it. The image illustrates my aim. The _yellow line_ hasn't changed, because it has no coincicent points with the _red line_. While the _red line_ as well as the _dark blue line_ is split into 4 different polylines, because both share a coincident point. ![enter image description here](http://i.stack.imgur.com/ekXvN.png) Thanks for any answer in advance!


title: What is the best open-source web mapping solution for handling 3000+ polygons

text: I have 3000+ polygons which represent watersheds (18 main watersheds). At the first stage all the 18 watershed will show and their list will be placed in the left sidebar. If any user mouse over the map only popup will show the watershed name. If any user click on any major watershed on map it will open pop up window with detail. At a larger zoom, all watersheds will show. I'm asking for your advice which open-source solution is the best for handling that large amount of polygons? I was thinking about Leaflet or OpenLayers but I don't know can they handle large amounts of polygons. Also, do you suggest do store them in mysql or? Thank you.


title: Is there a way to count the number of vertices of polygons and lines?

text: In QGIS, is there a way to count the number of vertices of polygons and lines in a given layer, and create a separate field containing these values? I suspect it's a function somewhere in the field calculator but I can't find it.


title: Should I implement TSP or Dijkstra?

text: I was asked to create a shortest path algorithm in java to use with OSRM. I want to create a route between some points (normally more than 3). The first one, will be the starting point and also the last point (close route). Other point have to be the last before the last one and rest of points have to be ordered taking the shortest path between the first and the last before the first again. I think this is no possible to do with TSP because TSP doesn't have a known first point and also it is not possible to determine if a point have to be visited in a concrete position. Is it possible? So I think I should use dijkstra for that and create a route between the starting point and the last before the starting point, after that, add the starting point to the route. Talking about Dijkstra with some friends, I was adviced about Dijkstra doesn't visit all points always, maybe some of them are left out the route.. .Is this true? Also I hear about Dijkstra doesn't create a route, it creates a tree... also, is it true? Which algotirhm should I implement? Thanks and sorry for brick ;)


title: Mapnik WKT rendering problem

text: I have a rendering problem with the mapnik GEOS plugin. I try to input a polygon WKT definition containing floats with around 16 digits. Mapnik however does not render these. When I decrease the number of digits to 5 it does work. Does anybody have a solution?               mmap.background = mapnik.Color('steelblue')     s = mapnik.Style()      r = mapnik.Rule()      polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#f2eff9'))     r.symbols.append(polygon_symbolizer)     [ more style here ]     mmap.append_style('My Style',s)           wkt_geom = 'POLYGON ((5.12345 51.6, 6.8 52.6, 7.8 51.6, 5.12345 51.6))' #DOES WORK     wkt_geom = 'POLYGON ((5.123456 51.6, 6.8 52.6, 7.8 51.6, 5.123456 51.6))' #DOES NOT WORK          layer = mapnik.Layer('wkt_layer_bla','+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs')     layer.datasource =  mapnik.Geos(wkt=wkt_geom)     layer.styles.append('My Style')     mmap.layers.append(layer)     mmap.zoom_all()     mapnik.render_to_file(mmap,'/foo/bar/bla.png', 'png')


title: Geocoder of google maps on openlayer map

text:
I'm trying to integrate the geocoder of google maps into my ol-map. There are a lot of similar questions asked. I tried out this: http://dev.openlayers.org/sandbox/topp/almanac/examples/geocoder.html example but it doesn't work. looking through that: How to use Google Maps geocoding with Openlayers? also didn't help me. (As I'm new here I cannot comment on peoples answers so I thought it might be ok to open a new question - tell me if I'm wrong) See code below for my attempt. Do you have any tips or know any working example??                    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />     <title>Google styled in OL mit kml</title>     <link rel="stylesheet" href="OpenLayers-2.12/theme/default/style.css" type="text/css">     <link rel="stylesheet" href="OpenLayers-2.12/theme/default/google.css" type="text/css">
OSM + Google Maps + KML Reprojection
      <input id="address" type="textbox" value="Sydney, NSW">      <input type="button" value="Geocode" onclick="search()">


title: Opening shapefiles with geotools

text: Trying to make a simple class that will read shapefiles and eventually be able to pull out the coordinates/attributes for each feature (working with points here). But can't quite figure out geotools.               File file = new File(filepath);     if(!file.exists() || !filepath.endsWith(".shp")) {         throw new Exception("Invalid shapefile filepath: " + filepath);     }     ShapefileDataStore dataStore = new ShapefileDataStore(file.toURL());     ContentFeatureSource featureSource = dataStore.getFeatureSource();     ContentFeatureCollection featureCollection = featureSource.getFeatures();          SimpleFeatureIterator iterator = featureCollection.features();     while(iterator.hasNext()) {         SimpleFeature feature = iterator.next();         List attributes = feature.getAttributes();         Collection properties = feature.getProperties();     }      Fails at `dataStore.getFeatureSource()` with the error: `java.lang.NoSuchMethodError: org.geotools.data.shapefile.ShapefileDataStore.getFeatureSource() Lorg/geotools/data/store/ContentFeatureSource;` The imports for every geotools class mentioned above are there. As some of you know, geotools comes with a crapton of jars. As I'm unfamiliar with Maven and this is an existing project that I don't know how to convert over, I'm trying to do the library imports the old fashioned way. Right now my imports are: `gt-api-11.1.jar`, `gt-data-11.1.jar`, `gt-main-11.1.jar`, `gt- opengis-11.1.jar`, and `gt-shapefile-11.1.jar`. I figure it may be a missing jar that could be the problem but I'm not sure which as every imported class is in one of these jars.


title: Problem with field calculator

text: New to QGIS and suspect I'm missing something basic but.... When I input the formula to create a new field, all the calculations are correct but when I press the "Save Edits" button, all the figures in the new column revert to "Null"


title: How to Set the width of a dropdown box?

text: I have created my own script in ArcGIS 10.1 and for one of my parameters I have a drop down list that the user has to select from. The longest entry in that list is only 6 characters, but the width of the drop down box is the total width of the window. Is there a way to set the drop down box width to something smaller.


title: How do I stack feature class annotation?

text: I have map annotation that is stored in a database (as a feature class), but that is not linked to the features. For the life of me, I can't figure out how to edit it so that I can take a label of say 3 words and stack it so that one word is on top of the other, like so: From:   How To Stack To:   How   To   Stack Any help would be greatly appreciated!


title: Cutting OSM roads in sections using PostGis

text: I'm a newby in postgis. I define a section by a part of a road that has no intersection. See image bellow: ![enter image description here](http://i.stack.imgur.com/imWJU.png) I'd like to extract all sections from all the roads of my postgis database. This gis database was created with osm2pgsl:                              List of relations      Schema |        Name        | Type  |  Owner     --------+--------------------+-------+----------      public | geography_columns  | view  | postgres      public | geometry_columns   | table | gis      public | planet_osm_line    | table | postgres      public | planet_osm_nodes   | table | postgres      public | planet_osm_point   | table | postgres      public | planet_osm_polygon | table | postgres      public | planet_osm_rels    | table | postgres      public | planet_osm_roads   | table | postgres      public | planet_osm_ways    | table | postgres      public | spatial_ref_sys    | table | gis     (10 rows)      I've read an interesting post about detecting the intersections but I cant' figure out how to create a table with sections. Any help appreciated.


title: cutting and extending one polygon layer in function to the other layer boundaries

text: I'd like a methodology that can edit my vector layer, extending and cutting the features in the polygon layer "cap_agricola" within the boundary layer "sp_edr", I would like to solve this problem in ArcGIS, but I can do without restrictions in GRASS or QGIS. ![Problem - Cut and Extending cap_agricola layer to the sp_edr boudaries layer](http://i.stack.imgur.com/XJxum.jpg)


title: Pyshp Editor error

text: Following my predecessors (here and here) I am trying to append new data to shapefile using Pyshp. The specifics of my problem is the need to save the file with the same name (without creating a new file). This simple example based on Pyshp documentation failes on the attempt of the edited file saving:               #!/usr/bin/python3.2          import shapefile          #create file     w = shapefile.Writer(shapefile.POINT)     w.field("Lon", "N", 20, 10)     w.field("Lat", "N", 20, 10)     w.field("Z", "N", 20, 10)     w.point(0.0, 0.0)     w.record(0.0, 0.0, 0.0)     w.save("test")          #append new record     e = shapefile.Editor(shapefile="test")            <-- Existing file is opened for editing     e.point(1.0, 1.0)     e.record(1.0, 1.0, 1.0)     e.save("test")                                    <-- Script falls out      * * * Here is my traceback:               ---------------------------------------------------------------------------     TypeError                                 Traceback (most recent call last)     <ipython-input-29-e3063921c803> in <module>()     ----> 1 e.save("out2")          /usr/local/lib/python3.2/dist-packages/pyshp-1.2.0-py3.2.egg/shapefile.py in save(self, target, shp, shx, dbf)        1026                 target = temp.name        1027                 generated = True     -> 1028             self.saveShp(target)        1029             self.shp.close()        1030             self.saveShx(target)          /usr/local/lib/python3.2/dist-packages/pyshp-1.2.0-py3.2.egg/shapefile.py in saveShp(self, target)         983             self.shapeType = self._shapes[0].shapeType         984         self.shp = self.__getFileObj(target)     --> 985         self.__shapefileHeader(self.shp, headerType='shp')         986         self.__shpRecords()         987           /usr/local/lib/python3.2/dist-packages/pyshp-1.2.0-py3.2.egg/shapefile.py in __shapefileHeader(self, fileObj, headerType)         713         z = self.zbox()         714         # Measure     --> 715         m = self.mbox()         716         try:         717             f.write(pack("<4d", z[0], z[1], m[0], m[1]))          /usr/local/lib/python3.2/dist-packages/pyshp-1.2.0-py3.2.egg/shapefile.py in mbox(self)         685     def mbox(self):         686         """Returns the current m extremes for the shapefile."""     --> 687         return self.__mbox(self._shapes)         688          689     def __shapefileHeader(self, fileObj, headerType='shp'):          /usr/local/lib/python3.2/dist-packages/pyshp-1.2.0-py3.2.egg/shapefile.py in __mbox(self, shapes, shapeTypes)         671             except IndexError:         672                 pass     --> 673         return [min(m), max(m)]         674          675     def bbox(self):          TypeError: unorderable types: NoneType() < int()      Thank you


title: Extract characters inside parentheses in QGIS

text: I have a shapefile of the buildings on my college campus. It has a "name" field that includes the name of the building and then in parentheses the building's campus address (i.e. Health and Counseling Services (71)). I want to create a new field, "number", that is only the number within the parentheses (i.e. '71'). I am working in QGIS with regular expressions but I can't quite get it right. Any help would be greatly appreciated!


title: Geoportal Server 1.1.1 search widget for Flex 2.5

text: Does anyone one know where I can get the above mentioned search widget? I have been tasked with installing this to an existing geopoprtal server and Flex viewer does anyone know how to do this. i have found some instructions but they seem to require reloading the viewer. Any Ideas or advice would be appreciated. Thanks


title: Export data from .mxd file to ms sql database

text: I am using ArcGIS product (v.10.1). I have also tempalte of Maps from ArcGIS. But i need to get data from maps (.mxd files) to store it into MS SQL database. Explanation: I have application that has map in it's window. When i click on some city on map i get name of the city and it's coordinates and save them into DB. Problem - i need to have list of cities in database before i click on map to make some custom autocomplete for textBox. * * * I read some articles about getting data using Python, but i don't know this language. Any advices how to export data from mxd files? * * * Note: I am using C# to work with data.


title: Import raster to postgis and add a foreign key constraint

text: I use raster2pgsql to import raster files into a table with this command line:               C:\Program Files\PostgreSQL\9.3\bin>raster2pgsql -s 4326 -M -F -a c:\temp\mun.tif public.basqui_raster | psql -U postgres -d basqui -h localhost -p 5433      However I would like to add a foreign key constraint to each row during the import. Is it possible to do that when importing raster with raster2pgsql? What could be the best way to do it? Edit: I would prefer not use the filename (-F) in a second step to create the foreign key. I would prefer to do it during the import.


title: Having 1 base layers hide when the other one is shown and visa versa

text: I have an open layers map that is using Mapquest as a base layers as well as my own base layer. I've seem to be having an issue with getting the one base layer to show while the other one hides. The one only needs to show from a certain zoom level to another and it goes the same with the other one. My code is as follows:                function initMap(smX, smY, zoom) {                 // OpenLayers                 OpenLayers.IMAGE_RELOAD_ATTEMPTS = 5;                      OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {                     defaultHandlerOptions: {                         'single': true,                         'double': false,                         'pixelTolerance': 20,                         'stopSingle': false,                         'stopDouble': false                     },                     initialize: function (options) {                         this.handlerOptions = OpenLayers.Util.extend(                                 {}, this.defaultHandlerOptions                             );                         OpenLayers.Control.prototype.initialize.apply(                                 this, arguments                             );                         this.handler = new OpenLayers.Handler.Click(                                 this, {                                     'click': this.onClick,                                     'dblclick': this.onDblclick                                 }, this.handlerOptions                             );                     },                          onClick: function (evt) {                     },                          onDblclick: function (evt) {                     }                      });                      var mapOptions = {                     numZoomLevels: 20,                     controls: [],                     projection: new OpenLayers.Projection("EPSG:900913"),                     units: 'm',                     restrictedExtent: new OpenLayers.Bounds(bbox = 1821416, -4151368, 3671526, -2516543)                 };                 map = new OpenLayers.Map('testMap', mapOptions);                mqBase = new OpenLayers.Layer.XYZ(                     "Map Quest - Streets",                     [                     "http://otile1.mqcdn.com/tiles/1.0.0/osm/${z}/${x}/${y}.png",                     "http://otile2.mqcdn.com/tiles/1.0.0/osm/${z}/${x}/${y}.png",                     "http://otile3.mqcdn.com/tiles/1.0.0/osm/${z}/${x}/${y}.png",                     "http://otile4.mqcdn.com/tiles/1.0.0/osm/${z}/${x}/${y}.png"                     ],                 {                     attribution: "Data, imagery and map information provided by MapQuest, Open Street Map and contributors, CC-BY-SA  ",                     transitionEffect: "resize"                 }                 );                mqBase.isBaseLayer = true;                //mqBase.setOpacity(1);                //mqBase.setVisibility(true);                //mqBase.displayOutsideMaxExtent = 0;                map.addLayers([mqBase]);                      cmaBase = new OpenLayers.Layer.WMS(                     "CMA Info - Streets",                     "http://192.168.1.40:8080/geoserver/CMAInfo/wms",                     //"http://gs.cmatile.co.za/geoserver/Cmainfo/gwc/service/wms",                     {                         layers: "Base", transparent: true, format: "image/png", buffer: 0,// isBaseLayer: true,                         //map: '/projects/ms4w/apps/tutorial/htdocs/base.map'}                         tiled: true                     }                      );                 cmaBase.isBaseLayer = true;                 cmaBase.setOpacity(1);                 cmaBase.setVisibility(false);                // cmaBase.displayOutsideMaxExtent = 0;                 map.addLayers([cmaBase]);                 //&& prevzoom < 15                 map.events.register("zoomend", this, function (e) {                     var curzoom = e.object.zoom;                     if (typeof prevzoom != "undefined") {                         if (curzoom > prevzoom && curzoom >= 15 ) {                             map.setBaseLayer(cmaBase);                             cmaBase.setVisibility(true);                             mqBase.setVisibility(false);                             }                     }                     //&& prevzoom >= 15                         else if (curzoom < prevzoom && curzoom < 15 ) {                             map.setBaseLayer(mqBase);                             mqBase.setVisibility(true);                             cmaBase.setVisibility(false);                       }                     prevzoom = curzoom;                 })                 map.setOptions(zoomLevel = 19);      I have tried several methods some of which it can be seen as comments as I was trying each one out. For some reason the setVisibilty does not seem to be working I have even used and If statement to check if it was a baselayer to try and set it that way and that also did not work. If someone can help to see if Im missing something or if Ive done something wrong it will be alot of help. I am using openlayers 2.13.1 and geoserver is my basemap for one.


title: ArcMap - WMS service layers not drawing

text: I'm using ArcObjects .NET to load WMS 1.3.0 service layers from our ArcGIS server into ArcMap 10.2.2. The process of connecting to the server and adding layers works fine, however none of the layers will draw - I simply have a blank canvas. I have tried two methods to add WMS layers to the map:   1. Using 'Add Data' button > GIS Server > Add WMS Server   2. With ArcObjects .NET (see code below)                  'setup connection properties         Dim pMap As IMap = pMxDoc.ActiveView.FocusMap         Dim pWMSLayer As IWMSGroupLayer = New WMSMapLayer         Dim pConnectionName As IWMSConnectionName = New WMSConnectionNameClass         Dim pPropSet As IPropertySet = New PropertySetClass         pPropSet.SetProperty("URL", txtURL.Text)         pConnectionName.ConnectionProperties = pPropSet         Dim pDataLayer As IDataLayer = CType(pWMSLayer, IDataLayer)              'connect to server         Try             pDataLayer.Connect(CType(pConnectionName, IName))         Catch ex As Exception             MsgBox("Error connecting to server: " & ex.Message)         End Try              'get top level group layer from server         Dim pServiceDesc As IWMSServiceDescription = pWMSLayer.WMSServiceDescription         Dim pLayerDesc As IWMSLayerDescription = pServiceDesc.LayerDescription(0)         _WMSGroupLayer = pWMSLayer.CreateWMSGroupLayers(pLayerDesc)              'add layers to map         For i As Integer = 0 To _WMSGroupLayer.Count - 1             Dim pLayer As ILayer = _WMSGroupLayer.Layer(i)             pMap.AddLayer(pLayer)         Next      Both methods add the layers to the map and they are visible in the table of contents with no errors or warnings. I can view layer source (seems correct), extent (also seems correct), some symbology properties, and I can also zoom to the extent of the layer (which changes map zoom), however nothing at all is visible on the map canvas. Scale range isn't the issue as I've tried zooming in far beyond the specified range. I've come across this page relating to WMS layers not drawing in ArcMap 9, but it isn't offering me any working solutions for 10.2: WMS layers won't draw in 9.3 I have tried changing coordinate system of data frame (and layers), no to effect. They are all in WGS1984 which is supported by the WMS service. At this point I don't know what else I can try besides contacting ESRI for assistance. Does anyone have any ideas or similar issues with WMS layers? Thanks. EDIT: also looks like this person had the same problem but no solution was posted.


title: Can I dewarp an image using GDAL?

text: I have a GeoTiff image mymap.tif https://www.dropbox.com/s/51chrzv9vp9olg8/mymap.tif ` $ gdalinfo mymap.tif Driver: GTiff/GeoTIFF Files: mymap.tif Size is 1085, 281 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4326"]] Origin = (-122.509475852140255,37.769141730303630) Pixel Size = (0.000010000000000,-0.000010000000000) Metadata: AREA_OR_POINT=Area Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=PIXEL Corner Coordinates: Upper Left (-122.5094759, 37.7691417) (122d30'34.11"W, 37d46' 8.91"N) Lower Left (-122.5094759, 37.7663317) (122d30'34.11"W, 37d45'58.79"N) Upper Right (-122.4986259, 37.7691417) (122d29'55.05"W, 37d46' 8.91"N) Lower Right (-122.4986259, 37.7663317) (122d29'55.05"W, 37d45'58.79"N) Center (-122.5040509, 37.7677367) (122d30'14.58"W, 37d46' 3.85"N) Band 1 Block=1085x32 Type=Byte, ColorInterp=Red Mask Flags: PER_DATASET ALPHA Band 2 Block=1085x32 Type=Byte, ColorInterp=Green Mask Flags: PER_DATASET ALPHA Band 3 Block=1085x32 Type=Byte, ColorInterp=Blue Mask Flags: PER_DATASET ALPHA Band 4 Block=1085x32 Type=Byte, ColorInterp=Alpha ` If I fetch the bounds of this image from Google Map I get this ground truth image https://www.dropbox.com/s/mrb3d1ho33mw9bw/truth.tif. The two images don't quite line up. Is it possible to dewarp my image so that when overlaid on the Google Maps image they line up correctly?


title: feature position on non-Geo-referenced Image

text: I use OpenLayers to show Map , when im creating the map, i want to show map image with an ImageLayer and VectorLayers, the image is non Geo-referenced. **The use case** : i receiving a map image url from the server and entities layers , the entities positions on the layers represented by pixels on the image that i'm holding. from the image url i'm creating a ImageLayer and then creating VectorLayer and want to put the features on the VectorLayer with the pixels positions. I tried to create the map in EPSG:4326 projection and then to convert the pixels point with getLonLatFromPixel(px,py), this is not helping because the pixels that the method receives are refers to the map container pixels and not to the image pixels. there is a way to put features on the map by using the pixels (x,y) of the image itself ? **Thanks.**


title: Select polygons that do not intersect any other polygon using QGIS, GRASS, SAGA or R

text: I would like to select the polygons of the layer below that do not intersect with any other polygon ( _within the same layer_ ). How can I do this using QGIS, GRASS, SAGA or R? ![Image](http://i.stack.imgur.com/wkHEw.png)


title: Mapnik XML Stylesheet in QGIS

text: Is it possible to use a Mapnik XML Stylesheet to style/render data in QGIS? I know there is Quantumnik to generate Mapnik XML Stylesheets. Yet I am wondering if it is possible to import a Mapnik XML Stylesheet into QGIS to display data?


title: Richmarkers not accurate on lower zoom levels

text: I am trying to enable my icon to rotate and found this link Show wind direction on Google Maps. Unfortunately when I use the richmarkers it can only show the accurate location when I zoom in else its not showing the accurate location any solution to this? Below are snippet of my codes.               var arrowImage           = new Image();                             arrowImage.src           = 'alertIcon/G.png';                                  // rotation in degree                             var directionDeg         = 300 ;                                  // create a container for the arrow                             var rotationElement      = document.createElement('div');                             var rotationStyles       = 'display:block;' +                                                        '-ms-transform:      rotate(%rotationdeg);' +                                                        '-o-transform:       rotate(%rotationdeg);' +                                                        '-moz-transform:     rotate(%rotationdeg);' +                                                        '-webkit-transform:  rotate(%rotationdeg);' ;                                  // replace %rotation with the value of directionDeg                             rotationStyles           = rotationStyles.split('%rotation').join(directionDeg);                                  rotationElement.setAttribute('style', rotationStyles);                                       // append image into the rotation container element                             rotationElement.appendChild(arrowImage);                                  // append rotation container into the richMarker content element                             richMarkerContent.appendChild(rotationElement);                                  // create a rich marker ("position" and "map" are google maps objects)                             new RichMarker(                                 {                                     position    : point,                                     map         : map,                                     draggable   : false,                                     flat        : true,                                     anchor      : RichMarkerPosition.MIDDLE,                                     content     : richMarkerContent.innerHTML                                 }                             );


title: How to get extent out of GeoTiff

text: How do I get the extent out of the GeoTiff format? I want something like this:               Extent(293518.1886150768,5680494.194041155,890338.5054657329,6267530.571271311)      AKA the xmin, xmax, ymin, ymax coordinates from the GeoTiff? Edit: I'm implementing my own reader so I will need to know how the file format tells me this information.


title: Is it possible to interrupt Select Layer By Location when first feature selected?

text: I´ve got a point feature class with some 7.000.000 points and I´m trying to check if any of these points is simply intersected by several polygon feature classes using 'for' loop and 'arcpy.SelectLayerByLocation'. Polygon feature classes contain just 1 polygon feature each and are saved in a folder as shapefiles. I can create a feature layer from points and then use Select Layer by Location with "INTERSECT" parameter but my problem is that the polygon feature can be intersected by e.g. 700.000 points and the process of selection takes ages. My question is: is there any way to break the Select Layer by Location process after it finds a first point that is intersected by the polygon? Or do I have to wait until the whole process finishes? The code is simple:               # Create a layer from feature class with 7.000.000 points      arcpy.MakeFeatureLayer_management(pointFC, "point_lyr")          # Builds a list of polygon feature classes within a workspace     fcs = arcpy.ListFeatureClasses("", "polygon")     for polygon in fcs:         arcpy.SelectLayerByLocation_management("point_lyr", "INTERSECT", polygon, "", "NEW_SELECTION")              # And now instead of:              # Count selected points after the process ends and do something if condition is met         points_count = int(arcpy.GetCount_management("point_lyr").getOutput(0))         if points_count != 0:             something              # I´d rather have              if a point intersects the polygon for the first time:             tell me that this happened, then break this loop and start another polygon loop      Is it possible to break a running arcpy process? Running on Win7, ArcGIS for Desktop Basic 10.1 + Python 2.7


title: Use Mapnik tiles in QGIS

text: I want to use the QGIS Map Composer with Python to automatically generate maps. As a base map of these maps I want to use tiles generated with Mapnik. Currently I am serving the tiles with mod_tile just like the regular osm tiles. I know I can simply use my url of my tiles in QGIS, but the quality of the tiles is made for web maps, and I need higher resolution. So my question is, if it is possible to somehow export a georeferenced image with Mapnik to use it in the QGIS Map Composer. Or if there is an other way to create some sort of high resolution tile with Mapnik and use it in QGIS?


title: Remove TileStache cache for a single feature

text: In a large parcel map, features are edited and need to be re-styled. I want to keep as much of the cache as possible, but invalidate the cache for features that are edited when the features are saved. In order to do this, I need some way of determining which tiles (at all zoom levels) contain the feature. Assuming my features are in PostGIS and served up in Leaflet, what are some ways to determine that intersection?


title: Combining Hover and Click events in OpenLayers

text: I can specify event callbacks in an OpenLayers layer:               var vectorA=new OpenLayers.Layer.Vector("A");     vectorA.events.on({       "featureselected": callback     });      and I can specify event callbacks in a control's constructor:               highlightControl = new OpenLayers.Control.SelectFeature(vectorA, {       renderIntent: "temporary",       hover: true,       highlightOnly: true,       eventListeners: {         featurehighlighted: hoverFeatureCallback,         featureunhighlighted: unhoverFeatureCallback       }     });      Can I mix them properly? Basically, one layer needs both select and hover events, so I'd have two controls, one that hovers, and give each their own callbacks (like the second example above.) AND I have multiple layers that need select events, so in that case I'd define the events in the layers themselves, like the first example. However, I can't use the first example with hover and select because, in OL's infinite wisdom, hover events are really just select events. How can I have two layers respond to hover events and one of them also respond to select events?


title: Point z value Return 1.#QNAN while set before

text: I am makin a python code for ArcGis 9.3. The code is as following :               rowsLine = gp.updatecursor(line_Feature_Class)     count = 0     rowLine = rowsLine.next()          rowsPoint = gp.InsertCursor(point_Feature_Class)     feat = rowsPoint.NewRow()               while rowLine:          #...          ptObj = gp.CreateObject('Point')         ptObj.x = x1         ptObj.y = y1         ptObj.z = z1         feat.Shape = ptObj         rowsPoint.InsertRow(feat)          del rowsPoint          desc = gp.Describe(point_Feature_Class)     shapefieldname = desc.ShapeFieldName          rowsAlti = gp.UpdateCursor(point_Feature_Class)     rowAlti = rowsAlti.Next()          while rowAlti :              featAlti = rowAlti.GetValue(shapefieldname)         pntAlti = featAlti.GetPart()         Alti = pntAlti.z         rowAlti.SetValue("Altitude", Alti)              rowsAlti.updateRow(rowAlti)         rowAlti = rowsAlti.Next()      The problem is that the return value of pntAlti.z is : 1.#QNAN. But if I replace pntAlti.z by pntAlti.x, then it works ! I saw that 1.QNAN means the value is not set. However I set the z value in my code : `ptObj.z = z1` Thanks for your help


title: Dissolving polygons closers than 50m into a single shape QGIS

text: I’m new with GIS (coming from a CAD background) and using QGIS 2.0.1. I have a very detailed map of urban fabric in a UK region, to the detail of individual buildings. I want to simplify that map by uniting regions where buildings are closer than 50m to each other into a single shape. Areas with more disperse urbanisation should not be considered - polygons here should be deleted. The result would be a land use map with polygons delimiting areas of dense building, rather than individual shapes for each unit. I suspect that the “Buffer” tool would do the job, if I choose 50m as buffer distance and combine the buffers into one region. Before I set out to do it, am I thinking right about this? Or is there other way to do it accurately? Thanks!


title: QGIS 2.0 Performance issue when selecting by location on a complex dataset

text: I have an issue when I use the select by location tool in QGIS 2.0 to select points within a very complex polygon. It is flood data, so has something like 350,000 nodes and is generated as a single polygon. When you explode the multipart polygon it seperates to 1500 individual features. Whether I do this within the single feature or the multiple features, QGIS hangs and crashes. The two datasets use the same coordinate systems. They are cached WFS layers. Does anyone know if I'm doing something wrong or if this is a logged bug? Is there a workaround?


title: Is there a way to modify an Excel sheet which is already open in ArcMap?

text: I have an excel sheet that is joined to a shapefile in ArcMap. I need to modify many elements of the table, but every time I need to close ArcMap in order to be able to do it. Is there any way to modify the Excel sheet realtime?? I mean without closing ArcMap? thanks!


title: 3 points meteorology stations

text: I want to extract meteorological data of a point which are 3 meteorological stations around it. I found that ArcGis has the capability to do it by interpolation. My essential data is wind speed , wind direction and cloud cover which interpolation method or any other method is useful for my case?


title: How to draw bus routes on OpenStreetMap with mapnik

text: I have successfully set up Postgres with PostGIS and MapNik to serve map tiles. The database is populated with data from OpenStreetMap project. I can see that bus routes are in the data-set, but I have hard time displaying them on the tiles. Can anybody provide some advice how to do it? Thank you.


title: How to clip an Openlayers OSM Background to a Vector Boundary?

text: I'm new to QGIS. In my project i have two layers: one is my vector (city boundaries) and the other is the OSM (obtained with the plugin). How can I cut the OSM along the edges of the vector? I want to obtain only a region of the map, to save it as an image.


title: How do I clip OSM basemap with a polygon?

text: I've created a map in qgis 2.0. The map uses the OpenLayers plug-in to add the OSM layer as a basemap. I would like to restrict the basemap to be visible only inside a given polygon defined in a shapefile layer. How do I achieve this?


title: How to calculate distance of gauging station to outlet point in watershed?

text: ![enter image description here](http://i.stack.imgur.com/g0a8K.jpg) I have a gauging station point with elevation. Also have a watershed generated for the same area. I just want to know how to calculate the distance along the flowpath from the station to the outlet(point) The green point being the gauging station with an elevation of 100m and the red one being the outlet.


title: Full Easting and Northing coordinates values cannot be encoded in GRASS r.water.outlet in SEXTANTE-QGIS

text: The Easting or Northing coordinate parameter accepts only 5 digits. The Easting and Northing values for my outlet are 6 and 7 digits long (integer). How can they be encoded in r.water.outlet? Can real numbers be encoded as well? Thanks.


title: How to return text from GP service?

text: So I've got a website that I want to extend by letting people look up postal/zip codes for an address they enter. I've already got all the addresses (with their postal codes) for the city in an ArcSDE geodatabase. What I'd like to happen is that the user will enter a civic number, select a street from a dropdown list and then upon clicking a button the postal code for the address will appear in it's text box. I've created a script tool that takes two arguments: civicNumber and streetName and then uses a search cursor to get the postal code for the address. If it finds the postal code the tool uses arcpy.AddMessage() to return the postal code. This actually works, I've turned the tool into a GP Service and tested it from the web and the postal code is indeed being returned in a `esriJobMessageTypeInformative` JSON object. My question is this: Is there any way for the Python script to return the Postal Code as a string without using the `arcpy.AddMessage()` function?


title: Creating "Inner Buffer" using ArcGIS for Desktop?

text: I am currently searching for how to create an "Inner Buffer" in ArcGIS for Desktop? I have a lot of polygons (houses) in one feature class, and one big polygon (in another feature class). I want so select all houses contained in that one big polygon within a distance to the outer border (e.g. 1000 meters). _Select Layer By Location_ could be one possibility to solve this but I can't figure out how.


title: gdal_grid unable to get nodata outside of original data area

text: I have a vector based csv file that I want to convert into a grid using gdal_grid, but I want to force null data values strictly around the area where there is data. I can get gdal_grid to work but it is always just smoothing the whole grid rather than giving me null data outside of the shape. Can someone tell me what settings to use to make sure that I get null data inside the grid where there is no data within the original csv. Many thanks


title: How can I give random color to each line in Geoserver layer that is served from MSSQL geometry table?

text: I have a table of routes and I want to give each route random color to distinguish them on map (Using OpenLayers). What will be the easiest way to do so ?


title: Delete records from table using ModelBuilder

text: I'm very new to ModelBuilder so forgive me if there is a relatively simple tool to perform this. I have looked into Delete Field/Table but it does not seem to achieve the goal. Goal: Create a tool that will update the master parcel feature class on one geoDB with combined data from a feature class and table in another location. How do I go about deleting the data that is in the master parcel feature class without deleting the headers? I basically want it to be a blank template with headers to later append/join the newest information with. Or would it be simpler to delete the file and just save the joined information in its place? I'm also not sure if I have permissions to delete files, only modify, so finding a way to go about this without deleting would be great. * * * edit: I realize that I may have have explained this poorly. Our master parcel feature class is the result of a manual combination of a geometry feature class and a parcel data table that is provided by our county assessor office. The goal of the tool would be to automate this process. My current workflow: BlountParcels - Master Parcel feature class (what needs to be updated/overwritten) Parcels - Geometry feature class camagisdata - table parcel data ![enter image description here](http://i.stack.imgur.com/9wb2w.jpg)


title: Importing a shapefile into Postgres: Can I use Multiple CPU Cores?

text: I'm importing some shapefiles into Postgres/PostGIS. I'm using Windows 8 and the shapefile importer/exporter GUI. It maxes out one CPU core. I'm wondering if there is a way I can use multiple cores to increase the speed?


title: converting random point order into polygons

text: We received a kml file with polygons but the points are not in a logical sequence and we end up with strange overlapping shapes that are supposed to be an open polygon. We did identify a point near the center and found the azimuth and then sorted in that order and it helped but still did not solve. The bounds of the polygon do contain convex sides. Is anyone aware of an approach to process the random points and get the best open polygon? We have considered calculating bearing and distance from each point and looking for non intersecting line segments but wanted to see if there is a better approach. We are using php to calculate and then pass into mapping api.. our database is mysql. Trying to do as minimal programming as possible.


title: Calculate the majority of one polygon in another (fishnet)

text: I am looking for a method to calculate the majority of one polygon file (a map of ecoregions) in a sampling grid (fishnet) in ArcGIS. The purpose is to assign each sampling grid (fishnet) to an ecoregion. The determining factor will be the majority of area. I had a decent method of working out the problem using Tabulate Area tool, but the numbers were almost all identical indicating something weird was going on. I didn't trust it. I was also playing around with the intersect tool then dissolve, but couldn't find a way keep the attribute data I needed (in this case, the ecoregion I want to label). Problem: Label each sampling grid by the majority area of ecoregion within the sampling grid. Can someone help me out?


title: Automatic batch update a file geodatabase from shapefiles

text: Every week I need to update the same shapefiles that are used to produce up to date maps. The files are downloaded from 3 different web-servers. The thing I want to do is to have a script that can search my folder where I put my downloaded shapefiles and replace the ones I've already got. This shouldn't be a problem when the shapefiles have the same name. But still, I can't seem to figure out how to solve this problem. What to use, model builder? Python script? I bet this is a no-brainer question, but I'm stuck here :)


title: How to remove black edges on georeferenced images with gdal

text: after an image has been georeferenced the transformation creates black edges. How could I set these to white instead of black using gdal? would I have to somehow create a colour table and set nodata to white? yours, Rob


title: For every point, distance to nearest polygon of type X?

text: In ArcGIS how might one calculate as an attribute value the distance to the nearest polygon of type X? Say we have a point layer of wildlife observations and we want to know how far each of them is from the nearest wetland. Please note if your answer is specific to arcgis v9 or 10 and the necessary license level.