title: How to automate creation of toe of slope from DEM?

text: I was looking for an automated process to "identify" the toe of a slope from a DEM and ran across this ESRI thread and technique: **I am at Step 7 but have no idea how this is accomplished. I get it is in Raster Calculator but how???** > Re: Finding toe-of-slope in a DEM > > Author William Huber > > Date Nov 29, 2003 > > Message There are a lot of approaches that might be effective. Here's one. > Implement it using the Raster Calculator. > > (1) Create an indicator grid of "flat" areas. In the example (attached), I > have begun with a DEM of part of Highland County, VA. The "flat" areas are > those with a slope of less than 6 percent. > > (2) Creat an indicator grid of "steep" areas. This is the complement of the > flat areas. > > (3), (4) Use the indicator grids to mask the original DEM, creating grids > showing elevations only where the slope is flat or steep. > > (5), (6) Compute neighborhood means of the steep elevations and flat > elevations. > > _**_(7)Create a "valley edges" grid at cells where the mean of the steep > elevations exceeds the mean of the flat elevations.**** > > (8), (9) Simplify this valley edge grid by performing a Boundary Clean > operation, then shrink it by a cell or two. Any help would be appreciated. Or maybe there is a different way.


title: Is it possible to use If/Then Statements with dynamic text in Data Driven Pages?

text: Is it possible to use If/Then statements with dynamic text in Data Driven pages labeling? I have been trying a bunch of different ways, but it has not been working, and I am fairly new at writing scripts. If I was to write out my statement in a label expression, my example would be:               Function FindLabel([Land_Owner_Zone1], [Land_Owner_Zone2], [Zone])     If [zone] = "Zone 1" Then     Findlabel = [Land_Owner_Zone1] & vbnewline & [Zone]     Else     Findlabel = [Land_Owner_Zone2] & vbnewline & [Zone]     End If     End Function      I am unsure how to write this in dynamic text format for data driven pages.


title: Filtering MSSQL layers in QGIS doesnt work

text: In my QGIS Plugin I am filtering layers with this Python code:               FeatureIterator = layer.getFeatures( QgsFeatureRequest().setFilterExpression ("Fieldname='SomeValue'") )     for feature in FeatureIterator :         attrs = feature.attributes()             print "Feature ID: " + str(feature.id()) + " " + str(attrs)      This works for most layer-types (shapefile, Postgis, personal GDB,etc). But when I tried it with MSSQL Spatial layers, no results are returned. The iterator is empty. Edit: It works if I use the same filter clause (Fieldname='SomeValue') as a Feature subset at Layer-Properties. Then the layer shows only the intended subset. I work with QGIS 2.2 and 2.4. MSSQL Server is 2008R2. My testdata was exactly the same in Postgis and MSSQL, so I know that the filterexpression is correct. Could this be a bug? Edit: I was able to track the sql querys that are executed on the Database. Instead of the filterexpression that I intended, QGIS adds an Intersect clause and an ID filter to the query: select .... from ... where [shape].STIntersects([geometry]::STGeomFromText('POLYGON((0.00000000 0.00000000, 0.00000000 0.00000000, 0.00000000 0.00000000, 0.00000000 0.00000000, 0.00000000 0.00000000))',25832)) = 1 and [objectid] = 30681189080760422 The Coordinates of the intersect-rectangle are allways 0 and the objectid- value doesnt exist.


title: Network distance computations using GRASS

text: I'm an ArcGIS user who's interested in learning to use GRASS. Can anyone point out resources that will help me in learning how to obtain network distances between (1) residential locations (points) and (2) service locations (points also)? While I do wish to obtain the shortest paths, I need multiple paths between each residential location and a given number of services (say, the shortest paths to the 5 closest stores for each residential locations). Help and suggestions will be much appreciated. MP


title: Imported points from csv file (latitude, latitude) don't show up in a post code shape file

text: I can produce a point vector layer by importing latitude and longitude info from a csv file in a new project. After saving this file as an ESRI shape file I can import it in a Postal Code shape file. However now, the points don't show up. The same happens if I import the csv file immediately in the project with the postal code shape file. I probably do something elementary wrong. Can some one help me out. My aim is to show locations of interest on this postalcode file. Romée


title: Geoserver Feature Versioning Support

text: I am creating a map feature editor that supports versioning & rollback using Geoserver as the WFS. From what I've read online, it seems Geoserver supports feature versioning. My question is how do I implement versioning on Geoserver? Do I need to make any sort of modification on the data store for this to work or will the versioning be handled automatically by Geoserver? Would be great if there was a tutorial on this. Thanks.


How do you compute vertical transformations in QGIS? Ideally, the procedure / plugin will allow the loading of sundry geoid models.


