title: parsing an osm-file with PHP simpleXML :: getting a few more tags

text: I am new to PHP's SimpleXML. **upate:** The original version of this question was derived from here: OSM Data parsing to get the nodes with child I am thankful that hakre offered a great example in the comments that makes a overwhelming starting point for my project. Below I have added my own answer to the question, how to refine the code to ad more tags. I want to filter the data to get the nodes with special category. Here is sample of the OSM data I want to get the whole schools within an area. The first script runs well - but now I want to refine the search and add more tags. I want to store all into MySQL. So we need to make some XML parsing with PHP: The following is a little OSM Overpass API example with PHP SimpleXM               <?php     /**      * OSM Overpass API with PHP SimpleXML / XPath      *      * PHP Version: 5.4 - Can be back-ported to 5.3 by using 5.3 Array-Syntax (not PHP 5.4's square brackets)      */      //     // 1.) Query an OSM Overpass API Endpoint     //          $query = 'node       ["addr:postcode"~"RM12"]       (51.5557914,0.2118915,51.5673083,0.2369398)->.point;       (       node       (around.point:100000)       ["amenity"~"school"];       way       (around.point:100000)       ["amenity"~"school"];       );       out;';          $context = stream_context_create(['http' => [         'method'  => 'POST',         'header' => ['Content-Type: application/x-www-form-urlencoded'],         'content' => 'data=' . urlencode($query),     ]]);          # please do not stress this service, this example is for demonstration purposes only.     $endpoint = 'http://overpass-api.de/api/interpreter';     libxml_set_streams_context($context);     $start = microtime(true);          $result = simplexml_load_file($endpoint);     printf("Query returned %2\$d node(s) and took %1\$.5f seconds.\n\n", microtime(true) - $start, count($result->node));          //     // 2.) Work with the XML Result     //          # get all school nodes with xpath     $xpath = '//node[tag[@k = "amenity" and @v = "school"]]';     $schools = $result->xpath($xpath);     printf("%d School(s) found:\n", count($schools));     foreach ($schools as $index => $school)     {         # Get the name of the school (if any), again with xpath         list($name) = $school->xpath('tag[@k = "name"]/@v') + ['(unnamed)'];         printf("#%02d: ID:%' -10s  [%s,%s]  %s\n", $index, $school['id'], $school['lat'], $school['lon'], $name);     }          ?>      how to get more out of it.. at least the adress and the website should this be added in this part!?               # get all school nodes with xpath     $xpath = '//node[tag[@k = "amenity" and @v = "school"]]';     $schools = $result->xpath($xpath);     printf("%d School(s) found:\n", count($schools));     foreach ($schools as $index => $school)     {         # Get the name of the school (if any), again with xpath         list($name) = $school->xpath('tag[@k = "name"]/@v') + ['(unnamed)'];         printf("#%02d: ID:%' -10s  [%s,%s]  %s\n", $index, $school['id'], $school['lat'], $school['lon'], $name);     }      since i am learning - i break down the code into pieces... For you the second part is more interesting. That is querying the XML data we have already. This is most easily done with xpath, the used PHP XML library is based on libxml which supports XPath 1.0 which covers the various querying needs very well. The following example lists all _schools_ and tries to obtain their names as well. I have not covered translations yet because my sample data didn't have those, but you can also look for all kind of names including translations and just prefer a specific one):               //     // 2.) Work with the XML Result     //          # get all school nodes with xpath     $xpath = '//node[tag[@k = "amenity" and @v = "school"]]';     $schools = $result->xpath($xpath);     printf("%d School(s) found:\n", count($schools));     foreach ($schools as $index => $school)     {         # Get the name of the school (if any), again with xpath         list($name) = $school->xpath('tag[@k = "name"]/@v') + ['(unnamed)'];         printf("#%02d: ID:%' -10s  [%s,%s]  %s\n", $index, $school['id'], $school['lat'], $school['lon'], $name);     }      The key point here are the xpath queries. Two are used, the first one to get the nodes that have certain tags. I think this is the most interesting one for you:               //node[tag[@k = "amenity" and @v = "school"]]      This line says: Give me all _node_ elements that have a _tag_ element inside which has the _k_ attribute value _"amenity"_ and the _v_ attribute value _"school"_. This is the condition you have to filter out those _nodes_ that are tagged with amenity school. Further on xpath is used again, now relative to those school nodes to see if there is a name and if so to fetch it:               tag[@k = "name"]/@v'      This line says: Relative to the current node, give me the _v_ attribute from a _tag_ element that as the _k_ attribute value _"name"_. As you can see, some parts are again similar to the line before. I think you can both adopt them to your needs. Because not all _school_ nodes have a name, a default string is provided for display purposes by adding it to the (then empty) result array:               list($name) = $school->xpath('tag[@k = "name"]/@v') + ['(unnamed)'];                                                         ^^^^^^^^^^^^^^^                                                     Provide Default Value      So here some of the results for that code-example:               Query returned 907 node(s) and took 1.10735 seconds.     more than 2000 School(s) found:     #00: ID:332534486   [39.5017565,16.2721899]  Scuola Primaria     #01: ID:1428094278  [39.3320912,16.1862820]  (unnamed)     #02: ID:1822746784  [38.9075566,16.5776597]  (unnamed)     #03: ID:1822755951  [38.9120272,16.5713431]  (unnamed)     #04: ID:1903859699  [38.6830409,16.5522243]  Liceo Scientifico Statale A. Guarasci     #05: ID:2002566438  [39.1347698,16.0736924]  (unnamed)     #06: ID:2056891127  [39.4106679,16.8254844]  (unnamed)     #07: ID:2056892999  [39.4124687,16.8286119]  (unnamed)     #08: ID:2272010226  [39.4481717,16.2894353]  SCUOLA DELL'INFANZIA SAN FRANCESCO     #09: ID:2272017152  [39.4502366,16.2807664]  SCUOLA MEDIA       and now i try to figure out how i can enter more xpath queries at the above mentioned code and get out even more important data - see here [url=http://wiki.openstreetmap.org/wiki/Key:contact]Key:contact \- OpenStreetMap Wiki[/url] > >     contact:phone >     contact:fax   >     contact:website >     contact:email >   i will digg into all documents and come back later the weekend... and report all the findings well - i think that i need to extend the xpath requests within the loop where xpath is used again, now relative to those school nodes to see if there is a name and if so to fetch it:               tag[@k = "name"]/@v'     tag[@k = "contact:website"]/@v'     tag[@k = "contact:email"]/@v'      What do you say...? many many thanks to you for all your help. You are very supportive. Keep up the superb work here.... look forward to any and all help


title: How to represent intensity of flow by multiple lines in QGIS Flowmapper

text: New to the software, it fits my needs perfectly. I am mapping the flow of a certain product through various stages of production and dissemination. However, though I do like representing volume by proportional width, I was wondering if I could represent this with a proportional number of flow lines instead. The reason is mainly aesthetic, as I say, I am content without it like this, but the fact we are dealing with multiple material products means that I feel multiple lines would be better than the thickness of lines. If anyone knows of any solution, please help, or otherwise let me down now so I don't waste too much effort on this. Lastly, any way to manipulate the flow so that it is curved rather than straight. Thanks in advance.


title: Points on map disappear when apply certain column

text: I am a linguist who studied some variables in New Zealand English. Now I wanted to map the proportion of certain variables onto a map. It has worked fine so war and I've been able to map quite a few, but now I seem to be stuck. I created a spreadsheet with the first column being the place names, the second and third columns are the coordinates, and not the forth column is the proportions of a certain variable. I saved it as a text file and then wanted to map it the usual way (add delimited text layer, then style, gratuated etc.. ) but when I want to apply it, the points on the map disappear.. why's that? I've also saved the text file as cvs and tried it that way, but the same problem arises.. Any ideas? I can't understand why it worked with all the other variables, but just this one doesn't.. I use qGIS 2.2.0-Valmiera Thanks a lot for your help


title: Grouping of NDVI values in erdas Imagine

text: I am trying to extract NDVI values from MODIS image i have image of my study area, when i look at the attribute values of image in Erdas i get the value ranging from -32768 to 9820, I want to create new groups like -32768 to -31768 -- 1(after recoding) Could somebody help in making me understand how to go do this in Erdas


title: DEM to point file & interpolation

text: I have a DEM which has noisy parts to it. I am trying to (eventually) calculate the volume. However when I do cut/fill, there are obviously white sections where it hasn't been able to calculate the noisy parts. To get around this I was going to interpolate my DEM with one of the interpolate tools (IDW/Kriging/NN/Spline etc) Problem: It has to be point file, so I was converting raster to point file when error 000870 appear stating that I already have an output file with the same name. I definitely don't and I have turned overwrite on just in case. I was then going to convert back convert back to raster to do the volume. If anybody can think of a simpler way or can solve my mysterious file problem it would be much appreciated. Using ArcGIS 10. All extensions possible, cannot download additional programs due to being on a university computer and not having admin rights. ![Noisy DEM](http://i.stack.imgur.com/Bl0lE.png)


title: How to extract impervious areas from Landsat data?

text: I'm working with multi-temporal datasets from Landsat (Landsat 8 & Landsat 7 & Landsat 5) and I tried to extract impervious areas, more exactly urban areas, but I can't obtained results very close to reality (accuracy was approximately 80 % by using reference point in Google Earth). I used supervised classification: SAM & Maximum Likelihood, then logical classification (between SAVI + MNDWI + NDBI) but I can't obtained very good results. I used bands: Red, Green, Blue, NIR, MIR, SWIR 2. How I can extract urban areas from different Landsat image? What methods would you me recommend ?


title: Date-format conversion

text: I'm using Time Manager and QGIS 1.8.0 and I am stuck at date-format conversion. Time Manager needs Gregorian format `YYYY-MM-DD HH:MM:SS.ssssss` in text field and I have thousands of records in "month/day/year hour:minute:second" format (Julian format). As I'm not a GIS expert I need some directions to make this conversion correctly and be happy with Time Manager presentations.


title: How to calculate distances from points to multiple features at the same time?

text: I have a point layer and I want to calculate the distance from each point to a number of line features (rivers, roads etc.). Using the "Near" tool I can do this for one feature at a time (eg rivers) and I get 2 new fields in the attribute table indicating the id of the nearest feature and the distance. Fine. However, trying to find distance to the next features (eg roads) simply overwrites the previous ones (for rivers). The problem is I want to get distances to multiple features and get them all added separately to the attribute table. I tried to change the field name to see if this helped but I can only change the "Alias" so it makes no difference. I am working in Arcmap 10.0. Thanks.


title: Tileserver with time dimension

text: We are preparing a Tileserver to serve tiles with a time dimension that updates every day. I was wondering what could be the best approach for that or the ideas we had are a good way. To illustrate the problem: We have observations for every full hour over 10 hours. That would be then 10 timesteps if implemented with a timeslider. The observations change every day. Our idea was to make 10 layers which are accessible over a different layername. The layername would reflect the time of the observation(e.g. ISO 8601 timestring) or could be a simple number like 1, 2, 3, ... for every timestep with a corresponding meta.xml containing the date of the first layer and the interval between the timesteps. Of course it would be the best if the meta.xml would be something like a WMS GetCapabilities. Different tile serving softwares as MapProxy or Mapcache serve some kind of metadata.xml for TMS and WMTS GetCapabilities for WMTS but seem unable to insert a time dimension. How could we implement this? There was a post asking for implementations of the WMS time dimension for a WMTS but maybe in the meantime there where some developments. There is also a blog post on how the guys at nasa added a Dimension block to the WMTS GetCapabilites but I don't know which software they used for it. In the end we would prefer to use some existing standards over implementing our own structure. Is this a good way to go? Do you think there is a better solution?


title: How to add lat/lon data in ArcGis 10.1?

text: How do you add xy data to a DEM Raster in ArcGIS 10.1? I have done this before, and it has worked just fine. However, when I recently tried adding a MS Access database which has lat/lon coordinates for 10 sites, it did not work (I am currently using Microsoft Office 2013, I am not sure whether this has anything to do with it but when it worked before I was using Microsoft office 2010). I have tried many ways, for instance: 1) Saving the Access table as a .mdb and then on ArcCatalog making the OLE DB connection using Microsoft Jet 4.0. OlE Provider. Apparently the connection was a 'success, so I proceeded to create a feature class 'from XY' Table. When I add this to the DEM raster in ArcMap, it does not display the site locations. 2) I have also tried saving the Access table as a .accdb, then adding the database connection using 'Microsoft Office 12.00 Access Database Engine'. I typed in the correct pathway and when I click on 'Test Connection' the program crashes. 3) Finally I have tried adding xy data in ArcMap, but encounter the same problems - either it doesn't display or the program automatically shuts down. If anyone can please let me know whether I am missing a crucial step or something is actually wrong with the software, or whether there is anything else I can do, it would be greatly appreciated!


