title: How to implement access control for WFS-T updates?

text: I'm implementing a system which has structure like in my question: Design question: what's the role of WFS-T in this kind of webapp? There's spatial data which is updated through WFS-T and non-spatial data which is updated through traditional (Java) web application. When user updates non-spatial data, it is easy to implement authorization rules such as **if current_user = sales_region.manager AND sales_region.status = ACTIVE then all spatial fields can be updated** But can I do the same for the WFS-T geometry updates? **if current_user = sales_region.manager AND sales_region.status = ACTIVE then sales_region.area geometry can be updated through WFS-T requests** Isn't the WFS-T interface just an open XML interface? Anyone who has http access to it can modify, add and remove geometries at will? How could I add business rules checks to protect the geometry updates? I know I could protect the whole WFS-T interface so that only some users can access it altogether, but that's not what I am after. I need to allow or deny some updates based on the data and some defined business rules. Are there any easy solutions for this?


title: GeoServer: Combining multiple raster stores into one wms layer

text: I would like to publish a `DTM` as WMS using GeoServer. For the visual part I have a _hillshaded_ version of the DTM and the original DTM contains the height as floating points (level above sea). Combining these two datasets I would like to create one WMS layer in which:   * `getMap` returns the requested BBOX from the hillshaded dataset.   * `getFeature` returns the height at the specified point from the original DTM. **Is this possible in GeoServer or are there alternatives which could solve this?**


title: Show alternating contour lines

text: Is it possible to show only every other counture line? And how do you do that? (Im new to QGIS Lisboa 1.8.0, but even newer to QGIS Dufour 2.0.1, but have used ArcGIS a bit a while ago, and also some ArcView 3.x earlier. But Im now using QGIS Dufour 2.0.1. for this).


title: Rearranging labels in screen to avoid overlapping

text: I have a map, with some icon over it. Every icon has a label next to it. When icons overlaps, labels are overlapped too. In this case, I want to move labels so they don't overlap anymore, with a line linking them with respective icons. Obviously I know size and position of labels and icon, related to screen coordinates. Which algorithm should I study in order to implement this (I can't use any GIS library)? I've searched for declutter on Internet (I don't speak English very well), but I've found only ideas for hiding/showing icons, not for label moving.


title: Find the co-ordinates of tangent points and intersection of a curved line

