title: Batch conversion of netCDF to Raster in ArcGIS using Python

text: I'm having trouble using the Dimension Values (Optional) argument in arcpy.MakeNetCDFRasterLayer_md. Whatever I seem to do the output raster is always the first layer which is the default value. Below is my attempt to export all layers to raster in the time dimension. Code:               def extractAllNetCDF():              variable = "RRt_10m"         x_dimension = "lon"         y_dimension = "lat"         band_dimension = ""         dimension = "time"         valueSelectionMethod = "BY_VALUE"              outLoc = "E:/New Folder/"         inNetCDF = "E:/netCDFFiles/RRt.nc"                  nc_FP = arcpy.NetCDFFileProperties(inNetCDF)            nc_Dim = nc_FP.getDimensions()              for dimension in nc_Dim:                  top = nc_FP.getDimensionSize(dimension)                  for i in range(0, top):                      if dimension == "time":                                     dimension_values = nc_FP.getDimensionValue(dimension, i)                     nowFile = str(dimension_values)                          arcpy.MakeNetCDFRasterLayer_md(inNetCDF, variable, x_dimension, y_dimension, nowFile, band_dimension, dimension_values, valueSelectionMethod)                     arcpy.CopyRaster_management(nowFile, outLoc + nowFile + ".img", "", "", "", "NONE", "NONE", "")                      print dimension_values, i      The print method at the end will show the dates as they should be and the index i is also moving alon so there is no reason to think that there are other problems with the code other than the Dimension Value being incorrect and reverting to the default. Does anyone have any idea how to get the subsequent layers to export? Is there any code online that has a specific example of this using the Dimension Value argument other than empty quotes?


title: creating unique ID-s for my file geodatabase

text: I use ArcEditor licence level, ArcGIS 10.0 version, and we're talking about a file geodatabase. I have several feature classes (electric energy domain) and I want to connect them using relationship classes, so I want to use their ID's as foreign keys. I want to automatically create a significant unique ID for each feature. It'a a text type ID, which would be composed of two letters and several numbers (ex. CH00123, CH00124). The ID would give me three info : what type of feature we're talking about (channel- CH), on which street it is (00), and it's unique code (123). I don't want to type this code for every feature I create, I need something to automatically assign an ID to each feature. I am not sure how to handle this situation. Do I need to use a data integration software? an engine to generate me special ID's? Thanks for suggestions.


title: ST_MakeLine equivalent for CIRCULARSTRING in PostGIS

text: I'd like to make a CIRCULARSTRING in PostGIS 2.0.0 out of points that are calculated as part of the query (rather than constants). For a LINESTRING I can do this using ST_MakeLine(). Is there an equivalent for CIRCULARSTRING or some other way to do this? I have 3 points that are calculated in some way (in the example let's just use ST_GeomFromText):               SELECT ST_GeomFromText('POINT(10 10)'), ST_GeomFromText('POINT(15 10)'), ST_GeomFromText('POINT(10 20)')      What I want to get is `CIRCULARSTRING(10 10, 15 10, 10 20)`. If I run:               ST_LineToCurve(ST_MakeLine(ARRAY[ST_GeomFromText('POINT(10 10)'), ST_GeomFromText('POINT(15 10)'), ST_GeomFromText('POINT(10 20)')]))      That gives me `LINESTRING(10 10,15 10,10 20)`, not a `CIRCULARSTRING`. Same question for CURVEPOLYGON - what can I use instead of ST_MakePolygon (which returns "Shell is not a line" if I pass in a CIRCULARSTRING)?


title: How to calculate the sum of all cells in a raster layer upstream of a specific point?

text: I have a rainfall map which covers large area. However, I'm interested in the total value of rainfall over a specific subcatchment. If I don't have the shapefile which defines the boundary of this subcatchment, is it possible to calculate the sum of the rainfall values upstream of a specific point in ArcGIS? or Do I have to delineate the catchment upstream of this point first?


title: How to show an array of 2d points/polygons in Polymaps (Cartesian coordinate system)

text: I need to display these 2d points/polygons (Cartesian coordinate system) on Polymaps with GeoJson layer. I do not use the other layers. Do I need to convert xy coordinates to the lat lon? If so, how? How to use Cartesian coordinate system in Polymaps?


title: How to convert Singapore Onemap Co-ordinate to Google maps co-ordinate system

text: I have to fetch few 3000 co-ordinates (I am Geo-coding ) and convert it into Google map co-ordinate system. Following is the process I am doing I am calling Mapone.sg API for geo-coding in return it gives me x.y values (I am unable to understand the co-ordinate system ) and which needs to be converted to Google map Co-ordinate system. If you can give me a right direction to work on it will be nice. Mehul


title: Oracle MapViewer / Map Builder - Tile algorithm

text: Oracle MapBuilder generates three types of layers - theme, base and tiles. I try to create algorithm for tiled layers. Have you any idea how to decode this? I mean algorithm for width and height of tiles. For a tiles is generated XML like this:               <map_tile_layer name="TL_TEST" image_format="PNG" http_header_expires="168.0" concurrent_fetching_threads="3">        <internal_map_source data_source="10.0.0.57_lrs" base_map="BM_AGC" bgcolor="none"/>        <tile_storage root_path="/temp"/>        <coordinate_system srid="2180" minX="-413374.575" minY="193748.248" maxX="1418971.699" maxY="766162.025"/>        <tile_image width="256" height="256"/>        <zoom_levels levels="11" min_scale="1.0" max_scale="1.0E7" min_tile_width="0.06773333184979038" min_tile_height="610782.0913333333">          <zoom_level level="0" name="level0" description="" scale="1.0E7" tile_width="610782.0913333333" tile_height="610782.0913333333">          </zoom_level>          <zoom_level level="1" name="level1" description="" scale="1000000.0" tile_width="65440.938357142855" tile_height="65440.938357142855">          </zoom_level>          <zoom_level level="2" name="level2" description="" scale="100000.0" tile_width="6761.425365313653" tile_height="6761.425365313653">          </zoom_level>          <zoom_level level="3" name="level3" description="" scale="10000.0" tile_width="677.1420081300813" tile_height="677.1420081300813">          </zoom_level>          <zoom_level level="4" name="level4" description="" scale="1000.0" tile_width="67.73172195320296" tile_height="67.73172195320296">          </zoom_level>          <zoom_level level="5" name="level5" description="" scale="100.0" tile_width="6.773322418713311" tile_height="6.773322418713311">          </zoom_level>          <zoom_level level="6" name="level6" description="" scale="10.0" tile_width="0.677333243384311" tile_height="0.677333243384311">          </zoom_level>          <zoom_level level="7" name="level7" description="" scale="5.0" tile_width="0.3386666216921555" tile_height="0.3386666216921555">          </zoom_level>          <zoom_level level="8" name="level8" description="" scale="4.0" tile_width="0.2709332973537244" tile_height="0.2709332973537244">          </zoom_level>          <zoom_level level="9" name="level9" description="" scale="2.0" tile_width="0.13546665869200755" tile_height="0.13546665869200755">          </zoom_level>          <zoom_level level="10" name="level10" description="" scale="1.0" tile_width="0.06773333184979038" tile_height="0.06773333184979038">          </zoom_level>        </zoom_levels>     </map_tile_layer>


title: Change projection system .tif file

text: I have a list of .tif files which have the wrong projection system (NAD_1983_UTM_Zone_13N). However, I want to have those files in (WGS_1984). I want to do this by a script that can transform tif files from one projection system to another. I have found many scripts online, but none seem to work for me. anyone a script or idea?


title: Arcgis Javascript latitude and longitude convert to WGS_1984_UTM_Zone_36N

text: I am using Arcgis Javascript api. I have two points comes from Google maps like this **x:** -111.67095840000002 **y:** 32.7558935 I wanna use this points on my Arcgis map. My map spatial reference is using wkt WGS_1984_UTM_Zone.... I mean it is not using wkid.               theX = results[0].geometry.location.lng();     theY = results[0].geometry.location.lat();          console.log("Longitude: " + theX + " Latitude: " + theY);          var pointESRI = new Point(Number(theX), Number(theY), map.spatialReference);     console.log(pointESRI);          var mp = WebMercatorUtils.geographicToWebMercator(pointESRI);     console.log(mp);          pointESRI = new Point({ "x": mp.x, "y": mp.y, "spatialReference": map.spatialReference });      ![enter image description here](http://i.stack.imgur.com/eKF40.png) WebMercatorUtils.geographicToWebMercator() methot not converting point as appropriate to WGS84 reference.


title: Looking for WPS Client

text: I would want to have a simple WPS client, that would work with GeoServer. QGis (1.8.0, Lisboa) WPS plugin (v1.1.0) doesn't work*. Any client that is SIMPLE to use, user friendly, will do (the one built in GeoServer as demo is perfect, but I'd like a third party client). It could either be a desktop client, or a web client. I tried using OpenLayers but I got a .js file that had to be incorporated to a web page and I am no programmer to know how to do it. A simple webpage with cute little buttons is what I want. If anybody has knowledge of an existing open-source program that does support OGC standards (WPS mainly), or simple tutorials on how to set up a web client, I would be grateful to read your answers. Thanks * * * Running on Win XP. * * * *I posted another question regarding this error.


title: How do I calculate the coordinates of the median point of a polyline?

text: I have a shapefile with some polylines. For each one of these polylines I'd like to calculate the median of the polyline (i.e. the point that's equidistant from each end) and then, calculate its coordinates. I don't know if this is feasible.   The only tools I've found regarding this matter are Mean Center and Median Center, however this is not what I'm looking for. Any help would be much appreciated.


title: Use bbox in geoserver sql view

text: I'm doing this sql query using geoservers sql views:               SELECT kmeans, count(*), ST_BUFFER(ST_CENTROID(ST_Collect(geom)), qrt(count(*))*1000  ) AS geom     FROM (        SELECT kmeans(ARRAY[ST_X(ST_CENTROID(geom)), ST_Y(ST_CENTROID(geom))], 10) OVER (), geom           FROM data WHERE foo = %FOO%     )     AS ksub     GROUP BY kmeans     ORDER BY kmeans      This always returns the same regardless which bounding box is requestet via the wms. But I want to use the bounding box in my WHERE clause. I know that the user gets a different result then if there are more tiles used to dislpay a map but that will be ok for me. Regards Daniel


title: QGIS distance from points to nearest line

text: How can I get the distances from each point at a shapefile to the nearest line from another? I've tried to follow the directions from Measuring the distance between lines and points in QGIS and Shortest distance from point to line (big dataset), but both methods didn't worked. I'm using QGis 2.4 (Chugiak) Thanks