title: PostGIS Create Table with Geometry

text: I have written a query that will create a table and perform a GROUP BY from an existing spatial table, but the table does have the geometry. I am wondering where, and what, in this query I put to make the new table have spatial geometry (points). Thanks...               SELECT st_id, name, COUNT(st_id)         INTO climate.test_stats     FROM climate.lugdc_stations     GROUP BY st_id, name     ORDER BY count DESC;          ALTER TABLE climate.test_stats     ADD COLUMN pk_id serial PRIMARY KEY;


title: Displaying geojson data in leaflet

text:
I am trying to access and display some data in a geojson layer and would appreciate if someone can take a look at my code below and suggest what I need to do to make it work. The result of the console.log is:   Object {latlng: o.LatLng, layerPoint: o.Point, containerPoint: o.Point, originalEvent: MouseEvent, type: "click"…}               var wrldbndryData;          function wrldbndryStyle(feature)      {         return {             weight: 2,             opacity: 0.1,             color: 'black',             fillOpacity: 0.1         };     }          function highlightFeature(e)      {         var layer = e.target;              layer.setStyle(         {             weight: 9,             color: '#666',             fillOpacity: 0.5         });              if (!L.Browser.ie && !L.Browser.opera)          {             layer.bringToFront();         }     }          function resetHighlight(e)      {         wrldbndryData.resetStyle(e.target);     }          function showFeature(feature)  <-- *** THIS IS WHERE I AM TRYING TO ACCESS THE GEOJSON DATA *****     {         console.log(feature);          if (feature.properties && feature.properties.NAME)          {             $('#sidecol').html(' COUNTRY INFO </center' + '
' + feature.properties.NAME + '');         }     }          function eachFeature(feature, layer)      {         layer.on(         {             click: showFeature,             mouseover: highlightFeature,             mouseout: resetHighlight         });     }          wrldbndryData = L.geoJson(wrldbndry,     {         style:wrldbndryStyle,         onEachFeature:eachFeature     });          wrldbndryData.addTo(map);


title: How to use right-click to finish drawing a Polygon in OpenLayers?

text: I'm using OpenLayers 2.13 and the DrawFeature control. By Default, the Polygon is completed, when you double click. Some of our clients are requesting that we allow the mouse right-click, to finish the polygon. I tried adding a rightclick function to the `OpenLayers.Handler.Polygon` class, but it is never called. I also confirmed that the right click event is being registered in the activate method of the `OpenLayers.Handler` class. How do I override this behavior to allow the mouse right-click to finish the polygon?


title: Moving Versioned Data between different geodatabases

text: Can you please let me know how to migrate feature class which is registered as version and archiving is enabled on it From: Database Environment: ArcSDE 9.3, RHEL 5, Oracle 11.2.0.2 To Database Environment: ArcSDE 10.1, RHEL 6, Oracle 11.2.0.3.0 I used import feature class tool and copy tool to migrate these feature classes, but they succeded only in feature class data and no versioned/archiving stuff was found in destination database. Please let me know the best practice for this....


title: Display MapBender map in webpage

text: ~~I 've built a webpage and I want to display in it, a map that I have in MapBender. Can someone tell me how to proceed?~~ I resolved this issue. All i had to do is put an iframe in the page and point it to the MapBender link. Now, I want to get the coordinates (and other information) from a point on the map when I click on it and put them in some text fields that I have. Can anyone help?


title: Offset WGS84 linestring based on pixel values

text: I would like to split up some roads based on their driving direction. For this, I need to offset the roads 2 pixels to the right, in the direction of the linestring. I have the coordinates in WGS84 and I have a value that mentions the amount of meters per pixel. Currently I have this code which I found online:                       public static ILineString ShiftLine(ILineString input, double distance)         {                  IPoint point0 = input.StartPoint;             IPoint point1 = input.EndPoint;                  // Calculate the angle of the line.             double avgAngle = Math.Atan2(point1.Y - point0.Y, point1.X - point0.X);             avgAngle -= Math.PI/2;                  // Now convert back to an x/y offset.             double yOffset = Math.Sin(avgAngle)*distance;             double xOffset = Math.Cos(avgAngle)*distance;                  Coordinate[] shiftedCoords = new Coordinate[input.NumPoints];             for (int i = 0; i < shiftedCoords.Length; i++)             {                 IPoint p = input.GetPointN(i);                 shiftedCoords[i] = new Coordinate(p.X + xOffset, p.Y + yOffset);             }             return new LineString(shiftedCoords);         }      The offset from meters to pixels, is done like this:                           var lengthInMetersOfAdegree = GlobalMercator.LengthOfDegreeForLatitude(entry.Y1);                 // 1 degree = lengthinmetersofdegree                 var offset = metersPerPixelResolution / lengthInMetersOfAdegree;      The offset is then passed to the ShiftLine function. I suppose something goes wrong with the offset code, but I'm not sure. Anyone can help me with pointing out the error, or supplying code sample that works? Basically all I want is WGS84 input coordinates, shift them 2 pixels. I'm not an expert on this field, any help is welcome.


title: GRASS v.to.rast.attr returning null values

text: I am having a problem converting vectors to rasters in QGIS. As there is no tool where you can directly set the pixel size (instead of the height and width of the new file) in QGIS 1.8 (the version I have at the moment, I am unable to update to 2.0 due to IT issues at work), I am now trying to use the GRASS plugin v.to.rast.attr. I would like to be able to set the pixel size to 10 m (the file and mapset area are projected in UTM). When I try and import the file and then use v.to.rast.attr to convert it from raster to vector, all pixels return the value of -1e+300, when the original values of the shpfile range between 5 & 50. Does anyone have solutions to fix this problem or an alternate option for converting vectors to rasters where the pixel size can be set?


title: CartoDB: using the layer selector tool

text:
I want to set up a layer selector with CartoDB like what is found in Usage Examples here. I viewed the source code and changed it for the table I want to use. When I open the html file and click on a different layer, I get a 404 error. The URL for the error states:               ({"errors":["Error: could not fetch source tables: column \"undefined\" does not exist"]});      But I can see and query my table in CartoDB just fine. Below is my createSelector function. I'm not trying to do anything more than what is in the example I linked above. Does anyone have an idea on what I am missing? To be clear, I'm not sure the error above is related to a problem with my code, but it could be! I have specified my CartoDB user name and my json file. Thank you.

Collisions

Pedestrian 2013

Pedestrian 2012

BIKE 2013

Bike 2012


title: How to merge two 'incompatible' polygon layers in ArcGIS?

text: I have two layers of polygons with administrative boundaries. Although they come from the same data provider, they seem to have small discrepancies and do not 'fit' each other. How could I merge these two layers in ArcGIS 9.3 removing overlaps and gaps between them? Would it be possible to prioritize one layer (grey one in my case) and 'fit' the other one (purple) accordingly? ![alt text](http://i.stack.imgur.com/pxbsf.png)


title: Forcing routes to visit closest facility in ArcGIS

text: My dataset contains three parts: A) Incident location B) AED locations C) locations of available "first responders"(don't think of these like an ambulance or a fire department, but rather normal civilians) Calculating the shortest path from C to A or C to B is doable with the Network Analyst tools, but is there a way for me to calculate two routes for each person in C? I.e. I'd like to have both a direct route, C to A as well as the route via the AED, C -> B -> A, where B in this case is the "closest facility"(closest AED) for the individual C. I tried to work with barriers with negative value, but that didn't help much. Any pointers on how to create the two routes?


title: In ArcPy, how can I save a clipped raster dataset to memory (not to file)?

text: I've been trying to clip a larger raster dataset by a polygon extent, where the resulting clipped raster dataset is saved in-memory. So far, my attempts produce vague errors.               >>> arcpy.Clip_management(in_raster="/path/to/raster_dataset.tif", rectangle="899329.78438381 1065232.74648145 912509.334774774 1080416.77223298", out_raster="in_memory/clipped", in_template_dataset="/path/to/boundary.shp", clipping_geometry="ClippingGeometry")     Runtime error <class 'arcgisscripting.ExecuteError'>: ERROR 999999: Error executing function. Failed to execute (Clip).      Is there a way I can store a raster dataset in-memory? Documentation for the VB.NET API seem to indicate it is possible beginning in Arc 9.2, and I'm using Arc 10 on my workstation. > ArcGIS 9.2 allows to create an in-memory RasterDataset where the pixels are > stored in memory, this type of raster dataset has the charastics of a file > RasterDataset.


title: Draw arc in Geo map WGS84 coordinate system

text: im trying to draw an arc in OpenLayers, in non Geo map it works well and here is my code:               function calculateArcPointsNonGeo(center, radius,startAngle,endAngle,segments)      {                 var pointList =  [],                     Angle,                     x,                     y,                     point,                     i;                      var dAngle = segments + 1;                 for (i = 0; i < dAngle; i++) {                     Angle = startAngle - (startAngle - endAngle) * i / (dAngle - 1);                      x = center.X + radius * Math.cos(Angle * Math.PI / 180);                      y = center.Y + radius * Math.sin(Angle * Math.PI / 180);                     point = new MapPoint(x,y);                     pointList.push(point);                 }                      return pointList;             }      in addition there is a great Example that demonstrates it , but when im trying to do this on Geo map (WGS84 coordinate system) its start to acting weird, i assume that i need to convert degrees to radinas or something like this. please help


title: Geoserver projection recognize problem

text: I'm publishing shape files in geoserver through curl,my problem is that geoserver didn't recognizing the projection which is in ESRI WKT format(layer was publishing).I searched for this problem i found that ESRI wkt format don't have the "EPSG:4326" . What i have to do to recognize the ESRI wkt format projection?


title: Does the psql prompt have a working directory in the file system?

text: I am learning postgis. This tutorial says "Navigate to where you downloaded the shapefile." Does that mean navigate to this location in the `#psql` prompt or navigate to this location with the `$propmt`? In other words: does the psql shell have a working directory in the file system that you can change with `cd`, like in *nix? I am confused because the tutorial says to run `shp2pgsql` after navigating to that location -- but `shp2pgsql` is a postgres command, not a bashcommand. But the sql prompt lists tables, records etc -- which does not have to do with the file system. So does the sql prompt know about the file system or not?


title: Export map to PDF using Esri JavaScript API

text: I've searched a lot and tried many ways to export a web map to PDF. I need to export a web map and all active layers. Can anyone help me?


title: Python problem visualization .tif image

text: I am trying to read a hyperspectral image in Uint16 .tif format and visualize it. I use GDAL in Python to read the data and everything works perfectly when I use matplotlib to show one band of the data (The first image). The image however becomes completely corrupted when I want to concatenate(stack) 3 different selected bands and visualize them. Here is the code:               import numpy as np     import struct     from osgeo import gdal     from osgeo import ogr     from osgeo import osr     from osgeo import gdal_array     from osgeo.gdalconst import *     import matplotlib.pyplot as plt          cube = gdal.Open(filename)     bnd1 = cube.GetRasterBand(29)     bnd2 = cube.GetRasterBand(19)     bnd3 = cube.GetRasterBand(9)      Now I turn each band into a ndarray:               img1 = bnd1.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)     img2 = bnd2.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)     img3 = bnd3.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)          and stack them to have a 3 band image               img = np.dstack((img1,img2,img3))          f = plt.figure()     plt.imshow(img)         plt.show()      I would appreciate any suggestion on it here are the results: ![visualization of one band](http://i.stack.imgur.com/o9U0b.png) ![Visualisation of 3 bands](http://i.stack.imgur.com/6yDeF.jpg)


title: Geocode a table, then append more rows to the table

text: I'm starting on a project with ArcGIS 10 where I'll be doing quite a bit of geocoding addresses. I understand the basics of how to do that. But a frustration on an earlier similar project is that when the boss would come down the hall with just a few more addresses for me to add to the map (over and over). So I would have to repeatedly re-geocode the whole table each time I add a few. Thinking "there must be a better way", I came across the option on the geocoding dialog to create a dynamic feature class related to a table. This promises to automatically update the feature class when the data in the table changes. So I have tried this. I built a file geodatabase for this project, and imported some test data into it, and then geocoded the addresses. In ArcCatalog, I see my data table, the feature class, and a relationship class that links them. However, when I edit an address in the table, I don't see any auto-updating going on. And more importantly, when I append more rows (via the Append toolbox item), they don't get auto-geocoded, either. Am I missing a step?


title: Python Error in recognising ENVI file

text: I am trying to read a .bip AVIRIS file using GDAL in Python. here is the code:               from osgeo import gdal     from osgeo import ogr     from osgeo import osr     from osgeo import gdal_array     from osgeo.gdalconst import *          data = gdal.Open(filename)      Here is the Error: "filename.bip does not exist in the file system, and is not recognised as a supported dataset name." This is while the same code can read another ENVI (.bil) file! I would appreciate any hint on it,


title: Is it possible to override autoincrementing keys when QGIS inserts into SQLite?

text: I have a table, without any geometry, loaded in QGIS 1.8.0 (from OSGeo4W) via the `Add Vector Layer` option. The table resides in a SQLite/SpatiaLite database, and has a field named "id" that was created with `INTEGER PRIMARY KEY` option. From SpatiaLite GUI or SQLite Manager in Firefox, I can add a new row using any unique 64-bit integer value for "id". However, an attempt to add the same row from QGIS results in "id" being autoincremented. My question is: Is there a way to override or fix this behavior? According to http://www.sqlite.org/autoinc.html it is valid to assign a value (within constraints) to an `INTEGER PRIMARY KEY` field.


title: How to batch merge rasters using FOSSGIS tools?

text: I'm working with QGIS 1.8 (gdal 1.9.2) and GRASS 6.4.3...but on Windows 7. I "just" want to merge my raster. I've got a folder with 50 GTiff raster from Lidar data. The aim is to merge raster before or after an import into GRASS in order to do some contours... Before GRASS, I tried gdal_merge with QGIS's gui, without any option. I had an error message : "permission denied". So I tried it in GRASS with r.patch but : First the gui don't allows you to select multiple raster at one time... Then I don't know how to batch with GRASS functions in Windows...I tried the GRASS command line but can't locate in the current mapset (didn't understand how). When I had the windows installer for QGIS, this commandline putted you in the current mapset but it's not the case with the basic installation of QGIS. I'm quite locked...Any help ?


title: Floating Point UniqueValueRenderer

text: I am trying to build a unique value renderer using ArcObjects .NET. All the samples that I have seen use the attribute table to build the UniqueValueRenderer. However, my raster does not have an attribute table and when I try to build the attribute table I get the following exception "Can't build VAT for multiband or floating-point raster dataset". Currently my renderer is defaulting to stretched. When I look at it in ArcMap - Properties - Symbology it is set to stretched (this is what we default rasters to in our code if there is not an attribute table). In ArcMap, I can click on Unique Values under Properties - Symbology. A message box pops up stating "Unique values do not exist. Do you want to compute unique values?" If I click yes, then my unique floating point values are there. How can I programmatically define these unique values and set this is as the default symbology / renderer for my floating point raster?


title: How to work around "OGR-failure while deleting field 13. Invalid field index"?

text: this post might be related to this one: QGIS Field Calculator Results Removed When Saving Edits * * * But there is   1. no suitable answer for me   2. It is not allowed to ask another time in the other post I guess so I created a new one. * * * I have got a Shapefile with all boroughs in France (projection = EPSG:2154 - RGF93 / Lambert-93) I like to add some columns: area, Population, average income. The answer to make the columns bigger for area calculation does not work because the maximum size seems to be 10 in my Qgis version (1.8 Lisboa). If I like to add simply numbers (population number for example that does not work either. It does not matter if I use "interger" "real" or "string"...while saving it tells me that saving did not work because it was not able to delete attributes. OGR-failure while deleting field 13. Invalid field index Afterwards the result seems to be different: sometimes it can save the result, sometimes there will be set "NULL". It can not be so difficult to add a column in a shapefile can it? To show you what I am working with 2 screenshots (I added the last 2 columns): ![failure](http://i.stack.imgur.com/7Lb6h.jpg) ![results](http://i.stack.imgur.com/IRP22.jpg)


title: Moving legend items up and down in QGIS

text: How can I change sort order in layer properties styles tab? The only way I've found is to change their labels and order by it.


title: Python: Grabbing X,Y coordinate from matching fields

text: ![enter image description here](http://i.stack.imgur.com/uDhaL.png) Here is my table. What I am trying to do is write a python script to go over the Zip_Crid values. If there are matching values, I need to grab x1, y1, x2, y2 and put it in a new field. Haven't really worked with python before, where should I begin?


title: JTS - Java Topology Suite - Calculate radius from a given point

text: From point A, I want to check if a point B is within a certain radius from A. How calculate it usindo JTS? I tried this:               Coordinate posicaoAtual = new Coordinate(-12.974900375893823, -38.48034163207694);      Localidade localidade = new Localidade(-12.976353604547635, -38.47492356986686);     Geometry ponto = new GeometryFactory().createPoint(localidade.getCoordenada());     Geometry pos = new GeometryFactory().createPoint(posicaoAtual);     Geometry raio = ponto.buffer(100.0, 8); // HERE     System.out.println(raio.contains(pos));      My doubt is to know which the distance unit from the first parameter of the buffer method?


title: "General function failure." creating Table with ESRI File Geodatabase API

text: **What I'm trying to do:** I am creating a C# .NET application for my organization to convert a folder full of shapefiles into a file geodatabase system. **What do I have?:** My application reads the folder of shapefiles, and puts them into memory. All shapes have been successfully converted to a list of custom objects I made. **My code:**               if (Directory.GetFiles(newFileGeoDBpath).Count() == 0) {         Geodatabase.Create(newFileGeoDBpath);     }          Geodatabase geodatabase = Geodatabase.Open(newFileGeoDBpath);      Table table = geodatabase.CreateTable("Anything!!!"); //crashes here     geodatabase.Close();      It crashes when trying a create a new table... I can't add shapes to a table, if I can't create one. Is there any better examples, besides the ones provided with the API of adding shapes to a fileGeoDB? **EDIT:** _Details as requested by Vince_ Windows 7 64 bit VS 2010 - Application is being run as 64 bit. Version: v 1.3 of the API I thought the XML portion was optional? I will make an attempt and update my question.


title: Snap Polygon vertices to Point features

text: I have a set of polygons drawn using an inaccurate GPS PDA. these polygons are sharing nodes. after post-processing GPS data to get accurate coordinates i have a set of all the nodes, then i have to update the vertices of already drawn polygons. the relation between old and new vertices is based on the nearest point. how could i modify polygon vertices in this case, using(arcengine , arcobjects) or just a GIS algorithm. May be the whole process is not good enough so any suggestion is welcome.


title: Looking for .NET Heat Map API / SDK

text: I have a requirement to add a 'heat map' style layer to a Javascript/jQuery based application that I am working on. So far, I have used heatmap.js with no problems at all. My only concern with this approach is with how well it will handle a large number of points. I have asked for some type of upper bound, and the answer I got was 'a-lot', most likely going to be in the thousands. So I am thinking I may try to create a WMS style web service - where a bounding box, image height & width are the request parameters, the server queries the points and does the rendering, and the resultant image sent back to client. Does anyone know any such open source apis/sdks that contain this functionality? I know such tools exist in qgis and arcgis desktop - but i am looking to do this in a service.


title: ArcSDE error: "Maximum number of connections to instance exceeded "

text: I have an ArcSDE error: " * _Failed to connect to database. Maximum number of connections to instance exceeded "_ Did what ESRI recommended :" > Check the value set for the CONNECTIONS parameter in the SERVER_CONFIG or > sde_server_config ArcSDE system table to be sure it is sufficient. If it is > not, increase the number of allowed connections using the sdeconfig command > with the alter operation" and still nothing works on my database server. (I can't create features, import, etc). We use a Microsoft SQL Server. And it's an enterprise geodatabase Thanks a lot. * * * Apparently it was a schema lock. I had a map service started which used that geodatabase. So how can I edit the geodatabase and have a started map service in the same time. I've already created several versions and tried to publish just the versions, but they still lock my gdb.


title: Elevation correction ASTER GDEM2 (+10m)

text: My ASTER GDEM2 has a average elevation error of -10m for my area. (I testet it using 90 defferential GPS reference points and compared it with old Topo25, SRTM, and GoogleMaps). I would like to correct it by +10m for every pixel. Any idea? Whats the reason for that "error"? It is measured at a different height reference? Thanks! Martin


title: How can I extend L.mapbox.legendControl?

text: I'm trying to override legendControl.onAdd, but I can't seem to do it early enough. I can override the other functions, since they're fired after map init, but I can't seem to change the prototype? How should I do it? This doesn't work:               L.mapbox.legendControl.__proto__.onAdd = function(map) {         console.log('map.legendControl._onAdd');              this._container = L.DomUtil.create('div', 'map-legends-wrap-container');         this._subContainer = L.DomUtil.create('div', 'map-legends wax-legends');         L.DomEvent.disableClickPropagation(this._container);              this._container.appendChild(this._subContainer);              this._update();              return this._container;     };


title: What should I expect in this case?

text: If while I'm executing gdalinfo against a raster layer I get a positive value for pixel size Y, what of these scenarios is possible?   * The raster's map coordinate system sets the origin of the data in the upper left corner.   * The raster is not georeferenced and its coordinates are interpreted as image pixels, instead of map units.   * Any other situation (corrupt file?, other?) I know this value is almost always negative. For example, in case of UTM coordinate system. But, how common is to have a positive value for Y coordinate? Many thanks in advance


title: QGIS Fill Ring tool problem

text: Using the new (in 2.2) 'Fill Ring' tool and it doesn't seem to work correctly in all cases. It seems that if I try and digitise the new ring close to the polygon boundary, TWO (or more) geometries will be created (and will get prompted to fill the attributes for), one for each adjacent polygon. If I digitize the ring around the polygon center everything works fine. Same behavior for shapefiles and Postgres layers that I tried. The screenshot below shows an example. The attributes form displayed is the SECOND one and refers to the attributes of the polygon on the LEFT of the polygon I added the ring to. ![enter image description here](http://i.stack.imgur.com/C25Hm.png) Anyone got any ideas? Is this a bug?


title: Shortest distance from pgrouting with cost not being length: optimal solution

text: I'm using pgrouting shortest_path_astar with a complicated cost function. It gives back the edges, vertices and costs. I would like to get the total distance of this shortest path (i.e., the sum of lengths of edges, not the sum of costs of edges). For now, I use:               SELECT SUM(t2.length) FROM (         SELECT edge_id FROM shortest_path_astar(..., TRUE, TRUE)     ) as t1,     routing.edges as t2 where t1.edge_id = t2.edge_id      Is it the most efficient way to get the distance? I need to run this query many times and I want to be sure to do it in the fastest way.


title: Is it possible to access the Catalog window programmatically in ArcMap?

text: Specifically, I would like to get the selected object, similar to the way you do in ArcCatalog through `IGxDialog.SelectedObject`. Is there a way to do this? So far I can get the Catalog DockableWindow but am stumped from there.


title: Applying layer file not formating labels in ARCGIS 10

text: I'm applying layer files in ARCGIS 10 using the ApplySymbologyFromLayer_management command via a python script tool. The symbology part works fine but the lables are not formatting. It is my understanding that applying a layer file should also format the lables as well. Has anyone else had this issue? Any ideas would be appreciated. Here is the code I am using. * * *               import arcpy, os     from arcpy import env          ## Create parameters for modelbuilder user specified fields.     symbolize_layer = arcpy.GetParameterAsText (0)     symbol_lyrFile = arcpy.GetParameterAsText (1)               arcpy.AddMessage("**** Step 1 ****")     #  apply symbology     arcpy.ApplySymbologyFromLayer_management (symbolize_layer, symbol_lyrFile)          arcpy.AddMessage("**** Step 2 ****")     # refresh the view     arcpy.RefreshActiveView()     arcpy.RefreshTOC()          arcpy.AddMessage("**** OK now were Done ****")


title: How to batch convert Shapefiles to GMLs?

text: I know how to do this within Quantum GIS but I have 30+ shapefiles that need to be converted to GML. Is there a Python script sample or plugin for Quantum GIS that I can use to convert my shapefiles to GML?


title: Lambert Conformal Conic Projection Inverse Transform after Fitting Ellipse to Noisy LLA Data

text: I'm using Python and Basemaps to convert some Lat Lon data to the LCC projection and back again. Since it is noisy data representing the footprint of a satellite, I am filtering it and fitting an ellipse to it. I'm noticing that the inverse transform from my fitted ellipse data is causing some strangeness ( See the first subplot in the attached image ). I don't understand why the southern most latitudes show this "jumping" on the inverse transform. Can somebody please explain it to me? Also, Wikipedia mentions that LCC is artificially cutoff after -30 degrees latitude. Is there a better projection to use for this or any other suggestions in general? Thanks! Plots: ![Left Most Plot shows the Black Line which is inverse of fitted Ellipse to Noisy Lat Lon Data](http://i.stack.imgur.com/eKCI0.png) Here is how I get my coordinates: > Get XY m = > Basemap(projection='lcc',lat_1=64.25,lat_2=65.75,lat_0=65,lon_0=-19.,width=12000,height=14000) > xb, yb = m(lons,lats) http://en.wikipedia.org/wiki/Lambert_conformal_conic_projection


title: "US NED" OnEarth JPL NASA Tiled WMS working?

text: I am working on developing a computer software application that displays a 3D terrain based on the 30m elevation data for USA. I am using the OnEarth JPL NASA Tiled WMS for fetching the elevation data for the "us_ned" layer using a series of WMS requests (e.g., http://wms.jpl.nasa.gov/wms.cgi?Service=WMS&VERSION=1.1.1&request_short_int&bbox=-105,35.33333,-103.6667,36.66667) based on a small local area displayed on the application screen. For example, I might zoom into Jackson, WY and want to show the elevation as a 3D terrain. Everything was working fine with the OnEarth JPL NASA elevation tile web service a couple of weeks back. However, the web service seems to have stopped working. It is no longer returning the tiles, but returning the following exception:               <?xml version='1.0' encoding="UTF-8" standalone="no" ?>     <!DOCTYPE ServiceExceptionReport SYSTEM "http://wms.jpl.nasa.gov/exception_1_1_1.dtd ">     <ServiceExceptionReport version="1.1.1">     <ServiceException>     This server no longer provides full WMS services!     </ServiceException>     </ServiceExceptionReport>      This same web service was working earlier with Global Mapper, but now it is not working as well. I know that the OnEarth JPL NASA for MODIS image tiles was being phased out, but that should have not impacted elevation tiles (I hope)? Can anyone please confirm whether the web service is down forever, or there is another way I should be calling it? I know that I could setup my own WMS tile server, pay ESRI a bunch of money, or figure out the Microsoft BING Elevation API. But this was just recently working and I don't want to start over on this project.


title: Tool for adding buildings to OpenStreetMaps by tracing their countours on satellite maps

text: There is an OSM map of some city, in which several buildings are missing (see the OSM map and the Yandex map of the same place below - there are lots of spaces, which look empty in the OSM map, but have buildings in the Yandex map). ![OSM map](http://i.stack.imgur.com/UViyz.png) ![Yandex map](http://i.stack.imgur.com/o86jB.jpg) I want to add missing buildings to the OSM map using a tool, which   1. displays the OSM objects over a satellite map (such that I can see what objects are present on the satellite map, but not in the OSM map),   2. allows to add new objects (ways) to the OSM map by clicking points on the satellite map and connecting them to a polygon and   3. allows to save my changes (coordinates of newly added buildings) locally in a file. What tool are there, which support these functions?


title: How to integrate Geoserver with Tomcat Server 7 in Eclipse

text: I have imported .war of geoserver into eclipse and synchronized with my application. As I created a new workspace and published everything fine. When the server is restarted I couldn't able to find my created workspace. And even, I couldn't able to find in the workspace folder. Please help me...


title: How to create one print layout with 3 different maps of the same area?

text: I have a problem with the Qgis print composer: I want to create a print layout with 3 different statuses of the same geographical area, stylizing with 3 different fields. So far everything works but when I do the export to PDF or image format all 3 different maps are correct at the last style used. I can not figure out how to solve the problem... (same problem here: How to avoid map refreshes in print composer? but this solution seems to me it does not work) I tested on the following systems:   * Windows XP + Qgis 1.8 (postgres/postgis data)   * Windows 7 + Qgis 1.8 (shape data)   * Ubuntu 13.04 + 1.9 Qgis nightly b. (postgres/postgis data) thank you very much


title: How to match google openlayer and shapefile data?

text: i've installed openlayers plugin and opened google satellite layer in qgis. for a point in beijing china. it's show 12961957,4852605 coordinate in blew status bar. and it should be 116.33710,39.83818 in longitude & latitude. and i have some shapefile which in lat & long coordinate for beijing area. so how can i overlaped the shapefile data on top of google layer in correct position ? anyone can help me. some hint. thanks.


title: Histogram Matching Python

text: I am trying to do histogram matching using Python to improve the mosaicking process of multiple overlapping rasters. I am basing my code on that found at: http://www.idlcoyote.com/ip_tips/histomatch.html To date, I have managed to clip the overlapping area of two adjacent rasters and flatten the array. so I have two 1 dimensional arrays of the same length. I have then written the following code based on that found at the above website. In the code shown I have substituted two very small datasets for the gd and bd images.               bins = range(0,100, 10)          gd_hist = [1,2,3,4,5,4,3,2,1]          bd_hist = [2,4,6,8,10,8,6,4,2]          nPixels = len(gd_hist)          # here we are creating the cumulative distribution frequency for the bad image     cdf_bd = []     for k in range(0, len(bins)-1):         b = sum(bd_hist[:k])          cdf_bd.append(float(b)/nPixels)          # here we are creating the cumulative distribution frequency for the good image     cdf_gd = []     for l in range(0, len(bins)-1):         g = sum(gd_hist[:l])         cdf_gd.append(float(g)/nPixels)                # we plot a histogram of the number of      plt.plot(bins[1:], gd_hist, 'g')     plt.plot(bins[1:], bd_hist, 'r--')     plt.show()                  # we plot the cumulative distribution frequencies of both images     plt.plot(bins[1:], cdf_gd, 'g')     plt.plot(bins[1:], cdf_bd, 'r--')     plt.show()          z = []     # loop through the bins     for m in range(0, len(bins)-1):              p = [cdf_bd.index(b) for b in cdf_bd if b < cdf_gd[m]]          if len(p) == 0:             z.append(0)         else:             # if p is not empty, find the last value in the list p             lastval = p[len(p)-1]                  # find the bin value at index 'lastval'             z.append(bins[lastval])          plt.plot(bins[1:], z, 'g')     plt.show()          # look into the 'bounds_error'     fi = interp1d(bins[1:], z, bounds_error=False, kind='cubic')       plt.plot(bins[1:], gd_hist, 'g')     plt.show     plt.plot(bins[1:], fi(bd_hist), 'r--')     plt.show()      My program plots the histograms and cumulative frequency distributions successfully...and I thought that I had the part of getting the transformation function 'z' correct....but then when I use the distribution function 'fi' on the 'bd_hist' to try to match it to the gd dataset it all goes pear-shaped. I am not a mathematician and it is highly likely I have overlooked something fairly obvious but if anyone can offer any help that would be absolutely amazing. Thank you for your time. Becky


title: Code to compute weights for Inverse Distance in C#

text: I want to test my data using IDW in C#. As I am new to C#, can someone please help me by providing the code to compute the weight for the point?


title: How to know Who is reconciling the Version to DEFAULT(Parent) in Geodatabase

text: We have Public DEFAULT version and working Versions(Child) have been created directly from DEFAULT. We have also provided Reconcile-Post Option to group of Users placed at separate locations. Hence sometimes when I go ahead with reconciling my version It says some one else is reconciling the version, Can I identify which USER is reconciling on Database at that moment.


title: Why does the map legend not update in print composer?

text: I'm having trouble with the QGIS map composer: If I change the map in QGIS and reopen the map composer and click on update, the legend elements DON'T update! For example: Before one certain country was filled with blue, now with violet. The little box in the legend showing the color of that map element doesn't change. It stays blue. How can I solve this problem?


title: Leaflet markers do not open popup on click (with markerclusterer)

text: I just started using Leaflet and Marker Clusterer to organize the markers. **Problem #1:** When an unclustered marker is clicked, no popup appears. **Problem #2:** When a cluster is clicked several times, all the markers within that cluster appears, and when one of this marker is clicked, its popup appears! However, after closing the popup by clicking on the map, clicking on any of these clustered markers again do not open any popups! If I only have 3 unclustered markers, the popup works fine. However, as more markers are added, once a cluster forms, clicking on the marker within any cluster will not cause the popup to open! **Initializing markerclusterer**               markers = new L.MarkerClusterGroup();     map.addLayer(markers);      **All markers added to markercluster`markers`** A loop calls on the `render` function to create the `marker` and add it to the markerclusterer's array `markers`. _(ignore the backbone.js code)_               ListingMarkerView = Backbone.View.extend({         template: _.template( $('#tpl_ListingMarkerView').html() ),              render: function() {                  // Create marker             var content = this.template( this.model.toJSON() );             var marker = new L.marker(                 [this.model.get('lat'), this.model.get('lng')],                 {content: content});             marker.bindPopup(content);                  // Add to markerclusterer             markers.addLayer(marker);         }     });      **Without markerclusterer** If I add the `marker` directly to `map` instead of the markerclusterer array `markers`, the popups work fine, so I guess the problem has something to do with markerclusterer. Did I do something wrong that resulted in such behavior of the popups? All help appreciated, thanks!


title: Deleteing the acutal shapefiles from the geoserver

text: I would like to delete the actual shapefiles that I have previously uploaded and published on a geoserver. Is this possible? I know how to delete the workspaces and data stores using REST commands but I don't know how to delete the actual shapefiles. Any ideas?


title: Possible to rollback from QGIS 2.0.1?

text: I'm running Ubuntu 12.04 LTS. Unfortunately Ubuntu has automatically upgraded my QGIS from 1.8 to 2.0.1. And 2.0.1 which creates GeoJSON files that OpenLayers cannot parse. Would be best for me to rollback to 1.8, so I can get on with my work for now. Can anyone help? thanks Mini


title: Creating New Modules in Saga GIS

text: I have recently started development with SAGA and have run into some difficulties. I have checked out the svn code base and have successfully compiled in VS-10 (except the project ViGrA). Now I just want to create a new project and begin my development. I already have a project with some modules. I have checked the entries for them in the required MLB_interface.cpp. The project was supposed to be ready for a batch build, so I just added it as an existing project to the full saga solution. When I build this new project I get errors like: > error LNK2001: unresolved external symbol "__declspec(dllimport) public: > static class CSG_String __cdecl CSG_String::Format(unsigned short const > *,...)" (_ _imp_?Format@CSG_String@@SA?AV1@PEBGZZ) for all of my module cpp/obj files. Please advice! Thanks


title: How can I create drawFeature control in leaflet?

text: I want to call a callback on leaflet after drawing a feature. In openlayers I do it like this :               customControl= new OpenLayers.Control.DrawFeature(                         tempLayer,                         OpenLayers.Handler.Point,                         {                             title: 'Custom Control',                             displayClass: "olControlClass",                                                     featureAdded: onCustomControlCallback                         });      How can I do the same in leaflet ?


title: "Hello World" introduction for ubuntuGIS?

text: Is the GIS equivalent of "Hello world!" = "You are here?" I installed ubuntuGIS recently, but I'm a little overwhelmed navigating the network of many packages to begin accomplishing simple things. I'd like to explore GIS more thoroughly (I have made maps using R and techniques like these). Are there any resources minimally illustrating going from... Point A: `sudo add-apt-repository ppa:ubuntugis/ppa && sudo apt-get install qgis` to Point B: Here is a simple world map with a marker for your location? Or could someone propose this "hello world" to get newbies started? I think a tiny, simple project like this---"print where you are in the world" analogous to "print words on a screen" would be a very good learning tool for providing some narrative, tutorial context---something like giving a newbie a "roadmap" to making maps. (Not required, but my personal preference on a thing like this would be that the output be d3 or kartograph, something in-browser.)


title: Index for data driven page created atlas

text: We have a parcel atlas created from a data driven page document where each page represents a particular 16 mile square township. I have created my own index of polygons instead of ArcGIS index creation because we are looking for a slight overlap. I would like to create or tweak a python script to clip the area for each page and pull out owner name information for each township; better yet, would be each 1 square mile section in order to create a index of names/ area. Anyone know of something similar or can point me in a direction of more information. I am relatively new with python but understand the basics pretty well. Thanks! Megan


title: Rectangle for Clip_management using arcpy

text: I want to draw an interactive rectangle on screen which should be used as the extent for clipping geometry for raster datasets. Any suggestion?


title: Where to find QGIS API documentation?

text: Where to find QGIS API documentation? Can someone please share a QGIS API documentation link, i couldn't find it anywhere on the site. please help thanks


title: Sextante (SAGA) "Normalisation" problem

text: I don't understand the SAGA-module "NORMALISATION". In SAGA Gis I just can define one MIN and one MAX value. In QGis I have to define 2 of both. How to handle that exactly?


title: PostGIS - Why are my ST_Union results (stored in a new table) outputting as points and not polygons?

text: I've decided to move to PostGIS to merge polygons rather than spend hours doing in in QGIS. However, I'm hitting an issue where the results of an ST_Union query are being stored, or at least outputted as points and not merged polygons. Here is my SQL query:               DROP TABLE opengeo.london_merged;     CREATE TABLE opengeo.london_merged (       the_geom geometry(MultiPolygon),       CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),       CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTIPOLYGON'::text OR the_geom IS NULL)     );          INSERT INTO opengeo.london_merged      SELECT ST_Union(the_geom) AS the_geom     FROM opengeo.tq_building      This is the output from the new opengeo.london_merged table that I'm getting: ![enter image description here](http://i.stack.imgur.com/TWkhf.png) However, if I just run a straight this same SQL on the orignal opengeo.tq_building table as an SQL view in GeoServer (so not creating a new table, just outputting the results):               SELECT ST_Union(the_geom) AS the_geom     FROM opengeo.tq_building      I get the following, correct output: ![enter image description here](http://i.stack.imgur.com/vpiza.png) There's obviously something going on with how the new table is created or how the data is stored / retrieved from it, I just don't yet know enough about PostGIS to put my finger on it. How can I store the ST_Union output in a new table so that the output is the expected polygons? **Edit** I ran the following SQL view on the opengeo.london_merged table in GeoServer and it gives me the correct output, suggesting that the data is stored correctly but it being outputted / given to me wrong when I view the table normally. Perhaps a type or SRID issue? I don't know.               SELECT the_geom FROM opengeo.london_merged      I'm still keen to solve this.


title: Problem while Reprojecting

text: When I try to transform geometries from SRID 2321 to 4326, the transformed coordinates are swapping. Why is this happening? Below are the details **Source** SRID - EPSG:2321 Source Data Projection Definition:   * WGS84 Bounds: 31.5000, 35.9000, 34.5000, 42.1000    * Projected Bounds:364573.0642, 3975550.1410, 635426.9358, 4663916.1654    * Area: Turkey - 31.5°E to 34.5°E below are the sample points :   * 470343.6,4402310.3   * 470105.0,4402410.8   * 470236.5,4402700.1 **Destination** SRID - EPSG:4326 Below are the sample points after conversion:   * 32.6539666652,39.753664831   * 32.6511781435,39.7545616174   * 32.6526994043,39.7571716999 When I try to project these point on a Google map they are referencing to Iraq when they should refer to Turkey. But when I tried with reversing the coordinates of the above points then they are referencing to correct locations on Google maps Why is this happening? Do I miss anything? **EDIT:** Software used for conversion : Spatialite v 4.1.1 and this conversion is confirmed using this link. Actually my source data is from ESRI personal geodatabase .mdb file. I used QGIS to import them into Spatialite and when I try to project them on google/Apple maps all data is referencing to a different location.


title: how to get all features from a layer in openlayers?

text: I deal with geoserver and openlayers,   i use WMS layer for displaying due to loading time is shorter than vector layer. My problem i couldn't get all features of a layer without event,   For sure i tray "getfeatureinfo" but it's return selected point not all And i try to use vector layer for getting the features, But it not work unless i add it to the map ,   which will take long long time in loading .! could any one help?


title: Open a local accdb file through ArcGIS API for Flex button

text: I have a web mapping application created with the ArcGIS Viewer for Flex application builder. Is there a way to create a button that will open a local accdb file? I dont want to connect to MS Access through my app. I just want to open the file. Thank you in advance Lia


title: Can't load Mxd file in MapControl?

text: I used below codes to save IMap into byte[] (to store in my database), but I can't load it again, are there anyone know what's wrong ?               public byte[] SaveMap(AxPageLayoutControl pageLayOut)     {         IGraphicsContainer graphicsContainer = pageLayOut.GraphicsContainer;         graphicsContainer.Reset();         IElement pElement = graphicsContainer.Next();         while (pElement != null)         {                if (pElement is IMapFrame)                {                    IMapDocument mapDocument = new MapDocumentClass();                    if (File.Exists(@"D:\test.mxd"))                              File.Delete (@"D:\test.mxd");                         mapDocument.New(@"D:\test.mxd");                    mapDocument.ReplaceContents((IMxdContents)pageLayoutTo2_3.ActiveView.FocusMap);                    mapDocument.Save(true, true);                         System.IO.FileStream objFileStream = new System.IO.FileStream(@"D:\test.mxd", FileMode.Open);                    byte[] result = new byte[objFileStream.Length];                    objFileStream.Read(result, 0, (int)objFileStream.Length);                         return result;                 }                 pElement = graphicsContainer.Next();           }     }          public void LoadMap(AxPageLayoutControl pageLayOut, byte[] result)     {         IGraphicsContainer graphicsContainer = pageLayoutTo2_3.GraphicsContainer;         graphicsContainer.Reset();         IElement pElement = graphicsContainer.Next();         while (pElement != null)         {                if (pElement is IMapFrame)                {                         System.IO.FileStream oFileStream = new System.IO.FileStream(@"D:\test.mxd", FileMode.Append);                    oFileStream.Write(result, 0, result.Length);                         IMapControl2 map = new MapControlClass();                    map.LoadMxFile(@"D:\test.mxd", Type.Missing, Type.Missing );                         IMapFrame frame = new MapFrameClass();                    frame.Map = map.Map;                    pElement = (IElement)frame;                         pageLayOut.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);                         return;                 }                 pElement = graphicsContainer.Next();           }     }


title: Scipy: Convert RGB TIFF to grayscale TIFF and output it on Matplotlib

text: I want to manipulate **RGB** bands in a TIFF file and output the **grayscale** map on matplotlib. So far I have this code, but I couldn't get it on grayscale:               import scipy as N     import gdal     import sys     import matplotlib.pyplot as pyplot          tif = gdal.Open('filename.tif')          band1 = tif.GetRasterBand(1)     band2 = tif.GetRasterBand(2)     band3 = tif.GetRasterBand(3)          red = band1.ReadAsArray()     green = band2.ReadAsArray()     blue = band3.ReadAsArray()          gray = (0.299*red + 0.587*green + 0.114*blue)          pyplot.figure()     pyplot.imshow(gray)     pylab.show()      And these are the arrays:               [[255 255 255 ..., 237 237 251]      [255 255 255 ..., 237 237 251]      [255 255 255 ..., 237 237 251]      ...,       [237 237 237 ..., 237 237 251]      [237 237 237 ..., 237 237 251]      [242 242 242 ..., 242 242 252]]          [[255 255 255 ..., 239 239 251]      [255 255 255 ..., 239 239 251]      [255 255 255 ..., 239 239 251]      ...,       [239 239 239 ..., 239 239 251]      [239 239 239 ..., 239 239 251]      [243 243 243 ..., 243 243 252]]          [[255 255 255 ..., 234 234 250]      [255 255 255 ..., 234 234 250]      [255 255 255 ..., 234 234 250]      ...,       [234 234 234 ..., 234 234 250]      [234 234 234 ..., 234 234 250]      [239 239 239 ..., 239 239 251]]      Any idea how can I fix this?


title: Openlayers toggle Marker controls and update information in external form fields

text: I am setting a form (for multiple addresses) where users can click on map by creating a marker for each address they want to input. Users can only create one marker per address and once the marker has been created I need to get nearest address and update the address field (form input) with that address. I am able to activate the marker drawing control but how can I make sure that user only creates one marker (not more) per address and then take that calculated address and update original text field (which fired the event) with new information. This is what I have so far:               <form name="xxxx">          POI 1: <input type="text">      click here to select address on map     POI 2: <input type="text">      click here to select address on map               </form>      Add more poi's       Here's the javascript code               poi_markers = new OpenLayers.Layer.Vector('markers',{styleMap:markerStyleMap});     poi_marker_control = new openLayers.Control.DrawFeature(poi_markers,OpenLayers.Handler.Point);                 function e_drawPOIMarker(evt){                  route_markers.events.register('featureadded',poi_markers,function(evt){              {              //get JSON response                          var place_name = getGeoName(evt.feature.geometry)                            poi_marker_control.deactivate();                                 //Somehow find the id of the text field which fired original event and update its contents with data returned in place name.             }         });             poi_marker_control.activate();           }      Basically I have to capture co-ods of locations clicked by users. Each field (in form) should only allow 1 marker. Suppose if I have unknown numbers of poi address fields (which users can add dynamically) whats the best way of capturing coods and then updating original fields with that information. At the moment I can activate and deactivate controls and allow users to create marker(s). but how can I know which marker belongs to which text field ? Thanks.


title: Problem setting IActiveView extent

text: I'm having a problem when setting the IactivewViews extent. I have done it the same way on other places and it does work. In here, it just throws me a LoaderLockException.                       // OLD CODE             private void airspaceElementsTreeView_NodeMouseDoubleClick(object sender, TreeNodeMouseClickEventArgs e)             {                 if (e.Node.Name.StartsWith("feature_"))                 {                     IFeatureWorkspace workspace = ((IDataset)_LoadedAirspace.UnderlyingObject.Table).Workspace as IFeatureWorkspace;                     ITable relatedTable = workspace.OpenTable(e.Node.Parent.Name);                          int oid = GetNodeObjectId(e.Node);                     IRow tempRow = relatedTable.GetRow(oid);                          IGeometry geom = null;                          if (tempRow is IFeature)                     {                         IFeature feat = tempRow as IFeature;                              geom = feat.ShapeCopy;                         ISpatialReference srCopy = ((IClone)geom.SpatialReference).Clone() as ISpatialReference;                              IEnvelope envelope = new EnvelopeClass();                         envelope.SpatialReference = srCopy;                         envelope.Union(geom.Envelope);                         envelope.Expand(.25, .25, false);                              // error here!                         IActiveView view = ((IMxDocument)m_application.Document).ActiveView;                         view.Extent = envelope;                     }                     else                         return;                          Marshal.ReleaseComObject(relatedTable);                 }                 else                     return;             }      By the way, this is inside a IDockableWindowDef. Any chance that THAT is ruining things for me? **edit** following the suggestions of my fellow friends, I did on two separate statements. First I'm trying to get IMxDocument from m_application and then casting it to IActiveView. It crashes when I'm trying to cast m_application.Document to IMxDocument (check new code below)                   private void airspaceElementsTreeView_NodeMouseDoubleClick(object sender, TreeNodeMouseClickEventArgs e)         {             if (e.Node.Name.StartsWith("feature_"))             {                 IFeatureWorkspace workspace = ((IDataset)_LoadedAirspace.UnderlyingObject.Table).Workspace as IFeatureWorkspace;                 ITable relatedTable = workspace.OpenTable(e.Node.Parent.Name);                      int oid = GetNodeObjectId(e.Node);                 IRow tempRow = relatedTable.GetRow(oid);                      IGeometry geom = null;                      if (tempRow is IFeature)                 {                     IFeature feat = tempRow as IFeature;                          geom = feat.ShapeCopy;                     ISpatialReference srCopy = ((IClone)geom.SpatialReference).Clone() as ISpatialReference;                          IEnvelope envelope = new EnvelopeClass();                     envelope.SpatialReference = srCopy;                     envelope.Union(geom.Envelope);                     envelope.Expand(.25, .25, false);                          // crash here                     IMxDocument document = m_application.Document as IMxDocument;                     IActiveView view = document.ActiveView;                     view.Extent = envelope;                 }                 else                     return;                      Marshal.ReleaseComObject(relatedTable);             }             else                 return;         }


title: How to handle list variable in a where clause

text: I am calling a list of names in from a text file and would like to use each of those to make new layers. I am running into an issue where I get "TypeError: cannot concatenate 'str' and 'list' objects" and am not sure how to get past this. I tried a couple of different where clauses without much luck. *EDIT: Working code below. The issue was I am calling the file directly in the walk and the list fields. so when I am appending the path I needed to remove the filename from os.path.join(dirpath, filename).                import arcpy, os     from arcpy import env          #setup input varible     ws = r"M:/"          cl = r"/Admin_1.shp"          dem = r"M:/SRTM/dem"          #set output as a variable     out = r"M:/Elevation_Masks/Masks"         with open('M:/Elevation_Masks/test_space/elevCountries_Test.txt', 'r') as f:         cb = [line.strip() for line in f]          CNTRY_NAME_FIELD = 'CNTRY_NAME'          for dirpath, dirnames, filenames in arcpy.da.Walk(cl, datatype="FeatureClass", type="ALL"):         for filename in filenames:             if CNTRY_NAME in [f.name for f in arcpy.ListFields(cl)]:                 for country_name in cb:                     print "Country:", country_name                     #setup where clause                     where = '"%s" = \'%s\'' % (CNTRY_NAME, country_name)                     arcpy.Select_analysis(os.path.join(dirpath), r"M:/Elevation_Masks/test_space/Country_Masks_Test/clipLayers_Test/clipLayers_Test.gdb/" + country_name + "_Boundary", where)


title: What is a best practice way to display/feed US census/tiger data into Google Maps API?

text: Let's say I want to display age, income, etc based on census block groups on Google maps. What is the best practice or recommended way to do so and also to allow to do simple queries?


title: Free software to publish tiles on Amazon S3?

text: Does anyone know of free software that would allow me to publish tiles on Amazon S3? I know of Arc2Earth, but it is not free.


title: Clipping a raster file into the multiple shapefiles

text: I splitted the single shapefile into the multiple shape file based upon their attributes. Now I want to clip the raster file into the each splitted shape file. How can I do it at once? I am using a arcgis version 10. I used spilt layers by attributes tools to split the shapefiles.


title: Null value error with Python in field calculator

text: I'm trying to calculate the angle between two line segments in ArcMap 10.0. I have a point shapefile, so I'm using the xy coordinates of 3 successive points to calculate the angle. I'm trying to use Python in the Field Calculator, because I'm not good with Python and I'm not sure how to use the Python window. I finally got my code to work, but encountered an error when two successive points have the exact same coordinates. This creates a null value, and I get the error "The field is not nullable". I tried to put in an if statement to change nulls to 0s, but I still get the same error. I'm not sure if I'm identifying null values correctly, though, because I've found conflicting ways to represent them, and have tried "None" and ' ' and "Null" to no avail. It seems like there should be some way to just replace the angle value with 0 if the calculation returns a null value. If anyone has any ideas on how to get past the null values, I'd really appreciate it. My code is below. Thank you!               Record = 0     def turnangle(X3, Y3):       import math       global X1, Y1, X2, Y2, Record       if Record == 0:  #first two rows aren't calculated         angle = 0         X1 = X3         Y1 = Y3       elif Record == 1:          angle = 0         X2 = X3         Y2 = Y3       else:               #calculating angle         BAx = X1 - X2         BAy = Y1 - Y2         BCx = X3 - X2         BCy = Y3 - Y2         dotprod = BAx * BCx + BAy * BCy         lenA = math.sqrt((X1-X2)**2+(Y1-Y2)**2)         lenB = math.sqrt((X2-X3)**2+(Y2-Y3)**2)         angle = math.acos((dotprod)/(lenA*lenB))*180/math.pi         X1 = X2       #advancing to next coordinates         Y1 = Y2         X2 = X3         Y2 = Y3       Record = Record + 1       if angle is not None:  #Here is where I'm trying to substitute 0 for null values         return angle       else:         return 0


title: PyScripter crashes at "import processing"

text: I set up PyScripter latest version, 2.7) to work with QGIS (2.2, OSGeo4W install) following http://chelahmy.blogspot.ca/2011/02/setting-up-pyscripter- for-quantum-gis.html. It generally works well. I'm running Windows 7. Here's my code:               import sys     import os     sys.path.append("C:\\OSGeo4W\\apps\\qgis\\python")     sys.path.append("C:\\OSGeo4W\\apps\\qgis\\python\\plugins")     from qgis.core import *     from qgis.utils import *     from qgis.gui import *     from PyQt4 import *     from PyQt4.QtCore import *     from PyQt4.QtGui import *          from processing.core.Processing import Processing      I need the sys.path.append to import qgis properly (rather than setting them in the PYTHONPATH). However, when I run the line               from processing.core.Processing import Processing      PyScripter simply shuts down, no error message, window just closes. Has anybody else experienced this? Same thing happens with other imports, like               import processing      I can import all the other plugin libraries (GdalTools, fTools, db_manager) with no problems - but I don't need those, I need the proceesing library to access all the nice OTB tools.


title: How to hyperlink an attribute in ArcGIS stand-alone tables?

text: I have arcgis 10.0 and have created a 1:m relationship in a File GDB. The 1:m relationship joins a table containing a URL to a point geometry. Some points have more than one file hence the 1:m join. After successfully creating the relationship, I wanted to use the field in the table to open the Image. It seems however, that ArcGIS stand-alone tables cannot use hyperlinks. Is there a workaround to this or can ArcGIS not use hyperlinks in stand-alone tables?


title: configuration de gestionnaire de license arcgis desktop 10 lors de l'installation

text: At the end of the installation of ArcGIS Desktop 10, we need to configure the license server administration. We find a problem when we write "localhost" as administrator of the license.                   Le nom du hote n'est pas valide


title: How to use Matplotlib in Python 2.6 installed by default ArcGIS 10.0 installation?

text: I am wondering if anyone has successfully run Matplotlib with the default installation of Python 2.6 that comes Arcgis 10.0. It does not seem to work with the python version installed by Arcgis 10.0. I was able to download and successfully run other packages such as Scipy. I have looked at many blogs where they talk about using other installations of python such as Python(x,y) etc. I prefer to run it without changing the version of python. Please post the details of the steps if anyone was successful in running Matplotlib. I have installed Matplotlib version 1.1.0. For example, the following command to check version works. > > > import matplotlib >>> >>> print matplotlib. **version** shows 1.1.0 I think the problem is with other things that comes with Matplotlib such as pylab, pyplot. Every time I try to import pylab or pyplot, the python IDLE environment is RESTARTED. Please give me the details of steps (from download, installation to running it) if anyone of your are successful in using with the default python installation of ArcGIS 10.0.


title: How to solve postgreSQL error when using QGIS Cloud plugin?

text: I created a qgis cloud account and installed the plugin. When I attempted to Upload Data using the plugin, I got the following error: postgreSQL exception: ERROR: "MD5" authentication with pgpool failed for user "database_name" Any ideas on how to fix this error so I can upload data? Thank you!


title: Is the Mercator projection or the Gall-Peters projection of the world more accurate?

text: I've been researching the concept of Eurocentrism, an ideology that views the world from a very European perspective. It's seen in literature, in academia, in world history, and perhaps even in cartography, as I've just found out. I found information on two types of world maps, the Mercator Projection map and the Gall-Peters Map. Supposedly, the Mercator Projection map is inaccurate, and an example of the political influences on map design. The Gall-Peters projection map as I've read, was a sort effort to counter that. The idea is that in the Mercator map, Europe is seen as a much larger land mass and appears reasonably large compared to the Gall-Peters projection, which shows Europe as a much smaller mass, almost insignificant, when compared to regions like Africa. So my curiosity brought me here. What is the accepted map projection of the world among Geographic Information Specialists and Cartographers and experts operating in related fields? Is the relation between the map projections and Eurocentric perspectives too much of a stretch? Is it neither? I'd love to hear from you all. Much thanks.


title: ArcObjects - Create a map in the background from layers and save as JPG

text: Given a series of map layers is it possible to arrange them in some order, add basic map features (north arrow, scale bar) and then save the result as a JPG. I'd like to accomplish this in the background as a button-triggered event in an ArcMap Add-On. I've seen ExportActiveView, but I think this will only work if I have the map layers already loaded into the ActiveView and on the users' screen. Is it possible to create this map behind the scenes?


title: Can GRASS read Spatialite geometry fields?

text: I've looked at the question: How to import SpatiaLite data into GRASS? already. The links there talk about two separate tasks: Passing SQL to a Spatialite DB from GRASS and connecting to Spatialite from GRASS. That is fine, but the point of Spatialite is to generate geometry metadata. So to get specific, let me explain what I am aiming to achieve. I create an SQLite database containing crime statistics with State Plane California Zone II (ft) projected coordinates (X_Coord and Y_Coord). All of that is done in R. I then run an SQL script that initializes the Spatialite script and adds geometry columns to the crime table, populates them, and then imports a couple of shapefiles. Now, my table still has the original X_Coord and Y_Coord fields, but it also has a SHAPE field that contains the geometry. I can import the data into GRASS because I have the coordinates, but it cannot recognize, as far as I can tell, the Point geometry contained in SHAPE. If I'm going through all of this to store geometry in Spatialite, it seems redundant (and time consuming) to then recreate geometry in GRASS. My aim is to use this database and have GRASS quickly extract the geometry and create a v.kernel raster for a density surface. It would also be nice to define the GRASS bbox based off the geometry stored in Spatialite (because it's easy to get that information from SHAPE instead of X_Coord and Y_Coord with Spatialite functions). But if GRASS currently cannot interface Spatialite geometry, there's an obvious disconnect. Is this just a feature GRASS is currently lacking or am I missing where on the web this is documented? I've spent the past couple of weeks looking into this and find the only solution is to redundantly create the geometry in GRASS based off X_Coord and Y_Coord, making the SHAPE field superfluous to this part of my analysis. (It's still useful for directly pulling into QGIS.)


title: Postgis - Converting 'geometry_dump' data type to a 'geometry' type

text: I am working with a dataset of property boundaries that has had ST_Union() run on it, thus all the polygons appear as a multi polygon in one row of a table. I would like to separate these out so I get a separate row/polygon for each shape. I have tried running ST_Dump() on this and it seems to separate out the polygons in the dataset but I can no longer view this in a GIS (I'm currently using QGIS). The data-type for this column says it's a 'geometry_dump' and I assume this is why I cannot visualize it. When I access my Postgres database from QGIS this new table is listed under the heading of datasets without a geometry. Does anyone know if there is something I can do to convert a 'geometry_dump' to a 'geometry'? Cheers, Roman


title: How to fix AttributeError: 'QgsVectorLayer' object has no attribute 'provider'?

text: When I try to start this algorithm I get an error: Traceback (most recent call last): File "C:/Users/trdina/.qgis//python/plugins\sextante\core\GeoAlgorithm.py", line 145, in execute self.processAlgorithm(progress) File "C:/Users/trdina/.qgis//python/plugins\sextante\algs\ftools\SumLines.py", line 85, in processAlgorithm spatialIndex = utils.createSpatialIndex(lineLayer) File "C:/Users/trdina/.qgis//python/plugins\sextante\algs\ftools\FToolsUtils.py", line 32, in createSpatialIndex provider = layer.provider() AttributeError: 'QgsVectorLayer' object has no attribute 'provider' I have no clue woh to solve this error. Any help is welcome. Best regards


title: Getting tiles from own Tilestache server to Openlayers application

text: The title sums up my current problem -   1. I have installed and set Tilestache in our webserver and the tiles are accessible in browser    2. My developer is working on an application to overlay some custom layers in kml format , over some basemap. We were using a basemap from Openlayers example, but now want to replace the url with tiles from the Tilestache. Which method we have to get the tiles ? XYZ ? UTF grids ? httprequest ? Experts please guide..


title: Shapefile postgis database to be queried by GeoDjango

text: I would like to create a Postgres database (postgreslq-9.1 and postgis-2.0.2) storing several shapefiles but I am uncertain on how to implement it so I'm searching for some advices. I have already done this last year but was pretty uncofortable beacause I had a table per shapefile, every table signed with its SRID (stored inside geometry_columns table - I was using postgis 1.5, so this was still a table and not a view as in postgis-2). I had to create a new table every time I converted a shapefile: pretty uncofortable because I had to query them using Django ORM, and create table mapping on the fly is not trivial. Is that possible to create one single table having an id indicating the original shapefile and the relative SRID **per geometry**? Otherwise have you got any better ideas? Thank you.


I am a true novice when it comes to ArcGIS and computer programs like this one. Am taking a class and cannot get the Project function to work under the tool box. The path I used was ArcTool box --> Data Management Tools --> Projections and Transformations--> Feature Path --> Project --> (This gives me a blank box and when I double click on the OK button I get the following message: Error 000735: Input Dataset or Feature Class: Value is required Error 000735: Output Dataset or Feature Class: Value is required Error 000735: Output Coordinate System: Value is required Can anyone tell me what I can do to remedy this. I would greatly appreciate it.


title: What are benefits of using ArcGIS with PostgreSQL over simply ArcSDE?

text: I understand that in order to benefit from ArcGIS spatial datatypes (ST_GEOMETRY, versioned geodatabases etc.) ArcSDE needs to be installed - and with the 10.1 client it is automatically installed (two-tier client>>DBMS, see **here**. However, I'm still trying to figure out the benefits of connecting an Esri front end to PostgreSQL (or Oracle). Is it purely back-end benefits, such as scalability /performance / redundancy / DBAability, or are there any extra front-end features offered by a PostgreSQL / PostGIS setup? The Esri **help** does not seem to elaborate FYI, I'm pretty much a database newbie who is designing a GIS dissertation aimed at illustrating why Postgresql / PostGIS / QGIS is more scalable / flexible (and cheaper) than ArcGIS / any database, so must also consider above scenario, which I can't afford to implement. Any links on this last point, especially regarding workflow (e.g. not having loads of shapefiles cluttering everything up) most welcome. d.


title: What is the formula to convert WGS84 coordinates to EPSG 2163?

text: I'm looking for the formula which converts wgs84 coordinates to epsg 2163. I know that QGIS can do the calculations but I need the exact formula. Cheers!


title: Spatial graphing, distribution of polygon area over latitudes

text: Given a polygon in a PostGIS database, what is the best way to create a graph of its distribution over a range of latitudes? I'm picturing a graph of a curve, where the proportion of the total area under the curve that lies between two points on the x-axis (latitude) represents the proportion of the polygon's area between those latitudes. The tools I have at my disposal are PostGIS and R. It seems like this should be a fairly routine exercise for which there are well established approaches that I simply haven't been able to find.


title: How do I add elevation data to a GeoJSON line from a GPX file?

text: I understand that elevation can be added as an optional third value in Coordinates (lat, long, elev) in a geojson file, but I can't find a way to actually do that from a gpx file. I've used ogr2ogr to convert gpx to geojson. I've also loaded the gpx file into QGIS and saved as a geojson file, but neither method brings in elevation when I do this to the gpx tracks. It will bring in elevation as a Property if I convert the points, but I need a line. Ultimately, I want to diplay these lines in Leaflet and build elevation profiles for them. Is there any easy way to do this other than parsing out the gpx files and building my own geojson files?


title: Creating a featureclass in memory and adding to the map error

text: I have the code below, which is supposed to create a new featureclass in memory then add the FC to the map, but I get a Automation error when I use the createfeatureclass method > Set pFC = pFeatWS.CreateFeatureClass("NewBus", pFields, Nothing, Nothing, > esriFTSimple, "Shape", " ") Here' s the full code               Public Sub InMemoryWorkspace()          'Create Workspace factory     Dim pWorkspaceFactory As IWorkspaceFactory2     Set pWorkspaceFactory = New InMemoryWorkspaceFactory          'Create workspace     Dim pWorkspaceName As IWorkspaceName2     Set pWorkspaceName = pWorkspaceFactory.Create("", "MyInMemoryworkspace", Nothing, 0)          'Get a handle on Name     Dim PName As IName     Set PName = pWorkspaceName          'Get handle on Featureworkspace     Dim pFeatWS As IFeatureWorkspace     Set pFeatWS = PName.Open          'Create field object and set number of fields     Dim pFields As IFields     Set pFields = New Fields     Dim pFieldsEdit As IFieldsEdit     Set pFieldsEdit = pFields     pFieldsEdit.FieldCount = 3          'Create objectID field     Dim pField As IField     Set pField = New Field     Dim pFieldEdit As IFieldEdit     Set pFieldEdit = pField     With pFieldEdit         .Name = "OBJECTID"         .Type = esriFieldTypeOID     End With     Set pFieldsEdit.Field(0) = pField          'Create postcode field          Set pField = New Field          Set pFieldEdit = pField     With pFieldEdit         .Name = "POSTCODE"         .Type = esriFieldTypeString     End With     Set pFieldsEdit.Field(1) = pField          'Create Shape field     Set pField = New Field     Set pFieldEdit = pField     With pFieldEdit         .Name = "SHAPE"         .Type = esriFieldTypeGeometry     End With     Set pFieldsEdit.Field(2) = pField               Dim pGeomDef As IGeometryDef     Dim pGeomDefEdit As IGeometryDefEdit     Set pGeomDef = New GeometryDef     Set pGeomDefEdit = pGeomDef     With pGeomDefEdit       .GeometryType = esriGeometryPolygon       Set .SpatialReference = New UnknownCoordinateSystem     End With     Set pFieldEdit.GeometryDef = pGeomDef     pFieldsEdit.AddField pField          ' Create featureclass      Dim pFC As IFeatureClass      Set pFC = pFeatWS.CreateFeatureClass("NewBus", pFields, Nothing, Nothing, esriFTSimple, "Shape", " ")           Dim pmemLayer As IFeatureLayer      Set pmemLayer.featureClass = pFC           Dim pmxdoc As IMxDocument      Set pmxdoc = ThisDocument      Dim pmap As IMap      Set pmap = pmxdoc.FocusMap           pmap.AddLayer pFC      Can someone explain where I am going wrong. Cheers


title: ArcGIS Routing on Android with Network Analyst module

text: I am currently developing an Android App for navigation purposes. It is based on the ArcGIS Android API. I have an ArcGIS server (with network analyst module) with all the needed data running. At the moment the App displays a map with the current position (point A) and has a search form, where you can search for a building (point B). The next step would be to display a route from "point A" to "point B". Does anybody know if there is a way to do that natively with the Android API? Also it would be possible to compute the routing information on a Java web service and then send it to the app, which only has to display it. I cannot find any sample code on the ArcGIS resource page.


title: IDRISI .vec file translator to ESRI shapefile format?

text: I've recently got hold of almost 20 years old IDRISI .vec files which I'd like to convert to ESRI shapefile. I've been struggling to find any translator for this file format so far. Does anybody know of a way to convert this type of file? Any help much appreciated! Many thanks. P.S. I'm going to post cross post on the http://forums.arcgis.com/


title: Sea routes network dataset

text: I am now working on a system to perform routing on a network of shipping lanes. See this dataset that I have as the base: http://geocommons.com/maps/109850 I have postgres and pgrouting up and running and I was able to run basic queries on the above dataset after importing it to postgres/postgis. However I have encountered some problems for which I seek advice from you.   1. A lot of lines have a duplicate, so say line A goes from source id 1 to source id 2 and line B goes from source id 2 to source id 1 and they overlap completely. This can be experienced a lot of times, and frankly, I would be happy to get rid of them so that finally I could have a network with only one polyline per direction.   2. I discovered that the distance value attributes of the referenced dataset are false. I noticed this when I checked the dataset with Global Mapper by using the Info tool. In the window that pops up, you can see the geometry textbox with bounds and length data. For a line I inspected that displayed a length value of 5300 some kilometres while the max value in the length column of the dataset was 110. I seek for a way to acquire those polyline length values, and assign them to a new attribute column. I have already found posts about how to do this in qgis and arcgis, but neither of them provided the results I wanted to get. The route dataset is in WGS84, latlon I hope the lines above make sense, and I appreciate any help! Thanks in advance! Tamas


title: Google maps POI download

text: Is there any legal way to search and download google maps POI , based on user location and search radius (50 km) + search criteria ( eg restaurent) ? Output shall be Kml/csv .


title: Why is the database connection lost on every query?

text: Every time I started a query the connection with the database is lost. Anyone familiar with this problem? This problem occurs with various query, for instance:               CREATE TABLE table     SELECT a.*, b.column, b.column      FROM table_a a JOIN table b      ON ST_intersects (a.geom, b.the_geom)     WHERE st_IsValid (a.geom);      I just updated from postgis 1.5.3 to 1.5.5 but that doesn't seem to help. Anyone an idea? * * * **UPDATE** I solved the problem. There was something wrong with the topology of one of the layers. After editing the query worked without disconnecting. I find it rather strange that the program and service crashes when the topology is not correct. Can someone shed a light on this issue.


title: uDig Connects to PostGIS but not QGIS

text: I recently installed QGIS 2.0.1 and uDig 1.40b on the same computer. PostGIS (from OpenGEO Suite) is running on my server. I'm able to connect to PostGIS and access layers with uDig but not with QGIS. Any thoughts why uDig works but not QGIS? I've tried previous versions of QGIS but without success. UPDATE: Thanks to all of you responding to my question. I have a few items to update my original question. Below is the parameters entered in QGIS to connect to PostGIS on a remote server (and, yes, the port for PostGIS is 54321). ![enter image description here](http://i.stack.imgur.com/BBi2d.png) The error message is below. ![enter image description here](http://i.stack.imgur.com/gGN78.png) I have recently connected to the PostGIS database with OpenJump without problems (as well as with uDig). Your thoughts are greatly appreciated. Scott UPDATE 2: In response to recent questions: (1) uDig, OpenJump and QGIS are all on the same client computer and OpenGEO is on a separate server, (2) pg_hba.conf has the following settings: # IPv4 local connections: host all all 0.0.0.1/0 trust, (3) postgresql.conf has the following settings: listen_addresses = '*' port = 54321, (4) tables and field names are all in lower case, and (5) same credentials used for QGIS, uDig and OpenJump. UPDATE 3: In place of the computer name (on the local network) I entered the IP address of the server and the connection was made between QGIS and PostGIS. It appears (in my situation) I cannot enter the hostname of the server to make the connection, only the IP. Anyone else encounter this phenomena? FINAL UPDATE: YES! The issue with the setting for IPv6 solved the problem. I'm able to connect to the server by its hostname as well as its IP address.