text: I have a line with curves in it. Is there an easy way to extract the X Y values for the tangent points and intersection points in ArcMap 10.2.1? Hopefully my simple paint diagram clarifies what I mean. ![enter image description here](http://i.stack.imgur.com/1LWOa.png) Thanks! EDIT: hm, fortunately I only have 1 line with about 8 of these rounded corners. Im going to create a new layer and draw in the straight lines leading up to the curve - snapping to the vertices but leave out the rounded corner. I will then extend the lines by whatever distance makes them all intersect and use that to give me the intersection coords. I can get the tangent ones easily enough. (My curves were originally just a lot of small connected straight lines but I have already manually edited them to arcs and simplified the vertices on straight sections - at most the line has moved 5m from the original which is fine for my uses). Please still suggest solutions =)


title: Very basic drawing on a map

text: I have lat/lng, angle with North and a length of a polygone. I want to draw the polygone on a map. Like in this figure: ![figure polygone](http://i.stack.imgur.com/2UzZo.png) And I want to adjust opacity of it in realtime. Which GIS or Map viewer I gonna use for this purpose. Sadly, I will use it without Internet. Any opensource-related advice or sample application will be best for me. Java is really okay for integrating my existing app. Will openstreetmap or google maps provide me these?


title: How to stack multiple graphics/marks/backgrounds on a single TextSymbolizer in Geoserver?

text: I tried simply adding two marks to a graphic, but only 1 got rendered. I tried adding two graphics to a TextSymbolizer, but only 1 got rendered. I tried adding two TextSymbolizers to the same rule, and it partially worked; they both got rendered, but they got rendered in different positions or stacked in the wrong order. To fix the wrong order, I tried changing their `priority`, but this caused them to never render on the same position. That's probably because of the `space-around` VendorOption, so I set that to -1000. This made it work in about 75% of the cases. I got them both to render, and in the correct order, but sometimes, I get one of the layers rendering in a different spot, or two different TextSymbolizers will overlap because of the negative `space-around`. The Text has varying length, so I can't use fixed sizes etc. My desired result can be seen in the image below; A mark (of type square) with a stroke stacked ontop of a bigger mark (also square). I tried creating a new mark (by implementing MarkFactory etc), but I only want the stroke on the inner mark, and I couldn't find any way to make an edge "unstrokeable". Any way to accomplish this in an easier way? Perhaps a way to get multiple strokes on a mark? ![Desired result](http://i.imgur.com/lEvs79d.png)


title: Geoserver REST api control

text: Recently I am using Geoserver grammatically through its REST api. When I request for a map it gives me a map image. **My request code is like this:**               public string GetMap(string name)                 {                     WebRequest req = WebRequest.Create("http://localhost:8888/geoserver/wms/reflect?layers=" + name);                          using (WebResponse resp = req.GetResponse())                     {                         FileStream imageStream;                         using (imageStream = new FileStream(@"E:\image.png", FileMode.OpenOrCreate))                         {                             var responseStream = resp.GetResponseStream();                             if (responseStream != null) responseStream.CopyTo(imageStream);                         }                     }                     return @"E:\image.png";                 }      It gives me an image like this: ![enter image description here](http://i.stack.imgur.com/4Bh5Y.png) But i need the whole control of Geoserver map like zooming or others: ![the whole map control](http://i.stack.imgur.com/yWNFo.png) How can i do that?


title: Generate rectangular polygons from an existing point dataset

text: I am working with ArcGIS 10.2 I have a point feature class, that looks like this ![enter image description here](http://i.stack.imgur.com/V2puP.png) and I want to generate polygon feature class from it , that is connecting the points like this ![enter image description here](http://i.stack.imgur.com/4XzQX.png) what are the geoprocessing tool that I can use?


title: Trouble importing Feature Classes into Feature Datasets

text: I'm sure there must be a simple reason for this but I can't seem to be able to import a shp file or a feature class from another file or sde geodatabase into a Feature Dataset without getting an error. The feature classes come in fine if I don't don't try to import them into a Feature Dataset. The feature dataset definitely has the same projection as the feature classes I'm trying to import so I don't think that's the problem. The Error I get is ERROR 000210 and just says Failed to convert and cannot create output and Failed to execute (CopyFeatures). Any ideas what I may be doing wrong here. Seems like something so simple. I'm trying to use the Feature Datasets to better organize my data which I think is there purpose. Update: So as you can see below, Jason helped me find the answer. As suggested by Matt below, I added a new Idea to the ArcGIS Ideas site. Here is the link to that, so please promote this up if you agree a better error message is needed here: http://ideas.arcgis.com/ideaView?id=08730000000brxl


title: Show only some buildings

text: How do I show only some buildings? I only need to show 3 cabins in the forrest, but not every house in the whole village around, and the city close by. Maybe some of the symbol buildings in the city area to, but thats less important.   Actually it would be best if the built area could just be gray (with the roads layer on top)(in other words less detailed outside the forrest). The buildings are made up of 2 layers, lines and areas (color inside the building lines). And one of those cabins need to be a different color than the two other, or it can be a symbol (I have an icon for it if needed, and its even preferable (a (DNT) tourist hut)). Im using QGIS Dofour 2.0.1 for this but also have QGIS Lisboa 1.8.0 (wich i used sligtly more), but I'm more used to ArcGIS (and the much older ArcView (3.x)(nope not the new one you might think of :-)).


title: Clip a raster using a shapefile into multiple rasters?

text: I have a raster file representing a large area and a shapefile of 30 polygons representing watersheds. I would like to clip the raster with the shapefile to create 30 rasters representing each watershed? I would like to automate the process as I will be repeating it for another 10 raster files representing large areas. I'm using ArcGIS 10.1 with an advanced arcinfo license. Can this be done without the use of python?


title: DTED Level 3, 4 or 5

text: Does the specification for DTED 0 - 2 also apply to 3 - 5? I know the finest resolution supported by spec MIL-PRF-89020B is 0.1 arc-seconds which I think would barely cover up through DTED level 4, right? I'm guessing then that my code for reading DTED 0 - 2 using just the header information should still work for the later DTED levels but I have no way of proving this. Anyone know where to download sample DTED Levels 3 and up? I could use that as well. Thanks! Matthew Hanna


title: Given latitude/longitude, how to determine whether point is in Japan or China?

text: For OpenStreetMap names transliteration to Latin characters, I need to know whether a point (latitude, longitude) should be transliterated to Japanese or Chinese. What formula `f(latitude, longitude)` would tell whether the point is in Japan or China? It does not need to be perfect, especially since there are islands whose status is disputed.   For transliteration purposes Taiwan is better included in the Chinese-speaking zone. ![China Japan map](http://i.stack.imgur.com/vCfU3.jpg) The goal is actually for the Android app OsmAnd to decide on-the-fly which language to show, so it must be extremely fast, and must be able to work offline.


title: OpenGeo Suite Dashboard not Showing Up on Windows 7 64-bit

text: I have installed OpenGeo Suite 3.1 on my machine (windows 7 64-bit) and I am encountering the following strange behaviour:   When I try to run the Dashboard, nothing appears! HOWEVER, the relative Process at the Task Manager, appears normally and is taking about 30MB of Memory!   Do you have any idea, what causes this strange behaviour?   Thanx in advance,   Petros Apotsos


How to use Google Maps V3 Geocode with Openlayers. Below is my code but the function does not execute when I try to integrate it with Openlayers. This works with my original google v3 file. Below is 2 sets of code one for Google maps with geocode and the other with openlayers                                                              <link rel="stylesheet" href="css/style.css" type="text/css">     <link href="css/style-map.css" rel="stylesheet" type="text/css" />

                 <input type="radio" name="type"                        value="none" id="noneToggle"                        onclick="polygonControl.deactivate()"                        checked="checked" />                 navigate

                 <input type="radio" name="type"                        value="polygon" id="polygonToggle"                        onclick="polygonControl.activate()" />                 draw polygon

     <input id="address" type="textbox" value="Sydney, NSW">          <input type="button" value="Geocode" onclick="codeAddress()">

                               The code that I am battling with is:                        var geocoder = new google.maps.Geocoder();     function codeAddress() {     var address = document.getElementById("address").value;     geocoder.geocode( { 'address': address}, function(results, status) {       if (status == google.maps.GeocoderStatus.OK) {         map.setCenter(results[0].geometry.location);         var marker = new google.maps.Marker({             map: map,              position: results[0].geometry.location         });       } else {         alert("Geocode was not successful for the following reason: " + status);       }     });      }      Below is the working code for the geocoder:                      <!DOCTYPE html>               <meta name="viewport" content="initial-scale=1.0, user-scalable=no"/>     <meta http-equiv="content-type" content="text/html; charset=UTF-8"/>     <title>Google Maps JavaScript API v3 Example: Geocoding Simple</title>     <link href="http://code.google.com/apis/maps/documentation/javascript/examples/default.css" rel="stylesheet" type="text/css" />
     <input id="address" type="textbox" value="Sydney, NSW">          <input type="button" value="Geocode" onclick="codeAddress()">


title: Failure to Delete Vector File from Personal Geodatabase - ArcGIS 9.3

text: I have a personal geodatabase - based on an MS Access .mdb file - from which I cannot delete an old feature class. The tables relating to the feature class have been deleted from the .mdb, but the feature class remains visible as part of the geodatabase in ArcCatalog. I cannot delete it through Arc Toolbox, nor through Catalog, either by right- clicking, or selecting the feature class and choosing the 'delete' button. I receive this error message: > 'Failed to delete selected object(s) The application is not licensed to > perform this operation The table was not found. The Microsoft Jet database > engine cannot find the input table or query 'postcodes_Project'. Make sure > it exists and that its name is spelled correctly.' I understand that the program is struggling as a result of the tables being missing from the database, but is there anything I can do to delete it permanently? Thanks in advance.


title: How do I disable Mouse Wheel zoom functionality in ArcGIS Silverlight?

text: I have an ArcGIS Silverlight app with a Map object. It containes an ArcGISTiledMapServiceLayer and a GraphicsLayer. We've written a custom measure tool for adding in a point at an offset on the graphic layer. I am trying to disable the mouse wheel zoom-in/zoom-out function when this tool is activated to prevent the user from being able to zooom while performing a measurement. I've played around with various combinations of setting the map.IsHitTestVisible = false, graphicLayer.IsHitTestVisible = false, overriding OnMouseWheel(), all with no sucess. It either had no effect, or it disabled ALL mouse events, including the mouse down and mouse move event which I need for the measure tool. I've read thru the documentation on the MouseWheel but I guess I'm not understanding it correctly. I just can't seem to get it to do what I want. And I can't find any specific examples that do what I want either. Any suggestions or places I can look to figure this out? I am using Silverlight 4 and ESRI Silverlight 2.4.


title: How can I optimize pgrouting for speed?

text: I am using pgrouting on a postgis database created through osm2pgrouting. It performs very good on a limited dataset (3.5k ways, all shortest path A* searches < 20 ms). However since I have imported a bigger bounding box (122k ways) from europe.osm the performance went down a lot (a shortest path costs around 900ms). I would think that using A* most of those edges will never be visited as they are out of the way. What I have done so far in an attempt to improve the speed:   * Put an index on the geometry column (no noticeable effect)   * Increased my memory from 8GB to 16GB   * Change the postgresql memory settings (shared_buffers, effective_cache_size) from (128MB, 128MB) to (1GB, 2GB) (no noticeable effect) I have a feeling that most of the work is being done in the C Boost library where the graph is being made so optimizing postgresql will not give me much better results. As I do minor changes to the set of rows I select for A* for every search I am a bit afraid that the boost library cannot cache my graph and has to rebuild all the 122k edges every time (even though it will only use a very limited subset every query). And I have no idea how much is spent doing that compared to the actual shortest path search. Does any of you use pgrouting on a 122k or greater OSM dataset? What performance should I expect? What settings affect the performance most?


title: RadMap for Silverlight reading in Shapefiles

text: I have created shapefiles now which is basically a couple of polygons representing rooms in a building. They appear fine in QGIS, but when I try and load it into the application (RadMap in Silverlight) the screen just goes blank. Anybody have any experience with this? EDIT: I'm using QGIS and have tried MapWindow too and the shapefile CRS is set to the recommended 'EPSG:4326 - WGS 84'. When I drew the polygon I made sure the scale was set to 1:1 in the hope of eliminating the possibility that the shape is just to small to see when it renders in RadMap/Silverlight. Thanks, Jacques


title: Intersecting a raster with a polygon using PostGIS - artefact error

text: I'm using PostGIS2.0 to do some raster/polygon intersections. I'm having difficulty understanding which operation I should use, and what the quickest way of performing this is. My problem is as follows:   * I have a polygon and a raster   * I want to find all the pixels that fall within the polygon, and get the sum of the pixel value   * And (updated problem): I'm getting massive values for some pixels that do not exist in the original raster when I perform the query I'm having difficulty understanding whether I should use `ST_Intersects()` or `ST_Intersection()`. I also don't know what the best approach for summing my pixels is. Here is the first approach that I have tried (#1):               SELECT          r.rast      FROM         raster as r,          polygon as p     WHERE         ST_Intersects(r.rast, p.geom)      This returns a list of `rast` values, which I'm not sure what to do with. I tried calculating the summary statistics using `ST_SummaryStats()` but I'm not certain if this is the weighted-sum of all the pixels that lie within the polygon.               SELECT               (result).count,             (result).sum         FROM (                SELECT                  ST_SummaryStats(r.rast) As result              FROM                 raster As r,                  polygon As p              WHERE                 ST_Intersects(r.rast, p.geom)             ) As tmp      The other approach which I have tried (#2) uses `ST_Intersection()`:                SELECT             (gv).geom,                      (gv).val      FROM       (         SELECT              ST_Intersection(r.rast, p.geom) AS gv         FROM              raster as r,              polygon as p                    WHERE              ST_Intersects(r.rast, p.geom)                ) as foo;      This returns a list of geometries that I analyze further, but I assume this is less efficient. I'm unclear on which is the quickest order of operation also. Should I always choose `raster, polygon` or `polygon, raster`, or convert the polygon into a raster so that it is `raster, raster` ? EDIT: I updated approach #2 with some details from `R.K.`'s link. Using approach #2, I have noticed the following error in the results which is part of the reason why I did not understand the output. Here is the image of my original raster, and an outline of the polygon that is being used to intersect it, overlaid on top: ![enter image description here](http://i.stack.imgur.com/pAlk0.jpg) And here is the result of the intersection using PostGIS: ![enter image description here](http://i.stack.imgur.com/xzXci.jpg) The problem with the result is that there are values of 21474836 being returned, which are not in the original raster. I've no idea why this is occurring. I suspect it is related to small numbers somewhere (dividing by almost 0), but it returns the wrong result.


title: Qgis nightly build - DEM error

text: I installed Qgis 1.8 several months ago. It's been working like a charm. Earlier today, I decided to check out the latest nightly build. Installed GDAL 1.10 (already had 1.9 installed) and SQLite 3.7.16 first, then the nightly build. It's been running fine, however, I just tried to to some DEM raster analysis...which does not work at all. I get an error message - "An error has occured while executing Python code". Any advice?


title: Postgis max contains query

text: I have two tables:               CREATE TABLE periphereies     (       gid serial NOT NULL,       per character varying(60),       geom geometry(MultiPolygon),       CONSTRAINT periphereies_pkey PRIMARY KEY (gid)     )     WITH (       OIDS=FALSE     );          CREATE TABLE archaia_theatra     (       gid serial NOT NULL,       "ΘΕΑΤΡΟ" character varying(254),       "X" numeric,       "Y" numeric,       geom geometry(Point),       CONSTRAINT archaia_theatra_pkey PRIMARY KEY (gid)     )     WITH (       OIDS=FALSE     );      First of all i am trying to write a query that return's me all the periphereies record's that contain at least one archaia_theatra.geom point record.. Then i want to select the periphereies which contains the maximum number of archaia_theatra (compared with all the other perifereies records) I am trying this:               CREATE OR REPLACE VIEW per_max_theatra AS     SELECT count(p.per) FROM periphereies AS p JOIN archaia_theatra AS a ON     ST_Contains(p.geom, a.geom) ORDER BY count(p.per) DESC LIMIT 1      but it's not working Any help? Thx in advance


title: Detection of geographical trend of a variable

text: What is the best method to detect and graphically display the spatial trend of a variable? I'm looking for something that tells me the probability that the trend is from "north-west to south-east" or something like that. Someone wrongly (my opinion) suggested me to use kriging interpolation but I think that it is not correct to infer spatial trend by using estimated variables. I can use ArcGis or R. Thanks


title: MapInfo "shape" file vs. naturalearthdata.com shape file

text: When i load this map shape file (exported from MapInfo as a shape file)               shp.WashingtonLand<-readOGR("C:\\...\\Maps","WLmap2_polyline")      R replies:               OGR data source with driver: ESRI Shapefile      Source: "C:\Users\Andreas\SkyDrive\Speciale\02. Data\02. Statistik", layer: "WLmap2_polyline"     with 2577 features and 23 fields     Feature type: wkbLineString with 3 dimensions     Warning message:     In readOGR("C:\\Users\\Andreas\\SkyDrive\\Speciale\\02. Data\\02. Statistik",  :       Z-dimension discarded      How do I see what this Z-dimension is? (this may be a separate question): Can this explain why the map won't show with points on, or as a layer on top of another shape file like the ones from Naturalearthdata.com?


title: Why are some columns not recognized as numeric even though I have a CSVT?

text: I have problems importing a .csv file. Some Values of the file are not supposed to be formated in text style, but in number style, thus i created an additional .csvt file like explained here. Both files with exactly the same name and the same location are created. The number of columns is the same and I add the .csv file with the "Add Delimited Text Layer" tool. Nevertheless I am not able to choose any columns for a graduated symbology. My .csv file includes 12 rows, of which the 5th, 11th, and 12th column consist numbers. My .csvt file therefore looks like that: _"String","String","String","String","Integer","String","String","String","String","String","Integer","Integer"_ I tried to create it with the standard editor from Windows and Notepad++. Are there any common mistakes done by newbies, like I am one? Thank you for your help!


title: Simulate a click event on FeatureLayer in ESRI JavaScript API

text: When my app changes context (user clicks on a tab) I want to programmatically "Select" the first of 6 points in my Feature Layer so I get the blue box around the point and activate the info window. I looked around the API Docs and found the FeatureLayer's selectFeatures method. This method takes a Query object, and a selection method (I'm creating a new selection so I'll use SELECTION_NEW). The function returns a promise, but do I really need any kind of callback? I'm just trying to select a feature. **Here is the definition of my feature layer:**               var projectsAGOL = new esri.layers.FeatureLayer("http://services.arcgis.com/QVENGdaPbd4LUkLV/arcgis/rest/services/rioLocoSites/FeatureServer/0",{          infoTemplate: infoTemplate,         visible: true,         outFields: ["*"]     });      **Click handler for the Tab:**               jQuery("#prjMap").click(function(){         var query = new esri.tasks.Query(); // Create a Query Object         query.where = "OBJECTID = 0"; // Select the first point from the Feature Layer         projectsAGOL.selectFeatures(query, esri.layers.FeatureLayer.SELECTION_NEW); // Apply selection     });      This code does not produce the desired outcome nor does it produce an error. Nothing happens. Am I even using the correct method to make a selection? I'm more familiar with the Leafletjs mapping library. **Update:** Provided jsFiddle


title: How to bring polygons closer to a point?

text: I have some scattered polygons with a lot of space between them.   I want to make a thematic map but all that white space doesn't look very good, so I want to bring the polygons closer to the center of the area.   Any idea on how to do that in a automatic way? It could be a raw algorithm, postgis query or python function.   Here is an example: Original:   ![enter image description here](http://i.stack.imgur.com/M51e2.jpg)   Final:   ![enter image description here](http://i.stack.imgur.com/cNzSM.jpg)


title: ModifiyFeature Control with an array of vectors

text: I want to put all the drawn vectors in an array and make editable only those which id match with a shapeId var. Do I have to add a new control in the for loop or can I setOptions of the feature to editable?               function initControls(){      colorPoly = document.getElementById('color_picker_'+shapeId).value;           var shapeStyle = {         strokeOpacity: 1,         strokeWidth: 1,         fillColor: colorPoly,         fillOpacity: 0.45,      }           var shapeStyle2 = {         strokeOpacity: 1,         strokeWidth: 1,         fillColor: colorPoly,         fillOpacity: 0.45,         pointRadius: 5,      }           var styleMap = new OpenLayers.StyleMap({         "default": shapeStyle,         "temporary": shapeStyle,         "select": shapeStyle,         "vertex": shapeStyle2      }, {extendDefault: false});           var newShape = new OpenLayers.Layer.Vector( "Editable",{             styleMap: styleMap         });                newShape.events.register("sketchcomplete", newShape, function(){             map.addLayer(newShape);             newShape.id = shapeId;             createdShapes.push(newShape);             alert(createdShapes.length)             toggle = 1;         });           editingToolBar = new OpenLayers.Control.Panel({             displayClass: 'customEditingToolbar',             allowDepress: true         });           var draw = new OpenLayers.Control.DrawFeature(         newShape, OpenLayers.Handler.Polygon,         {             title: "Draw Feature",             displayClass: "olControlDrawFeaturePolygon",             multi: true         });      for (var i = 0; i < createdShapes.length; i++){             if (createdShapes[i].id == shapeId){                 var edit = new OpenLayers.Control.ModifyFeature(createdShapes[i],                 {                     title: "Modify Feature",                     displayClass: "olControlModifyFeature",                     vertexRenderIntent: "vertex"                      });             }         }      editingToolBar.addControls([edit, draw]);      map.addControl(editingToolBar);     }      What is the right way to do it? Thank you very much for your help.


title: Compass error (can't find proper directories)

text: > Compass::SpriteException in Sessions#new > > Showing /home/cartodb/cartodb/app/views/layouts/front_layout.html.erb where > line #12 raised: > > No files were found in the load path matching "icons/*.png". Your current > load paths are: /home/cartodb/cartodb/app/assets/images (in > /home/cartodb/cartodb/app/assets/stylesheets/common/header.css.scss) I've googled and done this: "git submodule foreach git checkout master", and then "rvm use 1.9.2@cartodb --create && bundle install"…. I can get to the next error (and the next) by manually adding various png files but this is not ideal. Any suggestions?


title: Script to "extract by mask" a list of raster files

text: I'm trying to mask a group of raster files and save the new rasters in a different folder using the same name of the original files. It doesn't sound as a difficult task, but I'm having trouble saving-naming each masked raster inside the "for" loop. Here is what I have so far:               import arcpy     from arcpy import env     from arcpy.sa import *          env.workspace = "C:/rasters/threshold"     mask = "C:/GIS/mesoamerica.tif"          rasterlist = arcpy.ListDatasets("*", "Raster")         for i in rasterlist:             outExtractByMask = ExtractByMask(i, mask)             outExtractByMask.save("C:/SIG/MelasCA_30runs_avg/threshold/mesoamerica/" i)      I need to find a way to use the "i" variable in the outExtractByMask.save. Any idea? Thanks!


title: What is a SQL Query to all values in a column with Select by Attribute in ArcGIS?

text: I would like to select all values for an attribute using select layer by attribute to use in model builder and hopefully save the expression for automation purposes. If I have an attribute named COUNTY how do I select all values for this attribute? The end goal is to automate specific layer creation; while avoiding make feature layer and using the field info dialog box. Select layer by attribute(would always be all values in an attribute) Make Feature Layer from selected attributes Apply symbol Output database ![enter image description here](http://i.stack.imgur.com/sqUv4.jpg) I would like to use a SQL expression to select all rows in the county field. *Edit This is a work around to avoid having an array of check boxes in the field info from the make feature layer function. So I would like to select all values in the COUNTY field and make a feature layer from county. So my end attribute table would appear like: ![enter image description here](http://i.stack.imgur.com/SCDsQ.jpg)


title: How to turn this into a routable network that works with pgRouting?

text: I'm only asking again since I didn't get any answer and I'm still stuck here... I updated my last question with more info, but it wasn't answered. When I draw a line with snapping enabled and try to run driving_distance etc the network is broken. (When I zoom in enough I find broken lines here and there after importing it with DBManager. Also this question is new since I'm using the plugin Points to Path. Just trying to solve this for a school project.. Old question:   I have exported all buss stops on a route as points from google earth as a .kml file and converted it into a shape file. (GREEN on the attached screenshot) Using the plugin 'Points to Path' I've been able to draw a line between them, which I think is good. (No manual work so far) (RED on the attached screenshot) Now I need to divide that line into segments/nodes or what ever you call them, to be able to specify the cost for each buss journey (I have this in a excel file and think I can integrate them into the attribute table somehow) My problem is turning this into a routable network with costs for each part of the journey for pgRouting+QGIS+PostGIS+Postgresql. How would you do it? Marry x-mas to all gis.stack-exchangers! ![enter image description here](http://i.stack.imgur.com/hvNFN.png)


title: changing symbology

text: I have problem with my script. Layers are in one geodatabase, and have names like: orne1942, orne 2003, orne 2009 (the same forpart e.g. orne). Script should changing symbology layers wich have the same forpart (uniform it) in many map documents (.mxd). Mayby someone have idea what I am doing wrong ? (I am just begning my python adventure.... ) my script:               import arcpy     path = (r"C:\program2_D.Nowak\map\\")     arcpy.env.workspace=path     arcpy.env.overwriteOutput=True     orneWZOR=arcpy.mapping.Layer(r"C:\program2_D.Nowak\orneWZOR.lyr")     for mapDoc in arcpy.ListFiles("*.mxd"):        print mapDoc        mxd=arcpy.mapping.MapDocument(path + mapDoc)        print mxd        df=arcpy.mapping.ListDataFrames(mxd, "*")[0]        for df in arcpy.mapping.ListDataFrames(mxd):             orne=arcpy.mapping.ListLayers(mxd, "orne*", df)[0]             arcpy.mapping.UpdateLayer(df, orne, orneWZOR)        mxd.save()        print "done"     del mxd     del lasyWZOR, drogiWZOR, orneWZOR, odlogowaneWZOR, zabudowaneWZOR, laki_pastwiskaWZOR <


title: SQL server 08 Limitations

text: The Soil SURvey GeOgraphic (SSURGO) database is a set of very large and complicated shapefiles. We want to load all of the shapefiles for the US into a national SQL database. Are there import limitations for SQL 2008 for Spatial? Also, are there any websites that use a SQL database for interactive thematic mapping? Thanks Bernie


title: Understanding Geodatabase Archiving and H tables in ArcSDE?

text: We have ArcSDE 10.0 on MSSQL 2008 with Geodatabase archiving enabled. In the process of updating a topology I deleted the topology, disabled archiving, unregistered as versioned and then recreated the topology. It seems that when I did this a copy of the feature classes in the feature dataset were created at the root of the geodatabase suffixed with an "_H". When I deleted and created the topology again the features classes were copied to the root of the gdb again with an "_H1" suffixed. Does anyone know what's going on? I was told these were the Archiving tables. Please provide links if available.


title: In leaflet does the addControl automatically handle events via the leaflet API?

text: I need to add additional overlay layers to be controlled by the controls within the leaflet map boundaries, but the code utilizes a leaflet-custom.js file. I am assuming it handles the event and triggers creating DOM elements for the overlay but I cannot verify as the custom file is a minified file so I don't know what configurations were made. When I add my wms tile layers to the controls the wms tile layers do not display as a result of 403 errors with the Forbidden message, the urls are all modified, and the type changes from images/pngs/jpegs to text/html. Therefore is there any way to programmatically add geolocation with with the events?


title: how to convert lat and longitude into distance near poles?

text: I am trying to convert latitude and longitude into meters in polar regions. I am looking for higher accuracy upto 0.00 meters preferably. I tried the haversine formula but it doesnt seem to be accurate enough for small increments in latitude. does anyone know a better formula that can be used near poles?


title: Can't preview a shapefile in GeoServer 2.2.4

text: I have problem with previewing layers in GeoServer. I installed Geoserver 2.2.4 (following instructions in manual) and everything was going well. I managed to login to it and create my own new workspace. I've created a new store and published that layer but when I try to preview that layer with openLayers all I see is empty map with no features on it. I've also tried with example on http://docs.geoserver.org/2.2.4/user/gettingstarted/shapefile- quickstart/index.html to see if the problem is with my shapefile but the same thing happens. I really don't know why because all I did was the same as it is written in manual. I hope somebody can help me with that. Thanks


title: Can "Getting to know ArcGIS Desktop" for version 9 be used as a tutorial for version 10?

text: I am a GIS n00b. Our company has ArcGIS 10 desktop, and since it is multi-user license, I was allowed to install on my computer. I was given a copy of "Getting to know ArcGIS desktop" but it is for ArcGIS 9. Are the two versions so different that I should just get myself the latest copy for 10? ESRI sells it for about 70 dollars so I'd like to not spend the money, but I'd also like to learn the system as well as possible. Thanks!


title: How to get the Nth LineString from a MultiLineString in PostGIS?

text: I need a function that gets the Nth LineString from a MultiLineString in PostGIS, somthing like ST_PointN that gets the Nth point from a LineString.


title: mixing lon/lat and mercator layers in a MapServer .map file

text: I have a basic question about my MapServer .map file, which has two layers, each with a different coordinate system, longlat vs. mercator. See my .map file below. If I leave it as shown below, with the "version 2" UNITS/EXTENT active, it correctly draws my "MrSid" in the right place on the earth, but doesn't draw my "shape layer" layer. On the other hand, if I comment out the "version 2" UNITS/EXTENT and activate the "version 1" UNITS/EXTENT, it then draws "shape layer" but not "MrSid". So it's acting like the the UNITS/EXTENT has to match those of the layer for the layer to be drawn. Is that right? If so, how do I mix layers with different UNITS? (By the way, I also have a secondary question about this .map file. If I don't include any PROJECTION entry at the outer/MAP level, it will draw the "shape layer" layer but not the "MrSid" layer. But if I add a PROJECTION, regardless of whether the projection is `"proj=longlat" "ellps=WGS84" "datum=WGS84"` or `"proj=merc" "lon_0=0" "k=1" "x_0=0" "y_0=0" "a=6378137" "b=6378137" "towgs84=0,0,0,0,0,0,0" "units=m"`, it will draw either layer -- but not both at the same time, as per the main issue above. So my secondary question is why it seems to be happy using a default projection with my longlat shape file, but not with my mercator MrSid file.)               MAP       IMAGETYPE   PNG            # version 1       # UNITS dd       # EXTENT      -180 -90 180 90            # version 2       UNITS meters       EXTENT -20037508.34 -20037508.34 20037508.34 20037508.34            SIZE        550 450       INTERLACE OFF       PROJECTION         # "proj=longlat" "ellps=WGS84" "datum=WGS84"         "proj=merc" "lon_0=0" "k=1" "x_0=0" "y_0=0" "a=6378137" "b=6378137" "towgs84=0,0,0,0,0,0,0" "units=m"       END       SHAPEPATH   "../data"            OUTPUTFORMAT         NAME png         DRIVER "GD/PNG"         MIMETYPE "image/png"         IMAGEMODE RGB         FORMATOPTION  INTERLACE=OFF         TRANSPARENT OFF       END            # Start of LAYER DEFINITIONS-------------------------------       LAYER         NAME "shape layer"         DATA dourokukan         STATUS DEFAULT         TYPE LINE              CLASS           NAME "“shape class"           STYLE             COLOR 227 227 127           END         END       END        LAYER         NAME "MrSid"         DATA "c:\mrsid\adams\ortho_1-1_1n_s_nd001_2010_1.sid"         STATUS DEFAULT         TRANSPARENCY 100         TYPE RASTER         PROJECTION           "proj=utm" "ellps=GRS80" "datum=NAD83" "zone=13" "units=m"         END         LABELITEM ""         TEMPLATE VOID       END       # End of LAYER DEFINITIONS -------------------------------     END


title: ArcMap query layer with Oracle Locator displaying no results

text: I have an Oracle 11g database setup with Oracle Locator. I have created a table with the this query:               CREATE TABLE geotest (id NUMBER(38) NOT NULL, shape SDO_GEOMETRY);      Then I inserted one point with the this query:               INSERT INTO geotest VALUES(1,SDO_GEOMETRY(2001,3395,SDO_POINT_TYPE(-111.949439,40.722283,NULL),NULL,NULL));      The result is no data in the query layer in ArcMap. It allows me to create the layer, and validates. But then it makes me go through the advanced properties which I suspect means it can't see any of the rows to automatically pull the SRID. Any guidance would be appreciated! **Update** : So when I insert the metadata into the table using this:               insert into user_sdo_geom_metadata values ('geotest','shape',sdo_dim_array(sdo_dim_element('X',-180,180,0.005),sdo_dim_element('Y',-180,180,0.005)),4326);      And creating the spatial index, I can see the following point inserted into the database using this query:               insert into geotest values(1,'test',sdo_geometry(2001,4326,sdo_point_type(0,0,0),null,null));      So it would appear that a point with all 0s work, but anything beside that doesn't. I haven't found anything in the documentation that relates to this, or even any users complaining about that. Any either tricks that may be useful? **Update 2** : Alright it appears that after I drop the index, and create it again, the points then show up. I am not an Oracle expert, so this confuses me, anyone have an explanation for this behavior? I thought once an index was created, as soon as a point is inserted, the index updates.


title: User defined Coordinate Reference System (CRS) in QGIS is not located where expected?

text: I have received a shapefile with a user defined projection (i.e., +proj=krovak +lat_0=-32.923433 +lon_0=116.47485 +alpha=111.423433 +k=1.00004 +x_0=44328.529781 +y_0=1345943.211845 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs). The site is located in Australia but after automatic re-projection into GDA94 (the default CRS) it places the site in northern China. I have viewed this particular shape file on the client's GIS system (ArcGIS) and it seems to locate fine. Any suggestions as to what I need to do or am I missing something? I have attempted to manually change the CRS (did not work) and I also saved the layer to a new SHP file and manually chose the CRS but this also did not work.


title: Create filled contour map with GRASS

text: ## Background I'm trying to make contour/isoline/isochrone maps to show what areas are within walking distance to schools for my pet project to help with the Safe Routes to School program. I followed the pgRouting workshop and underdark's posts (all linked to from my project page) to import OpenStreetMap into a PostGIS database, then used the driving_distance() function to create a point layer that had a cost attribute, then ran the contour plugin of Quantum GIS to create contour lines, which I then styled in QGIS to create this image (pardon the spurious lines, not sure where they came from): ![enter image description here](http://i.stack.imgur.com/6uXTU.jpg) ## What I'd like to do **I'd like to replace the QGIS contour plugin and styling feature with GRASS** , as I'd like to create this image dynamically for a web service, and I don't think QGIS is meant for that (at least performance-wise, I think it may be possible with Python). So I'm trying to figure out how to do this same thing with GRASS, which appears to be used quite a bit in web services. From my reading, it appears I'll want to use something like `v.surf.rst`, though that seems like it would generate more of a heat map with a gradient, rather than the type of image shown above with solid colors for each band (in this case each band is 1/4 mile).


title: ArcGIS 10 creating a shape file from current dataframe extents in layout view

text: Where is the tool for this? Have looked around and the closet thing i can find is Toolbox's Grid/Strip Map Index tools under Data Drive Pages. I just want to be able to create a **single** polygon rectangle shp file based on the data frame (in layout view) for any given scale / page setup.


title: Is it possible to realign polylines from KMZ in Google Earth?

text: I have a street network that has been extracted from AutoCAD DWG's and converted to a KMZ file. The resulting polylines, when launched in Google Earth, do not visually align with the imagery. I am not sure what the source data is, and not really concerned with anything other than realignment (transformation?). I have ArcExplorer but no other GIS at the moment other than AutoCAD and Microstation, which aren't true GIS platforms. Apologies for the rambling, post, thanks.


title: How to echo gdal process to screen and pipe to text file in Windows?

text: gdal_warp and similar give a nice progress summary when run and I want this output to show in the osgeo shell and also pipe to a text file. The issue is when I add the pipe the output doesn't show on the screen. How can I do both?               gdal_translate -expand rgba -of GTiff -a_srs EPSG:3857 -co TILED=YES -co BIGTIFF=YES %1 %2\%3\%3_translated.tif > %2\%3\%3_translate_log.txt      puts the output to a text file but then doesn't show the progress on screen. I tried               Echo gdal_translate -expand rgba -of GTiff -a_srs EPSG:3857 -co TILED=YES -co BIGTIFF=YES %1 %2\%3\%3_translated.tif > %2\%3\%3_translate_log.txt      and               Echo.gdal_translate -expand rgba -of GTiff -a_srs EPSG:3857 -co TILED=YES -co BIGTIFF=YES %1 %2\%3\%3_translated.tif > %2\%3\%3_translate_log.txt      but neither works.


title: Auto-update the map service & REST in ArcGIS server 10.1

text: I have one map service which is updating per day so I need to manually overwrite the service in ArcGIS 10.1 In previous versions 10.0 its very easy..just update the .mxd then REST will be auto updated.. Is there any way to auto-update the map service & REST in ArcGIS Server 10.1 I tried the ADMIN API but its not proper way..also there is no settings for service definition file in 10.1.. Any help/idea/workflow will be highly appreciated.. thanks in advance : ) **UPDATE:** The problem is resolved by adding arcgis account.


title: Can markers changes undo through history/revert to original?

text: I'm a contributor to a food foraging map and found that many markers are in the ocean. I then also noticed that I had an undo at the top of gmaps and realised I'd moved 5 trees trying to get to my area. Is it possible to see marker/map history at maps.google.com? Is it possible to revert markers to original position? I'm completely open to programming this using the api if necessary - surely this information is somewhere...


title: What projection to use for Google Maps?

text: I'm loading a WMS layer which comes in WGS84 into Google Maps. From what I've read so far, I should reproject this layer to Google Pseudo Mercator (either 3857/900913). However, when I save the layer in either one of those two projections, it is not displayed in Google Maps, but I can view it in Data Preview of Geoserver. Others with WGS84 are displayed correctly but with a position shift. Do I really need to reproject the layer? I'm not sure whether it is not just a matter of some inaccuracy.


title: How to debug a QGis python plugin

text: I'm trying to write a `python` plugin for `QGis`. I'd like to debug some DB queries and I need a way to display the result. I thought about using a popup to display it, but I don't know how to show a **popup** from inside a `python` plugin. Can anyone suggest me how to do that?


title: Uninstalling ArcGIS10 License Manager?

text: I recently uninstalled the ArcGIS 10 from my computer, the software got uninstalled , but I am unable to uninstall the License manager. It is giving me an error "you don't have sufficient access to uninstall ArcGIS license manager 10"..and I am logged in as Admin. Due to this problem I am not able to install the ArcGIS 10.1.


title: Get layers and apply renderer to layer list from SDE ArcMap

text: I have the code that gets the layers from SDE, but when i add them to the map they have a different symbology. there is about 80 layers and i am accessing the SDE by getting the workspace from a selected layer in the TOC. The symbology renderer layer files are located in a directory. So all i want to do is loop through the directory and get the renderer for each layer in my list. But my code only gets the renderer for 1 file not for all? How can i compare feature class names so it van get the renderer for each layer by using the feature class name?? I appreciate your time and any ideas.               Public Sub AddSDECoverageLayer()               Dim pMxDoc As IMxDocument     Dim pMap As IMap               Set pMxDoc = ThisDocument                Dim intX As Integer          'Add each FeatureLayer to the map document     'This should assign layer names to the layer it gets from SDE     'But it only gives the first layer in the list a layer name.      'There is about 80 layers so cant do if statement 80 time     'What other way can i do this??           For intX = 0 To ListBox2.ListCount - 1                      Dim xsClass As String            xsClass = ListBox2.List(intX)            Dim xs As String            If xsClass = "PROPERTY_ADMIN.topo_ownership" Then            xs = "MM Ownership"            End If            If xs = "PROPERTY_ADMIN.MM_Rights_of_Way" Then            xs = "Rights of way"            End If                     Dim pFLayerXS As IFeatureLayer                Set pFLayerXS = importData(pFWorkspace, xsClass, xs)                pMxDoc.FocusMap.AddLayer pFLayerXS                pMxDoc.ActiveView.Refresh           pMxDoc.UpdateContents                    Next               End Sub      ## ---------------------------------------------------------------------------               Public Function importData(pWorkspace As IWorkspace, FClass As String, fcName As      String) As IFeatureLayer          'Function to create Feature layer from Feature class     Dim pFeatureWorkspace As IFeatureWorkspace     Dim pFClass As IFeatureClass     Dim pUnknown As IUnknown     Dim pMxDoc As IMxDocument     Dim pFeatureLayer As IFeatureLayer     Dim pDataSet As IDataset     Dim pFLayer As IFeatureLayer                     ' Get the selected layer     Set pMxDoc = ThisDocument     Set pUnknown = pMxDoc.SelectedLayer 'pMxDoc.FocusMap.Layer(0)      If Not TypeOf pUnknown Is IFeatureLayer Then      Exit Function     End If     Set pFLayer = pUnknown               ' Make sure the layer is from SDE and get the parent workspace     Set pDataSet = pFLayer.FeatureClass          Set pFWorkspace = pDataSet.workspace     Set pWorkspace = pFWorkspace     If Not pWorkspace.Type = esriRemoteDatabaseWorkspace Then     MsgBox "Not an SDE feature class"     Exit Function     End If     Set pDataSet = pWorkspace               Set pFeatureWorkspace = pWorkspace     Dim intY As Integer     For intY = 0 To ListBox2.ListCount - 1     Set pFClass = pFeatureWorkspace.OpenFeatureClass _       (ListBox2.List(intY))       Next                    Set pFLayer = New FeatureLayer     Set pFLayer.FeatureClass = pFClass          Dim pGFLayer As IGeoFeatureLayer     Set pGFLayer = pFLayer               Dim intz As Integer     For intz = 0 To ListBox2.ListCount - 1          'This gets the renderer but only for one layer,      'How can i make it get the renderer for all the layers in the list.     'There is 80 layers, in a directory. i just need to assign the renderer     'from this directoryto eac item in the list box          If ListBox2.List(intz) = "PROPERTY_ADMIN.topo_ownership" Then       Set pGFLayer.Renderer = GetRenderer("V:\Arc\Layers\MM Ownership.lyr")            End If       If ListBox2.List(intz) = "PROPERTY_ADMIN.MM_Rights_of_Way" Then       Set pGFLayer.Renderer = GetRenderer("V:\Arc\Layers\MM Rights of Way.lyr")     End If       Next                pFLayer.Name = fcName          Set importData = pFLayer          Exit Function          End Function


title: WMSgetFeatureInfo problem

text: I want to get attributes of feature using getFeatureInfo:               if(sss==null){        if (sld==null){          sld= new OpenLayers.Control.WMSGetFeatureInfo({         url: 'http://localhost:8080/geoserver/wms',          title: 'Identify',         layers: myLayer,         queryVisible: true          });        app.mapPanel.map.addControl(sld);        sss=app.mapPanel.map.getControlsBy('title', 'Identify');     }                                        }     for(var i=0;i<sss.length;i++){     alert("BROMACIA");      sss[i].activate();     }     sss[0].events.register("getfeatureinfo", app.mapPanel.map, function showInfo(evt) {     alert(evt.text)});      And nothing heppens.   When i try to do in firebug:               >>> sss[0].events.register("getfeatureinfo", app.mapPanel.map, function showInfo(evt) { alert(evt.text)});      firebug says `undifined`. Whats wrong?


title: Looking for bathymetry data for Lake Sakakawea in North Dakota

text: Navionics has a nice web map that has some incomplete coverage of Lake Sakakawea. Garmin's BlueChart looks to have something for about $160, but appears to be only compatible with GPS units. Don't need the entire lake just a few areas of it. Would like the data to be ArcGIS compatible. I have attempted to do raster to vector in ArcGIS taking screen shots of Navionics web map but haven't had much success. Any ideas / help would be most appreciated.


title: Slow performance on splitting a shapefile by attributes

text: I've tried to develop some workflows with the ArcPy module for splitting a shapefile based on the unique values in a specified field. The code works as expected except one thing I am not able to figure out.               import arcpy, os     from arcpy import env          #getting user inputs     in_fc = arcpy.GetParameterAsText(0) #'roads.shp'     outfolder = arcpy.GetParameterAsText(1) #r'C:\GIS\Temp\OutSHP'     fieldname = arcpy.GetParameterAsText(2) #"FEATURE"     createfolder = arcpy.GetParameterAsText(3) #true/false     env.workspace = outfolder          #user chosen NOT to create a folder with the name of feature class     #     if createfolder == "false":              with arcpy.da.SearchCursor(in_fc, ["{0}".format(fieldname)]) as cursor: #defining a search cursor             for row in cursor:                      if not arcpy.Exists(os.path.join(outfolder,row[0])): #checking if shapefile already exists                     #Analysis toolbox > Extract > Select GP tool to select features                     arcpy.Select_analysis(in_fc, os.path.join(outfolder,row[0]), ''' "FEATURE" = '{0}' '''.format(row[0]))                 else:                         print "The shapefile with this name already exists"               #user chosen to create a folder with the name of feature class     #     else:             with arcpy.da.SearchCursor(in_fc, ["{0}".format(fieldname)]) as cursor: #defining a search cursor             for row in cursor:                 directory = os.path.join(outfolder,row[0]) #create a path for directory that will be created                 if not os.path.exists(directory): #checking if directory exists                     os.makedirs(directory)                          if not arcpy.Exists(os.path.join(directory,row[0])): #checking if shapefile already exists                     #Analysis toolbox > Extract > Select GP tool to select features                         arcpy.Select_analysis(in_fc, os.path.join(directory, row[0]), ''' "FEATURE" = '{0}' '''.format(row[0]))                     else:                         print "The shapefile with this name already exists"      I've created a script tool in ArcGIS which takes an input shapefile, an output folder, a field name and whether user wants to create a folder for each unique value in the specified field before creating a shapefile there (the folder and the shapefile will have identical name except .shp in the end for shapefiles). The issue is that the tool gets the work done for 5 seconds when user chooses to create folder with the unique name for each feature class and for around 3 _minutes_ if the shapefiles are created directly in the output folder. I am having hard times trying to find out anything that could stipulate this behaviour in the code.


title: QGIS field calculator and Postgres permissions

text: I’ve given user x the following               GRANT ALL ON ALL TABLES IN SCHEMA y TO x;     GRANT ALL ON ALL SEQUENCES IN SCHEMA y TO x;     GRANT ALL ON ALL FUNCTIONS IN SCHEMA y TO x;      However in QGIS they dont have access to check/tick the 'update existing field' button. They can chose a field to update and write an expression for its update e.g $area but the OK button remains greyed out. If I as administrator load the layer, the check/tick box appears and I can write and run expressions. I'm sure its a permissions thing in PostGIS. Si ![update exisiting field](http://i.stack.imgur.com/QwRjJ.png)


title: Check if a land parcel (polygon) has access on a road or not

text: I have a feature layer with land parcels as polygons. Polygons (land parcels) are either adjacent in all its around sides by other polygons but some of them have a road a which is not actually a feature layer representing roads but is just an empty space between parcels. In this case, parcels have access to roads. The question is how can I automattically check if a polygon has access onto a road using VBA and AcrObjects. Thanks Deme


title: How to import ECW into my postgis database with raster2pgsql (windows)

text: I'm working on windows 7 64 bits, postgis 2.01. I have to import ECW files into my postgis database. Unfortunatly, my postgis seems to has been compiled without gdal ecw support.When i do a raster2pgsql -G, ECW does not appears. I have tried to follow instructions erected in this following post :   How to add ECW support to GDAL bundled with PostGIS Windows So I paste the raster2pgsql.exe into my C:\Program Files\PostgreSQL\9.1\gdal- data folder. And i paste gdal19.dll from C:\Program Files (x86)\Quantum GIS Lisboa\bin toward C:\Program Files\PostgreSQL\9.1\bin. I rename the dll like libgdal-1.dll. But it does not work and i can't import ecw with raster2pgsql.I have an error message saying : application error.For information/. Moreover, when i do a SELECT postgis_full_version(); "POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" And my qgis 's gdal has been also compiled with gdal 1.9.1. Moreover, for information, my GDAL_DATA (environment system variable) pointed to C:\Program Files\PostgreSQL\9.1\gdal-data Have you got any idea? How to use ECW with postgis? I have to use this ECW format. My client give to me data with this format and i did not convert raw data in ECW with a proprietary software.. Thank you!


title: How to print a WMS map with chosen scale and WMS zoom level?

text: I am trying to print a WMS map in QGIS. If, while in the print composer, I chose a "25000" scale, the composer chooses a WMS layer that actually shows a 1:50000 map. To get down deep at the 1:25000 WMS map in the print composer preview, I have to trick it by telling the composer to use a smaller scale, but then I don't get the scale I want. How does QGIS chose its WMS zoom level? I need to chose a scale AND a specific WMS zoom level. I tried the "Lock layers for map item" check box but it is not what I hoped for. Is there a way to accomplish this or is this feature missing? Any other open source alternative to get the same result would be fine if QGIS is not able to do that. Thank you very much.


title: How do you store brief cartographical tips, information, possible future designs that you find?

text: I find inspiring maps and/or brief, 1-line design tips (i.e. "when satellite image is used as a basemap, I like to desaturate the tarnation out of it" that I may want to consider later using in a project or as inspiration as a design for a future map, from a variety of sources (blogs, tweets, forums, here) during my daily internet browsing and online interaction. I'm curious to hear how each of you, especially cartographers and others specializing in map design store all of this information ? Right now, I store these bits of information across the sources that I find it on: favorited tweets, bookmarks in my browser (do you use tags ?), pinboard, favorited gis.se questions, and a text file (with the tidbit of information and the URL of where I found it). I occasionally notice that I ended up bookmarking it in a different place, or not knowing that I had that information already saved somewhere until after I needed it. Where do you find it efficient to store these tips ? Do you only look at these saved tips when you're looking for something specific or is it beneficial to occasionally review them without a specific goal in mind ?


title: QGIS 2.2 Crash dumped on Windows 8.1

text: After I installed QGIS 2.2 via OSGeo4W, as I click the icon of QGIS, then it shows this up, ![enter image description here](http://i.stack.imgur.com/hijmL.png) so I can't open QGIS 2.2 Here's some look like troublesome lines in my installation log,               2014/03/23 13:42:15 io_stream_cygfile: fopen(/etc/setup/setup.rc) failed 2 No such     file or directory     2014/03/23 13:42:15 Current Directory: C:\Users\Heinz\Installation files     2014/03/23 13:42:15 Root Directory: C:\OSGeo4W64     2014/03/23 13:42:15 Could not open service McShield for query, start and stop. McAfee may not be installed, or we don't have access.     2014/03/23 13:42:17 root: C:\OSGeo4W64     2014/03/23 13:42:17 Selected local directory: C:\Users\Heinz\Installation files     2014/03/23 13:42:17 Menu name: OSGeo4W     2014/03/23 13:42:17 net: Direct     2014/03/23 13:42:19 site: http://download.osgeo.org/osgeo4w/     2014/03/23 13:42:19 mbox note: Unable to get setup.ini from <http://download.osgeo.org/osgeo4w/>     2014/03/23 13:43:00 site: http://download.osgeo.org/osgeo4w/     2014/03/23 13:43:00 io_stream_cygfile: fopen(/etc/setup/installed.db) failed 2 No such file or directory     2014/03/23 13:43:01 io_stream_cygfile: fopen(/etc/setup/timestamp) failed 2 No such file or directory      How to fix this problem？ I am working under Windows 8.1 64-bit * * * **Update#1** I follow the suggestion of the answer in the post, How to debug QGIS installation problems on XP?, download Microsoft Visual C++ 2010 Redistributable and then install qgis 2.2 using OSGeo4W package, and here's the same part of the installation log,               2014/03/23 16:17:39 io_stream_cygfile: fopen(/etc/setup/setup.rc) failed 2 No such file or directory     2014/03/23 16:17:39 Current Directory: C:\Users\Heinz\Installation files     2014/03/23 16:17:39 Root Directory: C:\OSGeo4W64     2014/03/23 16:17:39 Could not open service McShield for query, start and stop. McAfee may not be installed, or we don't have access.     2014/03/23 16:17:40 root: C:\OSGeo4W64     2014/03/23 16:17:40 Selected local directory: C:\Users\Heinz\Installation files     2014/03/23 16:17:40 Menu name: OSGeo4W     2014/03/23 16:17:40 net: Direct     2014/03/23 16:17:43 site: http://download.osgeo.org/osgeo4w/     2014/03/23 16:17:43 io_stream_cygfile: fopen(/etc/setup/installed.db) failed 2 No such file or directory     2014/03/23 16:17:43 io_stream_cygfile: fopen(/etc/setup/timestamp) failed 2 No such file or directory      I can successfully launch Qgis browser, but I still have the same dumped message that I can't launch qgis desktop. * * * **Update#2** I have uninstalled the 64-bit version and installed the 32-bit OSGeo4W package. I have checked the installation log, just like the one I post in **Update#1** , and the rest part of it is fine. I still got the same error message, but unlike the 64-bit version message, the same mbox would appear after I click the confirm button, this time I just have to click **confirm** in the mbox one time. And I can still launch qgis browser. **Are there someone successfully launching qgis 2.2 on windows 8.1 x64 computer？**


title: How to create multi-ring buffers automatically in QGIS?

text: I'm wondering if anyone knows of a tool or plugin for QGIS that allows for multi-ring buffers. I've found a post about doing multiple buffers one by one, but was hoping there was a way to do this all in one place. I know in ArcMap that you can set the number of buffers you want and distance, so just hoping there's something out there in QGIS to accomplish this.


title: How do I change the ArcGIS Server Web Adaptor name?

text: If I decide I want to change the Web Adaptor name, is it as simple as going to the Web Adapter tab in the Server Manager and clicking unregister? And then would I just rerun the Web Adaptor install and choose a different name? Or what is the proper way of changing Web Adaptor names? Thanks.


title: Draw polygon shape via Qbutton

text: Ideally; a Ui Qbutton is needed that will activate map canvas editing to add a polygon via the classic mouse click rubber band method. Upon a concluding right click with the active tool, editing would be closed and the feature would be added to a designated layer. Novice python skills amplify my want to just 'drag and drop' the "Add Feature" resource from the tool bar into the developing plugin Ui (unfortunately not so easy). I'm wondering if there is a (relatively) simple means of connecting to existing code via imports for enabling use into a new plugin? Thanks for considering. [I hope to implement this tool connection and user QLineEdits from the Ui into " _create polygon reference, add feature polygon_ " sections of code in the following posting: Blank Shapefile Output : Solved!]


title: Styling based on hex color codes within table?

text: I have a table of polygons named 'testtable'. Within it is a column 'color', which has hex codes. I want to style the polygons based on the hexcode values.. is this possible? I was trying something like this, but it's obviously not working.. # testtable{ polygon-fill: {color}; } Thx,


title: Errors while running Region Group tool

text: I am getting an error while using Region Group on a large dataset. The errors are specified below: > ERROR 999999: Error executing function.   >  The table was not found. [VAT_final002]   >  A column was specified that does not exist.   >  A column was specified that does not exist.   >  ERROR 010005: Unable to allocate memory.   >  ERROR 010067: Error in executing grid expression.   >  Failed to execute (RegionGroup).   > I am currently running ArcGIS 10.1. I have 4 GB of total RAM. The raster I'm processing has 154242 columns and 97288 rows with 30x30 cell size (approximately 500 MB file size). All values = 1 or "nodata". Does anyone have a solution to this other than having to process in small sections?


title: Changing symbology/color in kml when exported from arcGIS/QGIS

text: I have marked some points and used symbology in ArcGIS/QGIS... when I export the file to kml, the color/symbology changes and only single color/symbol remains. Is there any possible way that the symbol/color that I use in Arc or qgis may retain in kml as well... NO KNOWLEDGE OF PROGRAMMING... :( any tool or software (easy) will be a great help... regards,


title: Spatial Join Python Command to summarize attributes

text: I am trying to run the Spatial Join function through the python shell in ArcGIS 10. When you do the point-and-click method of spatial join, you are given a tool window with an option that asks "How do you want the attributes to be summarized?" However, the python syntax doesn't show a parameter for this option. How do I control the attribute summarization through python? Why is there no parameter for this? Any suggestions are much appreciated, thanks, mike


title: How to copy a feature class within the same GDB

text: I'm trying to copy/clone a feature class within the same file gdb using C#. I've figured out that this won't work with IDataset or IFeatureDataset. So far I only found IGeoDBDataTransfer.Transfer but this is for copying between different workspaces. Copying and pasting geodatabase datasets ArcObjects 10.1 I know I could do the whole things with a featureBuffer, but that's a very long winded way. And I don't want to use geoprocessing like CopyFeatures. Yeah, I know I'm fussy ;-) If any of you got an idea how to get round this I would very much appreciate it! Many Thanks in advance!


title: Postgis 2.1 error after update

text: After a usual update progress - "apt-get upgrade" - I can not connect with qGIS to my PostGIS DB anymore. Error in the qGIS logfile:               Your database has no working PostGIS support.      Erroneous query: SELECT postgis_version() returned 7 [ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": liblwgeom-2.1.2.so: cannot open shared object file: No such file or directory]     Retrieval of postgis version failed      After asking for my version using `SELECT PostGIS_full_version();` This occured:               ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": liblwgeom-2.1.2.so: cannot open shared object file: No such file or directory      CONTEXT:  SQL statement "SELECT postgis_lib_version()"      PL/pgSQL function postgis_full_version() line 19 at SQL statement      How can I fix that? Thats my environment: xUbuntu 14.04, PostgreSQL 9.3, Postgis 2.1, QGIS 2.3,


title: Does Shooting star 1.05 not consider reverse_cost?

text: > **Possible Duplicate:**   >  PGRouting with 1.05 library in Postgres 9.1 Shooting star pgrouting lib 1.05 with postgres 9.1. Dijkstra and A-star working fine but "SHOOTING STAR" not considering "reverse_cost". Also if i try long route then it gives server disconnected error. Please help


title: Problem with showing Shape file in Geoserver

text: I have a number of Shape files from a customer. I can read them with both ogrinfo and QGIS. They don't have a prj-file, but I know that the projection is EPSG:3006. When I add them as layers in Geoserver, I declare the projection to be EPSG:3006 in the box "Declared SRS" by inputing "EPSG:3006". I set "SRS Handling" to be "Force declared". (It makes no difference if I add a .prj-file with a projection definition.) I can then calculate the bounds successfully. When I try to preview the layer in Openlayers, I get an empty map. When I try to load the layer as WFS, I get an empty result. There are no error messages in the logs as far as I can see. I have tried to change the log level. If I run ogr2ogr and save them as new shape files, they work in Geoserver. If I export them to a new shape file from QGIS, they work in Geoserver. Is there a way I can use ogr or QGIS to show my customer what the problem is?


title: Distance between DB field and Point using PostGIS

text: I'm using Postgresql 9.0 with PostGIS 2.0 extension. I've got the table and geographic field `location`:               =# SELECT * FROM geography_columns;      f_table_catalog | f_table_schema |   f_table_name   | f_geography_column | coord_dimension | srid | type       -----------------+----------------+------------------+--------------------+-----------------+------+-------               _dev   | public         |          objects | location           |               2 | 4326 | Point               =# \d objects                                            Table "public.objects"          Column     |            Type             |                           Modifiers                                ----------------+-----------------------------+---------------------------------------------------------------      ...      location       | geography(Point,4326)       |       And have got records with coordinates:               =# SELECT ST_X(location::geometry), ST_Y(location::geometry) FROM objects;       st_x  | st_y       --------+-------             |            -73.01 | 15.01      -73.02 | 15.02      -73.03 | 15.03     (4 rows)      I want to grab records within distance in meters from Point(-70 15). I've calculated this distances to each point (just for know it):               =# SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-70 15)'), ST_GeographyFromText('SRID=4326;POINT(-70.01 15.01)'));        st_distance         ------------------      1543.03968305103     (1 row)          =# SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-70 15)'), ST_GeographyFromText('SRID=4326;POINT(-70.02 15.02)'));        st_distance         ------------------      3086.04520065436     (1 row)          =# SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-70 15)'), ST_GeographyFromText('SRID=4326;POINT(-70.03 15.03)'));       st_distance        ----------------      4629.016529543      Now I see the closest point in 1543 meters and farthest point in 4629 meters. I've used very similar query to example from here: > Query and measurement functions use units of meters. So distance parameters > should be expressed in meters, and return values should be expected in > meters (or square meters for areas). >   >   >     -- Show a distance query and note, London is outside the 1000km > tolerance >       SELECT name FROM global_points WHERE ST_DWithin(location, > ST_GeographyFromText('SRID=4326;POINT(-110 29)'), 1000000); >                 =# SELECT location FROM objects WHERE ST_DWithin(location, ST_GeographyFromText('SRID=4326;POINT(-70 15)'), 5000);      location      ----------     (0 rows)      and I didn't get locations within 5000 meters. More than that I cannot calculate distances in meters from location to Point:               =# SELECT ST_Distance(location, ST_GeographyFromText('SRID=4326;POINT(-70 15)')) FROM objects;        st_distance         ------------------           323718.841976209      324792.390871594      325869.606072403     (4 rows)      There numbers are a lot bigger than numbers from distance in Point to Point example. I'm guessing something wrong with my column because I cannot even grab points with same coordinates:               =# SELECT location FROM objects WHERE ST_DWithin(location, ST_GeographyFromText('SRID=4326;POINT(-70.01 15.01)'), 5000);      location      ----------     (0 rows)


title: Multipoint routing logistics software solution?

text: I work for a state health program and am looking for the most efficient way to route inspectors from restaurant to restaurant. This is something like the travelling salesman problem, but with additional complexities. We have 75 inspectors, making less than 8 hour round trips from 17 different offices, inspecting ~18,000 businesses at 1.25 hours an inspection (not counting travel time). These businesses either need inspected annually, 3 times a year, or 5 times a year depending on their assigned risk level. Does anyone have any ideas/suggestions for this? I am talking to Telogis and Viamente as far as getting quotes for using their products but am interested in looking into all avenues. It's possible this is not on the right Stackexchange, if you know of a more appropriate venue I'd appreciate knowing that as well.


title: How to generalize a polygon file and maintain topology in QGIS?

text: my problem: I have a polygon-file in shape-format, which I need to generalize. How can "split" the polygons into lines, generalize and than "rebuild" the polygons, so that I can ensure, that there are no topological errors like "slivers". Please give a rookie some adive howto do this in quantumgis. Thank you very much Ps: I already tried: Vector -> Geometry Tools -> polygones to lines Vector -> Geometry Tools -> simplfy geometries here I can already see on the screen, that I am producing overlaps and gaps nevertheless now I try Vector -> Geometry Tools -> lines to polgones and here the program show an error and crahes ![enter image description here](http://i.stack.imgur.com/i6Wov.png)


title: Trying to convert to EPSG:27700 (Ordnance Survey) projection, unable to do so

text: I'm trying to code an interface that allows a user to draw a vector shape on openstreetmap, which is then converted to a WKT string, processed by the server and sent of as a WFS query to emapsite (which provide the Ordnance Survey data, including a veru detailed vector layer of the UK). However the issue I'm having is that when I draw a polygon around all of the UK (which should obviously return data), I'm not getting data returned (as in the server says there are no vectors there). However when I draw a polygon in the south- west corner of Africa (just off the coast of Nigeria), I'm getting a huge amount of data returned. Ordnance Survey uses EPSG:27700 (OSM uses EPSG:90013), which I don't really understand, but I imagine that I'm not doing projection conversions that are necessary. code below: http://pastebin.com/wi4CUAmc


title: WMS blank map returned

text: I am facing this problem 3th day and I really don't know where is the problem, so please help me. Here is my REQUEST: http://geoportal.cuzk.cz/WMS_SM5V_PUB/service.svc/get?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SM_STP_HRBITOV&STYLES=&CRS=EPSG%3A102067&BBOX=-895000,-865000,-623000,-620000&WIDTH=1600&HEIGHT=610&FORMAT=image/png& Here is Capabilities file: http://geoportal.cuzk.cz/WMS_SM5V_PUB/WMService.aspx?service=WMS&request=getCapabilities Can somebody please explain me why blank map is returned? Thank you


title: Obtaining *rendered* OSM data

text: I'd like to produce a simple walking map with some points of interest. Ideally, I'd like to use the rendered OSM data, as you can get from the page. Using the QGIS OSM plugin, I can only download a fraction of my area, and it is tricky to download a few PNGs at a high enough zoom level and stitch them together from OSM's "Export" tool. I only need to do this a few times, so I don't really see much point in downloading all the XML/Shapefile data, and re- rendering with Mapnik. Is there an easy way of getting hold of the already- rendered raster data fro a given region at a given zoom level in a format that I can easily use with QGIS/GDAL/OGR?


title: How to get the data type of each column from a postgis table

text: I need to get the column data types of all columns in a table, including the geometry types. What I wanted to know is if there is a function or SQL that gives something like this:               column_name | data_type     ------------+--------------     gid         | integer     descr       | character varying(32)     class       | character varying(10)     area        | double precision     geom        | geometry(Polygon,3763)      From a few answers on stackexchange and gis.stackexchange I know that I can get some of the information with the following query:               SELECT          g.column_name,         g.data_type,         g.character_maximum_length,         g.udt_name,         f.type,         f.srid     FROM           information_schema.columns as g JOIN          geometry_columns AS f               ON (g.table_schema = f.f_table_schema and g.table_name = f.f_table_name )     WHERE         table_schema = 'my_schema_name' and         table_name = 'my_table_name'      Result:               column_name | data_type         | character_maximum_length | udt_name | type    | srid     ------------+-------------------+--------------------------+----------+---------+------     gid         | integer           |                          |          |         |     descr       | character varying | 32                       |          |         |     class       | character varying | 10                       |          |         |     area        | double precision  |                          |     geom        | USER-DEFINED      |                          | geometry | Polygon | 3763      But, is there a proper more practical way to retrieve the information in the format I need? Or must I enter the "world" of `CASE WHEN` structures and string concatenation to gather all column attributes in one single column in that format? My fear is if a non expected data type surprises me by needing other attribute from the information_schema.columns table. I.e, in the example table before, I did not used any `numeric (15,2)` data type, that would need to use another attributes (numeric_precision and numeric_scale) to be parsed by a CASE WHEN.


title: How to get ArcMap reference without launching the application

text: I need to get hold of ArcMap properties in order to change the data source of the layers within it. I have developed a similar utility earlier utilizing ICommand, but I was curious if I could get a reference to ArcMap application outside of it and change the layer sources and display the new map. I am using the following code but it instantly launches ArcMap and after that I have no hold on the application i.e. if I terminate my application, ArcMap is still running in background. Is it possible to get a silent, invisible reference to ArcMap? Thanks in advance!               IApplication pApp;     IDocument pDoc;     IMxDocument pMxDoc;     string MXDPath = "C:\\NORMALTEMPLATE\\OLD_Normal.mxd";          pDoc = new MxDocumentClass(); //ArcMap launches here     pApp = pDoc.Parent;          pApp.OpenDocument(MXDPath);     MessageBox.Show("Doc opened");


title: How do you Select an existing vector Layer and Edit it, using python?

text: I'm developing a small plugin that needs to be able to select a vector layer programatically based on what layer the user has selected. However I'm not sure how to do so. I believe I found the two elements I need to make it work: setCurrentLayer and selectedLayer from the QgisLegendInterface, but I'm not certain how they work or how they are implemented, I cannot either find any examples of them other than in the API. Am I close? Or is there a better/proper form of doing this?


title: Select by attribute

text: In my shapefile, I have a bunch of features and one field is a integer numerical attribute. I'd like to select all of the features whose value is a multiple of another value. For instance, I have a field with typical values of 100, 101, 102, 103, 104, 105, 106, etc. (up to 1000) and I would like to select all features with a value that is a multiple of 5.


title: Mapserver: labels and tiles

text: I'm developing map with Mapserver to be rendered in tiles and testing output with Ka-map. Most of the time labels split on two tiles are matching, but sometimes they're not.   I configured Ka-Map like               $tileWidth = 256;     $tileHeight =256;     $metaWidth = 4;     $metaHeight = 4;     /* $metaBuffer = Buffer size in pixels to add around metatiles to avoid      * rendering issues along the edge of the map image      */     $metaBuffer = 256;      but even large buffer of 256 pixels doesn't protect label from occasional splitting.   What good practices do exist to make sure labels are placed deterministically?


title: 'Exceptions.NameError' > "Name 'sheet_count' is not defined" Error in Python Code

text: I know this is something simple but for some reason I am having trouble remembering. In my code I am getting a name is not define for sheet_count. This means I would have to define at top correct. For Example: Is it sheet_count = true? My code is below:               import arcpy, os, sys, string     from os import path as p     from arcpy import mapping as m          #Make parameters for people to choose mxd and folder of PDFs to compare          mxdList = string.split(arcpy.GetParameterAsText(0), ";")     pdf_path = arcpy.GetParameterAsText(1)          #Use Search Cursor to go through Attribuite Table to get Sheet number info          for mapDoc in mxdList:         arcpy.AddMessage(mapDoc)         mxd = arcpy.mapping.MapDocument(mapDoc)         for lyr in m.ListLayers(mxd, "*, Grid"):             if lyr.description == "*, Grid":                 max_list = []                 rows = arcpy.SearchCursor(lyr)                 for row in rows:                     max_list.append(row.Sheet_ID)                 sheet_count = max(max_list)                 arcpy.AddMessage('Sheet count: %s'%sheet_count)                 print 'Sheet count: %s'%sheet_count             else:                 arcpy.AddError('No Layers in %s match data source'%mapDoc)


title: How do I find the points of intersection between a small circle and a great circle?

text: I only actually need to intersect a great circle _arc_ i.e. the shortest path between two points, with a given small circle. I found this tutorial for intersecting two great circle arcs, but unfortunately it doesn't cover small circles.


title: How to pass a layer to a thread in PyQGIS?

text: In my plugin, a user can select one of the loaded layers from a combobox, and then proceed to do further calculations. These calculations are done in a separate thread, "WorkerThread". My question is that how can I pass the selected layer to my thread? Since this layer is variable and is based on the selected layer in combobox, I'm not able to initialize it like this:               def initGui(self):         ...         self.dlg.ui.find_cc_Button.clicked.connect(lambda: self.find_cc(layer)) # won't work, since layer is selected during "run".      "run" segment of my code, where user selects the layer from combobox is this:               def run(self):         # show the dialog         self.dlg.show()         self.dlg.ui.progress_bar.setValue(0)         self.dlg.ui.progressBar.setValue(0)         self.dlg.ui.layerCombo.clear()         self.dlg.ui.lts_combo.clear()              layers = QgsMapLayerRegistry.instance().mapLayers().values()         for layer in layers:             if layer.type() == QgsMapLayer.VectorLayer and layer.geometryType() == QGis.Line:                 self.dlg.ui.layerCombo.addItem( layer.name(), layer )       And the thread is:               def find_cc(self,layer):              self.WorkerThread.start() # I want to pass the layer to this thread          class WorkerThread(QThread):     """docstring for WorkerThread"""     def __init__(self):         super(WorkerThread, self).__init__()         # self.arg = arg          def make_column(self,layer,name):         index = layer.fieldNameIndex(str(name))         if index == -1: # field doesn't exist             caps = layer.dataProvider().capabilities()             if caps & QgsVectorDataProvider.AddAttributes:               res = layer.dataProvider().addAttributes( [ QgsField(str(name), \                 QVariant.Int) ] )             layer.updateFields()             return 0  # so I know if the column already existed or did I create it for the first time         return 1               def run(): #does the job         self.find_connected_components(layer)


title: Creating separate Polylines from Points using ModelBuilder and Python Script Tool?

text: I have a list of coordinate points (long and lat) of 4 locations (points 0, 1, 2, 3) saved as a feature class shapefile. I want to create a polyline that will connect points 0 and 3 and create a second polyline that will connect points 1 and 2. I've found several example scripts for creating points to polylines, but none that split them. Does anyone have any ideas that might work?


title: Google Map API 3 with OSM overlay

text: I'm using GeoXml3 to parse XML files in Google Map and would like to add OSM as custom maptype. I fiddled with this example but it hides completely the original Google layers. I've found a couple of other examples but unfortunately all of them are using the V2 API.


title: What are the pros and cons of using the ArcGIS Viewer for Flex vs. creating your own GIS application vs creating your own framework?

text: Simple question, but perhaps hard to answer. I have broken the title question into three distinct questions in order to clearly identify what we are considering for options each time we build an GIS web application in Flex.   1. What are the pros and cons of using the ArcGIS Viewer for Flex application as a starting point for GIS applications?    2. What are the pros and cons of creating your own application from scratch using the ArcGIS API for Flex? Note: This would be a focused application. In this case it is not intended to be a reusable application from which other applications could be configured / built.   3. What are the pros and cons of creating your own custom framework using the ArcGIS API for Flex from which multiple applications could be built? While the primary goal would be to build an application that meets the customers needs a secondary goal would be to have a reusable framework from which other applications can be configured / built, much like the Flex Viewer. Unlike many other decisions the choice to use one method over the other is not an easily reversible one. Assume you will have to live with this decision for a very long time (years).


title: ArcGIS Server 10.1 Integrated Windows Authentication with Token based

text: In ArcGIS Server 10.1 on Windows 7 and IIS, if wanting to secure services within an Intranet with Active Directory is that Integrated Windows Authentication http://resources.arcgis.com/en/help/main/10.1/index.html#//015400000517000000 ? If using Integrated Windows Authentication, does that dictate that web applications wanting to consume the services re-code their app if they do token based access? I'm involved with setting up AGS 10.1 and configuring and will have numerous 3rd party web app/desktop apps wanting to consume services. Base question is once configured are we setting the policy that any app that want gis services have to do it one way or the other?


title: how to find the point in a country furthest from a boundary

text: i'm trying to find the point in south africa furthest from the sea or a landward boundary. my method thus far:   1. i obtained a first guess at the centre of the inscribed circle and the three segments of boundary this circle might touch using a compass and 1:2 500 000 map.    2. i then clipped these three segments of boundary from a data set of boundary coordinates supplied by national geo-spatial information (WGS84; decimal degrees; 9 decimal places; roughly one point for every 40 metres of boundary).   3. i converted the boundary coordinates to radians and used cartesian equations (distance between two points and centre of a circle through three points) to iteratively improve the centre of my circle and the three boundary points it touches. (obviously, i'm ignoring curvature; i'm effectively finding the centre of the base of the spherical cap defined by the inscribed circle -- i.e. a point approximately 16km underground.) my questions:   1. is it reasonable to plug lat-long coordinates expressed in radians into cartesian equations to find the centre of the base of the spherical cap defined by the inscribed circle (given that the surface distance represented by one radian varies with latitude, and that the base of the spherical cap in question is roughly 800km in diameter and spans latitude 30S)? what order of error might be involved here?   2. is it reasonable to assume that the three boundary points tangential to the base of the spherical cap will likewise be the three boundary points closest to the actual point in the country furthest from a boundary?   3. i have the feeling i'm not even asking the right questions here. i guess there must be a better way to find the point in a country furthest from a boundary ... cheers, william


title: Is pre-1990 historical block-level census data available online?

text: I'm looking for historical block-level data for Chicago. According to this Census bureau document, historical block-level census data exists for major cities going back to 1960, and even back to 1940 (though only for housing statistics.) Only 1990 - 2010 is available on the census website. Do digital copies of historical census block boundaries and stats exist? Are they available anywhere? If not, where are the analogue versions of this data stored? And what is the most granular level of historical census data that is available digitally?


title: Help identifying a coordinate system

text: I have some "latitude" and "longitude" data that describes the Silverstone motor racing circuit in the UK. The coordinates look something like this:               latitude     longitude     +03124.27986 +00061.07652     +03124.27902 +00061.06950      +03124.27866 +00061.06632      +03124.27734 +00061.05942      +03124.27722 +00061.05852      +03124.27710 +00061.05780      These are clearly not "real" degrees & minutes latitudes and longitudes, since the latitudes aren't even valid, and the longitudes are clearly wrong for a place in the UK. Any help in identifying the coordinate system used here would be much appreciated! EDIT: Another set of coordinates in the same format, for a Singapore racetrack looks like               latitude     longitude     +00077.47468 -06231.29347      +00077.47699 -06231.29440      +00077.47890 -06231.29519      +00077.47957 -06231.29553      +00077.47959 -06231.29644


title: Hölzel projection - prj.file

text: Is there a prj-file-defintion available for the Hölzel-projection? (http://dl.maptools.org/dl/libproj4/manual.pdf#page=58, page 58)?


title: gdaltindex Limitations

text: I was wondering if there is a size, or count limitation with gdaltindex. I ran the operation from within OSGeo4W via a Python script, and from within QGIS 2.2 (Windows) and had issues processing a directory with 880 GeoTiffs. The file size of the tiffs are ~ 205MB each. However, when I processed a sub- sample (four tiffs) the operation ran perfect, from within QGIS and OSGeo4W. Any help would be appreciated.


title: Displaying only points within a polygon I created in QGIS?

text: Thank you in advance for any advice you might give me. I am new to qgis and believe what I'm asking is probably fairly basic. I am trying to map policing activity in NYC for two areas: one area in the Bronx and another in Manhattan. I want to display ALL of the activity but only WITHIN a couple of square blocks that I am interested in (not the whole precinct - which I have a shapefile for with boundaries already set). My assumption, is that I create a polygon outlining the block areas I want. But then, how do I attach those boundary coordinates to my csv database on police activity (that dataset has x y coordinates representing police activity)? And then, how do I only display the activity within the boundaries rather than all of the activity that exists throughout NYC? Thanks so much for any links, resources or explanations that might lead me to answering these questions


title: Can I avoid setting up a proxy for geoserver, apache?

text: I'm using Apache (port 80) and Geoserver (port 8080) on Windows 7. I guess I need to set up a proxy if I want to use WFS, since they're on different ports. My kind of silly question is, can I avoid it? If I re-set my website on Geoserver's www file, is everything going to work fine? Does Geoserver share the same characteristics with Apache? If Geoserver's www file is a bad idea can you please help me set up a proxy? I never done it before, so if you have a link with an analytical tutorial, it would be great... Thanks,


title: How can I set the base layer of map before initialization?

text: I have been using `map.setBaseLayer(layer);` to switch to my preferred map layer, using the localStorage variable `layer` that is stored in my web browser, upon initialization of the map. But I realize that this is creating unnecessary calls to load the default (first layer) map tiles that I may not need if the preferred layer is the second or third one. Is there a way that I can pre-select the base layer to be displayed (before map initialization), without changing the order of the map layers?


title: How can i merge lines with same name and intersect each other?

text: I have shp file with 50k lines (roads). Many roads have same name like 'high road' . some of them are the same road (they intersect each other) and some are located in different places of the city. I want to merge all lines with same name and that intersect to single road. How can I do that? I have available SQL 2012 spatial with geometries and geographic data of roads and QGIS with shp file, Someone told me in other post to use dissolve, but it looks like dissolve works on polygons only (i don't have selection for line vector layer) Thank you.


title: How to find coordinates for a lot of locations based on the town or village name?

text: I have got an Excel spreadsheet with roundabout **15,000 locations** , which are only identified by the towns or villages **name** and the country. I am intending to present these on a world map. Therefore, I need to know, whether there is a tool, which can **find the coordinates** for me based on the location's name. Pasting every name into google maps is no solution due to the high number of locations. Thank you for your help and I wish you a satisfying Christmas!!


title: How to install GDAL/OGR with GRASS vector support?

text: The GRASS 6.4.3Svn for Windows that i have in my machine does not have GDAL/OGR support. Is there any library (libgrass) for Windows to add into existing GRASS installation. If not, how to build or install GRASS with OSGEO4W installer including GDAL/OGR support?


title: How to minimize the number of nodes in a polygon?

text: I'm building a web application that will be transferring polygons from server to client. Initially, the shapefile is coming from an external source that I have no control over. Often, I find polygons with way too many points to describe a straight line. The result is, the size of the payload increases and therefore the fetch takes longer. This obviously isn't an issue with only a few polygons, but it will add up as more shapes are added to my dataset and I'd like to optimize bandwidth. ![too many points - qgis](http://i.stack.imgur.com/DqYHJ.gif) Is there a plugin, or native QGIS function that will take some kind of tolerance/variance as input, and produce (approximately) the same polygon with much fewer actual points?


title: Do I get the same result when merging and clipping DEM geotiff files using different steps?

text: I merged and clipped 12 GeoTif Files with GMTED2010 Data (15 Arcsec) obtained from Earthexplorer to get a large map covering Europa and the Near and Middle East. First I tried Variant A) Raster – Misc – Merge => Input whole Folder – Output as GeoTiff File – Load Into Canvas. After that I had planned to clip the resulting merged Layer using a Rectangle Polygon as Mask Layer. Sadly the program crashed during merging so I then tried this way: Variant B) Raster – Misc – Build Raster Catalogue => Input whole Folder – Output as vrt – Load in canvas I then clipped this vrt-Layer using my rectangle square and loaded the resulting geotiff File into the Canvas. Is there a difference between the resulting merged and clipped DEM geotiff File or do I get the exact same Geotiff File no matter which variant I use to produce it?