title: How to move an existing line to a parallel offset of another existing polyline

text: I'm currently trying to modify an existing line to be offset x distance from another existing polyline, is there any tool/addon that does that? I currently use Copy Parallel and create a dummy line that I move/snap my line to. I was just wondering if there was any way to directly specify that I want an existing line to be offset parallel from another line x units.


title: OGR CreateLayer returns TypeError

text: I'm trying to create a shapefile layer based on the spatial projection of an already existing raster. (Ultimately, this will be used in gdal.Polygonize to get a shape from my raster data.) However, I am encountering an error when trying to use CreateLayer.   * raster_uri- location of the raster file that I want to make into a shapefile. I know this exists.   * out_uri- Location I want the new shapefile to print to.   * layer_name- String layer name input.   * filed_name- String field name input. This is the code that I was using.               raster = gdal.Open(raster_uri)     driver = ogr.GetDriverByName("ESRI Shapefile")     ds = driver.CreateDataSource(out_uri)          spat_ref = osr.SpatialReference()     proj = raster.GetProjectionRef()     spat_ref.InportFromWkt(proj)          layer = ds.CreateLayer(layer_name, spat_ref, ogr.wkbPolygon)      It yields this error message:               TypeError: in method 'DataSource_CreateLayer', argument 2 of type 'char const *'      Outputing the type of spat_ref shows up as `<class 'osgeo.osr.SpatialReference'>` so I'm unsure what char const * it's complaining about. Interestingly, if I try to do this all from a python shell using the same file, layer creates just fine. Any ideas on what might be the issue?


title: A WMS layer fails to show at some scales

text: I am tring to overlay my WMS layer on Google but at a particular scale it doesn't show the data of my layer. What can be the problem?               <html xmlns="http://www.w3.org/1999/xhtml">     <head id="Head1" runat="server">         <title>Untitled Page</title>         <script type="text/javascript" src="http://maps.google.com/maps/api/js?v=2&sensor=false"></script>         <script type="text/javascript" src="http://openlayers.org/api/OpenLayers.js"></script>         <script type="text/javascript">             var map;             function init() {                 map = new OpenLayers.Map('map', { allOverlays: true });                 map.addControl(new OpenLayers.Control.LayerSwitcher());                 var gmap = new OpenLayers.Layer.Google("Google Streets");                 var wms = new OpenLayers.Layer.WMS("OpenLayers WMS", "http://localhost:5555/geoserver/cite/wms",                 {                     layers: 'cite:view_test',                     transparent: true,                     projection: "EPSG:4326"                 });                 map.addLayers([gmap, wms]);                 map.zoomToMaxExtent();              }     </script>     </head>     <body onload='init();'>         <div id='map' style='width: 500px; height: 500px; background-color:Gray'>     </div>     </body>     </html>      I am new to Openlayers: can anyone help me?


title: How to store and serve the time a GeoTIFF refers to

text: I'm writing software that creates temperature maps as GeoTIFF files. I need to store the time to which the GeoTIFF file refers; not the raster creation time, but the time during which the temperature was as the raster data says. What I did was store it in the GeoTIFF file as a GDAL metadata item; i.e. in a GDALMetadata private tag. The metadata item is called "TIMESTAMP". I now need to serve the GeoTIFF file as a WCS with MapServer.   1. Is storing the time as a GDAL metadata item a good idea? Is it standard practice?   2. How is this going to be accessed with WCS? Through `DescribeCoverage`?   3. How can MapServer be told to get the metadata item from the GeoTIFF and somehow serve it to the client?


title: Custom marker/icon based on attribute data in leaflet/geojson