title: Arcgis Server - Highlight Identified Features

text: I am identifying features from an ArcGISDynamicMapServiceLayer All the sublayers can be toggle on/off and when identifying, only the visible layers are queried/returned. Below are the initIdentify and executeIdentifyTask functions.                  dojo.connect(map, "onLoad", function(map) {initIdentify(map);});                  function initIdentify(map) {                              //create identify tasks and setup parameters                 identifyTask = new esri.tasks.IdentifyTask("http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer");                                 identifyParams = new esri.tasks.IdentifyParameters();                 identifyParams.tolerance = 3;                 identifyParams.returnGeometry = true;                 identifyParams.layerIds = wmsLayer.visibleLayers;                                identifyParams.layerOption = esri.tasks.IdentifyParameters.LAYER_OPTION_ALL;                 identifyParams.width = map.width;                 identifyParams.height = map.height;                              dojo.connect(map, "onClick", executeIdentifyTask);             }             function executeIdentifyTask(evt) {                 if (visibleLayers.length === 0){                                         return;                 } else {                     identifyParams.geometry = evt.mapPoint;                     identifyParams.mapExtent = map.extent;                     identifyParams.layerIds = wmsLayer.visibleLayers;                                    identifyTask.execute(identifyParams, function(idResults) { fillInfoWindow(idResults, evt); });                                    }       }      When the identify function is run, is there a way to highlight/select those features that are being identified?


title: Error 000732 - input raster does not exist when implementing script tool arcpy

text: I am trying to implement a script tool in ArcGIS 10 that completes a few equations with raster datasets. I seem to be getting the same error no matter what I try:               <type 'exceptions.RuntimeError'>: ERROR 000732: Input Raster: Dataset 5 does not exist or is not supported     Failed to execute (TWI).      TWI being the name of the script tool. My code is below:               import arcpy     from arcpy import env     import os     from arcpy.sa import *     arcpy.CheckOutExtension("spatial")          DEM = arcpy.GetParameterAsText(0)     cellsize = arcpy.GetParameterAsText(1)     FlowAcc = arcpy.GetParameterAsText(2)     K_Sat = arcpy.GetParameterAsText(3)     OutRaster = arcpy.GetParameterAsText(4)     OutSpace = arcpy.GetParameterAsText(5)          slope = (Slope(Raster(DEM)) * 1.570796 ) / 90     tan_slp = Con(slope > 0, Tan(slope), 0.001)     FlowAcc_scaled = (Raster(FlowAcc) + 1 ) * cellsize     TWI = Ln(FlowAcc_scaled / (Raster(K_Sat) * tan_slp))          arcpy.env.workspace = OutSpace          TWI.save(OutRaster)      I've got a feeling that the error is something to do with setting up the parameters correctly. It's something that I seem to be struggling to get my head around but as I understand it "DEM", "FlowAcc" and "K_Sat" are raster layers, "cellsize" is a long integer (also tried double), "OutSpace" is a workspace and "OutRaster" is a string. The error seems to be with the way is is saving/naming the raster. EDIT: Please see a screenshot of the parameter setup below (note I have tried the output workspace and output raster as both inputs and outputs and I get the same error) ![enter image description here](http://i.stack.imgur.com/MDabZ.png) **EDIT 2 - Answer** Ok the issue was that the second parameter 'cellsize' needs to be called with the 'arcpy.GetParameter' function rather than the 'arcpy.GetParameterAsText' function as it is an input value rather than a string.


title: ESRI Silverlight API 3.1 query, edit, update issues in the code behind

text: Using ESRI Silverlight API 3.1, .NET Framework 4.0, ArcGIS Server 10.0 Is there a way to update an attribute in the code behind without looping through the featureLayer.Graphics? Or a way to call FeatureLayer.Update without losing the edits made? We cannot hard code our layers into the xAML because we need them to be dynamic. Ideally we would like to run something similar to a QueryTask, returning a set of features that we could update the attributes of. To work around this we have attempted to use the FeatureLayer where and update functions. For this test there is only a single FeatureLayer in my map. The where and update successfully revises the layer to have a single entry -               FeatureLayer testLayer = (ESRI.ArcGIS.Client.FeatureLayer)featureLayer;     testLayer.UpdateCompleted += UpdateTask_ExecuteCompleted;     testLayer.Where = "TEST_ID=12345";     testLayer.Update();      // Update Event               private void UpdateTask_ExecuteCompleted(object sender, System.EventArgs args)     {       foreach (GraphicsLayer graphicsLayer in testEditor.GraphicsLayers)       {         if (graphicsLayer is FeatureLayer)         {           FeatureLayer featureLayer = graphicsLayer as FeatureLayer;                foreach (Graphic feature in featureLayer.Graphics)           {             if (Convert.ToInt32(feature.Attributes["TEST_ID"]) == 12345)              {                feature.Attributes["NAME"] = "Renamed";              }           }             }     }      This event handler will loop through my feature layer – now queried to 1 record, and update the name attribute. The issue at this point is that on my display I now only have the single record showing. If I reset the where clause to “1=1” and run an update I lose my edit. I can SaveEdits before I reset the layer which will save my edit. But I am hoping for a solution that will still allow the user to save all the edits at a single time when they are finished. Is there a way to update and maintain my edit? This all is a workaround due to not being able to query to get the feature to edit. If there is a way of doing that without iterating through every Graphic that would be the ideal solution to this. Are we missing the big picture here?


title: raster to polygon conversion

text: I 'm trying to automate the conversion of raster to polygon and the rasters files are in "TIF" format which is a "unsigned integer". I have the code below. However, I 'm getting an error while executing it. Could you have some suggestion please? Code:               import arcpy, os     from arcpy import env     from arcpy.sa import *          #To overwrite output     arcpy.env.overwriteOutput = True          #Set environment settings     env.workspace = "C:/Subhasis/Test/Neshanic_Python/extract"          outws="C:/Subhasis/Test/Neshanic_Python/extract/clip/"          #checkout ArcGIS spatial analyst extension license     arcpy.CheckOutExtension("Spatial")          inraster = arcpy.ListRasters("*", "TIF")          for i in inraster:         filename=os.path.splitext(i)[0]         outPolygons= outws + str(filename) + '.shape'         field = "VALUE"         arcpy.RasterToPolygon_conversion(inraster, outPolygons, "NO_SIMPLIFY", field)      > Error: Traceback (most recent call last): File > "C:\Subhasis\Test\Neshanic_Python\extract\Clip.py", line 23, in > arcpy.RasterToPolygon_conversion(inraster, outPolygons, "NO_SIMPLIFY", > field) File "C:\Program Files > (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\conversion.py", line 241, in > RasterToPolygon raise e


title: How to add multiple cartodb layers in viz.json to a leaflet layer control

text: I have create a visualisation in cartodb with 4 layers. The resulting viz.json is added to a leaflet map as follows:               cartodb.createLayer(cartomap, [url], {https:true})                     .addTo(cartomap)                     .on('done', function(layer) {                     // add individual layers     }      I want to accomplish that all 4 layers in the viz.json can be toggled separatly in the layer control (switching both the markers/features and the legend) The layer in on('done') is a valid leaflet layer and can be added to the Layer Control. However, it contains all sublayers which can then not be toggled individually. The layer also contains an array with all the "sub" layers in there, but these are not valid Leaflet layers (onAdd function is missing), so these can not be added to the layerControl individually. The layer.getSubLayer(i) also returns a not valid Leaflet layer. I have been trying for the past couple of days to get this working. Any help would be greatly appreciated.


title: On which TextBox event should I call ICompletionNotify.SetComplete?

text: I have a IToolControl derived control which contains a c# TextBox. In the IToolControl::OnSetFocus function documentation, it states that you need to call ICompletionNotify::SetComplete() when it's finished to let the application know that the control should lose focus. On another IToolControl ComboBox control, I called the SetComplete() function on the DropDownClosed event. On which TextBox event should I call ICompletionNotify.SetComplete()?


title: How do you hide a feature within a point layer?

text: I know this is basic, but for the life of me I can not remember how to hide/disable a specific feature point within a point layer. This seems like it should be a 2-3 click solution, but working with Arcmap I should always know better. One would think ESRI could simply make a right-click/hide point (polygon, or line), then have a "show all points in layer option" to bring it back. Thanks for the help!


title: Tilemill reversing inner and outer in multipolygon relations

text: Tilemill renders multipolygon relation (a middle part of a big lake (ID 2101875)) incorrectly: ![enter image description here](http://i.stack.imgur.com/Ww6fv.png) Multipolygon relations (inner/outer) seems to render mirror-like, compare: an island in OpenStreetMap but in Tilemill it looks like a lake: ![enter image description here](http://i.stack.imgur.com/tG3E7.png) Tilemill layer is a PostGIS source created by osm2pgsql. I read a similar post, downloaded the relation in JOSM and tried to find a gap but outer border of multipolygon consists of some line pieces connected to each other. If data would be wrong why does mapnik render this lake correctly in OpenStreetMap? Could someone help to find a mistake? Thank you an advance.


title: XY line tool will not execute

text: I wish to generate a road with approximately 1090 XY coordinates, my only problem is that my **XY To Line** tool will not execute, i get this message > **Error 999999** : error executing function. > > Failure to execute (XY To Line). what could the problem be? The model ran fine with small amounts of data but as soon as i tried the data I have for the road the error message occurred and now the tool won't even run small sums of data. Here is an image of my model and an example of how my input data is set up. .![enter image description here](http://i.stack.imgur.com/9yN6E.jpg)


title: Why do base and overlay only match at scale 1:1000 and 1:2000?

text: We are using GEOSERVER 2.1.4 to overlay WMS on GoogleMaps V3 with Openlayers . The problem we have is at different zoom levels the scale of the WMS overlay is different from the scale of the base layer. It is only at scale 1:1000 and 1:2000 that overlay and base layer coincide. Any suggestions? The part of the code:                       var map;             var untiled;             var tiled;             var pureCoverage = false;         OpenLayers.Lang.setCode('el');               // pink tile avoidance             OpenLayers.IMAGE_RELOAD_ATTEMPTS = 5;             // make OL compute scale according to WMS spec             OpenLayers.DOTS_PER_INCH = 25.4 / 0.28;                  function map_reg_init(){                 // if this is just a coverage or a group of them, disable a few items,                 // and default to jpeg format                 format = 'image/png';                 if(pureCoverage) {                     document.getElementById('filterType').disabled = true;                     document.getElementById('filter').disabled = true;                     document.getElementById('antialiasSelector').disabled = true;                     document.getElementById('updateFilterButton').disabled = true;                     document.getElementById('resetFilterButton').disabled = true;                     document.getElementById('jpeg').selected = true;                     format = "image/jpeg";                 }         var projEGSA=new OpenLayers.Projection("EPSG:2100");         var projLONLAT = new OpenLayers.Projection("EPSG:4326");             var projGOOGLE=new OpenLayers.Projection("EPSG:900913");                 var bounds = new OpenLayers.Bounds(                     2327363.1529655424, 4805801.76429202,                     2328463.0753220413, 4806793.842097636                 );         var scales = [500, 1000, 2000, 4000, 5000, 10000, 15000, 25000, 50000];         var resolutions = [];         for(var i = 0; i < scales.length; i++) {             resolutions.push(OpenLayers.Util.getResolutionFromScale(scales[i], 'm'));         }                      var options = {                     controls: [],                     maxExtent: bounds,             scales: scales,                          maxResolution: 156543.0339,                 numZoomLevels    : 22,                     allOverlays      : true,                         projection: "EPSG:900913",                     units: 'm'                 };         var mousePointerStyle = 'default';


title: How to obtain a DTM/DEM from DSM

text: I need to obtain a grid of bare ground elevation (removing all edifices, trees, roads and man-made elements) for quantitative geomorphological analysis from a high resolution (1 m) DSM obtained by orthoimages (not from LIDAR). Is there any suggestion you can provide me to filter my original grid?


title: OpenLayers.Layer.WMS and mapserver

text: I have a problem when I add a OpenLayers.Layer.WMS, I get a picture of my map but I can't do zoom on it. I read a lot of questions here but nothing help me. My mapfile is: **MAP ## CONFIG "MS_ERRORFILE" "/ms4w/tmp/ms_error.txt" ## DEBUG 2 ## CONFIG "CPL_DEBUG" "ON" NAME mapa2 SIZE 1400 900 IMAGETYPE png IMAGECOLOR 196 240 255 EXTENT -16.1737 27.6847 -13.1254 29.4506 SHAPEPATH "datos" FONTSET "C:/../fonts.txt" WEB TEMPLATE "/../index.html" IMAGEPATH "/../ms_tmp/" IMAGEURL "/ms_tmp/" METADATA "WMS_TITLE" mapa2 "WMS_SRS" "EPSG:4326 EPSG:4258" "WMS_ONLINERESOURCE" "../mapa2.map" "WMS_ENABLE_REQUEST" "*" END END PROJECTION "init=epsg:4258" END UNITS DD # FOTOS AEREAS LAYER NAME Foto1 TILEINDEX indice.shp TILEITEM location DATA Foto1 TYPE RASTER STATUS on OFFSITE 252 252 252 METADATA WMS_NAME "Foto1" WMS_FORMAT "IMAGE/PNG" WMS_SRS "EPSG:4258" END END LAYER NAME "MUNICIPIO" STATUS ON DATA MUNICIPIO LABELITEM "NOM_MUNICI" METADATA "WMS_TITLE" MUNICIPIO END PROJECTION "init=epsg:4258" END TYPE POLYGON CLASS STYLE COLOR 254 226 197 OUTLINECOLOR 255 0 0 END LABEL SIZE SMALL COLOR 0 0 0 ALIGN CENTER REPEATDISTANCE 9 POSITION CC END END END ** And I try to add the layer : **wms = new OpenLayers.Layer.WMS('Municipios', "../mapa2.map &service=wms&mode=map&request=getcapabilities",   {layers:"MUNICIPIO"}, {isBaseLayer:"true",srs: "EPSG:4258", format: 'image/png',visibility: "true", transparent: "true", singleTile: true, ratio: 1 }); ** I'm using mapserver by the way. Anybody can help me?


title: How to pick the coefficient in Inverse Distance Weighted method?

text: I'm using QGIS 1.8.0 to interpolate weather data. This is, I have a weather station located in say, A and C. With the Inverted Distance Weighting or Triangular Interpolation methods, I can know/interpolate information in B. When using the IDW method, I have to choose a coefficient. Does anyone know in which range the coefficients should be? Say from 1 to 10 or .1 to .99 or 100 to 1000? Any ideas? Any information would be really appreciated. Regards.


title: ArcGIS - How to copy and replace polygon (ONLY GEOMETRY)

text: I have two layer (A, B). I need copy polygon from layer A and replace polygon in layer B, but only geometry. I want to preserve original attribute table from layer A.


title: Can I use a cursor to select two rows and compare attributes

text: Can I use the cursor search function to look at the first two rows in a table and select them so I can compare attributes and possibly merge them or change attributes? I am using ArcGIS 10.0


title: How to use SLD to filter a WMS layer but without setting any styles?

text: I'm trying to filter (using OpenLayers) a WMS layer, served by MapServer. I've got this working by building the appropriate OpenLayers.Filter rule(s) with an empty Polygon symbolizer (so no specific styles are set). I then construct an SLD from the rule(s) and pass that in the SLD_BODY parameter. This has been working fine, however I've hit a snag. Originally, the map files for the various layers were setting `STYLEITEM "AUTO"`, which was pretty ugly, but it worked. I've now got rid of this and am manually setting the styles depending on the feature. e.g.               CLASSITEM "SomeField"     CLASSITEM "AnotherField"          CLASS         EXPRESSION ("[SomeField]" eq "SomeValue" && "[AnotherField]" IN "val1,val2")              STYLE             COLOR 255 255 255             OPACITY 25         END #STYLE              STYLE             SYMBOL "hatch"             COLOR 0 200 255             ANGLE 45             SIZE 2             WIDTH 3             GEOMTRANSFORM (buffer([shape],-1))         END #STYLE              STYLE             OUTLINECOLOR 0 200 255             WIDTH 1         END #STYLE     END #CLASS      ...Plus a good few more per layer. Again, this all works fine. The layers look as I expect them too...until a filter is applied. Then, the SLD_BODY gets generated with an appropriate filter and with no style data. However, my various styles defined (in the map files) then get ignored and every feature is just rendered in grey. What am I doing wrong? It seems like the classes (or the expressions in the classes) are being ignored when using a filter. The generated SLD document looks like:               <?xml version="1.0" encoding="UTF-8"?>     <sld:StyledLayerDescriptor xmlns:sld="http://www.opengis.net/sld" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">        <sld:NamedLayer>           <sld:Name>MyLayerName</sld:Name>           <sld:UserStyle>              <sld:FeatureTypeStyle>                 <sld:Rule>                    <ogc:Filter>                       <ogc:PropertyIsEqualTo matchCase="false">                          <ogc:PropertyName>SomeField</ogc:PropertyName>                          <ogc:Literal>SomeValue</ogc:Literal>                       </ogc:PropertyIsEqualTo>                    </ogc:Filter>                    <sld:PolygonSymbolizer>                       <sld:Fill />                       <sld:Stroke />                    </sld:PolygonSymbolizer>                 </sld:Rule>              </sld:FeatureTypeStyle>           </sld:UserStyle>        </sld:NamedLayer>     </sld:StyledLayerDescriptor>      And the appropriate part of the OpenLayers JS:               rule = new OpenLayers.Rule({         filter:     layerFilter,         symbolizer: {'Polygon': {         }}     });      Side note: If I manually change the `<sld:Fill />` to `<sld:Fill></sld:Fill>`, it works. :| Although I could just use JS to replace it, that seems like a really hacky, messy way of doing it and would be a real last resort. What's the proper way of doing this? Screenshots attached to show what I mean regarding styling: **Styles working (No filters)** ![Styles working \(No filters\)](http://i.stack.imgur.com/iNOeN.png) **Styles broken (Filter active for the main feature in the image)** ![Styles broken \(Filter active for the main feature in the image\)](http://i.stack.imgur.com/4GVNp.png)


title: Shapefile rotated after displayed on Bing Maps

text: I took shapefile of country here http://thematicmapping.org/downloads/world_borders.php Then opened in QGIS, after that save needed layer to KML file and took coorditaes for displaying polygons var country = new Microsoft.Maps.Polygon([[new MM.Location(35.360881805420263,46.286518096923942),... In QGIS shape is showing correctly but on Bing Maps it is rotated and in other place: http://clip2net.com/s/5SE7H7 Whats problem? Thank you for helping!


title: Is there a way to automatically calculate attribute values for new features?

text: I am creating a vector point layer which will be added to over time. I have been using the field calculator to populate fields such as X and Y coordinates. I'm wondering if there is a way to have QGIS automatically calculate values in the columns for which a new point is created. Eg. when I put in a new point it automatically fills out the X and Y coordinate fields. Thanks!


title: Google Earth import terrain

text: Currently I am able to import images to overlay on existing terrain in Google Earth. Is it possible to import custom terrain for viewing in Google Earth? There doesn't appear to be any options for doing so?


title: Draw regularPolygon from left to right in OpenLayers rather than center out

text: When I draw a Regular Polygon in Open Layers, it draws from the center out (like a circle). I would like it behave more like Handler.Box does (drawing from the UL corner to LR corner.) Is this possible or do I need to figure out how to get Handler.Box to work? Whenever I use Handler.Box I get this error:                TypeError: feature.geometry.getBounds is not a function      var bounds = feature.geometry.getBounds();      Which is an error in the OpenLayers.js file. I have all the code for Regular Polygon (all the after insert code and such) so I would prefer to get it to behave how I want. I couldn't find anything in the docs about it, or maybe I didn't understand it. Thanks.


title: OpenStreetMap content in commercial applications

text: I was considering the use of the OpenStreetMap API to retrieve vector features and combine them with other source data through a commercial tool. However, OpenStreetMap uses the Creative Commons Attribution-ShareAlike 2.0 (CC-BY-SA) license for its vector feature data. The OpenStreetMap copyright page mentions that: > OpenStreetMap is open data, licensed under the Creative Commons Attribution- > ShareAlike 2.0 licence (CC-BY-SA). > > You are free to copy, distribute, transmit and adapt our maps and data, as > long as you credit OpenStreetMap and its contributors. If you alter or build > upon our maps or data, you may distribute the result only under the same > licence. Does this mean that any item produced with OpenStreetMap data would need to be shared through CC-BY-SA? I do not want to place such restriction on our users, since they may combine OSM data with data that has distribution restrictions. I also wonder if people that use the OSM file import capabilities of several GIS tools are aware of the OpenStreetMap content license.


title: MapServer computes wrong map scales -- off by factor 0.77. Or is it just me?

text: I've noticed in a current project that MapServer seems to get the map scale wrong, or then I'm wrong, but don't see how. My mapfile contains the following definitions:               RESOLUTION 72   # this is the default, but I set it explicitly, just to make sure          PROJECTION         "init=epsg:21781"  # = CH1903 LV 03, which uses meters as coordinate units     END          LAYER         NAME "one"         MAXSCALEDENOM 7500         ...     END          LAYER         NAME "two"         MINSCALEDENOM 7502         ...     END      If I only display these two layers, I should definitely be able to get an empty map at scale 1:7,501. I have written a simple map viewer application to verify this. My viewer also assumes a screen resolution of 72 DPI, just like MapServer. However, it turns out that both layers are invisible around map scale 1:5,803. I would have to adjust my viewer to 93.062 DPI to get the empty map at 1:7,501. Why is this? Just for reference, this is how I calculate the map scale in my viewer:               double imageWidthInPixels = ...;   // passed to mapserv.exe (via 'map_size' arg.)     double extentWidthInMeters = ...;  // = extent's xmax-xmin (via 'mapext' arg.)          // constants:     const double dpi = 72.0;     const double inchesPerPixel = 1.0 / dpi;     const double inchesPerMeter = 1.0 / 0.0254;          double inchesOnScreen = inchesPerPixel * imageWidthInPixels;     double inchesInReality = inchesPerMeter * extentWidthInMeters;          double scaleDenominator = 1 / (inchesOnScreen / inchesInReality);      How could this go wrong? P.S.: I am quite certain that I have the right `extentWidthInMeters` since I can easily derive it from the extent. (I'm using the Swiss projected coordinate system CH1903 LV03, which already uses meters as coordinate units.) * * * **UPDATE:** It seems as if the inaccuracy was partly caused by differing aspect ratios of the map image control in my viewer and the extent sent to MapServer. I've now made sure that only map extents are passed to MapServer that have the same aspect ratio as my map image control. But some small error still remains: The map is now empty around scale 1:7'489 instead of 1:7'501.


title: IPoint -> Nearest Address

text: **Summary: How to use the MDB file (link below) as a locator to reverse geocode (point->address) with ArcMap?** I've implemented this sample: http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#/How_to_find_the_address_closest_to_a_point_using_reverse_geocoding/0001000001n9000000/ The sample works with a file ESRI provided, Atlanta.mdb, and with DataMaps10 US Streets locator. I am trying to get my own MDB files to work. I can do the following with them:               //Address -> Point      //Address -> Closest Address Match      //Point -> Nearest intersection      //Intersection -> Point       **I can't get Point -> Address to work.**   I don't know how to identify the problem. I project my IPoint coordinates to the projection of the MDB, but that doesn't work. I've tried everything I can think of. Can you give me some ideas of things to try to isolate the problem? I have a temporary license for ArcMap for another week, so I can use that to help. I'm using ArcEngine 10, C#, & VS 2010. **Edit:** I've uploaded a zip file here: http://www.filedropper.com/roads It has an MDB, a shapefile and a projection. I used the shapefile to create the MDB file. I need to know how to reverse geocode this MDB file. I need to just one more operation:               //Point -> Nearest Address       I would like to do it programmatically, but if you tell me the steps to do it in ArcMap, then I will have a decent idea how to implement it programatically, and that will be good enough. To summarize: How to use the MDB file as a locator to reverse geocode with ArcMap? **Edit #2:** When I implement the sample here: http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#/d/0001000001n9000000.htm There is this line of code:               IPropertySet addressProperties = reverseGeocoding.ReverseGeocode(point, false /*return intersection */);      If the second parameter, return intersection == true, then it works. But if the second parameter is equal to false, it throws a COM error, 0x80040217, which corresponds to GEOCODING_E_REVERSE_GEOCODE_NO_ADDRESS_FOUND . I've tried a lot of different coordinates and I have my search radius equal to 1 km, so something should be found. **Edit #3:** I found this http://resources.arcgis.com/content/kbase?fa=articleShow&d=38314 Which points to this: 9.3.1 Locator Styles http://resources.arcgis.com/gallery/file/geocoding/details?entryID=12D8D400-1422-2418-34B0-4FE1CC06C0ED I am hoping that I can get these 9.3.1 styles in the same way that the empirical programmer got it to work, and that might be a solution. If I can I'll post it, and give you the points empirical.


title: How to reverse geocode with ArcMap?

text: What steps to I have to take in ArcMap (ArcEditor) to input a lat / long and return nearest address? I found the Reverse Geocode tool but I can't figure out what to put in the "Input Features" and "Output Feature Class" fields or where to input my lat / long.


title: How to style ways with CartoCSS depending on a relation property (OSM)

text: I'm using CartoCSS together with an OSM Postgres DB (hstore on). I want to set the line color depending on the operator tag of relations (routes). How can I do that? **Update** Thanks to Paul Norman for showing how to include tags only in hstore. As Andre Joost points out, for the operator tag that's not necessary, because it's already included with the default osm2psql config. Doing some more research, I could dig deeper into the problem and with some surprise found out that there's a divergence between the documentation on http://wiki.openstreetmap.org/wiki/Osm2pgsql/schema and my DB instance. It says: _planet_osm_line: This table contains all non-closed ways which were imported. Additional rows are created for ways which are members of type=route relations. ... The osm_id is the negative ID of the relation_ While there are rows with negative osm_id's in the DB, I could not find such an entry for 2 randomly tested ways, e.g. for http://www.openstreetmap.org/way/280802857. As you can see, this way is a member of a relation (3731203) fulfilling the mentioned requirements. Now I'm not sure if that's the source of my problem. I don't know if that osm2pgsql doc is up to date and I don't know if mapnik uses this DB entries with negative osm_id to evaluate relation tags. But I would guess so. As a quick fix I manually added a column to the planet_osm_line table with the info I need. But I'd prefer to get rid of that hack asap. Any hints how to proceed?


title: How to automatically generate useful scales for data driven pages in ArcGIS desktop

text: I have been tasked with creating several hundred maps of properties that intersect a long pipeline. The maps will zoom to the part of the property which shows the intersection. The lengths of pipeline that intersect the properties is highly variable -from two metres to over 16,000 metres. Therefore, each map will be at very different scales. I am trying to work out a way to automate the updating of a SCALE field for use in data driven pages. I have found that in general, for my maps, there is between a 1:4 to 1:5 ratio between the length of pipe segment to the scale. That is, if the pipe segment is 1100 metres, then a useful scale would be about 1:5000. I would prefer to round the scale up rather than down so that I can see the full feature on the map. I want to avoid using scales that are not commonly used such as 1:7000 or 1:19,000, for example. I have come up with the following solution for use in the field calculator:               import bisect     def scale(n):       scaleMult=n*5       scaleList=[500,1000,1250,1500,2000,2500,4000,5000,7500,10000,12500,15000,20000,25000,30000,40000,50000,75000,100000,150000,200000]       roundedScale=scaleList[bisect.bisect(scaleList,scaleMult)]       return roundedScale      The above code multiplies the length by some number (I chose 5) then finds out where in my predefined list of scales it should go and returns the next one higher. The code works pretty well for about 90% of the properties. I played around with the scaleMult variable and multiplied n by a few different values between 4 and 5 but nothing really improved the maps. I was thinking that there might be a better solution. If the pipe is not very straight and doubles back on itself slightly or has a major bend in it, then the feature didn't fit well on the map. It was usually zoomed in too far. One idea I had was to use the features bounding box. Maybe there is a relationship between the feature's bounding box area and the scale. But, if the bounding box is long and thin, the scale could be different to a bounding box that was a square of the same area. Does anyone have any ideas on how to implement such a process? Are there any better solutions that I'm not thinking of? Looking for solutions in ArcGIS 10.1 and 10.2. I'll happily edit the question if more clarification is required.


title: GeoServer/GeoWebCache - showing raster layers depending on scale (zoom level)

text: I have a client application, "talking" to my GeoWebCache service (using TMS), showing a map. The map is built out or a few raster layers. I need to show only one of them a time - depending on the zoom level (scale), one of the layers should be shown - the rest not visible. The client can use only one "link" to access the map - I cannot switch the layers on the client. I created a group layer on GeoServer (2.1.4) built out of the two layers I am using (250k and StreetView from OS). I cannot see any way to specify that one of them should be seen when the scale is less then N and the other otherwise. I tried a style for raster layer, specifying opacity depending on the scale, but it does not seem to work. How is it done on GeoServer? If it is doable on GeoWebCache level, that would be fine too. I successfully added the layers to geowebcache.xml file (I need the epsg:27700 projection, that is why the config was needed), but cannot find a way of saying "use this wms url when scale is less than N". Any help will be much appreciated, thank you. ---- UPDATE after the answer - the solution ---- I have two styles created, called raster-0-50k and raster-50k-up. Here is the first one of them, which I applied to the layer at the bottom of the stack (the top one in the list of the group layer):                 <NamedLayer>         <Name>raster-0-50k</Name>         <UserStyle>         <!-- Styles can have names, titles and abstracts -->           <Title>Raster 0 to 50k</Title>           <Abstract>Raster, showing scale 0 to 50000</Abstract>           <!-- FeatureTypeStyles describe how to render different features -->           <!-- A FeatureTypeStyle for rendering rasters -->           <FeatureTypeStyle>             <Rule>               <Name>rule1</Name>               <Title>Opaque Raster</Title>               <Abstract>A raster with 100% opacity</Abstract>               <MaxScaleDenominator>50000</MaxScaleDenominator>               <RasterSymbolizer>                 <Opacity>1.0</Opacity>               </RasterSymbolizer>             </Rule>           </FeatureTypeStyle>         </UserStyle>       </NamedLayer>      the other style does exactly opposite - sets opacity to 1 for scale 50k and higher. In the group layer's list, I associated the two styles with my two layers. It does work now :) Thank you.


title: What are non-commercial sources of road network other than OSM

text: I want to get the graph (Road Network) of **Area of interest** : Lon/Lat - `Left/Bottom = -93.28587,44.9674` and `Right/Top=-93.25351,44.9855`   Minneapolis, MN Searching for a non-commercial source - aside from Open Street Map (OSM) - of a road dataset containing routing and network cost attributes.


title: Why can't I import my script from the python console?

text: i am working with the qgis python console. I usually store all my python scripts into ~/Scripts/python I do have set the PYTHONPATH to ~/Scripts/python in my .bashrc However, I can not manage to import my sripts from the console inside QGIS. Here is my error when I try to import a script called Create_Figures_2D               import Create_Figures_2D     Traceback (most recent call last):       File "<input>", line 1, in <module>       File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 309, in _import         mod = _builtin_import(name, globals, locals, fromlist, level)     ImportError: No module named Create_Figures_2D      How do I set up the PYTHONPATH for pygis ?


title: Join text with layer

text: In Mapinfo it´s easy to link a text-layer with a geometric layer by attributes with outcome as a new layer. Where and how to do this in Q-GIS?


title: User prompt not working

text: I created the following "trial.py":               import arcpy          #Get Name     Name = str(arcpy.GetParameterAsText(0))      Then I created a toolbox using ArcCatalog and then added a script, pointed the above script file. When I run from the ArcMap/Catalog by double clicking on the script name, it wont give a prompt/dialog to enter the user name. It just says "This tool has no parameters". What does "0" in arcpy.GetParameterAsText(0) mean? After getting the user input, I want to pass that parameter to sql querry on shape file, and zoom the map to that location. Would appreciate any ideas how to do that.   * Update1: Do I need to set parameter in the properties sheet: ![enter image description here](http://i.stack.imgur.com/xJ65A.jpg)   * Update2: How to focus the cursor on the text box so that I dont have to click on the text box using mouse. Is there anything in python like tab order in VB.


title: Qgis2.0.1 Zoom Problems

text: I'm working with project with EPSG:20790 in both basemap and SHP. i'm making a road in the middle of a farm between areas of vines, because of the precision that i want to have i'm working with some zoom in. The problem is, when i zoom out, it appears my problem, the road go to the vines area instead of staying in the place where i meant to be. Do you know the solution for this problem?


title: Map algebra for loop

text: I have never coded in Python before and I am trying to loop a simple map algebra equation across 365 rasters. I am trying to calculate total water flow over cells but in the raster I have a negative value to indicate the water table is below ground: for my purposes those cells should read 0. I have tried to write a for loop but it does not seem to be working. Python does not throw any errors, but the output data set still has negative numbers. Any suggestions on why the loop isn't working to remove the negative values?               arcpy.env.workspace="C:/..../2013_q1_tiff_v2prov"     Fcs= arcpy.ListRasters()      From arcpy.sa import *     from arcpy.sa import *     print Fcs          for Raster in Fcs:         Con(Raster <0, 0, Raster)          pleasework= CellStatistics(Fcs, "SUM")


title: Will I now need to have a paid account to add tools/scripts into the new and improved ArcGIS online cloud?

text: I've uploaded a few scripts into the ESRI galleries (which were an update to the legacy ArcScripts) in the past and I just wanted to upload an updated addin someone has asked about... I might be missing something but it seems to me that ESRI managed to make the already convoluted maze of their websites even more complex. It seems that there is a new structure in place whereby I need to create an ArcGIS account (separate from ESRI global account) in order to upload content. The problem is that although there is a 30-day free trial, it seems that I will only be able to upload and manage content by the means of "My Content" for which i need an ArcGIS online paid subscription. Is this correct? Will we now need to pay ESRI to upload our custom content and tools? ![enter image description here](http://i.stack.imgur.com/rSPlz.jpg)


title: Inserting union of two tables into another table Postgresql

text: I am working on a research project. I have two tables, one of a polygon that represents a severe weather risk area, one that represents the US. I want to merge them into one output polygon. I am stuck with the syntax of how to do this. `spc` is the table of the risk area, `simple_us` is the us. I want to output the union of these two table into the geometry column of a new table called `spc_clip`. My current command is:               INSERT INTO spc_clip (geom) VALUES (ST_Union(a.geom, b.geom)), spc as a, simple_us as b;      As it sits now I get a syntax error at spc when I run the command. Any ideas on how to proceed? Next I want to export the geometry of the union to a JSON. I am not sure how to do this. I cannot test this line until I get the above union done.               SELECT row_to_json(t) FROM (SELECT name,geom FROM spc_clip) t;      NOTE: name is already in the table and works fine. I am just struggling with the join and the json export. Thank you very much to anybody who responds!


title: Change the value of an attribute using QgsFeature object

text: I would like to update one attribute of a feature. However, I don't know to update it using the feature object. I have to use data provider to update it.               layers = QgsMapLayerRegistry.instance().mapLayersByName('my_line')      layer = layers[0]     dp = layer.dataProvider()     it = dp.getFeatures()          for i in range(0, dp.featureCount()):       feat = it.next()       attrs = { 2 : 30 }       layer.dataProvider().changeAttributeValues({ feat.id() : attrs })      Can I change the value of an attribute using QgsFeature object? Moreover, is it possible to loop using an iterator object?


title: How to use osm2po to solve TSP

text: I've got a TSP to solve. I created a routable street graph using osm2po and imported the result to PostGIS. On my PG db I installed the pgrouting extension (over PGAdmin3 GUI --> right click Extension-->pgrouting). I've got a second table with a visits_id (places to be visited) and the nearest node_id (node_id from osm2po). My question is how can I solve the TSP which seems to be build into osm2po already? a.) I couldn't find any document how to do it in TSP. b.) Do I have to install an additional function or something in my PostGIS? Or is installing pgrouting via unix-terminal and enabling it in PGAdmin3 is enough? c.) Do I solve the tsp via the osm2po inbuild tsp-algorithm or using pgrouting? Any help would be appreciated!


title: Speeding up OpenStreetMap PostGIS querying

text: I have OpenStreetMap data for the Netherlands loaded into a PostGIS database (PostgreSQL 8.3 / PostGIS 1.3.3) using the osmosis schema. This means all tags are stored in a hstore field. In addition to the GIST index that osmosis creates on the geometry field, I created an additional GIST index on the tags field. Trying to query using both a spatial constraint and a constraint on the tags field, I find that it is slower than I would like. A query like this one:               SELECT n.geom,n.tags,n.tstamp,u.name FROM nodes AS n        INNER JOIN users AS u ON n.user_id = u.id        WHERE tags->'man_made'='surveillance'        AND ST_Within(geom, ST_GeomFromText('POLYGON((4.0 52.0,5.0 52.0,5.0 53.0,4.0 53.0,4.0 52.0))',4326));      takes 22 seconds to return 78 records. There are around 53 million records in this table. Is there a way to significantly speed this up? I've heard that hstore is implemented significantly better in PostgreSQL 9, would upgrading help?


title: Classified Images of RandomForest-Classification look clustered

text: I performed a RF-classification. I have 6 bands, ~250 trainingpoints in a shapefile, OOB estimate of error rate: 1.56% Picture of take trainingpoints added. My code:               setwd("D:/BA-Workspace/sw_comp_ap_aug/")               library(sp)     library(rgdal)     library(raster)     library(randomForest)               # CREATE LIST OF RASTERS     rlist=list.files(getwd(), pattern="tif$", full.names=TRUE)           # CREATE RASTER STACK     xvars <- stack(rlist)            # CREATE RASTER STACK     xvars <- stack(rlist)       x <- coordinates(xvars)[, 1]     y <- coordinates(xvars)[, 2]          x_rst <- y_rst <- xvars[[1]]     x_rst[] <- x     y_rst[] <- y          xvars <- stack(x_rst, y_rst, xvars)     names(xvars) <- c("x", "y", "fitzefa1_rot1", "fitzefa1_gr1", "fitzefa1_bl1",     "fitzefa2_rot1", "fitzefa2_gr1", "fitzefa2_bl1")     # READ Raster TRAINING DATA     sdata <- readOGR(dsn=getwd(), layer="sw_trainshape")          # ASSIGN RASTER VALUES TO TRAINING DATA     v <- as.data.frame(extract(xvars, sdata))     sdata@data = data.frame(sdata@data, v[match(rownames(sdata@data), rownames(v)),])          sdata@data  <- sdata@data[-c(5,6)]           # RUN RF MODEL     rf.mdl <- randomForest(x=sdata@data[,3:ncol(sdata@data)],     y=as.factor(sdata@data[,"Class"]),                        ntree=501, importance=TRUE)          # CHECK ERROR CONVERGENCE     #plot(rf.mdl)     #sdata@data           # PLOT mean decrease in accuracy VARIABLE IMPORTANCE          #varImpPlot(rf.mdl, type=1)     #setOldClass(SpatialPointsDataFrame)     # PREDICT MODEL     predict(xvars, rf.mdl, filename="RfClassPred.img", type="response",          index=1, na.rm=TRUE, progress="window", overwrite=TRUE)      Here the classified .img ![enter image description here](http://i.stack.imgur.com/arzOD.jpg) The classes are obviously wrong. Could it be due to wrong training-points? I have 7 classes. (No excluded class for shadows or so, because I'm just trying to devide some tree-species) ![enter image description here](http://i.stack.imgur.com/nrNja.jpg) ![enter image description here](http://i.stack.imgur.com/TsqgH.jpg) ![enter image description here](http://i.stack.imgur.com/eOz46.jpg) ![enter image description here](http://i.stack.imgur.com/zIEAD.jpg)


title: tabbed popup window for Arcgis Javascript

text: Anyone know examples of a tabbed popup window for Arcgis for javascript, NOT an infowindow? If you could point me in the the right direction it would help.


title: OpenLayers data moves when panning map

text: I have an OpenLayers map. On to this I have loaded a GeoJSON file. When I pan the map the GeoJSON data moves relative to the basemap. Surely this isn't right! Can anyone shed some light no why this is happening?


title: WFS Spatial Operator Query - Crosses AND Within

text: I'm trying to query my wfs service with a polygon. I want the service to return both features that are inside the polygon and partially inside the polygon. However, when I use the "Within" spatial operator, I only get features that are completely inside. When I use the "Crosses" operator, I only get the features that are partially inside the polygon. I have tried using And/Or to query for features that are either one or the other, but WFS does not support multiple spatial operators in a query. It gives me this error: > "Unsupported filter - filter has more than one spatial operator." Does anyone have an idea of how to get around this without sending two separate ajax requests to the service? Here's an example that returns data with the "within" filter.               var xml = '<ogc:Filter xmlns="http://www.opengis.net/ogc">\n'                  +      '<ogc:Within>\n'                  +        '<ogc:PropertyName>Shape</ogc:PropertyName>\n'                  +        '<gml:Polygon>\n'                  +          '<gml:outerBoundaryIs>\n'                  +            '<gml:LinearRing>\n'                  +              '<gml:coordinates>' + queryObject.geometryString + '</gml:coordinates>\n'                  +            '</gml:LinearRing>\n'                  +          '</gml:outerBoundaryIs>\n'                  +        '</gml:Polygon>\n'                  +      '</ogc:Within>\n'                       +    '</ogc:Filter>\n'      Here's the capabilites of the wfs service:               <ogc:SpatialOperators>     <ogc:SpatialOperator name="BBOX"/>     <ogc:SpatialOperator name="Equals"/>     <ogc:SpatialOperator name="Disjoint"/>     <ogc:SpatialOperator name="Intersects"/>     <ogc:SpatialOperator name="Crosses"/>     <ogc:SpatialOperator name="Touches"/>     <ogc:SpatialOperator name="Within"/>     <ogc:SpatialOperator name="Contains"/>     <ogc:SpatialOperator name="Overlaps"/>     </ogc:SpatialOperators>


title: OpenStreetMap Neighborhood Boundaries

text: Is it possible to download neighborhood boundaries from OpenStreetMap? I see on maps that they often indicate neighborhood names; is there a way to retrieve their coordinates and boundary polygon coordinates, much like the Zillow Neighborhood Boundaries data set?


title: How to get GDAL to create statistics for GTiff in Python

text: I regularly create my own GeoTIFF rasters with GDAL in Python, e.g.:               from osgeo import gdal     from numpy import random     data = random.uniform(0, 10, (300, 200))     driver = gdal.GetDriverByName('GTiff')     ds = driver.Create('MyRaster.tif', 200, 300)     band = ds.GetRasterBand(1)     band.WriteArray(data)     ds = band = None # save, close      however when the result is viewed with ArcCatalog/ArcGIS, it looks either black or grey, since it has no statistics. This is solved either by right- clicking the raster and choosing "Calculate Statistics..." in ArcCatalog (there are several other ways to do this), or using gdalinfo in a command prompt:               gdalinfo -stats MyRaster.tif      will generate `MyRaster.tif.aux.xml`, which is used by ArcGIS to properly scale the raster. The PAM (Persistent Auxiliary Metadata) file contains the statistics, most notably the minimum and maximum values:               <PAMDataset>       <PAMRasterBand band="1">         <Metadata>           <MDI key="STATISTICS_MINIMUM">0</MDI>           <MDI key="STATISTICS_MAXIMUM">10</MDI>           <MDI key="STATISTICS_MEAN">5.0189833333333</MDI>           <MDI key="STATISTICS_STDDEV">2.9131294111984</MDI>         </Metadata>       </PAMRasterBand>     </PAMDataset>      My question: is there a built-in way of getting GDAL to create a statistics file (other than using the `gdalinfo -stats` command)? Or do I need to write my own?


title: QGIS 2.4: adding feature with setLineEdit() data not saving to database

text: A while ago, I struggled to get a combo box working correctly for QGIS 2.0. This involved transforming a combo box into a line edit if the item named "other" was chosen, enabling the user to add a value not already existing in the combo box. Now that i've updated to **QGIS 2.4 using Windows XP** , however, I am having other issues! The combo box still turns into a line edit, but the user-input does not get saved to the database. After the line edit appears, and a new value is added, "Autre" is still the value saved to the database.               #wiring up the signal to my function     Diametrefield.currentIndexChanged.connect(autre)          #changing the combo box to line edit     def autre():         global Diam_autrefield, Diametrefield         if Diametrefield.currentText() == '..Autre':             QMessageBox.warning(None,'Note de saisie',"Definir 'Autre' Diametre")             Diametrefield.setLineEdit(Diam_autrefield)          def validate():         global Diametrefield         if Diametrefield.currentText() == '..Autre':             Diametrefield.setStyleSheet("background-color: rgba(255, 107, 107, 150);")             QMessageBox.warning(None, 'Probleme de saisie',"Veuillez choisir 'Diametre'")         else:             # Return the form as accepted to QGIS.             myDialog.accept()


title: How are .mif files loaded to QGIS?

text: I have downloaded sources of Qgis 1.8. I would like to know how it accepts mif file. As I have learnt after choosing add VectorLayer. Method addVectorLayers is called. Then openvectorlayer dialog shows. Afterwards QgsVectorLayer is created. It is constructor               QgsVectorLayer::QgsVectorLayer( QString vectorLayerPath,                                     QString baseName,                                     QString providerKey,                                     bool loadDefaultStyleFlag )         : QgsMapLayer( VectorLayer, baseName, vectorLayerPath )         , mUpdateThreshold( 0 )     // XXX better default value?         , mDataProvider( NULL )         , mProviderKey( providerKey )         , mEditable( false )         , mReadOnly( false )         , mModified( false )         , mMaxUpdatedIndex( -1 )         , mActiveCommand( NULL )         , mRenderer( 0 )         , mRendererV2( NULL )         , mUsingRendererV2( false )         , mLabel( 0 )         , mLabelOn( false )         , mVertexMarkerOnlyForSelection( false )         , mFetching( false )         , mJoinBuffer( 0 )         , mDiagramRenderer( 0 )         , mDiagramLayerSettings( 0 )     {       mActions = new QgsAttributeAction( this );            // if we're given a provider type, try to create and bind one to this layer       if ( ! mProviderKey.isEmpty() )       {         setDataProvider( mProviderKey );       }       if ( mValid )       {         // Always set crs         setCoordinateSystem();              QSettings settings;         //Changed to default to true as of QGIS 1.7         //TODO: remove hack when http://hub.qgis.org/issues/5170 is fixed     #ifdef ANDROID         bool use_symbology_ng_default = false;     #else         bool use_symbology_ng_default = true;     #endif         if ( settings.value( "/qgis/use_symbology_ng", use_symbology_ng_default ).toBool() && hasGeometryType() )         {           // using symbology-ng!           setUsingRendererV2( true );         }              // check if there is a default style / propertysheet defined         // for this layer and if so apply it         bool defaultLoadedFlag = false;         if ( loadDefaultStyleFlag )         {           loadDefaultStyle( defaultLoadedFlag );         }              // if the default style failed to load or was disabled use some very basic defaults         if ( !defaultLoadedFlag && hasGeometryType() )         {           // add single symbol renderer           if ( mUsingRendererV2 )           {             setRendererV2( QgsFeatureRendererV2::defaultRenderer( geometryType() ) );           }           else           {             QgsSingleSymbolRenderer *renderer = new QgsSingleSymbolRenderer( geometryType() );             setRenderer( renderer );           }         }              mJoinBuffer = new QgsVectorLayerJoinBuffer();              connect( QgsMapLayerRegistry::instance(), SIGNAL( layerWillBeRemoved( QString ) ), this, SLOT( checkJoinLayerRemove( QString ) ) );         updateFieldMap();              // Get the update threshold from user settings. We         // do this only on construction to avoid the penality of         // fetching this each time the layer is drawn. If the user         // changes the threshold from the preferences dialog, it will         // have no effect on existing layers         // TODO: load this setting somewhere else [MD]         //QSettings settings;         //mUpdateThreshold = settings.readNumEntry("Map/updateThreshold", 1000);       }     } // QgsVectorLayer ctor      So for now I need to know where the data from .mif file is been parsed.


title: arcpy.Delete_management not deleting

text: I'm struggling with deleting an mxd in a for loop. My code is:                   import arcpy, glob, os, datetime, time         from time import strftime              for filename3 in glob.glob(os.path.join(folderPath, "*.mxd")):              fullpath = os.path.join(folderPath, filename3)         basename, filename2 = os.path.split(fullpath)              mxd = arcpy.mapping.MapDocument(filename3)              if os.path.isfile(fullpath):             basename, filename2 = os.path.split(fullpath)                       (shortname, extension) = os.path.splitext(filename2)              mxd.saveACopy(draftloc + "\\" + "SentToDraftOn(" + strftime("%d%b%y %Hh%Mm%Ss") + ")"+ "_" + shortname + extension)         arcpy.Delete_management(mxd)      I'm trying to do some processing on each mxd in a directory and then I want to saveACopy to a new location and then delete the old mxd. I keep getting the following traceback error: > Object: Error in executing tool I think it's an indentation thing, but I can't seem to figure it out. Can anyone suggest what I am missing? Thanks, Mike


title: How to calculate raster statistics for polygons?

text: how can I calculate sums, averages etc of raster-points (multi-band) per polygon of a vector-layer. I was told that this is called "zonal statistics". I tried that with QGIS first. There is a way to do it but that is much too slow (convert raster to vector, intersect with second vector-layer, calculate geometry, export numbers, calculate statistics with spreadsheet or other program, re-import the results, takes forever for my 350.000 raster-points). I was also given the hint to use saga-gis. That has "zonal statistics" but those are based on categories from a raster-band, not on polygons from a vector-layer. So to use this I would have to convert my vector-layer to raster and then calculate the statistics. This seems to be the wrong way to solve this. There would be no way to account for raster-points that belong to 2 or more polygons because they are intersected by the polygon-boundary. I assume that polygon-based statistics should be able to handle this so I also assume that I haven't found the correct module yet. Saga-gis has really many modules. Please let me know which one is the right one for this application. THX, stn


title: Can't Seem to Replicate the KML Layer in jsFiddle

text: I can't see why I am unable to replicate the KML Layer Example: http://openlayers.org/dev/examples/sundials.html in jsFiddle: http://jsfiddle.net/JPN4C/ I'm assuming that it's the map protocol or something else near the url, but can't see what that is.


title: How to load WFS in jsFiddle?

text: I'm trying to implement a WFS service in jsFiddle for an OpenLayers map. It works fine here: http://openlayers.org/dev/examples/cql-format.html My example of it is (the map loads but the filter doesn't): http://jsfiddle.net/L7sFY/ But I can't seem to get the filter to work properly. WFS: http://en.wikipedia.org/wiki/Web_Feature_Service XML WFS file: http://demo.opengeo.org/geoserver/wfs?service=WFS&version=1.0.0&request=GetFeature&typename=topp:states&outputFormat=json&format_options=callback:loadFeatures


title: javascript function to search all OpenLayers vector features covering a pixel

text: I have an OpenLayers map that draws features in a vector layer. The features are selectable and have a popup on select. Unfortunately in a lot of cases the features overlap so it can be impossible to select some features. I think that what I need to do to fix this is change my select control so that it uses a click handler and searches the map for features at this point. What sort of function do i need to write? are there any examples of this being implemented before? This is how the features are drawn:               var vector_Layer = new OpenLayers.Layer.Vector();          function GetFeaturesFromKMLString (strKML) {         var format = new OpenLayers.Format.KML({         'internalProjection': new OpenLayers.Projection("EPSG:900913"),         'extranalProjection': new OpenLayers.Projection("EPSG:4326")             });         return format.read(strKML);     };          vector_Layer.addFeatures(GetFeaturesFromKMLString('$newkml'));      And this is how Layers are currently selected:               var select = new OpenLayers.Control.SelectFeature(vector_Layer, {clickout: true});              vector_Layer.events.on({                 "featureselected": onFeatureSelect,                 "featureunselected": onFeatureUnselect});                             map.addControl(select);             select.activate();                  select.handlers['feature'].stopDown = false;              select.handlers['feature'].stopUp = false;      here is the click event listener I was planning on using:               OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {                                 defaultHandlerOptions: {                     'single': true,                     'double': false,                     'pixelTolerance': 0,                     'stopSingle': false,                     'stopDouble': false                 },                      initialize: function(options) {                     this.handlerOptions = OpenLayers.Util.extend(                         {}, this.defaultHandlerOptions                     );                     OpenLayers.Control.prototype.initialize.apply(                         this, arguments                     );                      this.handler = new OpenLayers.Handler.Click(                         this, {                             'click': this.onClick,                             'dblclick': this.onDblclick                          }, this.handlerOptions                     );                 },                       onClick: function(evt) {                     //function that seachers for and selectas features at this point                 },


title: Map algebra expression with dynamic raster layers

text: Attached is the model that i am working on. ![Model1](http://i.stack.imgur.com/iaBFm.gif) The "To Value" and "From Value" are model parameters and hence subject to change each time according to requirements, This leads to varying number of files generated each time. This variation in the number and names of output files in a workspace is creating the problem in hardcoding the following calculator expression in a subsequent model directly. ![Model2](http://i.stack.imgur.com/KFVR0.gif) This is the expression in the second map algebra tool that i want to use. The input rasters are those generated from the model shown above. Raster Calculator (output1+(output2/Output1)+(output3/Output2)+(output4/Output3)+(output5/Output4)+....) I saw some effort to a similar question here, Raster Calculator dynamic map algebra expression, but I didn't find it conclusive. (I didn't understand the advice to use the FOR loop there) Help and advice please! I am using ArcGIS-10.1


title: SPOT NDVI time series

text: I have 423 spot.vgt ndvi from first April 1998 till 31 Dec 2009. I want to know how to stack these NDVI to have one image for time series analysis? Or better to ask how is it possible to analyze a pixle change during these period?


title: Batch Aggregate Raster using Python

text: Trying to write a python script to batch Aggregate a bunch of different rasters. get a Parsing Error Syntax error for line 11.               # Batch Aggregate raster tool     import arcpy     from arcpy import env       from arcpy.sa import *          env.workspace = r"C:\rasters"     out_workspace = r"C:\rasters\agg"          rasters = arcpy.ListRasters()          for raster in rasters         outAggreg = Aggregate(raster, 10, "MEAN")         outAggreg.save(out_workspace, raster + "_agg")


title: Make InfoTemplate moveable but still highlight attributes

text: I’m using the esri/InfoTemplate in the Javascript API to identify features. However, it gets anchored to the mouse click and cannot be moved around. Based on this answer I used dojo/dnd/Moveable to make it moveable. However, now I cannot highlight or copy the attributes inside of it. It’s like Moveable turns it into a graphic or something. Is there a way to make the InfoTemplate floating but still be able to highlight and copy the attributes? Thanks in advance for your help. (The other option I’m considering is using the jQuery Dialog Widget, but I cannot get it to highlight the geometry for multiple features. Any suggestions there would be great too. Whichever is easier.)


title: Looking for GIS Data for Curitiba, Brazil?

text: I am wondering if anyone knows where to find building footprints, road network, transit, or environment feature data for Curitiba, Brazil. Thanks


title: Drawing Graphic at street layer/level only

text: I am new to ArcGIS JavaScript API and would like some help. I have a scenario where a user has to draw a graphic, mainly a polyline and some points on a map. To avoid confusion I want to limit the user to zoom in to street level prior to start drawing. What I am trying to achieve here is when the user clicks on the locality he/she wants to draw on, on the click even the map zooms in to street layers. I have been looking around but can't find a way to do it.


title: Why do I get a ". not nullable field error" when attempting to split a feature?

text: I'm trying to split feature, by point, into 2 new features.                var pFeatureEdit = lineToSplit.Feature as IFeatureEdit;      var set = pFeatureEdit.Split(point);      When this feature has configured relation, I got exception. Got next COM error: The field is not nullable [FIELDID] Where FIELDID - not nullable foreign key in database Anybody know how to solve this issue? UPD. When I split this feature using standard Split Tool in ArcMap - it splits successfully. This uses another method to split features?


title: Qgis Plugin development with postgres data

text: can someone help me with initial developing on qgis plugin that displays postgres data table on a `QTTableView` or `QTTableWidget`? So far:   * I can connect to database using psycopg2 or using QtSql(tested inside QGIS console)   * Using idle i can't because it does not load the QPSQL drivers...although i have them on my qt binary install.   * I'm using pluginbuilder for initial setup, but i don't know how to get a reference to my tablewidget, to pass data from sql query.   * Using IDLE i can draw a form with a tablewidget that is filled with data from sel query. I already read this question :(but still have doubts...) I need help putting things together. Thanks for reading. _UPDATE:_ In here if i run from terminal i can fill a `qtablewidget` with sql data: http://pastebin.com/Lm6GsUMk I want to draw a widget like this inside a qgisplugin. `My Plugin(TestPlugin)` was buit with plugin builder, asnd inside my `ui_testplugin.ui`, i have a `QtTableWidget` draw with `QtDesigner` and compiled using a generated makefile. I know that i have to insert logic inside `testplugin.py`, inside run block, but i don't know to get a reference to `QtTableWidget`. Hope this makes it clear. Thanks.


title: Qgis map composer shows a huge red cross

text: I can show this best by a picture: ![Figure](http://i.stack.imgur.com/twNl3.png) This cross also shows when I export it as raster or pdf. What is it and how can I get rid of it? I am using Qgis 2.4 and the map was done with the help of the Grass plugin `AniMove` on Linux Mint.


title: OpenLayers Map with custom SRS (without EPSG code)

text:
Is it possible to show geodata by OpenLayers in custom projection? I want to show my geodata in OpenLayers interface. Source layer(s): shape-file(s), standard epsg:4326 lon/lat projection; processed by Mapserver (WMS Layers). Result map: Albers Conic Equal Area projection, with specific central meridian (It has no EPSG code). Definition in PROJ.4 format:               +proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=93 +x_0=16500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs      As to Mapserver Layers – they are all well, I can show geodata through standard Mapserver CGI template-based interface, in this custom Albers projection… For OpenLayers I tried to define Proj4js.defs –                         <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>     <link rel="stylesheet" href="ol/style.css" type="text/css" />
OpenLayers Test

                It seems that the presented coordinate transformation works correctly (displaying the appropriate numbers when I move the mouse), but the map does not appear… How can I do this? # ==================== **Discussion of the problem & Answer to the question --**   It looks like I found the answer myself...   (The above source code has been corrected, and already has a working version) I changed WMS-layer upon layer of MapServer type Original (non-working) version was:               ...         var rus_general = new OpenLayers.Layer.WMS( 'Russian Regions',              'http://localhost/cgi-bin/mapserv.exe?map=C:/OSGeo4W/apache/htdocs/test/z-albers.map',          {layers: 'rus_general'}             );     ...      In this case, the bottleneck - WMS-service: I do not know how to specify a custom (non-EPSG) projection in this case ...   This is not a problem for the Mapserver layer type. Mapserver generates a map projection with PROJECTION section in the MAP-file. There are no restrictions such as "EPSG Only"... Thank you all.


title: QGIS python plugin, segmentation fault, what can I do?

text: I'm trying to write a python plugin to create an input file for an urban climate model out of three shape files. Unfortunatelly QGIS closes when I try to load the plugin. When I start QGIS from the Terminal (Ubuntu), I get "segmentation fault" when it happens. Thats all. Is there some way to prevent QGIS from crashing or another way to track down the error?


title: ArcGIS Server won’t release SDE table locks

text: My ArcGIS Server services are locking tables in SDE. This is a problem if I need to edit or overwrite them. Unfortunately the locks don't go away upon stopping the services. The only solution I have found is to restart the Linux `arcsde` service, which has many implications. Is there a better way to lose these phantom locks?


title: Unable to inspect GML file

text: I use FME to make GML file based on SQLServer data. When I create GML, I cann't see anything in FME Data Inspector or Universal Viewer. Like something is missing in file. I have schema, and XMLValidator validates file against schema, but I still cann't see graphical data. Link to Result.gml Link to schema. Does anyone knows what is the problem and what is missing in schema? **EDIT:** As requested! ![FME workspace](http://i.stack.imgur.com/xT1e0.jpg) This is just a part, but main part. All features come into XMLTemplater, like on the picture. Some of them (with non-point geometry) goes first in GeometryExtractor and then into XMLTemplater. These StringReplacers works fine. XMLValidator validates document successfully.


title: Website, tool or app to produce colored map of the United States?

text: Does anyone know of a website, Windows program or iOS app that has a clickable map of the United States and colors in the states that are clicked on? Looking to produce a map of states I have visited or lived in. Ideally it would be nice to be able to save and print this image as well. Thanks.


title: Where is Sextante in QGIS?

text: I am using QGIS 1.7.3 (can't get 1.8 to work on my WIN 7 64 OS) and I cannot find the Sextante plugin. I have found other plugins with no problems, accessing all the various repositories, I have a list of many to choose from several repositories but there is no Sextante although the Sextante website says it should be downloadable from the QGIS main repository. What am I missing, doing wrong?


title: When iterating over a vector layer with pyqgis, how do I check whether a feature is selected?

text: While iterating over a vector layer using the following code (summerized from the example in the pyqgis cookbook), is there a way to check whether a feature is selected?               provider = vlayer.dataProvider()     feat = QgsFeature()     allAttrs = provider.attributeIndexes()     provider.select(allAttrs)     while provider.nextFeature(feat):         geom = feat.geometry()         attrs = feat.attributeMap()         for (k,attr) in attrs.iteritems():             print "%d: %s" % (k, attr.toString())      Alternatively, I could create a list of selected features using `vlayer.selectedFeatures()`, but I'm hoping there is a way to check each feature directly.


title: How to write XML attributes from feature with ArcView license?

text: I need to generate xml of attribute data of a feature, but I have only Arcview license so, I can't use Export XML tool. What are my options? What is wrong with this code (after applying search cursor writing the attributes of node):                for row in rows:       order.attrib['name'] = row.getValue('Id')       order.attrib['UTM'] = row.getValue('UTM')      I am not able to see any addition to order tag. missing something... Thanks Ok - just realized there was some issue with the feature: Here is the code and error message after the change (So something is broken fundamentally):               >>> import xml.etree.ElementTree as xml     >>> import arcpy     >>> fclass= arcpy.SearchCursor(r"..\SiteMarker_Project")     >>> order = xml.Element('Order')     >>> site = xml.Element('Site')     >>> order.append(site)     >>> rows = fclass     >>> for row in rows:       order.attrib['name'] = row.getValue('Id')       order.attrib['UTM'] = row.getValue('UTM')          Traceback (most recent call last):       File "<pyshell#8>", line 3, in <module>         order.attrib['UTM'] = row.getValue('UTM')       File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 1011, in getValue         **return convertArcObjectToPythonObject(self._arc_object.GetValue(*gp_fixargs(args)))     RuntimeError: ERROR 999999: Error executing function.**      Thanks


title: Estimating Landsat 7 azimuth (heading)

text: I am trying to understand if the Landsat 7 satellite ETM+ orbit has any significant azimuth or if it could be assumed as 0° (as it is not a pointable sensor). I searched in the web but could not find any good answer so far. Zenith is supposed to be equal zero as Landsat images are always nadir. What about the azimuth? Would knowing the scene center coordinates (they're 46.036323 N and 10.456876 E) be of any help? The image looks like a sort of isosceles trapezoid, with the major basis towards the North (as I would expect for the Northern Hemisphere). It is also tilted differently from south to North (as I also would expect). I would appreciate any help, thanks in advance. **UPDATES** One question arised while I was trying to solve my task. The metadata file of Landsat provides coordinates of the four image corners. Anyway they refer to the black background. The azimuth is intended to be calculated through these corners or the corners of the "real"-non-background image? And why do I have to use latitude rather than longitude? Is it not better to use longitude? I thought to use (for instance) the right corners, subctracting the longitudes of them (e.g. upper right long-lower right long) and obtain what I was looking for in terms of degrees, does it make sense for you? Here is an image showing a part of Landsat paths in a Worldwide Reference System (descending orbit) and my Landsat scene. RS is geographic, WGS84. According to this what I should looking at is the non-background image. ![Landsat in the WRS path](http://i.stack.imgur.com/sp8Hu.jpg)


title: Why is a Shape polygon file not editable?

text: I received a .shp polygon file. When I choose the edit toggle tool in QGIS the little red X's appear as expected. Then when I try to edit the polygons such as add a polygon all seems to go well, when finishing with a new polygon the dialog box for the attributes appears which allows entries and then upon closing the new polygon disappears. Also when trying to use the node tool it looks like a node can be moved but not really. Or any other editing such as editing the attribute by selecting a polygon with the information tool, the usual table appears and not the dialog table which allows entries into the attribute variables for the selected polygon. The file is not read only, other files are in the same directory which are editable. Any ideas what the difficulty is?


title: insert cursor "Row contains a bad value"

text: I'm having a problem copying data from a shapefile to a feature class using an insert cursor. Here is the part that is causing trouble:               newRows = arcpy.InsertCursor(Target)     sourceRows = arcpy.SearchCursor(Christiansburg)          for source in sourceRows:         neue = newRows.newRow()         neue.setValue(TName, source.getValue(CName))         newRows.insertRow(neue)      It takes data from a field in the source row, and then copies it to the new row. The trouble is, when the source row has a forward slash in it, it fails when I try to do the insert, saying "The row contains a bad value.". If I edit the forward slash out of the source data, it works fine. If I replace the variable with a hardwired string that has a forward slash in it, it also works fine, so I am not 100% sure that it is the forward slash causing the problem. I have tried using the str function on the result of getValue, and it still doesn't work. The field being copied from is street names, and the row causing trouble is a street that changes name. If worst comes to worst, I could simply edit out the forward slash permanently, but I would rather not if I can help it. Does anyone know what I am doing wrong here? EDIT: Derped on copy/paste of code snippet. Should be fixed now.


title: Color gradient in symbology in ArcGIS 10.1 for interpolation raster

text: I have measured noise and want to create with this measured points a noise map with interpolation....is there any way in symbology to create a flowing transition between the different colors for the created interpolation raster? I used color ramp, so the colors have a gradiation but the borders between these colors are not fluent or in simple words: i want to have a result like a rainbow, where the colors are fluent but without recognizable borders thanks in advance!


title: Basic arcpy script with cursors running SLOW!

text: I've got 2 polygons layers (no gaps): IDUs and crops. What I want is to identify the 3 largest crops by % area for each IDU. This seems common enough but I couldn't find a way to do it all using tools, so I decided to try and learn basic arcpy to finish the job. I'm in ArcGIS 10.1 SP 1 (Advanced) and have done lots of PHP but this is my first try at python. I found the Tablulate Intersections tool which works great and gives a table with all the overlapping crops with area and %area fields, but I want to add attributes to the IDU layer for the 3 largest crop codes - eg. LU_1=123, LU_2=201, LU_3=132. Here's the script I wrote. It works, but it's taking 100sec to process 160 test IDU records, and eventually it needs to process 100,000. I've got a 2013 Retina MacBookPro which flies through everything else I've thrown in in Arc so far, so something is obviously poorly done in my code.               import arcpy     idusFC = "c:/Envision/StudyAreas/PEI-Dev/NewIDUs/PEI.gdb/IDUs_Test"     cropTB = "c:/Envision/StudyAreas/PEI-Dev/NewIDUs/PEI.gdb/IDU_CROP2011_table_Sort1"     iduFields = ["CROP2011_A_Code", "CROP2011_B_Code", "CROP2011_C_Code", "IDU_ID"]     cropsFields = ["CROP_2011_GRIDCODE"]     print "Starting processing"     start = time.time()     with arcpy.da.UpdateCursor(idusFC, iduFields) as idus:         for idu in idus:             itr = 0             iduID = idu[3]             whereClause = arcpy.AddFieldDelimiters(cropTB, "IDU_ID") + " = " + str(iduID)             with arcpy.da.SearchCursor(cropTB, cropsFields, where_clause=whereClause) as  iducrops:                 for iducrop in iducrops:                     idu[itr] = iducrop[0]                     itr += 1                     if itr > 2:                         break             idus.updateRow(idu)     end = time.time()     processTime = end - start     print "Processing Completed - {:.1f} seconds".format(processTime)      Is there another way I should do this so that it runs reasonably fast? Thanks!


title: Adding layer into ArcMap via Python - .pyt versus .py

text: In the final step of a script I've built, I'm trying to add a layer into an empty blank map document. This is my block of python code below that adds the layer. Everything works fine if its run from a script tool pointing at a .py from ArcToolbox, but if I run the same block of code inside a .pyt toolbox, the layer flashes and then goes away. Is there a parameter to set somewhere that I'm missing? I would like to use the .pyt tool because it makes it easier to distribute to users. -"layername" comes from previous part of script               outlayer = layername[:-4]           arcpy.MakeFeatureLayer_management(layer, outlayer)     lf = arcpy.mapping.Layer(outlayer)     mxd = arcpy.mapping.MapDocument("CURRENT")     df = arcpy.mapping.ListDataFrames(mxd, "*")[0]     arcpy.mapping.AddLayer(df,lf,"TOP")


title: Custom C# tool populate field combobox parameter from layer parameter

text: I'm creating a custom tool for arcmap toolbox where user select FeatureLayer in one combobox, and in next combobox user can select a field from selected FeatureLayer. FeatureLayer combobox fills automaticaly from:               ...     inputParameter.DataType = (IGPDataType)new GPFeatureLayerTypeClass();     inputParameter.Value = (IGPValue)new GPFeatureLayer();      The problem arises when I need to fill field combobox. I'm not realy sure how to set up its parameters. right now i'm trying something like this:               inputParameter = new GPParameterClass();     inputParameter.Enabled = true;     inputParameter.Name = "FieldTrg";     inputParameter.DisplayName = "Feature Class Fields";     inputParameter.Direction = esriGPParameterDirection.esriGPParameterDirectionInput;     inputParameter.ParameterType = esriGPParameterType.esriGPParameterTypeRequired;     inputParameter.DataType = (IGPDataType)new GPStringType();     inputParameter.Value = (IGPValue)new GPString();     //inputParameter.AddDependency("FeatureClassTrg");     //IGPFeatureSchema fgs = new GPFeatureSchemaClass();     //fgs.FieldsRule = esriGPSchemaFieldsType.esriGPSchemaFieldsAll;     //inputParameter.Schema = (IGPSchema)fgs      In UpdateParameters function I were able to get fields in selected FeatureLayer/FeatureClass, but I had no idea where and how should I put them back in combobox. Also currently the field 'combobox' isn't really a combobox, but just a string field, so probably GPString isn't corect thing to use as its DataType/Value. Any help regarding this problem would be appreciated.


title: How do I write a script to list those layers with broken data sources?

text: In my SHP to GDB toolbox, I create a list of layers in MXD, then I check for those not broken I move them. Despite using ListBrokenDataSources I have errors like this: > Runtime error : [Error 3] Impossible trovare il percorso specificato: > u'J:\20120208_Backup_disco_E\Archivio\ATLArete\Shape Atlarete > 1-1-2008\miste.shp' **Translate:** > Runtime error : [Error 3] Unable to find the specified path: > u'J:\20120208_Backup_disco_E\Archivio\ATLArete\Shape Atlarete > 1-1-2008\miste.shp' Here is my code:               # Create a list of SHP in mxd          list = []         df_list = arcpy.mapping.ListDataFrames (mxd)              for df in df_list:              for fc in arcpy.mapping.ListLayers (mxd, "", df):                 if (fc.supports("DATASOURCE")) and (fc.dataSource.endswith(".shp")):                             list.append(fc)         print("List created")              # Check for borken datasource         for fc in list:             if arcpy.mapping.ListBrokenDataSources(fc):                 list.remove(fc)         print ("Broken data removed from list")              # Move fc to gdb         add_list = []         for fc in list:             fc_date = time.strftime("%Y%m%d", time.gmtime(os.path.getctime(fc.dataSource)))             try:                 dsc = arcpy.Describe(fc)                 out_temp = '_' + fc_date + '_' + fc.name + '_' + fc_prj                 out_name = arcpy.ValidateTableName(out_temp)                 print ('Created ' + out_name + ' to ' + gdb_full_path)             except:                 out_temp = '_' + fc_date + '_' + fc.name + '_' + 'UNKNOWN'                 out_name = arcpy.ValidateTableName(out_temp)                 print ('Created ' + out_name + ' to ' + gdb_full_path)             arcpy.FeatureClassToFeatureClass_conversion(fc, fd_path, out_name)             replace_name = fd_path + '\\' + out_name             fc.replaceDataSource (gdb_full_path, "FILEGDB_WORKSPACE", out_name)             print ('Replaced data source ' + out_name + ' in ' + gdb_full_path)                   # Create a new MXD         new_mxd = os.path.splitext(gdb_full_path)[0] + '\\' + 'MXD' + '\\' +  time.strftime("%Y%m%d", mxd_date )+ '_' + (mxd_name)         mxd.saveACopy(new_mxd)         print ('Created new MXD ' + new_mxd + ' in ' + ws)