title: How to connect two line ends in QGIS without resorting to other software?

text: I have two adjacent (ends meet to make one line) polylines in QGIS. I need to properly connect the two ends of the two sections. I've tried the obvious (and less obvious too) but QGIS won't forget that they started out as two lines. On this occasion my overall aim is to create a GPX track using a set of currently unconnected lines with meeting ends (i.e. unconnected in the data, but end to end geometrically). I guess that QGIS needs to understand how each line is meant to connect together - after all it would still be a valid GPX track with one section of line in the track, followed by another random section, and then a third, and so on. Typically the issue is that QGIS does just this - so that the GPX track isn't one continuous line (I get separate 'trkseg' sections of track). To (try to) make one track line I have:   * used a python script to reverse any lines drawn in the wrong direction so that they all point the same way;   * run the Topology Checker plugin to look for 'dangling ends'; and with snapping turned on dragged the offending ends apart and then brought them back together.   * used 'merge features' to make the two (or more) lines into one. and then saved as a GPX track. I've done these things repeatedly and in differing orders and with more than one set of data. It rarely works. What does work is having taken the above steps - to save the data as a MapInfo TAB file, load it in MapInfo, use MapInfo's 'combine' command, save, re-open in QGIS, save back into a SHP file... then create the GPX track (which often works). My interpretation of this is that QGIS is never correctly combining the lines into one - even though the endpoints of the lines are correctly orientated (one in the same place as the other). Am I missing something? Is there some hidden "permanently combine" feature which is stronger than "merge features"? I've read (I think) that it might be possible to use GRASS to do this - but as I already have a solution using a different package I haven't pursued this. I was hoping for a simple solution.


title: Error:getFeatureFromEvent is not a function - Openlayers Multiple layer selection

