title: How to add wms layer in OpenLayers from a geoserver?

text: How to overlay a wms layer from geoserver (localhost)in Open layers? What is the basic syntax of Geoserver WMS Layer? Details are as follow ! workspace name: india   layer name:hospitals.                <script defer="defer" type="text/javascript">      var map = new OpenLayers.Map('map');       var wms = new OpenLayers.Layer.WMS( "OpenLayers WMS", "labs.metacarta.com/wms/vmap0";, {layers: 'basic'} );      var twms = new OpenLayers.Layer.WMS( "World Map", "localhost:8080/geoserver/wms? service=wms& version=1.1.1& request=GetCapabilities",       {           map: 'hospitals',           transparent: 'true',           layers: 'hospitals'      } );       map.addLayers([wms, twms]);      map.zoomToMaxExtent();       </script>       plz solve this!! =====   satya


title: Looking for (free) WMS of surface waters of Europe

text: Do you know any WMS with hydrographic map of Europe? Especially I am looking for map with surface waters of Europe


title: Map projection lat, lon to pixel?

text: I'm trying to find the concave hull of a lat, lon pointset but when I project my points on my map it doesn't look like expected. I've tried a lineare projection and a spherical projection. My data is about 10000-70000 lat lon pairs. How is this possible?               $mapWidth    = 2000;     $mapHeight   = 2000;     $mapLonDelta = $mapLonRight - $mapLonLeft;     $mapLatDelta = $mapLatTop - $mapLatBottom;          $worldMapWidth=(($mapWidth/$mapLonDelta)*360)/(2*M_PI);     $f=min(max(sin($mapLatBottom*(M_PI/180)),-0.9999),0.9999);     $mapOffsetY=$worldMapWidth/2 * log((1+$f)/(1-$f));     $f=min(max(sin($mapLatTop*(M_PI/180)),-0.9999),0.9999);     $mapOffsetTopY=$worldMapWidth/2 * log((1+$f)/(1-$f));     $mapHeightNew=$mapOffsetTopY-$mapOffsetY;     $mapRatioH=$mapHeight/$mapHeightNew;     $newWidth=$mapWidth*($mapHeightNew/$mapHeight);     $mapRatioW=$mapWidth/$newWidth;          $tx = ($lon - $mapLonLeft) * ($newWidth/$mapLonDelta)*$mapRatioW;     $f = sin($lat*M_PI/180);     $ty = ($mapHeightNew-(($worldMapWidth/2 * log((1+$f)/(1-$f)))-$mapOffsetY));      This should look like Germany. What's wrong?: ![Map](http://i.stack.imgur.com/6hqMW.png) Update: It seems to work. This is a map of Bonn(Germany): ![enter image description here](http://i.stack.imgur.com/HIEVW.png)


title: JS API - onDrawEnd firing twice

text: I have a simple editing application with the following code:               startEditing : function(template) {              var drawingTool = template.template.drawingTool;              switch(drawingTool) {             case FeatureTemplate.TOOL_POINT:                 drawingTool = Draw.POINT;                 break;             }              on(this.drawingToolbar, "draw-end", lang.hitch(this, this.createFeature, template));              this.drawingToolbar.activate(drawingTool);     }      and               createFeature : function(template, evt) {              var featureLayer = template.featureLayer;         template = template.template;              var prototype = template.prototype;              var geometry = evt.geometry;              var graphic = new Graphic(prototype.toJson());         graphic.setGeometry(geometry);              this.initAttributes(graphic, featureLayer).then(function() {                  var features = [graphic];                  featureLayer.applyEdits(features).then(function(addResults) {                 var objectIds = array.map(addResults, function(addResult) {                     return addResult.objectId;                 });                      var q = new Query();                 q.objectIds = objectIds;                      featureLayer.selectFeatures(q).then(function(features) {                     main.openForm(features);                 });             });              });          }      Using a point tool with the drawing toolbar, I am finding that when I click a point on the map, the createFeature function is called twice resulting in 2 features being created with the same geometry. How can I stop the draw toolbar from creating 2 features?


title: How to generate pgrouting source and target with OSM?

text: Is there anyone here can explain me more details for pgrouting's source and target? I'm wondering how to generate with the given lonlat in OSM PostGIS entities. In pgrouting's website, they are telling that source and target are identifier of the vertex but not telling if this an `osm_id` for OSM database or any. My only given to generate is lonlat of the potential source/target. I'm using osm2po to create routable database. Anyone can give me a clue, so I can make my map more dynamic? Thank you in advance,


title: Android native app - Geoserver

text: I am involved in a bigger university mobile application project and my task is to design and develop a mobile module regarding map and location-based services in a native Android application. Shortly my part of app should display different kind of info on a georeferenced map (Google Map or Openstreetmap) based on Oracle data:   * Shortest path to reach a certain event from a certain location   * Different layers with domain information to the event (like parking places around it, bus lines ecc. )   * GetFeatureInfo option I am new in this kind of topic and I would like to ask you an architectural question regarding it. I thought about using Geoserver to define layers, styles, datasets etc., and WMS protocol to generate the map over a website using HTML 5, OpenLayers, and jQuery-mobile. But how can I accomplish this in a native android app? Building a webpage where I display the map using HTML 5, OpenLayers and jQuery-mobile to interact with Geoserver and then using the Webviews in Android to display the webpage in my app? Would be this a possible solution? Possible drawbacks? What are the alternatives using Geoserver?


title: How do I open mxd backup file?

text: How do I open a mxdbk file ? My ArcMap frozen and save the mxdbk files. I would like to know how to open them up rather having to re-do the work. Thanks !


title: And/or query on WFS service?

text: Is it possible to use "and" & "or" queries on a WFS service through the URL? Here's my current URL:               http://mywebsite.com/county_roads_wms_only/MapServer/WFSServer?request=GetFeature&service=WFS&typename=county_roads_wms_only:Roads&maxfeatures=2&filter=<ogc:Filter><ogc:PropertyIsLike escape="\" singleChar="_" wildCard="%"><ogc:PropertyName>Road_Label</ogc:PropertyName><ogc:Literal>%mcculloch%</ogc:Literal></ogc:PropertyIsLike></ogc:Filter>'      Is there a way to form it so that it would be "%mcculloch% AND %st%" or "%mcculloch% OR %st% in the URL?


title: Issue with Buffer with projection

text: I am having an issue with buffer. My data is in Geographic coordinate system. I am trying to do a buffer of 600m and i am finding different results for each case, have tried thse so far: Map in web maercator map in UTM zone 12 map in geogaphic coordinate I am missing some data points due to this and wanted to check how I can resolve this issue. I know that UTM is the most accurate projection for distance but I wanted to know if by using WebMercator map if there is a way to get correct buffer results , any ideas? Thanks jay


title: help in field calculator

text: In Arcmap (10.1 version) I create a new field call "gender", and i want the content of "gender" will be changed according to others fields. For example, if a word(boy) is shown in a field call "person" in one row, then the content of "gender" will also change to boy (replace the original content). Here are my code(VB script):               IF [personA] = "boy" and [personB] = "boy" THEN         [gender] = "boy,boy"     EndIF          IF  [personA] <> "boy"  THEN         [gender] = "Non-boy,boy"     EndIF          IF  [personB] <> "boy"  THEN         [gender] = "boy,Non-boy"     EndIF      however, it cannot run and a message was shown.               A field name was not found or there were unbalanced quotation marks.     ERROR 999999: Error executing function.     Item not found in this collection.     Failed to execute (Calculate Field).


title: Syntax to add attributes in a line in the Expression String Builder

text: I have several columns with certain values, here's an example:                 | Apple | Banana| Orange | Pear       | Expected result       | 1     | 1     | 1      | 1          | red       | 1     | 2     | 2      | 1          | green       | 2     | 1     | 1      | 1          | red       | 2     | 2     | 2      | 1          | green      I would like to add **ALL** columns **ONLY** when the values are 2. I'm using this as a style filter so there will be colour ranges (i.e. `Apple` \+ `Banana` \+ `Orange` \+ `Pear`is less than 4 then red; equal or greater than 4 then green). Can this be done? The following is, I think, close to what I want but it gives `red` for all results instead of the expected results as shown above:               case when "Apple" IN ('2') + "Banana" IN ('2') + "Orange" IN ('2') + "Pear" IN ('2')     >=4 then 'green'     else 'red'     end


title: From where to download standard shapefiles

text: Is there any standard site for downloading legitimate shapefiles for US and Canada, Alaska?


title: Map Rendering Issues

text: I have tried simplification of geometry using Postgis function (ST_Simplify) in a map, for map rendering issues. I am still not able to get a good result, are there any other steps to follow?


title: StArea() in sqlserver is in which unit?

text: I have a problem with saving `STArea` in sqlserver. Here is a query where I get the area of a polygon:               select Geometry::STGeomFromText('POLYGON ((51.3386535086000000 31.7951253522000000, 51.3388695102000000 31.7952389270000000, 51.3388979249000000 31.7952464707000000, 51.3390467037000000 31.7953440361000000, 51.3392599393000000 31.7954571079000000, 51.3394563273000000 31.7955843452000000, 51.3396358676000000 31.7956905439000000, 51.3398606703000000 31.7958423402000000, 51.3399567269000000 31.7958984151000000, 51.3402088545000000 31.7960463557000000, 51.3404525165000000 31.7961833999000000, 51.3405501656000000 31.7962423246000000, 51.3405803405000000 31.7962540593000000, 51.3407558575000000 31.7963579111000000, 51.3407710288000000 31.7963713221000000, 51.3408634812000000 31.7964024190000000, 51.3410357293000000 31.7965043429000000, 51.3411784731000000 31.7965789419000000, 51.3412407506000000 31.7965510301000000, 51.3414824009000000 31.7962486111000000, 51.3417441677000000 31.7959164362000000, 51.3419324253000000 31.7956736963000000, 51.3419411425000000 31.7956510652000000, 51.3417580817000000 31.7955185473000000, 51.3413181156000000 31.7952600494000000, 51.3409249205000000 31.7950360849000000, 51.3404942583000000 31.7947866395000000, 51.3400418032000000 31.7945303209000000, 51.3396468479000000 31.7942999862000000, 51.3395502046000000 31.7943521217000000, 51.3392463606000000 31.7945348471000000, 51.3389606215000000 31.7947612423000000, 51.3386898860000000 31.7950401921000000, 51.3386374991000000 31.7951089237000000, 51.3386374153000000 31.7951182276000000, 51.3386535086000000 31.7951253522000000))',32639).MakeValid().STArea()      It gives me `3.49106150461542E-06`. How I can find that this area is right? In which unit is the `STArea()` measured in? And which code must I use instead 32639 for the points that came from the gps device? **Edit,This is how I read points **                var xmlDocument = XDocument.Load(xmlpath);                      XNamespace ns = xmlDocument.Root.GetDefaultNamespace();                 var values = xmlDocument.Root                               .Descendants(ns + "trkpt")                               .Select(x => new                               {                                   lat = (string)x.Attribute("lat"),                                   lon = (string)x.Attribute("lon"),                                   ele = (string)x.Element(ns + "ele")                                    })                               .ToList();                 //List<IPoint> ipoint = new List<IPoint>();                 string poinnts = "";                 foreach (var val in values)                 {                     string[] cordinates = { val.lat, val.lon, val.ele };                          poinnts += "," + val.lon.ToString() + " " + val.lat.ToString()+" " + val.ele.ToString();                      }                 poinnts=poinnts.Remove(0, 1);                      poinnts += ","+poinnts.Split(',')[0];                      GeomConverter gc = new GeomConverter();                 Polygon poly = gc.CreatePolygonFromPoints(poinnts);                           FeatureType res = poly.FeatureType;                 Fset = new DotSpatial.Data.FeatureSet(res);                      Fset.Projection = DotSpatial.Projections.ProjectionInfo.FromEpsgCode(Properties.Settings.Default.UTMCode);                      Fset.AddFeature(poly);                 area = poly.Area;                 length = poly.Length;                 this.status = true;                      return Fset;                       }      Here Is `Create PolygonFromPoint()`                   public Polygon CreatePolygonFromPoints(string feature)         {                       string[] coordinates = feature.Split(',');                       double[] xy = new double[coordinates.Length*2];             double[] z = new double[coordinates.Length];             for (int i = 0; i < coordinates.Length; i++)             {                 string[] cord = coordinates[i].Split(' ');                 double lon = double.Parse(cord[0]);                 double lat = double.Parse(cord[1]);                 double ele = double.Parse(cord[1]);                      xy[i*2] = lon;                 xy[i * 2 + 1] = lat;                 z[i] = ele;             }                  ProjectionInfo pStart = KnownCoordinateSystems.Geographic.World.WGS1984;                       ProjectionInfo pEnd = KnownCoordinateSystems.Projected.UtmWgs1984.WGS1984UTMZone39N;                       Reproject.ReprojectPoints(xy, z, pStart, pEnd, 0, coordinates.Length);                       List<DotSpatial.Topology.Coordinate> co = new List<DotSpatial.Topology.Coordinate>();             for (int i = 0; i < coordinates.Length; i++)             {                 co.Add(new DotSpatial.Topology.Coordinate(xy[i*2],xy[i*2+1],z[i]));             }                  Polygon polygon = new Polygon(co);             return polygon;                   }


title: Getting Annotation Text Geometry as Polygon Using Arcobjects

text: How can i get the rectangle boundary of an annotation object using ArcObjects? Once we select an object on the map, there is a boundary as shown in the first picture. Also is there an option to get the minimum rectangle that covers the text as shown in second picture? ![enter image description here](http://i.stack.imgur.com/3YqYf.png) I tried                   IAnnotationFeature annoFeature;         annoFeature.Annotation.Geometry;      but this only gives me a polyline geometry which seems to be the baseline of text


title: Is there a way to use gdal functions from OSGeo4W out of their shell?

text: I would like to use some gdal functions from OSGeo4W (gdal2xyz, and gdal_merge) for my work. I can use them in the OsGeo4W shell provided after installing the package. But I want to use them in Pythonout of the OSGeo4W shell. I have Python2.7, gdal, gdal-python binding installed (all for 64 bit). However, I don't have a clue on how to use those gdal functions out of their shell. The reason I want to do this is that I have a module created from ArcGIS ModelBuilder which I exported as a Python script. With this script, I want to use gdal2xyz, gdal_merge, etc. to continue data process. Any suggestion is very welcomed! Thank you!


title: How to change the symbology on many image files within a group?

text: I have a group of raster images from lidar data (so a DEM). When loaded in, they all have different symbology - so their color scheme and color ramp are all different. I would like to create a group with images that all have the same symbology, so that I can visually compare topography between the adjacent images. There are approximately 600 images within this group, so I don't want to change the color scheme for each image individually. I've tried changing the symbology of one image, and using the "apply symbology to layer" function, but this does not work because it seems my layer or group of images is not considered a raster layer... I also tried first making this layer a raster layer with "make raster layer", but the layer or group that I have is not considered a valid type of input. I'm using arcmap 10.0. Any advice or suggestions on how to do this would be most helpful.


title: zonal statistics in Python crashes

text: I have a python script that works (it's worked for several other projects in the past and it worked for my current project initially). But whenever I run it, Python crashes sporadically. It would go through some iterations (I'm calculating zonal stats as table for 540 files for 10 countries, so over 5000 files in all, it's worked fine for the first country and now it's not working). I'm not getting an error message or anything, it's just stopping. I mean I would get the: > > > ========= RESTART ================= but that's it. Then I would click back to the script, click F5 and restart it, it would stop randomly and I'll restart it again. I have no idea why it's doing this and there's no way that I can babysit the program. HELP?!


title: Maxent output clips to minimum extent of all input variables

text: I am experiencing an issue with the spatial extent of my output ASC files. I am using the UK Land Cover Map (LCM) 2007 at 25 metres resolution as a categorical variable in my Maxent modelling, alongside select continuous Bioclim variables interpolated (spline) to 25 metres and to the same extent. However, as the Bioclim files are originally at 1 km resolution, in some areas along the coastline they do not fully cover the sand dunes and mudflats covered the LCM which I am trying to model for. The 'extent' of the Bioclim layer is the same as for LCM but these missing cells are given a 'No Data' value in the input ASCs. The output Maxent ASC extents are clipped to that of the Bioclim layers but I want them to extend to the extent of the LCM (which is included in the analysis and is in fact the 'most important' variable). Coastal sand dunes are an important habitat for my species and omission of these areas will negatively affect the accuracy of further connectivity analyses based on my Maxent outputs. Is there a setting in Maxent to resolve this, or would I need to interpolate the input Bioclim layers outside of their original extent? If so, would this method have merit (bearing in mind that Bioclim is already an interpolated surface)? ![Maxent logistic output \(red to green gradient\) and underlying LCM 2007 \(monochrome\).](http://i.stack.imgur.com/Jw0SQ.jpg)


title: Filtering out a subpart of a polygon shape file

text: I have ArcGIS 10.1 for Desktop software. I am working with a polygon shape file (having n rows) in ArcMap. The project demands splitting this file into two polygon shape files having n1 and n2 rows such that n1+n2=n. I can filter out the data using a where statement on the attribute table, but the resultant data does not have the field called "Shape", due to which I cannot plot the polygons. I am new to ArcMap software.


title: ArgGis seems to snap but doesn't

text: We are building a project in ArgGis 10.2 were we create buildings using rectangles. Some of them need to be attached at one or more sides (they have common walls), so we would use snapping to snap vertex to edge to create a new building next to another. It seems to snap, but when zoomed in (1:0.3 scale) we can see that they are not snapped. at 1:1 looks good, but if we look even closer we see that the lines or point are not joined.


title: unable to open the QGIS Tab

text: I have download the QGIS Softarwe of “gis-osgeo4w-1.7.3-00624b3” and “QGIS- OSGeo4W-1.8.0-2-Setup” both from net and installed both in my system, My system is running window07 professional Service Pack1, 32 bit operating system. This is the new install in my system. Now when I am doing double click over given display icon , it is not responding no error is geting. Please suggest where I am doing mistake, my other friends are using this.


title: using WMS feature with OpenLayers object

text: I am sending query to a WMS address URL and getting results like following.               data: {        features: [...],        totalFeatures: 2,        type: "FeatureCollection"     }      I want to use features to calculate their lengths or area. I am creating a Vector and adding features in it.               (new OpenLayers.Feature.Vector(data.features[0].geometry)).geometry.getArea();     (new OpenLayers.Feature.Vector(data.features[0].geometry)).geometry.getLength();     (new OpenLayers.Feature.Vector(data.features[0].geometry)).geometry.getCentroid();      But not working. Is there a practical way?


title: Issues creating slope raster layer from a point layer

text: I'm trying to create a slope raster layer using QGIS 2.0 on Windows XP. I started with a points layer with the Z elevation attribute. ![points](http://i.stack.imgur.com/0U9zH.jpg) I then used the interpolation plugin to create a DEM. ![enter image description here](http://i.stack.imgur.com/hhYHp.jpg) But, when I used the Raster Terrain Analysis plugin, it gave me a strange output. What can I be doing wrong? ![enter image description here](http://i.stack.imgur.com/CkWaL.jpg)


title: can I get map's bounding box at any time, after user clicks button

text: I use openlayers. I am creating a form that the user can insert a bounding box to the database. In my page I have the form and a map, so the user can preview what he entered. So I thought it would be nice if there was also the other way around : User hits a button like "Get current bounding box", and that gets the bounding box from the map at any given time and transport it to the form. The numbers of the bounding box that is. I'm talking about the bounding box that the users see, not the whole map's box. After that, user can hit "Save" and submit the form, or can clear the form, or edit it, or zoom/pan the map and hit "Get current bounding box" again. I think you got the idea... Is that possible? How can I get map's bounding box at any given time? Thanks in advance


title: Update PostGIS after updating GDAL?

text: I am using PostGIS 2.0 and I want to use the raster functionality through mapserver. It seems as if the latest official GDAL release 1.9.1 has an unstable PostGIS raster driver, as pointed out in a previous question. After posting a ticket, was advised on github to update GDAL from trunk for this to work. I did that but I still get the same error. Thus, my question is: Is it necessary to build PostGIS again from source after updating GDAL for that update to become effective in PostGIS?


title: How to lock map labels for inset map in Print Composer?

text: OK, so I am making a print layout which has two maps, one an inset map with just a skeleton of a line layer and the other is the detail map at higher resolution. In the detail map, I want to show the labels for each arm of the skeleton. For the inset map, I want the labels to be off. Even when I 'lock layers for map item', the labels still update in the inset map when the map is exported to PDF. The update preview is set to Cache for both, but I don't think this has an impact of what gets exported.


title: Couldn't load plugin 'processing' on windows 7 64bit

text: I installed QGIS (2.0.1-Dufour Dufour, d94c044) using both the stand alone installer and the OSgeo4w installer and seem to have an issue with the processing plugin. At first I thought that there was a python path issue, but now I am wondering if everything installed correctly. I don't know where the libraries live, so I would like some help debuging this issue. Here is what I know:   1. I get the error message below when I start up QGIS   2. The following folders do not exist on my computer:      * C:\Users\gstein/.qgis2/python     * C:\Users\gstein/.qgis2/python/plugins     * C:\PROGRA~1\QGISDU~1\bin\python27.zip     * C:\PROGRA~1\QGISDU~1\apps\Python27\lib\plat-win     * C:\PROGRA~1\QGISDU~1\apps\Python27\lib\lib-tk Any thoughts about what may be wrong? **Error Message** > Couldn't load plugin 'processing' from > ['C:/PROGRA~1/QGISDU~1/apps/qgis/./python', 'C:\Users\gstein/.qgis2/python', > 'C:\Users\gstein/.qgis2/python/plugins', > 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\pip-1.2.1-py2.7.egg', > 'C:\Python27\ArcGIS10.1\Lib\site-packages', 'C:\Program Files > (x86)\ArcGIS\Desktop10.1\bin', 'C:\Program Files > (x86)\ArcGIS\Desktop10.1\arcpy', 'C:\Program Files > (x86)\ArcGIS\Desktop10.1\ArcToolbox\Scripts', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\win32', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\win32\lib', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\Pythonwin', > 'C:\PROGRA~1\QGISDU~1\bin\python27.zip', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\DLLs', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\plat-win', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISDU~1\bin', > 'C:\PROGRA~1\QGISDU~1\apps\Python27', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\PIL', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32\lib', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Pythonwin', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site- > packages\Shapely-1.2.18-py2.7-win-amd64.egg', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg', > 'C:\PROGRA~1\QGISDU~1\apps\qgis\python\plugins\fTools\tools'] > > Traceback (most recent call last): File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 182, in > loadPlugin **import** (packageName) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing__init__.py", > line 20, in from processing.tools.general import runalg, runandload, > alghelp, alglist, algoptions, load, extent, getobject File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing\tools\general.py", > line 29, in from processing.core.Processing import Processing File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing\core\Processing.py", > line 43, in from processing.algs.QGISAlgorithmProvider import > QGISAlgorithmProvider File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing\algs\QGISAlgorithmProvider.py", > line 86, in from processing.algs.CreateConstantRaster import > CreateConstantRaster File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing\algs\CreateConstantRaster.py", > line 1, in from processing.core.RasterWriter import RasterWriter File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing\core\RasterWriter.py", > line 28, in import numpy File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:\Python27\ArcGIS10.1\Lib\site-packages\numpy__init__.py", line 137, in > import add_newdocs File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:\Python27\ArcGIS10.1\Lib\site-packages\numpy\add_newdocs.py", line 9, in > from numpy.lib import add_newdoc File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:\Python27\ArcGIS10.1\Lib\site-packages\numpy\lib__init__.py", line 4, in > from type_check import * File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:\Python27\ArcGIS10.1\Lib\site-packages\numpy\lib\type_check.py", line 8, > in import numpy.core.numeric as _nx File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) File > "C:\Python27\ArcGIS10.1\Lib\site-packages\numpy\core__init__.py", line 5, in > import multiarray File > "C:/PROGRA~1/QGISDU~1/apps/qgis/./python\qgis\utils.py", line 453, in > _import mod = _builtin_import(name, globals, locals, fromlist, level) > ImportError: DLL load failed: %1 is not a valid Win32 application. > > Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit > (AMD64)] > > QGIS version: 2.0.1-Dufour Dufour, d94c044 > > Python path: ['C:/PROGRA~1/QGISDU~1/apps/qgis/./python', > 'C:\Users\gstein/.qgis2/python', 'C:\Users\gstein/.qgis2/python/plugins', > 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\pip-1.2.1-py2.7.egg', > 'C:\Python27\ArcGIS10.1\Lib\site-packages', 'C:\Program Files > (x86)\ArcGIS\Desktop10.1\bin', 'C:\Program Files > (x86)\ArcGIS\Desktop10.1\arcpy', 'C:\Program Files > (x86)\ArcGIS\Desktop10.1\ArcToolbox\Scripts', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\win32', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\win32\lib', > 'C:\Python27\ArcGIS10.1\Lib\site-packages\Pythonwin', > 'C:\PROGRA~1\QGISDU~1\bin\python27.zip', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\DLLs', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\plat-win', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISDU~1\bin', > 'C:\PROGRA~1\QGISDU~1\apps\Python27', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\PIL', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32\lib', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Pythonwin', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site- > packages\Shapely-1.2.18-py2.7-win-amd64.egg', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', > 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg', > 'C:\PROGRA~1\QGISDU~1\apps\qgis\python\plugins\fTools\tools']


title: Looking for data sources for real road speeds

text: A data question really, though there may be a useful geo-processing solution. Context: I am calculating "drivetime" catchment areas using PgRouting and OSM / OSM2PO. If my employer/client is to use this as information for decisions, then I think the base data needs to be good. Problem: The fact is that in a large city say, London speed limits exceed traffic speeds due to congestion. Therefore, my polygons are not really valid (except at times of low activity). Question: Thus, I'd greatly value any help (all ideas good) with how to update my OSM speeds on an area by area or road by road basis. (E.g: Is there crowd- sourced data (from mobiles)? Can I get access to traffic flow data? Has somebody else tackled this?). Many thanks, John


title: How to find polygons inside polygon with PostGIS 2.0

text: I have a relation representing all the counties in my country and the other one containing all the districts - e. g. each and every county is composed of several districts. How do I select them with spatial query? I tried with               select nazorp from kraje as k, orp_wgs as o WHERE ST_Within(o.geom,k.geom) AND     k.nazev = 'Liberecký'      but had wrong results returned (well, they are probably not wrong, they are just not what I expected them to be). ST_Within only returns districts that don't share a boundary with the county, but I need to get all the districts within the county. Is that possible? I haven't found any built-in function suitable for my needs yet.


title: Convert polygons in shape file to a numpy array?

text: I have a shape file containing a number of polygons and all attributes being text. I'm trying to import the polygons into numpy as an array where each polygon is represented as unique values. I approach this by using gdal_rasterize to generate a geotiff, which I then can convert to an array:               gdal_rasterize -a provcode -l longhurst PG:'host=localhost dbname=biomes' -tr 0.25 0.25 out.tif      and               tif = gdal.Open('out.tif')     tifArray = tif.ReadAsArray()      My problem is that all polygons get the same values since provcode is of type string. How can I make gdal_rasterize burn different values to represent different polygons? Alternatively, is there a better way to do this conversion?


title: Aligning layers with questionable metadata in SW Ontario

text: After making my best efforts to read/try everything I can on this forum and elsewhere, I unfortunately have to contribute to the wide literature of "my layers aren't aligning" questions. I am working in SW Ontario & have some somewhat reliable reference data in NAD 1983 CSRS (Canadian Spatial Reference System). Then I have two shapefiles of points: 1). One that I created from a raw spreadsheet of lat-longs. The source reports (somewhat half-heartedly) that it should be in NAD 1983, so that's the geographic coordinate system I've defined for it. 2). One that was sent to me already as a shapefile, projected into NAD 1983 Great Lakes Basin Albers. I think it's quite possible this was not the data's "original" coordinate system, but I don't know. I've gotten the original raw data and can try to plot it myself. Ultimately I want everything to be in Great Lakes Basin Albers. Between these 2 datasets, certain points are supposed to be redundancies. When I hone in on such points, they have the same lat-longs (in the fields that were used to plot them), but the point from #2 is always 13-23 meters NE of the point from #1. I say my reference data (streams) is "somewhat reliable" b/c, while I know it's projected right, the points aren't necessarily going to be right on top of the stream lines, though the points from #1 often are, while #2 is offset. Thus, I've assumed my streams and #1 are okay. This kind of looks like a central N. America NAD83 / NAD27 offset to me, but when I reproject #2 to a geographic coordinate system of NAD27, it only gets further away to the northeast. Throw in the datum transformations needed when the CSRS version of NAD comes into play, and I'm quite unsure how to diagnose my problem in a systematic way or fix it, i.e. to get these points with identical lat-longs to actually be on top of each other, preferably at the location of #1. EDIT: The #2 points appear to have an attribute for "Geodetic Datum" specifying either NAD 1927, NAD 1983, or WGS84 (& very often blank). They're museum records covering roughly 1930-2010, so this kind of makes sense. Does this mean their lat/long attributes are all based on different reference points & the shapefile can't be fixed in a uniform way? ![enter image description here](http://i.stack.imgur.com/dfEvT.jpg)