title: Add Shapefile now fails :(

text: This code used to work...               if (File.Exists(@"C:\Roads.shp"))     {       mapControl.AddShapeFile(@"C:\", "Roads.shp");     }      I have Roads.shp on the C:\ drive I also have Roads dot {sbn, sbx, shp.gcd, shp.xml, shx} in the same place, c:\ The program throws an "The specified filename is invalid" exception with no Inner Exception. Could it be related to installing ArcEngine 10 SP1? That's the only thing I can think of... Really appreciate any ideas, thanks in advance. Using ArcEngine 10 with SP1 with C# VS2010 * * * I guess I was making a Noob mistake it looks like I forgot to copy Roads.dbf to the target directory. Too bad the error message isn't slightly more helpful.


title: Resolving ImportError: No module named arcgisscripting?

text: I am trying to run an ArcGIS python script from Eclipse and get the following message: > ImportError: No module named arcgisscripting I can run this in IDLE... any ideas?


title: Arcgis SDK for Android - Query multiple feature layers

text: I'm developing a android application with arcgis sdk. I have a basemap and then 4 different feature layers on top of it. I began with the arcgis sample "attribute editor" and removed the edit functioncs because I just want to view the features of the clicked point. I had a little issue with querying a point or a 1px line on the touchscreen but I solved it with creating an envelope around the clicked point and then query the envelope. So far so good! But my problem now is, that I can only query one feature layer at a time. I want to click on the map and then it queryies the layer 1, layer 2, layer 3, and finally layer 4. How can you do this? I was just able to query layer 1 OR 2 OR 3 OR 4... I have this code sections:               featureLayerDefault.selectFeatures(query,     ArcGISFeatureLayer.SELECTION_METHOD.NEW,     new CallbackListener<FeatureSet>() {         // handle any errors         public void onError(Throwable e) {         Log.d(TAG, "Select Features Error" + e.getLocalizedMessage());         }          if (queryResults.getGraphics().length > 0) {     Log.d(TAG, "Feature found id=" + queryResults.getGraphics()[0].getAttributeValue(featureLayerDefault.getObjectIdField()));     // set new data and notify adapter that data has changed     listAdapter.setFeatureSet(queryResults);     listAdapter.notifyDataSetChanged();


title: How to remove small "salt & pepper" polygons from a layer?

text: I have a map made up of multiple polygons, with a few hundred thousand ~~points~~ very small polygons ~~on~~ within it. I need the larger polygons, but not the ~~points~~ small polygons. How do I select and remove all the ~~points~~ small polygons, while preserving the larger polygons? (I obviously do not want to do this manually - there are a few too many points for that!) I'm totally new to QGIS, and I've Googled extensivey with no luck. Probably because I'm unsure what the correct terminology to use in a Google search is! The files I am using are a .dbf, .prj, .qpj, .shp and .shx. I'm simply dragging the .shp file on to QGIS, and the map appears as one layer. Basically a map of the UK divided into regions (which I want) and covered in points (which do not want). ![map view](http://i.stack.imgur.com/IUyjM.jpg) ![table view](http://i.stack.imgur.com/BZe0i.jpg) This is what it looks like up close. I did a dissolve to see if they'd go in to the surrounding layer, with no luck: ![close-up view of polygons](http://i.stack.imgur.com/4h2WE.jpg) Thanks!


title: How to create raster attribute table to a raster data using gdal python?

text: I classified a remote sensing data and saved it as a raster using Gdal python. I want to create a raster attribute table for this classified result.. my classified raster has five classes with 1,2,3,4,5 representing classes A,B,C,D,E respectively... Is it possible to assign color to each class.. since it is a single band raster image..


title: code for finding gaps between polygon parcels

text: i am trying to write code for finding gaps between polygons, the idea is to convert polygons to polylines and checking the attribute table, each selected segment that shows two records in the tables represent adjacent edge between polygons, if the segment shows only one record,it then can be considered as part of gap contour. hope someone can help me to convert this procedure to VBA code.


title: How to get a QGIS useable style file from OS .xml and .xsl files?

text: Ordnance survey provide some style files to support their vectormaplocal product. The style files appear to come in pairs, similarly named but one with a .xml extension and the other .xls The two files together seem to be designed to define the style but my problem is how do I get the two into a format that QGIS can handle?


title: How to set/save the extent of a map in QGIS?

text: I'm creating a bunch of maps with TimeManager in QGIS in order to create an animation. In order to overcome some of the limitations of TimeManager for dealing with raster files, I'm creating animations 30 frames at a time, which I can export and later combine the frames to make a longer animation. The problem with this approach is that it becomes important that the extent and scale of the map is identical for each of the 30 frame projects. So far, the only way I know to set the extent of a map in QGIS is to select the hand tool and drag it around, setting the scale numerically at the bottom of the map window. It would be far preferable if I could set the extent numerically as well, so that I could ensure that it is the same each time. Is there a way to do this?


title: Mysterious Problems with .IGE file

text: I'm trying to manipulate an .IGE file, which I understand is a raster data format developed by Erdas Imagine. The data I've downloaded is the National Land Cover Database (NLCD) 2006 Percent Developed Imperviousness data, available here: http://www.mrlc.gov/nlcd06_data.php I've tried using winzip and winrar to unzip the file, both of which bring up error messages during the unzipping process. Files do appear, however, in the folder I have designated for extracting the files, including an .IMG file. I can load this .IMG file into ArcMAP 10.1 and QGIS 1.8.0. In either of those programs, I can view the image, but I can't clip it to a shapefile, lest error messages appear. I have tried loading this .IMG file into Opticks 4.10.0 32-bit as well, but an error message comes up indicating the opticks cannot find the file it needs to import. I wonder if I didn't unzip the compressed file properly? Any suggestions will be greatly appreciated, thanks!


title: Display Labels in Google Maps from Spreadsheet Mapper

text: How would I go about displaying labels on placemarks on the map canvas where the source kml comes from the Spreadsheet Mapper? This particular use case is where we want to share and control the source spreadsheet in Google Docs, and have a display interface that is fairly ubiquitous, i.e. on thin blade client browsers, android browsers. Not either of these display devices will allow the use of Google Earth (thin blades don't have it installed cause of those well meaning network admins.... and droid GE is not able to display limited share files AFAIK). It is highly desired to be able to edit the source data in a GDocs SS, since the placemark locations will rarely change, but the labels and comments should pretty straightforward edits.


title: Wrap GEOS functionality into python

text: I've been looking on how to get GEOSExtractLine available in python. I looked into Shapely source to if it was available for it, but it's not. I'm trying to wrap this funcionality into Django, much like other linear referencing methods (https://code.djangoproject.com/ticket/11948). I changed the code properly but I'm getting a "GEOSExtractLine" not found error when trying to load the function. Here's the traceback:                 File "C:\Python27\lib\site-packages\django\contrib\gis\geos\__init__.py", line 6, in <module>         from django.contrib.gis.geos.geometry import GEOSGeometry, wkt_regex, hex_regex   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\geometry.py", line 15, in <module>         from django.contrib.gis.geos.coordseq import GEOSCoordSeq   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\coordseq.py", line 10, in <module>         from django.contrib.gis.geos import prototypes as capi   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\prototypes\__init__.py", line 13, in <module>         from django.contrib.gis.geos.prototypes.geom import from_hex, from_wkb, from_wkt, \   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\prototypes\geom.py", line 134, in <module>         geos_extract_line = geom_output(GEOSFunc('GEOSExtractLine'),[c_double,c_double])   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\prototypes\threadsafe.py", line 39, in __init__         self.cfunc = getattr(lgeos, func_name)   File "C:\Python27\lib\ctypes\__init__.py", line 366, in __getattr__         func = self.__getitem__(name)   File "C:\Python27\lib\ctypes\__init__.py", line 371, in __getitem__         func = self._FuncPtr((name_or_ordinal, self)) AttributeError: function 'GEOSExtractLine' not found      How do I find out the correct name that GEOS publishes for this method? I tried to follow the convention used by Shapely and Django's GEOS wrapper, found here http://geos.osgeo.org/doxygen/classgeos_1_1linearref_1_1LengthIndexedLine.html#a20e3c3e6faa7183e23e3f9a312893aca What I'm truly after is this:               def test_extract_line(self):         line = LineString((0,0),(10,0))         extract = line.extract_line(.5,1)              from_point = Point(5,0)         to_point = Point(10,0)              self.assertEqual(extract[0][0],from_point[0])         self.assertEqual(extract[0][1],from_point[1])         self.assertEqual(extract[1][0],to_point[0])         self.assertEqual(extract[1][1],to_point[1])      Any tips? Sean Gillies come to rescue? Help! :D EDIT: based on this ticket https://code.djangoproject.com/ticket/11948 and the patches it describes I've added the following code to the respective files: geometry.py (the actual wrapper):               def extract_line(self,from_distance,to_distance):         if not isinstance(self,(LineString,MultiLineString)):             raise TypeError('extract line only works on LineString and MultiLineString geometries')         return capi.geos_extract_line(from_distance,to_distance)      geom.py (where I point to geos ctype library)               geos_extract_line = geom_output(GEOSFunc('GEOSExtractLine'),[c_double,c_double])      __init__.py I've added the import instruction for geos_extract_line. When I do this, I get the traceback written above "function GEOSEXtractLine not found". Digging a bit further on GEOS docs, it looks like the correct function I'm after it's called Extract Line By Location (http://geos.osgeo.org/doxygen/classgeos_1_1linearref_1_1ExtractLineByLocation.html) Lastest EDIT: My guess is that this function is not published by GEOS capi, because:               >>> from django.contrib.gis.geos import libgeos     >>> libgeos.lgeos     <CDLL 'C:\OSGeo4W\bin\geos_c.dll', handle 4200000 at 2ad9870>     >>> dir(libgeos.lgeos)     ['GEOSArea_r', 'GEOSBoundary_r', 'GEOSBuffer_r', 'GEOSContains_r', 'GEOSConvexHull_r', 'GEOSCoordSeq_clone_r', 'GEOSCoordSeq_create_r', 'GEOSCoordSeq_getDimensions_r', 'GEOSCoordSeq_getOrdinate_r', 'GEOSCoordSeq_getSize_r', 'GEOSCoordSeq_getX_r', 'GEOSCoordSeq_getY_r', 'GEOSCoordSeq_getZ_r', 'GEOSCoordSeq_setOrdinate_r', 'GEOSCoordSeq_setX_r', 'GEOSCoordSeq_setY_r', 'GEOSCoordSeq_setZ_r', 'GEOSCrosses_r', 'GEOSDifference_r', 'GEOSDisjoint_r', 'GEOSDistance_r', 'GEOSEnvelope_r', 'GEOSEqualsExact_r', 'GEOSEquals_r', 'GEOSFree_r', 'GEOSGeomFromHEX_buf_r', 'GEOSGeomFromWKB_buf_r', 'GEOSGeomFromWKT_r', 'GEOSGeomToHEX_buf_r', 'GEOSGeomToWKB_buf_r', 'GEOSGeomToWKT_r', 'GEOSGeomTypeId_r', 'GEOSGeomType_r', 'GEOSGeom_clone_r', 'GEOSGeom_createCollection_r', 'GEOSGeom_createLineString_r', 'GEOSGeom_createLinearRing_r', 'GEOSGeom_createPoint_r', 'GEOSGeom_createPolygon_r', 'GEOSGeom_destroy_r', 'GEOSGeom_getCoordSeq_r', 'GEOSGeom_getDimensions_r', 'GEOSGetCentroid_r', 'GEOSGetExteriorRing_r', 'GEOSGetGeometryN_r', 'GEOSGetInteriorRingN_r', 'GEOSGetNumCoordinates_r', 'GEOSGetNumGeometries_r', 'GEOSGetNumInteriorRings_r', 'GEOSGetSRID_r', 'GEOSHasZ_r', 'GEOSInterpolateNormalized_r', 'GEOSInterpolate_r', 'GEOSIntersection_r', 'GEOSIntersects_r', 'GEOSLength_r', 'GEOSLineMerge_r', 'GEOSNormalize_r', 'GEOSOverlaps_r', 'GEOSPointOnSurface_r', 'GEOSPrepare_r', 'GEOSPreparedContainsProperly_r', 'GEOSPreparedContains_r', 'GEOSPreparedCovers_r', 'GEOSPreparedGeom_destroy_r', 'GEOSPreparedIntersects_r', 'GEOSProjectNormalized_r', 'GEOSProject_r', 'GEOSRelatePattern_r', 'GEOSRelate_r', 'GEOSSetSRID_r', 'GEOSSimplify_r', 'GEOSSymDifference_r', 'GEOSTopologyPreserveSimplify_r', 'GEOSTouches_r', 'GEOSUnionCascaded_r', 'GEOSUnion_r', 'GEOSWKBReader_create_r', 'GEOSWKBReader_destroy_r', 'GEOSWKBReader_readHEX_r', 'GEOSWKBReader_read_r', 'GEOSWKBWriter_create_r', 'GEOSWKBWriter_destroy_r', 'GEOSWKBWriter_getByteOrder_r', 'GEOSWKBWriter_getIncludeSRID_r', 'GEOSWKBWriter_getOutputDimension_r', 'GEOSWKBWriter_setByteOrder_r', 'GEOSWKBWriter_setIncludeSRID_r', 'GEOSWKBWriter_setOutputDimension_r', 'GEOSWKBWriter_writeHEX_r', 'GEOSWKBWriter_write_r', 'GEOSWKTReader_create_r', 'GEOSWKTReader_destroy_r', 'GEOSWKTReader_read_r', 'GEOSWKTWriter_create_r', 'GEOSWKTWriter_destroy_r', 'GEOSWKTWriter_write_r', 'GEOSWithin_r', 'GEOSisEmpty_r', 'GEOSisRing_r', 'GEOSisSimple_r', 'GEOSisValidReason_r', 'GEOSisValid_r', 'GEOSversion', '_FuncPtr', '__class__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattr__', '__getattribute__', '__getitem__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_func_flags_', '_func_restype_', '_handle', '_name', 'finishGEOS_r', 'initGEOS_r']      With this I've got a glimpse of all the functions published by GEOS API and there is no Extract, ExtractLine or extractByIndex, or anything. In fact, if you look at CAPI (the published API) there is truly no reference of these functions: http://trac.osgeo.org/geos/browser/trunk/capi/geos_c.h.in Still, if anyone can help, I appreciate. Another edit: PostGIS has the funcionality. How did they do it? http://postgis.net/docs/ST_Line_Substring.html Thanks!


title: ArcObjects - Listening for Map Click in ArcMap

text: I have a seemingly simple requirement that is causing me all types of headaches - listening for a map click event in ArcMap via ArcObjects. Simple as that, just tell me when the map is clicked (and x,y location obviously). After some research, it appears that creating a custom tool that implements ITool was the way to go. Well it turns out that approach gets ugly due to some 'bugs' in ArcGIS Desktop (v10 at least). I don't want to add my custom tool to a toolbar, as it doesn't do anything and is irrelevant for a traditional ArcMap user. So, I thought I would just create a 'hidden' toolbar containing my custom tool. Well it turns out hiding a toolbar doesn't work so well (see this link). The only other thing I can think of to do is add my custom tool to a 'standard' toolbar. The problem with this approach is that if the user then hides or removes whatever standard toolbar I put the custom tool in from ArcMap, the tool is no longer accessible. Can someone please tell me that I am missing some very easy way to listen for a map click event without all the aforementioned complications?


title: How to use ArcGIS JavaScript API for ArcGIS Server service, where to save code?

text: I have found some sample codes online for implementing traditional map components such as legends and scale bars to an ArcGIS Server service. https://developers.arcgis.com/javascript/jssamples/index.html#map Could someone explain how this is implemented into the service? Where to save the code?


title: How can I show points on the same coordinate as e.g. a bigger symbol in QGIS?

text: I use QGIS 1.8 for mapping archaeological findings. Some of the findings within one layer are located on the same site and therefore have the same coordinates. Therefore you just see one point instead of the actually overlapping points. There is the tool - in German its called "Punktverdrängung" (maybe point displacement in english?), with which you can circle overlapping around a central point, but that does not really fit my needs. Is there the possibility of showing more than one point on the same coordinates as a bigger symbol or anything like that? Thanks for your help!


title: Empirical Bayesian kriging - does data need to be normal?

text: With ebk, does your data need to be normally distributed? If I apply a normal score transformation, my data is still not normal. Is ebk kriging thus an invalid interpolation technique for my data, or can I get around this?


title: How to join 2D polylines and DEM into 3D polylines?

text: I got a simple question but could not find a satisfactory answer yet: I have a road network in 2D/PostGIS and a DEM. I want to turn the road network into 3D by retrieving the Z-value from the DEM. Both files are quiet large so an efficient process is needed. I'd prefer to use QGIS or any other open source tool. thanks very much


title: Our best table-analysis with GIS?

text: Sometimes I think about possible analysis using GIS models (I use QGIS 1.8), and I always think about few types of visualization techniques: demographics, traffic flows, distance ones... Thanks to you and reading those pages I discover many possible criteria (and I begin to use it) analysis and I like to collect the most important one in a single item to help the consultation. Is this a good idea or not? I try with some examples of a project where I worked on ![generated traffic by zone](http://i.stack.imgur.com/GMtOR.jpg) generated traffic by zone by magnitude ![traffic by zone](http://i.stack.imgur.com/Xwz7m.jpg) traffic IN/OUT zone by zone by bar diagram ![journey time on different routes](http://i.stack.imgur.com/fc3wz.jpg) journey time on different routes ![OD from a zone](http://i.stack.imgur.com/QtdLd.jpg) OD from a zone (thanks to Flow Mapper plugin my question on it ) ![quantity and quality parks on links](http://i.stack.imgur.com/lONln.jpg) quantity and quality parks on links


title: Using all dimensions with postgis functions

text: I have points as geometry format in table that I use for some analysis. All points are saved in four dimensions eg. x, y, z and m - lon, lat, alt and time. The issue is that for the analysis are used functions like ST_GeomFromEWKT and ST_MakeLine, which included only two dimensions form the four present. My specific questions is can I have the four dimensions using these functions? I tried simple example               SELECT ST_AsText(ST_MakeLine(ST_MakePoint(1, 1, 0, 1339590840), ST_MakePoint(2, 2, 0, 1339590840)));      Which returns "LINESTRING(1 1,2 2)", without z and time dimensions. Thanks in advance!


title: How do I process xy data to be compatible with census block data (so I can include them in my GWR model)? (ArcGIS 10)

text: I'm pretty new to GIS, so apologies if my problem is rather simplistic I'm planning to do regression analysis (GWR) on several datasets, where all but one is based on census areas (UK OS lower level super output areas, lsoa) for which I have a shape file that describes the polygons. Regression with the census-based variables is easy to set up and execute because both the dependent and explanatory variables are held within the same Input feature. But I'm finding difficulty with my final explanatory variable - concentrations at x,y coordinates (UK OS co-ordinates). I believe I need to 'join' or 'relate' with my lsoa shapefile (which is what I've done with the other variables); I've tried all sorts of things such as kriging and extracting to vector, 'relates', 'joins', adding new fields, calculating new fields, but to no avail. I'm sure there is someone out there who has done this before - I'd really appreciate your advice. thanks Jane


title: OpenLayers: How to create a rectangle?

text: I want to use OpenLayers to draw a simple rectangle (based on 4 coordinates) on my map. I would also like to fill the rectangle with a color. For example (EPSG:4326):   * 13.40, 52.50   * 13.50, 52.50   * 13.50, 52.60   * 13.40, 52.60 The documentation tells me to use this:               var poly = new OpenLayers.Bounds(0,0,10,10).toGeometry();      How I am supposed to convert this?


title: How to test list elements in FME

text: My question is how to test each element in one list with each element in another list? I have lot of lines segments and using lengthCalculator to find all segments less than some specific length. Further with help of TopologyBuilder I get 'to_node' attribute related to LINE output port and 'node_number' related to NODE output port. I have to detect all dangles. So I am trying to test if there is a "node_number =! to_node" If so, then I know at this segment is dangle. Agree? for this I used transform Tester again, but this is not working. My next idea is to create list of all nodes and another list of all segments < specific_length, and then check (test) each node element in one list with to_node element in teh second list. Will this work? Not good idea? This question is linked to: Find dangling lines . Tnx


title: How to allow for multiple lines in legend?

text: How can I force the description text in my legend to occupy more than one line? This is basically the same question as asked before. However, I want to use _Categories/Unique values_ and not _Features_ as my Symbology style, in which case there's no Description button. Even when I change to _Features_ , it seems I can't use multiple lines (I suppose that there's been some changes to this functionality). How would I go about getting the description/header text across multiple lines?


title: How to fix "Raster layer: Provider is not valid" error when loading PostGIS rasters?

text: I compiled and installed the current QGIS master as usual and I imported some raster data into my postgres 9.3 database with postgis 2.0.2 extension. Whenever I try to import the raster data into qgis using the database manager, I get this message:               GDAL provider: Cannot get GDAL raster band:      Raster layer: Provider is not valid (provider: gdal, URI: PG: dbname=gis host=127.0.0.1 user=xxx password=xxx port=5432 mode=2 schema=public table=elevation      In the command line, I get               ERROR 1: Computed PostGIS Raster dimension is invalid. You've probably specified unappropriate resolution.      I have actually no idea where to search for the error, so any help is highly appreciated. * * * PS: Just moved to Postgis 2.1.2 - same problem there...


title: OpenLayers.Format.SLD - Ignores 'matchCase'?

text: I've trying to use OpenLayers to construct an SLD element to include in a WMS layer's parameters. This SLD document can contain a number of filters. This works fine and resultant SLD does work to filter the layer. However, the `matchCase=false` property of the filter(s) are being stripped out on formatting the SLD. Example:               var filter    = new OpenLayers.Filter.Comparison(     {       type:       OpenLayers.Filter.Comparison.EQUAL_TO,       property:   'my_property',       value:      'A Value',       matchCase:  false     });          var rule = new OpenLayers.Rule({       filter: filter,       symbolizer: {"Polygon": {}}     });          var style = new OpenLayers.Style();          style.addRules([rule]);          // Create the SLD     var sld_body = new OpenLayers.Format.SLD().write({       namedLayers: [         {           name: 'myLayer',           userStyles: [style]         }       ]     });      I think I know what the problem is, but I have no idea how to fix it. If I check the output of `OpenLayers.Format.Filter({version: "1.1.0"}).write( filter );` then I see correctly formatted XML, which includes my `matchedCase` property. If I do the same with version 1.0 (the default), I see the same XML but **without** the `matchedCase` parameter. The documentation for `OpenLayers.Format.SLD` states that it inherits from `OpenLayers.Format.Filter.v1_0_0`, which explains why it won't read the `matchedCase` parameter. I'd hoped using `OpenLayers.Format.SLD({version: "1.1.0}).write()` would work, but it just throws an error about a driver not being available for that version. So, how should I deal with this? I need to generate an SLD document and it needs to acknowledge the matchedCase property. Thanks, Tom


title: GDAL and Python: how to create new files for (e.g.) clipping with gdalwarp

text: I'd like to learn using GDAL, and I want to use it in Python. Following the tutorial (http://www.gdal.org/gdal_tutorial.html), I succeded in opening vector and raster files using ogr.Open and gdal.Open () calls respectively. I then tried to understand how to create new files, but I didn't figure it out so far. For instance: if I want to use "gdalwarp" function to (let's say) clip a raster with a polygon, it does its job perfectly and clearly within QGIS, but using GDAL in Python how do I: 1) specify the output folder? 2) create the new file (I think it's done using the Create or CreateCopy () methods, is it?)? I hope someone could help me, thanks in advance to the community. UPDATE I came up with creating a file in a specific folder. Anyway, it is just a file to be filled, namely, now I want to know how can I do my calculation and output to this "empty" file. Here is my code to let you understand me:               from subprocess import call     call(["ls", "-l"])     import gdal, ogr     from gdalconst import *          input= 'C:/Users/unimi/Desktop/sept2000_B5_refl.tif'     clip_area= 'C:/Users/unimi/Desktop/study_area_foscagno.shp'     output= 'C:/Users/unimi/Desktop/clip_sept2000_B5_refl.tif'          shp= ogr.Open (clip_area)     if shp is None:         raise SystemExit("The shapefile could not opened")     raster= gdal.Open (input, GA_ReadOnly)     if raster is None:         raise SystemExit("The raster could not opened")     print raster.GetProjection ()     format= 'GTiff'     driver= gdal.GetDriverByName (format)     metadata= driver.GetMetadata ()     dst_ds= driver.Create (output, 2133, 1273, 1, GDT_Float32)          ClipByPolygon= 'gdalwarp -dstnodata 0 -q -cutline %s -crop_to_cutline -of GTiff %s %s' % (shp, raster, dst_ds)     call (ClipByPolygon)          dst_ds= None      The code doesn't raise any error. However, the "empty" file from the "Create" call is produced, but gdalwarp has no effects afterwards.


title: PostGIS extension not available in Postgresql list

text: I have installed Postgresql 9.2 from graphical installer in Ubuntu 13.10 and Install PostGIS 2.0.3 through StackBuilder. All things went smooth. But it is not possible to create PostGIS database as there is no extension present in Extension. As far I know two more extension is needed to be there named 'postgis' and 'postgis_topology'. I have tried `CREATE EXTENSION postgis` but it returns there is no postgis.control file within the extension folder. What should I do now?


title: Transit Directions: Preferred Vehicle Type?

text: I would like to be able to get transit directions through the JS api, that is limited to a single vehicle type (e.g. Bus). This functionality is available on maps.google.com, but how do i submit the vehicle type as part of the DirectionsRequest?. Using the code below will provide routes via Bus, Train               google.maps.DirectionsTravelMode.TRANSIT      But if i want to restrict only to Train, how do i do it?


title: Coordinate System Google Maps API

text: I am making an application in java using GMaps4JSF where the person will enter the address on the map I'll take this coordinate and consult on my postgis base regions in which region this point is inserted and return this information on the screen. It's very simple. Google Maps Api uses the Mercator projection (EPSG: 3857), the coordinate of the queried map will be returned in this projection? I must redesign my data to it does not need to do conversions.


title: GeoMedia WebMap with WFS 1.1.0 Filter Examples?

text: I'm having problems forming a GET query for WebMap spatial data using WFS 1.1.0. My query looks like this:               http://USGVDCALIX2/WFS/Request.aspx?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&MAXFEATURES=500&FILTER=<Filter xmlns="http://www.opengis.net/ogc"><BBOX><PropertyName>Geometry</PropertyName><Box srsName="EPSG:4326" xmlns="http://www.opengis.net/gml"><Envelope srsName="EPSG:4326"><lowerCorner>-72.742129821777 43.944567565918</lowerCorner><upperCorner>72.545288391113 44.148836975098</upperCorner></Envelope></Box></BBOX></Filter>&TYPENAME=Counties      But the result is:               <ExceptionReport xmlns="http://www.opengis.net/ows" version="1.1.0"><Exception exceptionCode="GWMWFS.GeometryUtility" locator="GWMWFS.GeometryUtility"><ExceptionText>Unsupported geometry type.</ExceptionText></Exception></ExceptionReport>      According to this GeoMedia WebMap Professional 6.1 is WFS 1.1.0 compliant. And according to this my query looks valid. So does this error suggest that WebMap is unable to do a BBOX query against a polygon? Should the query be formed differently? Links to WebMap WFS 1.1.0 filter examples would be welcome!


title: Thematic mapping

text: What is the best way to develop thematic maps? I would like to generate pie charts and fill polygons with different colors depending on a parameter value. I'm working with ExtJS 4.2 and OpenLayers.


title: How to show two different labels of a particular vector layer with the help of QGIS?

text: I am working under a government project in Kolkata with QGIS. I want to show two different labels of a particular vector layer. just like we can show two different labels in Arc GIS with the help of Append option. Suppose I want to show Male population Numbers and Female population numbers on a District vector layer, so how it can be possible? please tell me the option of this Software (QGIS). Actually i am able to show one label on a layer but how to show two labels?


title: Sample FindNearest ArcGIS REST URL

text: Could anyone provide me with a sample FindNearest URL using the REST API for Arc GIS? It should look something like this: server.arcgis.com/FindNearest/submitJob/, but I have not been able to get a working link. The documentation provides no sample, so I am wondering if anyone here could provide me with a link, or explain how it is structured. I would greatly appreciate any help.


title: hypsometric colors alà IMHOF, SYDOW or PEUCKNER

text: I need a hypsometric color-map to add it to my multishade-modell. (High Mountain between 800 and 4000m) There are a lot of differend typs explained in IMHOFS "Cartographic Relief Presentation". I want to do it traditionally... :-) So, where can I find "prebuild" colormaps alà PEUCKNER, IMHOF or SYDOW? (I dont like the "Terrain" or "SRTM" ones from Grass) Thanks so far! Martin


title: Passing Python List into Layer Definition Query

text: I am trying to cursor search through items in a field in a feature class. I then append them to a list called 'distList'. I want to pass the list into layers definition query as follows:               "Field1" in distList      or in long hand, it would look like this:               "Field1" in ('ds(c1)', 'ds(b1)', 'ds(c2)', 'ds(g1)')      The list looks like this:               ['ds(c1)', 'ds(b1)', 'ds(c2)', 'ds(g1)']      The problem I am having is, the definition query above doesn't like to look 'in' a python list. It doesn't recognize the square brackets. It prefers round brackets, so I thought about using a tuple. The problem I am facing is I don't know how to build a tuple after cursor searching the items from the feature class. here's my code example:               distList = []     for lyr in arcpy.mapping.ListLayers(mxd, "", df):         if lyr.name == "Disturbance":                for row in arcpy.SearchCursor(lyr):                     if "ds" in row.Field1:                          distList.append(row.Field1)         lyr.definitionQuery = '"Field1"' + "in " + distList      Can anyone suggest a way to get my list into a tuple or maybe just a better way to get my items into a format where they have round brackets instead of square ones? As a work around, I tried converting the list to string `str(distList)` and then replacing the brackets. That works fine, but it seems cumbersome and I'm sure there's a better way. Thanks, Mike


title: How to implement Landmarks in Network Analyst of ArcGIS 10.1?

text: I have been trying to implement the landmark functionality into route solving in Network Analyst using ArcGIS 10.1. I have added a point class as turn landmarks (as well as playing around with confirmation landmarks as well), verified they have a valid name field, and adjusted the tolerance, yet the landmarks are not showing up in the directions. Am I missing something? Thanks, Eric


title: Is it possible to create a new ArcGIS toolbar with new tools using Python?

text: I have written a number of fairly-simple python functions using `arcpy` which help a user significantly when performing a certain workflow using various tools (both ArcGIS toolbox tools and external software). I'd really like to bind some of these functions to buttons on a new ArcGIS toolbar, as they are not really geoprocessing tools, but more like 'Workflow Management Tools', as they do things like show and hide various layers, change symbologies etc, and generally make it easier for the user to perform specific tasks. Is there a way to do this? I have looked in the documentation and have found various things about how to create toolbar buttons using ArcObjects, and I have found this question about accessing ArcObjects from Python, but I'm not sure how easily these can be combined. Does anyone know if this can be done, and if so, how?


title: Translating all features in a shape file

text: I have a friend with a bit of a problem. He has 1000+ shape files mapping paddocks and being used for precision agriculture. He's recently discovered the base station positions used to create these maps were incorrect. The base stations are being updated, which means each shape file needs an x and y offset applied to all features. The offsets vary but are typically small (~1-2 meters). The correction needs to be applied with sub-centimeter accuracy. I've been brought on board to automate the process. I've written a C# application that can accurately calculate the offsets to apply to each file. What I'm looking for is something to apply these offsets for me. Either a C# library or an external command line application would be ideal. I could implement something myself but I don't want to re-invent the wheel. Does anyone have a suggestion for achieving this? I've read that I can use ogr2ogr and specify a different false easting/northing for the source srs, although I'm not sure if this will provide the accuracy I require. I also don't know what projection the files are using. Any help would be greatly appreciated!


title: spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelWithin; not working

text: I have a task of getting the total area of all the polygon features within a boundary. So I got the envelope of the boundary and did a spatial query to get the features inside that element but my problem is I am able to get the elements that overlaps the envelope but not exactly inside it.. my code goes like this               IMxDocument mxDoc = (IMxDocument)GetMxDocumentFromArcMap(m_application);             IEnvelope gridEnvelope = GetGridExtents("Z29");             IFeatureClass featureClass = GetFeatureClassFromShapefileOnDisk("D:\\Projects\\Working\\BMRDA\\data", "Hoskote_Satellite_Polygon");                  // Create the spatial filter and set its spatial constraints.             ISpatialFilter spatialFilter = new SpatialFilterClass();             spatialFilter.Geometry = gridEnvelope.Envelope;             spatialFilter.GeometryField = featureClass.ShapeFieldName;             spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelWithin;             spatialFilter.set_OutputSpatialReference(featureClass.ShapeFieldName, mxDoc.FocusMap.SpatialReference);                  // Set the attribute constraints and subfields.             // We want to exclude ramps, highways and interstates.             spatialFilter.WhereClause = "LEVEL3 = 'AGRICULTURE LAND'";                  // Execute the query.             IFeatureCursor featureCursor = featureClass.Search(spatialFilter, true);             ESRI.ArcGIS.Geodatabase.IFeature feature;  // Automatically initialized to null. Used to test existence of a feature in the featureCursor             double dbArea = 0.0;             int intCount = 0;             while ((feature = featureCursor.NextFeature()) != null)             {                 IArea area = (IArea)feature.Shape;                 dbArea = dbArea + area.Area;                 intCount++;             }          the spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelWithin; does not work it returns me null I want the polygon to be clipped and get the area of it Am i missing something here? Please guide..


title: Change the image size in ARCGIS

text: I did interpolation operation on temperature and precipitation data. Output image size is 3300 (col) and 2550 (row). I have some NDVI images with size of 15456 (col) and 5074 (row). Generally, I want the interpolated images to be the same size with NDVI images. After doing interpolation, I wanted to do resampling on temp/prec images in order to not loss any information on image (3300: 2550 to 15456:5074 size). But resampling method only works on raster image. Is there any tool or procedure in ARCGIS to set a image size when I'm interpolating temperature and precipitation data. Do you have any suggestion to achieve this? Thanks a lot


I wanting to reproject the vertices of my polygon. I have been able to achieve this for the bounds but having a bit of trouble with the vertices.                          function onFeatureModified(event) {                var bounds = event.feature.geometry.getBounds();                var answer = "bottom: " + bounds.bottom + "\n";                answer += "left: " + bounds.left + "\n";                answer += "right: " + bounds.right + "\n";                answer += "top: " + bounds.top + "\n";                alert(answer);                     var mapProjv =  map.getProjectionObject();                var epsg4326v = new OpenLayers.Projection("EPSG:4326");            var verticesNative = event.features[0].geometry.getVertices();            var verticesLatLon = verticesNative.transform( mapProjv , epsg4326v  );            var vertices = verticesNative;                alert(verticesNative);                }      Below code that works:                          var vertices = event.features[0].geometry.getVertices();                alert(vertices);


title: What are the differences between the OGC and ISO WKT/WKB specification?

text: Are the OGC-Simple feature access and the ISO (ISO/IEC 13249-3:2011 "Information technology - Database languages - SQL multimedia and application packages - Part 3: Spatial") significantly different ( _different enought to justify purchasing the ISO/IEC 13249-3:2011 PDF_ )?


title: Error while trying to upload a .bsq file to the geoserver

text: I am trying to upload a .bsq file (ENVIHdr Coverage Format) to the geoserver. While doing this I got an error message. I am using the "GDAL" Extension for the geoserver. The error message is: "Could not list layers for this store, an error occurred retrieving them: Incorrect input type!"      Do I use the wrong extension or is there an other mechanism to upload a .bsq file to the geoserver ? After the upload I want to publish this raster data as a WCS service.


title: How to convert GRIB-1 to GRIB-2?

text: I need to get data for my database from some GRIB files (binary), I downloaded wgrib2, which can export GRIB-2 files to .csv and csv is very suitable for me, but some GRIB files are GRIB-1, so I used wgrib2 on them and program wgrib2 told me: "grib1 message ignored (use wgrib)", so I used wgrib, but wgrib can't export files to different format. So I tried to install recommended **cnvgrib** to convert GRIB-1 to GRIB-2, but still getting some errors during installation due different library versions, 32/64 bit problem, different directories etc., it's written in Fortran, that's also the problem I can't compile it easily... So is there another way to convert GRIB-1 files to GRIB-2 (or to some non- binary format)? I can google many ways to do it, but I need help from some1 who have already converted GRIB-1 to another format, because there can be some easy way how to do that, but I have so much pain fixing installation problems, wrong configurations, Fortran version incompatibilities etc during installation of cvngrib, which I haven't even fixed and gave up... EDIT: OK, my co-employee fixed degrib (it linked own libraries with newer system libraries), so I can run degrib now. So I run "degrib "GRIB file" -C -msg all -nMet -Csv" on GRIB-1 file and it failed the conversion to .csv with:               ERROR: In call to Grib2Convert.      ERROR: In call to ReadGrib2Record.      Problems with ReadGrib1Record called by ReadGrib2Record      Inside ReadGrib1Record      Grid projection number is 4      Don't know how to handle this grid projection.  I try to fix it soon, but meanwhile you can help me again!


title: GeoJSON support in ArcObjects 10 .NET

text: i am currently trying to build a geoprocessing application using OpenLayers and ArcObjects, I have some features that user draws on openlayers and i need to send these features to arcobjects for further processing.so i need a common feature format between these two systems. ## solution 1 first i have tried getting esri-json support on openlayers with agsjsadapter but i got some unhandled exceptions - still trying to succeed, any hints will be really appreciated. ## solution 2 currently i am trying to convert a geojson format feature to Igeometry or Ipolygon on arcobjects, but i couldn't find any support for this in arcobjects. i guess i need to solve this issue with one of the solutions above. and i am open to any kind of solution thanks already


title: "Field Mapping" in ArcGIS 10 - ArcPy

text: I've wrote a **Python** script that does a spatial join and some simple calculations. My problem is with setting the merge rule for one specific field, and leaving the rest of the fields as is. For example, I have a population field that when joined by spatial location uses the merge rule _"First"_ which grabs the first occurring of a Population count. I want to be able to set the merge rule to _"Sum"_ to **add up the population values between all of the polygons found in the spatial extent of another polygon**. I have done some intense tinkering with field maps and field mapping objects, but can't seem to get it working properly. Specifically I have tried the method: _popFieldMap.mergeRule = 'Sum'_ to set the mergeRule, but it always reverts to "First". Any ideas how I can programmatically change the merge rule for one field in a Spatial Join? Thanks! **Here is my code (keep in mind that it is pretty specific to my data and contains lines to test certain stages of the script):**               import arcpy,sys,os          #Get the Files involved, set some variables.     SectorTable = sys.argv[1]     SectorShape = sys.argv[2]     MaxDev = sys.argv[3]     PopulationFC = sys.argv[4]     OutputFC = sys.argv[5]     DeviationField ="Angle_Deviation"     ID = "SectorID"     newID = "BP_ID"     mxd = arcpy.mapping.MapDocument('CURRENT')     df = arcpy.mapping.ListDataFrames(mxd)[0]          #Check to see if ID field types and name match     try:         SectorShapeFields = arcpy.ListFields (SectorShape,ID)         SectorTableFields = arcpy.ListFields (SectorTable,ID)         arcpy.AddMessage("Finished Listing Fields")         nameCheck = SectorShapeFields[0].name == SectorTableFields[0].name         typeCheck = SectorShapeFields[0].type == SectorTableFields[0].type     except:         arcpy.AddMessage("Failed to List Fields")          #If they do not match, add new fields to correct.     if not nameCheck:         arcpy.AddMessage("Field Names do not match!  Adding new field to circumvent...")     if not typeCheck:         arcpy.AddMessage("Field Types do not match!  Adding new field to circumvent...")         if SectorShapeFields[0].type != SectorTableFields[0].type:             try:                 arcpy.AddField_management(SectorShape, newID, SectorTableFields[0].type,10)                 arcpy.CalculateField_management(SectorShape, newID,"!"+ID+"!", "PYTHON")                 arcpy.RefreshTOC()             except:                 arcpy.AddMessage("Error in Creating Field. Does it already exist?")               #Join the two together     arcpy.AddMessage("Performing Join")     arcpy.AddJoin_management( SectorShape, newID, SectorTable, ID)     arcpy.SelectLayerByAttribute_management (SectorShape,"NEW_SELECTION","Angle_Deviation>"+str(MaxDev))     df.zoomToSelectedFeatures()          #Field Mapping ...     myMap = arcpy.FieldMappings()     myMap.addTable(PopulationFC)     myMap.addTable(SectorShape)          #Verify the field merge rule for the pop10 field.     fIndex = myMap.findFieldMapIndex("POP10")     arcpy.AddMessage(str(fIndex))     popFieldMap = myMap.getFieldMap(fIndex)     arcpy.AddMessage(str(popFieldMap.mergeRule))     popFieldMap.mergeRule = 'Sum'     arcpy.AddMessage(str(popFieldMap.mergeRule))     popFieldMap2 = popFieldMap          ##Test     fIndex = myMap.findFieldMapIndex("POP10")     arcpy.AddMessage(str(fIndex))     popFieldMap = myMap.getFieldMap(fIndex)     arcpy.AddMessage(str(popFieldMap.mergeRule))          #Perform Spatial Join     arcpy.AddMessage("Performing Spatial Join")     arcpy.SpatialJoin_analysis(SectorShape, PopulationFC, OutputFC,"JOIN_ONE_TO_ONE","",myMap,"INTERSECT")          #Add Result and Symbolize     arcpy.mapping.AddLayer(df,arcpy.mapping.Layer(OutputFC))     translayer = arcpy.mapping.ListLayers(mxd,"",df)[0]     translayer.transparency = 50          arcpy.RefreshActiveView()      **EDIT - Below is the code with the solution implemented!**               import arcpy,sys,os          #Get the Files involved, set some variables.     SectorTable = sys.argv[1]     SectorShape = sys.argv[2]     MaxDev = sys.argv[3]     PopulationFC = sys.argv[4]     OutputFC = sys.argv[5]     DeviationField ="Angle_Deviation"     ID = "SectorID"     newID = "BP_ID"     mxd = arcpy.mapping.MapDocument('CURRENT')     df = arcpy.mapping.ListDataFrames(mxd)[0]          #Check to see if ID field types and name match     try:         SectorShapeFields = arcpy.ListFields (SectorShape,ID)         SectorTableFields = arcpy.ListFields (SectorTable,ID)         arcpy.AddMessage("Finished Listing Fields")         nameCheck = SectorShapeFields[0].name == SectorTableFields[0].name         typeCheck = SectorShapeFields[0].type == SectorTableFields[0].type     except:         arcpy.AddMessage("Failed to List Fields")          #If they do not match, add new fields to correct.     if not nameCheck:         arcpy.AddMessage("Field Names do not match!  Adding new field to circumvent...")     if not typeCheck:         arcpy.AddMessage("Field Types do not match!  Adding new field to circumvent...")         if SectorShapeFields[0].type != SectorTableFields[0].type:             try:                 arcpy.AddField_management(SectorShape, newID, SectorTableFields[0].type,10)                 arcpy.CalculateField_management(SectorShape, newID,"!"+ID+"!", "PYTHON")                 arcpy.RefreshTOC()             except:                 arcpy.AddMessage("Error in Creating Field. Does it already exist?")               #Join the two together     arcpy.AddMessage("Performing Join")     arcpy.AddJoin_management( SectorShape, newID, SectorTable, ID)     arcpy.SelectLayerByAttribute_management (SectorShape,"NEW_SELECTION","Angle_Deviation>"+str(MaxDev))     df.zoomToSelectedFeatures()          #Field Mapping ...     myMap = arcpy.FieldMappings()     myMap.addTable(PopulationFC)     myMap.addTable(SectorShape)          #Verify the field merge rule for the pop10 field.     fIndex = myMap.findFieldMapIndex("POP10")     arcpy.AddMessage(str(fIndex))     popFieldMap = myMap.getFieldMap(fIndex)     arcpy.AddMessage(str(popFieldMap.mergeRule))     popFieldMap.mergeRule = 'Sum'     arcpy.AddMessage(str(popFieldMap.mergeRule))          myMap.replaceFieldMap(fIndex,popFieldMap)          ##Test     fIndex = myMap.findFieldMapIndex("POP10")     arcpy.AddMessage(str(fIndex))     popFieldMap = myMap.getFieldMap(fIndex)     arcpy.AddMessage(str(popFieldMap.mergeRule))          #Perform Spatial Join     arcpy.AddMessage("Performing Spatial Join")     arcpy.SpatialJoin_analysis(SectorShape, PopulationFC, OutputFC,"JOIN_ONE_TO_ONE","",myMap,"INTERSECT")          #Add Result and Symbolize     arcpy.mapping.AddLayer(df,arcpy.mapping.Layer(OutputFC))     translayer = arcpy.mapping.ListLayers(mxd,"",df)[0]     translayer.transparency = 50          arcpy.RefreshActiveView()


title: Filter Points Layer by Polygon Layer

text: I come from the Java world and go straight to C# and ArcGis 10 (I have to make a simple task for the university). So this means I don't really know much of C# but I can live with it and don't know ArcGis API at all. All I have is an example from the classes. It is about how to filter the cities by coordinates X and Y and a `Buffer`. So we have layers `City`(Point) and `States`(Polygon). My task is to make other functionality, and by functionality I mean method, that takes the name of 1 state and get all the cities from state's area.   In the example code that we have some `ISpatialFilter` that we use with the buffer like this:               public List<City> SearchCities(...)     {         IFeatureClass citiesClass = MyHelper.GetLayer(Layers.CITIES);              ISpatialFilter spatialFilter = getSpatialFilter(centerX, centerY, distance);              IFeatureCursor citiesCursor = citiesClass.Search((QueryFilter)spatialFilter, true);         .....         .....     }          private static ISpatialFilter getSpatialFilter(double centerX, double centerY, double distance)     {         IPoint center = getCenter(centerX, centerY);         IGeometry buffer = getBuffer(distance, center);              ISpatialFilter spatialFilter = new SpatialFilterClass();         spatialFilter.Geometry = buffer;         spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;              return spatialFilter;     }               private static IGeometry getBuffer(double distance, IPoint center)     {         return (center as ITopologicalOperator).Buffer(distance);     }          private static IPoint getCenter(double centerX, double centerY)     {         IPoint center = new Point();         center.PutCoords(centerX, centerY);         return center;     }      So I decided to reuse some of the code and only change the way the data is filtered. So, I came up with something like this:               private IQueryFilter GetQueryFilter(string stateAbbr)     {         IGeometry stateGeometry = getGeometry(stateAbbr);              ISpatialFilter spatialFilter = new SpatialFilterClass();         spatialFilter.Geometry = stateGeometry;         spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;              return (QueryFilter)spatialFilter;     }          // TODO fix it     private IGeometry getGeometry(string stateAbbr)     {         IFeature state = GetStateByAbbr(stateAbbr);              return (state as ITopologicalOperator).Boundary;     }          private IFeature GetStateByAbbr(string stateAbbr)     {         IFeatureClass statesClass = MyHelper.GetLayer(Layers.STATES);         IQueryFilter queryFilter = GetStateQueryFilter(stateAbbr);         IFeatureCursor stateCursor = statesClass.Search(queryFilter, false);              return stateCursor.NextFeature();     }          private IQueryFilter GetStateQueryFilter(string stateAbbr)     {         ISpatialFilter spatialFilter = new SpatialFilterClass();         spatialFilter.WhereClause = StatesWhereBuilder(stateAbbr);              return (QueryFilter)spatialFilter;     }          private string StatesWhereBuilder(string stateAbbr)     {         return (MyHelper.isEmpty(stateAbbr)) ? "" : MyHelper.CreateEqualStatement(StateConst.STATE_ABBR, stateAbbr, true);     }      I know it is probably very amusing for some of you, but this is all I came up with :] The problem is that I get `NullReferenceException` after calling `(state as ITopologicalOperator).Boundary` but the debug shows that `state` is not null. So after the `as` something happens. I might be wrong with the `ITopologicalOperator`. If so, what should I use instead? Thanks in advance. EDIT: Please tell me what do you think about this question (how it is exposed, etc.) I changed the GetGeometry() method to:               private IGeometry getGeometry(string stateAbbr)     {         IFeature state = GetStateByAbbr(stateAbbr);         IPolygon statePolygon = state.Shape as IPolygon;              return (statePolygon as ITopologicalOperator).Boundary;     }      and now I don't have NullPointerException, but I don't get any result.


title: Calculating Topographic Exposure with Raster Calculator

text: I need you help trying to calculate **Topographic Exposure** (TOPEX) on a DEM in ArcGIS. I am trying to figure out **wind exposure** on a DEM with 25m cells. There is a great blog post where someone explains how to do it but he uses GRASS. Jamie Popkin's GIS Blog I thought I would be able to put the same equation into ArcGIS Raster Calculator but it says the equation is invalid. I've tried playing with it as much as I can but I can't seem to make it work. His first equation is: r.mapcalc incl4north = "atan(((dem - dem[4,0])) / ((25 * 4)))" And then he goes on to string a bunch of them together ![enter image description here](http://i.stack.imgur.com/N8ZvN.png) Does anyone know how the equation should appear in Raster Calculator in either ArcGIS or QGIS? I would be super grateful for any assistance. Thank you!


title: BuildCircleMbr with 1km of accuracy - spatialite

text: I want to create a Circle MBR in spatialite of 1km for example from lat/lon in WGS84 I use `BuildCircleMBR` but i am not able to generate this well. My query `SELECT BuildCircleMbr(40.0,-4.0,1, 1000)` I have read the third argument (1000) should be in degrees and not in meters but i don't know how to do this. Any idea? Thanks


title: Eye-friendly colours for long periods of work

text: Apologies if this question doesn't fit the rules here. I am digitising a large amount of remote footways and I will be panning, zooming a lot and I've done 2 days on this so far and always get quite a bad headache. Is it the colours I'm using? I have a pre-existing footway/road layer that I am adding features to. A polygon of footway which I am creating new features over, a buffer around this footway and a plain polygon layer that I use to cover areas I have completed. Can somebody recommend some alternative colours to use to minimise eye strain and headaches? ![enter image description here](http://i.stack.imgur.com/csoDV.png)


title: Smoothly interpolating Polygons from set of coordinates?

text: I'm still quite new to QGIS and I'm using it to display and analyze my soil data. I have point information on the distribution of the soils in a specific area and I need to create a soil map from that. I'm using QGIS 2.2.0 Valmiera and I used the voronoi polygons function to interpolate the areas between the points. Now these are very edgy but in soil science it is custom to have smooth bent lines in the soil maps because that rather represents the nature of soil distribution then straight lines. Is there another function or algorithm that interpolates smooth polygons, or is there a way to smooth the voronoi polygons? I have seen that there a some threats on smoothing polygons already but none refer to exactly this problem. Below is a picture of the polygons I've got so far. ![Below is a picture of the polygons I've got so far:](http://i.stack.imgur.com/MOPy8.jpg)


title: Convert from gps coordinates to EPSG:32640

text: I need to convert latitude and longitude witch i get from gps to EPSG:32640 (in C#). What is best and fastest method for that? Thanx


title: gdal_contour transparent contour as png

text: I started using the GDAL Library in a c++ project. It is a great library, but I am a bit stuck right now with different aspects of it. What is the best way to use gdal_contour in my c++ programm? Should i use the system() function or simply grab the source and use it? I tried to get the contours of a hightmap with gdal_contour as a PNG or SVG(or different vector format, like eps or so) but had no success so far. Is there a list wich output formats are supported. Until now I only got shp and dxf files. Do I have to burn the shp file to a png with gdal_translate? I tried this and had some troubles with the alpha channels. Maybe someone has the commands for me to get a transparent contour as png. Sometimes the contour lines look a bit to pixelated, is there a way to smoothen things up with more antialiasing or so? Is there a way to get a colored contour automatically, or do i have to use a color table txt file? I tried the following where output.tif was a transparent tif. but had no success.               gdal_contour -i 50 <DEM> contour     gdal_rasterize -burn 0 -l contour contour output.tif


title: LZW compression differences between Global Mapper and gdal_translate

text: I compressed an GEOTIFF image with both gdal_translate and Global Mapper with LZW compression. The one Global Mapper compressed is much smaller, is almost 30% of the one gdal_translate compressed. gdal_translate -co COMPRESS=LZW -co TILED=YES -co PREDICTOR=2 1.tif lzw1.tif How can I compress it further like Global Mapper? Thanks!


title: Arcpy how to conduct addition in a script, function needed?

text: I'm trying to work on my first arcpy script and am running into some troubles. The objective is to create a new field, check several existing fields (Cnt1,2,3 and Pop1,2,3) for missing data (-99 value) and if find missing data populate new field with a 0, if not then populate the field with the result from a simple calculation. From a separate related post I got some feedback that the problem was that I am concatenating the values, not adding them and this if-elif statement might need to be within a function and then called. Any advice greatly appreciated! ![enter image description here](http://i.stack.imgur.com/6Z4Pz.jpg)


title: Calculating sinuosity of a path with GME in ArcGIS 10?

text: Does anyone have any experience with pathmetrics in GME? I would ideally like to use the old Hawths's tool "line metrics" for calculating sinuosity of a daily path but I have ArcGIS 10. I downloaded GME & R and have run some things successfully (e.g. min convex polygon), but have had lots of problems trying to run "path metrics" in GME. It basically won't run and I don't know why. I get an output that says 0 files created. Also, the format of my 'time' column in the original xls file doesn't seem to get copied when I create a polyline - it ends up with no time codes (shows 0:00:00 in every row). this means I can't easily calculate the time intervals (and so the rate of path movement). Is there a specific time format that is needed? Lastly, imagining that I do get it to run, how do I calculate the path sinuosity from the segment lengths, turn angles etc..? Thanks!!


title: What could cause garbage hillshade output

text: I've used Hillshade successfully in the past, but for this DEM I keep getting a weird output. It's a 30-meter TIFF downloaded from the National Elevation Dataset. ![normal looking DEM](http://i.stack.imgur.com/zSwQl.jpg) And this is what the output looks like: ![crazy looking hillshade](http://i.stack.imgur.com/PHbsW.jpg) Kind of cool, but not at all correct. I've done other Surface analyses with the same DEM (e.g. Viewshed) and they ran without any obvious problem. Any ideas what may be causing the problem, and/or how to fix it? **UPDATE:** Data came in UTM, so I re-projected it to the appropriate South Carolina State Plane. (I also copied it from the original TIFF into a geodatabase, ESRI raster.) It took a lot longer to process the same area, which got my hopes up, but the **results are a similar "mesh"** of gray squares/lines. I'm going to try hillshading another tile from the same source and see if it has the same problem. **MORE UPDATE:** The tile to the south hillshades just fine (without reprojecting or anything), and processes very quickly, so this seems to indicate there's something specific to the tile I want that is corrupted (or otherwise borked). It was downloaded from the NRCS Geospatial Gateway -- South Carolina, Richland County, 30m DEM, and it's tile `ned30m33080.tif`. (This FTP link should work for the next few days.)


title: Trying to identify the program which created tab file in question

text: i am trying to identify format or the program which created following geo tiff and tab files. Basically I am trying to update map images in our custom program which uses following tab format and tiff files and no way to contact original developer as he is no where to be found. ## Following is content of one of the tab files:               !table     !version 300      !charset WindowsLatin1               Definition Table        File "1000001.TIF"        Type "RASTER"        (524448,182304) (0,6912) Label "Pt 1",      (524448,192672) (0,0) Label "Pt 2",         (534816,192672) (6912,0) Label "Pt 3",       (534816,182304) (6912,6912) Label "Pt 4"       CoordSys Earth Projection 8, 79, "m", -2, 49, 0.9996012717, 400000, -100000       Units "m"       * * * If anyone can shed any light I would be very grateful. Thanks in advance.


title: How is Precise Point Positioning (PPP) exactly done?

text: I am looking for a positioning accuracy that is higher than normal consumer- grade GNSS receiver (mostly GPS/GLONASS hybrids with +- 5m precision). I came across RTK with dual receiver setup (base station + rover) which is supposed to do cm-level precision although it has the "disadvantage" while it needs a base station or another source of correction data that one may not be able to supply. I am interested in getting the most of a single receiver setup. I read Precise Point Positioning (PPP) is able to achieve sub-meter (or maybe even 20cm- level) accuracy. I think PPP needs raw data from sensor which I can do. I will do readings from a moving vehicle in my application so I can not go with stationary "averaging" over time for getting the accuracy. However, I do not need a real-time solution as I will only collect mobile data and evaluate it later so I can do positioning in post-processing. U-blox also has a page on PPP where I get the impression their solution is just an ordinary raw data chip with a "better firmware" that can do some additional calculations (quote: "... by stabilizing measurements of the distance between GNSS satellites and the receiver (pseudo-ranges) using carrier phase tracking") so it has better precision than ordinary chip. I read the Wikipedia article but I do not really understand how PPP works. So how exactly is this PPP correction technically done? Where do the corrections come from? I read RTKLIB can do PPP too. My location is Czech Republic.


title: OpenLayers popup - random losing click event, have to refresh page

text:
I have a OpenLayers, Mapfish, GeoServer map on a webpage. I am using `OpenLayers.Popup` and `getFeaturesInfo` to get the attributes from GeoServer (2.1.3). I have several wms layers with markers and text labels. `The problem is that suddenly the the markers do not react to click.` There is no action in FireBug at all. If I refresh the page and start all over again I can get click events on the markers. This is very unstable and since I get no js error in Firebug it is not easy to tell what happens. It just seems to lose the click event. Is there anything here in the code that could be done better? I am using a ExtJS lightbox plugin inside the popup.                        info = new OpenLayers.Control.WMSGetFeatureInfo({                     url: 'http://kart.naturkart.no/geoserver/wms',                     title: 'Identify features by clicking',                     queryVisible: true,                     infoFormat:'application/vnd.ogc.gml',                     eventListeners: {                         getfeatureinfo: function(event) {                             if (popup) {                                 map.removePopup(popup);                             }                             var contentHtml = '';                             // Manage the features                             if (event.features.length > 0) {                                 for (var i = 0; i < event.features.length; i++) {                                     var feature = event.features[i];                                     // Identify the type                                     if (feature.gml.featureType == 'Architect_1' || feature.gml.featureType == 'Architect_2') {                                         // Create HTML content for this feature type                                         var gardsnavn = feature.attributes['gardsnavn'];                                         var navn = feature.attributes['navn'];                                         var epost = feature.attributes['epost'];                                         var lokalitetsid = feature.attributes['lokalitetsid'];                                         var x = feature.attributes['eu8933x'];                                         var y = feature.attributes['eu8933y'];                                          if (gardsnavn === null) {                                             gardsnavn = '';                                         }                                         if (navn === null) {                                             navn = '';                                         }                                         contentHtml = contentHtml + '
' + feature.attributes['art'] + '
';                                         contentHtml = contentHtml + gardsnavn + ' | ' + navn + '
';                                         if (feature.attributes['bilde']>0){                                         var imgcount = feature.attributes['bilde'];                                         var lokalitetsid = feature.attributes['lokalitetsid'];                                         var i=0;                                         for (i=1;i<=imgcount;i++){                                         contentHtml = contentHtml + ' ';                                         }                                         };                                         contentHtml = contentHtml + 'Les mer om dette minne.';                                         contentHtml = contentHtml + 'Direkte link til dette minne.';                                         contentHtml = contentHtml + 'Meld inn opplysninger om dette eller andre minner.
';}                                      else if (feature.gml.featureType == 'gamlenavn_archi') {                                             var navn = feature.attributes['navn'];                                             var kommentar = (feature.attributes['kommentar'] === null?'':feature.attributes['kommentar']);                                             var x = feature.attributes['eu8933x'];                                             var y = feature.attributes['eu8933y'];                                             contentHtml = '
'+ navn + '
';                                              contentHtml = contentHtml + kommentar + '
';                        contentHtml = contentHtml + 'Direkte link til dette stedsnavn.';                                                                  }                                 }                             } else {                                 // Doesn't show any popup if no features.                                 return;                             }                             var popup = new Popup(                                     "chicken",                                     map.getLonLatFromPixel(event.xy),                                     null,                                     contentHtml,                                     null,                                     true                                     );                             popup.autoSize = true;                             map.addPopup(popup, true);                             Ext.ux.Lightbox.register("a[rel^=lightbox]",true);                         }                     }                 });                 map.addControl(info);                 info.activate();


title: How to use Heatmap plugin and its output?

text: I have a layer of points which I would like to depicted as a heatmap concentration of points for visualization. In QGIS 1.8 I tried the built-in heat map plugin, leaving the default values as is, and got this grey box. What's the next step after getting this layer produced by the heatmap plugin? I'm hoping to produce something similar to this output from here. http://www.sethoscope.net/heatmap/ P.S. I actually tried this and the output was beautiful but couldn't get the image as a layer into qgis. ![enter image description here](http://i.stack.imgur.com/35gfe.jpg)


title: Accessing coordinates of polygon

text: I have designed a polygon in a web page and the coordinates are stored in the database of the type polygon. However when I check the database, I find GEOMETRY - 120B under polygon field. What does it mean and how can I know the coordinates? In the end, I want to know whether a point is inside or outside the polygon?


title: pgRouting and Openlayers error when insert controls variable

text: I've got another problem here :( I'm following the foss4g_w12 tutorial for building a pgRouting WebApp using openlayers, everything goes well until I added the following statement to my code:               controls = {               start: new OpenLayers.Control.DrawFeature(start, SinglePoint),               stop: new OpenLayers.Control.DrawFeature(stop, SinglePoint)             }      the map suddenly showed up like this: ![enter image description here](http://i.stack.imgur.com/paQuT.jpg) I already refresh the browser and it didn't change, after I delete the code it goes normal again. here is my entire javascript code:                   var SinglePoint = OpenLayers.Class.create();         SinglePoint.prototype = OpenLayers.Class.inherit(OpenLayers.Handler.Point, {     createFeature: function(evt) {         this.control.layer.removeFeatures(this.control.layer.features);         OpenLayers.Handler.Point.prototype.createFeature.apply(this, arguments);         }         });              var start_style = OpenLayers.Util.applyDefaults({             externalGraphic: "start.png",             graphicWidth: 18,             graphicHeight: 26,             graphicYOffset: -26,             graphicOpacity: 1         }, OpenLayers.Feature.Vector.style['default']);              var stop_style = OpenLayers.Util.applyDefaults({             externalGraphic: "stop.png",             graphicWidth: 18,             graphicHeight: 26,             graphicYOffset: -26,             graphicOpacity: 1         }, OpenLayers.Feature.Vector.style['default']);               var result_style = OpenLayers.Util.applyDefaults({              strokeWidth: 3,              strokeColor: "#ff0000",              fillOpacity: 0          }, OpenLayers.Feature.Vector.style['default']);               var map, parser, start, stop, downtown, result, controls;          function init() {     map = new OpenLayers.Map({         div: "map",         projection: "EPSG:900913",         displayProjection: "EPSG:4326",         numZoomLevels: 18,         // approximately match Google's zoom animation         zoomDuration: 10     });          // create Google Mercator layers     var gphy = new OpenLayers.Layer.Google(         "Google Physical",         {type: google.maps.MapTypeId.TERRAIN}     );     var gmap = new OpenLayers.Layer.Google(         "Google Streets", // the default         {numZoomLevels: 20}     );     var ghyb = new OpenLayers.Layer.Google(         "Google Hybrid",         {type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20}     );     var gsat = new OpenLayers.Layer.Google(         "Google Satellite",         {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22}     );          // create Bing layers          // API key for http://openlayers.org. Please get your own at     // http://bingmapsportal.com/ and use that instead.     var apiKey = "AqTGBsziZHIJYYxgivLBf0hVdrAk9mWO5cQcb8Yux8sW5M8c8opEC2lZqKR1ZZXf";          var veroad = new OpenLayers.Layer.Bing({         key: apiKey,         type: "Road",         wrapDateLine: true     });     var veaer = new OpenLayers.Layer.Bing({         key: apiKey,         type: "Aerial",         wrapDateLine: true     });     var vehyb = new OpenLayers.Layer.Bing({         key: apiKey,         type: "AerialWithLabels",         wrapDateLine: true     });          // create OSM layers     var mapnik = new OpenLayers.Layer.OSM();          // layer WMS Lokal     var admin1 = new OpenLayers.Layer.WMS("Administrasi Kota Tangerang Selatan",                 "http://localhost:9090/geoserver/wms",                 {layers: "rute:admin1", format: "image/png", transparent: true},                 {maxScale: 150}             );     var jalan1 = new OpenLayers.Layer.WMS("Jalan Kota Tangerang Selatan",                 "http://localhost:9090/geoserver/wms",                 {layers: "rute:jalan1", format: "image/png", transparent: true},                 {maxScale: 150}             );     start = new OpenLayers.Layer.Vector("Start point", {style: start_style});     stop = new OpenLayers.Layer.Vector("End point", {style: stop_style});     downtown = new OpenLayers.Layer.Vector("Downtown data area", {style: result_style});     result = new OpenLayers.Layer.Vector("Routing results", {style: result_style});          // create a vector layer for drawing     //var vector = new OpenLayers.Layer.Vector("Editable Vectors");          map.addLayers([         gphy, gmap, gsat, ghyb, veroad, veaer, vehyb, mapnik, admin1, jalan1, start, stop, downtown, result     ]);     parser = new OpenLayers.Format.WKT();     var feature = new OpenLayers.Feature.Vector();         downtown.addFeatures([feature]);         controls = {               start: new OpenLayers.Control.DrawFeature(start, SinglePoint),               stop: new OpenLayers.Control.DrawFeature(stop, SinglePoint)             }     map.setCenter(new OpenLayers.LonLat(11879804.55, -700076.0324),12);     map.addControl(new OpenLayers.Control.LayerSwitcher());     //map.addControl(new OpenLayers.Control.EditingToolbar(vector));     map.addControl(new OpenLayers.Control.Permalink());     map.addControl(new OpenLayers.Control.Scale());     map.addControl(new OpenLayers.Control.MousePosition());     //map.zoomToMaxExtent();          }      please help me, because my routing function can't started yet regards


title: Explode/Decompose Polyline Network with MapInfo

text: I have a network of polylines which possibly include multi polylines as well. I wish to decompose these polylines into their fundemental line segments so I can produce the following table:               Link_ID | BegNodX | BegNodeY | EndNodeX | EndNodeY | Distance     --------------------------------------------------------------     0001    | 554612  |  456375  | 557864   |  455487  |  0.348     0002    | 554632  |  456355  | 557841   |  456522  |  3.454     ....      I've noticed there is a Mapinfo tools such as Polybrk.mbx by Mr. Dressel but this tool only explodes one polyline selection. I need a whole network exploding and I do not wish to do each one manually.


title: Query Huge data passing in POST method

text: i have a feature layer in arcGIS, which contains huge number of features, from the Java desktop application i need to make a call to open a web browser and pass the id's of the features, which has to be shown on map. how can i pass the id's from java application to browser where number of id's are 20K. where we can not use GET method to call, and if i pass the data in POST, How can i get the POST data using dojo is any other way to do this?


title: Mapserver - label at beginning of line string (WMS)

text: I am creating a WMS service made up of linestrings pulled from a POSTGIS database I want to label the start of each line with a sequence number. At present I can label the centre of each with:-                PROCESSING 'LABEL_NO_CLIP=1'  (to prevent repetition across tiles)      ....      LABEL           COLOR 132 31 31           TYPE TRUETYPE           FONT "arialbd"           SIZE 12           POSITION AUTO      END      This places the label in the middle of the line string. Any idea how to place it at the start? Note that GEOMTRANSFORM "start" does not work with labels


title: Unselect all layers and clear the legend in Openlayers

text: Is the there a way of clearing the legend in Openlayers by unselecting all the layers? We have a map that starts with several layers. I was thinking that we might be able to add a button to the map that would do this.


title: Creating a radius of influence around points

text: I have a vector file of points. What I need to do is create a circle with a radius of 12km round each point. Would it be possible to then find out what percentage of my total map is covered with circles and what area is uncovered?


title: How to get LonLat Points between a LineString or 2 LonLat Points

text: i have a marker and i want to move him on a LineString. So i know start LonLat and end LonLat and also the time between. I tried to calculate it my self with "getLayerPxFromLonLat" and then count x+1 y+1 and after that i switch them back with "getLonLatFromPixel". So i can calculate it my self (e.g. 10 Seconds = 5 Points between start and end) but when i the map moves the function got wrong data. I hope u have a idea for that Problem or an Idea how to get Points on a LineString.


title: Large dataset processing via cloud computing

text: I need to process a large dataset. Aside from processing locally on my computer, can I install, for instance, arcgis10 on a cloud service like amazon. I don't have esri server, only arcgis10. I'm not interested in storage or distribution, just processing. I have ample bandwidth so upload/download isn't an issue. Might seem like a silly question, but thought it was worth asking.


title: How to export png to file system with python?

text: I use python to make a plugin that contains shape file,color_bar,time_bar 3 items above now I want to export my result to a png image After I find the method to do it I find a similar way to do it http://www.qgis.org/pyqgis-cookbook/composer.html?highlight=layout But according to the site ,it's can only print the shape file in qgis It's can't make additional items(in my situation is time bar and color bar)export as well Is there other way to do it with python?


title: Map project a raster having separate latitude and longitude raster bands

text: I have a satellite image where one of the channels is a raster having the longitude values, another channel a raster containing the latitude values. What would be the best way for map projecting such files? Is there some standard routines or tools (I use Python / gdal preferrably)? My only thought so far is to loop through the pixels, convert the lat/long value of a given pixel and drop it into the closest matching pixel of a raster already in the projection I want to have (where it seems not that straightforward to actually find the closest matching value) -- but somehow I wonder if there is a better and quicker way?


title: Which data formats are best for Open Data distribution?

text: What are the pros and cons of different data formats (performance, file size, etc.) when considering open data distribution? Our organisation wants to publish data as open data. However, there is no clear idea on which data formats to use. Ofcourse, the more 'open' a data format is, the easier it is to use. Which data formats are the most 'open' and therefore most usable for the distribution of Open Data when taking the following types in consideration?:   * raster data (I'm thinking: GeoTIFF, Erdas Imagine IMG?)   * vector data (I'm thinking: GML, CSV, ESRI Shapefile, DXF?)   * tabular data (I'm thinking: CSV?)   * 3D data (I'm thinking: CityGML?)   * 3D point coulds / LIDAR (I'm thinking: LAS?)   * am I forgetting something here? Also, if there is documentation about open data formats I'm very interested if you would like to share. Cheers! Mark


title: How to properly import QImage from PyQt4.QtGui?

text: I am starting with QGIS python console. Naturally, I came across the PyQGIS Developer Cookbook. While trying to cut my way through it, I encountered a lot of problems with "undefined attributes and names". Some of them are caused by not importing some library. I had the same problem as described here. > NameError: name 'QImage' is not defined Even though I imported the `PyQt4.QtGui` library. At least, I thought:               Python 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)] on N-PC     ## Type help(iface) for more info and list of methods.     help(QImage)     Traceback (most recent call last):       File "<input>", line 1, in <module>     NameError: name 'QImage' is not defined     import PyQt4.QtGui     help(QImage)     Traceback (most recent call last):       File "<input>", line 1, in <module>     NameError: name 'QImage' is not defined     from PyQt4.QtGui import *     help(QImage)     Help on class QImage in module PyQt4.QtGui:          class QImage(QPaintDevice)      |  QImage()      |  QImage(QSize, QImage.Format)      |  QImage(int, int, QImage.Format)      And the help goes on. Why did QGIS fail to load `QImage` by the 1st `import` command? Or what actually happened here? I am using QGIS 2.4.0, and I get the same issue in 2.0.1.


title: Selecting correct ArcGisFeatureLayer to query onSingleTap event?

text: My problem: I'm working on a Android (ArcGis) Attribute Editor following the sample code from their tutorials. The problem occurs when I have more than 1 ArcGisDynamicLayer and, more important, more than 1 ArcGisFeatureLayer. My question: What is the (correct) approach to IDENTIFY the right Feature Layer to query when I tap on a "field" displayed on the map? I read a similar question (Arcgis SDK for Android - Query multiple feature layers), but I didn't understand how to use IdentifyTask using only a single tap.


title: QGIS - Auto-Update Color-Ramp to Current Display Extent

text: When displaying raster files in Qgis with the color-ramp stretched to the min/max values of the current display extent, the min/max values don't update automatically when zooming or panning to a new area. The only way i know around it is to manually recalculate the values each time, which can be quite a pain. I know this functionality to be a basic feature in ArcMap, so maybe I'm just missing something in the display options. So if anyone knows how to get the min/max extent to auto-update, please let me know how to do it. Thanks! BTW: I'm currently using Qgis 2.4.0, but I encountered the same issue in all earlier versions I worked with.


title: Watershed Analysis

text: I am performing watershed analysis using Archydro tools in ArcGIS10. After calculating all flow direction, accumulation, drainage network when I perform watershed delineation, it gives me the same boundaries of watershed as of adjoint catchment. So my question is: can I use the adjoint catchment grid as basin boundaries? Plus, these boundaries don't cover the entire area to be mapped (shown in attached picture). Editing the polygons will change the basin area. What can I do to form a continuous map of river basins in my area of interest? Sample snapshot of my work is attached for further clarification. ![Sample](http://i.stack.imgur.com/H4BSX.jpg)


title: Get Features from OpenLayers Layer.Vector not working

text: polygon_layer = new OpenLayers.Layer.Vector("Polygone", {                 projection : projLonLat,                 strategies : [ new OpenLayers.Strategy.Fixed() ],                 protocol : new OpenLayers.Protocol.HTTP({                     url : 'polygons.geojson',                     format : new OpenLayers.Format.GeoJSON()                 })             });                       getFeatures(polygon_layer);              function getFeatures(maps_layer) {              alert(maps_layer.features.length);              var features = new Array();         for (var i = 0; i < maps_layer.features.length; i++) {                  if (maps_layer.features[i].cluster) {                 for (var j = 0; j < maps_layer.features[i].cluster.length; j++)                 features.push(maps_layer.features[i].cluster[j]);             } else {                 features.push(maps_layer.features[i]);             }              }     }      Output for alert = "0"; What am i doing wrong getting the features? There should be 6 Features, which are also correctly appearing in the application. The (valid) *.geoJSON looks like this:                   {       "type": "FeatureCollection",       "features": [         {           "type": "Feature",           "properties": {},           "geometry": {             "type": "Polygon",             "coordinates": [               [                 [...]]]}....


title: How to use the Maplex abbreviation dictionary only if space is not available for the full label?

text: Maplex has a nice feature which can abbreviate labels based on a dictionary that you configure. eg, replace STREET with ST. However, it seems to be an all-or-nothing setting, with Maplex always replacing the words in the dictionary, even when the abbreviation is not necessary. In the case below, Maplex has abbreviated STREET even though the word would easily fit in all cases: ![with abbreviation dictionary on](http://i.stack.imgur.com/3iBDg.png) ![with abbreviation dictionary off](http://i.stack.imgur.com/RCdrZ.png) In these examples, I want SKITTLE LN and TEMPERANCE LN to be abbreviated, but not the STREET names. Edit: According to the Help File entry mentioned by ChrisW: > the Maplex Label Engine first attempts to place the full text of the label, > then abbreviates any words found in the dictionary for labels that could not > be placed. Why then is Maplex abbreviating the word STREET in these examples below? ![enter image description here](http://i.stack.imgur.com/AcvzC.png) I've simplified the options so that the only Fitting Strategy used is Abbreviation Dictionary. The style is Centered Straight with Regular placement. Using a line featureclass in an fGDB. The selected line is a single feature, so it's clear there is ample room for the full label. There are no other layers drawn. It sure looks like a bug to me....


title: Create a fishnet in a 'for' loop with Arcpy

text: When running the following script, I get the error saying that I need a value for “origin_coord” and for “y_axis_coord”. So is it not enough to define a template dataset? But how can I fill the “origin_coord” and “y_axis_coord” parameters with Arcpy working with ArcMap 10.0? Because in the loop are a lot of shapefiles from different cities, so I need a loop.               import arcpy     import os     from arcpy import env          env.overwriteOutput = True     env.workspace = r"D:\Users\julia\out_09_09\urbanA"     env.qualifiedFieldNames = False     #make a list with input cities as shapefiles.          fcList = arcpy.ListFeatureClasses()          for shpFile in fcList:         shpFileName= os.path.splitext (shpFile) [0]         print shpFileName # works         # Process: Make Feature Layer         arcpy.MakeFeatureLayer_management(env.workspace + "\\" + shpFile, "shpFile_ly")         # join the city with the reclasstable         jointable = r"D:\Users\julia\out_09_09\reclass.dbf"         arcpy.AddJoin_management("shpFile_ly", "CODE", jointable, "joinCODE", "KEEP_ALL")         # Process: Copy Features         outFeatureClass = shpFileName + "_join.shp"         arcpy.CopyFeatures_management("shpFile_ly", outFeatureClass)         # Process: Polygon to Raster         outraster = shpFileName + "_join.img"         cellSize = 200         arcpy.PolygonToRaster_conversion(outFeatureClass, "CODE", outraster, "CELL_CENTER", "NONE", cellSize)         # Process: Make Raster Layer         arcpy.MakeRasterLayer_management(outraster, "outraster_ly")#works         # Process: Add Join          arcpy.AddJoin_management("outraster_ly", "Value", jointable, "code", "KEEP_ALL")         # Process: Copy Raster         outraster2 = shpFileName + "_join2.img"         arcpy.CopyRaster_management("outraster_ly",outraster2) #works         # Process: Raster to Polygon         outPolygons = shpFileName + "reclass.shp"         arcpy.RasterToPolygon_conversion(outraster2, outPolygons, "NO_SIMPLIFY", "Value")#works              # Process: Create Fishnet         cellSizeWidth = "200"         cellSizeHeight = "200"         outFishnet = shpFileName + "_net.shp"         arcpy.CreateFishnet_management(outFishnet, "", "", cellSizeWidth, cellSizeHeight, "0", "0", "", "NO_LABELS", outPolygons, "POLYGON")      ![enter image description here](http://i.stack.imgur.com/rhnqx.png) Is there someone who can help me? I’m trying it without the Iteration but it is not working even if I add a layer before? Why?               import arcpy     from arcpy import env     env.overwriteOutput = True     env.workspace = r"D:\Users\julia\erste_aufg"          # Process: Make Feature Layer     #arcpy.MakeFeatureLayer_management(r"D:\Users\julia\erste_aufg\de013l_hannover    \de013l_hannover\de013l_hannover.shp", "hannover_ly")          #Process: Create Fishnet     outFeatureClass = r"D:\Users\ju\ers\result\hannover_fisch.shp"     cellSizeWidth = '200'     cellSizeHeight = '200'     templateExtent = r"D:\Users\ju\ers\de013l_hannover\de013l_hannover\de013l_hannover.shp"          arcpy.CreateFishnet_management(outFeatureClass, "", "", cellSizeWidth, cellSizeHeight, '0', '0', "", "NO_LABELS", templateExtent, "POLYGON")


title: psql schema error while importing a table

text: I am trying to create a table from a shape file and import it into my postgis database. I am using pgadmin3 and postgresql 9.3 on Ubuntu 14.04. I have a database `postgis_test` and enabled postgis with `CREATE EXTENSION postgis;`. `SELECT postgis_full_version()` has the result "POSTGIS="2.1.3 r12547" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="UNKNOWN" RASTER" To create and import a table I use `shp2pgsql -s -1 -W "latin1" test_geom.shp test_geom.dbf | psql -h localhost -d postgis_test -U postgre` At the execution following errors are displayed               Shapefile type: PolygonM     Postgis type: MULTIPOLYGONM[3]     SET     SET     BEGIN     ERROR:  schema "test_geom" does not exist     ERROR:  current transaction is aborted, commands ignored until end of transaction block     ERROR:  current transaction is aborted, commands ignored until end of transaction block     ERROR:  current transaction is aborted, commands ignored until end of transaction block     ERROR:  current transaction is aborted, commands ignored until end of transaction block     ROLLBACK      I don't understand why it is creating the schema test_geom. Shouldn't it use the username or public or something? Thanks


title: How to identify "true arcs" in Arcmap?

text: This is really annoying: the only good way of creating circles in ArcMap is using true arc, (compound curves) but I need to convert those to vertex (densifying) because of a compatibility problem. I need info for geodatabase and SDO_geometry (arcsde) and yes there is a way. I use: `UPDATE layer1 a SET arctype = 'compound' WHERE has_compound_curves(a.shape) <> 0;` and then using EG Geotools to densify layer by layer.... But there should be easier way of doing this straightforward in Oracle or ArcMap, and I just want to densify the features that have true arc. I don't want any unnecessary vertices on straight lines and so on. So how can I identify true arcs in ArcMap? And how can I densify them as fast and easy as possible? Same for Oracle SQL. anyone? :)


title: book/tutorial on using rpy2 in python for kriging?

text: I am proficient in Python, but have no knowledge of R. Can you suggest a good tutorial in using **rpy2** or any other R library to perform **Block/Cokriging "Kriging"**


title: how to aggregate raster values for each grid in a fishnet-like map

text: here is my problem: I have a raster layer containing the elevation info. and another gridded layer (shapefile type). What I want to do is to calculate the mean elevation within each grid to have a single representative value. Does anyone know how to handle this? Thx in advance for any suggestion .