title: How to convert MultiLinestring to Linestring?

text: I have a multiline string and want to convert it to a linestring in order to use some functions that only work with linestrings. I used `ST_Dump()` and got a collection of my geometries. If I merge them back together, I still get the multilinestring. Should I get the Npoints and combine them together or something, in order to end up with a simple linestring? `ST_LineMerge()` also does not work on my multilinestring, it returns the same. my example, MULTILINESTRING((-3.16420835153456 55.9269166007097,-3.164222 55.926918),(-3.1642070167833 55.9269296196706,-3.16421351659546 55.9268662214904),(-3.16421351659546 55.9268662214904,-3.16421636372824 55.9268384509897),(-3.16421636372824 55.9268384509897,-3.16422182573761 55.9267851753802),(-3.16422182573761 55.9267851753802,-3.16422870102352 55.926718114886),(-3.16422870102352 55.926718114886,-3.16423309121073 55.926675293667),(-3.16423309121073 55.926675293667,-3.16423565148822 55.9266503211093),(-3.16423565148822 55.9266503211093,-3.16424103159897 55.9265978443265),(-3.16424103159897 55.9265978443265,-3.16424680776317 55.9265415044985),(-3.16424680776317 55.9265415044985,-3.16425267254583 55.9264843002995),(-3.16425267254583 55.9264843002995,-3.16425541048045 55.9264575949012),(-3.16425541048045 55.9264575949012,-3.16426111146586 55.9264019883556),(-3.16426111146586 55.9264019883556,-3.1642667032531 55.9263474469124),(-3.1642667032531 55.9263474469124,-3.16426957768543 55.9263194101362),(-3.16426957768543 55.9263194101362,-3.16427488261739 55.9262676666359),(-3.16427488261739 55.9262676666359,-3.16428009893088 55.9262167875066),(-3.16428009893088 55.9262167875066,-3.164282741107 55.9261910161221),(-3.1642875546472 55.9261440655823,-3.164282741107 55.9261910161221),(-3.1642875546472 55.9261440655823,-3.16429466890915 55.9260746741522),(-3.16429466890915 55.9260746741522,-3.16430092974527 55.9260136069079),(-3.16430092974527 55.9260136069079,-3.16430822838418 55.9259424170929),(-3.16430822838418 55.9259424170929,-3.16431547242401 55.925871759829),(-3.16431547242401 55.925871759829,-3.16431448732505 55.9258328901507),(-3.16431770120536 55.9257327846001,-3.16431547242401 55.925871759829),(-3.164339 55.925777,-3.16431770120536 55.9257327846001))


title: Enable Add-In Button if One Feature is Selected from Specific Layer

text: I am trying to enable an add-in button on my toolbar if only one feature is selected from a specifc Layer. For example, when ArMap loads, the button should be disabled. If the user selects one feature from the AJ_RecommendedImprovementsPoints layer, then the button enables. Here's what I have, but it's not working at all. suggestions??               Protected Overrides Sub OnUpdate()         'Get Transportation Improvement Projects         Dim pFeatLyr As IFeatureLayer         pFeatLyr = GetLayerByTOC("RV_ImprovementPoints")         'Enable Button         If GetFeatSelCount(pFeatLyr) = 1 Then             Me.Enabled = True         Else             Me.Enabled = False         End If     End Sub               Public Function GetFeatSelCount(ByVal pFeatLyr As IFeatureLayer) As Long         'Returns # of selected featurues of feature layer as long integer         'Initialize the required variables         Dim pMxDoc As IMxDocument         pMxDoc = My.ArcMap.Document         Dim pMap As IMap         pMap = pMxDoc.FocusMap         Dim pFeatSel As IFeatureSelection         pFeatSel = pFeatLyr         Dim pFeatSelectionSet As ISelectionSet         pFeatSelectionSet = pFeatSel.SelectionSet         GetFeatSelCount = pFeatSelectionSet.Count     End Function               Public Function GetLayerByTOC(ByVal lyrName As String) As IFeatureLayer         'This function finds a feature layer based on its TOC Name         Try             'Initialize the required variables             Dim pMxDoc As IMxDocument             pMxDoc = My.ArcMap.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         Catch ex As Exception             'Layer not found. Show a message             MsgBox("Layer " & lyrName & " not found !", vbCritical, "Error")         End Try     End Function


title: How to convert shp to csv

text: I am trying to convert my shape files to csv through a batch processing in Java using Geotools. The problem is I have more than 500 shape files (point features). The code is working fine for a small number of shp files. But when I am iterating the process over a large dataset (number of files), some of the files are not written and of 0 KB. I have no clue why this is happening. Any help would be highly appreciated. I want all the files to be converted from shp to csv. Here is my java code-               // shp2csv conversion     import java.io.BufferedWriter;     import java.io.File;     import java.io.FileWriter;     import java.io.IOException;     import org.geotools.data.FileDataStore;     import org.geotools.data.FileDataStoreFinder;     import org.geotools.data.shapefile.ShapefileDataStore;     import org.geotools.data.simple.SimpleFeatureCollection;     import org.geotools.data.simple.SimpleFeatureIterator;     import org.geotools.data.simple.SimpleFeatureSource;     import org.opengis.feature.simple.SimpleFeature;          public class GeoOpen     {               private static SimpleFeatureIterator simpleFeatureIterator;     static FileDataStore store ;     static String ID="";                    public  static boolean openShapeFile(File srcfname) throws Exception                    {         String filename=srcfname.getAbsolutePath();          File dataFile = new File(filename);         dataFile.setReadOnly();         store = FileDataStoreFinder.getDataStore(dataFile);        // ShapefileDataStore store = new ShapefileDataStore(dataFile.toURL());              SimpleFeatureSource source = store.getFeatureSource();         SimpleFeatureCollection featureCollection = source.getFeatures();         simpleFeatureIterator = featureCollection.features();         return true;     }     public  static boolean iterate(File srcfname, String dest_path)     {         File dest_file=srcfname;              String dest_filename=GeoOpen.getFileNameWithoutExtension(dest_file)+".csv";              try{             BufferedWriter  bw=new BufferedWriter (new FileWriter(dest_path+dest_filename));             bw.write("ID"+","+ "date"+","+"tstamp"+","+"X_prj"+","+"Y_prj"+","+"NEAR_FID"+","+"NEAR_DIST");         while(simpleFeatureIterator.hasNext())         {                  SimpleFeature f = simpleFeatureIterator.next();               //  System.out.println(""+f.getID()+" ,"+f.getAttribute(1)+", "+f.getAttribute(2)); //ID, field 1: lat, field 2: long             //writing in a CSV file            bw.write("\n");                 int index = f.getID().lastIndexOf('.');            ID=f.getID().substring(index+1); //get only the ID //by default ID comes with filename.ID format when getID() method called           bw.write(Integer.parseInt(ID)+" ,"+f.getAttribute(1).toString()+", "+f.getAttribute(2)+","+f.getAttribute(3)+","+f.getAttribute(4)+","+f.getAttribute(5)+","+f.getAttribute(6));              }              }         catch (IOException e)         {             e.printStackTrace();         }         catch(Exception ex){             ex.printStackTrace();              }         finally{             simpleFeatureIterator.close();             store.dispose();           }              return true;     }          public static String getFileNameWithoutExtension(File f)     {         String s="";         int index = f.getName().lastIndexOf('.');         //System.out.println (index);       if (index>0&& index <= f.getName().length() - 2 )        {           s= f.getName().substring(0, index);       }       return s;     }     public static void main(String[] args) {         try {                String dest_path="C:\\ArcGIS test data\\commonfiles csv\\001\\nr_train\\";                File src_dir= new File ("C:\\ArcGIS test data\\commonfiles prj\\001\\nearTrain");                File [] file_nr = src_dir.listFiles();  //array containing files from nearRoad             for (int i=0; i<file_nr.length; i++)  //near road             {                 if (file_nr[i].getName().endsWith(".shp"))                 {                     GeoOpen.openShapeFile(file_nr[i]);                     GeoOpen.iterate(file_nr[i], dest_path );                 }             }              } catch (Exception e) {                  e.printStackTrace();         }     }     }


title: Exporting all attributes from a GPX file