title: Read and Export OSM Road Data as a Network of Road Segments [Adjacency Matrices]

text: I wish to extract from an OSM file a road network, and export it as an array (adjacency matrix) where the nodes are the road intersections and the edges represent the road segments. I don't want to just visualise it, but also export the resulting adjacency matrix as a file with a network data structure in at least one of the following formats: numpy file (.npy), pajek file (.net), gephi file (.gexf), csv file (.csv). Alternatively, a list of nodes and edges with geometric info would do. A solution using either Python or QGIS would be highly appreciated. Edit: I was hoping for something like this Matlab package, especially the 3rd item in the description. Unfortunately, I don't have Matlab installed on my machine.


title: How do I calculate the distance among all points in set in QGIS?

text: I am a very new user of any GIS. In QGIS I have input my UTM data as a CSV file. Under Vector, Analysis Tools, Distance Matrix I have been able to create a distance matrix but it just takes point 1 and calculates the distance to each other point. I Want the distance between all possible points. Is there a simple way to do this?


title: Folium GeoJSON Overlay Not Appearing

text: I am in the midst of a task which involves the clustering of assessment neighborhoods into larger neighborhoods (data link). The goal is build these neighborhoods in a manner consistent with major economic corridors in DC. I'd like to provide a dynamic method of selection for others in the office, and my chosen approach is the overlay of assessment polys on a map tile via Folium. Here is the set up:               #Read in data     aneigh=gp.GeoDataFrame.from_file(data_dir+'ASNBHDPLY.shp')          #Transform to WGS84     aneigh_ll=aneigh.to_crs(crs={'proj':'longlat',                                  'datum':'WGS84'}).rename(columns={'geometry':'geom_temp'})          #Define function to swap lat-long coords in Polygon objects     def swap_coord(poly):         '''Mappable function that switches coordinates in polygon objects'''         #If the object is not a Polygon...         if poly.geom_type != 'Polygon':             #...take the first component and switch the coordinates             new_poly=Polygon([(x[1],x[0]) for x in poly[0].exterior.coords])         #...otherwise...         else:             #...swap the coordinates of the Polygon             new_poly=Polygon([(x[1],x[0]) for x in poly.exterior.coords])         return new_poly          #Swap coordinates     aneigh_ll['geometry']=aneigh_ll['geom_temp'].apply(lambda x: swap_coord(x))          #Drop geom_temp     an_ll=aneigh_ll[[x for x in aneigh_ll.columns if x != 'geom_temp']]          #Write transformed data to disk     an_ll.to_file('asmnt_neigh_ll2.json',driver='GeoJSON')      First, I know the transformation above works because 1) I can plot one coordinate of one of the constituent polygons, and 2) I used it before. I should also note that the `to_file()` method that writes the GeoJSON to disk is from the geopandas package (which relies on fiona for I/O operations). Where I am running into trouble is getting the overlay to work. I know overlays are possible with my setup because I can recreate the example provided here. In fact, I wrote the set to disk to follow as closely as I could to the example. However, when I proceed with the following, I get only the base map.               wd='/home/choct155/dissertation/MiscProj/'          dcmap2=folium.Map(location=[38.9072309,-77.0364641])     dcmap2.geo_json(geo_path=wd+'asmnt_neigh_ll2.json')     dcmap2.create_map(path='dcmap.html')     showmap(dcmap2,'dcmap.html')      Any thoughts on why this might be the case? I am wondering if there is a fundamental difference between the GeoJSON files fiona has generated and those used on the folium site.


title: Con() does not recognize raster values

text: I have a problem when using the `con()` argument. What I want to do is to remap all raster of a raster file ("Red_SL") cells with values larger than 40 to the new value 5. All other values should be replaced by another raster file called "remap2". "remap" 2 contains only 4 definite values (1,2,3,4), and "Red SL" has continuous values between 0 and 90. My syntax is the following:               import arcpy     from arcpy.sa import *          arcpy.env.overwriteOutput = True     arcpy.env.workspace = r"D:\\Projekt\\Testprojekt\\Daten\\Ergebnisse\\"     arcpy.CheckOutExtension("Spatial")     remaped2 = "remaped2.img"     slope = "Red_SL"          CF= Con((slope < 40), 5, remaped2)     CF.save("Red_SL_Neigung_int.tif")     arcpy.CheckInExtension("Spatial")      now the problem is: I do not get an error message but the produced file is nonsense! Python is somehow not able to recognize the slope value in: CF= Con((slope < 40), 5, total) This means it doesent matter what value I give python (here I set it to 40, but I will get the exact same result if I give values of 1 or 3000), the result will always be the a new file which is exactly the same as "remaped2". Python therefore thinks all values are below 40 (or whatever number I give), which is not true. The exact opposite happens if I change the "<" to ">" `(CF= Con((slope > 40), 5, total))`. Now I will get i new file, which has only 5 in every raster cell. Here again the value whch I type in (instead of 40) does not matter, any value between 1 and 3000 will produce the exact same result. I double checked, but the file "Red_SL" has all values between 0 and 90, therfore the result I get from python is nonsense. I already tried (without sucess): -Change "Red_SL" from floating to integer -try differend formats (tif, grid, img) for "Red_SL" \- replace comma in floating raster with points (german vs. englisch delimiters) Can anyone see what is wrong with my use of the Con arguement?


title: Select and export by attribute iteration

text: I'm trying to do a loop for select and export layers by attribute and name them based on the field value. I've written the script below and it creates the shapefiles which picks up field values as shapefiles' name but they are empty output shapefiles. looks like the query in the select process isn't run. Any help would be appreciated.               arcpy.env.workspace = r'C:\temp'          # Create the search cursor     cursor = arcpy.SearchCursor('fc', '')          # Iterate through rows in cursor     for row in cursor:         outfc=r'C:\temp\test'+"\\"+str("{0}".format(row.Project))         arcpy.Select_analysis("fc",outfc,'"Project"=\'"{0}".format(row.Project)\'')          ### Project is the field name which is a String.


title: How to make a geo reference image using GeoTool based on a matrix mapped data?

text: My ultimate objective is to make a geo reference image and make it as a raster layer to overlay the image above my openlayers map. I have a set of geo coordinate data, and each coordinate is assigned to a value. The value is used to represent a color code.               cell           coordinate           value        Color     [x0 , y0]      [125.0 , 36.0]         1          Yellow     [x0 , y1]      [125.0 , 36.1]         1          Yellow     .     .     [x1, y0]       [125.1 , 36.0]         1          Yellow     [x1, y1]       [125.1 , 36.1]         2          Orange     .     .     [x2, y0]       [125.2 , 36.0]         1          Yellow     [x2, y1]       [125.2 , 36.1]         2          Orange     [x2, y2]       [125.2 , 36.2]         3          Red     [x2, y2]       [125.2 , 36.3]         3          Red      Hence, I will have a matrix mapped data which looks like below: ![matrix mapped data](http://i.stack.imgur.com/jOhij.png) Now I want to create a geo reference image with inner projection : `EPSG 4326` and outter projection : `EPSG 900913` **Can this be done by using GeoTool? If can, then how?** Help is needed.


title: Free and global basemaps for download

text: I'm pretty sure that a lot of people might have faced the situation where they search the internet for geographical basemaps on a global scale but can't hardly find the needed data in a short time and with sufficient quality and/or documentation. I think especially for beginners it is hard to know where to download good data because their are just too many websites and data-sources out there of which some are for free and others are not, some are of good quality and others not...and often all of them lack proper documentation. So my question is: Is there any specific website that bundles information on where to find good basemaps (such as administrational boundaries, cities, physical maps, relief maps, waterbodies, vegetation, soil and others) **on a global scale that also include proper documentation and reliable data- quality**? For me a good starting point is Worldgrids.org. check out here: http://worldgrids.org/doku.php people interestet in soils should consider: http://www.globalsoilmap.net/ some layers from the FAO: http://www.fao.org/geonetwork/srv/en/main.home (If you go to the portal make sure to do the advanced search and click on downloadable data as a filter option)


title: Openlayers KML script loads slow

text: Hy, I made a map with multiple KML layers with large datasets (total 500.000 markers). It works fine but as expected it takes a very long time for the script to load al markers. Is there a way to reduce this load to the script so it will only "load" the markers within the current view instead of loading the complete KML file? Exporting the KML to other formats is possible but i don't have a clue to what. I read allot about BBOX etc but i just can't get to wrap my head around this.


title: OpenLayers drawing features on clustered layers

text: I'm having problems adding and moving features on a clustered (wfs) layer. When I add a new feature all the clustered ones disappear from view. Plus much hacking needed to actually save the new feature. Bit of code:                   map = new OpenLayers.Map('map', {             projection: "EPSG:900913",             displayProjection: "EPSG:4326",     });     var base = new OpenLayers.Layer.OSM();          filterStrategy = new OpenLayers.Strategy.Filter();          cluster = new OpenLayers.Strategy.Cluster();          protocol = new OpenLayers.Protocol.WFS({         version: "1.1.0",         url: "/geoserver/test/wfs",         featurePrefix: "test",         featureNS :  "test",         featureType: "XXXX",         geometryName: 'geom',     });          wfs = new OpenLayers.Layer.Vector("Vector", {         projection: "EPSG:4326",         strategies: [              new OpenLayers.Strategy.Fixed(),             filterStrategy,             cluster         ],         protocol: protocol,     });          map.addLayers([ base, wfs ]);          drawControl = new OpenLayers.Control.DrawFeature(         wfs,  OpenLayers.Handler.Point,                           {               featureAdded: featureAdded,         }     );          map.addControl( drawControl );          function featureAdded(feature) {         alert("feature added!");     }      I have similar problems if I want to move a feature. Any help appreciated Mini


title: Is there a way to download and import labels for the CanVec Named Features (Code 1580010) layer (Natural Resources Canada)?

text: I used the CanVec symbology toolkit (described here) to symbolize all the downloaded data but now I am wanting to downland and import labels or a procedure that I should follow for creating all necessary labels. I have an add-on that I can use to import labels from another layer so ideally I would like to obtain a layer file that already has all the labels classes defined and properly symbolized.


title: PostGIS: Approximating a house number from address range

text: I'm developing a reverse-geocoder for Canada. So far, given a lat/lng, I can find the nearest street (line segment), which includes line segment direction and address ranges for both sides of the street. I'm now trying to figure out the best way to programmatically approximate the nearest house number to the given lat/lng point. Here's an example of a row containing the street data:               -[ RECORD 1 ]-     [...]     l_adddirfg | Same Direction     l_hnumf    | 3219     l_hnuml    | 3235     l_stname_c | Breen Road North-west     r_adddirfg | Same Direction     r_hnumf    | 3224     r_hnuml    | 3236     r_stname_c | Breen Road North-west     the_geom   | 0105000020E610000001000000010200000002000000B0F6990E78885CC088DF2B5F3C8C49400875B39A89885CC0A0BCA6AC4B8C4940      So, given a lat/lng coordinate that lies near the "the_geom" line segment, a person could tell visually which side of the street the point is on (left or right side), and how far along the segment it is -- thereby approximating a house number. For example, if the point lies on the right side, three-quarters down the street, I would use the fields r_hnumf (right side, first number) and r_hnuml (right side, last number)... The street address is probably close to:               3232 Breen Road North-west      What I'm looking for is a best practice in either computing/approximating this in PostGIS (which I'm new to), or in the application layer once the row is fetched. Thanks!


title: QGIS Line End Caps

text: When styling a road network where all sections of the network between junctions are separate paths, I style the ends of the lines to rounded to remove any gaps that might appear where the joins are at angles. The only problem is that you then have rounded ends to cul-de-sacs or dead end side roads. ArcGIS has a tool that allows you to specify the end caps to either square or butt, which removes the rounded ends. Is there a method of doing this in QGIS?


title: Unit of Google Maps API geometry services

text: When i call the method               google.maps.geometry.spherical.computeDistanceBetween(from, to);      it returns a number, is this distance in KiloMeters or Miles?? How can we find this and doc on this API?


title: Doubts with Ripley's K function Interpretation

text: I'm a self-taught person, I'm studying the K Ripley's function and I would need some person who give me some advices. First, I have written a little script of R in QGIS processor:               ##Point pattern analysis=group     ##Layer=vector     ##Titulo=string     ##showplots     library("maptools")     library("spatstat")         K <- Kest(as.ppp(Layer))     plot(K, main=Titulo)      Then, to test the different results I have prepared 3 layers about the same extension. One regular, other random, and the last clustered, with 100 points each one. This is the random graphic ploted: ![enter image description here](http://i.stack.imgur.com/UplbK.png) This is my question: The random graph should be equal than the poisson line, Why the observed K is higher than poisson one? (I see the lines shapes are equal) Is this due to the little number of points? Here are the files to reproduce it (included script for QGIS). Zip file


title: Drag and drop the chart widget to new browser window in flex viewer

text: I am trying to drag and drop the chart widget to new browser window. I have using DragManager to drag and drop the chart widget. But I does not drag the chart widget to another browser window. How can i achieve it, Is there any possibilities in flex.


title: How to have all nodes of my search in google earth as kml file

text: I need data for my research, and i am trying to get data from google earth to find densities of the facilities. I am using google earh desktop program and i activate 7 days trial of pro version. Yet, i have problem to create my search as kml file. When i search a keyword, ii find plenty of nodes relevant to my query, but i can see 20 pages of the result (every page has 10 point), thus, i can reach only 200 nodes and it is not enough to make predictions. When i click to page 21, or 22 or 23 i see the page 20 and i cannot reach other pages. The screen shot where i stuck can be seen as follow : ![enter image description here](http://i.stack.imgur.com/2Lnjd.jpg) How can i reach all the results as kml file?


title: Retrive data from database in openlayer GIS

text: i am using this code but it shows nothing, whats wrong with this code.               <!DOCTYPE html>                 <?php        $conn_string = "host=localhost  dbname=postgis_21 user=postgres      password=1234";        $dbconn4 = pg_connect($conn_string);        //$result = pg_exec();        $rs = pg_query("select ST_AsText(geom) from pak_adm2_pco_20110324");       $row = pg_fetch_array($rs, 0, PGSQL_NUM);       print_r($row[0]);        ?>           <meta charset='utf-8' />           <?php /*      */?>


title: GUI tools for PostgreSQL/PostGIS?

text: as a cartographer I work with TileMill to visualize geodata. I installed a local PostgreSQL server with PostGIS and osm2pgsql. Everything works fine. Anyhow, I am new to working with servers and the command line approach. So, in learning the basics of database management, some GUI tool would be of great help. Could I for example safely run the PgAdmin app next to the PostgrSQL installation (http://postgresapp.com/) on my Mac OS X Mavericks? Any suggestions most welcome.


title: Arcgis 10.2 - Table was not found error in a personal geodatabase

text: I'm working with a geometric network within a personal geodatabase. I've imported a feature class the was topologically fixed from one PGDB to another. Than created a geometric network with weights. When I've tried to add new weights I've received an error (don't have a print screen). I've restarted ArcMap, but the .mxd couldn't locate the feature classes. When I've located it for ArcMap I've received the following error. (It is the second time that it happens - exporting the feature class before forming a geometric network didn't work). What have I done wrong? ![enter image description here](http://i.stack.imgur.com/BC0D5.png) Thanks


title: How to create table containing contents of a directory with FME 2014?

text: I have to match the contents of a directory with the attributes of a FeatureClass. The FeatureClass field contains the names (without file endings) of documents which are either .pdf or .jpg files. I have two directories - one for PDFs and the other for JPGs. I need to match the PDF file names with the FeatureClass Attribute and if they match rename them with the PDF file ending. The ones that are not matched are then obviously JPG. How can I create the list of files using FME2014?


title: Copy parts from several attribute tables in shapefiles?

text: I have several shapefiles from which I wish to take parts from their attributes tables to make a new attributes table. Is this possible without destroying the original attributes table? I have features connected to them. If I were to use the join tool, how do I join from several layers? How can I do this?


title: Using arcgis.rand() with variables from code block of Calculate Value tool?

text: Within a larger ArcGIS ModelBuilder model I am trying to populate a field in an attribute table with a random number, using Calculate Value. The overall distribution of random numbers needs to follow a Normal distribution. The parameters of that normal distribution (mean and standard deviation) need to be calculated from model variables specified by the user. Is this possible? The arcgis.rand() function can do what I want and works perfectly if I specify the mean and stdev directly in the Calculate Value expression:               arcgis.rand("Normal 6 2.5")       If I use model variables directly, it also works:               arcgis.rand("Normal %mean% %stdev%")      However, in some instances I need to do some simple calculations to obtain the parameters from variables the user has entered. For example, they have already entered variables A, B and C, and I need the mean to be A*B and the stdev to be B/C. I have attempted to do this via the code block:               import random     MN = %A%*%B%     SD = %B%/%C%      and then the expression               arcgis.rand("Normal MN SD")      but this always results in values following the default distribution (mean 0, stdev 1). Troubleshooting shows this happens even if I define the code block variables as simple numbers, suggesting the problem seems to be using code block variables in the expression, not in the definition or calculation of the variables themselves. My question is: how can I do what I need using arcgis.rand()? I am open to suggestions of other methods/workflows for achieving what I need, but I am a novice at Python, in case that is not already obvious! I am using ArcGIS Desktop Advanced 10.2.2.


title: Raster to Vector

text: I want to map out moisture availability using satellite image but I do not know how to convert the raster image to a thematic map with limiting boundaries like state boundaries. I am using ArcGIS 9.3 and ENVI 4.3. An example of expected result is shown in a link provided below. I have satellite image who follows the same color gradient and I want the image to be displayed in shapefiles according to certain boundary shapefile. http://www.drought.unl.edu/dm/monitor.html ## EDIT I have posted a sample MODIS image here and a sample expected result of a drought map accessed from UNL.edu. ![](http://i.stack.imgur.com/W2bC4.jpg) ![](http://i.stack.imgur.com/rBUMg.gif) So, I am not sure about the classification. Please do enlighten me on this matter.


title: Changing attribute values in a memory layer

text: I'm trying to change attribute values in a memory layer. My memory layer is defined in this way:               pro.addAttributes([QgsField("nom", QVariant.String), QgsField("valeur", QVariant.Int)])      I made several tries:               attrs = { 0 : "hello"}     bFlag3=vl.dataProvider().changeAttributeValues({ fet.id() : attrs })      Or in this way :               bFlag3=vl.changeAttributeValue(fet.id(),0,'hello')      Each time my boolean flag returns **True** but I never see any change in my memory table. I make a               vl.commitChanges()      in both cases. When I try these two pieces of code on the two other layers (which aren't memory layers...) I have in my project, I see the attributes change... Can anyone help me troubleshoot this?


title: GeoExt 2 ComboBox customization

text: First question so will try to be as thorough as possible. I'm trying to customize the GeoExt.form.field.GeocoderComboBox to use Here maps geocoder. I can get the label field to display in the combo box and also convert the position into the format I need however I am unable to actually select the value in the combobox. It looks like the 'select' event is not firing and I have no idea why. It looks like this does fire when I go with the default queryParam "q" although obviously it won't return what I actually want. So it would seem that me changing the queryParam to the one required by Here maps for some reason is stopping me selecting the address I put in. Any help would be much appreciated. Here's my code:               Ext.define('HQMap.view.tools.HereGeoCode', {     extend: 'GeoExt.form.field.GeocoderComboBox',     requires: [         "Ext.data.proxy.JsonP",         "Ext.data.JsonStore"     ],     initComponent: function() {         var me = this;         var hereUrl = "http://geocoder.api.here.com/6.2/search.json?app_id=&app_code=&gen5&country=GBR";         var hereStore = Ext.create("Ext.data.JsonStore", {             root: null,             fields: [                 {name: "bounds", convert: function(v, rec) {                     var bbox = rec.raw.Response.View[0].Result[0].Location.MapView;                         return [bbox.TopLeft.Latitude, bbox.BottomRight.Longitude, bbox.BottomRight.Latitude, bbox.TopLeft.Longitude];                         }                 },                 {name: "lonlat", convert: function(v, rec) {                     return [rec.raw.Response.View[0].Result[0].Location.DisplayPosition.Longitude, rec.raw.Response.View[0].Result[0].Location.DisplayPosition.Latitude];                     }                 },                 {name: "name", mapping: "Response.View[0].Result[0].Location.Address.Label"}             ],             proxy: Ext.create ("Ext.data.proxy.JsonP", {                 url: hereUrl,                 callbackKey: "jsoncallback"                 })         });              Ext.apply(me, {             srs: "EPSG:900913",             emptyText: "Search Address",             url: hereUrl,             locationField: "lonlat",             displayField: "name",             valueField: "bounds",             queryParam: "searchtext",             forceSelection: true,             store: hereStore         });              me.callParent(arguments);     }      });


title: How to use Leaflet Marker Cluster on markers from Leaflet-layerJSON (Overpass API)

text:
I'm having troubles getting the Leaflet Marker Cluster working with Leaflet- layerJSON. I'm trying to make a Leaflet map with data from OpenStreetMap. The markers are showing on the map, but they are not clustering. It would be a great help, if you could point out my mistakes in the following code:                    <meta charset="UTF-8"/>          <link href='https://api.tiles.mapbox.com/mapbox.js/v2.0.0/mapbox.css' rel='stylesheet' />                    <link href='https://api.tiles.mapbox.com/mapbox.js/plugins/leaflet-markercluster/v0.4.0/MarkerCluster.css' rel='stylesheet' />     <link href='https://api.tiles.mapbox.com/mapbox.js/plugins/leaflet-markercluster/v0.4.0/MarkerCluster.Default.css' rel='stylesheet' />

                     I have tried to find answers to my question in other posts, but without any luck. I'm new to webmapping and stackexchange and apologize if my question is unclear.


title: QGIS: Plot a layer given with Longitudes [0 360] to [-180 180]

text: I have a world map in correct Longitudes from -180° to 180°. Additionally, I have a shapefile of data from South America, which shows at 290° instead of -70°. How can I reproject it using QGIS? ![World map and data plotted at the wrong place](http://i.stack.imgur.com/EhoWE.png)


title: Can I convert ESRI Grid (DEM) raster to USGS DEM?

text: There is a tool for DEM to Raster in ArcGIS but I need to convert my ArcGIS DEM rasters (mostly LiDAR) into USGS DEM. The idea is to import this DEM into Sketchup and create a "sandbox" TIN surface and for this SketchUp needs the data to be in this particular format. In the past I mostly used the Sketchup Plugin to export TIN from ArcGIS to Sketchup which **worked great** but since the plugin no longer works in version 10. I am trying to create a new workflow but I am not having much success. The "Collada method" does not work well for complex surfaces. It requires to first create a multipatch and there is the first issue as described here. Besides, it is painfully slow if the conversion even succeeds at all. There is another option described here but I have not tried it. Can anyone comment on how well the Cloud v8 plugin mentioned there works?


title: Create Thiessen Polygons (Analysis) by defining an external feature boundary

text: Although this issue has been addressed in this post some months ago: Thiessen Polygon delineation within a Feature, the appropriate answer has not been given. So, new ideas are welcome. My question is how can we use the Create Thiessen Polygons (Analysis) tool available with ArcInfo license by defining an external feature boundary (e.g. study area boundaries defined by a polygon) for the created thiessen polygons. As I have seen this is not possible via the Environments...options in that tool. Is there any available code in VBA for doing this? or How can I run this tool within a VBA code? Thanks Demetris EDIT The comments indicate that ArcGIS 9.2 and 9.3 is being used, the tag indicates version 10, the first post is for that version not previous ones.


title: Auto Balancing in EarthWork

text: I am not sure whether this question belongs here. I think it belongs more to Geo-technical, but anyway, I just give it a shot. During earthwork development, earth terrain has to be cut ( if the terrain is higher than platform) or fill ( if the terrain is lower than platform) in order for the construction of **platforms** such as roads, housing area, pond and whatnot. Usually, one would have to make sure that the cut and fill within a development area is roughly the same, in order to reduce the import or export of the soil. The process of calculating the height of the platforms with relative to the terrain so that the cut/fill is roughly the same is called **_auto- balancing_**. How do engineers ( or earthwork software) do the auto balancing in practice?


title: In OpenLayers, is there any way to convert polygons to tiles?

text: I have a OpenLayer map which by enabling a feature by user, a very huge table of polygons (lots of small areas, with information and metadata for each) will appear on the map, which rendering this big amount of data is very slow and zooming in and out is kind of miserable. Do you know any solution to this problem? I was thinking of converting these polygons to tiles, which of course if possible, will not show any information for areas. But when the user reaches a certain zoom level, map starts showing information for areas, in case that area is being clicked by user. What I mean is behaving like tiles before a certain zoom level, and after that, behaving like polygons. I will appreciate any suggestion to the problem itself, or to my (maybe crazy) proposed idea.


title: Is there a standard vertical datum for meters below the WGS84 ellipsoid, or below MSL?

text: I have data I am publishing a WMS layer representing data under the surface of the water, say temperature at 100m depth. What is the accepted way to represent this in WMS GetCapabilities document? Do I define the dimension like:               <Dimension name="elevation"                 units="EPSG:5030"                 unitSymbol="m">0,50,100</Dimension>      or               <Dimension name="elevation"                 units="EPSG:5030"                 unitSymbol="m">0,-50,-100</Dimension>      or something else or should a custom dimension be used in this case? ~~EPSG:5030 is a length parameter, it doesn't indicate that it's related to any projection or that it might indicate distance from mean sea level for that matter. seehttp://trac.osgeo.org/geotiff/wiki/VerticalCS~~   EPSG:5030 represents height _above_ the WGS85 ellipsoid. I found a reference to EPSG:5715, but it doesn't say much other than it's 'msl depth' - does this mean values below msl should be positive?


title: which Coordinate system should I use for 3D analysis?

text: I am creating 3D surface below the earth's surface using ArcGIS and then I need to find distances to the plane from the surface. Which coordinate system shall I use? 3D interfaces like IProximity3D and IArea3D look like they work only with Projected coordinate systems which are 2-dimensional. I am new to this area and need direction to decide and use.


title: Is there an opensource software to route openstreetMap data based on database?

text: I have a project that consists in aggregate traffic information on OpenStreetMap for routing proposes. My first idea is to use the OSM data from the region on a database, so I can aggregate the traffic information on the table, probably adjusting the cost column. The problem is that all the routing software that I have already downloaded and tested creates a bin or obj file from the OSM instead of working with databases, so does anyone know an opensource software that creates and draw routes on OpenStreemaps based on database?


title: How to speed up the process of vector data visualization for a web application

text: I have a GIS web application and vector data that is going from the model to the web application. For developing the User Interface I am using JavaScript and OpenLayers. But sometimes the amount of geospatial data that I need to simultaneously visualize on a map is very big (10 000+ objects) and the browser cannot respond to all user requests. So, how can I speed up the process of visualization **(both model and UI part)** of large amounts of geodata?


title: Can we access the ArcMap Image Analysis tools using Python?

text: Esri recently added a new Image Analysis window to ArcMap - it's found under > Windows in the main toolbar: ![enter image description here](http://i.stack.imgur.com/ByDdO.png) > The Image Analysis window's Processing section simplifies the experience of > applying complex processing and analysis techniques to layers of image and > raster data in ArcMap. The tools provide one-click options to apply > processes, such as clipping, normalized difference vegetation index (NDVI) > creation, mosaicking, and exporting ![enter image description here](http://i.stack.imgur.com/hN6P9.png) **Can we access these tools using Python?** Specifically, I want to automate the process of applying a Clip, then adding Functions to the clip, before exporting to a new layer.


title: Formatting tags doesn't seem to work?

text: I am trying to do a simple format of a text box in Arcmap. I went into the 'About Formating Text' and copied a tag to change color and when I apply the setting, the text doesn't change color and the tags appear in the text like this: ![enter image description here](http://i.stack.imgur.com/eM6UA.jpg) Can anyone suggest what I am doing wrong? Thanks, Mike


title: location of a dragged and dropped marker using leaflet

text: I am using the leaflet API with OSM. I need to implement a functionality where on dragging the marker and dropping it on user preferred location, I need to get its latitude and longitude and save its location to database. For now on double clicking the marker I am able to drag the marker. This line helped me:               draggable:true      Now when I release the marker I want to find the latitude and longitude of the location or invoking a function on dragging the marker, returning the `lat,lng` values.


title: Finding all neighbors for all polygons in Spatialite

text: Excuse me if this is a basic question but I'm new to Spatialite... I have a dataset with zip code polygons and I would like to create a list where for each zip code we have listed all the neighboring (i.e. adjacent) zip codes. If I do the following it works but only for 1 zip code. What I would like to do is to repeat this automatically for all the zip codes in the data (since I have thousands of codes).               SELECT p1.ROWID, p1.ZCTA5CE10 FROM "US_zcta5_2010" p1, "US_zcta5_2010" p2     WHERE p2.ROWID=1 AND     Touches (p1.Geometry, p2.Geometry)=1;      I'm using the spatialite-gui on a Mac.


title: arcpy points in polygon check

text: I have a list of points and a list of polygons both in a separate shape files. I am trying to write a ArcPy script to loop through the points, and then the polygons respectively and tell me when a point is within a polygon. I have tried both the contains and touches methods on polygons but haven't had any luck. Any help would be appreciated. So far I have the following:                   pointFeatureClass = "C:\\Data\\POINTS2.dbf"         pointFieldList = arcpy.ListFields(pointFeatureClass)         for field in pointFieldList:             curField = field.name              arcpy.AddMessage(curField)              ptRows = arcpy.SearchCursor("C:\\Data\\POINTS2.dbf", "", "", "CID;Shape","CID;Shape")          for ptRow in ptRows:              #Now loop through poly's testing for intersect             shpRows = arcpy.SearchCursor("C:\\Data\\counties.shp", "", "", "NAME;Shape","NAME;Shape")             for shpRow in shpRows:                  #arcpy.AddMessage("Checking for point "+ str(ptRow.CID) +" in: " + shpRow.NAME)                 if shpRow.Shape.touches(ptRow.Shape):                     arcpy.AddMessage("Got it: " + shpRow.NAME)


title: Importing external bookmark (.dat) files into many MXDs using arcpy?

text: Does anyone know any python script to use for arcview to import a bookmark (saved as a .dat) that was created and saved externally using ArcMap? I tried with this script but without a result:               import arcpy     from arcpy import env          env.workspace = r"F:\Projects\Roni\abaloa"     for mxdname in arcpy.ListFiles("*.mxd"):         print mxdname         mxd = arcpy.mapping.MapDocument(r"F:\Projects\Roni\abaloa\\" + mxdname)         df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]         arcpy.LoadSettings("F:\Projects\Roni\abaloa\gis\layers\b1.dat")         for bkmk in arcpy.mapping.ListBookmarks(mxd, 'Bookmark1', df):             df.extent = bkmk.extent      mxd.save()      del mxd        I've also searched all documentation and forums with no luck.


title: Link GPS Data Logger to a website

text: Is there a way of linking GPS data logger to a self made website for tracking purposes? I want to create a tracking website which uses the data from GPS logger and connect to Google maps to show the track through the website I made.


title: 95 percent volume contour

text: I generated, from point data, a kernel density map using GRASS, and I would like to identify the 95% volume contour in it. Is it possible in GRASS? Thank you


title: How do I modify the "Playlist" Storymap to include a basemap switching widget?

text: I'm fairly new to HTML and the ArcGIS Javascript API. I would like to modify the default "playlist" storymap template to include a basemap gallery widget. Can someone help break down the steps to do this? This is what the default template looks like: http://storymaps.esri.com/stories/2013/20towns/ Here is an example of the modified playlist template with the basemap gallery widget. http://fishandgame.idaho.gov/ifwis/maps/wma/


title: How to add hyperlinks & embedded images to TileMill / Mapbox Tooltips?

text: I have a shapefile I am styling in TileMill; which will be imported later into a Mapbox webmap as MBTiles. There is an attribute field for each feature in the shapefile that contains a link to an online image. Can I embed this image and add a clickable hyperlink within a tooltip? The formatting I have now in the tooltip dialog is: "Link: {{{hyperlink}}}" where {{{hyperlink}}} is the unique url of an image for each feature. The link will display in the tooltip, but it isn't 'clickable.' I would also like to embed the image within the tooltip if that is possible as well. Thanks


title: Projecting GCS Australian 1984 to WGS 1984 UTM (or similar) to get units of meters?

text: I am currently trying to run a kernel density analysis and need to change the co-ordinate system to a UTM measurement in order to set some meaningful biological constraints to my seabird data (eg, Search Area and cell size output). The best out come i am wanting is to set these in meters. Currently my data is point data (XY or .CSV) and set in GCS_Australian_1984. I have been trying to change GCS_Australian_1984 to WGS 1984 UTM zone 12s (not sure if this an appropriate conversion in terms of file comparability) and when i run the "project" tool in data management i keep getting an error that says "table name is invalid" and "The project method cannot do a datum transformation". My data is located in an area of W: 110 N: -6 S: -21.5 E: 131 (Timor Sea(Southern Hemishpere)) Is there any way or a recommendation, that i can transform my data into a format where kernel density analysis will give meters and not SQUARE_MAP_UNITS, this is not very meaningful to me. Ps This is a good summary of what i am trying to achieve (http://www.youtube.com/watch?v=P3GP-IiQ-0Q) at 5 minutes and 16 sec and using arcmap 10.1 SOLUTION For any one having similar issues with getting your kernel density values (Search area or cell size out put) into biological terms then here is how to do it. If you have point data add it as X and Y data under your "normal" co-ordinate system. File > add data > add XY data. Then right click on the layer that was just imported. A tab will appear, select data > Export data > a box will appear select the file directory and export as a shape file. then import this shape file back into arcmap 10.1 and change the co-ordinate system to UTM and save. this can be done in the data management tab in the arctool box and is under projections and transformations, select the project tool. then run your kernel density with the transformed shape file. The option area units should be "ungreyed" select your units and done. NOTE: XY data can not be transformed into UTMs it has to be changed into a shape file and reimported.


title: NAs not allowed error message in attempt to plot raster object in R

text: I have imported a clipped raster (performed in ArcGIS 10.1) from the National Landcover Dataset into R, but when I try to plot I get the following error message. Error in .readCellsGDAL(x, uniquecells, layers) : NAs are not allowed in subscripted assignments In addition: Warning message: In .readCells(x, cells, 1) : NAs introduced by coercion I am relatively new to using rasters in R, but I cannot locate the NAs. levels(nlcd) yields the following:               ID  objectID    Value   Count   Red Green   Blue    Land_Cover_Class    opacity     0   1   0   4   0   0   0   Unclassified    1     1   12  11  584837  0.2784314   0.4196078   0.627451    Open Water  1     2   22  21  2247379 0.8666667   0.7882353   0.7882353   Developed, Open Space   1     3   23  22  954664  0.8470588   0.5764706   0.5098039   Developed, Low Intensity    1     4   24  23  647735  0.9294118   0   0   Developed, Medium Intensity 1     5   25  24  333254  0.6666667   0   0   Developed, High Intensity   1     6   32  31  12173780    0.6980392   0.6784314   0.6392157   Barren Land 1     7   42  41  5640937 0.4078431   0.6666667   0.3882353   Deciduous Forest    1     8   43  42  6987282 0.1098039   0.3882353   0.1882353   Evergreen Forest    1     9   44  43  957021  0.7098039   0.7882353   0.5568627   Mixed Forest    1     10  53  52  56354805    0.8 0.7294118   0.4862745   Shrub/Scrub 1     11  72  71  9100560 0.8862745   0.8862745   0.7568627   Herbaceuous 1     12  82  81  477862  0.8588235   0.8470588   0.2352941   Hay/Pasture 1     13  83  82  453351  0.6666667   0.4392157   0.1568627   Cultivated Crops    1     14  91  90  2671979 0.7294118   0.8470588   0.9176471   Woody Wetlands  1     15  96  95  495808  0.4392157   0.6392157   0.7294118   Herbaceuous wetlands    1      The only issue I can see being the problem is the first row. Below is the code that I used.               library(rgdal) ## spatial     library(raster) ## spatial     library (rgeos) ## spatial     nlcd = raster("C:/Extract_img11.tif")     plot(nlcd)      I can't post Extract_img11.tif because it is too large. Sorry for the vagueness of this question. My guess is that the clipping created the problem, but I haven't ever seen this before. Thanks for any information.


title: OpenLayers: adding a loading indicator for WMS requests

text: I'm new to Openlayers and GIS in general. I have a WMS request followed by the addLayer method, and surrounded by the addition and deletion of a 'loading' class to an HTML element:               $(toggle).parent('div').first().addClass('loading');     Map.currentOverlay = new OpenLayers.Layer.WMS(name, Constants.webServices.wms, {         format: 'image/png',         tiled: false,         transparent:true       }, {         singleTile:true,         ratio:1,         getURL: Projection.getDefaultOverlayUrl        });     }     if (Map.currentOverlay != null) {       Map.currentOverlay.initialZoom = map.zoom;       Map.currentOverlay.isBaseLayer = false;       Map.currentOverlay.visibility = true;       Map.currentOverlay.events.on({moveend: Map.__onZoomUpdate});       map.addLayer(Map.currentOverlay);       $(toggle).parent('div').first().removeClass('loading');      My goal is to add a loading gif in a dialog window (with CSS) while the layer is not displayed yet. According to this post: http://stackoverflow.com/questions/9569847/when-does-openlayers-wms- physically-call-the-remote-server the request is done in the addLayer method. But it seems asynchronous: the class 'loading' is added and removed straight away, while the map is not yet displayed. Can I add a callback or something similar to addLayer ? If not, how do I do? NB: I didn't produce the above code, I just added the "loading" behaviour. Do not hesitate to tell me if something seems wrong with it.


title: Problems setting up eclipse (with PyDev) for ArcGIS 10

text: I'm trying to use eclipse with PyDev to program in python. In Windows 7 x64, I've added a PATH environmental variable (PATH C:\Python26\ArcGIS10.0\python.exe;) and what I thought were the relevant folders to PYTHONPATH (C:\Dropbox\CODE\Python\; C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy) ![alt text](http://i.stack.imgur.com/NNd7H.png) However, when I type the command:               import arcpy      in eclipse I get an error message. However, when I try importing ArcPy via IDLE I have no problems. Is there some way that I should be including arcpy in the build path that I am missing?


title: Getting ArcGIS.com-created map to show up in sample app

text: I've been experimenting with the My First iOS app example, trying to get it to display a Bing map layer I created on ArcGIS.com. However, the layer doesn't display. In the sample app's code, when it says:               NSURL *mapUrl = [NSURL URLWithString:@"http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"];      AGSTiledMapServiceLayer *tiledLyr = [AGSTiledMapServiceLayer tiledMapServiceLayerWithURL:mapUrl];      [self.mapView addMapLayer:tiledLyr withName:@"Tiled Layer"];      what am I supposed to replace the mapUrl from? The bit.ly link created through the share function on ArcGIS.com doesn't work, since it links directly to a map viewer displaying the map. I'd like to get the map directly to put into my sample iOS app.


title: Why is Layer missing from Export to PDF?

text: I am using the arcpy module to loop through a number of map features (community areas), zoom to them, title the map with the name of the community area, and export to pdf. The trouble is that while the title, legend, and line features export, a symbolized polygon layer doesn't. It still shows up in the legend, but the map doesn't show it. Trying a couple of sanity checks, I have   1. confirmed that that layer in question is visible in the ArcMap map that I'm starting with,    2. tried making sure that that layer is on top,   3. instructed my code to make that layer the only visible one before exporting, and   4. had Python list the set of layers that are visible, right before exporting. None of those tests give me anything better, or give me any insight to the problem. My code is below. Since I'm also a relative newbie to coding with arcpy, any suggestions that you have for solutions or more elegant coding would be very welcome! -Nick My code is:               import arcpy, os     myPath    = r"C:\Users\<myname>\Documents\Heat Maps\\"     myMapPath = r"C:\Users\<myname>\Documents\Heat Maps\Map Output\\"          arcpy.env.workspace = myPath     myMap = arcpy.mapping.MapDocument(myPath + "Heat Maps of Need - v4-0 - Added Streets for Zoomed in Maps.mxd")          AllLayers = arcpy.mapping.ListLayers(myMap)     df = arcpy.mapping.ListDataFrames(myMap)[0]     lyrList = ["Needs Index"]     nComm = 77          # Clear visibility of all layers except for the community areas layer          for lyr in AllLayers:         if lyr.name = "Needs Index":             print "Turning on Layer " + lyr.name             lyr.visible = True         else:             lyr.visible = False          # Loop through specified layers and export to pdf          if os.path.exists(myPath + "Heat Maps of Need Indicators.pdf"):         os.remove(myPath + "Heat Maps of Need Indicators.pdf")     PDFdoc = arcpy.mapping.PDFDocumentCreate(myPath + "Heat Maps of Need Indicators.pdf")          for l in lyrList:              # Turn on proper heat layer         for lyr in AllLayers:             if lyr.name == l:                 lyr.visible = True         arcpy.RefreshActiveView()              # Check which layers are visible         for lyr in AllLayers:             if True == lyr.visible:                 print "Layer " + lyr.name + " is visible"              # Loop through all community areas         for c in range(1, nComm + 1): # Add 1, since the top value of the specified range is not included                  # Run tasks specific to the community area in question             for lyr in AllLayers:                 if lyr.name == "Community Areas":                          # Get name of community area                     SqlWhere = " \"AREA_NUMBE\" = \'" + str(c) + "\'"                          CommunityRowObj = arcpy.SearchCursor(lyr, SqlWhere)                     for row in arcpy.SearchCursor(lyr, SqlWhere):                         CommunityName = row.getValue("COMMUNITY")                         CommunityName = CommunityName.title()                          # Zoom in to community area                     arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", SqlWhere)                     df.zoomToSelectedFeatures()                     arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")                  for x in arcpy.mapping.ListLayoutElements(myMap):                 if x.name == "Title":                     x.text = "Heat Map of Needs Indicators\n" + l + " - " + CommunityName                  arcpy.RefreshActiveView()                  # Export the layer to PDF             arcpy.mapping.ExportToPDF(myMap, myMapPath + "Needs Heat Map - " + l + " - " + CommunityName + ".pdf")             PDFdoc.appendPages(myMapPath + "Needs Heat Map - " + l  + " - " + CommunityName + ".pdf")              # Turn off current Needs Layer to clear it from the next one         for lyr in AllLayers:             if lyr.name == l:                 lyr.visible = False          arcpy.RefreshActiveView()     arcpy.RefreshTOC()     PDFdoc.saveAndClose()          del PDFdoc     del AllLayers


title: Can you point me to an example of a modern-day John Snow?

text: It's became pretty much a custom to start spatial epidemiology or medical geography textbook / lecture with an example of John Snow's cholera investigation in 1854 Soho. ![alt text](http://i.stack.imgur.com/p6qTf.jpg) He's definitely still remembered in the current literature and even appeared on top answer to a recent question on this site. Could you give some recent examples where mapping, GIS or sptaial analysis provided significant contribution to understanding and explaining an epidemiological issue?


title: What tools for regionalization/redistricting can you recommend?

text: I am looking for tools for regionalization / districting in GIS. So far I have been able to identify   * Districting ArcGIS extension (now also online)   * REDCAP (currently limited to 50k input polygons)   * iRedistrict (not yet released)   * R method Is there anything else available out there? Especially something that will have no problems in dealing with quite large dataset [~ 1.3 million polygons]?


title: What type of calculations must be done for traffic emission for buffers in land use regression?

text: I have a raster layer that shows traffic emission distribution by kernel density method. Traffic emission will be used as independent variable in land use regression. What type of calculations must be done for traffic emission for buffers in land use regression like land use areas and road lengths. How can we build a model for this using ModelBuilder?


title: Interpolate/disaggregate point values to an existing raster grid

text: I have several thousands of ASCII files whose names have the following format: data_"lat"_"long".mon (e.g. data_41.21875_-116.21875.mon). Each file represents one point in WGS84 projection and the resolution of the gridded data is 1/16th degree (~ 6-7 km) Sample data file:                   Year Month Data1 Data2 Data3 Data4         1915 10 25.07710   41.88800   69.79910    0.00000         1915 11 20.32660   41.14250   66.66640    0.00000         1915 12 29.64930   65.21980   63.92670    2.29000         1916 1  24.38410   78.68940   73.11800   21.30610         1916 2  21.93770   72.59410   78.67960   77.18490      I got another raster grid (Lambert Conformal Conic projection) which has the cell size of about 400m (1320 ft) ![enter image description here](http://i.stack.imgur.com/NJ8s0.jpg) **This is what I'd like to do**   1. Read the names of all data files in the folder and extract the coordinates of all the points   2. Read the Year, Month & Data3 (column 5)   3. Loop through all the timesteps, Interpolate/disaggregate the point values of Data3 to the raster grid   4. Write all the grids to a file one after another (Raster to ASCII) The first two steps I can do it relatively easy in R or Python. However from the 3rd step, I can't find any tool that allows me to do interpolation from points to existing raster grid in ArcGIS. **Could anyone suggest any R or Python codes/packages that could deal with iteration of interpolation/disaggreation?** I have some experience with R but pretty novice to Python (but I'd love to gain more experience with it). Another concern is that, in terms of accuracy, **I'm not sure which projection (WGS84 or Lambert Conformal Conic) I should use for interpolation** Thank you very much for looking! Any suggestion would be greatly appreciated!


title: How to export figure like Post In Action demonstrate?

text: I am reading the book "PostGIS In Action" and now I follow the listing 2.4 here               SELECT AddGeometryColumn ('public','my_geometries',     'my_circular_strings',-1,'CIRCULARSTRING',2);     INSERT INTO my_geometries(name,my_circular_strings)      VALUES ('Circle',     ST_GeomFromText('CIRCULARSTRING(0 0,2 0, 2 2, 0 2, 0 0)')),      ('Half Circle',     ST_GeomFromText('CIRCULARSTRING(2.5 2.5,4.5 2.5, 4.5 4.5)')),      ('Several Arcs',     ST_GeomFromText('CIRCULARSTRING(5 5,6 6,4 8, 7 9, 9.5 9.5,     11 12, 12 12)'));      And out put as               SELECT name,ST_NPoints(my_circular_strings) As cnpoints,      ST_NPoints(ST_CurveToLine(my_circular_strings)) As lnpoints      FROM my_geometries      WHERE my_circular_strings IS NOT NULL;      and postgresql give me a table                                       name                        | cnpoints | lnpoints      ----------------------------------------------------+----------+----------      Circle                                             |        5 |      129      Half_Circle                                        |        3 |       65      Several Arcs                                       |        7 |      113     (3 rows)      My question is the book also give me a figure as result like below ![book's example](http://i.stack.imgur.com/zpu9s.jpg) I want to make a picture like that ,how to do it??


title: How to use GeoServer securely in production with Tomcat Policies?

text: I want to use GeoServer (first testing, then production) but we have to use Tomcat policies for security reasons. I have already tried to deploy GeoServer with the security manager enabled, but GeoServer will only work if granting explicit access to the root (/) of the server. Is this so? I didn't find anything in the docs. Does anyone know how to get GeoServer working with Tomcat policies?


title: Adding geometry to rows with missing X,Y inArcGIS

text: I have a table with X, Y data that contains many rows with still missing coordinates. Now I would like to edit it in ArcGIS 10 and simply add the geometry for selected lines after identifying the correct location on a base map. This seems so simple and basic, but I can't get it to work. Always, adding a point in ArcGIS will just create a new row. With more than 1000 yet unidentified spots, this is really cumbersome.


title: If i use Cartaro to manage my geospatial data handling website, will it be possible to access other Drupal distributions?

text: If i use **Cartaro** to manage my geospatial data handling website, will it be possible to use or integrate other Drupal distributions or modules to customize other sections of my site?


title: Why does QGIS appear to ignore the bbox addition to my WFS uri?

text: I'm pulling a WFS layer from the LINZ data service (1.0.0 and 1.1.0 complient) using a python script as follows               from PyQt4.QtCore import *     from PyQt4.QtGui import *     from qgis.core import *     sourceUri = "http://wfs.data.linz.govt.nz/<MYPRIVETAPIKEY>/wfs?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=v:x245&BBOX=1656381,5303502,1657111,5304135"     vlayer = QgsVectorLayer(sourceUri, "bbox test", "WFS")     QgsMapLayerRegistry.instance().addMapLayer(vlayer)      Except instead of pulling only two objects QGIS imports the whole layer. It's as though its dropping the bbox clause that I appended. When I copy that uri into a simple browser I DO get the bbox clipped layer (two objects) as expected. Altering the uri for other WFS versions still give problems               http://wfs.data.linz.govt.nz/<MYPRIVETAPIKEY>/wfs?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&SRSNAME=EPSG:4326&TYPENAME=v:x308&BBOX=1656381,5303502,1657111,5304135,EPSG:2193      Also gives the correctly clipped layer in a browser but gives an empty set in QGIS. NB no errors are logged in either case


title: Find area that it's in one GeoJSON file but not the other

text: I have two GeoJSON files, A and B. The features in A contain everything in B, and more (they're contour lines, and A is a lower elevation). How can I get the difference between the two files? As in, I want a file with a feature or features that covers all the land that is in A but is NOT in B. It doesn't really matter whether the result is just one feature, or a hodgepodge of features, as long as it covers the correct area. But I have no idea how to calculate that difference in QGIS or by any other method.


title: How to use relative path in mac for open internal and external file?

text: I would like that my project would be portable, I mean that I can share it without lose all the path I have in my project. This seems not so easy and many times I have to re-set all the path when I give the project, with all the folders with rasters and shape file into to a colleague. ![enter image description here](http://i.stack.imgur.com/xKZRh.png) So I would like to know:   1. can I set by default only relative path to the folder in which my project is (for opening raster and shape file that QGIS can handle) and this remain relative even if I share the whole project to a colleague with a different operating system.   2. For opening an external file with `action` in `properties->Action` I use `Open file` with this action `[% "Ref" %]` where `Ref` is the attribute where I write the file path. I supposed that writing the path with `./` before the path (e.g. `./RASTER/test.jpg`) could be use to create a relative path but this seems not working. Why? And this however would work only in Unix like system! i'm using QGIS 1.8. Thanks a lot!


title: ListFields (arcpy) not reading from layer

text: I have a very simple python script that opens an mxd from a folder and then loops through each layer, and finds a ObjectID field and changes the alias name.               #setup     import arcpy     from arcpy import env          #Go into folder where mxd's are     env.workspace = r"D:\Temp\Test_Mxd"     #for each mxd     for mxd in arcpy.ListFiles("*.mxd"):         print "Working on map document: " + mxd         mapdoc = arcpy.mapping.MapDocument("D:\\Temp\\Test_Mxd\\" + mxd)         #   for each dataframe         df = arcpy.mapping.ListDataFrames(mapdoc,"Layers")[0]         ##layers = arcpy.mapping.ListLayers(mapdoc,"*",df):         layers = arcpy.mapping.ListLayers(mapdoc,"*",df)              #   for each layer         for layer in layers:             print "Working on layer: " + layer.name             **fieldList = arcpy.ListFields(layer)**             for field in fieldList:                 if field.name == "OBJECTID":             # Set the new alias name                     field.aliasName == "IGNORE::OBJECTID"         mapdoc.save()     del mxd         my problem is on the line `fieldList = arcpy.ListFields(layer)`, I get a error in the interpreter (pythonwin) which says that the layer does not exist. Has anyone got any idea on where I am going wrong? Thanks.


title: How do you manage existing GIS data and allow for integration with an Asset management system?

text: I am using ESRI based GIS software, Postgresql/PostGIS/ArcSDE DB and we have a Mincom Ellipse asset management system. Currently our all our spatial information regarding asset sites have been recorded as points, this has been versatile for the fact that it accommodates mapping at various scales. Now that we are integrating our Asset Management System with our GIS database the asset management guys want the GIS features to reflect the structure e.g A building footprint as a polygon instead of a point. My question in terms of the spatial data management. Should I be maintaining two sets of data? One for the asset representations and then one for various mapping tasks? Thanks DB


title: In Postgis 1.5, how can I split lines in points while retaining the id?

text: # Hypothesis: Having 2 linestring tables:               pipelines (id integer, street varchar(50), shape geometry)      streets (id integer, name varchar(50), shape geometry)      Each pipeline passes/intersects several streets. ## Need I need to:   1. Fill the street field in the pipelines automatically and   2. I need to know to which street does the pipeline "belong the most". I thought about densifiying the vertices of the each pipeline at a maximum 5 meter distance from each other, and do the analysis on them and then statistically find out the answer. **How In Postgis can I create a table that presents the geometry for each point of the pipelines along with their corresponding id?**


title: Is there an existing way to do linear referencing in ArcGIS JavaScript API 2.1?

text: First of all, let me specify that this is for the JavaScript API 2.1 or later only, not ADF/Google/Virtual Earth/Silverlight/Flex or any of the other extensions. We've been creating geometry from a road geometry service to display travel routes on our map. We also have a layer of known routes. We're working on trying to speed things up by using the known routes in a JSAPI map, and calling the road geometry service only when we need to update or add routes. Using the road geometry service, we can draw the route so that when a vehicle is shown on the route, the road already traveled can be red, and the remaining route as green. Has anybody tried to render lines with different symbols on different parts of the line?