text:
I would like to add different icons to my leaflet map based on attribute data. My layer is a geojson layer, it is coming from geoserver with jsnop. I can see my data with the default marker, but when I try to change the code to add custom style I get this error: > Cannot read property 'labelAnchor' of undefined It would be great if I can get some help.               var WFSLayer = null;     var ajax = $.ajax({         url : URL,         dataType : 'jsonp',         jsonpCallback : 'getJson',         success : function (response) {            WFSLayer = L.geoJson(response, {                     pointToLayer: function(feature, latlng) {                         var smallIcon = L.Icon.extend({                             options: {                                 iconSize: [27, 27],                                  iconAnchor: [13, 27],                                  popupAnchor:  [1, -24],                                 iconUrl: 'leaflet/icons/' + feature.properties.pcp + '.png'                             }                         });                         var myIcon = new smallIcon();                         return L.marker(latlng, {icon: smallIcon});                     },                 onEachFeature: function (feature, layer) {                     popupOptions = {maxWidth: 600};                     layer.bindLabel('
'+feature.properties.musno+'
');                     //sidebar.setContent('
'+feature.properties.musno+'

'+'
'+feature.properties.exchange_name+'

'+feature.properties.pcp, popupOptions);                     layer.bindPopup('
'+feature.properties.musno+'

'+'
'+feature.properties.exchange_name+'

'+feature.properties.pcp                         ,popupOptions);                 }             }).addTo(markers.addTo(map));         }     });


title: Calculating the distance between points and multiple lines

text: I have a point shapefile of about 100 points and a polyline shapefile of about 200 lines and I am interested in finding the distance between every point and each individual polyline. My polyline shapefile represents individual fault lines in a region. So for each polyline (or fault line) I want to calculate the distance to each of my points. I know that I can calculate the distance between the points and the closest polyline to each point but I'm not sure how to go about calculating the distance between my points and each individual line. Any ideas?


title: Can I tie-in a data collection form with an online map service to crowdsource values for map features?

text: Apologies if this has been covered before. I've done some research but unfortunately my Google Fu has failed me. I want to be able to collect some simple numeric data from multiple people over the web. I have a 10 x 10 grid of polygons (covering a total of 25 Hectares) and each value that I collect should pertain to one of the 100 cells within that grid. Example of the grid: ![](http://i.stack.imgur.com/v3C5U.jpg) I've managed to work around doing this in the past as I've only had four features over a similar area. All I did then was to export the features as a KML and import them into Google Maps Engine. I then ran a Google maps iframe alongside a google forms iframe and used a simple labelling system to refer the Google form questions to the map. Example: ![](http://i.stack.imgur.com/V02zs.jpg) Of course doing this for 100 features in place of 4 could become a little unwieldy. My question therefore, is whether anyone is aware of a relatively straightforward means of acquiring feature specific data from multiple people?


title: Populating a Value List in one Parameter, from Values Derived from a Second Parameter, Python

text: Using the Validator in the ArcGIS 10.1 tool GUI, can I loop through and List feature classes from a workspace type parameter (ie a file geodatabase) and then go on and list the values of a field in a feature class from that list? I'm trying to do this:               if self.params[0].value and arcpy.Exists(self.params[0].value):         for fc in arcpy.ListFeatureClasses(self.params[0].value):           if fc == "TestFC":             value_set = set(row.getValue("SITE")                           for row in arcpy.SearchCursor(str(self.params[0].value)))             self.params[1].filter.list = sorted(value_set)      where `param[0]` is a workspace type parameter where the user chooses a geodatabase. `param[1]` will be a string type parameter set as a value list. I want to read in a setted value list from all the values in the 'SITE' field so the user can choose which site the tool will further process on. Maybe in other words:   1. User opens tool   2. User chooses a file geodatabase as the first parameter   3. User will then choose a SITE from the list derived from the first parameters testfc feature class, 'SITE' field Thanks, Mike


title: ArcObjects - find nearest feature to a point - has to be FAST

text: I've seen lots of discussion on this, especially using **IIndexQuery.NearestFeature()** on a **FeatureIndex** instance. Something like this:               private IdDistancePair FindNearest(IGeometry g, IFeatureClass fc, string where = "")     {     IFeatureIndex fi = new FeatureIndexClass();     fi.FeatureClass = fc;          //you can optionally set a FeatureCursor on the IFeatureIndex      //to reduce the number of features indexed.     //also, limiting subfields may speed up the indexing operation     IQueryFilter qfilter = null;     if (!string.IsNullOrEmpty(where))     {       qfilter = new QueryFilterClass();       qfilter.WhereClause = where;       IFeatureCursor fcursor = fc.Search(qfilter, true);       fi.FeatureCursor = fcursor;     }          //create the index     fi.Index(null, ((IGeoDataset)fc).Extent);          int fid;     double distance;     ((IIndexQuery)fi).NearestFeature(g, out fid, out distance);          return new IdDistancePair(fid, distance);     }      However, it appears that the FeatureIndex is an in-memory thing unrelated to any existing spatial index on the query feature class. It takes a long time (dozens of seconds) to build. (and the above doesn't work for me anyway, always returning a FID of -1 after a 45s wait). I only need to do one query per layer, repeated for ~5 layers. The docs on the FeatureIndex class suggest that the above approach is expensive in such a circumstance. I need this to complete in say 10 seconds or less. I CAN limit the search to a radius, say a few thousand meters. But there are still likely to be thousands of query features inside that radius, so I'd prefer to avoid any brute-force loop/measure the distance solution in favor of one that uses a spatial index, if possible. Any suggestions?


title: access attributes of one class for use in another class in python add-in

text: I want to get combobox attribute 'dem' and use it for further processing in class ExtentData. But I am not sure how to get attributes from one class to another for processing in pythonaddin. In following code, I want attribute 'dem' from class ComboBoxFrame right after if statement in class ExtentData. Any suggestion would be appreciative.               import arcpy     import pythonaddins     import wx          class ExtentData(object):     """Implementation for ExtentData_addin.tool (Tool)"""         def __init__(self):             self.enabled = True             self.checked = False             self.cursor = 5             self.shape = 'Rectangle'              def onRectangle(self, rectangle_geometry):         """"""              dlg = None         extent = rectangle_geometry              if extent != None:             self.dlg = ComboBoxFrame()               class ComboBoxFrame(wx.Frame):         def __init__(self):             wx.Frame.__init__(self, None, -1, 'Input Data Selection', size=(350, 300))             self.Bind(wx.EVT_CLOSE, self.OnClose)             panel = wx.Panel(self, -1)                  DemList = ['DEM1', 'DEM2', 'DEM3']             self.Dem = wx.ComboBox(panel, -1, "", (15, 30), wx.DefaultSize,DemList, wx.CB_DROPDOWN)             self.Bind(wx.EVT_COMBOBOX, self.OnSelect, self.Dem)                  self.btnApply = wx.Button(panel, label="Apply", pos=(65, 200))             self.Bind(wx.EVT_BUTTON, self.OnClose, id=self.btnApply.GetId())             self.Show(True)              def OnClose(self, event):             self.Show(False)              def OnSelect(self, event):             dem = str(event.GetValue())          app = wx.PySimpleApp()     app.MainLoop()


title: How to reference a number variable in r.mapcalc within QGIS graphical modeler

text: I need to run the following formula in r.mapcalc: (xor(abs(A-B),abs(A-C))>120)*A I would like to however use a user defined number input variable instead of the 120. I already have this variable loaded up in the graphical modeler but Im not sure how to reference it within the formula....


title: Georeferencing a jpeg in QGIS but the scale is way off

text: I've gone through the really helpful walk-through(s) of how to georeference a jpeg in order to make it a layer, but after inputting the coordinates (either by typing in actual coordinates or clicking "From map canvas" and aligning points) and loading the new layer, the scale is way off (and I think the jpeg isn't even remotely in the right place--hard to tell because of how off the scale is). Any ideas?


title: Calculate actual Landsat image corner coordinates to derive azimuth (heading)

text: This is the continuation of an unsolved problem I had (see this thread if you like). I have a Landsat image and I would like to estimate the satellite azimuth for my location. As the satellite follows a near-polar orbit (see Landsat Handboook), the satellite must have an azimuth angle different than 0°, getting bigger as it moves poleward. If I have the coordinates of the corners of my image, I can estimate how the image is tilted against North direction (estimate the azimuth). So, the idea is:   1. to choose one of the two sides of my image (wether the left or right) to get the angle of the normalized direction (see also How to calculate orientation of line segments using open source GIS?)   2. to use the coordinates of the chosen corners (e.g. lower and upper left, as `a` and `b` variables, see below) in the following MATLAB function (inspired gist.github.com/604912)              function d = line_dir( ax, ay, bx, by )          dx= bx - ax;     dy= by - ay;     m= sqrt((dx^2)+(dy^2));     r= atan2(dx/m, dy/m);     d= (-r*180)/(pi);          if d<0         d=-d     end     end      However, the coordinates of the corners in the metadata file refers to the background image, not the real image (see the figure below to understand), so can I find now the real corners of the image, to then find the azimuth angle? ![enter image description here](http://i.stack.imgur.com/xMsbj.jpg)


title: ArcObjects SQL Server Connection not recognizing instance

text: I am trying to make a connection to a SQL Server database with SdeWorkspaceFactory. This is from my Windows 7 x64 Desktop to a Network database. If I use this:               SDEWorkspace = NewObj(esriDataSourcesGDB.SdeWorkspaceFactory, esriGeodatabase.IWorkspaceFactory2)     connectionString = "server=SQLServer;" + \                         "instance=hbffcacqsql01;" + \                         "database=ACQ_HBED_PROD;" + \                         "authentication_mode=OSA"          result = SDEWorkspace.OpenFromString(connectionString, 0)      or this:               connectionProperties = NewObj(esriSystem.PropertySet, esriSystem.IPropertySet)     connectionProperties.SetProperty("SERVER", "SQLServer")     connectionProperties.SetProperty("INSTANCE", "hbffcacqsql01")     connectionProperties.SetProperty("DATABASE", "ACQ_HBED_PROD")     connectionProperties.SetProperty("AUTHENTICATION_MODE", "OSA")          SDEWorkspace = NewObj(esriDataSourcesGDB.SdeWorkspaceFactory, esriGeodatabase.IWorkspaceFactory2)          result = SDEWorkspace.Open(connectionProperties, 0)      I get this error:               Traceback (most recent call last):       File "C:/Users/tstanley/Downloads/ArcMap_Python/FirstTry.py", line 88, in <module>         result = SDEWorkspace.OpenFromString(connectionString, 0)     _ctypes.COMError: (-2147216021, None, (u'Entry for SDE instance not found in services file', u'esriDataSourcesGDB.SqlWorkspace.1', u'esri_csGeoDatabase.hlp', 0, None))      If I leave out the instance parameter, the connection window pops up and if I put the same value into the window, I get a connection. ![enter image description here](http://i.stack.imgur.com/Gcx94.jpg) Can someone tell me what I'm missing? I've read a bunch of other posts of people with this error and they mention replacing the instance with a service number, 5151. It seems that they are trying to connect to an SDE database, which this isn't. I'm really new to ArcObjects so I have just about everything still to learn. Any help is appreciated.


title: How do I use variables inside a arcpy toolbox function?

text: I have a point file and I need to make a lot of buffers in specific increments but I can't figure out how to interate through a list inside the buffer function. Here is what I have:               import arcpy       city_pts = "C:\Urban_Land_Project\Cities\WI_cities.shp"     buffer_sizes = [1,2,4,7,10,15,20,30]     for x in buffer_sizes:         arcpy.Buffer_analysis(city_pts, "C:\Urban_Land_Project\City_buffers\WI_city_buffers*km.shp","2 Kilometers" ,"FULL","ROUND","NONE"," #")      I need to have each buffered file to be named "WI_city_buffers*km.shp" with the * being the buffer size. How do I get the arcpy.Buffer_analysis to iterate through my list and apply each buffer distance and proper name?


title: R plot raster with real colors

text: I'have a raster image from Scran25 IGN. In Qgis, it opens perfectly but in R (with rgdal or raster packages), I do not know how to keep the original colors like this exemple? ![IGN scan25 from Banyuls](http://i.stack.imgur.com/CcBl7.jpg) I have tried               library('rgdal')     library('raster')          scan25<-readGDAL("./Banyuls/MyCreate/etudehydoAOC/etudehydoAOC/Scan25/vrt_4communes.vrt")     communes<-readOGR(dsn = "./Banyuls//MyCreate/Surface_4communes/",layer = "4communesAOC")          scan25<-raster(scan25)     test<-aggregate(scan25,10,mean)          plot(communes)     plot(test,add=T)      And the result is ![IGN scan25 import in R with bad colors](http://i.stack.imgur.com/Jd7Wh.png)


title: Erdas grids names to field names in tables renamed when saved

text: I am trying to use zonal stat on several imagine rasters and save values to a table with several automatically added fields that will have the same name of the corresponding raster, I can do it but my problem is when the name of imagine raster is added to the table as anew field name the .img is converted to _img so the name of the field is no longer the same and calculate field value breaks since it gets the field name from NAME in the iterate raster tool which will have the .img, if i have my input grids in esri grid format it should work but I want to keep it in Erdas format , is there anyway to do that, thanks. raster


title: Creating model that simulates water flowing down crack, then splitting at 'Y' joint?

text: Maybe this is a simple problem, but I've been given a task of simulating how water flows in one line, along a crack, and then hits a 'Y' junction, where half the water flows in each of the split directions. I'm having trouble imagining how to do this. I had thought about just creating 2 routes, that diverge at the 'Y' joint. Is there a better way to do this?


title: Geodjango- get the distance between 2 points

text: Sorry to be asking yet another geodjango newbie question. What I want to do is quite simple, I have 2 locations defined by gps coords, lat / long like returned by google maps: http://gmaps-samples.googlecode.com/svn/trunk/geocoder/singlegeocode.html I need to calculate the distance between them. I know I could use the google API but I'll be processing bulk queries so I'd rather do it on my own server. I've spent a few hours with the docs, installed geodjango OK, but I can't seem to find an example of this. Everything in the literature is way more complicated than I need. Thanks for any feedback


title: Implementing OpenLayers3 basemap from tiles created using gdal2tiles

text: I recenty discovered the gdal2tiles.py tool which is working great with OpenLayers 2.x, but I wanted to use the new release of OL and I'm having some troubles. I have googled for a solution but i did not found a way to go. Here is my code, inspired from the solutions offered here: https://groups.google.com/forum/#!searchin/ol3-dev/TMS/ol3-dev/JLn0lfHvI7w/q7C3O8_b0-UJ (still referring to the dev version of OL3, different it seems to the release).               var tileUrlFunction = function(tileCoord, pixelRatio, projection) {       if (!tileCoord) { return ""; }            var x = tileCoord[1];       var y = tileCoord[2];       var z = tileCoord[0];            if (x < 0 || y < 0) { return "";}            return 'http://localhost/map-tiles/'+ z.toString()+'/'+ x.toString() +'/'+y.toString() +'.png';     }          var layers = [       new ol.layer.Tile({         source: new ol.source.XYZ({           tileUrlFunction: tileUrlFunction         })       })     ];          var map = new ol.Map({       target: 'map',       layers: layers,       renderer: 'canvas',       view: new ol.View({         center: [2900,-3900],         zoom: 3       })     });      I've tried several combinations of the `x`, `y` and `z` assignment without success. Tiles on the browser are completely wrong and or missing. Any help would be much appreciated.


title: Whitebox is not finding my java

text: I was recently presented to the software Whitebox GAT 'Iguazu'. Unfortunately, I tried to use it in two different machines: one with Windows 7 (64) + JRE 7 (65) and other with W7 (64) + JRE 7 (67) . In both computers I could not start the program. I tried double click and also open specifying the program path. Nothing happens. Am I doing something wrong?


title: Save map canvas as image with drawings

text: Is there any way to save the map canvas as an image, with the drawings on map? I draw a line on map, and saved it as image, but only the map saved(without line). I used the following code.. iface.mapCanvas().saveAsImage('/home/kasun/output_canvas.png', None, 'PNG') Please help me to find a solution to this problem. Thank You


title: Chloropeth Map in R, Data on map not represented as in listed in the file

text: Below is the code for creating the Chloropeth Maps, once created, the actual problem is that the heat intensity indicated by the map is not according to the dataset column called **value** mentioned in the **fill** property of the MAP code. I do not understand why that is so! I am attaching the Map for better understanding, the scale is actually not correct and neither the map. Scale and Mapping are opposite to the data. note: This is the link from where i downloaded the shape file zip file (http://www.gadm.org/country)               #Data     1. The Shapefile has all columns such as **lat, long,order,group, id** etc     2. There is a **indiamap csv** file listing state against every id     3. I merged to make a new file (**merge.map**) with some of my values and **indiamap csv**.     4. Finally as the code mentioned I merged **merge.map** and **india.map** to get a final file.               # Reading Shapefile         india.shp<- readShapePoly("IND_adm/IND_adm1.shp")         draw.india<- fortify(india.shp)              #Reading data for the corresponding states.         india.map<-read.csv("indiamap.csv", header = TRUE,sep=",")              #merging the state file with main file with lat and long         indiamap1<-merge(x=merge.map,y=india.map,by="id")              #Sorting the dataframe         indiamap1<-arrange(indiamap1,group,order)              #Creating the MAP         ggplot(indiamap1, aes(x=long, y=lat, group=group)) +                 geom_polygon(aes(fill=value),colour="white") +                 scale_fill_gradient2(low="337", mid="87", high="83",                                      midpoint=median(indiamap1$value)) +                 coord_map("polyconic")      ![enter image description here](http://i.stack.imgur.com/rnTN3.jpg)


title: Does QGIS attribute edit form support "relationship" data?

text: Consider I have a geometry table having one-to-many relationship with a no- geometry table. When I create a feature (geometry), I would like to be able to add the relative no-geometry table data in the corresponding feature attribute edit form as well. Is this possible?


title: OpenLayers: WMSGetFeatureInfo shouldn't return feature when a layer is hidden

text: In OpenLayers, I have a map with several layers. I created a control to query the layers and return feature info.                var info= new OpenLayers.Control.WMSGetFeatureInfo({          url:'http://onesharedocean.org/geoserver/ocean/wms',          title:'select polygon',          output:'features',          infoFormat:'application/vnd.ogc.gml',          format: new OpenLayers.Format.GML,          eventListeners: {             getfeatureinfo: function(event) {                console.log(event);                if (typeof(event.features[0])=='undefined'){                   console.log('nada');                   return;                };             roiName=event.features[0].attributes.name;             if (typeof(roiName)=='undefined') {roiName=event.features[0].attributes.Name}             thisPolygon=event.features[0].gml.featureType;             document.getElementById('infoblock').innerHTML='fid: '+thisPolygon+'  Selection: '+roiName;          }          }       });      Is there an option to say that I don't want any information to be returned if the user set the layer visibility to off (from the layerSwitcher), or should I query the layerSwitcher and check if the corresponding layer is off or on (and how can I do that)?


title: how can we convert vector data to raster data in geotools

text: In my project iam able to open and save shape file. I want to convert this vector data to image format as a tiffimage plese help me the procedure to do the same.


title: Accuracy problems with Openlayers\Openstreetmap\AnimatedCluster

text: I am working on a OSM map that clusters several points with a modification of AnimatedCluster (https://github.com/acanimal/AnimatedCluster). However, somewhere the accuracy seems to lack. Can someone help me figure out at which point it might go wrong? The image below displays how a cluster ends up in the middle of the sea, as well as the cluster's longitude/latitude coordinates. The points are still in the sea when I zoom in ![Wrong clusters and their coordinates](http://i.stack.imgur.com/9Fup9.jpg) The code I used to display them goes as follows:               var coor=clusters[i].geometry.getBounds().getCenterLonLat().transform(map.getProjectionObject(),fromProjection);     coords = coords +coor.lat+","+coor.lon +"; ";      If I type in the coordinates manually into Google Maps or Openstreetmap, none of the coordinates are in the sea. However, they are all in different places using Google and OSM. Could it be OSM and/or Google can't handle more than 1-2 decimal places? Or could it be an issue with the AnimatedCluster? I couldn't find anything while debugging there. Parsing between geometry and long/lat seems to cause inaccuracy (the original latitude was 21) on the 13th decimal place, but that should not be visible on the map. This is the code for adding the features:               map = new OpenLayers.Map('geoMap');     osm = new OpenLayers.Layer.OSM();     map.addLayer(osm);          vector = new OpenLayers.Layer.Vector("Features", {     renderers: ['Canvas','SVG'],     strategies: [         new OpenLayers.Strategy.AnimatedCluster({             distance: 50,             animationMethod: OpenLayers.Easing.Expo.easeOut,             animationDuration: 10,             clusterType: clustertype         })         ],         styleMap:  new OpenLayers.StyleMap(style)     });          map.addLayer(vector);          for (var i = 0; i < dataselect.length; i++) {     var lon = dataselect[i][5];     var lat = dataselect[i][6];     lon = Math.round(lon * 100000) / 100000;     lat = Math.round(lat * 100000) / 100000;     var select = getRandomArbitary (0, 4)        var lonlat = new OpenLayers.Geometry.Point(lon,lat);     lonlat.transform(fromProjection, map.getProjectionObject());          var f = new OpenLayers.Feature.Vector(lonlat, null, null);     features.push(f);        }     vector.addFeatures(features);      I'd appreciate any help whatsoever :)


title: ST_Intersects incorrectly returning false in PostGIS

text: I have a shapefile created from ArcMap that contains about 400 geometry's. I've added this shapefile into PostGIS using the shp2pgsql tool. I am trying to use the ST_Intersects function to determine if a given point, using Lat and Long, is within one of the 400 geometry's. I have a point that I know to be within one of the shapes, but ST_Intersects is returning false for this. The record where cocnum equals "MA-500" contains the point Lat: 42.3581, Long: -71.0636.                                   SELECT                              cocname,                             cocnum,                              geom,                                        st_intersects(geom, point(-71.0636,42.3581)::geometry)                         FROM                              public.coc                         WHERE cocnum like 'MA-500'      However, I get the following result: ![Query Result](http://i.stack.imgur.com/nAR7J.png) When I create a KML of the geometry using `st_asKml(geom)` and then insert that into google maps, along with the point `google.maps.LatLng(42.3581, -71.0636)` I get the following result: ![map](http://i.stack.imgur.com/iIOf8.png) Any ideas on why ST_Intersects is returning false for this? Any guidance given is greatly appreciated, Thanks! Ed


title: How to extract pixel value counts from a raster in QGIS?

text: How does one go about obtaining a summary of the counts of cell values within a raster in QGIS? I know QGIS is doing this, in order to provide the summary statistics and histogram in the properties of the layer. I do not need the summary statistics nor a histogram image, but rather counts of individual values. It looks like there was at one point an attempt to add a special feature for this (Feature #4321, http://hub.qgis.org/issues/4321). Am I missing something really obvious?


title: Shapefile for India which includes Telangana (state)?

text: I need a shapefile for **India** which is up-to-date to display the regional borders. By up-to-date, I mean that it **includes Telangana** as a state. Many shapefiles that I downloaded do not have Telangana as a state as it was partitioned a month ago, it is displayed as part of Andhrapradesh. So I need the latest shapefile.


title: Converting shapefile to KML using JavaScript or PHP?

text: My client needs to generate maps from shapefiles. Through my initial search, I came to know that we need to convert shapefiles to .kml then we can generate maps from google maps. So I searched for a script but could not find anywhere regarding that. So please help me with this problem. I am new to this, tell me am I going in right way? One more thing, can we generate shapefile lines on a jpg image on the right place if I take an image of a state?


title: Sending e-mail on feature update in ArcGIS JavaScript application?

text: We have an ArcGIS feature service that users can edit to change a "project status" attribute. An example would be a change from "pending" to "in- progress". Certain groups need to be notified when that attributes new value impacts them. For instance, a change to "in-progress" would let accounting know that charges may start coming in to the job. I'd like to e-mail the appropriate parties using the user's default email client but I'm not sure how to handle the attribute change event. I am new to the JavaScript API and at this point any information will be useful. I've investigated the ArcGIS JavaScript API guide related to event handling but would benefit from a more specific example. If my question needs clarification please let me know.


title: Postgres geometry in a pdf file with php

text: Is there a way to draw a postgres geometry in a pdf file using a php script? Thank you.


title: How to implement a spatial query on polygon data for a specific point?

text: I have vector data fed into a database which returns a report at province level upon querying. This is implemented in a web tool. I want to add an extra functionality at point level, where I can input lat/lon and run a spatial query to calculate some attribute, say _distance to coast_ , based on that point. Let us assume I have data which contains the distance from centroid of the province. Please suggest me an open source approach to implement this. How feasible would it be to divide the province into grids and dynamically take up the values?


title: How to crop Shapefiles in QGIS?

text: I need to clip a large and very detailed .shp file in QGIS showing environmental data similar to a temperature map. Is there something like the very simple crop tool alike in photoshop where you can drag a rectangle over an area of interest and crop the shapefile to this area only? I do have an additional shape file of neighborhoods/districts that I could also use for cropping but I do not know how to specify "extract a shp/kml only within the borders of district X, row xxxxx in the attribute table of districts.shp" How can I extract the data primitively within a rectangle of choice or per district?! I need smaller pieces from the large map that has been provided by an environmental geo-institute to be implemented into an mobile Application (Android) to view the quite heavy data in a lighter way. I made some screenshots of the data and posted them here: http://merglindev.wordpress.com/2012/10/22/extracting-districts-from-big- shapefile-data/. You can see the districts and the data layers. I need to have every district as a separate file / separately export every dataset per district to kml.


title: ArcObjects/ArcServer 10.1 - How to programatically publish a map service?

text: I am migrating some ArcObjects C# code to 10.1 and am having some trouble programatically publishing an MSD. We are using a multi-site configuration (not sure if that matters). Here is the code I am using to publish. I have added some comments to help clarify.               IAGSServerConnectionFactory3 pConnectionFactory = new AGSServerConnectionFactoryClass();          //admin connection file works fine, tested in ArcCatalog on the same server     IAGSServerConnection4 server = pConnectionFactory.OpenFromFile(@"e:\admin-connection.ags", 0) as IAGSServerConnection4;          IAGSServerConnectionAdmin pAGSServerConnectionAdmin = server as IAGSServerConnectionAdmin;     IServerObjectAdmin pServerObjectAdmin = pAGSServerConnectionAdmin.ServerObjectAdmin;     IServerObjectConfiguration4 pConfiguration = (IServerObjectConfiguration4)pServerObjectAdmin.CreateConfiguration();          //Set the general configuration settings     pConfiguration.Name = "output_56";     pConfiguration.TypeName = "MapServer";     pConfiguration.Description = "Map Output"; // required field in 10.1? Read that somewhere and added     pConfiguration.IsPooled = true;     pConfiguration.MinInstances = 1;     pConfiguration.MaxInstances = 15;     pConfiguration.WaitTimeout = 60;     pConfiguration.UsageTimeout = 600;     pConfiguration.IdleTimeout = 1800;     pConfiguration.StartupType = esriStartupType.esriSTAutomatic;     pConfiguration.IsolationLevel = esriServerIsolationLevel.esriServerIsolationAny;          //Set the configuration properties of the MapServer     IPropertySet pProps = pConfiguration.Properties;          //I have confirmed that the MSD exists and the FilePath is correct     pProps.SetProperty("FilePath", mapDocument);                    //output directory is at network path on primary server on site \\10.7.1.XX\arcgisserver\directories\arcgisoutput\UserSims     // it should be noted that an empty directory is created here before the program crashes     pProps.SetProperty("OutputDir", outputDir);           //serverHost = "localhost", which is a secondary ArcServer on the site     // I commented this out because it looks like the arcgisoutput directory does not exist anymore.     // I'm not sure if this is required to publish, since we set OutputDir     //pProps.SetProperty("VirtualOutputDir", "http://" + serverHost + "/arcgisoutput");          pProps.SetProperty("MaxImageHeight", "2048");     pProps.SetProperty("MaxRecordCount", "1000");     pProps.SetProperty("MaxBufferCount", "100");     pProps.SetProperty("MaxImageWidth", "2048");     pConfiguration.Properties = pProps;          //Set the info segment properties     IPropertySet info = pConfiguration.Info;     info.SetProperty("WebEnabled", "true");     info.SetProperty("WebCapabilities", "Map,Query,Data");     pConfiguration.Info = info;          //Set the recycle properties of the MapServer object     IPropertySet recycle = pConfiguration.RecycleProperties;     recycle.SetProperty("StartTime", "1:00 AM");     recycle.SetProperty("Interval", "86400");     pConfiguration.RecycleProperties = recycle;          //Add the configuration to the server     pServerObjectAdmin.AddConfiguration(pConfiguration); //crashes here     pServerObjectAdmin.StartConfiguration(service, "MapServer");      Here is the error I am getting.               Unhandled Exception: System.Runtime.InteropServices.COMException (0x80040208): E     RROR: service failed to start, Cannot connect to this server.        at ESRI.ArcGIS.Server.IServerObjectAdmin.AddConfiguration(IServerObjectConfig     uration config)        at Matric.Modsim.SimulationCruncher.Publish.AddMapService(Geoprocessor gp, IA     GSServerConnection4 server, String serverHost, String outputDir, String service,      String mapDocument) in C:\Users\matt.moyles\Documents\Visual Studio 2010\Projec     ts\simulation2\simulation\Publish.cs:line 111      The error does not make sense to me, because I have confirmed the connection file is valid. If anyone could provide any insight or suggestions I would be very grateful! EDIT: I also found this relevant entry in the arcgisserver logs:               <Msg time="2012-06-20T15:13:27,824" type="SEVERE" code="8274" source="Server" process="632" thread="1" methodName="" machine="ROGER-WATERS.TECH.LOCAL" user="" elapsed="">Failed to construct instance of service 'UserSims/output_56.MapServer'. </Msg>


title: Anyway to remove points outlying shape in ggplot2?

text: I am to map a plot of an area and then plot points over it (given by Easting and Northing). However, I found out that there are certain mistakes in the data entered for Easting and Northing resulting in some outlying points, and `ggplot` (which I am using) tends to reshape the map in order to fit the outlying points. Is there any quick way to fix this issue? I do not know of any except to identify and correct the points manually. I am using shape files from the Ordnance survey. Here is the screen-shot. The image does not appear to be skewed. However, upon comparison with other images, the differences become apparent. ** ## EDIT ** The code I'm using is as follows.               countyRegion<- readShapePoly(file.choose())     norfolkCounty<- countyRegion[countyRegion$NAME=="Norfolk County",]     #convert shape file into data that can be plotted on graph     gpclibPermit()     norfolk<- fortify(norfolkCounty,region="NAME")     Norfolk<- merge(norfolk, norfolkCounty@data,by.x="id",by.y="NAME")          #eit data     jul2012<- read.table(file.choose(),header=TRUE,as.is=TRUE,blank.lines.skip=FALSE,sep=",")     jul2012$Crime.type<- factor(jul2012$Crime.type)     jul2012<- jul2012[jul2012$Crime.type!="",]     jul2012$Crime.type<- factor(jul2012$Crime.type)     levels(jul2012$Crime.type)<- c("Anti-social behaviour","Theft/burglary","Criminal damage/arson","Drugs","other crime","Theft/burglary","Public disorder and weapons","Theft/burglary","Vehicle crime","Violence and Sexual Offences")          #plot     jul12Map<- ggplot(data=Norfolk,aes(long,lat)) + geom_polygon() + geom_point(data=jul2012,aes(Easting,Northing,colour=as.factor(Crime.type)),alpha=0.6) + scale_colour_brewer(palette="Set1",name="Category") + theme_bw() + scale_x_continuous(breaks=0,labels="") + scale_y_continuous(breaks=0,labels="") + theme(axis.ticks=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank())      The shape file and the crime stats file are linked, as well as the dbf, shx, and prj files. ![enter image description here](http://i.stack.imgur.com/GDdzJ.png) ![enter image description here](http://i.stack.imgur.com/peMDM.jpg) Many thanks.


title: overlay wfs google in openlayers

text: ![enter image description here](http://i.stack.imgur.com/yA1eb.png)                   map = new OpenLayers.Map('map', options);         // BASELAYER                   var google_maps = new OpenLayers.Layer.Google(             "Google Maps", {                 numZoomLevels: 20             }         );              var google_satellite = new OpenLayers.Layer.Google(             "Google Satellite", {                 type: google.maps.MapTypeId.SATELLITE,                 numZoomLevels: 20                  }         );         // OVERLAYS         tiled = new OpenLayers.Layer.WMS(         "Kaski",         wmsUrl,         {             LAYERS: 'KaskiGrup',             transparent: true,             visibility: true,                  },         {             transitionEffect: 'none',             isBaseLayer: false,                   }     );                                       // NOKTASAL LAYERS         wfs_layer = new OpenLayers.Layer.Vector("Point",             {                 /*style: {                     externalGraphic: 'http://www.openlayers.org/dev/img/marker.png',                     graphicWidth: 21,                     graphicHeight: 25,                     graphicYOffset: -24                 },*/                 styleMap: new OpenLayers.StyleMap({                     "default": new OpenLayers.Style({                         pointRadius: 10, // sized according to type attribute                         fillColor: "#00FFFF",                         strokeColor: "#33CCCC",                         strokeWidth: 3,                         graphicZIndex: 1                     }),                 }),                      strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],                      protocol: new OpenLayers.Protocol.WFS({                     version: version,                     // loading data through localhost url path                     url: wfsUrl,                     featureNS: featureNS,                     visibility: true,                          // layer name                     featureType: featureType,                     // geometry column name                     internalProjection: WGS84,                     externalProjection: WGS84_google_mercator,                     geometryName: geometryName,                     schema: schema,                      })             });      I want to add google map on the background and overlay my wms and wfs layers but my wfs layers is not at correct location, also wms layers is at correct location. As you can see, i add google prejection on map option where is my mistake ? Thanks


title: Find Nearest Line Segments to Point

text: **Background** From a known Point, I require to establish the nearest surrounding "visible perimeter" against a table of MultiLineStrings, as shown on the diagram. I've searched this site with a number of terms (e.g. minimum edge, minimum perimeter, nearest neighbour, clip, containing polygon, visibility, snap, cut nodes, ray-trace, flood fill, inner boundary, routing, concave hull) but cannot find any previous question which seems to match this scenario. **Diagram**   * The green circle is the known Point.   * The black lines are the known MultiLineStrings.   * The grey lines are an indication of a radial sweep from the known Point.   * The red points are the nearest intersection of the radial sweep and the MultiLineStrings. ![enter image description here](http://i.stack.imgur.com/OOALL.jpg) **Parameters**   * The Point will never intersect the MultiLineStrings.   * The Point will always be nominally centred within the MultiLineStrings.   * The MultiLineStrings will never fully enclose the Point, therefore the perimeter will be a MultiLineString.   * There will be a table containing approximately 1,000 MultiLineStrings (normally containing a single line of about 100 points). **Methodology Considered**   * Undertake a radial sweep by constructing a series of lines from the known Point (at, say, 1 degree increments).   * Establish the nearest intersection point of each radial sweep line with the MultiLineStrings.   * When one of the radial sweep lines doesn't intersect with any of the MultiLineStrings, this would indicate a gap in the perimeter which would be accommodated in the perimeter MultiLineString construction. **Summary** Whilst this technique will find the nearest intersections, it won't necessarily find all the closest perimeter node points, dependent on the resolution of the radial sweep. Can anyone recommend an alternative method to establish all the perimeter points or supplement the radial sweep technique with some form of buffering, sectoring or offsetting? **Software** My preference is to use SpatiaLite and/or Shapely for the solution but would welcome any suggestions which could be implemented using open source software. **Edit: Working Solution (based on answer by @gene)**               from shapely.geometry import Point, LineString, mapping, shape     from shapely.ops import cascaded_union     from shapely import affinity     import fiona          sweep_res = 10  # sweep resolution (degrees)     focal_pt = Point(0, 0)  # radial sweep centre point     sweep_radius = 100.0  # sweep radius          # create the radial sweep lines     line = LineString([(focal_pt.x,focal_pt.y), \                        (focal_pt.x, focal_pt.y + sweep_radius)])          sweep_lines = [affinity.rotate(line, i, (focal_pt.x, focal_pt.y)) \                    for i in range(0, 360, sweep_res)]          radial_sweep = cascaded_union(sweep_lines)          # load the input lines and combine them into one geometry     input_lines = fiona.open("input_lines.shp")     input_shapes = [shape(f['geometry']) for f in input_lines]     all_input_lines = cascaded_union(input_shapes)          perimeter = []     # traverse each radial sweep line and check for intersection with input lines     for radial_line in radial_sweep:         inter = radial_line.intersection(all_input_lines)              if inter.type == "MultiPoint":            # radial line intersects at multiple points            inter_dict = {}            for inter_pt in inter:                inter_dict[focal_pt.distance(inter_pt)] = inter_pt            # save the nearest intersected point to the sweep centre point            perimeter.append(inter_dict[min(inter_dict.keys())])              if inter.type == "Point":            # radial line intersects at one point only            perimeter.append(inter)              if inter.type == "GeometryCollection":            # radial line doesn't intersect, so skip            pass          # combine the nearest perimeter points into one geometry     solution = cascaded_union(perimeter)          # save the perimeter geometry     schema = {'geometry': 'MultiPoint', 'properties': {'test': 'int'}}     with fiona.open('perimeter.shp', 'w', 'ESRI Shapefile', schema) as e:          e.write({'geometry':mapping(solution), 'properties':{'test':1}})


title: Is it possible to insert note to remember the info of a project?

text: Sometimes I need to open an old project, and mostly I remember how it works (i.e. the structure, the formatting and so on) but sometimes not. At other times I open a project composed by a colleague of mine and it is hard to understand how the project was conceived. Is it possible to insert in the project a non-printable note to help to understand the project more quickly? edit: I use QGIS 1.8 (2.0 soon) on windows 7


title: Python code to determine the position on a map using Geo-Tagging

text: What is a proper python code, which enables to scan photographs (when captured from a series of snapshots) and match it with a geo-tagged image database. Based on the extent of similarities existing between captured photo and geo- tagged images, it must help in finding my position on a map. So far, I am aware of finding latitude and longitude of an image file using following python code               from osgeo import gdal          geotiff= gdal.Open('C:\Users\nikhil\Desktop\north_rgb.tif')          width,height=geotiff.RasterXSize,geotiff.RasterYSize          gt = geotiff.GetGeoTransform()          minx = gt[0]          miny = gt[3] + width*gt[4] + height*gt[5]          maxx = gt[0] + width*gt[1] + height*gt[2]          maxy = gt[3]          print(minx,miny,maxx,maxy)      I need some extension or modification of this code as my requirement stated above.


title: Calculating distance between two sets of Decimal Degree coordinates?

text: I have coordinates in decimal degrees e.g: point A:Lat,Lng(5.317101667,-2.023303333) and point B:lat,lng(5.317093333,-2.023318333). The altitude in Meters at Point A is 11m and the altitude and Point B is 33m, How can i accurately calculate the distance between these points? Note I know about Haversine, howeve,r these points are close to each other and I want an accurate way to calculate the distance between these two, would finding the distance between the points with d= _sqrt_ (x2-x1)^2+(y2-y1)^2+(Z2-Z1)^2, be accurate? if so what units would the answer be in? I am looking to get the answer in metres.


title: Devising measure of vertical distance from spring line?

text: I am building a model (using Logistic Regression) to explore the relation between a binary categorical dependent variable (parcels agricultural production quality) and several predictors, among which slope, curvature, distance from the coastline, etc. One of the predictors I wish to include in the model (the significance of its contribution pending) is distance from spring line. I have digitized the spring line as a polyline. The spring line is usually rather high in elevation (in the territory under study) and it lies where two geological layers meet up. Now, I would like to devise a distance measure which allows to discriminate between areas below the spring line and above it. The rationale of it would be that, due to gravity, areas (and related parcels) above the spring line could possibly be more difficult to irrigate, whereas areas (and related parcels) below the spring line would be easier to irrigate. I hope my poor English manages to make my point.


title: Offline entire planet in Java Desktop app

text: Sorry if this question has already been asked, but I need your help. In program that we are developing we need to integrate maps. Here are the requirements that are put for integration:   1. Maps should be off-line, no access to online   2. Must be present all world, so the user can go anywhere as on website http://www.openstreetmap.org/   3. Maps must to be integrated(embedded) into the program(no third-party systems install)   4. Maps data should be in reasonable aisles   5. Maps need to work on the user's computer(not server configuration) I'll add that we at this moment interested in the rapid integration of(a deadline) and again, our program is developed under the Java for desktop. Also we interested in free solutions. As I understand it we need three solutions:   1. Database that can contain info about entire planet(from planet.osm file)   2. Lib that can render(on the fly) a piece of info from database into the tiles   3. Lib that shows renders results to user Please help in choosing solutions/libraries and specify the way in which these things have already been made.


title: Migrating from GeoMedia to ArcGIS for Desktop?

text: I have just started a new job and one of my first tasks is to swap my section from GeoMedia to ArcGIS for Desktop. Has anyone completed this task before? Are you able to share advise, strategic plan, roadmap etc???


title: MODIS calculations with non cloud contaminated pixels

text: I have 11 files showing 8day daytime LST temps in .tif format from MODIS. Is there any way to calculate the average 8day daytime temperature for march 2013 for example but using in the calculations only pixels which are not cloud-contaminated in any of the 11 files that participate in this process?


title: Storing a SketchUp file in PostGIS

text: I am sorry if this is a relatively dumb question, but I was wondering if there is a way to store a SketchUp file (.skp) in a PostGIS database? I was thinking about as a blob, but I wasn't sure if anyone has had success with this. The record would be a point location. Right now I am doing a file path link to the file on the directory, but I was hoping there was a way to store the file in the database. I would like to use this for data storage, but also as a means to retreive and view the data as well. Thanks...


title: Is Hillshade option in ArcMap rubbish or is it me?

text: Here is a pretty pic of the raster I am using - with no effects applied - just the color ramp and historgram adjusted to show out all the relief - all be it in a very 2D way! ![enter image description here](http://i.stack.imgur.com/YKejh.png) I have a raster layer that I want to display with hillshade as a default. I can do this in code and in effect it is the same as selecting the Hillshade option in ArcMap (see below): ![enter image description here](http://i.stack.imgur.com/gKjds.png) When you do this however - it looks rubbish: ![enter image description here](http://i.stack.imgur.com/bHyA2.png) Especially if you zoom into the raster you see all kinds of gridded articfacts - which does not look nice at all! ![enter image description here](http://i.stack.imgur.com/Wbhej.png) However if I create a proper Hillshade layer using the raster with Spatial Analyst and overlay it with transparancy of 85% it looks pretty good and very different!: ![enter image description here](http://i.stack.imgur.com/4EbRG.png) My question is in two parts really: 1) How can I acheive this effect in ArcMap and override the default functionality of the Hillshade option (I supsect this is impossible, or is it some setting I have missed out that is caused the raster to be displayed poorly). If it is possible to override or select the correct settings - is it possible to do in ArcObjects. 2) If option 1 is not possible, can I create my effect with creating a Hillshade layer (without using Spatial Analyst) and effectively fuse it to the orginal raster - i.e. I don't want people to see two layers in the table of contents in ArcMap. I just want one layer but with a default hillshade rendering of my choosing/configuration. I need to do this programmatically too i.e ArcObjects. **UPDATE 13 Jan 2012** OK - so thanks to everybody on here with their thoughts, some nice suggestions too and things I never knew existed. One thing I did learn is that you must have your DEM in the correct Spatial Reference and measurement units or the resulting hillshade looks too black and over exaggereated - check out this video to find out more - its a very good explanation. So I think I now know what I need to do in regards to producing some high quality hillshading, but I think it's going to take some low level tweaking ArcObjects code. I am thinking I need to come up with my own graphics rendering and plug this into ArcMap (I have a dll that renders bitmaps of grids in a cool way). I think a good place to show this would be to override the Draw() event on something like the IRasterRenderer, and insert my own image rendering. Of course it might get a bit complicated with reprojections etc. Hmmmmm. Has anyone done any kind of this low level graphics rendering in ArcMap - if so can you offer any guidance or things I should be aware of?


title: How to measure area correctly?

text: I have now used my entire evening, but I can't figure out how to make areal calculations (eg. in km2). I have tried the 'field calculator', I have tried the 'measure area' function on the "attributes" toolbar, and i have tried fTools in the "Vector"-menu.      How I do:   1\. Open Qgis   2\. Add Google Hybrid Layer (OpenLayers Plugin)   3\. Find tenniscourt (should be 200-260m2)   4\. Measure   \- attribute toolbar->measure area - result: 65.060 m² hmm?   \- layer->new shapefile layer->polygon, save, "add feature" and drawing polygon. Then field calculator $area - result: 118.48828125   \- repeat with new shapefile - result: 4.05503328693158e-08   \- vector->geometric tools->export/add geometric colums = another random result   I can't figure out what I am doing wrong! For me it seems like there is some error in the program? I have tried saving the shapefile and opening it again in a "blank" project - same result - and new result everytime?      I do not understand how this can be so difficult, for me is seems like a pretty basic function, but I am completely lost here.      looking forward to hear from you!


title: Creating dynamic GIS charts in ArcGIS 10.1

text: I have generated a table which has road segments which get inundated during different flood events. I want to be able to add a pie chart to the map which shows that during each flood ARI (1y, 2y, 5y, 10y, 20y, 50y and 100y) x% of the road network in the AOI gets inundated. I tried the following - ![enter image description here](http://i.stack.imgur.com/laGH0.png) and get ![enter image description here](http://i.stack.imgur.com/dOkJf.png) The table looks like ![enter image description here](http://i.stack.imgur.com/lsV3f.png) Is it possible for me to add python code for the chart to add all the 1y lengths and display that number only, 2y total lengths etc only so that you get a pie chart with just the 7 categories. I need to automate this as it needs to be run in several areas and some of the input data changes. best,


title: In osm2pgsql how is the planet_osm_roads table populated?

text: I'm using osm2pgsql to populate a North America OSM database from planet.osm. I like how simple the resulting PostGIS schema is, but I'm a bit confused by its contents. I have the two tables `planet_osm_roads` and `planet_osm_line`. I assumed that planet_osm_roads contained all the highway data, but then discovered that it only contains a fraction of the roads in some of the 'highway' types. For example, planet_osm_line has over 12,000,000 lines with a 'highway' value of 'residential'. However, planet_osm_roads has only 2,000. By comparison both tables have the same number of 'motorway' features. Looking at the data I can't tell why some features would be in the one table but not the other. I'd really like to see the rules the osm2pgsql uses to direct line features into each table. I tried looking through the osm2pgsql source code to shed light on the mystery but so far I've not found what I'm after. Any suggestions very welcome.


title: Merging many map sheets together into one large sheet for later georeferencing?

text: I have hundreds of images, each one is a scan of a map sheet, and I want to georeference them. The maps cover a large area (a whole country), they just split the area into many smaller cells. Each map sheet aligns with (but doesn't overlap) another sheet. The sheets are predicably named, so I can "go north" of the current sheet etc. Each scanned sheet has a border (the map image wasn't pushed up against the edges of the scanner). Each sheet also has an internal paper border. Each scanned sheet was not aligned to the scanner, so each sheet will be skewed/tilted a little bit. I know I have to manually crop all of them. I know how to georeference using QGis etc., but it's just the task of doing ~500 individual georeferencing that's daunting. Is there any way I can, once I have the cropped skewed images, 'merge' the sheets together to produce a small number of larger images? The top left corner of the cropped image A, should match up with the bottom left corner of image B, etc., so is there a way I can automatically 'glue' image A and B together so I only have to georectify the final image? I'm not a GIS professional, so apologies if this is a hard to understand question. I'm using Linux, so has to work there.


title: Hand fit curve and then apply to time series of rasters, make a map. How?

text: Maybe a data management problem, not too sure, but likely a challenge which someone has dealt with before. I have a time series of rasters. I can summarise the information by taking averages of each raster and then plotting the averages with time (which is what is shown in the image). I also want to show a spatial representation and produce a map of the interpolated history on the polygon area of a site. ![Fittedcurve](http://i.stack.imgur.com/O2fuM.png) The green data points are the averages from each of the rasters. The fine dashed line is the expected/modelled performance. The red line with square data points is the hand fitted curve of what I "know to be true", something between the model and the measured values, an interpolation. Unfortunately the relationship with time is not amenable to regression, at least not with any of the simplistic tools I use. Note that there is tail on the curve which the model shows but the lack of measurements at that time is ignorant of. I would like the hand fitted curve/model to be modified by each of the rasters so that a final map can be created. The extent of variation is about 0.9 - 1.35 on the max values, this substantial when the whole lot are summed. Experienced workers will recognise the feature being plotted, thats great, but I suspect it would be a common problem in many other disciplines as well. I am looking for solutions I could implement in QGIS and GRASS. cheers.


title: Hotlink vector/raster files to open in current workspace in ArcGIS 10 similar to Mapinfo

text: I want to add datasets to ArcMap and place them in a particular group layer by simply clicking on hotlinked polygons in a shape file. The same process is easy in MapInfo and I want to implement the same for ArcGIS 10. I found the following post http://forums.arcgis.com/threads/8236-Hyperlinks- in-ArcGIS-10 that describes the process in 9.3.1 On properties-->html popup you can load url and any file will open in an external program but I want it to add to the current workspace. A solution similar to mapinfo--layer control--hotlinks will be ideal but a python based solution will be acceptable as well if it can be added to the relevant vector layers within an mxd. Thanks,


title: How to emulate this grayscale basemap look+feel in Google Maps

text: I love the design on the Development Seed homepage, where they've used a grayscale basemap with bold 3D-style graphic markers. ![enter image description here](http://i.stack.imgur.com/ZeowL.jpg) It appears that they've used the ModestMaps API to create a tiled basemap at various resolutions, with the text and graphics burned into the basemap. I'd like to emulate this look+feel (for a personal, just-for-fun project) using freely available tools such as the Google Maps API. I need the ability to change the marker overlays rapidly, so they would need to draw dynamically over the basemap. **Is it possible to achieve a similar grayscale, terrain-based backdrop using Google Maps?** The Google Maps API Styling documentation shows how to reduce the saturation of the topo map, but I can't see an option to do it with the Terrain map. I think I can emulate the bold graphics using the complex Icons sample as a starting point but would welcome any suggestions. Thanks


title: Improving ArcGIS for Desktop geocoding process?

text: I would like to improve my ArcGIS for Desktop geocoding process comparing the city names between two shapefiles . First one Store_city_name (point) and the second is city_name_boundaries (polygons). My method : After using spatial join tool in Arcgis, If the city names from store shapefile are not the same as polygon shapefile, I’m trying to improve by identifying the distance between the store location and the city XY coordinate where it should be. Everything is working well while calculating the distance between the store location and the centroide of the right city. But sometimes this distance is so high specially if the store was geocoded, for exemple, 20 meters outside the edge of the right city. I would like to find the distance and the XY coordinates from the store location until the edge of the right polygon (city feature) - where the store should be. Can anyone help me ?


title: A count of points within a set distance using updatecursor

text: I have a script that calculates the average distance between points. What I want is to rewrite it so I can input the point layer, input a distance, and the output field will tell me a count of the number of points at that distance. Can anyone help? so basically I just want a script that will tell me a count of the number of points from a point at a certain distance.               #Import standard library modules     import sys, string, arcpy              inputFC = arcpy.GetParameterAsText(0)# Input Point Layer         inputdistance = arcpy.GetParameterAsText(1) #input distance         outputfield = arcpy.GetParameterAsText(1) #output field          # Read the point data into list.       #         xs = []         ys = []      rowCursor = arcpy.SearchCursor(inputFC)      for row in rowCursor:         geom = row.getValue (properties.shapeFieldName)         cen = geom.centroid         xs.append(cen.X)         ys.append(cen.Y)      del rowCursor           # compute distance       sumdist = 0.0      count = 0.0      for i in range(0,len(xs)):         for j in range(i,len(xs)):             if (i <> j):                dist = ((xs[i] -xs[j])**2 + (ys[i] - ys[j])**2)**0.5                sumdist = sumdist + dist                count = count + 1       print "The average inter-point distance is " + str(sumdist/count)       arcpy.AddMessage("The average inter-point distance is " + str(sumdist/count))      With the following script. Ive tried to run it as Distance.add_nearby_points_count_field('gauges.shp', 100)                import sys      import arcpy      import math      def add_nearby_points_count_field(inputFC, inputdistance):               # Read the point data into list.     coord_pairs = []     desc = arcpy.Describe(inputFC)     rowCursor = arcpy.SearchCursor(inputFC)     for row in rowCursor:         geom = row.getValue(desc.shapeFieldName)         coord_pairs.append(geom.centroid)          # Add a field called 'cnt' and calculate nearby points     arcpy.AddField_management(inputFC, 'cnt', 'LONG') #i also tried cutting   out this part and creating the "cnt" field manually     rowCursor = arcpy.UpdateCursor(inputFC)     for row in rowCursor:         geom = row.getValue(desc.shapeFieldName)         from_point = geom.centroid         near_pts = 0         for to_point in coord_pairs:             distance = math.sqrt(pow((to_point.X - from_point.X), 2) + pow((to_point.Y - from_point.Y), 2))             if distance <= inputdistance:                 near_pts += 1         row.cnt = near_pts - 1 # Subtract 1 to remove the measurement to itself         rowCursor.updateRow(row)     del rowCursor          if __name__ == '__main__':         inputFC ='guages.shp' #originally its a featureclass but I converted it to a shapefile to try both methods         inputdistance = 100         add_nearby_points_count_field(inputFC, inputdistance)


title: Can you set custom fields for a [Feature] Layer?

text: Is there some way I can "tag" a [Feature] Layer so that I can grab the layer from the service and see in some attribute LocationField = I suppose I could put that information in the description and just parse it out, but I'm looking for a way to see what layers (if any) map to some "object" in my SQL database. I would much rather not have to put that kind of config info in my DB.


title: PostGIS: ERROR: lwpoly_from_lwlines: shell must have at least 4 points

text: I have achieved to close lines in order to build a polygon with:               st_addpoint(st_makeline(geom ORDER BY .cd_noseq), st_pointn(st_makeline(geom ORDER BY cd_noseq), 1))      but then, when I try to make a polygon with st_makepolygon like this:               st_makepolygon(st_addpoint(st_makeline(geom ORDER BY cd_noseq), st_pointn(st_makeline(geom ORDER BY cd_noseq), 1))) as geom      I get the following error:               ERROR:  lwpoly_from_lwlines: shell must have at least 4 points     ********** Error **********          ERROR: lwpoly_from_lwlines: shell must have at least 4 points     SQL state: XX000      Technology: "PostgreSQL 9.2.7 on x86_64-unknown-linux-gnu, compiled by gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-4), 64-bit" "POSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.6" LIBJSON="UNKNOWN" (core procs from "2.1.0 r11822" need upgrade) TOPOLOGY (topology procs from "2.1.0 r11822" need upg (…)" Any ideas will be wellcome. Many thanks in advance. Best regards,


title: Check if an object is empty in mapbasic

text: I have mapbasic code that manipulates the obj-column in a point-table. Some rows doesn't have any geometry so the obj-column is empty. My question is how to check if there is no object? All the methods i found in the mapbasic-reference that givs information about the obj-column, like ObjectInfo(), try to access the object which givs an error since there is none. Is there any way to use SQL to select the rows with non empty obj-columns? Regards


title: How to determine if a mouse click falls within a polygon

text: I have a polygon drawn on top of a layer using the graphicsContainer of the current activeView, then I want to be able to show a custom context menu if the user right clicks INSIDE the polygon. So I have to figure out a way to determine if the mouse click falls within the polygon. I have tried to use the IRelationalOperator2 as suggested in this question Does anyone know a function to determine if a Point exists within a polygon?, but I can't seem to get it working properly. Here is what I have so far:               System.Drawing.Point screenMousePoint = System.Windows.Forms.Form.MousePosition;     IPoint screenPoint = new Point();     screenPoint.X = screenMousePoint.X;     screenPoint.Y = screenMousePoint.Y;          IPoint mapPoint = this.controller.GetMapCoordinatesFromScreenCoordinates(screenPoint, activeView);          // this.currentDefinedArea is a polygon     IRelationalOperator2 relation = this.currentDefinedArea as IRelationalOperator2;          if (relation.Contains(mapPoint))     {         // Show context menu     }      UPDATE: The "GetMapCoordinatesFromScreenCoordinates" method looks like this, and it's one of the snippets that is included in arcobjects:               public IPoint GetMapCoordinatesFromScreenCoordinates(IPoint screenPoint, IActiveView activeView)         {             if (screenPoint == null || screenPoint.IsEmpty || activeView == null)             {                 return null;             }                  IScreenDisplay screenDisplay = activeView.ScreenDisplay;             IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;                  return displayTransformation.ToMapPoint((System.Int32)screenPoint.X, (System.Int32)screenPoint.Y);         }      UPDATE 2: OK, found out what was causing the above code to fail. Using:               System.Drawing.Point screenMousePoint = System.Windows.Forms.Form.MousePosition;      ...gives a different set of X, Y coordinates than using the MouseEventArgs provided by the OnMouseDown method. Don't know why I wasn't using them in the first place...here's how my working code looks:               protected override void OnMouseDown(ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs arg)     {         if (arg.Button == System.Windows.Forms.MouseButtons.Right)         {             if (this.controller.CurrentDefinedArea != null)             {                 IPoint screenPoint = new Point();                 screenPoint.X = arg.X;                 screenPoint.Y = arg.Y;                      IPoint mapPoint = this.controller.GetMapCoordinatesFromScreenCoordinates(screenPoint, activeView);                      IRelationalOperator2 relation = this.controller.CurrentDefinedArea as IRelationalOperator2;                      if (relation.Contains(mapPoint))                 {                     this.controller.CreateContextMenu(ArcMap.Application);                 }             }         }     }


title: What UTM like coordinate system is this?

text: I have this set of coordinates: `4622580.00000, 503749.00000` They should roughly correlate with `51°01'39.4"N 13°44'44.2"E` in Dresden, Germany. At first I thought they might be UTM, but that puts me off the coast of Somalia. Does anybody have an idea? Or maybe they are UTM but with special _options_? I am unfortunately not very knowledgeable with coordinate systems. ![enter image description here](http://i.stack.imgur.com/v19uG.png)


title: What web mapping libraries are suitable for indoor maps?

text: I am looking for a custom mapping library but not to display an actual geographical area but to display a custom floor-plan. I have been looking for a while but most of the solutions that I found are for geographical areas. I would like to customise a few zoom levels (probably 3 levels) and ideally I'm looking for a JavaScript, HTML5 (canvas, SVG) solution. Does anyone have any suggestions?


title: convert from C# to VBA

text: I have a code that wrote with C#.But I need to convert it to VBA and then use from it into GIS. how can I convert it be VBA?


title: How can I iterate Selection by Attributes?

text: I have a large data shapefile for the state of New York. Each shape has a code corresponding to each county, so I want to iterate "Selection by Attribute" for each code, create a layer from all the parcels with the same county code, and export that into a separate shapefile. How would I go about doing this? I am using Arcmap 10.1


title: how to tidy up line intersections

text: I am very new to QGIS so please bear with me. I am using version 2.0. I am trying to tidy up lines that I have collected using a vertex function on my gps. When I zoom in to the map in QGIS it shows node points all over the place where the lines intersect. Is there a way to delete the unwanted nodes and just leave one node at the point where the lines intersect?


title: How can I strip out Z values in PostGIS?

text: How do I get rid of the Z dimension in PostGIS? Specifically, I have some `MultiPolygonZ`s that I want to copy into another table as `MultiPolygon`s. I don't care about the Z values - it's okay to just throw them away. There must be a simple way, but I've been unable to find any functions that do this... Any tips are greatly appreciated!


title: coordinate system for calculating distance at the nation scale

text: I want to use ArcGIS 10.1 to calculate the straight (euclidean) distances from zip codes to their nearest city centers in about 300 U.S. metropolitan areas. I have two point shpfiles: one includes the centroids of all zip codes and the other includes the points that represent city centers. Both of the two shpfiles are in the GCS_North_American_1983, which prohibits me from doing the distance calculation. I have converted the GCS into North America Lambert Conformal Conic to do the calculation. It worked. But I am not sure whether this is the best way to do it. Any advice is appreciated!


I use Google maps in my application but because this service has limitation (2500 request per day for Geocoding) I want to use OSM instead of Google maps. But OSM doesn't cover Iran as well as Google maps. My question is: Can I use Google maps tiles in OSM for Iran? You can compare Google maps tiles with OSM in Iran on the following pictures: First is OSM and second is Google maps. ![enter image description here](http://i.stack.imgur.com/yPgzs.jpg) ![enter image description here](http://i.stack.imgur.com/6dzrL.jpg)


title: Setting Layer Properties Using ArcPy

text: I am writing an ArcPy script that creates a Point Density Raster file, and then changes the following properties in the resulting layer:   1. Symbology Classified Classes = 10 Classification Method = Geometrical Interval Excluded Value = 0   2. Display Resample during display using: bilinear interpolation I know I can set `lyr.symbology.numClasses = 10`. and `lyr.symbology.excludedValues = 0`, and even create a dummy layer with the above properties (symbologyLayer) and `arcpy.ApplySymbologyFromLayer_management(lyr, symbologyLayer)`. Is there a way to do this simply by setting attributes of `lyr.symbology`, do I have to do it in one fell swoop with `ApplySymbologyFromLayer_management`, or is there a better way?


title: Primary key field hidden on attribute table added to QGIS

text: I'm having the following issue with QGIS. Adding a sqlite text table (Layer > Add Vector Table) QGIS hides the primary key field (having name `id`). I understand also QGIS uses this field as unique id, but I need to have it in the attribute view because all my links and joins use this field. I'm not able actually to perform a join with my vector tables. I can neither change the column names, nor create views.


title: Using ST_Affine() and trans_param() functions on vector layer

text: After reading this post How to georeference a vector layer with control points? I decided to give it a try to move/rotation/scale a polygon using Postgis. I have two polygons, `polygon_from` and `polygon_to` Below are their 4 vertices' longtitude-latitude. I intend to move `polygon_from` to `polygon_to` position.                gid |      from_x      |      from_y      |       to_x       |       to_y     -----+------------------+------------------+------------------+------------------        1 | 103.771789705736 | 1.29874218723039 | 103.771893494233 | 1.29880598084823        2 | 103.7718553679   | 1.2987732903611  | 103.771964782989 | 1.29879195095702        3 | 103.771887622999 | 1.29870762819626 | 103.771951679355 | 1.29871997730869        4 | 103.771820808866 | 1.29867710104945 | 103.771879806787 | 1.29873515522057      I copy and pasted a function (below) from this website http://casoilresource.lawr.ucdavis.edu/drupal/book/export/html/532               -- Function: trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision)          -- DROP FUNCTION trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision);          CREATE OR REPLACE FUNCTION trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision) AS $BODY$ DECLARE             cc_row record;                  cc_det float;                  inv_cc00 float;             inv_cc01 float;             inv_cc02 float;             inv_cc11 float;             inv_cc12 float;             inv_cc22 float;          BEGIN             EXECUTE 'SELECT              count(a.to_x) as cc00,             sum(a.to_x) as cc01,             sum(a.to_y) as cc02,             sum(a.to_x * b.to_x) as cc11,             sum(a.to_x * b.to_y) as cc12,             sum(a.to_y * b.to_y) as cc22,             sum(a.from_y) as aa0,             sum(a.from_y * b.to_x) as aa1,             sum(a.from_y * b.to_y) as aa2,             sum(a.from_x) as bb0,             sum(a.from_x * b.to_x) as bb1,             sum(a.from_x * b.to_y) as bb2             from (' || sql || ') a,(' || sql|| ') b WHERE a.gid = b.gid' INTO cc_row;                       SELECT INTO cc_det                      cc_row.cc00*cc_row.cc11*cc_row.cc22 +                      cc_row.cc01*cc_row.cc12*cc_row.cc02 +                      cc_row.cc02*cc_row.cc01*cc_row.cc12 -                      cc_row.cc00*cc_row.cc12*cc_row.cc12 -                      cc_row.cc01*cc_row.cc01*cc_row.cc22 -                      cc_row.cc02*cc_row.cc11*cc_row.cc02;                  SELECT INTO inv_cc00 (cc_row.cc11*cc_row.cc22-cc_row.cc12*cc_row.cc12)/cc_det;             SELECT INTO inv_cc01 (cc_row.cc12*cc_row.cc02-cc_row.cc01*cc_row.cc22)/cc_det;             SELECT INTO inv_cc02 (cc_row.cc01*cc_row.cc12-cc_row.cc11*cc_row.cc02)/cc_det;             SELECT INTO inv_cc11 (cc_row.cc00*cc_row.cc22-cc_row.cc02*cc_row.cc02)/cc_det;             SELECT INTO inv_cc12 (cc_row.cc01*cc_row.cc02-cc_row.cc00*cc_row.cc12)/cc_det;             SELECT INTO inv_cc22 (cc_row.cc00*cc_row.cc11-cc_row.cc01*cc_row.cc01)/cc_det;                  SELECT INTO a cc_row.bb0*inv_cc01+cc_row.bb1*inv_cc11+cc_row.bb2*inv_cc12;             SELECT INTO b cc_row.bb0*inv_cc02+cc_row.bb1*inv_cc12+cc_row.bb2*inv_cc22;             SELECT INTO d cc_row.aa0*inv_cc01+cc_row.aa1*inv_cc11+cc_row.aa2*inv_cc12;             SELECT INTO e cc_row.aa0*inv_cc02+cc_row.aa1*inv_cc12+cc_row.aa2*inv_cc22;             SELECT INTO xoff cc_row.bb0*inv_cc00+cc_row.bb1*inv_cc01+cc_row.bb2*inv_cc02;             SELECT INTO yoff cc_row.aa0*inv_cc00+cc_row.aa1*inv_cc01+cc_row.aa2*inv_cc02;          END;          $BODY$   LANGUAGE 'plpgsql' VOLATILE; ALTER FUNCTION trans_param(IN sql text, OUT a double precision, OUT b double precision, OUT d double precision, OUT e double precision, OUT xoff double precision, OUT yoff double precision) OWNER TO postgres;          SELECT trans_param('select * from above_table');      After applying the PostGIS function, I obtained `a,b,d,e,x_off,y_off` parameters as below               (0.0824142911101262,47.5134245217341,0.00103114836995566,0.594639499427103,-9.93126928180573,-0.124261944539057)      However when I apply those parameters to the `polygon_from`, it goes to a very strange lat/long position and becomes an irregular shape (not what I expected).               Update polygon_table Set the_geom=(Select st_affine(the_geom,0.0824142911101262,47.5134245217341,0.00103114836995566,0.594639499427103,-9.93126928180573,-0.124261944539057) Where gid=polygon_from_gid;      I am wondering how to use the above **trans_param() function** and those **ST_Affine** parameters to shift/rotate/scale a polygon?


title: Shifting image in a vrt

text: I have a set of geotif images that I want shifted by an integer number of pixels in X and Y. However, those images are quite large, so I don't want to duplicate them (and I need to keep the original unmodified), so modifying the geotransform is not a solution for me. Is it possible to create a set of vrts (with different shift values) without modifying or duplicating the original images?


title: OpenLayers 3: "goog is not defined"

text:
First time using OpenLayers. I am getting an error "goog is not defined". Here are the steps I performed:   1. Downloaded and unzipped v3.0.0-gamma.4   2. It created the following directory structure:     * apidoc     * build     * closure-library     * css     * doc     * examples     * ol     * resources   3. I proceed to create a file called ol3.html and place it in the examples folder shown above. Here is the code of that file: `<!doctype html>   <link rel="stylesheet" href="../css/ol.css" type="text/css">   <title>OpenLayers 3 example</title>
My Map

   `   4. When I load up the page, I get the error "goog is not defined". I can see that it finds the local **ol.js** file but it give the error "goog is not defined".   5. If I switch this line: `` To this: `` Then everything is fine and I can see the map. What am I missing?


title: Grid-like artifact visible in contour layers?

text: I am creating contours from an Aster DEM in ArcMap 10.1 I see a grid pattern appearing, particularly in areas of dense contours: ![enter image description here](http://i.stack.imgur.com/29Rfq.jpg) I see the same artifact if I create a Hillshade layer so presumably it is a problem in the DEM. ![enter image description here](http://i.stack.imgur.com/S6Btt.jpg) The grid artifact doesn't correspond to the sample size of the DEM. I am using the original DEM - no resampling or focal statistics. The DEM has not been reprojected (displayed in the WGS84 UTM projection it was delivered in). Any ideas about the cause of the artifact? If so could anyone suggest a workflow (other than resampling/ focal statistics) for removing it?


title: Better way to combine 2 rasters

text: I have two raster layers: One with 56 m resolution and the other has 30 m resolution. Usually, I just use the 'Combine' command in ArcGIS to combine them, getting a single value for each unique combination of values of the 2 rasters. ArcGIS internally controls the resampling, and the resulting raster is 56 m in resolution. Alternatively, one can split the 56m raster into two (getting two 28m rasters as result), resample the 28m rasters into 30 m rasters and then combine with the other 30 m raster. Which way is better in terms of being more accurate? Any citations to this extent would be super helpful.


title: ArcPy inconsistently accepts SQL queries for MakeFeatureLayer or SelectLayerByAttribute

text: I need to buffer a subset of features in a particular feature class. To do this, I am first creating a feature layer (using a SQL query to select only what I want), and then running Buffer. So my ArcPy script includes this little function to buffer my handful of feature classes:               def funcBuffer(inFeature, outFeature, bufDist, querySQL):         arcpy.MakeFeatureLayer_management(inFeature, "tempLayer", querySQL)         arcpy.Buffer_analysis("tempLayer", outFeature, bufDist, "FULL", "ROUND", "ALL", "")      This works perfectly on _most_ of my feature classes. However, there is one that fails. The Python error is:   `ERROR 000230: Failed selecting with "OH_UG" = 'UG'` The help page for that error suggests that the SQL query is invalid. However,   * the field `OH_UG` exists, and there are some features with the value "UG"   * it successfully validates in Query Builder   * I can use it to Select by Attribute with no problem (python snippet from within ArcMap: `arcpy.SelectLayerByAttribute_management("TransmissionLines","NEW_SELECTION","OH_UG = 'UG'")`)   * similarly structured queries work on other feature classes I've tried different ways of passing the SQL query and all fail: `'"OH_UH" = \'UG\''`, `'OH_UH = \'UG\''`, `""" "OH_UG" = 'UG' """` I revised the script to do the selection after making the feature layer:               def funcBuffer(inFeature, outFeature, bufDist, querySQL):         arcpy.MakeFeatureLayer_management(inFeature, "tempLayer")         arcpy.SelectLayerByAttribute_management("tempLayer", "NEW_SELECTION", querySQL)         arcpy.Buffer_analysis("tempLayer", outFeature, bufDist, "FULL", "ROUND", "ALL", "")      and now `MakeFeatureLayer` works, but `SelectLayerByAttribute` throws `ERROR 000358: Invalid expression` \-- still a complaint regarding invalid SQL. I attempted to select using a different (numeric) field in the table, and am still getting the same `000358` error. * * * So I am not sure what's going wrong -- is it the SQL, is it the particular feature class, is it gremlins? How can I work around the problem in ArcPy?


title: Water Flow Analysis

text: If I know that 10,000+ gallons of water were going to be dropped at a single point how could I model where it would flow. I would prefer something that worked in ArcGIS and could create a layer. Thanks!


title: ArcGIS Server - host ArcGIS Services through fully-qualified domain name and servername

text: I’m in need of some direction… We currently host an ArcGIS Server 9.3.1 SP2 .NET on our Intranet. A situation has come up that requires me to have the ArcGIS server respond to both the servername and through a fully qualified domain URL to access GIS services. How do I configure the server to respond to the fully-qualified domain name and servername for ArcGIS requests?               http://servername/ArcGIS/rest/Services          http://sub.domain.foo/ArcGIS/rest/services      Is there a whitepaper, directions, something that explains how I can fulfill this scenario? * * * mwalker - Thank you for the response. Reverse proxy would be ideal and actually after reading through this setup I will be proposing we setup our environment as such going forward. In the mean time I need to have this particular Silverlight site up and running by Monday morning within our Intranet for remote office access (Secure VPN), using a single domain name where the GIS services can be referenced through the domain name itself. Now I know I could bind the domain name to the Default Web Site within IIS, setup the Silverlight project, and GIS services and the site should respond as expected. This is the quick and dirty answer to my dilemma for a Monday morning availability. But for conversation and learning more about the ArcGIS server flexibility - are there any other options to consider?


title: ST_Collect runs out of memory

text: I am trying to calculate the convexhull for an entire table containing 3 200 000 records. the_geom is of type POINT:               SELECT ST_ConvexHull(ST_Collect(the_geom)) FROM mytable;      After 1-2 seconds it fails with the following error:               ERROR:  out of memory     DETAIL:  Failed on request of size 134217728.          ********** Error **********          ERROR: out of memory     SQL state: 53200     Detail: Failed on request of size 134217728.      I have tried removing the ST_Convexhull call and just called ST_Collect with the same result. I am currently running this on PostGIS 1.5, PostgreSQL 9.1 on Windows 8. Is it likely running on Linux will somehow fix the problem? I have tried to set the following values in postgresql.conf, but nothing seem to help:               shared_buffers = 1024MB     temp_buffers = 1024MB     work_mem = 1024MB     maintenance_work_mem = 1024MB      I have also tried using ST_MemCollect but it runs extremely slow.


title: How to print two pagelayout on two sides at the same time ?

text: I am having a PageLayoutControl with name: PageControl1. I used below codes for printing it:               public override void OnClick()      {                     try                     {                         PrintDialog printDiag = new PrintDialog();                         DialogResult res = printDiag.ShowDialog();                         if (res == DialogResult.Cancel)                             return;                              ESRI.ArcGIS.Output.IPaper paper = new ESRI.ArcGIS.Output.PaperClass();                         paper.Attach(printDiag.PrinterSettings.GetHdevmode().ToInt32(), printDiag.PrinterSettings.GetHdevnames().ToInt32());                         paper.Orientation = PageControl1.PageLayout .Page .Orientation ;                         PageControl1.Printer.Paper = paper;                         PageControl1.PrintPageLayout(1, 1, 0);                          }                     catch (Exception ex)                     {                         MessageBox.Show(ex.Message , null, MessageBoxButtons .OK , MessageBoxIcon.Exclamation );                     }           }      And I have another PageLayoutControl with name: PageControl2. Do you know how to print these PageLayoutControl on two sides of one paper (only one click printing button) ? (PageControl1 on the front side, PageControl2 on the back side) ? Thanks and regards,


title: How to know if flow accumulation map reflects true catchment areas?

text: I have calculated flow accumulation (img file) using this DEM (img file), which is 0.5 degree resolution, using arcmap. I am wondering if the cell accumulation value can directly be translated into drainage areas in square kilometers or not, and if these values reflect true catchment areas? For example, if the flow accumulation value is 140 and my resolution is 0.5 * 0.5 degree = 56km^2, the upstream area is 140*56 = 7840km^2. Any help please


title: How can I display multiple regions shapefiles on an iOS device

text: How can I display multiple shapefiles on an iOS device? My client is asking to display all the states in a different color, and again cities in a different color and all the fields in a different color with lat/longs. I was not able to proceed further with such a huge requirement so how should I start? What are the requirements for moving forward? Are there any GiS specialists needed for me to know completely about the shapefiles? Please give me a suggestion which would be very helpful. Thanks, Rayapudi.


title: How to modify gml using openlayers?

text: i have to edit GML file in openlayers but i am not able to do same. please guide me with some code if possible. i also have to save edited data to hard disk using openlayers Thanks for any help....


title: Value Map Attributes to Multiple Files in QGIS

text: I am hoping someone can share a way that I can some what automatically (none manual) associate a pre-defined CSV file for use with multiple shapefiles as an attribute Value Map. I have a number (in the hundreds) of shapefiles that have the same attributes (Type1 and Type2). Both Type1 and Type2 are numeric codes and I will like to apply text descriptions to them. I have generated a CSV file that has the value and description. What I am wondering is, is there a way to apply this CSV to all my shapefiles automatically, or will I have to open each file and associate the CSV file through the Properties dialog. If there is an alternative method that would provide the same results I would be willing to give it a try. Thanks...


title: How to calculate the distance to multiple locations from each cell of my raster?

text: I'm working with Arcgis 10.1 and Python 2.7.2 How can I calculate the distance to multiple locations from each cell of my raster?   * I tried euclidean distance but it only provides the distance to the **closest** location. I need to know the second closest, third closest and so on until the 150th closest.   * I tried generating a near table but I'm working with high resolution data. The lowest I could go for this dataset are 1,9 million cells. A near table would generare 150 times 1,9 million rows. Looping through these rows takes up a tremendous amount of time. Here is an example (scratch of course) of what I'm trying to accomplish. ![enter image description here](http://i.stack.imgur.com/yzFLZ.jpg) _Red dots are locations, grid represents (49) cells, green lines are to show the distance._ I thought maybe something with an alternated euclidean distance function or an addition of euclidean distances to each single location. I tried combining euclidean rasters, but it didn't work. Can anyone help me or point me in the right direction? A solution in python would be appreciated but I'm glad for any idea , workaround, anything...


title: Openlayers: How to make layer come to top after choosing it in layer switcher?

text:
I have a base layer and 3 image layers . I also have setup a layerswitcher so that user can choose which layer is to be overlayed on the base map and I have 3 sliders such that user can choose the opacity of each layer. The problem as of now is that because of layer indexing, for the user to see the image layer right above the base layer other layers have to be switched off. I want it such that when any of the layers is switched on it automatically comes to the top. This is a **link to the site**. And my code is below:                       var map = new OpenLayers.Map('map');                  var NaturalEarth = new OpenLayers.Layer.WMS(             "Satellite","http://arcims.csl.gov.uk/wmsconnector/com.esri.wms.Esrimap/naturalEarth?",         {layers: 'naturalEarth'}         );             map.addLayer(NaturalEarth);             map.zoomToMaxExtent();                  var imgLayer1 = new OpenLayers.Layer.Image(                 'T2',                 "http://wrf1.geology.um.maine.edu/bipush/OpenLayers-2.13.1/OpenLayers-2.13.1/kml/T2.png",                 new OpenLayers.Bounds(-180, -90, 180, 90),                 new OpenLayers.Size(512, 256)                       //{numZoomLevels: 3}          );                        var imgLayer2 = new OpenLayers.Layer.Image(                 'PWTR',                 "http://wrf1.geology.um.maine.edu/bipush/OpenLayers-2.13.1/OpenLayers-2.13.1/kml/PWTR.png",                 new OpenLayers.Bounds(-180, -90, 180, 90),                 new OpenLayers.Size(512, 256)                      //{numZoomLevels: 3}          );          var imgLayer3 = new OpenLayers.Layer.Image(                 'MSLP',                 "http://wrf1.geology.um.maine.edu/bipush/OpenLayers-2.13.1/OpenLayers-2.13.1/kml/MSLP.png",                 new OpenLayers.Bounds(-180, -90, 180, 90),                 new OpenLayers.Size(512, 256)                      //{numZoomLevels: 3}          );              imgLayer1.setIsBaseLayer(false);         imgLayer2.setIsBaseLayer(false);         imgLayer3.setIsBaseLayer(false);         map.addLayers([imgLayer1,imgLayer2,imgLayer3]);         map.addControl(new OpenLayers.Control.LayerSwitcher());         map.zoomToMaxExtent();                 </script>
             0% <input id = "slide" type = "range"             min ="0" max="100" step = "1" value ="100"                     onchange = "Layer1OpacityUpdate(this.value)" />100%

             0% <input id = "slide" type = "range"             min ="0" max="100" step = "1" value ="100"               onchange = "Layer2OpacityUpdate(this.value)" />100%

             0% <input id = "slide" type = "range"             min ="0" max="100" step = "1" value ="100"               onchange = "Layer3OpacityUpdate(this.value)" />100%
                Thanks in advance