text: I have a GPX file that i have imported into Google Earth, it looks like this. ![enter image description here](http://i.stack.imgur.com/f6gAn.jpg) When I try and use the GPX to features tool in ArcGIS, it does not export out all the attributes.. ![enter image description here](http://i.stack.imgur.com/x0ksj.jpg) It seems that GE has access to more data. What do I need to do to export all the attributes shown in the Google Earth popup into my feature class. Thanks


title: Move Duplicate Nodes in the Same Layer

text: In ArcMap 10.1 can I move duplicate nodes in the same layer like I can in MapInfo? I want to click on the place where two or more lines converge and move the endpoints of all the lines to the location where I release the mouse. In MapInfo it is Options > Preferences > Map Window > Editing > Move Duplicate Nodes in the Same Layer > OK.


title: Python Script examples for geoprocessing shapefiles without using arcpy

text: I would like to use a Python script that is not based on arcpy to do things like query a shapefile by attributes, create new layer from selection, and calculate areas of a polygon and convert polygons to points. Anyone have any code examples of using other Python modules or libraries? I am able to do this easily using arcpy but i wanted to explore other options.


title: Can QGIS create a "virtual view" from a given position?

text: I would like to be able to choose a position on a map (by clicking, or by entering an address, or other means -- please suggest) and have an artificial picture created, like what I would see if I stood in real world in that spot. What would a photograph from that spot show? I think I would need layers of elevation, and maybe buildings ("infrastructure"?) and maybe land cover (vegetation type). What do you call what I am trying to do, and how do you do it with QGIS? Please if you answer remember I am quite a beginner -- I do not understand 1/4 of the answers I have read on this site! Thank you.


title: MultiPatch to Point Shapefile

text: I have a File GeoDatabase that has a Feature Dataset with MultiPatches. The MultiPatches are tree models. I would like to convert the MultiPatch into a point layer. Is there a tool, plugin, or script that will convert the MultiPatch layer into a point feature class. Thanks...


title: Feature Service from plain Postgres/PostGIS in ArcGIS Server 10.2

text: Using 10.2 Arcmap/ArcGIS Server on Windows I've successfully published some map services from my Postgres/PostGIS data and now I want to enable editing the features through a JavaScript web app. This is NOT through SDE tables (i.e., creating a geodatabase). Yes, I am using Postgres 9.3.1 and Postgis 2.1... but everything else seems to work fine but this PK error. How do I get ArcGIS Server to trust me that the incrementing PK I've created is OK? > "00166: Layer's ObjectID is not maintained by the database and the feature > service is being published with Create capability enabled" So, following the JS API examples I am enabling both the geometry service and the Feature Access service with all the Operations and Capabilities enabled. I began using the serial data type for the OBJECTID primary key which appears to be a shortcut for creating an integer, auto-incrementing key as follows:               CREATE TABLE test_schema.volmonwgs84e     (       objectid serial NOT NULL,       volmon24_i integer,       site_id integer,       waterbody character varying(30),       geom geometry(MultiPoint,4326),       CONSTRAINT volmonwgs84e_pkey PRIMARY KEY (objectid)     )     WITH (           OIDS=FALSE     );     ALTER TABLE test_schema.volmonwgs84e       OWNER TO postgres;               CREATE INDEX volmonwgs84e_geom_gist       ON test_schema.volmonwgs84e       USING gist       (geom);      However, when I go to publish the service (from ArcMap) the Analyzer shows the error above. Right-clicking the error kindly offers to "Add an auto-incrementing ID" for me, to which it later responds with the following. AFAIK, I ALREADY have an auto-incrementing field anyhow.               "Executing: AddIncrementingIDField testing.test_schema.volmonwgs84e BOO     Start Time: Tue Dec 10 15:55:06 2013     ERROR 000152: You cannot add an incrementing ID field to a table registered with the geodatabase.     Failed to execute (AddIncrementingIDField).     Failed at Tue Dec 10 15:55:07 2013 (Elapsed Time: 0.70 seconds)[/CODE]      So, what data type and rules need to be on this auto-incrementing field, and what is AddIncrementingIDField trying to do for me? Thanks! UPDATE: Thinking that my issue was the "new" PostGIS/Postgres approaches to things I've also tried creating a table with these older approach as follows:   1. manually setting the sequence on integer PK column instead of using the SERIAL type,    2. using "SELECT AddGeometryColumn" to create the geometry column, and therein using the old constraints approach with the final argument = false (http://postgis.net/docs/AddGeometryColumn.html)  > >     CREATE SEQUENCE test_schema.mypoint_const_b_seq >     INCREMENT 1 >     MINVALUE 1 >     MAXVALUE 9223372036854775807 >      START 1 >     CACHE 1; >                  ALTER TABLE test_schema.mypoint_const_b_seq      OWNER TO postgres;          CREATE TABLE test_schema.mypoint_const_b     (       objectid integer NOT NULL DEFAULT nextval('test_schema.mypoint_const_b_seq'::regclass),       volmon24_i integer,       site_id integer,       waterbody character varying(30),       CONSTRAINT mypoint_const_b_pkey PRIMARY KEY (objectid)     )     WITH (       OIDS=FALSE     );          SELECT AddGeometryColumn ('test_schema','mypoint_const_b','geom',4326,'MULTIPOINT',2,     false);          ALTER TABLE test_schema.mypoint_const_b       OWNER TO postgres;          CREATE INDEX mypoint_const_b_geom_gist       ON test_schema.mypoint_const_b       USING gist       (geom);      But I still get a fatal error when trying to create a service with Feature Access "Create" capability enabled with the same message: [CODE]"00166: Layer's ObjectID is not maintained by the database and the feature service is being published with Create capability enabled" [/CODE] I know that the database IS maintaining the PK because I can add features using QGIS just fine. AND, ArcMap's properties for the indexes on the layer shows objectid is a unique index. PS, buried in this excellent link (http://proceedings.esri.com/library/userconf/proc13/tech- workshops/tw_192.pdf) I see that Arc doesn't like fields with quotes (which is how you force CAPS on OBJECTID. Lowercase objectid _appears_ to be just fine to ArcMap, but perhaps not ArcServer.


title: How do I perform a raster overlay analysis in ArcObjects?

text: Using the "IMultivariateOp.MLClassify" Method, it will create a field call "Count",which can be used for calculating area, as shown in the figure below: ![MLClassify](http://i.stack.imgur.com/FZwK9.png) But now I want to calculate the increase area of two rasters, for example a raster in "2000" year and "2013" year. So I have to analyze the overlay, as shown in the figure below: ![The algorithm or the increasement](http://i.stack.imgur.com/Xc06E.png) How do I execute this analysis in ArcObjects using C#? Thanks very much.


title: How to create shapefiles with z coordinate with QGIS

text: I'm working with QGIS, and I want to create a shapefile with the elevation coordinate, I can't find something similar to what I was used to do in Arcgis, can you help me? Thanks for your help Luís


title: Python problems with select by attributes and user input

text: I am attempting to create a tool where the user can input a census tract number within the city of Plano (there are 77 tracts) and it ultimately creates a map layout showing the population density. I am having problems using the select tool to copy the inputed tract layer into a new feature class. It is executing successfully but it keeps copying the original plano_tract layer into a new feature class instead of copying the selected tract that was inputed. Not sure if I needed a where clause in the script or if I could simply just do a select by attribute function. I was using this thread as reference but it I still getting the same result. Python Script for Select by Attributes taking user input I am pretty new to Python and this is a beginner exercise for me. Here is my script:               import arcpy     import os     from arcpy import env          #set workspace     Workspace = arcpy.GetParameterAsText(0)          env.workspace = r"C:\Users\dallascowboy83\Desktop\project 2\project\project.gdb"     env.overwriteOutput = 1          #set local variables     tracts = arcpy.GetParameterAsText(1)          #ask user for tract to search     tract = raw_input("Please type Tract Number: ")          plano_tract = r"C:\Users\dallascowboy83\Desktop\project 2\project\project.gdb\plano_tract.shp"          plano_tract_lyr = "plano_tract_lyr"          #user selects tract, Select by Attribute tool runs     arcpy.MakeFeatureLayer_management("plano_tract", "plano_tract_lyr")          arcpy.SelectLayerByAttribute_management("plano_tract_lyr", "NEW_SELECTION", "\"Tract_Number\" = "+ tract)          arcpy.CopyFeatures_management("plano_tract", "plano_tract_feat")          #add field     arcpy.AddField_management("plano_tract_feat", "density", "DOUBLE")          #calculate field     arcpy.CalculateField_management("plano_tract_select", "density", "[DP0010001] / [sq_mile_1]", "PYTHON", "")          print "executed successfully"


title: Command Line FGDC/ISO 19115 Metadata Editor

text: Are there any metadata editors compliant with ISO 19115 and/or FGDC metadata standards which can work with XML formatted metadata and run solely on the command line? Often times I'm working remotely or otherwise desire a headless environment and would like to do editing of metadata from the shell. Most of the metadata editors I've seen either run in a GUI (i.e. CatMDEdit) or are web based like GeoNetwork. I've looked at this list of editors for now, and haven't really found anything that fits the bill. http://www.fgdc.gov/metadata/geospatial-metadata-tools Although I know the 'mp' program can validate metadata, I haven't come across an actual command line metadata _editor_ yet. Any suggestions would be appreciated.


title: Port ArcGIS Engine to ArcGIS Runtime SDK for WPF?

text: We have a WPF application, which has arcengine controls embedded in it. We are looking at options to replace ArcGIS Engine with someother component which is 64 bit and we are evaluating ArcGIS Runtime SDK for WPF (possibly .net as well). We do not need an ArcGIS for Server or ArcGIS Online based solution, we dont use web maps or tiles services. We have file geodatabases and need to create, edit and modify feature classes within the gdb. 1**. Can I use ArcGIS Runtime SDKs to build a non-cloud / non - server WPF Solution.i.e. The map control in the WPF sdk will need to access only a local fgdb (not hosted as a service).**


title: Adding values to Fishnet Grid?

text: I have developed a fishnet grid for my study area. I have various datasets for the area and want to allocate a value (from one of the pre-determined datasets) to each polygon within the Fishnet. How do i display that the data from one source falls into the various polygons? I am using ArcGIS, after completing the Spatial Join the values i wished to display within the Fishnet Polygons all came up as Null? Any solutions to this?


title: Randomly select a record to update using arcpy.da.UpdateCursor

text: This question is related to my question here. What I would like to do is use a cursor with a where clause to select a number of records. For each row, the values in two fields, `TYPE_A` and `TYPE_B` should be compared. If those two values are not the same, the field `NUM_ID` should be updated with the current value of the counter. Once that is done, _all but one_ of the rest of the records in the selection should be updated with the value of the incremented counter. Some example data:               FEATURE 1 :NUM_ID = 4   TYPE_A = NEW   TYPE_B = NEW     FEATURE 2 :NUM_ID = 4   TYPE_A = NEW   TYPE_B = EXISTING     FEATURE 3 :NUM_ID = 4   TYPE_A = NEW   TYPE_B= NEW     FEATURE 4 :NUM_ID = 4   TYPE_A = NEW   TYPE_B= NEW     FEATURE 5 :NUM_ID = 11  TYPE_A = NEW   TYPE_B = EXISTING     FEATURE 6 :NUM_ID = 14  TYPE_A = NEW   TYPE_B = NEW     FEATURE 7 :NUM_ID = 11  TYPE_A = NEW   TYPE_B = EXISTING      The 4 features with `NUM_ID = 4`would be selected and Feature 2 will be given the NUM_ID 21, as its values for TYPE_A and TYPE_B are different. The next step would be to update _all but one_ of the other 3 records by incrementing the counter. The next selection would be the 2 features with `NUM_ID = 11`. As both the records have different values for both fields, they would both be updated. The last selection will be for NUM_ID = 4. As both its values for the fields are the same, and it is the only record in its selection set, no update is needed. The result:               FEATURE 1 :NUM_ID = 23   TYPE_A = NEW   TYPE_B = NEW     FEATURE 2 :NUM_ID = 21   TYPE_A = NEW   TYPE_B = EXISTING     FEATURE 3 :NUM_ID = 22   TYPE_A = NEW   TYPE_B= NEW     FEATURE 4 :NUM_ID = 4   TYPE_A = NEW   TYPE_B= NEW     FEATURE 5 :NUM_ID = 11  TYPE_A = NEW   TYPE_B = EXISTING     FEATURE 6 :NUM_ID = 14  TYPE_A = NEW   TYPE_B = NEW     FEATURE 7 :NUM_ID = 25  TYPE_A = NEW   TYPE_B = EXISTING      The code I have:               getNum("%Selected Features%")      * * *               def getNum(ftr):       lst = [4, 11, 14, 15] #values for SQL expression       x = 20       fields = ('TYPE_A', 'TYPE_B', 'NUM_ID')       assert "NUM_ID" in [f.name for f in arcpy.ListFields(ftr)], "NUM_ID field no in feature class"       for l in lst:         y = x + 1        whereclause = "%(sql_name)s = '%(l)s' """ {"sql_name":arcpy.AddFieldDelimiters(ftr, "NUM_ID"), "l":l}         with arcpy.da.UpdateCursor(ftr, fields, whereclause) as cursor:           for row in cursor:             if row[0] ! = row[1]:               row[2] = y             #some condition here             cursor.updateRow(row)         y += 1      `TYPE_A` only has one value in it for all features i.e `NEW`. `TYPE_B` could have any value. It does not matter which feature gets updated amongst those which have the same values for both fields (`NEW`), as long as one of them keeps the original value for NUM_ID. While this code updates the row correctly where the two fields don't match, I do not know how to let one of the records keep the original value, while the rest should be updated.


title: How to preserve original resolution of scanned 1-bit PDF raster with gdal_translate?

text: I am using `gdal_translate` to convert non-spatial PDFs to non-spatial TIFFs (intended for subsequent manual georeferencing). This is my `gdal_translate` command:               gdal_translate -of GTiff -sds -co NBITS=1 --config GDAL_PDF_DPI 300 %input_pdf% %output_tif%      Here is what it looks like at ~8x magnification (note the "nibbling" around the edges of text and the complete loss of many "speckles"): Input (rendered by Adobe Reader, disregard the slight blurring of edges): ![Input](http://i.stack.imgur.com/CP7P1.png) Output (rendered by ArcMap): ![Output](http://i.stack.imgur.com/gfqvx.png) If I do not specify the output bit depth (by omitting the `NBITS` option), I can tell that GDAL is re-rasterizing and resampling the input 1-bit image into a blurred grayscale image: ![Resampled](http://i.stack.imgur.com/Io63s.png) Ideally it would just copy each pixel in the input straight to the output. How can I retain the original resolution without information loss in the output? Should I use something else, like ImageMagick, Poppler or GhostScript? All this re-rasterization GDAL is doing is painfully slow, I might add, so I think I will try one of the alternatives for that reason alone.


title: How to style a point based on Actual Attribute Values?

text: I have published a point layer in Geoserver, which is the data from a GPS Survey. In addition to the geometry, my layer, also contains an Attribute for Accuracy, which indicates the Accuracy in meters of the reading. I have been trying to style the layer, such that a point is show as a point surrounded by a circle whose radius is equivalent to the Accuracy in map units. (Similar to how most mobile apps show the GPS points) All of the Attribute based symbolization that I have seen is Rule based (i.e. Classification). How do you Symbolize a point, such that its radius is equal to its accuracy?


title: Can FieldMappings be used to drop all attributes?

text: We have a script that regularly copies a bunch of feature classes out to shapefiles, keeping only the most important attributes (e.g. customer category) and dropping useless/private information (e.g. customer name). Right now, it repeats the same set of steps without much variation. I am attempting to significantly simplify it by defining a function to do all those steps (the entire purpose of functions, of course.) Also, by introducing the arcpy.FieldMappings objects instead of verbose field mapping strings, it's possible to quickly understand what fields are coming from which data sources. This works great for most of my feature classes.               def generate_field_mapping(inTable, inFieldList):         fm = arcpy.FieldMappings()         for field in inFieldList:             vars()[field] = arcpy.FieldMap()             vars()[field].addInputField(inTable, field)             fm.addFieldMap(vars()[field])         return fm          def subset_attributes(origFC, fieldList, outFC):         fStart = time.time()         name = origFC.split('.')[-1]         # make feature class with only desired fields         stringFieldMappings = generate_field_mapping(origFC, fieldList)         tempfc = os.path.join(scratchWS, "temp_{}".format(name))         arcpy.FeatureClassToFeatureClass_conversion(origFC, scratchWS, outFC, "", stringFieldMappings)         print "  > done with {} ({} seconds)".format(name, round((time.time() - fStart), 2))      But, I have one feature class for which we want _all_ of the fields dropped and I am not sure if I can accomplish this in the existing function workflow.   * If I pass an empty list of fields, it breaks: `FieldMap: Error in adding input field to field map`   * If I pass an empty (`""`) value, and put in a conditional to skip fieldmapping if there's no list, then all the original fields are passed.   * If I pass `OBJECTID` (logic being, that would definitely be in the output), it breaks: `Field mapping error: merge rule not valid for output field OBJECTID` I'm happy to just bypass the fieldmapping part of the function with a conditional, but I am not sure how to drop all the resulting attributes. (I don't want to DeleteField_management each one, that would be clumsy and I switched to FieldMappings to get away from that.) **Any ideas for a workaround?**


title: Alternate option for a Data Player plugin for GPS in qgis

text: I found a plugin for QGIS to play a GPS Track. (Data player plugin).The Problem is I can not read .gpx File with this plugin. Is there another option?


title: How do I maintain a resposive GUI using QThread with PyQGIS

text: I have been developing some batch processing tools as python plugins for QGIS 1.8. I have found that while my tools are running the GUI becomes non-responsive. The general wisdom is that the work should be done on a worker thread, with the status/completion information passed back to the GUI as signals. I have read through the riverbank docs, and studied the source of doGeometry.py (a working implementation from ftools). Using these sources I have tried to build a simple implementation in order to explore this functionality before making changes to an established code base. The overall structure is an entry in the plugins menu, that lanuches a dialog with start and stop buttons. The buttons control a thread that counts to 100, sending a signal back to the GUI for each number. The GUI receives each signal and sends a string containing the number both the message log, and the window title. The code of this implementation is here:               from PyQt4.QtCore import *     from PyQt4.QtGui import *     from qgis.core import *          class ThreadTest:              def __init__(self, iface):             self.iface = iface              def initGui(self):             self.action = QAction( u"ThreadTest", self.iface.mainWindow())             self.action.triggered.connect(self.run)             self.iface.addPluginToMenu(u"&ThreadTest", self.action)              def unload(self):             self.iface.removePluginMenu(u"&ThreadTest",self.action)              def run(self):             BusyDialog(self.iface.mainWindow())          class BusyDialog(QDialog):         def __init__(self, parent):             QDialog.__init__(self, parent)             self.parent = parent             self.setLayout(QVBoxLayout())             self.startButton = QPushButton("Start", self)             self.startButton.clicked.connect(self.startButtonHandler)             self.layout().addWidget(self.startButton)             self.stopButton=QPushButton("Stop", self)             self.stopButton.clicked.connect(self.stopButtonHandler)             self.layout().addWidget(self.stopButton)             self.show()              def startButtonHandler(self, toggle):             self.workerThread = WorkerThread(self.parent)             QObject.connect( self.workerThread, SIGNAL( "killThread(PyQt_PyObject)" ), \                                                     self.killThread )             QObject.connect( self.workerThread, SIGNAL( "echoText(PyQt_PyObject)" ), \                                                     self.setText)             self.workerThread.start(QThread.LowestPriority)             QgsMessageLog.logMessage("end: startButtonHandler")              def stopButtonHandler(self, toggle):             self.killThread()              def setText(self, text):             QgsMessageLog.logMessage(str(text))             self.setWindowTitle(text)              def killThread(self):             if self.workerThread.isRunning():                 self.workerThread.exit(0)               class WorkerThread(QThread):         def __init__(self, parent):             QThread.__init__(self,parent)              def run(self):             self.emit( SIGNAL( "echoText(PyQt_PyObject)" ), "Emit: starting work" )             self.doLotsOfWork()             self.emit( SIGNAL( "echoText(PyQt_PyObject)" ), "Emit: finshed work" )             self.emit( SIGNAL( "killThread(PyQt_PyObject)"), "OK")              def doLotsOfWork(self):             count=0             while count < 100:                 self.emit( SIGNAL( "echoText(PyQt_PyObject)" ), "Emit: " + str(count) )                 count += 1     #           if self.msleep(10):     #               return     #          QThread.yieldCurrentThread()      Unfortunately it's not quiet working as I hoped:   * The window title is updating "live" with the counter but if I click on the dialog, it is non-responsive.   * The message log is inactive until the counter ends, then presents all of the messages at once. These messages are tagged with a time stamp by QgsMessageLog and these time stamps indicate they were received "live" with the counter i.e. they are not being queued by either the worker thread, or the dialog.   * The order of messages in the log (exerpt follows) indicates that startButtonHandler completes execution before the worker thread comences work i.e. the thread behaves as a thread.              end: startButtonHandler     Emit: starting work     Emit: 0     ...     Emit: 99     Emit: finshed work        * It seems the worker thread just isn't sharing any resources with the GUI thread. There are a couple of commented out lines at the end of the above source where I tried calling msleep() and yieldCurrentThread(), but neither seemed to help. Is anyone with any experience with this able to spot my error? I'm hoping it is a simple but fundamental mistake that is easy to correct once it is identified.


title: Is PgQuery for QGIS still available?

text: I was wondering if PgQuery for QGIS is still available. I am trying to install it on some new machines with QGIS 1.8.0 and I can not find it in the Python Plugin Installer. Does anyone know of a place where I can download it directly? Do I need to add a specific repository? I currently only have the Official Repository loaded. Thanks...


title: GDAL/OGR Determine if dataset is raster or vector

text: I am writing a script using the GDAL/OGR Python bindings and it is possible that the user may provide a raster dataset or a vector dataset as input. If the user provides a vector dataset, the script will convert it to raster. Is there a way to determine whether the dataset is raster or vector before opening them? I could just try to open the dataset with OGR first and if it returns None I could try to open with GDAL, but I am looking for a more simple way like `gdal.IsRaster(filepath)` if it exists.


title: GeoServer FeatureID Filter Limits

text: we've made an application that uses GeoServer WMS tiles on an OpenLayers map. We have a layer of Bus Stops which we filter based on search results from previous search i.e. User searches for 'The Square', we find 20 stops that say 'The Square' and filter them in the WMS requests using FeatureID like so. > > http://wms.geowessex.com/geoserver/wms?LAYERS=dev%3ADFT_BUS_STOPS&CACHE=true&TRANSPARENT=true&FORMAT=image%2Fpng8&EXCEPTIONS=application%2Fvnd.ogc.se_xml&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FEATUREID=DFT_BUS_STOPS.2106%2CDFT_BUS_STOPS.239%2CDFT_BUS_STOPS.357%2CDFT_BUS_STOPS.1217%2CDFT_BUS_STOPS.1185%2CDFT_BUS_STOPS.245%2CDFT_BUS_STOPS.1175%2CDFT_BUS_STOPS.2261%2CDFT_BUS_STOPS.1199%2CDFT_BUS_STOPS.1665%2CDFT_BUS_STOPS.1664%2CDFT_BUS_STOPS.1214%2CDFT_BUS_STOPS.1181%2CDFT_BUS_STOPS.1712%2CDFT_BUS_STOPS.1258%2CDFT_BUS_STOPS.2078%2CDFT_BUS_STOPS.1741%2CDFT_BUS_STOPS.1176%2CDFT_BUS_STOPS.1165%2CDFT_BUS_STOPS.1137%2CDFT_BUS_STOPS.1068%2CDFT_BUS_STOPS.2084%2CDFT_BUS_STOPS.353%2CDFT_BUS_STOPS.1139%2CDFT_BUS_STOPS.1248%2CDFT_BUS_STOPS.203%2CDFT_BUS_STOPS.1731%2CDFT_BUS_STOPS.205%2CDFT_BUS_STOPS.1239%2CDFT_BUS_STOPS.2453%2CDFT_BUS_STOPS.1679%2CDFT_BUS_STOPS.1678%2CDFT_BUS_STOPS.356%2CDFT_BUS_STOPS.1701%2CDFT_BUS_STOPS.1715%2CDFT_BUS_STOPS.250%2CDFT_BUS_STOPS.1704%2CDFT_BUS_STOPS.358&SRS=EPSG%3A900913&BBOX=-275784.79880351,6550959.0701793,-275173.30257734,6551570.5664055&WIDTH=256&HEIGHT=256 (this particular example only shows two of the filtered stops, but you get the picture) This works fine, but with some of the search results the amount of stops can be significantly higher, so we end up with a request like this. > > http://wms.geowessex.com/geoserver/wms?LAYERS=dev%3ADFT_BUS_STOPS&CACHE=true&TRANSPARENT=true&FORMAT=image%2Fpng8&EXCEPTIONS=application%2Fvnd.ogc.se_xml&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FEATUREID=DFT_BUS_STOPS.2106,DFT_BUS_STOPS.239,DFT_BUS_STOPS.357,DFT_BUS_STOPS.1217,DFT_BUS_STOPS.1185,DFT_BUS_STOPS.245,DFT_BUS_STOPS.1175,DFT_BUS_STOPS.2261,DFT_BUS_STOPS.1199,DFT_BUS_STOPS.1665,DFT_BUS_STOPS.1664,DFT_BUS_STOPS.1214,DFT_BUS_STOPS.1181,DFT_BUS_STOPS.1712,DFT_BUS_STOPS.1258,DFT_BUS_STOPS.2078,DFT_BUS_STOPS.1741,DFT_BUS_STOPS.1176,DFT_BUS_STOPS.1165,DFT_BUS_STOPS.1137,DFT_BUS_STOPS.1068,DFT_BUS_STOPS.2084,DFT_BUS_STOPS.353,DFT_BUS_STOPS.1139,DFT_BUS_STOPS.1248,DFT_BUS_STOPS.203,DFT_BUS_STOPS.1731,DFT_BUS_STOPS.205,DFT_BUS_STOPS.1239,DFT_BUS_STOPS.2453,DFT_BUS_STOPS.1679,DFT_BUS_STOPS.1678,DFT_BUS_STOPS.356,DFT_BUS_STOPS.1701,DFT_BUS_STOPS.1715,DFT_BUS_STOPS.250,DFT_BUS_STOPS.1704,DFT_BUS_STOPS.358,DFT_BUS_STOPS.204,DFT_BUS_STOPS.1198,DFT_BUS_STOPS.2363,DFT_BUS_STOPS.2449,DFT_BUS_STOPS.2450,DFT_BUS_STOPS.1260,DFT_BUS_STOPS.1707,DFT_BUS_STOPS.1699,DFT_BUS_STOPS.1700,DFT_BUS_STOPS.228,DFT_BUS_STOPS.1675,DFT_BUS_STOPS.1711,DFT_BUS_STOPS.1710,DFT_BUS_STOPS.1528,DFT_BUS_STOPS.1706,DFT_BUS_STOPS.1705,DFT_BUS_STOPS.243,DFT_BUS_STOPS.265,DFT_BUS_STOPS.1160,DFT_BUS_STOPS.1261,DFT_BUS_STOPS.227,DFT_BUS_STOPS.1164,DFT_BUS_STOPS.1177,DFT_BUS_STOPS.1736,DFT_BUS_STOPS.985,DFT_BUS_STOPS.263,DFT_BUS_STOPS.1748,DFT_BUS_STOPS.1266,DFT_BUS_STOPS.1178,DFT_BUS_STOPS.1163,DFT_BUS_STOPS.2125,DFT_BUS_STOPS.1061,DFT_BUS_STOPS.1780,DFT_BUS_STOPS.1268,DFT_BUS_STOPS.1747,DFT_BUS_STOPS.1152,DFT_BUS_STOPS.1738,DFT_BUS_STOPS.2531,DFT_BUS_STOPS.251,DFT_BUS_STOPS.1155,DFT_BUS_STOPS.2530,DFT_BUS_STOPS.1151,DFT_BUS_STOPS.1156,DFT_BUS_STOPS.1132,DFT_BUS_STOPS.1073,DFT_BUS_STOPS.1154,DFT_BUS_STOPS.252,DFT_BUS_STOPS.1722,DFT_BUS_STOPS.1723,DFT_BUS_STOPS.1744,DFT_BUS_STOPS.2104,DFT_BUS_STOPS.236,DFT_BUS_STOPS.1721,DFT_BUS_STOPS.1720,DFT_BUS_STOPS.1064,DFT_BUS_STOPS.195,DFT_BUS_STOPS.1745,DFT_BUS_STOPS.369,DFT_BUS_STOPS.1170,DFT_BUS_STOPS.858,DFT_BUS_STOPS.1070,DFT_BUS_STOPS.1193,DFT_BUS_STOPS.1192,DFT_BUS_STOPS.2420,DFT_BUS_STOPS.1670,DFT_BUS_STOPS.1671,DFT_BUS_STOPS.1245,DFT_BUS_STOPS.2412,DFT_BUS_STOPS.2427,DFT_BUS_STOPS.206,DFT_BUS_STOPS.1255,DFT_BUS_STOPS.1256,DFT_BUS_STOPS.1742,DFT_BUS_STOPS.1072,DFT_BUS_STOPS.1133,DFT_BUS_STOPS.1734,DFT_BUS_STOPS.1733,DFT_BUS_STOPS.293,DFT_BUS_STOPS.2447,DFT_BUS_STOPS.2448,DFT_BUS_STOPS.1144,DFT_BUS_STOPS.2347,DFT_BUS_STOPS.249,DFT_BUS_STOPS.1693,DFT_BUS_STOPS.1692,DFT_BUS_STOPS.197,DFT_BUS_STOPS.1142,DFT_BUS_STOPS.1784,DFT_BUS_STOPS.1783,DFT_BUS_STOPS.968,DFT_BUS_STOPS.1785,DFT_BUS_STOPS.1143,DFT_BUS_STOPS.1781,DFT_BUS_STOPS.1782,DFT_BUS_STOPS.202,DFT_BUS_STOPS.986,DFT_BUS_STOPS.862,DFT_BUS_STOPS.2107,DFT_BUS_STOPS.247,DFT_BUS_STOPS.898,DFT_BUS_STOPS.237,DFT_BUS_STOPS.248,DFT_BUS_STOPS.1191,DFT_BUS_STOPS.1265,DFT_BUS_STOPS.927,DFT_BUS_STOPS.1737,DFT_BUS_STOPS.261,DFT_BUS_STOPS.262,DFT_BUS_STOPS.1189&SRS=EPSG%3A900913&BBOX=-275784.79880351,6550959.0701793,-275173.30257734,6551570.5664055&WIDTH=256&HEIGHT=256 Now obviously I have hit some sort of limit here, my main question is, is the limit based on the length of the url, or is it based on the amount of FeatureIDs I am filtering by. I think ultimately I am going to have to put a limit in place on the filter, so that if there are more than 120 results (my request failed at 147 FeatureID filters, this gives some wiggle room) the results do not get filtered, unless anyone has any other ideas? Thanks!


title: How to join two table based on one unique column in PostGIS

text: I wonder if there any simple query that perform just like Join features of ArcGIS Desktop.


title: get valules for all fields from selected features

text: The following code works for retrieving values within a particular field(Name1), but I want to retrieve mulitple fields. Am I to create multiple variables for each field?               IMxDocument pMxDocument = ArcMap.Application.Document as IMxDocument;             IMap pMap = pMxDocument.FocusMap;             ILayer pLayer = pMap.get_Layer(0);             IFeatureLayer pFeatureLayer = (IFeatureLayer)pLayer;             IFeatureClass pFeatureClass = pFeatureLayer.FeatureClass;                  ITable pTable = (ITable)pFeatureClass;             ICursor pCursor = pTable.Search(null, false);             int fldIndex = pFeatureClass.Fields.FindField("Name1");             IRow pRow = pCursor.NextRow();             while (pRow != null)             {                 MessageBox.Show(pRow.get_Value(fldIndex).ToString());                 pRow = pCursor.NextRow();                  }      Edit: What I'm trying to do is get the selected features of a feature class, get the values for fields like Name, address, city, etc and then pass those values onto a document. I thought I would give some background to what I am doing since I may be on the wrong track to start off with.


title: AddData tool stop working after Interoperability connection has been used?

text: Bug in ArcGIS Engine? Has anyone seen this behavior? We get an exception from our ArcGIS Engine test projects after we had been testing the Interoperability connection in ArcEditor. We had to enable unmanaged code debuging in VS because the test application crashed without telling what happened, we get a stacktrace that look like this:               'TestMapWFSLayer.exe': Loaded 'C:\Test\TestMapWFSLayer\TestMapWFSLayer\bin\Debug     ...     'TestMapWFSLayer.exe': Loaded 'C:\Windows\SysWOW64\rasadhlp.dll'     First-chance exception at 0x7680b727 in TestMapWFSLayer.exe: Microsoft C++ exception: FNPNS::TSM::CDoesNotExistException at memory location 0x003799c4..     ...          'TestMapWFSLayer.exe': Loaded 'C:\Windows\SysWOW64\wship6.dll'     'TestMapWFSLayer.exe': Loaded 'C:\Windows\SysWOW64\WSHTCPIP.DLL'     The thread 'Win32 Thread' (0xd90) has exited with code 0 (0x0).     The thread 'Win32 Thread' (0xfd8) has exited with code 0 (0x0).     First-chance exception at 0x7680b727 in TestMapWFSLayer.exe: Microsoft C++ exception: FNPNS::TSM::CDoesNotExistException at memory location 0x00379974..     ...     'TestMapWFSLayer.exe': Loaded 'C:\Windows\SysWOW64\icmp.dll'     First-chance exception at 0x7680b727 in TestMapWFSLayer.exe: Microsoft C++ exception: FNPNS::TSM::CDoesNotExistException at memory location 0x00379f20..     First-chance exception at 0x7680b727 in TestMapWFSLayer.exe: Microsoft C++ exception:      ...     First-chance exception at 0x7680b727 in TestMapWFSLayer.exe: 0xE0000001: 0xe0000001.     ...           ...     First-chance exception at 0x7680b727 in TestMapWFSLayer.exe: 0xE0000001: 0xe0000001.     ...     ...     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\DIDatasource.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ArcGIS\Engine10.0\bin\fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Windows\SysWOW64\wsock32.dll'     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\xerces-c_2_8_fme.dll'     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\icuuc38_fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\icudt38_fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\tcl_fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\threadlockmgr.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Windows\SysWOW64\shfolder.dll'     'TestMapWFSLayer.exe': Loaded 'C:\Windows\winsxs\x86_microsoft.vc90.openmp_1fc8b3b9a1e18e3b_9.0.30729.4974_none_80bb23871e9ab973\vcomp90.dll'     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\xerces-depdom_2_8_fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\icuin38_fme.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\fme_stringconverter.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\fme_sqlite3.dll', Binary was not built with debug information.     'TestMapWFSLayer.exe': Loaded 'C:\Program Files (x86)\ESRI\Data Interoperability Extension\res\fmearcgisext_res.dll', Binary was not built with debug information.     Critical error detected c0000374     Windows has triggered a breakpoint in TestMapWFSLayer.exe.          This may be due to a corruption of the heap, which indicates a bug in TestMapWFSLayer.exe or any of the DLLs it has loaded.          This may also be due to the user pressing F12 while TestMapWFSLayer.exe has focus.          The output window may have more diagnostic information.      Update: Tried uninstalling the Data Interoperability Extension (DIE): When I tried using the AddData tool after uninstalling DIE, it was i bit slow but after a while it gor visible and I was able to use it. Installed the DIE again and the AddData dialog workes until I hit the geodatabases tab.![AddData dialog](http://i.stack.imgur.com/QymoA.png)


title: It seems there is a bug in St_Clip function of PostGIS

text: The code I am using for clip feature is :               raster ST_Clip(raster rast, geometry geom, double precision[] nodataval=NULL, boolean crop=true);      I tried clip a shape depends on the polygon from the Raster Data. But the PostGIS would clip the raster depend on the Bounding Box of the polygon, and put the Non-data value as 0. Even I tried to modify that with code:               SELECT ST_Clip(rast,geometry, ARRAY[-9999],true) FROM dtm, polygons WHERE polygons.polygonid = 1;      After I dump the data as a ASC-II file from the PostGIS, using gdal_translate with another option -a_nodata -9999, I find out it only changes the Header of the ASC-II file not the exact bounding box value. ![enter image description here](http://i.stack.imgur.com/kCCTE.png) Is there anyway to change the value from 0 to -9999, or someway to clip the raster depending on the polygon rather than the Bounding Box of the polygon. THanks.


title: DEM data - Central Europe - noise in heights

text: I have downloaded DEM data from this page: http://www.viewfinderpanoramas.org/dem3.html. I have written simple program in C++ that takes "bricks" and build resulting height map. Hovewer, I have problem with result. I am using Mercator projection and take every "gps pixel" from DEM data and recalculate its position to real pixels (same for border). It works fine (or at least it seams to). As you can see on image, there is something weird. That "mess" in lower part is incorrect and I dont know, if data are incorrect, if I am using them incorrect (wrong bricks ?) or something in my program is wrong. ![Height map of Czech Republic](http://i.stack.imgur.com/R8fkI.png) EDIT: Here is algorithm for reading "bricks" heights:               //iterate all bricks and set its starting coordinates     double startLon = 19;     double startLat = 48;          ///-------------------------     int tileSize = 1201;     double step = 1.0 / tileSize;          for (int y = 0; y < tileSize; y++)     {               for (int x = 0; x < tileSize; x++)        {              short height = brick[x + (tileSize - 1 - y) * tileSize];             //tileSize - 1 - y -> reverse brick coordinates              height = swap_endian<short>(height);              if (height == -32768)         {                            height = 0;         }              lon = (startLon) + x * step;         lat = (startLat - 1) + y * step;              int pixelX = 0;         int pixelY = 0;         this->helper->ConvertCoordinates(lon, lat, &pixelX, &pixelY);              if (pixelX < 0) {  continue;}         if (pixelY < 0) {  continue;}         if (pixelX >= this->width) {  continue;}         if (pixelY >= this->height) {  continue;}                  //write it to map             //if multiple height in same pixel, take only last one         this->heightMap[pixelX + pixelY * this->width] = height;                     }     }


title: How do I build a custom menu to change preferences for a MapBasic application?

text: I would like to be able to input a few different tables for an application. I want to be able to set the path before the app runs. How do I create a preference box and what code do I use to call the tables that I set in the box? I want it to resemble this:


title: Software for analyzing GPS tracklogs to extract human-readable placenames

text: I have a pile of GPX-formatted tracklogs for many countries where I've travelled over the years (the GPS-device has been continuously on and recording for most of the trips). I would now like to analyse the logs automatically to extract a list human-readable names for the places that I've visited. For example:   * detect and name border crossings   * detect and name inhabited places where I've stayed longer Does there exist ready-made software (preferably free and open source) that can be used to perform such analysis? (By "tracklog" I mean: a sequence of timestamped longitude-latitude pairs.)


title: How to get a string from a Double with a fixed number of decimals?

text: I have a featureclass (in a File Geodatabase), with a particular field of type Double. I want to copy it over to a text field, so that only 3 decimal places are shown. For example: `34.78902345 becomes 34.789` I have gone through this question: How to convert a double field to a string?, but that does not solve my issue. I have even tried to format the double field, but even then the text field contains all the decimals. How do I do this in ArcGIS 10.0?


title: Which flow routing algorithms does ArcGis support?

text: Other than standard D8 are there support for multiple flow algorithms like FD8, DEMON, Dinfinity etc?


title: How to easily label point values in ArcScene?

text: I would like to be able to quickly label point values in my scene. I am aware of the fact that there is no out-of-the-box solution to label features in ArcScene similar to the one in ArcMap. (No doubt the 3D environment is presenting a lot of complexity to the programmers.) I've worked with some (expensive) programs that were able to successfully and cleanly label feature in 3D environment but redrawing/refreshing was required to update the labels. (Gemcom GEMS) I do not need my points to rotate correctly with the scene or anything as fancy as that however I am a little tired of placing text manually in the scene or in post-processing. Does anyone know of a working solution?


title: What generation of the MrSID format is supported by ArcGIS 9.3?

text: This is probably old news to some, but I just learned that MrSID (.sid) files have different 'generations' of file formats - so there are MrSID Generation 2, Generation 3, and Generation 4 files. Does anyone know what generation of .sid is supported in ArcGIS 9.3? Thanks!


title: Polygon gradient fills / tint bands in QGIS

text: Is it possible to use a gradient fill on a polygon in QGIS? I would like to have my polygons with a colour at the edges, fading to clear/white as it moves away from the edge. The majority of the centre of the polygon will be white or clear. It is something people I know who use Illustrator can do. I wondered if there was a way to emulate it in QGIS. Here is an example. I would like something like the purple shading. ![enter image description here](http://i.stack.imgur.com/b9MRt.jpg) Update: This method suggested below does work, however it produces some annoying extra bits for my polygons. For example, see the image below: ![enter image description here](http://i.stack.imgur.com/kWGqo.jpg) I am unsure how to remove these. I cannot use a mask as I need to be able to see the base mapping and other data outside of this polygon on layers that are below this layer. Using a mask would obscure all of these layers. I have included a link to Underdarks blog post about this, which gives more detail, and contains the code to set this up. http://underdark.wordpress.com/2011/08/08/creating-a-gradient-fill-for- polygons-in-qgis/


title: Geomajas: How do I add context menu on a map marker

text: I'm evaluation Geomajas. Apparently, there's very little info on this framework across the Internet. I've managed to create a context menu for the map (the menu that appears when you right click on the map). I want to do the same for a map point marker (I may need to extend this to other geometry types, i.e., Polygons, Multipolygons, Lines, etc). Anyone with an idea on how to do this, please assist. Additionally, when I add context menu, the "About Geomajas" menu item hides my first menu item. Is there a way that I can add menu items from below this menu, or make it the last item in the menu? I'm adding context menu by adding the line:               map.getContextMenu().addChild(MyMenus.getMyContextMenu());


title: How to delete records from multiple, similar feature classes by attribute values using an SQL statement?

text: I have several hundred geodatabase feature classes which all have the exact same fields. I need to go through each of those feature classes and delete any records which match an SQL query (the same SQL query can be applied to all the feature classes). Alternatively, I could export them all to new feature classes, keeping the records I want. In ArcCatalog I can export them singly using an SQL query to specify which records get exported, but the SQL statement option isn't available when you want to export multiple layers simultaneously (using "Feature Class to Geodatabase (multiple)").


title: Performing Geoprocessing SQL Queries in QGIS from a PostGIS Database

text: This may be a dumb question, but is there way to run a PostGIS geoprocessing SQL query in QGIS? Specifically, I am thinking of running SQLs such as ST_Union, ST_Split etc etc. I have tried the following SQL in the QGIS PostGIS Query Editor without any luck, while the query worked in PostGIS.               SELECT covtype, ST_Multi(ST_Union(f.the_geom)) as sGeom     INTO testLCC     FROM lcc2000_001k AS f     GROUP BY covtype;          SELECT populate_geometry_columns     ('testLCC'::regclass);          ALTER TABLE testLCC ADD CONSTRAINT pkid PRIMARY KEY(covtype);      Thanks in advance...


title: ArcGIS 10.0/arcpy: How to create a polygon geometry from inner and outer ring point arrays?

text: I have two arrays, each containing a closed ring of points, one for an outer ring and the other for an inner ring, the latter running counterclockwise, of course. I am using the arcpy.Polygon() factory function to build the polygon. **My question: how do I input the points to Polygon()?** Documentation gives the general syntax as Polygon(inputs, spatial reference, has_Z, has_M). The inputs part is what is sketchy. The inputs supposedly can be point objects or arrays. My problem is that I have a polygon with a doughnut hole. I can build a simple polygon with a single point array as input ... no problem; however, in what way do I add the _two_ sets of points? In the "Reading Geometries" help page it says that parts (in this case rings) are separated by a null point object. For polygons with holes, the outer ring goes first. So, I tried adding the outer array then the inner array to a container array. That didn't work: no hole. I tried the same thing with a null point object between them. No joy: same result. Then, I tried putting both arrays into a Python list as the input. That threw an error. The sample scripts in the help files never deal with _writing_ polygons with holes. So, I really don't know what format to use for the inputs. Does anyone have any insight into this one? I would surely appreciate the help. Thanks, all.


title: How to grid unevenly sampled categorical data?

text: I am looking for a method for gridding categorical data. I have extracted from nautical charts and field sheets a collection of points specifying the nature of the surface of the bottom of the ocean. These data are categorical rather than numerical and they are not regularly nor even randomly sampled. Nautical charts are created to aid in navigation and anchoring; they are not created to map habitat. As such, more soundings are made close to shore where relatively shallow depths can pose a hazard to navigation and where ships to tend to anchor. Further from shore, where the depths are more than adequate for navigation and anchoring is impractical, soundings are made much less frequently. **Has anyone else tried to create gridded substrate maps from nautical charts?** I looked at Thiessen (Vornoi) polygons but the concentration of soundings along shores leads to a fine 'honeycomb' along the shore, large polygons offshore and in between long pie-shaped polygons stretching way offshore. Gridding using nearest neighbour yields pretty much the same result. I need a way to limit the influence of the shallow, near shore points -- a way to limit those long pie-shaped polygons. In deeper waters I do not expect that the nature of the bottom will be a continuation of the near shore bottom. I started thinking along two lines -- both using depth. **One** is weighting the choice of the 'nearest' neighbour using the difference in depth between a grid cell and neighbouring points. **Another** is de-selecting neighbouring points which are more than some specified tolerance different in depth. Or, perhaps rather than a pre-specified tolerance, I could bin depth ranges and then limit the choice of neighbouring points to those in the same depth range or bin. **Any thoughts on how to implement either of these two options?** Since talking with colleagues on other forums, I have been looking at a couple of other approaches. The first involves using a barrier -- the 100m depth contour -- to limit the influence of near shore data. The challenge with this approach is that any of the ESRI interpolation routines that can use barriers are designed to work with continuous data rather than discontinuous data. I could use the barriers to break up the points into shallow near shore and deeper points before creating Thiessen polygons. However, I anticipate rampant edge effects since ArcGIS createsThiessen polygons for rectangular areas not for complex areas. The second approach -- suggested by several colleagues -- was kriging. I had initially dismissed kriging out of hand because I have only ever considered it for continuous data. The challenge with kriging is that it too is not designed for categorical data. Now, I am looking at cokriging with depth and nature of surface but, any type of kriging is going to involve using integer numerical codes for nature of surface. Afterward the resulting floating point numerical codes will have to be reduced back to the original integer coding. Not pretty. Can any one suggest other lines to follow? (One could, perhaps, use terrain analysis. For example, slopes steeper than the angle of repose could not be sediments. I am looking for something simpler and, at any rate, I do not have data at sufficient spatial resolution.) Regards,


title: Problem installing gdal - works in build directory, not in /usr/bin

text: I've been trying to install gdal with the FileGDB driver. I installed it fine on a linux mint machine. I am having problems installing on a CentOS machine though.   It builds fine, I can run the apps in the build directory but when I run make install, it runs without error.   The machine did have gdal 1.9 installed. After we had tried several times installing I used rpm to uninstall the gdal packages shown as installed. Reran the make install, reran ldconfig just in case. No difference. For now I'm hard coding the path to the build directory, but obviously I'd rather have it run properly.               ./configure LDFLAGS=-L/lib64 --with-fgdb=/data/HOMEDIRS/SAI/dev/src/FileGDB_API      --without-python --with-perl --with-pg=/usr/bin/pg_config --prefix=/usr     /data/HOMEDIRS/SAI/dev/gdal-1.11.0/apps/ogrinfo --formats|grep File      -> "MapInfo File" (read/write)      -> "OpenFileGDB" (readonly)      -> "FileGDB" (read/write)     /usr/bin/ogrinfo --formats|grep File      -> "MapInfo File" (read/write)      -> "OpenFileGDB" (readonly)      Interestingly, if I delete the file in /usr/bin then reinstall, the install does something to the binary (link in libraries?):                ls -l /data/HOMEDIRS/SAI/dev/gdal-1.11.0/apps/ogrinfo /usr/bin/ogrinfo      -rwxr-xr-x 1 nrgadmin root  7671 Sep  4 20:00 /data/HOMEDIRS/SAI/dev/gdal-1.11.0/apps/ogrinfo      -rwxr-xr-x 1 root     root 78727 Sep  4 20:17 /usr/bin/ogrinfo


title: how to make a join with sql server table

text: I have a question. I'd like to join .shp to table by using AddJoin_management. But I'd like to make join to sqlserver table (not to *.dbf). Is it possible to do? it should look something like this:               AddJoin_management ("my_layer.shp", "ID", SQL-SERVER-TABLE, "ID")


title: QGis and developpement version

text: On various forums, when I asked some question, I have got more and more answer like that: "you can do this with the development version of QGis 1.9... but you can't with the current version..." It is quite frustrating! Why QGis don't do like some other opensource software, propose the both version on their website (the current and the development one with a message which inform that it is not a stable version)? I know that I can found this development version, but I find it quite a shame to have to looking for this version and not to find it on the main website, especially for software that tries to be very friendly user. User of several opensource software's, I find it quite a shame. Just to share my point of view, don't worry, I'm not angry against QGIS, which I find excellent, it's just a reaction to the answers on many forums that refer to this version. Cheers.


title: How can I convert files from a county government to shapefiles?

text: I work for a company that creates software for multiple listing services and real estate agents. We're expanding our coverage area into a new market, and I've been given some files for two different counties which they are calling their "shapefiles". One file is an access database, and the other is a fixed width column text file. Is there some way for me to convert these files to actual ESRI shapefiles? Here are some examples - **Here's what I'm given in the fixed width file for the "vector" column (from Palm Beach county, Florida):** A0CL10D10L12D40R24U40L2U10 A4L30D2L2D6R2D14L2D6R2D7D11L7D14L14D20R26R8U6L8D6U6R8D11R17U14R5U6L5D6U6R5D6R8D6R14U38L14D12L8CL5D7R5U7 A1L24D12L5D10R5D28R12U4CR22U20L10D10L12D10 I have no idea what the above could be. If you want to see the file here's a small sample: sample.txt **The Access file (from Broward County, Florida):** The access file just says "Long binary data" for the geometry columns. It has a lot of tables named GDB_xxxx (GDB_SpatialRefs, GDB_GeomColumns, etc). I'm guessing this has something to do with Microsoft's spatial format. Any help you guys can give would be greatly appreciated!


title: Solution for slow display of postgis-tables ? (qgis 1.8 + postgis 2.x)

text: I have the problem that my postgres-tables are displayed very slow when I scroll them in QGIS. I have no problem with other tested GIS (e.g. UDig, GVSig, ArcGIS). The number of objects isn't large (5000-10000) and the rendering is no problem. I've observed that every action like scrolling causes queries in the database, looks like there's no cache - Is that true??? 'Can't believe that. The connection to the database server runs through a fast Internet (and it works with other GIS). Does anyone have an idea how to solve the problem? I would love to work with QGIS, but the problem makes it impossible so far. I appreciate your advices, regards Bob


title: Digital Number (DN) on satellite imagery?

text: How would I know if any of the satellite imagery (Landsat, ASTER, Quickbird and so on) have DN or Digital Numbers such as integer or uninteger on its pixel. I download from the USGS EarthExplorer and have seen some in 0 -255 or into floating numbers. My question do they at USGS re-calibrated or left untouched them before they put in into their website for the user to download. Does this means that each different imagery has its own different digital numbers from the sensor ? I have ArcGIS 10.1, ENVI 5.0 and ERDAS 2013


title: Persistent datum conflict error. Datums are the same

text: I'm clipping a number of features using other feature layers, and was careful to get everything into WGS84 beforehand, using what I thought were the proper datum transformations. Now I keep getting Error 001003, Datum conflict between input and output. I'm not sure if this is referring to my input shapefiles, or if it's a problem with output coordinate environment settings, or none of the above. ![Feature 1 Coordinate System Info](http://i.stack.imgur.com/WAf8O.jpg) ![Feature 2 Coordinate System Info](http://i.stack.imgur.com/b3hFl.jpg) Contents of the .prj file of the input feature: > > GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],VERTCS["Unknown > VCS",VDATUM["Unknown"],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]] Contents of the .prj file of the clip feature: > > GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] Any ideas?


title: How do I create an index layer for multiple maps?

text: I am trying to create an index layer of just the area within the Orange lines. I am having difficulty creating this layer due to it not being a simple box area. The goal is to create a index layer, with grids, so that the Data Driven Pages would be used to generate it as a mapbook for the response area for my firehouse. Attached is an image of what I am speaking about. Any ideas? I am using ArcGIS but should I be QGIS? I am fairly a beginner but pick up quickly thanks in advance. ![enter image description here](http://i.stack.imgur.com/L5XJw.png)


title: Drawing tools to add circles, arrows, lines for QGIS?

text: I'm pretty new in QGIS and i'm searching a drawing tool like it is possible to find in microsoft word or excel. I search something very user friendly (for dummies!) to add circles, arrows, lines... with no need to create a new layer!! Can you help me?


title: Why does polygon.intersection(complex_linestring) split the linestring and mess up Z?

text: I have a complex (self crossing) line with Z values which I'd like to truncate to the area inside a box. The problem is Shapely splits the line into simple (non-crossing) sections and the Z value at the intersection of the loop is incorrect. My high-level goal is taking a pile of GPS tracks, trimming them to only include the area in a box and calculate speed within the box. I'm using Z to represent seconds. I suspect this problem is deep in geos. Can you suggest a better approach? I made a pretty picture of this example at http://i.stack.imgur.com/pQxtH.png but it may not work since I'm new here.               from shapely.geometry import LineString     import shapely.geometry     box = shapely.geometry.box(1, 1, 4, 5)     line = LineString([(2, 0, 0), (2, 4, 4), (3, 4, 5), (3, 3, 6), (0, 3, 9)])     [[p for p in l.coords] for l in line.intersection(box)]      results in:               [[(2.0, 1.0, 1.0), (2.0, 3.0, 5.0)],      [(2.0, 3.0, 5.0), (2.0, 4.0, 4.0), (3.0, 4.0, 5.0), (3.0, 3.0, 6.0), (2.0, 3.0, 5.0)],      [(2.0, 3.0, 5.0), (1.0, 3.0, 8.0)]]      but I was expecting only this segment:               [(2,1,1), (2,4,4), (3,4,5), (3,3,6), (1,3,8)]               shapely.__version__     '1.2.17'     shapely.geos.geos_version_string     '3.2.2-CAPI-1.6.2'


title: How can I use Android to add some additional data to geotiff images?

text: I would like to create an application which will use a geotiff image and add some additional data on it. All these with the Android Platform. I've done some research and found GDAL but I don't know if there is something else which exists? Thanks Martin


title: Choosing between Dynamic layers and Feature layers

text: I'm building a client using the ArcGIS WPF libraries, that's displaying my entities on a world map. I'm using an `ArcGISDynamicMapServiceLayer` to display the geographic and topographic base map, and I'm now wondering what's the best way to show my custom entity layer. At first, I used a `FeatureLayer` object, linked to the layer's published URL, to display the entities on the map. This worked great in that it showed the entity symbology and allowed me to interact with it (catch right-clicks, select/deselect), but it doesn't show the labels associated with the entity, just the symbology. So I tried using a second ArcGISDynamicMapServiceLayer to display the entities, which gives me data labels, but without the interactivity. So now, since these are both client requirements, I find myself using BOTH, and synchronizing visibility between the two instances of the layer, once as a FeatureLayer, and once under a dynamic layer. This is not a good practice, and bugs rear their head constantly. So my question is - what's the best way to achieve those two goals? Is there a simple way to add my own data labels on the map for each entity? There could be hundreds of them on 5-6 of these entity layers, and this feels like replicating functionality that ArcGIS already offers. On the other hand, I can simulate the right-click-on-entity behavior using a DynamicMap, catching the MouseClick event, and running a Query action on the map to find the entities clicked, but this, again, is reinventing parts of the wheel, and leaves some unsolved issues. So what would you do? How do I get those two requirements together without ugly hacks?


title: How to get pgrouting source and target value using openlayers 2/3?

text: Is there anyone can teach me how to get the pgrouting's source and target parameters in openlayers 2/3? I plan to get those parameters by click the OSM map layer. Currently, I can only get the lonlat x and y values, but my this is not exactly what I need. For those willing to help, thank you in advanced.


title: How to control line width in legend in QGIS?

text: I have a map in which I am defining the sizes of various features in map units. When I come to generate a pdf of the map, using the Print Composer, the map itself is fine, but when I include a legend then the symbol elements of the lines are represented at a totally inappropriate size (e.g. a 'Main road' symbol swamps the whole legend). I think it is the use of map units that is the cause of the problem as switching the line width to absolute units (e.g. mm) on the original map produces an appropriately sized line symbol in the legend. It seems like the legend that appears in the Print Composer is directly reflecting the way it appears in the Layer TOC - where, again, using map units produces oversized symbols relative to the map. However this is only true for lines. Point symbols rendered in map units appear exactly as I would expect in the Print Composer legend. I would be very grateful if anyone can enlighten me as to whether this is a known problem, my misunderstanding of something, or there is a solution/setting/trick that I have missed. I'm using QGIS 1.7.4 on Mac OSX 10.6. Thanks for any help. Phil


title: Go to X,Y dialog box

text: If someone gave you a x,y coordinate numbers he/she got from. I put the numbers in the X,Y events on my ArcMap and they are in Meters numbers. We entered the number into the x,y dialog box. The Arcmap display in the wrong place. I checked the bottom on the right corner and says they are in meters. I am thinking maybe that person could have gotten the numbers from Garmin GPS, or probably from the internet of that site where they went to. Is it possible that ArcMap gave the wrong events if I entered the UTM numbers ? Is it a bug or what? I would like to know from anyone who have been using it ? Thanks ! Using ArcMap 9.3.1


title: How to edit an ArcGIS Server WFS-T service in uDig?

text: Now a straight forward workflow to do this is present in the below link. http://udig.refractions.net/files/docs/uDigWalkthrough2.pdf Point 3 of page 8 states to change to the ‘Edit Geometry’ tool from the tool bar. ![enter image description here](http://i.stack.imgur.com/moJET.png) But I can’t find this anywhere in the uDig interface. ![enter image description here](http://i.stack.imgur.com/rbNBF.png) What could I be missing? Do note that the service added to uDig is a WFS-T service from ArcGIS Server and not GeoServer.


title: Equally sized countries?

text: Is there a map projection (or a more appropriate term) which allows you to see every country as an equally sized land area, but positioned approximately in the correct geographic position? Water could still be incorporated as it shows continental separation, etc.


title: Arcpy raster calculation producing random no data cell in output

text: I've been trying to get the raster calculation below to work using python:               import arcpy     from arcpy import env     import os     from arcpy.sa import *     arcpy.CheckOutExtension("spatial")          KSat_mm = "C:\\GIS\\Sat_Cond_V_mm.tif"     m = "C:\\GIS\\vG_m.tif"     env.workspace = "C:\\GIS\\WC_Dimensionless"     for raster in arcpy.ListRasters():       WC_unsat_workspace = "C:\\GIS\\WC_Unsat"       dim_m = (1/Raster(m))       WC_unsat_out = (Raster(KSat_mm) * (Raster(raster)**0.5)) * (1 - (1 - Raster(raster)**dim_m)**Raster(m))**2       WC_unsat_out.save(os.path.join(WC_unsat_workspace, raster.replace("WC_dim", "WC_unsat")))      The issue I'm having is that the output is giving me quite a few cells that are NoData when none of the input data has NoData cells in the same location. It might just be something I'm missing in the above code but any input would be helpful thank you. I'm using ArcGIS 10. All of the inputs are .tif files in a double precision format. Let me know if you need any more information. UPDATE: Ok the issue seems to be with negative values in the input datasets not sitting with the equation rather than the code itself.


title: How to derive Z-score values for a given field

text: I have a shapefile with an attribute field with a range of numerical values. I would like to derive z-score values for this field within ArcGIS. Any tips?


title: Exporting .tif from referenced mosaic dataset does not use functions in raster dataset

text: Very frustrating little issue in ArcGIS 10.1 -- I have a referenced mosaic dataset of terrain data that has a series of functions inserted (slope and several iterations of a statistics function to simplify and smooth the resulting mosaic). When I export the project to a .ai file (OR even the data as a .tif), the resulting .ai or .tif file produces my slope image WITHOUT the functions (i.e. unsmoothed). Anybody think of anything i might not have tried to make it work? Why would it display the mosaic with the functions running on the screen, but not the graphics without it?


title: How to get shapefile of river from OpenStreetMap?

text: I dont want to download whole map of an area. I just want rivers and other waterbodies present in an area in shapefiles so that I can use them in Qgis. I tried geofabrik but it is giving whole map and that also of full country while i need data of a city. Right now I am trying JOSM. will update you when it will works. A noob here.


title: How do I calculate acreage grouped by owner, masked by a poly?

text: I have statewide land ownership data and have created a mask to project ownership for only one county. I need to calculate the acreage of each owner within that county. Is there a quick way to achieve this? I'm using Arc10. Thanks for any help!


title: Creating donut holes without losing existing polygons using ArcGIS for Desktop?

text: Using ArcGIS 10.1 for Desktop is there a way to do a similar operation to creating a donut hole or island in a polygon without losing the existing features? I have mapped small lakes as several different polygons based on vegetation type. However, I want to map the vegetation surrounding the lakes, too - but if I create a circular polygon over the top of the lake and the surrounding area, it includes the area I've already mapped and given attributes to. So I want to 'cut out' a hole in the middle of the new polygon without losing the polygons that are now covered by it.


title: How show filter image out of feature in other color?

text: I am using openlayer js to filter the features display. I am trying to show the filtered imaged in other color keeping rest of feature also. Following is the code i am using for filtering:                  var a="73.84";        var b="73.89"        var filter= new OpenLayers.Filter.Logical({          type: OpenLayers.Filter.Logical.AND,        filters: [         new OpenLayers.Filter.Comparison({            type: OpenLayers.Filter.Comparison.BETWEEN,              property: "x",              lowerBoundary:a,      upperBoundary:b                 }),                 new OpenLayers.Filter.Comparison({                     type: OpenLayers.Filter.Comparison.BETWEEN,                     property: "y",                     lowerBoundary:"18.51",                     upperBoundary:"18.56"                 })]         });                   var parser = new OpenLayers.Format.Filter.v1_1_0();         var filterAsXml = parser.write(filter);         var xml = new OpenLayers.Format.XML();         var filterAsString = xml.write(filterAsXml);              allwards.params["FILTER"] = filterAsString;              allwards.redraw();      Using this code my feature is filtered but original feature are removed and only filted feature are show. I want to show original features along with filtered feature in different color. Please help me in this case.


title: How to use Raster Calculator in ArcGIS 9.3?

text: Raster Calculator tool is not available in ArcGIS 9.3 , so how to do in 9.3.??? thanks in advance


title: What is the problem with the polygon? Getting a *"found non-noded intersection"* exception error

text: I am working on a GIS program in which I needed to find the difference between the two polygons so, I decided to used NetTopologySuite (NTS) to find the difference. It seems to be working on some of the polygons and on some it is not working at all, it gives an _"found non-noded intersection"_. So, I drilled down more and found that it is not the difference function but it was the polygon which was screwing everything up. Here is the polygon which was giving exception _"found non-noded intersection"_ : > POLYGON ((51020.548917151937 215, 51025.909273225909 215, 51011.728078042863 > 215, 51009.760352395453 215, 51000.668102852258 215, 50988.318928099558 215, > 50979.498088990484 215, 50974.5448485677 215, 50969.59160814491 215, > 50963.0777577259 215, 50943.129090817689 215, 50940.482839084973 215, > 50940.482839084973 215, 50885.820734580833 215, 50862.136328915352 215, > 50840.782769521713 215, 50836.271454156857 215, 50829.654858288406 215, > 50821.008170505767 215, 50814.767517584383 215, 50803.113286225176 215, > 50779.203314791455 215, 50763.940031140366 215, 50751.308348118779 215, > 50723.714135803755 215, 50702.886896536016 215, 50677.02202177752 215, > 50659.72864621225 215, 50653.1120503438 215, 50649.954129588405 215, > 50646.495454475349 215, 50620.630579716853 215, 50588.073920500494 215, > 50574.840728763593 215, 50574.840728763593 215, 50554.915752568821 215, > 50548.223968110957 215, 50481.230934942883 215, 50474.614339074433 215, > 50425.365813008117 215, 50420.177800338533 215, 50405.440836813345 215, > 50400.252824143761 215, 50348.523074626777 215, 50333.259790975688 215, > 50317.9965073246 215, 50225.965673881583 215, 50209.424184210453 215, > 50191.228545572209 215, 50174.762244490492 215, 50126.1904157289 215, > 50126.1904157289 -100, 51020.548917151937 -100, 51020.548917151937 215, > 51020.548917151937 215)) and following is the exception which I was getting after creating the polygon in NTS: > found non-noded intersection between LINESTRING(51020.5489171519 215, > 51025.9092732259 215) and LINESTRING(51025.9092732259 215, 51011.7280780429 > 215) [ (51020.5489171519, 215, 53973.4510575744) ] Please help me figure out what is the problem with the polygon? Thank You


title: Installing knn on PostgreSQL for windows

text: I am trying to install the extension KNN for clustering as an extension on my PostgreSQL database. I have downloaded the file > http://www.pgxn.org/tag/analysis I found instructions on how to install the package on windows here > http://www.postgresonline.com/journal/archives/216-PostgreSQL-9.1-Exploring- > Extensions.html However the file i downloaded does not contain the file make or make install. Can anyone help. It doesn't need to be for this particular extension as im sure its the same for all extensions and it would be a nice bit of knowledge to have KNN is a clustering algorithm which groups data points based on their attributes into clusters or groups which are similar. My data set is a geographic data set with point information. Classic KNN doesn't treat information like latitude and longitude as spatial elements but as regular attributes (in this case numbers) The extension i am trying to install uses the logic of KNN but opens it up so that it also incorporates the spatial information i.e. where the data points are on a map of Europe to cluster the information. There are methods for doing this using Python and R but it is easiest to implement in PostgreSQL and then visualize in QGIS. An example was outlined in the thread **_Spatial clustering with PostGIS_** in these forums previously but it was for linux/unix systems and i wish to modify it so that i can incorporate it into my windows PC for my particular geographic dataset Thanks for your time


title: Mapinfo 10.5 and 8.5 compatibility? Unable to open TAB file

text: I have recently submitted some data to a client as a TAB file. They are using v8.5 and I am using v10.5, the data is just one .TAB (with the associated other files). And all the field types are character, so I thought this wouldn't be an issue. Are there any other things I may have overlooked when saving my file? Regards, Alan


title: How can I mosaic images in ArcMap without it affecting the colour?

text: This is 4 images (NE, SE, NW, SW) sections of imagery. They have some overlap between them.   * Working with 1m 3-band orthoimagery.    * Different parts are sometime flown at different times/sources hence the difference in times of the year.    * The images are in .tif format.    * Imagery properties says it's 8 bit pixel Depth and the pixel Type is unsigned integer ![enter image description here](http://i.stack.imgur.com/UrdZ6.jpg) When I mosaic them the colours change. **Looks the same as below** Number of bands = 3 Mosaic Operator = Last Mosaic Colormap Mode = Last **Looks the same as below** Mosaic Operator = Last Mosaic Colormap Mode = Match ![enter image description here](http://i.stack.imgur.com/zyeuX.jpg) Mosaic Operator = Last Mosaic Colormap Mode = First This image looks little better, features are better defined. Less Haze to. ![enter image description here](http://i.stack.imgur.com/KdzG0.jpg) How can I mosaic these images without it affecting the colour? The brightness of the image has increased.. Is there alternative software to Arcmap that merges raster's better for free? **EDIT** Images are 8 Colour bit. Updated better Mosaic image. Still using Raster Mosaic still alters the image colour. I used the Mosaic to New Raster tool. Is there a way to remove the black parts?


title: Polygons of Land Use Change

text: I am doing a project where I am digitizing a map of Haiti from the 70s. Through R2V I digitized the lines and then I exported it into ArcGIS 10. I am just interested in the land use change, therefore, I have deleted the other vectors. I want all of the land use change areas to be complete polygons. I performed through data management in the toolbox a feature class to polygon operation. This was somewhat successful. However, it did not identify all polygons. When I have complicated and nested polygons how do I get ArcGIS 10 to identify them as such? Your help would be much appreciated! Meg!


title: Highlighting Link in OpenLayers Map dynamically

text: I'm having a map which is rendered through OpenLayers using GeoServer. The map contains various attribute data field like link_no, X_coord, Y_coord etc. I want to highlight the link using link_no. I internally use WMS layer to render the map. I got the following code from this wiki for highlighting purpose:               // define a blank Vector layer and add it to the map     var highlight_style = { fillColor:'#99CCFF', strokeColor:'#3399FF', fillOpacity:0.7 };     hilites = new OpenLayers.Layer.Vector("Highlighted",                 {isBaseLayer:false, features:[], visibility:true, style:highlight_style}             );     map.addLayer(hilites);          // define the WFS server which will fill requests     var wfs_url = '/cgi-bin/mapserv?map=/maps/spraywatch2/wms/mapfile.map&SERVICE=WFS&VERSION=1.0.0';     /* highlightFeatures() makes the WFS request, then highlight_them() callback does the real work     * Args for WFS filter: typename, attribute, value     */     function highlightFeatures(typename,attribute,value) {     var wfsurl = wfs_url + '&REQUEST=getfeature&typename=' + typename +            '&Filter=<Filter><PropertyIsEqualTo><PropertyName>'+attribute+'</PropertyName><Literal>'+value+'</Literal></PropertyIsEqualTo></Filter>';     OpenLayers.loadURL(wfsurl,'',null,highlight_them);     }     function highlight_them(response) {                           // use the GML parser to turn the XML into a list of Feature objects                           var features = new OpenLayers.Format.GML().read(response.responseText);                           // have the Vector layer purge its feature list, replace them with the new ones                           hilites.destroyFeatures();                           hilites.addFeatures(features);                           hilites.setVisibility(true);                        }      Now what exactly is the typename, attribute and value parameters in the function highlightFeatures ? My aim is to highlight a particular link using a function call like this as explained in the wiki:               highlightFeatures('parcels','parcel_id','123456')      The above function call was an example given in the wiki link.


title: Using ESRI geoprocessing services from a non-ESRI client (Python)?

text: I'm thinking it should be possible to consume an ESRI Geoprocessing service from a python script without invoking ArcPy or needing an ESRI installation on the client machine. For example I create a geoprocessing service that takes a series of point location and returns a csv file with the values under them as a file. I would like to include this service call in a python client application that will be running on machines that don't have any ESRI installation. Couldn't you create the rest call to kick off the service and consume the result without arcpy? Does anyone have an example of how you might do this? Thanks, Colin


title: Should I use ISO 3166-1 Numeric codes as Primary Keys for Countries in a database table?

text: My main concern is, will one country ever adopt a numeric code that has been used for another in the past, such as can happen with alphanumeric 3166-1 codes? If not, then it seems like a perfect identifier for a primary key in a database...


title: So I am a programmer. What do I do to become a mapper?

text: I have now been rejected on 5 jobs asking for mappers on the pure basis that I am a programmer. I don't get interviews and when I ask I get response of "You are a programmer and we wanted a mapper with programming experience". I have my basics, I can use ArcMap, configure ArcSDE, understand coordinate systems. I have experience writing ESRI add-in code in ArcMap and direct code accessing the server without ArcMap. I understand how the spatial type is implemented and know how to create a matching UDT in accordance to OGC. I am familiar with the basic spatial operations. I know that I am missing some things (possibly a lot of things), but nobody wants to tell me what they are. So please tell me. Please assume some generality, I want to know what most programmers lack when compared to mappers, just listed my case to prove motive. Note: Open source GIS is coming next.


title: How to fix 'GEOSDifference: TopologyException'?

text: Why does Postgis throws an error in this SQL? Error:               ERROR:  GEOSDifference: TopologyException: found non-noded intersection between LINESTRING (-5.67027 43.5374, -5.67085 43.5376) and LINESTRING (-5.67111 43.5377, -5.6706 43.5375) at -5.6708422794605911 43.537631835802237      Query:               select st_astext(         st_difference(             st_geomfromtext('POLYGON((-5.66666 43.53301,-5.66704 43.53335,-5.66931 43.53374,-5.66975299919511 43.533380063154,-5.66979 43.53335,-5.66971 43.53478,-5.66969091189271 43.5348058613067,-5.67026768864541 43.5374211525033,-5.67084666354158 43.5376334432986,-5.67111 43.53773,-5.67059550578171 43.53754135212,-5.66854 43.53971,-5.66854 43.53971,-5.66854 43.53971,-5.66758 43.53976,-5.6665 43.5402,-5.6665 43.5402,-5.66609712188377 43.539781395014,-5.66401 43.54013,-5.66401 43.54013,-5.66351576708437 43.5392869195307,-5.66281 43.53909,-5.6626 43.53851,-5.66159162498075 43.5369268527413,-5.66131 43.53695,-5.66196 43.53617,-5.66272791899155 43.5356642562857,-5.66378432333548 43.5344801289574,-5.66384 43.53414,-5.66452 43.53344,-5.6656 43.53337,-5.66611 43.5331,-5.66666 43.53301))', 4326),             st_geomfromtext('POLYGON((-5.66577 43.53447,-5.66612 43.53455,-5.66675167677363 43.5349178897119,-5.66735 43.53529,-5.66811 43.5355,-5.66786 43.53678,-5.66816 43.53697,-5.66801 43.53717,-5.66729498033211 43.5375449497042,-5.6665 43.53782,-5.66634200990626 43.5386936730934,-5.66634 43.53876,-5.66554366174751 43.5384386684555,-5.66521 43.53856,-5.66457 43.53781,-5.6640744056249 43.5375662650614,-5.66341231504476 43.5372589777288,-5.66324 43.53718,-5.66334987968078 43.536344363177,-5.66304550704754 43.5357063513112,-5.66299 43.53559,-5.66358 43.53544,-5.66413 43.53529,-5.66521 43.535,-5.66577 43.53447,-5.66577 43.53447))', 4326)         )     );      This is a example of shapes with this problem. All share a common pattern. Some parts of edge are near to be parallel. ![example 1](http://i.stack.imgur.com/MIfpo.png) ![example 2](http://i.stack.imgur.com/7BWzU.png) ![example 3](http://i.stack.imgur.com/90XIZ.png) ![example 4](http://i.stack.imgur.com/bFoml.png) The point reported in the error is pointing to this location. In example 4, I marked this point with a circle. Versions:   * PostgreSQL 9.1.6 / x86_64   * POSTGIS 2.0.1 r9979   * GEOS 3.3.5-CAPI-1.7.5   * PROJ Rel. 4.7.1, 23 September 2009   * LIBXML 2.7.8   * Ubuntu Server 12.04.1 LTS x64 Thanks!


title: Viewshed Analysis With GDAL

text: I am looking to perform a simple (one observer location) viewshed analysis on a standard DEM using GDAL (either python or .NET bindings). I have searched, but haven't found any readily available implementation. Can anyone point me in the right direction?


title: What are Definition, Algorithms and Practical Solutions for Concave Hull?

text: ### Convex Hull A convex hull of a shape is defined as: > In mathematics, the convex hull or convex envelope for a set of points X in > a real vector space V is the minimal convex set containing X (Wikipedia) Wikipedia visualizes it nicely using a rubber band analogy (image below), and there are some good algorithms to compute it. ![alt text](http://i.imgur.com/XEelc.png) ### Concave Hull A concave hull is visualized using the red line in the image below (the blue line visualizes the convex hull). Intuitively, it is a polygon which embraces all the points, but has less (minimal?) area compared to the convex hull. As a result, the polygon's boundary length is longer. A concave hull may be the solution for some real-world problems (e.g. finding the reasonable boundary of a city). ![enter image description here](http://i.stack.imgur.com/hwxSW.jpg) I have failed to find a proper definition, algorithm and practical solution for the notion of a Concave Hull. The Grass Wiki has some descriptions and images, and there is a commercial solution in concavehull.com. Any ideas, algorithms and links will be very much appreciated. Adam


title: How to draw polygons around known clusters

text: I need to split up an area into non-overlapping, but contiguous territories. I have n points on a map, which are partitioned into k sets. The partition is only partially geographic, meaning that there is no reason to believe that the points are mixed in with one another. No two points overlap each other. The problem is to draw polygons such that each polygon contains all points from one set, and excludes all points from all other sets, such that the polygons have the following properties: No two polygons overlap each other. Polygons must be contiguous. A further "nice to have" would be that the polygons minimise the proportion of the area they cover which is not within r distance units of any of the points in the set, subject to the condition that no two lines of the polygon are closer than d to each other. I can't post an image of what I'm talking about until I have reputation 10. Please let me know if what I'm asking is unclear. My question is this: Has anything like this been done before? Are there any papers I should be reading, or known algorithms I can make use of? I'm new to GIS in general, are there any suggestions on search terms I could be using to find similiar concepts? ![http://i.stack.imgur.com/2kEEG.png](http://i.stack.imgur.com/XD1dN.png)


title: How to extract breaklines from a given DEM?

text: Breaklines are very useful for describing a terrain - for example to build an exact TIN. My question is, how to extract these breaklines from a given DEM? I calculated the curvature (vertical) - raster (using arcgis 10), where you can see the breaklines (min/max-curvature), but I did not figure out how to extract them to get polylines in the end. What is the workflow to get these breaklines from a DEM? Does anyone have an idea of a workflow, or a "tool", or skript, which can help? I think I am not the first/last man who wants to perform this task:) ...in the meanwhile i have found **a paper on this topic** that seems to be a good solution, but the big question remains: how to reduce that to practice? Or has anyone a different (easier?) approach with some software?


title: How to speed up band math in R 'raster' package?

text: I have some GeoTiff rasters with over 100 bands in total, but they are not huge: resolution is about 1200x600 and total size is less than 7 mb. I'm trying to do a simple NDVI calculation:               raster <- stack('path_to_raster_1', 'path_to_raster_2')     red <- 60 # red band number     NIR <- 75 # NIR band number     f_NDVI <- function(x){       (x[NIR]-x[red])/(x[NIR]+x[red])     }     NDVI <- calc(raster, f_NDVI)      The calculation runs **extremely** slow and takes about 10 minuets. I wonder if there is a way to optimise it? And why calculation is so slow in my case?


title: Shooting star algorithm in postgres 8.4 and postgis1.4 showing fatal error and remaining are working fine

text: When I try to run this shooting star algorithm, the postgis suddenly loosing connection to postgres and showing error no connection to server               SELECT * FROM shortest_path_shooting_star('        SELECT gid AS id,            start_id::int4 AS source,            end_id::int4 AS target,            shape_leng::float8 AS cost,           x1, y1, x2, y2,           rule, to_cost        FROM viz_network',        173,        97,        false,        false);      And digkstra algorithm is working fine :               SELECT * FROM shortest_path('       SELECT gid AS id,            start_id::int4 AS source,            end_id::int4 AS target,            shape_leng::float8 AS cost       FROM viz_network',       1,       310,       false,       false);


title: add features to a featurelayer and save changes in the server

text: I am developping a web application with arcgis javascript api. I am using arcgis server 10.2. I declared a featureLayer, the URL of this featureLayer is a feature service on my server. Using a geoprocessing service, I am able to create features from the results of the service. I am able to add the features to the featureLayer I declared and show them in the map, but the problem is the changes made are not saved on the server. this is the code sample:                function displayResults(results, messages) {         var flayer = new esri.layers.FeatureLayer("http://hamri-pc:6080/arcgis/rest/services/edt/FeatureServer/0", {           mode: esri.layers.FeatureLayer.MODE_ONDEMAND,           outFields: ['*']         });           var features = results[0].value.features;         var theExtent = null;         for (var f = 0, fl = features.length; f < fl; f++) {             var feature = features[f];             var attribs = feature.attributes;                 feature.setSymbol(normalpictureMarkerSymbol);                 feature.setInfoTemplate(infoTemplate);                                 flayer.add(feature);                           flayer.applyEdits(feature,null,null);                                 map.addLayers([flayer]);             }         }      Does anyone have an idea of what should I do to save the features in the server?


title: Resample, ArcGIS: How do I avoid resample taking into acount NoData values?

text: I am trying to **Resample** a raster dataset to match the resolution of another dataset. I have fine resolution data on small islands (surrounded by NoData in the sea). When I re-sample the data to a resolution that is coarser than the island size, I lose this data. I need a way to re-sample to a coarse resolution, that takes into account the data on small islands, and not the NoData which takes up the majority of the new pixel size. Any help? David