text: OpenLayers.Control.SelectFeature works well when I try to select single layer from map. But it throws error as **getFeatureFromEvent is not a function** when I pass a collection of layers to it. Following is my detailed code.               var TapestryLayers = [];          fontColor = "#FF0000", fontSize = "13px", strokeColor = "#FF0000";     for(var i=0;i<data.tapestryShapes.length;i++){      var geom = new penLayers.Geometry.fromWKT(data.tapestryShapes[i].the_geom).transform(geographicProj, sphericalMercatorProj);      var fencefeature = new OpenLayers.Feature.Vector(geom, {}, {                                                 "strokeColor": strokeColor,                                                 "fillOpacity": 0.04,                                                 "strokeWidth": 1.5,                                                 "strokeOpacity": 0.3,                                                 "pointerEvents": "visiblePainted",                                                 'label': data.tapestryShapes[i].name,                                                 'fontColor': fontColor,                                                 //'fontSize': fontSize,                                                 'fontFamily': "arial, monospace",                                                 'fontWeight': "bolder",                                             });           fencefeature.attributes = {                                  name: 'tapestryDataLayer',                                  type: 'tapestry'                                 };      var tapestryDataLayer = new OpenLayers.Layer.Vector(data.tapestryShapes[i].name,                                                     {'displayInLayerSwitcher':false}                                                     );      map.addLayer(tapestryDataLayer);      map.getLayersByName(data.tapestryShapes[i].name)[0].redraw();      tapestryDataLayer.addFeatures([fencefeature]);      TapestryLayers.push(map.getLayersByName(data.tapestryShapes[i].name));       }     //select = new OpenLayers.Control.SelectFeature(tapestryDataLayer,{clickout: true, multiple:true}); //But only for single layer as expected!      select = new OpenLayers.Control.SelectFeature(TapestryLayers,{clickout: true, multiple:true}); //Get error getFeatureFromEvent is not a function (but it still a vector layer!)       map.addControl(select);       select.activate();            tapestryDataLayer.events.on({                'featureselected': function(feature) {                                  alert("selected");                                           },                 'featureunselected': function(feature) {                                   alert("unselected");                                           }                                       });      Any suggestions as why its giving that error. I have searched for same error and same functionality on stackexchange but that couldn't help. Thanks


title: How to calculate percentage of property's perimeter affected by flooding (excluding internal walls)?

text: I am trying to determine the percentage of a polygon perimeter that is intersected by another. In short, I am carrying out some flooding analysis in ArcGIS 10 and am trying to determine properties that   1. are directly affected by the flood mapping and   2. the percentage of a property's perimeter that is affected excluding any internal walls.  I'm using Ordnance survey master map data for building polygons. Any guidance on how to achieve this would be much appreciated.


title: Open Source .Net GIS software end-to-end

text: I need to provide a client with an end-to-end web-based GIS solution running on IIS. I will admit to not knowing much about the subject, other than having developed a Google Earth enabled version of the project. However, I need to develop a version that works offline, in the sense that I have an internal network without internet connection. My problem is twofold:   1. I don't know what I need, or rather I don't know how many separate layers make up a GIS server. For instance, I thought QGIS Server by itself would be sufficient, but I still need to get my hands on a spatial data store (or do I?) and it seems that QGIS will only run on a Windows environment if IIS is disabled.   2. I don't know where to get my hands on the various components I don't know I need. Has anyone implemented a .NET web-based GIS service utilizing Open Source technology? UPDATE: The features I need are:   1. Javascript API   2. Ability to easily import custom tiles, e.g. for places like Angola that don't have low-level tiles   3. Ability to add and remove place markers with balloons(like google maps pins with info window popup) using javascrip API


title: mapbasic, create polygon, insert node

text: Hi I try to create a parsel table from coordinates. My scripts are below. Why it does not show the results on the map? Thank you in advance!               Create Table "parsel" (Parsel_Adi Char(10)) file "D:\data\parsel\parsel.tab" TYPE NATIVE Charset "WindowsLatin1"     Create Map For parsel CoordSys Earth Projection 8, 28,"m", 33, 0, 0.9996, 500000, 0     Map from parsel     Set Map Layer parsel Editable On     Create Region Into Variable obj_parsel 0     Alter Object obj_parsel Node Add ( 417000, 4401000)     Alter Object obj_parsel Node Add ( 418000, 4401000)     Alter Object obj_parsel Node Add ( 417000, 4403000)     Alter Object obj_parsel Node Add ( 418000, 4400300)     Insert Into parsel(Object, Parsel_Adi) Values (obj_parsel, "SEMIH")     Commit Table parsel Interactive


title: unable to upgrade gdal-bin from 1.9 to 1.10 on ubuntu 12.04

text: I am trying to install import files to cartodb installed on local machine. But import fails with these errors like               Error 2001 unable to load data     Error 3009 SHP normalization error      Also the ruby console says the ogr2ogr error. So tried to update **gdal-bin 1.10.x**. Current version of gdal-bin is **1.9.2**. But during the sudo-apt get install gdal-bin, these error appears               gdal-bin : Depends: libgdal1h (>= 1.10.0-1~precise1) but it is not going to be installed      Tried to add gdal 1.10 ppa to system. Also removed older version. But not successful.


title: How to apply the four colors theorem in a polygon map in ArcGIS/ArcToolBox automatcally?

text: I need to apply the four colors theorem in a polygonal shape in a way that I do not need to choose manually each color to put in each region. I wish to knows if there is any extension, plug-in, script or data base that may be used with ArcGIS and ArcToolBox to do it mathematically or programatically, so I could use it for now on with every map I come to create. ![enter image description here](http://i.stack.imgur.com/ZPCtw.png) ![enter image description here](http://i.stack.imgur.com/x2rFQ.png) ![enter image description here](http://i.stack.imgur.com/kmBZx.png)


title: Extract lat/long from point type column

text: In my table I have a column "location" that is a POINT. I tried this:               SELECT ST_X(location), ST_Y(location) FROM locations;      But I still get the error message: ERROR: function st_x(point) does not exist


title: How to add attribute in shapefile guidedly

text: I have a point layer (31 by 31 grid = 961 points) created by `fishnet tool > feature vertices to point tool` (but with duplicate points). `NOx Concentration` data are in an Excel file. I want to populate an `NOx Concentration` field in the point layer using the _corresponding_ values in the Excel file automatically. How can this be done? ## Excel File ![My excel file screenshot](http://i.stack.imgur.com/fkkd9.jpg) ## Map ![Arc Map Image](http://i.stack.imgur.com/a0ssh.jpg) I am writing this question in connection with my previous post.


title: MapInfo Region Selection

text: We have large amount of data so to filter it region wise and collect info in that specific region we use region filtering in Actix. So the only data is filter and displayed in that filtered region which is unique and all the statistic is based on that filtering. Can we do that in Mapinfo if yes please answer. Thanks for your time.


title: QGIS - Is it possible to increase Output Field Name?

text: When creating a new Column, the default value for the length of the name is 10. Using the Field Calculator via the Model Builder, I want to include the names of regional areas into the columns which has more than 10 characters. Can the length of the **Output Field Name** be increased?


title: Is there an easy way to refresh / update a raster tif file in GIS?

text: Apologies for the newbie / basic question, but it seems to be a fairly simple concept which I can't get to work properly, and I'd be really grateful for a pointer or two! I have a load of different georeferenced raster tif maps loaded into ArcGIS (9.2). If I update the content of one of these maps in a program outside of GIS, (i.e. photoshop or equivalent,) I want the most up to date map to show in ArcMap the next time I open it. It's linking to the file in windows, so surely it's straightforward enough to refresh the link? Right now, I don't know how to do this without either removing the map from the GIS file and re-adding it, or navigating to the file in windows and deleting the .aux and .rrd files, which is annoying. Surely there is a setting or an easy function I can use to do this? (Please say yes?!) I would be much obliged if anyone could offer any advice. Thanks, Lisa P.S. I am just learning to GIS, and have no Python experience as yet, so please be gentle!


title: Nested loops for raster calculation arcpy

text: I am trying to figure out a way to input corresponding raster datasets from two separate folders into an equation using arcpy. I think it requires some sort of nested loop but I can't figure out how to do this so that only corresponding rasters are inputted into the equation. For the sake of this example lets say I have two folders with rasters listed in them as follows:               Folder 1        Folder 2     WC_top_2008     WC_int_2008     WC_top_2009     WC_int_2009     WC_top_2010     WC_int_2010     WC_top_2011     WC_int_2011     WC_top_2012     WC_int_2012      If I wanted to multiply the corresponding years in each dataset by one another (i.e. WC_top_2008 * WC_int_2008, WC_top_2009 * WC_int_2009 etc.) how would I do this in arcpy? I am using ArcGIS 10 with spatial analyst. Thank you.


title: How to Retrieve field name from Uploaded Excel file using ArcGIS Javascript?

text: I published GP Service which involves uploading data to the server using Upload feature of 10.1. I wrote the code for getting the itemID of the uploaded file. I wanted to have user select the X and Y field from the excel file. How can I access the field names from the uploaded excel file before running GP Task? Any pointers or code samples will help. Thanks in advance Ann


title: ArcPy PointToRaster memory locking error solution?

text: RasterToPolygon ArcPy tool fails on generated raster (PointToRaster) with large(?) dataset. Not sure the size threshold, but fails on 78559 x 33544 raster, completes on smaller datasets so I know syntax, variables, and logic are sound. I feel it is some kind of size triggered memory lock, because if I stop the interpreter and run the tool on the raster dataset on disk it also completes. Image for further reference... ![enter image description here](http://i.stack.imgur.com/pjdaQ.png)


title: Iterate zonal statistics table model builder - add field and name each row by raster

text: I have a model builder model to run zonal statistics as a table on multiple rasters and merge the output table into one table However, I would like to add a field to the merged output table that identifies each each record by it's raster name. Since I'm new to figuring this out, can someone tell me what steps I need to add in order to make that happen? See image attached. Thanks, ![enter image description here](http://i.stack.imgur.com/pwks5.png)


title: Copy ActiveView between AxMapControl

text: i'm starting to work with GIS engine, ESRI v9, and i have troubles with axMapControls. axMap is the main AxMapControl that is shown. I have a static function that scales an image and needs axMap.ActiveView but this can't be declared as static, so (i guess) i should copy current axMap in a new AxMapControl that could be static in order to use it in my func; something like this:               static void foo(){         //do stuff         GetScreenCoordinatesFromMapCoorindates(position, aux_axMap.ActiveView);      So i set a timer with the following:               aux_axMap = axMap;      And in the function "foo":               GetScreenCoordinatesFromMapCoorindates(position, aux_axMap.ActiveView);      But this is not working because aux_axMap is not getting axMap.ActiveView cause when i zoom up/down aux_axMap doesn't update and the function `GetScreen....()` is returning same values. How should i do this?


title: How to format an attribute of type esriFieldTypeDate from a FeatureSet Graphic

text: Hoping someone can tell me if I'm going about this wrong, but I'm executing a QueryTask / Query on a map service and getting back a FeatureSet as expected. Within the attributes of each Graphic, I have two fields that are of type esriFieldTypeDate and their values are the long integer value representing the date. My problem is I'm not seeing it formatted in any way when it's displayed in my dynamically created DataGridColumn. I'd like to see it in MM/DD/YYYY format. I was hoping to check the attribute's data type and if it was equal to 'esriFieldTypeDate' - run some formatting / conversion method on it. But, I can't even get the data type of the attribute. There is no spatial querying going on here, this is purely based on attribute data entered on a form in the application. The last if statement of the sample below is an error - just showing my thought process. Also, the "resultFields" variable is an array of field names. I use that to allow the field order to configurable. Should I not be using a QueryTask for this???               function onResult(featureSet:FeatureSet, token:Object = null):void     {     ...         for each (var gr:Graphic in featureSet.features)         {             for (var i:int = 0; i < resultFields.length; i++)             {                 var hg:HGroup = new HGroup();                 var lab:Label = new Label();                 lab.width = 200;                 var dgc:DataGridColumn = resultFields[i] as DataGridColumn;                      lab.text = dgc.headerText + ": ";                 hg.addElement(lab);                      var valLabel:Label = new Label();                 if (gr.attributes[dgc.dataField].type != "esriFieldTypeDate")                 {                     //Not a date - display text value                     valLabel.text = gr.attributes[dgc.dataField].type;                 }                                                else                 {                     //Attribute is a date - format it to MM/DD/YYYY                     valLabel.text = gr.attributes[dgc.dataField];                    }                      // ... add value label to group, add group to column, add column to grid etc.


title: How to upload a list of 5-digit zip codes to CartoDB?

text: I'm playing around with the free version of CartoDB right now. Trying to upload both a CSV and excel document with a list of 5-digit zip codes. I have received a "there was a problem georeferencing your data". Any idea why this might be? Thanks.


title: Best practice for loading shapefiles into ArcSDE geodatabase using web application

text: We are developing a web GIS using Flexviewer and we want to allow users to upload a shp file and insert it into a ArcSDE geodatabase. I think there are two steps at least.   1. upload the shp file (zip) to a directory on the server   2. call some function to load it to the ArcSDE geodatabase I am not sure what is the best practice for this requirement. using ArcObjects or Geoprocessing service (ArcPy)? Any advise is appreciated.


title: Problems installing libgdal-ecw package

text: I have problems to visualize ecw files for qgis 2.3.0 using ubuntu 12.04 LTS. From terminal, following the instructions for qgis 1.8 found at ECW for QGIS 1.8 in Ubuntu 12.04, because launching the command sudo gdal-ecw-build /usr/local/. I added the repository ppa:ubuntugis/ubuntugis-unstable and then did these steps: wget http://de-mirror.org/gentoo/distfiles/libecwj2-3.3-2006-09-06.zip unzip libecwj2-3.3-2006-09-06.zip cd libecwj2-3.3 ./configure make sudo make install sudo gdal-ecw-build /usr/local/ Launching the last command I receive some warnings like this: " ecwdataset.cpp:1501:15: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]" and then an error message: " /usr/local/include/gdal_priv.h:225:19: error: ‘char** GDALOpenInfo::papszSiblingFiles’ is private ecwdataset.cpp:2526:48: error: within this context" finally (in italian because I have an italian version of Ubuntu 12.04 lts): make: _*_ [ecwdataset.o] Errore 1 ...anybody can help me? Thanks in advance Paolo


title: How do you show the "Select By Location" dialog window when it is open but behind other ArcMap windows?

text: My "Select By Location" dialog window keeps falling behind other windows(such as my table window) that I have on my second monitor. So far, the only way to fix it is to move the table window to the first monitor which allows me to get the the "Select By Location" dialog window, and then window juggle. The "Select By Location" dialog window doesn't create an item in the taskbar which would be useful (I could right-click and move). Or even just having the window raised to the top when I select the "Select By Location" menu-item again would help. Are there any clever workarounds or maybe a shortcut key to deal with this issue?


title: Hyperlinking multiple PDF's to multiple points

text: Hi I am kinda new to all of this and appreciate any help. (no critisicm please, questions welcomed) I would like to hyperlink a pdf to a point. Each point will have a link to a different pdf. So far it seems that I can do it one by one. Is there a quicker way of doing it? Thanks for your patience with the Newby We're using ESRI ArcMap


title: Passing options to ogr2ogr using python subprocess.Popen

text: I'm trying to programatically import some shpefiles to postgis, since I need to catch the exit code of the process I decided to use subprocess.Popen. I've managed to import the files using:               cmd = ["ogr2ogr","-f", "PostgreSQL", "PG:dbname=idegeo host=localhost port=5432      user=postgres password=postgres",             "-nlt", "PROMOTE_TO_MULTI","-nln","tablename",             "-unsetFieldWidth","shape.shp"]     process = subprocess.Popen(cmd,stdout=subprocess.PIPE)     output = process.communicate()[0]      Now, the only thing I'm having problems is in passing the argument `--config PG_USE_COPY YES` to increase performance (as suggested here). I've tried using `[someargs,"--config", "PG_USE_COPY YES"]` but Popen refuses to recognize the value passed to --config Any help is appreciated


title: Joining Textfile to Shapefile Arc 10 and 10.2 Differences

text: I am trying to join a textfile generated in stata that includes caoname, caoid, and fips. The fips code should match to a shapefile county map that matches to the cntyidfp category.   In ArcGIS 10.0, doing a straightforward join based on these attributes works. However, in 10.2, ArcGIS only allows for the joining of the cntyidfp to the caoname, which is completely useless.   What's the difference between the two versions that makes this difficult?   I have tried exporting to an excel document first, ensuring both the fips and cntyidfp data have the same number of characters, and creating a new table of the textfile within arcgis, but none of these solutions seem to work.   What am I missing?


title: Creating tile packages (.tpk files) in ArcGIS for Desktop using multiple cores?

text: It takes hours or days to create tile package files on my workstation, which has 4 cores, 32G ram and a SSD drive. I noticed only one core is in active use while it is working. Is there a way to tell ArcGIS 10.2 for Desktop to use multiple cores? Also, it would be nice to see what's going on - a preview panel would be nice? Is there any way to keep an eye on what is going on?


title: User experience and design pattern resources for web mapping interfaces

text: What are your personal best practices with designing web mapping interfaces for laypeople? Do you have any resources or libraries for good web mapping design patterns and user experience practices? Edit: Bonus points for answers that include information on UX testing and website analytics!


title: In ArcGIS, how to automatically recalculate a calculation if I modify the Excel data source?

text: There's a database with buildings. Each of them has a few parameters, for example building type codes. In a distinct table each building type codes has 4 attributes. During the Join, the datas will be connected to the buildings, according to the code. The problem is if I modify something in the linked Excel table (the secondly mentioned table), the computation should be automatically refreshed based on the new values. (just like in the case of Excel functions) How can I do this?


title: What are the differences between formats for storing styles?

text: I'm looking for a format for store information about data styles for my web application. There are different formats for describing vector data styles. For example, SLD, Mapnik XML, and there may be others. What are the strengths and weaknesses of these formats? What are the other alternatives?


title: pgRouting pgr_tsp() not returning results when used with SQL parameter

text: I’m trying to solve the travelling salesman problem with locations pulled from a database, e.g.:               mydb=# SELECT id::integer, location_lng::float AS x, location_lat::float AS y FROM addresses where region_id = 76 ORDER BY id;       id   |      x      |     y          --------+-------------+-----------        169 | -114.130995 | 51.039266       1946 | -114.122258 | 51.026517       2693 | -114.135513 | 51.031368       3080 |  -114.12952 | 51.037816       3601 | -114.138326 | 51.037731       3819 | -114.135513 | 51.031326       3831 | -114.128151 | 51.033531       4001 | -114.139744 | 51.031526       5918 | -114.134081 | 51.031903      64222 | -114.122329 | 51.030329      74248 |  -114.13445 | 51.025564      99867 | -114.145568 | 51.033872     104685 | -114.134065 | 51.036855     131470 | -114.138351 | 51.036562     131504 | -114.125255 | 51.031422     231023 |  -114.12083 | 51.038339     398043 | -114.125218 | 51.031422     414612 | -114.122333 | 51.033205     535239 | -114.145643 | 51.034149     (19 rows)      And then I use the function as documented:               mydb=# SELECT seq, id1, id2, round(cost::numeric, 2) AS cost     mydb-# FROM pgr_tsp     mydb-# ('SELECT id::integer, location_lng::float AS x, location_lat::float AS y FROM addresses where region_id = 76 ORDER BY id', 1);     seq | id1 | id2 | cost      -----+-----+-----+------     (0 rows)      Any ideas why I'm not getting any results? Thanks! ...René


title: How to delete duplicate nodes/points in QGIS?

text: I'm new to this site and most impressed with the help available. I have a question regarding removal of duplicate points which has been asked before. However the QGIS solutions suggested do not appear to work for me. I'm using QGIS 2.0.1. I have a file of about 100 attributed features which I have brought over from MapInfo into QGIS. I've run the Vector\Geometry Tools\Check geometry validity which gives about 40 duplicate points/nodes. To remove duplicates I have run the MMQGIS plugin which doesn't remove any and also the Topology Checker plugin using "must not have duplicates" but which tells me it can't find any errors. Other suggestions have been using Generaliser plugin (which I cannot find in the repository)and also the Delete Identical (does not appear in the Data Management tools). Other suggestions are GRASS which I’m not familiar with yet. Any suggestions from the experts as to the best/easiest way of doing this, specifically in QGIS? Thanks Update. The following are the steps: 1) Run Vector \ Geometry Tools\ Check geometry validity 2) Save a file with only these duplicate points 3) Select the layer with duplicates 4) Run MMQGIS\ Modify \ Delete duplicate geometries specifying the selected layer as source and a new file name for the "cleaned" layer 5) Run the geometry validity tool again on the “cleaned “ layer and get any captured duplicate points 6) Compare with original duplicates, no change As regards an image, the only thing I get is about 40 points displaying where the duplicates are.


title: Can QGIS resample Rasters during display?

text: When I was using ArcGIS a user could resample a raster during the display (I think via layer properties under the display tab). This would smooth the appearance of the raster so it wouldn't appear pixelated when the user zoomed in too far, I found this useful when using hillshades. I haven't found a way to do this in QGIS, but can it be done? If not is there a GDAL tool I can use to smooth the appearance of my hillshades so they don't look pixelated when they are getting near the end of their optimum scale. Thanks Ando


title: Is it possible to override default QGIS error messages?

text: From my QGIS plugin I would like to override/replace one of the the default QGIS error message ("The layer layer_url is not a valid layer and can not be added to the map") with my own message. I have looked at the API cant find anything along these lines. Is it possible to achieve this?


title: How can I get a join for pre-existing fields?

text: I am joining some values from a table to a shapefile. The table and shapefile have the same column names and types, but when I join, Arcmap creates a new field, and I just want the join to write the data from the table into the fields that already exist. How would I do that?


title: Are there any examples of creating custom openlayers geometry types?

text: I am working on a project that uses openlayers where IE8 compatability is a requirement. Unfortunately, there seems to be some performance issues with OpenLayers in IE8 when parsing a large number of polygons / points from server side data. I have already developed a custom data exchange format based on JSON (After experimenting with both GeoJSON and WKT) which has yielded a large performance boost, but I still think I need more. One current bottleneck in the application seems to be in initializing points - I have a test case where I am initializing geometries with over 85000 points and it takes a few seconds on a virtualized IE8 box. I will be implementing strategies of lazy initialization, but it also occurred to me that in server side libraries like JTS, multi-point, line string, and linear ring instances may be backed by a packed coordinate sequence - a single array containing all required ordinates http://www.vividsolutions.com/jts/javadoc/com/vividsolutions/jts/geom/impl/PackedCoordinateSequence.Double.html Has anybody seen any examples of using this pattern with OpenLayers, or can anybody provide any insight or advice in implementing it? Aside from creating the geometry type objects themselves, I would need to add support for these types to the renderers. Are there any other places where I would need to add support for them? In short, what are the steps required in order to add support for custom geometry type to OpenLayers? Thanks


title: python script to convert shapefiles in a folder to .kmz's with the SAME shapefile name

text: Pretty new at working with python and am Stuck...I have a "for" Loop that looks in a folder and, based on a criteria (Homes valued over 500k), creates .kmz's with symbology for each of the .shp files in the folder. The script effectively creates the layers based on criteria and applies symbology. However I need help carrying the Name (fNAME) into the kmz's. Can anyone help me complete this script? Regards Script: ~~~~~~~~~~~~~~~~~               import arcpy          from arcpy import env          # Variables     env.workspace = r"C:\Work\TEST\HOMES"          fcs = arcpy.ListFeatureClasses("D*", "Point")          HOMES_Lyr = r"C:\Work\TEST\HOMES\LAYER\HOMES.lyr"          arcpy.env.overwriteOutput=True          for fc in fcs:              fNAME = arcpy.Describe(fc).basename+"_Green"         arcpy.MakeFeatureLayer_management(fc,fNAME,"\"PRICE\" > 500000","","")             arcpy.arcpy.ApplySymbologyFromLayer_management(fNAME,HOMES_Lyr)         arcpy.LayerToKML_conversion(fNAME, ????.kmz, "12000", "false", "DEFAULT", "1024", "96", "CLAMPED_TO_GROUND")      ~~~~~~~~~~~~~


title: ImportError: cannot import name QgsVectorFieldSymbolLayer

text: I am using QGIS 1.7.4., downloaded April 2012, but receive an import error loading the class QgsVectorFieldSymbolLayer. As i understand this class has been added in the Python bindings on 2011-11-15. Does anybody know whether this class has already been included and if so, what else can be wrong? No, while executing a Python script. First line:               from qgis.core import (QgsSymbolLayerV2Registry, QgsDataSourceURI, QgsVectorLayer, QgsVectorFieldSymbolLayer, QgsSymbolLayerV2AbstractMetadata, QgsCoordinateReferenceSystem, QgsMapLayer, QgsField, QgsComposition, QgsComposerLegend, QgsComposerLabel, QgsRectangle, QgsFeature, QgsMapLayerRegistry, QgsGraduatedSymbolRendererV2, QgsSymbolV2, QgsRendererRangeV2, QgsMarkerSymbolLayerV2)      All classes are imported okay, except QgsVectorFieldSymbolLayer.


title: Raster data array output flipped on x-axis using python/gdal?

text: I'm trying to create a raster using the python gdal libraries and I've gotten to the point where data is being output, but the output data is flipped on the _origin point_ 's x-axis. I know I must be overlooking something, but I can't figure out where I'm going wrong. Any ideas? When creating the raster I set the upper-left x/y values, and the array appears to be indexed from the upper-left and continue down to lower-right. In the code below I'm filling the array with the value of the row. When printing the array out it looks like this:               [[  1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.         1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.         1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.   1.         1.   1.   1.   1.   1.   1.]      [  2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.         2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.         2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.   2.         2.   2.   2.   2.   2.   2.]      [  3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.         3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.         3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.   3.         3.   3.   3.   3.   3.   3.]      [  4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.         4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.         4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.   4.         4.   4.   4.   4.   4.   4.]     ...      And this data writes successfully to the raster band. However when viewed in MapWindow GIS, the data appears to go in the opposite direction with the originally set _upper-left_ origin point, appearing as the _lower-left_ value. In other words, the data is flipped on the _origin point_ 's x-axis.               import gdal     import osr     import numpy          OUTPUT_FORMAT = "GTiff"     def create_raster(filename="test.tif"):         driver = gdal.GetDriverByName(OUTPUT_FORMAT)         band_type = gdal.GDT_Byte         number_of_bands = 1              x_rotation = 0 # not supported         y_rotation = 0 # not supported         cell_width_meters = 50         cell_height_meters = 50              (min_lon, min_lat, max_lon, max_lat) = _get_point_bounds() # retrieve bounds for point data                 srs = osr.SpatialReference()         srs.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon                 srs.SetUTM( 54, True) # Set projected coordinate system  to handle meters                      # create transforms for point conversion         wgs84_coordinate_system = srs.CloneGeogCS() # clone only the geographic coordinate system         wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, srs)              # convert to UTM         top_left_x, top_left_y, z = wgs84_to_utm_transform.TransformPoint(min_lon, max_lat, 0)              lower_right_x, lower_right_y, z = wgs84_to_utm_transform.TransformPoint(max_lon, min_lat, 0)               cols, rows = _get_raster_size(top_left_x, lower_right_y, lower_right_x, top_left_y, cell_width_meters, cell_height_meters)         dataset = driver.Create(filename, cols, rows, number_of_bands, band_type) #              # GeoTransform parameters         # --> need to know the area that will be covered to define the geo tranform         # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution         geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, cell_height_meters ]         dataset.SetGeoTransform(geo_transform)         dataset.SetProjection(srs.ExportToWkt())              dataset_band = dataset.GetRasterBand(1)         data = dataset_band.ReadAsArray(0, 0, cols, rows).astype(numpy.float32) # returns empty array               for row in xrange(rows):             for col in xrange(cols):                 data[row][ col] = row + 1              dataset_band.WriteArray(data, 0, 0)         dataset_band.SetNoDataValue(0.0)         dataset_band.FlushCache()         dataset = None # Close file      I've also noticed when I calculate the pixel position for a given lat/lon the _y-value_ results in a negative index, which seems sort-of correct considering that the array is from _upper-left_ to _lower-right_.               inverse_geo_transform = gdal.InvGeoTransform(self.geo_transform)[1] # for mapping lat/lon to pixel     pixel_x, pixel_y = gdal.ApplyGeoTransform(self.inverse_geo_transform, utm_x, utm_y)


title: How to batch interpolate many shp´s to DEM

text: Here´s the problem: I have about 40 shapefiles of contourlines (type line) and I need to interpolate them to create DEM's (tiff, or better some grid format). To do so, some batch scripting would be nice but I just cant figure out a good way. I could do it with the batch tools of QGIS (grass v.bspline would do the job), but I would need to set the parameters for EACH single run. Not very slick at all. So I tried to do it with SAGA CMD in a bash script, but the problem is that the interpolation function needs the input of the area extension and I dont really know how to code that again. It´s my first real contact with coding for GIS or Python but I´ve programmed before in Matlab and know some linux bash stuff so that should not be too much of a problem. So can anyone help me with a simple way to batch script basically what the Interpolation Plugin of QGIS does? Any help much apreciated, Martin


title: Export to new projection in QGIS 1.8

text: I would like to use QGIS to convert a CSV file from UTM to lat/long. I imported the CSV file with the "Add Delimited Text Layer" plugin. This works fine, the data looks like it is projected properly on my map. However, when I try Layer> Save as... and set the Selected CRS to NAD83, it does not convert my UTMs to lat/long. When I go Open Attribute Table, the x and y coordinates are still in meters, not degrees. Any thoughts why this isn't working for me?


title: Label inside top left corner of polygons

text: Does anyone know of a tool or procedure to automatically place a label in the corner of a polygon, rather than at the centroid? I am using Mapinfo and know that I can manually place the label there, but I have a lot of polygons and would prefer to do it automatically, especially as the polygons change every so often. Even better would be a way to do this in ArcGIS, but it doesn't seem to have as much control over labelling.


title: Interpolate current Course value of GPS using previous Course values

text: I have a set of GPS data with Lat, lon and Course. I need to interpolate current value of course, using previous values of Course. The formula I am using is not accurate in curves.In C language I was trying the implementation. Formula I have used is : Current_Course = previous_Course + ((previous_course - previous_to_previous_course)* Ratio) where Ratio is ratio of time difference and GPS update time. Does anyone knows a better formula to interpolate course value. Thanks in advance.


title: How to change line thickness of multiple symbols without altering other properties?

text: I want to change all of my line thicknesses in of my symbology to a 2 within the same feature class. I have about 30 lines. Right now, I have to change it separately with each line. Is there a way to select them all to change them?


title: How to find the distance within which a set percentage of features fall?

text: Someone asked a question earlier about finding the number of features within a certain distance of another feature in ArcView. I have the opposite question: is there any way, in ArcGIS 10, to find the distance within which a set percentage of other features fall? I guess I could create a distance raster and then extract the values to points (the other features) and consult the attribute table. However, I will later need to export the distance to another tool somehow, so the distance would need to be part of the tool's output. Is there any way to find such a distance? In other words, is there a tool that works like the density tools in reverse?


title: How to Add New Freehand Shapes and Data in QGIS?

text: I'm new to QGIS. Can anyone either direct me to a good tutorial or explain? I want to draw freehand new shape layers on to Google Earth to map vegetation types and embedded data for each feature. I am able to draw a shapes as a layer on the map, but I don't understand if or how to customise the attribute table for that layer and every time I try to draw a new shape for that layer it comes up with window saying id. or should I be using one of the different database functions? Thanks for your time. paddy


title: WMS and Raster Data

text: I am calling in some WMS layers using openlayers. One of the layers is a raster. I have it working fine. But how does this fit into the spatial database model? I thought a WMS is a protocol for handling vector data. How raster data stored? Are they within PostGIS? Does Geoserver handle raster data? Thanks so much!


title: How to add a raster legend with min/max values in QGIS 1.6?

text: I calculated/made NDVI with landsat data using ERDAS. Results of NDVI is in black and white map but I need it in colour to present the scale for NDVI with value of pixels in Map Composition. How to make the NDVI map in colour in ERDAS? I made the coloured map of NDVI in QGIS but in Map Composition I can not put the legend of the map with value of pixels. Is there any way to put the legend with values of pixels in QGIS? Thanks a lot in advance. ![NDVI I made](http://i.stack.imgur.com/GEnCH.jpg) ![enter image description here](http://i.stack.imgur.com/oSEBW.jpg)


title: Why does qgis become unresponsive when adding a .sid raster layer?

text: Can someone provide insight into why qgis becomes unresponsive when i try to add a raster layer? Specifically a .sid raster from NAIP imagery. I can no longer add these raster layers to projects, even to projects which already have raster layers from the same source (NAIP *.sid)added previously. The same .sid raster files that will not load now have been successfully added to qgis projects in the past. What I have tried: Uninstalled qgis 1.7 and replaced with qgis 1.7.4 using the OSGe04W installer. Moved raster files from external to internal hard drive. Freed up internal hard drive space (>25% of capacity). Tried multiple raster files. I only have 1Gb of RAM, but this has been sufficient in the past. The raster file sizes are on average around 750Mb. Thanks


title: error printinf pdf in opengeo suite client sdk app

text: I am using OpenGeo Suite 3.0.2 I created an app using Client SDK. Stored the raster and vector data on postgreSQL/PostGIS print tool is not working, giving error i.e. "Java.io.IOException: the byte array is not recognized imageformat."   * @require plugins/Print.js      * @require GeoExt.ux/PrintPreview.js     * @require GeoExt/widgets/PrintMapPanel.js     * @require GeoExt/plugins/PrintProviderField.js     * @require GeoExt/plugins/PrintPageField.js     * @require OpenLayers/Control/ScaleLine.js **These are the dependencies, I included.** { ptype: "gxp_print", customParams: {outputFilename: 'DAA Map-PDF'}, printService: "/geoserver/pdf/", actionTarget: "map.tbar", showButtonText: true } **This is the code I included for print service.**


title: How to convert decimal degrees units to km in R?

text: I am trying to figure out how I can transform my unit of length from decimal degree to km when I import the following simple ArcGIS polygon shape. Ultimately, I would like to have this information available when I assign this spatial polygon to my window in spatstat [i.e., Window area = xx square unit & Unit of length: xx km]. I did came across this post: R-sig-geo This is what I've done so far, but I don't understand why/if I need to do so. Doesn't my div0A SpatialPolygonsDataFrame already has my CRS argument defined? Why do I need to spTransform then?               library(spatstat)     library(maptools)     library(lattice)     library(sp)     library(foreign)     require("rgdal")     div0A <- readOGR(dsn="Archive", layer="Projections")     > summary(div0A)     Object of class SpatialPolygonsDataFrame     Coordinates:          min       max     x -80.00 -57.63430     y  66.25  78.16666     Is projected: FALSE     proj4string :     [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]     Data attributes:            Id      NPAzimutha         UTM20            UTM19           AlberEA            Min.   :0   Min.   :348877   Min.   :349232   Min.   :349162   Min.   :348656        1st Qu.:0   1st Qu.:348877   1st Qu.:349232   1st Qu.:349162   1st Qu.:348656        Median :0   Median :348877   Median :349232   Median :349162   Median :348656        Mean   :0   Mean   :348877   Mean   :349232   Mean   :349162   Mean   :348656        3rd Qu.:0   3rd Qu.:348877   3rd Qu.:349232   3rd Qu.:349162   3rd Qu.:348656        Max.   :0   Max.   :348877   Max.   :349232   Max.   :349162   Max.   :348656           proj4string(div0A) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")     xxlcc <- spTransform(div0A, CRS("+init=epsg:4326"))     res1 <- lapply(slot(xxlcc, "polygons"), function(x)       sapply(slot(x, "Polygons"), slot, "area"))          > summary(unlist(res1))        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.       89.72   89.72   89.72   89.72   89.72   89.72          newW<-as(xxlcc, "owin")     > summary(newW)     Window: polygonal boundary     single connected closed polygon with 29 vertices     enclosing rectangle: [-80, -57.6343]x[66.25, 78.16666]units     Window area =  89.7204 square units      Does this look correct? How can I specify/know my unit of length in km? Moreover, this window is ~Baffin Bay, west of Greenland and I am wondering if I should be using a different datum? E.g. EPSG Projection 4747 - GR96? Although, I will have to create ppp object, for which the coordinates have been taken with GPS satellite navigation system. I would appreciate if someone can walk me through this or suggest an online example or further readings. I read chapter 4: Spatial Data Import and Export in Applied Spatial Data Analysis with R (2008), but I'm obviously very new to all of this. Thank you very much in advance, Aurelie


title: ArcGIS Online:Unable to launch application: Unable to access application's configuration

text: I have a custom configuration for an ArcGIS Online Web Application, and have followed all instructions. I have changed the Web App settings and receive this message. **Unable to launch application: Unable to access application's configuration** My JSON file is configured correctly. What could be the cause of this?


title: Converting kml to shape file via ogr2ogr including all attributes

text: Here is the example kml file I am using:               <kml xmlns="http://www.opengis.net/kml/2.2">     <Document>       <name>My Golf Course Example</name>       <Placemark>         <name>Club house3</name>         <visibility>1</visibility>         <phoneNumber>554545454</phoneNumber>         <Point>           <coordinates>-111.956,33.5043</coordinates>         </Point>       </Placemark>     </Document>     </kml>      I have named it `input3.kml`. I put it, and ogr2ogr in the same folder, then I run this command:               ogr2ogr -f "ESRI Shapefile" output3.shp input3.kml      After that, I open up Quantum GIS and add vector layer, then add `output3.shp`. Then I find the `output3 layer`, right click and open attribute table. How come only two attributes show up?(Name and Descriptio) Reza


title: How to find the service pack level in ArcGIS Addin?

text: I can query the current version from RuntimeInfo.Version. But I would also like to know what service pack is also installed. Do I need to just check the registery key HKEY_LOCAL_MACHINE\SOFTWARE\ESRI\ArcGIS\RealVersion. Or is there a proper ArcObject that allows such a thing?


title: Add Attritube Field to Combo Box on Form

text: I am trying to create a user box that:   1. When the button on the toolbar is clicked, a user form appears with a combo box filled in with a feature's location (Taken from the attribute table - such as I-10 and I-15)   2. Once the user selects a location from the combo box, additional text boxes on the form are populated with additional attribute information of that selected feature (i.e., description, costs, etc).   3. The user can change the additional attribute information and then click update to save the changes.  I have this running in ArcGIS 9.3 using VBA ArcObjects and am now trying to convert it to VB.Net using Visual Basic 2008 Express. I am having issues with my AddItem Sub Routine- it says that my List Count, List, and AddItem are not members of the System.Windows.Form.ComboBox. Suggestions on what to change? **Errors:**               Public Sub AddItem(ByVal sItem As String, ByVal cmbBox As System.Windows.Forms.ComboBox)         ' adds item retaining alphabetic order, unique values         Dim l As Long         For l = 0 To cmbBox.ListCount - 1             Select Case StrComp(cmbBox.List(l), sItem, vbTextCompare)                 Case 0 : Exit Sub ' already in the combobox                 Case 1                     cmbBox.AddItem(sItem, l)                     Exit Sub             End Select         Next l         cmbBox.AddItem(sItem)              End Sub      **Entire Code:**               Public Class TransportationImprovementLinesAttributesForm              Public pMxDoc As IMxDocument              Private Sub cmdCancel_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles cmdCancel.Click             Dim pMxDoc As IMxDocument             pMxDoc = My.ArcMap.Document             Dim pMap As IMap             pMap = pMxDoc.FocusMap                  'Close Form when button is pushed             Me.Close()                  'Clear Selection and Refresh View             pMap.ClearSelection()             pMxDoc.ActiveView.Refresh()              End Sub              Private Sub TransportationImprovementLinesAttributesForm_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load             'Get Improvement Projects FeatureLayer             Dim pFeatLyr As IFeatureLayer             pFeatLyr = GetLayerByTOC("RV_ImprovementLines")                  'Sort Location combobox             Call SortCboFromFeatLyr(pFeatLyr, "LOCATION", cboLocation)         End Sub         Private Sub cmdUpdate_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles cmdUpdate.Click             Try                 Dim pMxDoc As IMxDocument                 pMxDoc = My.ArcMap.Document                 Dim pMap As IMap                 pMap = pMxDoc.FocusMap                      'Get a reference to the editor.                 Dim editor As IEditor                 editor = My.ArcMap.Application.FindExtensionByName("esriEditor.Editor")                      'Get Crossings FeatureLayer                 Dim pFeatLyr As IFeatureLayer                 pFeatLyr = GetLayerByTOC("RV_ImprovementLines")                      'Get Dataset from FeatureLayer                 Dim pEditDataset As IDataset                 pEditDataset = pFeatLyr                      'If not editing, start editing dataset's Workspace                 If editor.EditState <> esriEditState.esriStateEditing Then                     editor.StartEditing(pEditDataset.Workspace)                     'editor.StartOperation()                 End If                           'Get Feature Reference                 Dim pQueryFilter As IQueryFilter                 pQueryFilter = New QueryFilter                      ' Set the where clause                 pQueryFilter.WhereClause = "LOCATION ='" & cboLocation.Text & "'"                      ' Execute the query filter                 Dim pFeatureCursor As IFeatureCursor                 pFeatureCursor = pFeatLyr.Search(pQueryFilter, False)                      'Initialize first (only) feature                 Dim pFeat As IFeature                 pFeat = pFeatureCursor.NextFeature                      'Update attributes with user input                 pFeat.Value(pFeat.Fields.FindField("Location")) = cboLocation.Text                 pFeat.Value(pFeat.Fields.FindField("ID")) = txtProjectID.Text                 pFeat.Value(pFeat.Fields.FindField("Phase")) = txtPhase.Text                 pFeat.Value(pFeat.Fields.FindField("Descriptio")) = txtDescription.Text                 pFeat.Value(pFeat.Fields.FindField("Agency")) = txtAgency.Text                 pFeat.Value(pFeat.Fields.FindField("Cost")) = txtCost.Text                 pFeat.Value(pFeat.Fields.FindField("DateEdited")) = txtDate.Text                 pFeat.Value(pFeat.Fields.FindField("UserEdited")) = txtEdit.Text                 'pFeatureCursor.UpdateFeature(pFeat)                 ' pFeat = pFeatureCursor.NextFeature                      'Stop editing and save edits                 editor.StopEditing(True)                      ' Flag the new selection                 My.ArcMap.Document.ActivatedView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, Nothing, Nothing)                      'Close form                 Me.Close()                  Catch ex As Exception                 Windows.Forms.MessageBox.Show(ex.ToString, "error")                       End Try         End Sub              Public Sub SortCboFromFeatLyr(ByVal pFeatLyr As IFeatureLayer, ByVal strField As String, ByVal CmboBox As System.Windows.Forms.ComboBox)                  'Pass FeatureLayer to Table             Dim pTable As ITable             pTable = pFeatLyr                  'Create TableSort object             Dim pTS As ITableSort             pTS = New TableSort                  'Specify TableSort properties             With pTS                 .Fields = strField  'Field to sort (just one for now)                 .Ascending(strField) = True  'Ascending or Descending                 .CaseSensitive(strField) = False                 .Table = pTable  'Pass Table in                 .Sort(Nothing)  'Perform the sort             End With                  'Get index of Field's position in FeatureLayer table             Dim lngFldIdx As Long             lngFldIdx = pFeatLyr.FeatureClass.Fields.FindFieldByAliasName(strField)                  'Pass sorted Rows to Cursor             Dim pCursor As ICursor             pCursor = pTS.Rows                  'Get first Row in Cursor             Dim pRow As IRow             pRow = pCursor.NextRow                  'Loop through Cursor to add items             Do While Not pRow Is Nothing                 AddItem(CStr(pRow.Value(lngFldIdx)), CmboBox)                 pRow = pCursor.NextRow             Loop              End Sub              Public Sub AddItem(ByVal sItem As String, ByVal cmbBox As System.Windows.Forms.ComboBox)                  ' adds item retaining alphabetic order, unique values             Dim l As Long             For l = 0 To cmbBox.ListCount - 1                 Select Case StrComp(cmbBox.List(l), sItem, vbTextCompare)                     Case 0 : Exit Sub ' already in the combobox                     Case 1                         cmbBox.AddItem(sItem, l)                         Exit Sub                 End Select             Next l             cmbBox.AddItem(sItem)              End Sub         Private Sub cboLocation_Change()             GetAttributesQuery()             SelectMapFeatures()         End Sub              Sub GetAttributesQuery()             Dim pMap As IMap             Dim pDoc As IMxDocument             pDoc = My.ArcMap.Document             pMap = pDoc.FocusMap             Dim pLayer As IFeatureLayer             pLayer = GetLayerByTOC("RV_ImprovementLines")             Dim pFL As IFeatureClass             pFL = pLayer.FeatureClass             Dim pTable As ITable             pTable = pLayer.FeatureClass             Dim pQueryFilter As IQueryFilter             Dim pCursor As ICursor             ' create a query filter             pQueryFilter = New QueryFilter                  ' create the where statement             pQueryFilter.WhereClause = "LOCATION = '" & cboLocation.Text(cboLocation.Text) & "'"                  ' query the table passed into the function and use a cursor to hold the results             pCursor = pTable.Search(pQueryFilter, False)             Dim pRow As IRow             pRow = pCursor.NextRow             txtProjectID.Text = pRow.Value(pRow.Fields.FindField("ID"))             txtPhase.Text = pRow.Value(pRow.Fields.FindField("Phase"))             txtDescription.Text = pRow.Value(pRow.Fields.FindField("Descriptio"))             txtAgency.Text = pRow.Value(pRow.Fields.FindField("Agency"))             txtCost.Text = pRow.Value(pRow.Fields.FindField("Cost"))             txtDate.Text = pRow.Value(pRow.Fields.FindField("DateEdited"))             txtEdit.Text = pRow.Value(pRow.Fields.FindField("UserEdited"))              End Sub              Public Sub SelectMapFeatures()             Dim pMxDoc As IMxDocument             Dim pMap As IMap             Dim pActiveView As IActiveView             Dim pFeatureLayer As IFeatureLayer             Dim pFeatureSelection As IFeatureSelection             Dim pQueryFilter As IQueryFilter                  pMxDoc = My.ArcMap.Document             pMap = pMxDoc.FocusMap             pActiveView = pMap                  'For simplicity sake let's use the first layer in the map             pFeatureLayer = GetLayerByTOC("TransportationImprovementProjects_Lines")             pFeatureSelection = pFeatureLayer 'QI                  'Create the query filter             pQueryFilter = New QueryFilter             pQueryFilter.WhereClause = "LOCATION = '" & cboLocation.Text(cboLocation.Text) & "'"                  'Invalidate only the selection cache             'Flag the original selection             pActiveView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, Nothing, Nothing)             'Perform the selection             pFeatureSelection.SelectFeatures(pQueryFilter, ESRI.ArcGIS.Carto.esriSelectionResultEnum.esriSelectionResultNew, False)             'Flag the new selection             pActiveView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, Nothing, Nothing)              End Sub         Public Function GetSelection(ByVal pFeatLyr As IFeatureLayer) As ISelectionSet                  'Gets selection of feature layer and returns selection set             'Initialize the required variables             pMxDoc = My.ArcMap.Application.Document             Dim pFeatSel As IFeatureSelection             pFeatSel = pFeatLyr                  GetSelection = pFeatSel.SelectionSet              End Function              Public Function GetLayerByTOC(ByVal lyrName As String) As IFeatureLayer                  'This function finds a feature layer based on its TOC Name             'Initalize global variables             pMxDoc = My.ArcMap.Application.Document             Dim pMap As IMap             pMap = pMxDoc.FocusMap                  Dim lyrCntr As Integer             Dim pFeatLyr As IFeatureLayer                  For lyrCntr = 0 To pMap.LayerCount - 1                      'Ensure that the layer is valid                 If pMap.Layer(lyrCntr).Valid = True Then                          'Ensure that the layer is a feature layer                     If TypeOf pMap.Layer(lyrCntr) Is IFeatureLayer Then                         pFeatLyr = pMap.Layer(lyrCntr)                              If UCase(pMap.Layer(lyrCntr).Name) = UCase(lyrName) Then                             GetLayerByTOC = pFeatLyr                             Exit Function                         End If                     End If                 End If             Next lyrCntr                  Return Nothing             'Layer not found. Show a message             MsgBox("Layer " & lyrName & " not found !", vbCritical, "Error")              End Function              Public Function GetFeatSelCount(ByVal pFeatLyr As IFeatureLayer) As Long                  'Returns # of selected featurues of feature layer as long integer                  Dim pFeatSel As IFeatureSelection             pFeatSel = pFeatLyr                  Dim pFeatSelectionSet As ISelectionSet             pFeatSelectionSet = pFeatSel.SelectionSet                  GetFeatSelCount = pFeatSelectionSet.Count         End Function          End Class


title: How to infer a relationship between two spatial points using R

text: I am quite new to the area of spatial statistics, but I'm very interested. For learning and demo purposes, I've created three datsets.   **Dataset - Persons:** This describes individuals at a certain location with a few variables. Please note, that the persons are located in the provided cities. A short explanation:   * POINT_X: X-coordinate of city.    * POINT_Y: Y-coordinate of city.   * city: The name of the city, where they live.   * ill: "1" states that they are ill. For learning purposes, all persons are ill.   * job: If they have a job or not. "1" means: they have one, "0" means they haven't got one.   * disnw: The distance to the nearest waterpoint.   * wID: not relevant. **Dataset - City:** This describes a number of cities including some variables. A short explanation of these:   * city: The name of the city.   * population: The population of the city.   * POINT_X: X-coordinate of city.    * POINT_Y: Y-coordinate of city.    * ill: Number of ill persons in the city.    * notill: Number of healthy persons in the city.    * disnw: The distance (in km) to the nearest waterfeature.    * wID: _not relevant_   * rate_ill: The rate of ill persons in the city.    * rate_notill: The rate of healthy persons in the city.  **Dataset - Waterfeatures:** . Please note that the viallages are on the same location as persons. This is a collection of spatial points, which describes waterfeatures.   * POINT_X: X-coordinate of a waterfeature.    * POINT_Y: Y-coordinate of a waterfeature. geographic overview about the setting (red are persons, blue are waterfeatures, yellow are cities) ![Plot](http://i.stack.imgur.com/hhajK.png) Now I want to check the hypothesis that cities, which are nearer to waterfeatures (so where the variable _disnw_ is lower), have a higher number of ill persons. So is there a correlation between the number of ill persons(`no_ill`)/rate of ill (`rate_ill`) and the proximity to water features (`disnw`). I know, that the datasets are possibly not representative or suitable for my hyptothesis, but for now this fact shouldn't matter.   I've already looked at some functions and packages, but I'm very unsure about a suitable method. Methods, which might be useful (at least from my point of view): semivariogram, variogram, Ripley's K function, G-Function, correlation coefficient. To give you a better overview, I've created example datasets. You can find these here:               persons = read.csv("http://pastebin.com/raw.php?i=3aMGi9Ax", header = TRUE, stringsAsFactors=FALSE)     city = read.csv("http://pastebin.com/raw.php?i=Lk3KXLQT", header = TRUE, stringsAsFactors=FALSE)     water = read.csv("http://pastebin.com/raw.php?i=hQRvMZwE", header = TRUE, stringsAsFactors=FALSE)      It would be awesome to get some input from your side. Maybe you have a tip, how to perform this kind of analysis.   Thanks in advance!


title: ArcGIS Server and SQL Spatial. Examples?

text: In-house we have an ArcGIS Server license that I am somewhat familiar with, and a license for SQL Server. I'm looking for some examples of existing web maps that use ArcGIS Server for the display of geographic data and SQL Server (Spatial) for storage. I've looked into this but haven't found any good/recent examples of this combination, especially which utilize an editing interface for modifying the stored geometries. I'm looking for a few examples and links to documentation that will help me build a good understanding of the capabilities and requirements of these software packages before embarking on the project.


title: Extrapolate data and district names from map

text: My data: soil samples in each district. My question: I want to extrapolate from the soil sample's map the names of the districts in which the samples have been collected. I have already found the centroids coordinates of each district and I have separate data with the names of the districts and their centroids coords. Ideally, I would like to get a table which says _which soil samples are in which district_. Could anybody help me please? Thank you very much!


title: How to achieve equivalent to Dissolve from ArcGIS for Desktop using ST_Union?

text: I am trying to achieve an equivalent to the Dissolve tool from ArcGIS for Desktop using ST_Union from PostGIS but it seems I am not getting expected result. I have one table which has certain attributes with the Polygon Geometry. (like FID, Locstat, Loccnt, Shape) Here is my query :               SELECT c.fid, ST_Union(c.boundaryshape) FROM c Group by c.fid,c.boundaryshape;      Thanks in advance.


title: How to truncate the %Name% string for inline variable substitution

text: I'm using model builder to iterate through feature classes in datasets. With each sucessive process step (each model) I'm using the %Name% variable to name new feature classes in a new dataset. I found that I can't just use %Name% without adding a little to the name of the new feature class; the model didn't seem to like working with two feature classes with the same name. To write the new feature classes to disk, I'm using something like %Name%_AppendedBit. The problem is that with sucessive steps I keep adding more and more bits to the name of the feature class. I tried using Calculate Value with something like Replace(%Name%,_AppendedBit,_JustIdentifyThisStep) to make a variable with which to name the new feature class, but %output_value% doesn't appear to have any value (just "1"). I tried this as a variant and as a string. I was using vb Replace in the "expression" box, not messing with the "code block" box. I could use %n% to name the feature classes, but then I lose useful detail from the original feature class name. This could be worked around by doing a spatial join later. So, the long and the short of it is I'd like to use some kind of string truncation with inline variable substitution.


title: QGIS Python Composer Attribute Table has no content

text: I post 2 Month ago the same question. So I try it a second time with more code. My Problem still is that I add at the composer a attribute table and in this table there is no content. My Code:                   mapRenderer = self.iface.mapCanvas().mapRenderer()         c = QgsComposition(mapRenderer)         c.setPlotStyle(QgsComposition.Print)         x, y = 0, 0         w, h = c.paperWidth(), c.paperHeight()         composerMap = QgsComposerMap(c, x,y,w,h)         c.addItem(composerMap)              composerLabel = QgsComposerLabel(c)         composerLabel.setText("Hello world")         composerLabel.adjustSizeToText()         c.addItem(composerLabel)              tableGes = QgsComposerAttributeTable(c)         tableGes.setVectorLayer(vlayer)         tableGes.setComposerMap(composerMap)         tableGes.setDisplayOnlyVisibleFeatures(True)         c.addItem(tableGes)         ## Here i try this Methode also:         ##c.addComposerTable(tableGes)         ## Same result!              printer = QPrinter()         printer.setOutputFormat(QPrinter.PdfFormat)         printer.setOutputFileName("E:\***\out.pdf")         printer.setPaperSize(QSizeF(c.paperWidth(), c.paperHeight()),         QPrinter.Millimeter)         printer.setFullPage(True)         printer.setColorMode(QPrinter.Color)         printer.setResolution(150)              pdfPainter = QPainter(printer)         c.renderPage(pdfPainter, 0)         pdfPainter.end()      Additional I post the out.pdf as a image. ![enter image description here](http://i.stack.imgur.com/NjfV1.png) I hope someone can help me. Thanx. Greetings Andreas _EDIT:_ Can someone try the code in a own plugin? If there is the same result i can make a bug report.... If not ... i need to search the bug in my plugin....


title: How to render a raster multiple times across the date line?

text: Using EPSG:4326 and WGS84, I have a GeoTIFF with longitude extents (167, 208) - so it crosses the date line. I'd like mapserver to display the image regardless of whether the user's requesting bbox comes from the west or the east, however:               [...]bbox=150,0,250,90[...]      displays the image, while               [...]bbox=-210,0,-110,90[...]      displays a blank image. Everything else unchanged, I'd like to see the same image returned in both cases. Some research brought up this 3+ year old mail from FrankW: http://lists.osgeo.org/pipermail/mapserver-dev/2010-February/009713.html Points #3 & 4 imply that adding PROCESSING "RESAMPLE=AVERAGE" should solve this problem, but I've added that to my mapfile with no luck. Has this changed in the past few years? Is there another way I can get the desired behavior? Is there a setting I could check for that could be preventing the image from mirroring every +/- 360 degrees longitude? Thanks!


title: Can't find MapFrame of GraphicsContainer in PageLayoutControl (C#) ?

text: I am getting a problem: I added element (text, line) in a ActiveView.FocusMap of a PageLayoutControl, then all elements are showed in the MapFrame of PageLayoutControl.GraphicsContainer. If I don't select elements (in MapFrame) I can find the MapFrame by the loop:               IGraphicsContainer graphicsContainer = myPageLayout.ActiveView.GraphicsContainer;     graphicsContainer.Reset();     IElement pElement = graphicsContainer.Next();      while (pElement != null)                 {                     if (pElement is IMapFrame)                     {                           // do something                      }                 }          But, when one element is selected by using "Select Elements" command, the program won't enter the "if" function (can't find MapFrame). Are there anyone know how to get The MapFrame in this case ? It's very urgent Thank you very much !


title: How to read .adf (TIN) map files?

text: I have a set of maps in few formats :   * TIN > >     ├── [4.0K]  info >     │   └── [2.6K]  arc.dir >     ├── [ 708]  log >     ├── [374K]  M-33-45-D-d-4.asc >     └── [4.0K]  m-33-45-d-d-4tin >         ├── [  96]  tdenv.adf >         ├── [1.1M]  tedg.adf >         ├── [ 332]  thul.adf >         ├── [8.2K]  tmsk.adf >         ├── [ 116]  tmsx.adf >         ├── [1.1M]  tnod.adf >         ├── [751K]  tnxy.adf >         └── [188K]  tnz.adf >     * ASCII  > >     ├── [ 12K]  M-33-45-D-d-4_c.asc >     ├── [ 12K]  M-33-45-D-d-4_j.asc >     ├── [ 365]  M-33-45-D-d-4_k.asc >     ├── [1.6K]  M-33-45-D-d-4_o.asc >     ├── [518K]  M-33-45-D-d-4_p.asc >     ├── [114K]  M-33-45-D-d-4_pz.asc >     ├── [109K]  M-33-45-D-d-4_s.asc >     ├── [164K]  M-33-45-D-d-4_sz.asc >     └── [317K]  M-33-45-D-d-4_z.asc >     * TTN  > >     └── [3.2M]  M-33-45-D-d-4.ttn >   I would like to somehow interpret those maps (preferably TIN) for my visualization application but I am quite confused where to start. I know about GDAL library but I am not able to use it on those files to convert them to anything else. Can anyone point me the direction how can I read those .adf files?


title: Generating polygonal blocks using growing lines

text: We first generate a set of **locations** (i.e., points(x,y)), randomly. Then we assign for each point a number showing the **angle** of a line being centred at that point. _The criteria_ for the length of lines is that, the length of lines (all together) grows slowly until a line intersects another line (as shown below). Then that end-point of growing line will be frozen (no further growth) for that line. The final result will be a set of closed areas. **We are looking for algorithms and ideas regarding the implementation of a solution for the mentioned problem.** ![enter image description here](http://i.stack.imgur.com/VcuLc.png) **Update:**   Here we add the result of our implementation in **Python** of the **R** code developed below by `whuber`, the accepted answer. His mathematical solution is very nice, such that even its Python edition took only `186`s (~`3`m) for `500` points, Wow! This is awesome. Note that in this solution there is no further requirement for refining the resulting lines. **Final Note:**   In order to get lines instead of rays you can replace (in _R_ ):               rays <- array(as.vector(cbind(pts,pts+limits*vec)),dim=c(n.points, 2, 2))      with something like (in _Python_ ):               pts += limits*vec     n = len(pts)//2     lines = hstack(pts[:n],pts[n:])      **Really Final Note**   To mention that to make `blocks` from `touching lines` (the product of the post) if you are interested in using of `DCEL` there is a simple (with some bugs) implementation at {here} for **Python** lovers. Also you may follow {this post} for other methods in generating blocks (polygons) from touching (intersecting) lines. ![enter image description here](http://i.stack.imgur.com/scZ39.jpg)


title: Generate polygons from a set of intersecting lines

text: This is a simple and quite common question which has already been asked for different purposes (see this link and this too, for example), here, however, we are looking for _not a software package_ but **algorithms** that we could try to implement say in **Python**. So, as shown below a set of lines are mapped (they are already trimmed, BTW).   **_Algorithms/ideas to generate polygons (as red ones shown)_**? ![enter image description here](http://i.stack.imgur.com/XyJxG.png)


title: How to assign value to field of layer based on spatial query with another layer?

text: I want to assign a value to a field of a polygon layer based on the value of a field on another polygon layer, using a spatial query. My first polygon layer contains buildings and the second contains administrative boundaries. I want to assign the administrative area code based on which administrative area polygon the building is completly contained in. So, on my building layer I've a field named "AreaCode" that I want to update with the Area code value from the Administrative Boundaries layer. I want to do it using the program interface if possible and not by coding a script or macro. Otherwise, does exist an external tool/plugin for ArcGIS that can do this?


title: Styled Layer Descriptor importing Problem

text: I Published a shape file in geoserver and changed that shape file style using Udig via WFS and then SLD file exported to some other place . My problem is while I'm trying to upload that SLD file in Geoserver i'm getting following error is "Oops something went wrong".. Please Clearfy.. Thanks... Satya


title: *Value relation* widget in QGIS

text: I cannot understand what the **Value Relation** widget of the Fields menu in QGIS actually does. For example: I have a point layer with some columns (ID, value1...) and a table (so just a dbf file) with other values (ID, valueA,...). If I open the **Fields** menu of the point layer I can set some options for the selected column: Layer (set to table.dbf), Key column (set to ID) and Value column (set to valueA). Now, if I create a new point, it appears in the row of the widget column a dropdown menu where I can choose the valueA of the table.dbf. And that's ok. But, what is the meaning of the **Key column**? Does someone have an explanation for that? Am I missing something?


title: Shortest path for a given waypoints

text: I have an application where the user selects some "markers" from a google maps (all those coordinates' markers are stored into a postgresql database). Right now, I can draw the route between all selected markers (using DirectionService) but when I have over 10 waypoints the route seems to be drawn in the user selected order (not in the "shortest route" order - I also used optimizeWaypoints param to true and split the route into subroutes of 8 points + start and stop points without success). What I`m trying to do now: when the user selects a marker, I save his selection (marker coordinates) and then all selected coordinates are sent (via ajax) to the server. The received coordinates are in a format something like this: array((lat1,lng1), (lat2,lng2), ...(latn,lngn)) From here, I want to return back to the user all those coordinates ordered by "shortest path" and I dont know how to achieve this. Any idea will be great. Thanks.


title: RESTful Geoprocessing Service Request returns empty Result object without errors

text: This is a bit difficult for anyone here to test as our ArcGIS Server 10.2 instance is behind a firewall. I'll try to present all the necessary information, but feel free to ask for more and I'll adjust the question accordingly. I have a Geoprocessing Service hosted on ArcGIS Server. The service is a very simple Python script (below) that takes in a string of Hydrologic Unit Codes (in this case sub-watersheds [HUC8 level]). The script takes this string, splits it and compares against a table of species occurrence by watershed. This data is compiled into a Python dictionary and returned as JSON with the json module. I'm accessing this service via an ajax request from a web mapping app, but for simplicity's sake I'll leave the JavaScript side out unless it's really necessary. **Python**               import arcpy, json          # files     hucs = r'S:\ES\ES 2007 Efiles\PROP-500 Maps and Photography\PROP-540 Geographical Information System (GIS) Data files\MegaPetition.gdb\HUCs'     species = r'S:\ES\ES 2007 Efiles\PROP-500 Maps and Photography\PROP-540 Geographical Information System (GIS) Data files\MegaPetition.gdb\CandidatesHUC'     data = {}     Watersheds = arcpy.GetParameterAsText(0)     Watersheds = Watersheds.split(",")          cursor = arcpy.SearchCursor(hucs)          for row in cursor:         for watershed in Watersheds:             if row.HUC_CODE == watershed:                 currentWatershed = {}                 speciesList = []                 arcpy.AddMessage("Working on " + row.HUC_NAME + "...")                      whereClause = "HUC8_CD = '" + watershed + "'"                 speciesCursor = arcpy.SearchCursor(species, whereClause)                 arcpy.AddMessage("    Contains " + str(row.NumPetitionedSp) + " At-Risk Species.")                      for animal in speciesCursor:                     speciesList.append(animal.GCOMNAME)                 currentWatershed = {'attributes': [watershed, row.NumPetitionedSp],                                     'species': speciesList}                 data[row.HUC_NAME] = currentWatershed          print json.dumps(data)     del row, cursor      Perhaps I don't understand how the script gets the string of watersheds from the request so that I can use it in my script. With a Script Tool I use `arcpy.GetParameterAsText(0)`, but that doesn't seem to apply here. I've tried to execute this script through the ArcGIS Server GP Service Page, but I still get an empty result: ![RESTful Request](http://i.stack.imgur.com/HyBga.png) Generates the following URL: `http://server- address/arcgis/rest/services/NatureServe/AtRiskGPService/GPServer/NatureServe%20At%20Risk%20GP%20Service/execute?Watersheds=06040003%2C06040002%2C06030003%2C06030002%2C06030005&env%3AoutSR=&env%3AprocessSR=&returnZ=false&returnM=false&f=pjson` With the following Response: ![Empty Response](http://i.stack.imgur.com/2H6bl.png) **Update!** I changed the messaging option for the service from _Error_ to _Info_ , so I can see that the script is taking in the watershed parameter and processing it. It looks like it's just not sending the response back as part of the Result object. The Result object is created by running a GP tool (ex. `result = arcpy.GetCount_management(in_table)`). Is there any way I can add my own data to the Result object if I'm not actually running any GP tools? If not I guess I need to send the JSON through a message instead.


title: How to calculate pixel per meter ration according to zoomlevel in google map(bing map)?

text: I am new here. I need some help. Right now, I got an image and detect two points in the image. I want to calculate the distance between these two points. I got:   * 1 . x,y pixel value for both two points.    * 2 . zoomlevel 19 or something else   * 3 . Using google map and bing map What I need is just pixel per meters(ratio).


title: Ontario Topographic and Property Data

text: Looking for Ontario topographic data and property line data; willing to purchase the latter. Shapefiles or file geodatabase files would be optimal. In terms of topographic data, I am looking for: basic outline of Ontario, county or municipal outlines, contour data, water (lines, poly, etc), buildings, and so forth.


title: How to create a new column to store the length of a linestring in postgis?

text: I am working with pgRouting, following this nice guide：A Beginner’s Guide to pgRouting I am now in this step：Calculating shortest routes, with these codes below               SELECT * FROM shortest_path('     SELECT gid AS id,            start_id::int4 AS source,            end_id::int4 AS target,            shape_leng::float8 AS cost     FROM network',     1,     5110,     false,     false);      and my network table looks like this ![enter image description here](http://i.stack.imgur.com/zCjRg.png) apparently I have no column "shape_leng", so I need to create a new column in table network, and it is the length of every row. But I don't know how I should type the ST_length code. Can someone give me some advice？


title: Help with WPS.client

text: Hy! I am trying to make an application which will test if any of the segments (multi line layer acquired from GeoServer) intersects with polygon acquired using OpenLayers.Protocol.WFS.fromWMSLayer, and if it does add it to the map. After the polygon has been selected it's added to an empty vector layer which is then used as an input in WPS.Client. The code looks like this:               int = client.getProcess('server', 'JTS:intersection');             for (var i=0; i<v_layer.features.length+1; i++){             int.execute({                      inputs: {                               a: selected.features[0].geometry,                               b: segment_layer.features[i].geometry                              },                      success: function(outputs) {                                                  map.layers[3].addFeatures(outputs.result);                                                 }        });                  }      This works, when I tested it on just a few segments, but since this layer contains about 200 of them it crashes the browser. I couldn't find any other WPS process that i could use other than JTS:intersection. Any suggestions? Thank you for your time.


title: NVIZ runs from GRASS but not via QGIS

text: When I try to run NVIZ from QGIS I get an error which states 'Can't find usable init.tcl in the following directories" and then goes on to list a number of directories. (I did have a screenshot here but couldn't post it as a newbie.) However when I try to run NVIZ directly from GRASS it runs properly. I have searched for the init file which the QGIS error report says cannot be found and it is actually in "C:\Program Files (x86)\Quantum GIS Lisboa\lib\tcl8.5". According to the QGIS error report this is not one of the directories actually searched. I do have a workaround running NVIZ from GRASS but that is clearly not as satisfactory as running from within QGIS. I have noted the reply to an earlier question on this topic which recommends installing tcl/tk via OSGEO4W. However I do not believe that these components are actually missing since I've found the init file. Can anyone suggest a solution, please? I am a new user of QGIS and have downloaded 1.8.0 Lisboa, running under Windows 7.


title: Java shared hosting with PostGIS

text: I've been developing a disaster risk management application that makes use of Java/Smart GWT, PostgreSQL/PostGIS, and Tomcat 7. I'm looking for a shared hosting environment for the beta testing period (about 6 months). My current hosting provider does not provide PostGIS, or PostgreSQL greater than version 8.1. Any pointers or advice will be appreciated.


title: Geoserver 2.5.1 Image Mosaic Plugin

text: I have geoserver 2.5.1 running on ubuntu 14.04 server. I tried to compose an image mosaic with ImageMosaicPlugin but I always get a message **Could not list layers for this store, an error occurred retrieving them: Failed to create reader from file:data/dalekovod and hints nul**. Images are in the right file format because geoserver recognizes them, and I can load each at one of them separately but I have 120 images. Any idea what could be the problem?


title: Transect with square polygons

text: would like to create with QGis 2.2 (german version) several square polygons (1km ²) that seamlessly connect to each other, so that a strip ( transect ) of 9km times 1km length results . ( I once attached a picture of how it should look ) -  http://www.myimg.de/?img=TransektArcMap75b51.jpg Although I have installed the plugin " Rectangles ovals digitizing ", but can not choose here the edge length . In addition, the following square is to connect seamlessly . In ArcMap the cursor shows me when digitizing contacts with neighboring polygons or their vertexes or lines on , so you could work " seamlessly " . Furthermore, each of the squares , partially each quadrant , then an attribute (here, an air quality index) to be assigned , in accordance with which the square to be dyed. Are different colors within a shape based on various values (how it works in ArcMap ) possible? About tips on how I can tackle this project , I'd be very happy . PS: For information on tools , etc. ..... I use the German user interface. Best regards and thank you in advance Daniel


title: Can't start ArcSDE (esri_sde) service

text: I am getting an error message when atempting to create a connection in ArcCatalog (I'm working under the 10.0 version). The message is : Failed to Connect to the specified server. No ArcSDE server licence found. Then I Check the Service esri_sde and it has not started, so I try to start it then I get this error message : The ArcSde Service(esri_sde) service on Local Computer started and then stopped. Some services stop automatically if they have no work to do. I've tried deleting the service then creating it all over again, checked the registry, launched CCleaner.. Rebooted the computer, still the same error When I check the log files, I got this in the giomgr_esri_sde.log : Sat Jun 14 16:15:23 2014 - ****** Spatial Database I/O Manager ****** Sat Jun 14 16:15:23 2014 - SDE System Startup Initiated . . . Sat Jun 14 16:15:23 2014 - Validated System Paths. Sat Jun 14 16:15:23 2014 - SDE IOMGR going into background . . . Sat Jun 14 16:15:23 2014 - System Parameter table 'SERVER_CONFIG' Read. Sat Jun 14 16:15:24 2014 - Initialized Shared Memory Segment. Sat Jun 14 16:15:24 2014 - Error (-327), No ArcSDE server license found. And in the sde_esri_sde.log : * * * ## ArcSDE 10.0 for PostgreSQL Build 685 Fri May 14 12:05:43 2010 Connected to instance . . . DBMS Connection established... Error (-327), No ArcSDE server license found. I don't understand this licence problem! I have already registered and I do have the licence. What could be wrong ?


title: Why do TIF loaded into PostGIS lose their colormap?

text: I have successfully loaded a TIF image (using raster2pgsql) into my PostGIS 2.0.1 database. I can load the raster image into QGIS and it displays in the correct space. However, it has lost all colour and appears as a greyscale image. (Layer Properties > Style > Single band gray, colormap = grayscale) If I load the same image directly from file into QGIS it displays colours correctly.(Layer Properties > Style > Single band gray, colormap = colormap) If I convert the TIF to JPG and load the JPG to PostGIS and then display in QGIS the colours render correctly. (Layer Properties > Style > Three band color) I can save the colourmap from an image and apply it to the greyscale images but the colours aren't quite right. What do I need to do to get the TIF images to retain their colourmap when I load to PostGIS? Windows 7 32-bit PostgreSQL 9.1.4 PostGIS 2.0.1 QGIS 1.8.0 raster2pgsql -s 27700 -I -C -M *.tif -F -t 200x200 gis_rasters.osvmd > vmd.sql Thanks in advance Ross


title: Merging raster data with different resolutions

text: I need to merge different raster datasets which have one layer and contain discrete values to a layer stack. The resolutions of the data can have any value so not only multiples of the resolutions from the other input rasters thus there are problems with alignment. Also the corner coordinates are not necessarily identical. Has anybody got a creative idea of what is the best way to do the merging without losing information as a layer stack can only have layers with the same resolution? I have used gdal_merge.py and of course there is quite some loss of information through the nearest neighbour resampling (though I must use nearest neighbour as I want to keep discrete values). One thing I thought of was to polygonize each raster and then union them and rasterize the resuling polygon with a appropriate resolution. That would keep most information but here I am afraid that it will not work properly with large datasets. Perhaps somebody has experience with this? Thanks for any help! Kind Regards DomR


title: pgrouting crashes and restarts the postgresql server when searching for endpoints in 2 disjoints graphs

text: I have               >apt-show-versions postgresql-9.1-pgrouting     postgresql-9.1-pgrouting/precise uptodate 2.0.0-release-ppa1~precise2      Whenever I do a search http://docs.pgrouting.org/dev/src/astar/doc/index.html between two end points in disjoint graphs, I get the error error sending I/O to backend and my postgresql server crashes and restarts. Why doesnt this simply return a no route found?


title: Alternative indexing methods for point set operations

text: It is common to use a bounding box spatial index to improve performance when working with a large number of features. Where operations are performed against individual geometries with a large number of vertices do any similar optimization strategies exist? For example, do any data structures exist that can speed up point in polygon, or union operations?


title: How do I extract text from dwg?

text: I'm trying to convert data from an AUtoCad Civil 3d drawing to a text format to use in ArcGIS. Each dwg has coordinate and text information that I would like to use for point location and attributes. The text looks to be a table but is in fact just mtext and line objects. I have found LISP code that will export all the mtext to a text file but I am wondering if there is a way to place all the text found at the same xy location in the dwg to be placed on one line in the text file? Thanks!


title: PostGIS: ST_Transform function conversion problem

text: I'm inserting map data from a GML file with srid EPSG:27700 into a PostGIS database table with srid EPSG:4326. I'm using the ST_Transform function as in this snippet:                ST_Transform(ST_GeomFromGML('<gml:Point srsName="EPSG:27700">         <gml:coordinates>527347.0,181346.0</gml:coordinates></gml:Point>'), 4326);      After the data was inserted in the database I tried to find if its geometry intersects with an existing table's geometry using the ST_Intersects function. The function returned FALSE in all cases, but I knew that the two geometry layers _should_ intersect. So, I loaded the layers in QuantumGIS to see what's happening. The PostGIS layers indeed don't intersect, but when I loaded the GML data directly from the GML file using the original srid (EPSG:27700) the map data correctly appeared on top of the pre-existing layer with srid 4326, with which I knew it should intersect. The problem is that the ST_Intersects function (as well as all the others of the same kind) works only for geometries with the same srid, so I need to solve the transformation issue. Apparently there's either a problem with the ST_Transform function, or I'm missing something that I should know. I use PostGIS 1.5.3. Any help will be greatly appreciated.


title: Reinventing the wheel: Public Transit System

text: I'm building a public transit service in South Africa, and I initially started the project on PHP and MySQL. I've managed to get far with the simple stack, with:   * Routing algorithm in PHP   * Data storage as coordinates in MySQL   * front-end in WordPress   * JS for Google Maps and a REST API for mobile apps (incomplete) _My question based on the above is: I want to move away from using PHP in the back-end, and I'd like to know what alternatives there are._ I'm fairly inexperienced, and used PHP as it's the only server-side language that I currently know. I've seen that a lot of map/routing applications use Java and Python, and am thinking of moving in that direction in future. I'm finding myself having to do a lot of geographic calculations, and as my data grows I find myself reinventing the wheel. **Would a PostGIS + PgRouting setup benefit me?** A working demo of the project can be found at http://rwt.to. **UPDATE:** The project has subsequently went live. Thanks for all the suggestions.


title: How to draw (or add) a line indicating a specific latitude/longitude?

text: I'm trying to add a single line that show a specific latitude (e.g. the Tropic of Cancer) I don't want a gird but only a single line even in a vector layer.How can I do that I've try with pyQGIS console but seems I can't... I would like to measure angle with that line so I'd like to have it in a layer. However I found useful a method to use it in print composer, but my main goal is to have it as line in QGIS. For the moment I'm using Google Satellite open layers with a EPSG:3857.


title: Issue on IFeatureDataConverter.ConvertFeatureClass

text: Does anyone know if it's possible to reuse this method with a "feature dataset" already present ? I'm trying to import multiple shapefiles under a featuredataset inside a fgeodb. Everything it's working fine the first time when I create a new featuredataset but if I try to load a second shapefile inside it I get an error. I know that a work around it's to use the geoprocessor but there are some interesting thinks in dataconverter that I need. Tnx


title: ArcGIS Server: Better to have one massive map service or multiple smaller services?

text: I have a Silverlight application that consumes around 30 map service layers at different types. The service layers are categorized based on type and only applicable to certain screens of the application. For organization's sake, I split each of the layer categories into separate map services, but I'm considering whether or not that makes sense. What is the best practice in this situation? Should I just lump all of the map layers into a massive map service or keep them separate? I'm curious about server load vs application performance.


title: CORINE LandCover: what about the accuracy?

text: How accurate is CLC? Do anyone know any scientific papers about that question? I tested the new CLC2006 for a region in turkey and its totally rubbish!


title: Using web mapping and geocoding to track genealogical data

text: Over the last year I have gotten involved in researching my family tree. My paternal grandfather came from a remote island in the Atlantic, along with hundreds of others of his generation, and many of them settled within the same area. As the island is very small, they shared a small number of family names (about 30 surnames, with some variations of the same name). I would really like to use my GIS knowledge to firstly gather data on where the descendents are living now (the recently established Society has compiled a small database), and then to analyse the data to map out the various family names/bloodlines/races (the islanders are a melting pot of races, absolutely fantastic). Before I offer my services to the Society, I'd like to know if anyone has undertaken something similar, and what would be the best way to approach it? As mentioned in the question title, I want to use ArcGIS Online so that others will be able to add their info and I would manage the data. Some light googling indicates that geocoding surnames might be part of what I am looking for, but as I have no experience with it I'm not sure.


title: Loading a multipoint shapefile in R

text: This issue has only been raised in the community once, two years ago, to the best of my knowledge. I think it's high time to raise the issue again! `readOGR`, from the rgdal package, cannot handle multipoint shapefiles, even though they load fine in other GIS packages like QGIS. Of course, you can save the file as a single part object (see image below) in an external program, but ideally it would also be possible to do it in R. Please see worked example below. ![multipart shp file in qgis](http://i.stack.imgur.com/UPRoR.png)               download.file("http://www.personal.leeds.ac.uk/~georl/egs/lnd-stations-multi.zip",                    "lnd-stations.zip") # download multi part polygon     unzip("lnd-stations.zip") # unzip     lndS <- readOGR(".", "lnd-stations", p4s = "+init=epsg:27700") # load          OGR data source with driver: ESRI Shapefile      Source: ".", layer: "lnd-stations"     with 2532 features and 27 fields     Feature type: wkbMultiPoint with 2 dimensions     Error in readOGR(".", "lnd-stations", p4s = "+init=epsg:27700") :        Incompatible geometry: 4      After performing "Multipart to singleparts" in the QGIS Vector menu (illustrated above) and saving the file, it loads fine in R. Try it:               download.file("http://www.personal.leeds.ac.uk/~georl/egs/lnd-stns.zip",                    "lnd-stns.zip")     unzip("lnd-stns.zip")     lndS <- readOGR(".", "lnd-stns", p4s = "+init=epsg:27700")     plot(lndS)


title: PostGIS/QGIS and granting permissions

text: I've got a PostGIS databases, with lots of schemas. I can easily create a new user that will allow them to update and delete features. However when I get the user to create a new feature I get the following error in QGIS:-                 PostGIS error while adding features: ERROR:  new row for relation "conservation_areas" violates check constraint "enforce_geotype_geom      I've used the following SQL to GRANT the user (anne) privilages that will allow her to update and delete, but I guess its the insert bit that is letting me down.               GRANT ALL ON SCHEMA conservation TO anne;     GRANT SELECT ON ALL TABLES IN SCHEMA conservation TO anne;     GRANT UPDATE ON conservation.conservation_areas_gid_seq To anne;     GRANT UPDATE ON conservation.conservation_areas TO anne;      The postgres user has no such issues when edited data, therefore I'm missing something somewhere when i'm granting the user permissions. Is there a simple SQL that will allow a user to do everything that the Postgres user can do but to a specific tables (and not all tables) within a schema? Med