title: Alternative to Snap (Editing) in QGIS/MapInfo?

text: Is there any alternative to Snap (Editing) in either QGIS or MapInfo (user- made tools), or another way to do this in ArcGIS without an advanced license?![enter image description here](http://i.stack.imgur.com/EPulk.png)


title: Where is my Polygon?

text: I created this script to merge all of the polygons in a shapefile into a single polygon. I am using Python 2.7.5, Shapely 1.3.0, and Fiona 1.1.5.               def merge(path_in, path_out):         #Import Libraries         import fiona         from shapely.geometry import mapping, shape, MultiPolygon, Polygon         from shapely.ops import cascaded_union         from fiona import collection              #Import shapefile         with collection(path_in, "r") as input:             input_driver = input.driver             input_crs = input.crs             input_schema = { 'geometry': 'Polygon', 'properties': { 'INT_DATE': 'int' } }              #Open the file to be written             with collection(path_out, "w", driver=input_driver, crs=input_crs, schema=input_schema) as output:                 shapes = []                      #Append all of the polygons to a list                 for f in input:                     shapes.append(shape(f['geometry']))                      #Merge all the polygons                 merged = MultiPolygon(cascaded_union(shapes))                      #Write the new shapefile                 output.write({'geometry': mapping(merged), 'properties':{'INT_DATE':'4'}})          if __name__ == "__main__":         import sys         merge(sys.argv[1],sys.argv[2])      When I run this program it returns no error message and outputs a shapefile with the attribute table I expect. However, the polygon is not there. Does anyone have any idea what is causing this? I have relied heavily on the Fiona users manual and "GIS with Python, Shapely, and Fiona" by Tom MacWright. Is there a better way to use Shapely than using Fiona?


title: View for SDE Permissions for SQLServer

text: Does anyone know the correct tables and joins to write a view to see the permissions on objects in SDE? Here is what I had in 9.3 and it doesn't seem to work (ie., all Permissions are Edit b/c the actadd are all NULL):               SELECT dbo.sysobjects.name, dbo.sysusers.name      AS [User/Role], CASE [actadd]      WHEN 1 THEN 'View' ELSE 'Edit'      END AS Permission, dbo.sysobjects.xtype, dbo.sysobjects.status, dbo.sysobjects.category, dbo.sysobjects.sysstat, dbo.syspermissions.actadd          FROM dbo.sysusers      INNER JOIN dbo.syspermissions ON dbo.sysusers.uid = dbo.syspermissions.grantee      RIGHT OUTER JOIN dbo.sysobjects ON dbo.syspermissions.id = dbo.sysobjects.id      Likewise I look in syspermission table and all fields are NULL besides id, grantee and grantor.


title: How "complex" should a JavaScript map become? When does a plugin make more sense?

text: Background: A client is looking to upgrade an ancient ArcIMS site, which is quite complex. It has a table of contents, the ability to change symbology on-the-fly, print the map, use spatial bookmarks, identify the "active" layer, export data, etc. It's almost like an online version of ArcView 3. They wish to upgrade to ArcGIS Server 10, and are seeking recommendations on which of the ArcGIS Server APIs to use. I'm most familiar with the ArcGIS Server JavaScript API, so I don't know whether the Flex or Silverlight versions would offer advantages over the JS API. 1) Is there a point of complexity after which the JS API no longer makes sense, and a plugin would offer easier development and better performance? 2) Is it easier to develop a complex application (as described above) using Flex or Silverlight? 3) Are there any online examples of very complex ("ArcView-like") sites built on the JavaScript API? Thanks, Steve


title: How to move point in Spatialite Layer through manual entry of X/Y coordinates?

text: My Spatialite Database (Point Layer) looks like this (arch_locations.sqlite)               Fields:     pkuid     short_ID     name     geometry     X     Y      The field "X" and"Y" I added later (after first creation of DB) using the editing functions within the attribute table. They are of type "decimal number double" and I assigned the expressions $x and $y. Now they show X/Y Coordinates for viewing in an Access Database I have linked to the Spatialite Database through ODBC. Is it possible for me to move a Point on the layer by manually changing the values in the X and Y column? I will have exact coordinates for some points while others will derive their location through digitalizing georeferenced scanned maps. Devdatta Tengshe hinted that this might be possible in my other question. Thanks for reading! Edit to clear things up: I created the Database in QGIS through "Layer - new - spatialite Layer". I set it to Point and geometry and digitalized some points on a georeferenced map. Then I wanted to have two columns which show me X/Y Coordinates to be viewed in Access. I asked here how to do it: How to add columns with x/y data to Spatialite Layer? Following this advice I went to the attribute table of the layer, enabled editing and added two columns. Then I used the field calculator to put in the two new columns the expressions $x and $y. Now I only want to move some of the points I digitalized to an exact position by entering the x/y coordinates in those two columns. It doesn´t matter if I have to add the coordinated somewhere else I just want some points to be moved slightly.


title: Batch Clip in ModelBuilder of ArcGIS for Desktop?

text: I have multiple rasters and feature classes that I want to clip to a study area. I am pretty new to ModelBuilder. I know you can right click the clip tool to batch clip outside of ModelBuilder but how do I bring that into ModelBuilder at ArcGIS 10.1 for Desktop?


title: Layer overlay problem in Openlayers /Geoserver

text:
Hallo I am trying to create a Web Map servise using Openlayers and Geoserver. I want Open Street Map (OSM) as back ground and that both vagar_alla and adresspt layers to be tranparent, so that it will be possible both layers at the same time. I have the following code sofar. The problem is I OSM layer appears only at the begining zoom-in it disappears AND the top most layer obscures the other layers beneath. Any help?               <?xml version="1.0" encoding="UTF-8"?>     <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">                                                                                Lund Map

location

                     --> Click on the map to get feature info


title: how to add HTML popup in exported PDF map?

text: I am producing some maps of a point layer using arcpy.mapping.ExportToPDF. Is there a possibility of creating HTML popup for the final PDF based on this link? http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Using_HTML_pop_ups_for_feature_layers/00s500000018000000/


title: How to write to shape file in MapGuide

text: Here the code I have written to a Shape file in MapGuide Open Source (in PHP running on IIS 7.5). The code is actually the modified class MarkupEditor (markupeditor.php) from the General Tracks (Gt) Example.               class MarkupEditor     {         private $args = null;         private $site = null;              function __construct($args){             $this->args = $args;             $this->site = new MgSiteConnection();             $this->site->Open(new MgUserInformation($args['SESSION']));         }         //         // Rest of the original code in the class         //         function AddPointToShp($featureID){             $geometryFactory = new MgGeometryFactory();             $agfWriter = new MgAgfReaderWriter();             $vertices = explode(",", $this->args['GEOMETRY']);             $coord = $geometryFactory->CreateCoordinateXY((double) $vertices[0], (double) $vertices[1]);             $point = $geometryFactory->CreatePoint($coord);             $point = $point->Transform($this->GetTransform());             $byteReader = $agfWriter->Write($point);             $propertyValues = new MgPropertyCollection();             $propertyValues->Add(new MgInt32Property("FeatId", $featureID));             $propertyValues->Add(new MgStringProperty("TYPE", "Breader"));             $propertyValues->Add(new MgStringProperty("DESC", "OCB'));             $propertyValues->Add(new MgGeometryProperty("Geometry", $byteReader));             $this->InsertMarkupFeatureToShp($propertyValues);         }              function InsertMarkupFeatureToShp($propertyValues){             try{                 $featureService = $this->site->CreateService(MgServiceType::FeatureService);                 $featureSourceId = $this->args['MARKUPLAYER'];                 $commands = new MgFeatureCommandCollection();                 $commands->Add(new MgInsertFeatures($this->args['OPENMAKUP'], $propertyValues));                 $featureService->UpdateFeatures($featureSourceId, $commands, false);             }             catch (MgFeatureServiceException $ex){                 echo '<b>MgFeatureServiceException : </b>'.$ex->GetExceptionMessage();             }             catch (MgInvalidArgumentException $ex){                 echo '<b>MgInvalidArgumentException : </b>'.$ex->GetExceptionMessage();             }             catch (MgInvalidOperationException $ex){                 echo '<b>MgInvalidOperationException : </b>'.$ex->GetExceptionMessage();             }             catch (MgFdoException $ex){                 echo '<b>MgFdoException : </b>'.$ex->GetExceptionMessage();             }         }     }      The members of $args array passed in the constructor are:               Array     (         [MARKUPLAYER] => Library://DataStore/Data/Electrical/Breaker.FeatureSource         [SESSION] => 970bab06-4bd2-11e1-8000-000000000000_en_7F0000010AFC0AFB0AFA         [MAPNAME] => DataStore         [OPENMARKUP] => Default:breaker         [TGT] => 1         [POPUP] => 0         [GEOMETRY] => 1.5041909054501184,0.39827301781943014     )      **`$featureID`** in function **`AddPointToShp()`** is a numeric value. The code always stops execution at this line in function **InsertMarkupFeatureToShp()** :               $featureService->UpdateFeatures($featureSourceId, $commands, false);      It does not even generate an exception so that I may be able to find out what is wrong. Some posts on **osgeo.org** say that GT Example is used to write to SDF only. Where as, I am able to read any field from a desired SHP file using some other modified functions in the same file. Any help please? **UPDATE** I have converted the target Coordinate System to the one as in existing SHP file. Now I get the **MgInvalidArgumentException** :               Invalid argument(s): [1] = " " The string cannot be empty.      Also I have tested **`$propertyValues`** [`$propertyValues->GetCount()`]. It returns the count as 4.


title: Comparing polygon attributes that share the same parcel?

text: I'm new to ArcGIS, or at least anything beyond the basics, especially ArcToolBox. so the answer might be extremely simple. Basically I have one layer of building footprints, and another layer of parcel data. Many residential parcels contain multiple buildings, one of which is the house (coded in field BTYPE as 401) and one or more accessory buildings (BTYPE 999). each building footprint has a height attribute. What I'd like to do is highlight parcels that contain a 401 footprint and a 999 footprint, but only the parcels where the height of the 999 footprint is higher than the height of the 401 footprint. I have no idea how to do this.. I have tried a spatial join between the two layers but all that gets me is a count of the footprints on the parcel. I don't need the exact answer but a rough outline of steps that I need to take to accomplish this would be welcome.


title: Google Map Maker: How do I prevent users from re-adding streets from outdate satelite imagery

text: I live in a street that used to have a row of parking spaces in the middle. In Google Maps, this was represented as two one way lanes. Now, all traffic has been merged onto one of the lanes because they the building on the other side was torn down and a larger building is currently under construction. This can not be seen on the satellite imagery yet, but I was still able to correct the map pretty accurately - and this took me quite a while with all these approval roundrips. Today, I saw that a Google Maps Editor just re-added the street, undoing all my work I have done in that part of the map. This is really frustrating. Any ideas to prevent this?


title: Does merging segments in Google Map Maker ever work?

text: When editing line segments in Google Map Maker, the number of crossings one needs to touch for a change is pretty relevant. From my experience, every move of a crossing needs to be reviewed, so the number of crossings needed to be moved is the number of round trips one needs for a change. So, one way to improve this would be to remove crossings which anyway only connect two segments first. However, I've never been able to do this. I always get an error message saying that the segments to be merged have different attributes. I tried to find the differences, but without success. But let's start with a general question: Did anyone ever successfully remove a crossing in Google Map Maker? * * * This is the full error message (in German): > Segmente mit ungleichen Attributen können nicht zusammengeführt und die > entsprechenden Kreuzungen können nicht gelöscht werden. Die Segmente, die an > dieser Kreuzung enden, weisen Attribute auf, die nicht zueinander passen. > Bitte korrigieren Sie dies und versuchen Sie dann erneut, die Kreuzung zu > löschen.


title: How to compare two spatial point patterns?

text: If I have two point pattern distributions within the same geographic region, how would I go about visually and quantitatively comparing those two distributions? Also assume I have many points within a smaller region, so simply displaying a pin map is uninformative. Any examples or resources on calculations are always appreciated.


title: Better way to get raster values for a large number of point features?

text: I have several projects with tree point features (SpecialStatusSpeciesPoint) that contain up to 2,000,000 records. I am trying to record tree height (Hgt) and base elevation (BaseElev) values, based on DEMs representing vegetation height (veg_HGT) and bare earth (BEm). I am using ArcMap 10.1 with Spatial Analyst. Features are in a File GDB. DEMs are initially IMGs, but moved into a scratch file gdb as a GRID. I initially tried Sample but was left with inconsistently null values in the attribute table. The same was true for Extract Multi Values to Points as well. Ultimately, I found Zonal Statistics as Table to work, but at a painstaking slow speed. Sometimes it would crash or not complete all records. I thought if I broke the points up into chunks, it would ensure all records were copied and perhaps in less time. I wrote the following code (with the aid of snippets from you guys) as part of a larger script. All necessary modules and variables are taken care of in the code. For the sake of brevity, I only included the following. The processing time for only 10,000 records is on average 55 mins, 83 mins for 100,000 records, and 88 mins for 250,000 records (the code below).               # Define workspaces     arcpy.env.workspace = sdsfie_db     arcpy.env.scratchWorkspace = scratchGDB     arcpy.env.overwriteOutput = True     print("Adding attributes to SpecialStatusSpeciesPoint...")          # Add Field Id     if (not FieldExist(SpecialStatusSpeciesPoint, "Id")):         arcpy.AddField_management(SpecialStatusSpeciesPoint, "Id", "LONG", "", "", "", "Id", "NULLABLE", "NON_REQUIRED", "")          # Use updateCursor to populate Id field     print("      Copying Id field...")     rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint)     for row in rows:         OID = row.getValue("OBJECTID")         row.setValue("Id", OID)         rows.updateRow(row)     del row, rows          # Count number of records     recCount = int(arcpy.GetCount_management(SpecialStatusSpeciesPoint).getOutput(0))          # Decide whether to break data up into chunks     if recCount > 250000:         # Break data into chunks         print("   Breaking data into chunks of 250,000...")         inDesc = arcpy.Describe(SpecialStatusSpeciesPoint)         oidName = arcpy.AddFieldDelimiters(SpecialStatusSpeciesPoint, inDesc.oidFieldName)         sql = '%s = (select min(%s) from %s)' % (oidName,oidName,os.path.basename(SpecialStatusSpeciesPoint))         cur = arcpy.da.SearchCursor(SpecialStatusSpeciesPoint,[inDesc.oidFieldName],sql)         minOID = cur.next()[0]         del cur, sql         sql = '%s = (select max(%s) from %s)' % (oidName,oidName,os.path.basename(SpecialStatusSpeciesPoint))         cur = arcpy.da.SearchCursor(SpecialStatusSpeciesPoint,[inDesc.oidFieldName],sql)         maxOID = cur.next()[0]         del cur, sql         breaks = range(minOID,maxOID)[0:-1:250000] # makes slices of 250000         breaks.append(maxOID+1)         ranges = []         for b in range(len(breaks)-1):             ranges.append([breaks[b], breaks[b+1]-1])         numChunks = str(len(ranges))         print("      Broken into " + numChunks + " chunks...")              # Begin timer         timerStart = time.time()              tempOut = scratchGDB + "\\tempOut"         counter = 0              # Iterate through chunks         for i,j in ranges:             counter += 1             lapStart = time.time()             print("   Processing chunk " + str(counter) + " of " + numChunks + ". ID " + str(i) + " - " + str(j) + "...")                  # Build sql where_clause independent of GDB type             fieldDelimited = arcpy.AddFieldDelimiters(SpecialStatusSpeciesPoint, "OBJECTID")             fieldType = arcpy.ListFields(SpecialStatusSpeciesPoint, "OBJECTID")[0].type             if str(fieldType) == 'String':                 i = "'%s'" % i                 j = "'%s'" % j             whereClause = "%s BETWEEN %s AND %s" % (fieldDelimited, i, j)                  arcpy.env.workspace = sdsfie_db             arcpy.env.scratchWorkspace = scratchGDB             arcpy.env.overwriteOutput = True                  # Create layer of chunk             print("      Creating temp layer...")             # Create temp layer slice             arcpy.Select_analysis(SpecialStatusSpeciesPoint, tempOut, whereClause)             # Count features in temp layer             tempOutCount = arcpy.GetCount_management(tempOut).getOutput(0)             # Zonal Statistics as Table             print("      Creating zoneTable7...")             arcpy.gp.ZonalStatisticsAsTable_sa(tempOut, "Id", veg_HGT, zoneTable7, "DATA", "MAXIMUM")             # Join Field             print("      Joining MAX field...")             arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable7, "Id", "MAX")             # Zonal Statistics as Table             print("      Creating zoneTable8...")             arcpy.gp.ZonalStatisticsAsTable_sa(tempOut, "Id", BEm, zoneTable8, "DATA", "MINIMUM")             # Join Field             print("      Joining MIN field...")             arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable8, "Id", "MIN")             # Use updateCursor instead of calculate field             print("      Copying MAX and MIN fields...")             rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint, whereClause)             for row in rows:                 veg_HGT_max = row.getValue("MAX")                 BEm_min = row.getValue("MIN")                 row.setValue("Hgt", veg_HGT_max)                 row.setValue("BaseElev", BEm_min)                 rows.updateRow(row)             del row, rows             # Clean up for next iteration             arcpy.DeleteField_management(SpecialStatusSpeciesPoint, "MAX;MIN")             arcpy.Delete_management("in_memory")             # Print lap time             lapStop = time.time()             lapTime = lapStop - lapStart             print("   Chunk run time: " + str(lapTime/60) + " minutes...")     else:         rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint)         for row in rows:             OID = row.getValue("OBJECTID")             row.setValue("Id", OID)             rows.updateRow(row)         del row, rows         # Process: Zonal Statistics as Table         arcpy.gp.ZonalStatisticsAsTable_sa(SpecialStatusSpeciesPoint, "Id", veg_HGT, zoneTable7, "DATA", "MAXIMUM")         # Process: Join Field         arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable7, "Id", "MAX")         # Process: Zonal Statistics as Table         arcpy.gp.ZonalStatisticsAsTable_sa(SpecialStatusSpeciesPoint, "Id", BEm, zoneTable8, "DATA", "MINIMUM")         # Process: Join Field         arcpy.JoinField_management(SpecialStatusSpeciesPoint, "Id", zoneTable8, "Id", "MIN")         rows = arcpy.UpdateCursor(SpecialStatusSpeciesPoint)         for row in rows:             veg_HGT_max = row.getValue("MAX")             BEm_min = row.getValue("MIN")             row.setValue("Hgt", veg_HGT_max)             row.setValue("BaseElev", BEm_min)             rows.updateRow(row)         del row, rows          # Clean up     arcpy.DeleteField_management(SpecialStatusSpeciesPoint, "Id;MAX;MIN")     arcpy.Delete_management("in_memory")     arcpy.Delete_management(scratchGDB)          # Stop timer and print results     timerStop = time.time()     elapsed = timerStop - timerStart     hrs, remain = divmod(elapsed, 3600)     mins, secs = divmod(remain, 60)     print("   SpecialStatusSpeciesPoint processing time: %d:%d:%f" % (hrs, mins, secs))      What might be the bottleneck? Is there something terribly wrong with the code? (I began learning Python 4 months ago) Has anyone who has had similar issues with Sample/Extract Multi Values to Points/zonal Statistics as Table found a more efficient workaround for null values or missing records? Thank you in advance for any help you can offer.


title: Is it possible to create a tooltip with a map service?

text: I am trying to modify the ESRI Parks Finder application to include a tooltip that shows the name of a park when you mouse hover over the icon. It is a map service. Most of the documentation I have seen shows doing this with either a feature service or graphic. In fact, there are other areas in the application that utilize tooltips but they are based on graphics not the actual symbology of the data. Any ideas or help on how I can go about this?


title: MSHP file from SAGA tough to delete

text: After running a SAGA line dissolve process in QGIS, an MSHP file is created (amongst other things). Fine by me, though I'm not sure what it's used for (metadata in XML format?). Problem is that it can't be deleted without closing QGIS. This is a problem, as I'm running the process repeatedly in a Python script. Here is the relevant code:               lines_layer = some_vector_polygon_layer     dissolved_lines_path = "SOME SHAPEFILE PATH HERE"     unused_field = lines_layer.dataProvider().fields()[0].name() # Name of a field, not used     processing.runalg("saga:linedissolve", lines_layer, unused_field, unused_field, unused_field, 1, dissolved_lines_path)      I then load the resulting layer and add it to the map registry & table of contents. Shortly after, I delete the layer from the map registry (& table of contents). I can delete the resulting SHP quite fine, but I am unable to delete the MSHP. There are no errors in the resulting shapefile or in the processing log. Any ideas?


title: getting a merged color satellite image from red,green, blue, pan, and nir tifs

text: I am using Tiff images from Digital Globe and trying to get a color satellite image. I have merged the red, green, blue, pan and nir tiff files into one file. I then changed the band properties to a max of 500 and contrasted enhancement to stretch and clip to the min and max values. Doing all this produced a merged satellite image but the image has no color and is being rendered as single band gray. How can I create a colored image? I'm assuming I need to somehow create a tif such that each red, blue, green, nir band is a band for the tiff?


title: Dynamically generating data for map tiles from Python function

text: I'm brand-new to GIS and relatively new to web services. I've found a ton of great resources on the web, but I'm having trouble addressing one specific issue: One of the data layers on my planned map is dynamically generated in Python from some scattered and dynamic data (it is the "results" output of a computational weather-related model). I can generate this layer in chunks of arbitrary size and location, and output in any standard format (a numpy array, a .png, a geotiff (using GDAL), etc.), and I can expose a function such as:               def create_tile(width, height, ulc_lat, ulc_lon, lrc_lat, lrc_lon):      to hide the implementation details. But, how do I get this layer to display in the map? The options that I can see right now are:   * Serve as a WMS by creating a new Layer Type in a Python-based WMS server such as TileCache or TileLite ... this is probably what I'll try unless someone has a better idea.   * Implement as a PythonDataSource in Mapnik. _Problem:_ It's a new feature, requires custom compilation, and causes Mapnik to segfault even on the provided examples.   * Somehow map filenames on disk (my_tile_z0_x3_y4) so that their contents is generated on-demand by the create_tile function above... I have no idea how this would be done, though.   * Create the layer on disk ahead of time, and re-generate it often to keep the appearance of being "dynamic". (lots of disk space, needless CPU use, etc) Does anybody have any suggestions? It seems like this must be a common problem, but I can't find any standard solution... Thanks! Rob


title: Aggregating a raster to fit a polygon grid

text: I am working at the scale of 10 square km grid cells for the UK, and have a polygon layer of the Ordnance Survey 10k grid cells in ArcGIS. I have another layer of the human footprint index, also for the entire UK, that I am using to stratify my data. However, the resolution of this human footprint index is much greater than that of the 10k grid cells, so I want to aggregate the cells, such that I have a human footprint index value for each 10 sqkm cell. When I try to do this, however, the squares and grid do not line up - as can be seen in the pictures. They are in the same projection. Is there any way to aggregate a layer to the confines of a polygon grid? Any help would be very much appreciated! Images here: ![Before aggregation](http://i.stack.imgur.com/KOH9C.png)   ![After aggregation](http://i.stack.imgur.com/WsC15.png)


title: How to get SAGA working with QGIS Sextante?

text: I just created my account on this forum to find an answer to a stupid (?) problem I can't deal with.   First, let's precise I just (some months ago) learned on my own what GIS was, and gave it a try using QGIS. I'm working on Windows XP. I apologize in advance for my not so good English, and my poor skills (and maybe bad vocabulary) in GIS and kriging. My problem is: I need to use kriging to interpolate punctual data from a vector layer (points); this layer is created by me from fieldwork measurements and is in .shp format (try.shp).   I had a look on the internet and saw there were 2 possibilities : use SDA4PP plugin, or Sextante toolbox with SAGA. I wasn't able to have any of them working.   Apparently, SAGA is much easier, thus I would prefer to work with SAGA. What I have done so far: I downloaded the last version of QGIS via OSGeo4W (all packages selected).   I downloaded SAGA 2.0.8 from Sourceforge.net (version *_bin_msw_win32), activated it Sextante toolbox, entered the right folder path (which contains spaces) in Sextante configuration. When I want to perform ordinary kriging, I select the layer, the attribute (which is "id"), ask for a variance grid, target grid is "[0] user defined" (no other choice), I choose an exponentiel variogram model, no block kriging and no logarithmique transformation, and let the default values for the other parameters (especially, grid and variance rasters are saved to temporary files). Then I click "OK" and get this message error:               Could not layer: C:\DOCUME~1\user\LOCALS~1\Temp\sextante\sagaordinarykriging6fd1782d35df....tif. Check the Sextante log to look for errors.      Here is the error log :               Could not load model watersheds.model     ALGORITHM:ftools:export/addgeometrycolumns      and the algorithm log :               sextante.runalg("saga:ordinarykriging","C:/Documents and Settings/user/Bureau/Kriging/try.shp","id",True,0,1,False,100,False,0,10,100,1,0.1,1,0.5,1000,4,20,1,True,"-55.4887564869,-55.4887477484,-8.61474380806e-07,4.59453003047e-06",None,None)      Any idea of the reason why it doesn't work?   Thanks. * * * I tried to uninstall-reinstall OSGeo4W using the advanced options. I selected SAGA as a package when reinstalling, but SAGA was not really installed: the tar.gz was downloaded and that's all. I unzipped it, and then I had to enter manually the path, which I chose to be without spaces (C:\OSGeo4W\saga, after having cut-pasted the SAGA unzipped folder at this very place). And ... same as before ... Same error, same behaviour, all the same. Any clue ? -- Edit: The installation worked well, no errors. I precise I removed all the files concerning SAGA and QGIS before reinstalling (with Add/Remove Programs and Search). * * * It worked ! I really can't explain why, but it worked. Until now, I was working on a self-made point shapefile, with random points with "id" within the range 0-150. Recently I used a file containing real data from fieldwork, and it worked. I did exactly as for the first file, no clue why it suddenly worked. One more thing. When using SAGA software, I can choose to select a variogram or not; it seems I can't select a variogram using SAGA via the Sextante toolbox, even if I want to. Why ? Is there a way to select a variogram using SAGA via the Sextante toolbox ? Besides, I guess it's always better to select a variogram (I'm a beginner in kriging). What do you think ? Cheers


title: How to compute area for polygons in an EPSG:4326 shapefile?

text: I have a shapefile covering the whole globe, in EPSG:4326 (geodetic WGS84). I would need to compute the polygons areas in a metric such as km2. How can I do that in QGIS or with gdal/ogr and having a good precision?


title: Stop rounding values in Qgis attribute table

text: I am calculating the areas of my polygons in Qgis using $area and everything works fine until I click on "save". All my values are rounded to an integer. I tired setting the output field type to decimal but still see the same problem? Any solution to that? I am using Qgis 1.8


title: What does `findCloseEdges` do in osm2po?

text: I have something like this:               int[] edges = graph.findCloseEdges(lat, lon, 1); // returns only one element     RoutingResultSegment rrs = graph.lookupSegment(edges[0]);     LatLon[] latlons = rrs.getLatLons()      For the red dot in the image (with `lat`/`lon`), array `latlons` contains the blue points. I would instead expect "close edge" to be a segment of the A3/E70 road. Here is the image (I don't have enough reputation to post images apparently). What is the third integer argument in `findCloseGraph` and what does this method return? I need to do simple "snap to road".


title: GeoJSON in Leaflet

text: I'm just starting to work with leaflet and have a basic question about showing a geoJSON file. First, I found this excellent article with example code: http://bl.ocks.org/lxbarth/4019660 However I have a more basic/beginner question, in that example the set the JSON file through this url: http://polymaps.appspot.com/county/{z}/{x}/{y}.json However, the file I converted (from shapefile to geojson using quantum gis) is just a .geojson extension. I'm clearly missing something here, perhaps a way to convert the geojson file into a tile format? Can someone elaborate as to how I'm supposed to do this or what my missing step is? Do I need to be hosting the geoJSON file on a server in order to do this? If my file is called polygons.geojson, how would I go about using it in the above reference code? Thanks!


title: Google base layer opacity change in OpenLayers 2.13

text: I have recently swapped from using OpenLayers 2.12 to 2.13.1 as it fixes a bug. However in 2.12 I was using a base layer opacity of 0.7 to 'knock back' the map so that the vector layers stood out more, but this doesn't work in 2.13. In 2.13 OpenLayers the base layer opacity on a Google map is applied to all layers, so this trick doesn't work. I believe this is because 2.13 changes how Google base layers are handled (see release notes). Anyone have a suggested work round that will work in OpenLayers 2.13.1? i.e. I want to make the map layer more transparent/less contrast while keeping the vector layers at full contrast so they stand out. Also I am testing with Google but may go with a different map provider in the final application, due to licencing costs, so a generic answer for any base map type would be helpful. Thanks in advance for your help.


title: How do I find location of minimum value within a buffer (around a point)

text: 1. I have a digital elevation model (DEM) with a point location defined.    2. I want to locate the lowest height values at a set distance from this point.    3. I envisaged creating a buffer around the point and then locating the lowest height values; firstly is this the best way of achieving this, if not can you recommend another option ?    4. If it is the best way, once the buffer is created around the point, how do I extract DEM cells with the lowest height values? I dont want to know what their values are, I want to be able to identify the cells and know where the cells are.  Thanks for any assistance you can provide.


title: Where can I find GIS data to test GIS coordinate operations?

text: I am trying to find some data that can be used to test GIS coordinate operations from one CRS to another CRS. So far all I have found is the Gold Data set and also various example points in the EPSG G7-2. The gold data is OK for testing some projections from WGE and the few points in the G7-2 document are good for testing individual operations. What I need now though are some points that can test implicit as well as explicit concatenated coordidnate operations such as EPSG:4277 to EPSG:4230 for example, which I had found a sample for in http://www.epsg.org/Exchange/1065.pdf. There has to be a good source out there with more than one point and CRSs using multiple datums. Fly, Google monkeys, fly! **EDIT:** I also found something called GIGS


title: How to create an oriented buffer using arcpy?

text: I would like to create an oriented buffer for every polygon in my shapefile using arcpy. By oriented I mean that I have two angles a1 and a2 that constraint the direction of the buffer. This is represented in the graph below: ![enter image description here](http://i.stack.imgur.com/tjCrb.jpg) Any ideas?


title: How to create heatmap (concentration of points) in QGIS?

text: There was already a question about it in 2010, where the answer was basically "use QGIS Heatmap plugin". I've also read on QGIS webpage that it is now "implemented as a native C++ plugin" My question is probably stupid, but I haven't found the answer: what is the best way to generate a heatmap from a list of observations, in QGIS?


title: How to determine the longest line segment within a polygon that passes through its centroid?

text: A PhD student approached me recently asking how one would determine the longest straight line that passes through the centroid of a polygon, the output being a polyline rather than just a table of numbers. The polyline would be within the polygon starting at the edge at one end, passing through the centroid, and then end at the opposite edge. Does anyone know how to calculate this? I'm surprised (but I guess I'm showing my ignorance) that this is not some interface in ArcObjects. I would have posted an image to explain what I am asking (I hope it is clear) but as a newbie to this site I was not allowed to upload images! The polygon represents the crater edge of a volcano with no islands so the polygon can have an irregular shape.


title: Problems with spatial autocorrelation in Geographically Weighted Regression (GWR) analysis

text: I have conducted an **Ordinary Least Squares** (OLS) and **Geographically Weighted Regression** (GWR) analysis in ArcGIS where I have tried to predict biodiversity patterns in an area. When I examine the results with **Spatial Autocorrelation** the residuals tends to be clustered with a fairly high Z-value. The result from the GWR, however, tends to be the expected when I visualize it. How do I handle the clustering? What problems does it creates, IS there any problems?


title: Make selected polygon a different color without covering everything

text: Working in ArcGIS 9.3 I'm selecting polygons by attributes and trying to have the selection show up as a light yellow. I set the selection color to the light yellow and it's the 'lowest/last' layer but it covers all of the other layers. Is there a way to change it? Thanks!


title: Why does only one of the following queries throw an error?

text: There are 2 SQL queries:               SELECT wkb_geometry FROM building_polygon WHERE ST_Contains(wkb_geometry, ST_GeomFromText('POINT (29.4 60.1)'));     SELECT wkb_geometry FROM building_polygon WHERE ST_Contains(wkb_geometry, ST_GeomFromText('POINT (30.32670 59.92843)'));      The first query works fine but returns 0 rows. But the second query throws an error `Operation on mixed SRID geometries`. It is not problem to point SRID for second query but it is not clear for me why the first query doesn't returns the same error. Is it expected behavior?


title: What's the purpose of wms_metadata.xml and admin.sld?

text: I've installed QGIS Server and it works fine. But in /usr/lib/cgi-bin/ folder I've found 2 additional files - `wms_metadata.xml` and `admin.sld`. What is the purpose of these file? I've never used them.


title: How to perform repetitive tasks in QGIS?

text: I am trying to do a treatment on many feature files, so I would like to automate it. In fact, I have one shapefile with the spatial distribution of some species, and one with vegetation type. I would like to select (by attribute) one species in the Species Shapefile, then to select (by locality) all vegetation areas which intersect with their distribution area. Finally, I would like to have a shapefile with, as name, the species name, and with attributes and forms of the frequented vegetation types. And I would like to repeat this for all species (more than 100), and if possible do this in an easy manner (so it could be done by another person). I have already tried this task using Sextante plugin, but I can't have the species name as a shapefile name at the end. Can someone suggest a method for this?


title: How to declare folder name as global?

text: I have asked similar question at stackoverflow but posting here as well since it involves ArcGIS and python Add-In. Following code shows first step where i have to declare output folder as global so that later outputs can be saved in it as well. Right now I am getting an error at output folder string r'optfile/ras1'. Any help how to correctly store files in output folder and declare it as global would be appreciative.               import arcpy     import os     import pythonaddins          from datetime import datetime          now = datetime.now()     month = now.month     year = now.year          optfile = "C:/temp/"+str(year)+"_"+str(month)          class DrawRectangle(object):     """Implementation for rectangle_addin.tool (Tool)"""         def __init__(self):             self.enabled = True             self.cursor = 1             self.shape = 'Rectangle'             os.makedirs(optfile)                  def onRectangle(self, rectangle_geometry):         """Occurs when the rectangle is drawn and the mouse button is released.         The rectangle is a extent object."""              extent = rectangle_geometry         arcpy.Clip_management(r'D:/test',                                "%f %f %f %f" % (extent.XMin, extent.YMin, extent.XMax, extent.YMax),                                r'optfile/ras1', "#", "#", "NONE")         arcpy.RefreshActiveView()


title: How to calculate orientation of line segments using open source GIS?

text: I'd like to calculate the orientation of line segments relative to north direction using open source tools. Which tools or functions would you recommend?


title: How to calculate polygon length and width using open source GIS?

text: How can I calculate the length and with of polygons using open source tools? By "polygon length" I mean the length of the longest line within the polygon (update: to be correct, I just need the longest line within the polygon's convex hull) and the width as the longest measure perpendicular to the length measurement. **Update:** I'm trying to create a non-axis aligned minimum bounding rectangle for a polygon. The length of this minimum bounding rectangle would be the polygon's diameter.


title: How to write a script that will automatically draw lines from one point to multiple destination points? (ArcMap 9.3)

text: I need to create lines from an airport to numerous destination locations. The destinations change often, so drawing them manually becomes tedious. I want to create a script that will draw them automatically. Does anyone know how to do this? I've written a number of scripts for ArcGIS 10 using arcpy, but nothing for 9.3. If you could point me in the right direction that would be great.


title: How to improve speed of rendering polygons from Postgis in Leaflet?

text: I want to improve speed of rendering polygons in my Leaflet map. Polygons are stored in Postgis, and are rendered with select query with bbox. I wonder can I cache somehow those polygons? I read about Tilestache a bit, but I'm still not sure if that would help in speed performance. Also, one of my problems is that I need to relate postgis tables with non- spatial tables in my postgresql. So, some of the attributes I need to show in the popup are in some other non-spatial tables in Postgresql. Just to be more clear, that map is already developed and all works but I want to improve speed of polygons rendering. Can you give me some advices?


title: How to transform shp file into raster file in Quantum GIS?

text: I'd like to transform a shape file into raster file in Quantum GIS... and the opposite way back also (raster to shp). Can you help me?


title: MapServer MapCache WMS

text: I just installed mapserver and mapcache on a Debian server. I struggled to put those two to work :D. Now when I call this url (by the way I don't know if it's the good one) //myserver/cgi- bin/mapserv/service/wms?map=/opt/basemaps/osm- google.map&mode=browse&template=openlayers, I get the map pat it is just a single image ! no tiles. If the map appears this means my mapcache works fine ? Is there anything i missed.               ./mapserv -v     MapServer version 6.5-dev OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WFS_SERVER SUPPORTS=WCS_SERVER SUPPORTS=FASTCGI SUPPORTS=GEOS INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE      Thank you Reagards


title: unionize intersections of merged polygons

text: [Updated quesion] I am using the following code to merged polygons from shapefiles for the purpose of calculating the surface area;               library(maptools)     library(rgdal)     library(rgeos)          #enter slots, e.g. area     shp_list[[3]]@polygons[[1]]@area          setwd("~/type1/speciesX1/")     # Select species     species <- 'speciesX1'     print(species)          # List shapefile-dirs for that species      shp_dir <- list.files(path=species, pattern="*.shp", full.names = TRUE, recursive = TRUE,   include.dirs = FALSE)          # Read all the shapefiles and place them in a list     shp_list <- lapply(shp_dir, function(x) {readOGR(dsn=x, layer=ogrListLayers(x), verbose=T)})          # Check how many shapefiles have been loaded     number_of_loaded_shapefiles <- length(shp_list)     print(paste("Number of loaded shapefiles: ", number_of_loaded_shapefiles))          # Set the number of fields to 3, i.e. name, id, mrgid and set id and iso_3digit to id if present     set.seed(1234)     for(i in seq_along(shp_list)){     # EEZ       if('objectid' %in% colnames(shp_list[[i]]@data)==T){        shp_list[[i]]@data <- subset(shp_list[[i]]@data, select=-c(objectid, country, sovereign, remarks, sov_id, eez_id, iso_3digit, date_chang  ,area_m2, longitude ,latitude))        shp_list[[i]]@data['id']<-paste0(shp_list[[i]]@data$id,shp_list[[i]]@data$iso_3digit)        colnames(shp_list[[i]]@data)[colnames(shp_list[[i]]@data) == 'eez'] <- 'name'       }else{     # IHO + EEZ (e.g. Mexican part of the North Pacific Ocean, 25569)         if('iho_id' %in% colnames(shp_list[[i]]@data)==T){          shp_list[[i]]@data<-subset(shp_list[[i]]@data, select=-c(objectid_1, iho_sea, iho_id, iho_mrgid, country, sovereign, sov_id, eez, eez_mrgid, longitude, latitude, area_m2))          colnames(shp_list[[i]]@data)[colnames(shp_list[[i]]@data) == 'marregion'] <- 'name'          colnames(shp_list[[i]]@data)[colnames(shp_list[[i]]@data) == 'eez_id'] <- 'id'          shp_list[[i]]@data['id']<-paste0(shp_list[[i]]@data$id, paste0(sprintf("%1s", sample(LETTERS,1)),sprintf("%1s", sample(LETTERS,1))))          shp_list[[i]]@data <- shp_list[[i]]@data[,c(1,3,2)]         }else{     # MEOW (ecoregion, e.g. Bahamian 21980)           if('eco_code' %in% colnames(shp_list[[i]]@data)==T){            shp_list[[i]]@data<-subset(shp_list[[i]]@data, select=-c(lat, long, placetype))            colnames(shp_list[[i]]@data)[colnames(shp_list[[i]]@data) == 'ecoregion'] <- 'name'            colnames(shp_list[[i]]@data)[colnames(shp_list[[i]]@data) == 'eco_code'] <- 'id'            shp_list[[i]]@data['id']<-paste0(shp_list[[i]]@data$id, paste0(sprintf("%1s", sample(LETTERS,1)),sprintf("%1s", sample(LETTERS,1))))            shp_list[[i]]@data <- shp_list[[i]]@data[,c(2,1,3)]                }         }       }     }          # See the row names (IDs)     lapply(shp_list, function(x) row.names(x@data))          # Generate unique id string on the fly     set.seed(1234)     shp_list <- lapply(shp_list, function(x) {     new.id <- paste0(x@data$id, sprintf("%03d", round(runif(1, 1, 999)),0))     spChFIDs(x, new.id)     })          # Print mrgids     for(i in seq_along(shp_list)){       print(shp_list[[i]]@data$mrgid)     }          # Merge the SPDF      merged_shp_list <- do.call("rbind", shp_list)          # Some of the polygons (areas) are the same in the two original SPDFs.     # Check the attribute data     attribute_data_in_merged_shp_list <- merged_shp_list@data          # E.g. Alboran sea (mrgid=3324) is present in both. Make a spatial union based on a common ID-value (mrgid)     union_polygons <- unionSpatialPolygons(merged_shp_list, IDs=merged_shp_list@data$mrgid)          # Check how many polygons are present in the original merge        number_of_polygons_before_merge <- nrow(merged_shp_list)     print(paste('Number of polygons before merge: ', number_of_polygons_before_merge))          # How many polygons are present after merge     number_of_polygons_after_merge <- length(union_polygons)     print(paste('Number of polygons after merge: ', number_of_polygons_after_merge))          # unionSpatialPolygons returns a SpatialPolygons object, but we want to retain the original attribute data as well. Get only the unique rows from the original attribute data     attribute_data <- unique(merged_shp_list@data)          # Row names (IDs) need to match the IDs in union_polygons     rownames(attribute_data) <- attribute_data$mrgid          # Create a new SPDF     final_shp_list <- SpatialPolygonsDataFrame(Sr=union_polygons, data=attribute_data)          # Checks     number_of_merged_polygons_final <- nrow(final_shp_list)     print(paste('Number of merged polygons: ', number_of_merged_polygons_final))          # Calculate area     area <- gArea(SpatialPolygons(final_shp_list@polygons))     print(paste('Area (',species,')',':',area))          # Plot marged polygons     final_shp_plot <- spplot(final_shp_list, "mrgid", colorkey=F, fill='cyan')     print(final_shp_plot)      The shapefiles I'm using are from Marineregions.org In some cases, the merged polygons might overlap, depending on what type of shapefile have been used (e.g. EEZ, IHO, MEOW). I want to calculate the surface area of the merged polygons, but figured that as it currently is the calculated area includes potential intersections, so I wonder how to unionize these from a set of merged polygons. I've looked into the `rgeos package; gUnion, gIntersection, gDifference etc` but these functions only works on a pair. I'm new to this so please bare with me. Below are links to two shapefiles (choose shapefile in the Download drop- menu). For the code, place the files in `~/type1/speciesX1/` and name e.g. the folders 1 and 2 and the files within correspondingly. Sorry, for this part. **Example Data**   * Comoran part of the Mozambique Channel   * Mozambique Channel Comoran part of the Mozambique Channel is part of Mozambique Channel, thus we can see if the intersection is calculated twice for the merged polygon. **Output of Example Data** Plot of the merged polygons ![Plot of merged polygons](http://i.stack.imgur.com/Z3mNf.jpg) **Test**               area      125.6221 (as its calculated from above code)          # Area of the two polygons     A1 <- gArea(SpatialPolygons(shp_list[[1]]@polygons))     A2 <- gArea(SpatialPolygons(shp_list[[2]]@polygons))          A1     6.807844     A2     118.8142          # Area of the union     u <- gUnion(shp_list[[1]], shp_list[[2]])     uA <- gArea(SpatialPolygons(u@polygons))          uA     118.8142          # Area of the intersection     i <- gIntersection(shp_list[[1]],shp_list[[2]])     iA <- gArea(SpatialPolygons(i@polygons))          iA     6.807844 (same as A1)          uA + iA     125.6221 (same as calculated by the above code)      So the area calculated by the code does in fact include the intersection twice, i.e. no union is created, since the `unionSpatialPolygons` function creates unions based on `IDs`. Is there any function e.g. like `gUnion` (or how can I use gUnion?) for multiple (>2) spgeom objects, which I can use and add to my workflow after creating unions based on IDs? Any pointers would be highly appreciated. (If anything is unclear, please tell me so I can modify the question)


title: Symbolize layer and group names in TOC?

text: I regularly have quite large mxd projects, with lots of feature layers referenced in Table of Contents (TOC). For simplicity I group them together. The problem is that I can't see the difference between a feature layer or a group (or a group of groups). It is not only me using these documents, and it would be really useful to be able to colorize the name of groups and/or feature layers or have any other way of separate them from each other (other than having to press the +). At the very least I'd like to have one color for groups and another for feature layers, but being able to choose freely would be a big improvement. Then I can have blue for water related stuff, red for administrative, green for... you get the picture. Is there any such functionality in Arcmap? Or has anyone made a neat add-in that could do this? ![TOC](http://i.stack.imgur.com/zo19L.png) _Arcgis 10.1 on Windows 7_


title: Label halos in OpenLayers

text: I would like to add halos on my labels on an Openlayers map. I have tried adding the following to my code but no halos display.               labelHaloColor: "#FFFFFF",     labelHaloWidth: "2",      I got that from here: http://trac.osgeo.org/openlayers/ticket/2965 Do I need to use a newer version of OpenLayers to get this working? I am currently using the dev version.


title: Polygons for postcodes in Spain

text: Is there a source (free or paid) for polygons of border coordinates for spanish postcodes? Google show this borders in google maps, but doesn't provide them.


title: Toggle Layers from Multiple Services

text:
I'm trying to adapt my code so that I can create a layer list of all layers in multiple services. I'd like them to be able to be turned on/off. My current code builds a dynamic layer list from one service, however, I've recently been asked to add another service to the app and append the layers to the same div. I've found this question from a while back and am running into the same issue the OP described, where I can toggle one service off, ONCE, but never back on, and I can't toggle the second service at all. The code I have modified is below. Can anybody notice any glaring mistakes or things that I might be missing? I've attempted to contact the OP in the post I linked, but thought this might help get a fresh look.                     dojo.require("esri.map");                 var layer, parcellayer, map, visible = [];             var startExtent            function init() {                  //Set the initial map extent                   startExtent = new esri.geometry.Extent({               "xmin": -9199621.530456403,                   "ymin": 3365235.5214724857,                   "xmax": -9077475.159256855,                   "ymax": 3443965.6606061114,                   "spatialReference": {                   "wkid": 102100               }           });         map = new esri.Map("map", {               extent: startExtent,               logo: false           });              var basemap = new esri.layers.ArcGISTiledMapServiceLayer("http://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer");         map.addLayer(basemap);                   //Use the ImageParameters to set the visible layers in the map service during ArcGISDynamicMapServiceLayer construction.         var imageParameters = new esri.layers.ImageParameters();         imageParameters.layerIds = [0,1,2,3,4];         imageParameters.layerOption = esri.layers.ImageParameters.LAYER_OPTION_SHOW;         //can also be: LAYER_OPTION_EXCLUDE, LAYER_OPTION_HIDE, LAYER_OPTION_INCLUDE              layer = new esri.layers.ArcGISDynamicMapServiceLayer("http://gisprod2/ArcGIS/rest/services/Dynamic/Zoning/MapServer", {"imageParameters":imageParameters});         map.addLayer(layer);              //Use the ImageParameters to set the visible layers in the FEMA map service during ArcGISDynamicMapServiceLayer construction.         var FEMAimageParameters = new esri.layers.ImageParameters();         FEMAimageParameters.layerIds = [0];         FEMAimageParameters.layerOption = esri.layers.ImageParameters.LAYER_OPTION_SHOW;         //can also be: LAYER_OPTION_EXCLUDE, LAYER_OPTION_HIDE, LAYER_OPTION_INCLUDE              FEMAlayer = new esri.layers.ArcGISDynamicMapServiceLayer("http://gisprod2/arcgis/rest/services/Static/EnvironmentallySensitiveOverlayZones/MapServer", {"FEMAimageParameters":FEMAimageParameters});         map.addLayer(FEMAlayer);       }     //Sets up the ability for dynamic layers to be toggled.       function updateLayerVisibility() {         var inputs = dojo.query(".list_item"), input;         //in this application layer 2 is always on.         visible = [];         for (var i=0, il=inputs.length; i<il; i++) {           if (inputs[i].checked) {             visible.push(inputs[i].id);           }         }         //if there aren't any layers visible set the array value to = -1         if(visible.length === 0){         visible.push(-1);         }              layer.setVisibleLayers(visible);       }               //Sets up the ability for FEMA layers to be toggled.     function updateFEMALayerVisibility(){         var inputs = dojo.query(".FEMA_item"), input;         //Placing layer numbers below will enforce them to be visible.         visible = [];         for (var i = 0, il = inputs.length; i < il; i++) {             if (inputs[i].checked) {                 visible.push(inputs[i].id);             }         }         //if there aren't any layers visible set the array value to = -1         if(visible.length === 0){         visible.push(-1);         }              FEMAlayer.setVisibleLayers(visible);          }            dojo.addOnLoad(init);     </script>                          This sample loads an ArcGISDynamicMapServiceLayer and presents check boxes for only the layers that should be toggled on and off by users.
           FEMA:
         <input type='checkbox' class='FEMA_item' checked='checked' id='0' value='0' onclick='updateFEMALayerVisibility();' /> FEMA Layers

           Zoning:
         <input type='checkbox'  class='list_item' checked='checked' id='0' value=0 onclick='updateLayerVisibility();'/> Zoning Classifications


title: What software would be useful to create General Transit Feed Specification (GTFS)?

text: I am doing a project in which i am required to create GTFS. I look at transit google developer forum which provide a format to create feed but that completely have to done by hand (we have to write text files), i was looking for a software that can help in this matter.


title: OGIS Ole Db provider sample with read/write capability

text: I'm trying to determine what level of effort would be required to develop an OGIS Ole db provider for ArcGIS. There is a sample that is read-only. I would like to be able to edit featureclasses in the ole db workspace using the out-of-the-box editing tools in ArcMap. Plus, I would prefer to develop in C# instead of C++. I see that there was an effort to do this with zigGIS, but it never was completed. Does anyone know of sample code showing how to do this?


title: Is there a way to traverse and select all vectors connected to a starting vector?

text: Preferably using QGis I'd like to remove errors from a road network .shp file and was wondering if there was a away to perform a flood fill or traversal on a vector layer to select all connected vertexes from the starting vertex to do basic validation and remove unconnected roads. Thanks.


title: ArcGIS Calling Wrong fme.exe When Python Script Converted To An ArcGIS Script Tool

text: I'm using this this bit of code to execute a FME Workbench file through a Python script ...               import os     os.system("fme.exe P:\\Mapping\\Scripts\\FME\\Easements_ROW_Extract_From_CAD.fmw")      However, when I try using it as an ArcGIS script tool, ArcGIS tries to execute it using the Data Interop's fme.exe. I only have the stand alone version of FME. I tried using the full path to the fme.exe i'm trying to use, but that's not working ( the script runs without error, but nothing is produced). Any ideas on how to get this script to work, using the stand alone version of FME, as an ArcGIS script tool?


title: Spatio-temporal block kriging with R package gstat?

text: In the documentation of `sp` I saw that the `coords` argument of the `SpatialPolygons` function only accepts 2 columns. Is it not possible to create 3-dimensional polygons? I am doing spatio-temporal kriging using package `gstat` by adding time as a third spatial dimension of a `SpatialPointsDataFrame`. I would like to extend that to block kriging with polygons as blocks. Is there any way to do that?


title: How to restrict access to GeoServer webapp to Localhost in Tomcat 7

text: I am currently running Tomcat 7 on a Windows Server 2008 machine. I have GeoServer 4.2.1 installed as a `war`. What I need to do is restrict access to GeoServer Manager only to Localhost. As it currently stands you can access the Manager anywhere by going to               https://www.mydomain.com/geoserver      I would like to only give access to Geoserver Manager by means of               https://127.0.0.1/geoserver      I have updated by `context.xml` file to limit access to Tomcat Manager via localhost with this:               <Valve className="org.apache.catalina.valves.RemoteAddrValve"allow="127.0.0.1" deny=""/>      Would I do something similar? Any help is greatly appreciated.


title: Georeferencing PNG against Google maps with openlayers plugin

text: I am trying to georeference blueprints to the footprint of a building. I'm using the Google Satellite image from the plugin to do this. The blueprints are rotated about 30 degrees. When I run the georeferencer after placing GCP, the blueprints load in vertically, but the rotation doesn't occur... Why is this? I figured it would line up the GCPs properly... Thanks


title: Applying pixel reliability mask to time series MOD13Q1

text: I am working with a 23-band GEOTiff in which each band is a 16-day MVC MODIS NDVI image. The MODIS .hdfs include a "pixel reliability" image with the following values:                    [-1]: Fill/No Data-Not Processed.           [0]:  Good data      - Use with confidence           [1]:  Marginal data  - Useful, but look at other QA information           [2]:  Snow/Ice       - Target covered with snow/ice           [3]:  Cloudy         - Target not visible, covered with cloud      I want to apply a quality mask using ENVI to elimate pixels with values of -1, 2 and 3. I know how to do this for a single image using the "Build Mask" function in ENVI. If I use the 23-band NDVI.tiff file as my input and use the 23-band pixel reliability file as my mask, will this mask each NDVI image with its respective pixel reliability (i.e. January 1st NDVI will be masked with the January 1st pixel reliability mask), or will the "Build Mask" perform some other operation? If I am thinking about this in the wrong way, is there a better way to apply the quality mask to the time-series stack? Many thanks! Emily


title: How to Share Open Source Python Code for ArcGIS?

text: My company just purchased ArcGIS software. I'll be writing scripts for ArcGIS in Python. Is there a community or open source project where this code can be shared?


title: Python: Change SDE data source for mxds in specific folder

text: I've searched all the forums and help and have made progess but recently got an error that has thrown me for a loop - also some mxds are not being changed at all! What I'd like my code to do: Iterate through mxds in a specific folder, change the SDE data source in each layer within an mxd and change it to the new SDE source (we just upgraded our server). Then I'd like to either save over that mxd, or save to a new folder (whatever works!). Here's the code:               import arcpy     import os          folderPath = "G:\GIS\Services\Test"     for fileName in os.listdir (folderPath):         fullPath = os.path.join(folderPath, fileName)         if os.path.isfile(fullPath):             basename, extension = os.path.splitext(fullPath)         if extension.lower() == ".mxd":             mxd = arcpy.mapping.MapDocument(fullPath)             print "MXD: " + fileName             arcpy.env.workspace = fullPath                         mxd.replaceWorkspaces("", "NONE", r"Database Connections\Connection to gisserver.sde","SDE_WORKSPACE")                  print "successfully changed data sources"     mxd.save()     del mxd      This worked sucessfully when tested with one mxd in a folder, so i thought it would apply to a folder with more mxds, however   The error I received while attempting to iterate through the folder: > Traceback (most recent call last): File "C:/Users/lyee/Test2", line 10, in > mxd = arcpy.mapping.MapDocument(fullPath) File "C:\Program Files > (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\mixins.py", line 608, in > **init** assert (os.path.isfile(mxd) or (mxd.lower() == "current")), > gp.getIDMessage(89004, "Invalid MXD filename") AssertionError: Invalid MXD > filename.


title: Unknown gridlines on old Greenland base map!

text: Can anyone help me out please? I'm using an old 1:10,000 base map of Greenland for a mapping project with nice 10cm gridlines already drawn on. However it uses a coordinate system which is unfamiliar to me, it has values of 92 up to 112 on y axis and 124 to 138 on x axis, which increase by increments of 2. My mapping area is around 61 N, 45.15 W, can anyone please tell me what the grid system my map uses represents? Thanks Eleni


title: How to view DGN file in ArcMap?

text: I have a DGN file made by Workshop Paint Pro and would like to view it in ArcMap 10 so that I can create a risk map. What are the procedures I need to follow?


title: Why doesn't the Kernel Density Estimate function appear in AniMove plugin?

text: can anyone tell me why the Kernel Density Estimate function does not appear under my AniMove plugin in QGIS? I can only perform the MCP. I run Windows OS.


title: Adding a GIS feature to existing database entry with QGIS

text: I've built a database with postgreSQL with tables, which include geometry columns. When I add some data to my non-geometry columns, I can see them in QGIS when I connect to my database using postgis. Now, I want to add polygons to the existing entries with the editing tool in QGIS. That means, filling up the geometry columns in QGIS. How can I do this? I've only heard about adding whole shapefiles to the database but not to add single features to a single column. An example... I have the following table (municipality):               id smallint NOT NULL,     name varchar(50) NOT NULL,     mun_geom geometry(Polygon, 31468)      I can fill the columns id and name in postgresql when I'm connected to my database (like (1,'municipality1')). After adding municipality as a postgis layer I want to complement the entry with adding a polygon feature.


title: No module named PyQt4

text: I'm totally new to python and am getting the above message when I run my script. Can't find any basic step-by-step instructions on how to import this module. Just wondering if anyone could give me some hints? Thanks!


title: How to Resolve Inconsistent OpenLayers vector position Over Google map

text: I have an OpenLayers map with a Google base layer and an OpenLayers vector layer on which people can draw their own bounding boxes using the DrawFeature control. I have found that, as you zoom in, the vector is not plotted covering the same co-ordinates at each zoom level. Rather it seems to zoom at a different ratio to the underlying layer. Strangely, if I resize the browser window, the vector then stays at a consistent position through all the zoom levels (i.e. becoming larger or smaller but covering the same geographic range, as one would want it to), but slightly offset from the position it was initially drawn at. The remains the case even if I make the browser window full size again. See example below:   1. I draw a vector touching the Belgian border. ![first example](http://i.imgur.com/qRams1v.jpg)   2. I zoom out and the vector is now well within Germany. ![second example](http://i.imgur.com/QFTktl9.jpg)   1. I resize the browser window and the box anchors to a consistent set of co-ordinates but offset from its initial position, now overlapping the Belgian border ![third example](http://i.imgur.com/Rxu0z1r.jpg) Please note that this is not the same issue as the vectors scrolling more slowly than the map when dragging the map. The vectors seem to stick to the same area when dragged, but zooming causes them to move relative to the underlying Google layer. Is this a known issue and can anyone help resolve? In answer to the question about my initialisation of the layers, my code is as follows:               var miny=-2000000; var minx=6000000; var maxy=600000; var maxx=10003750;              // allow testing of specific renderers via "?renderer=Canvas", etc     var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;     renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;                    //define the layers needed         g_phys  =   new OpenLayers.Layer.Google("Google Physical",{type: google.maps.MapTypeId.TERRAIN, numZoomLevels:20}, {isBaseLayer:true});         g_sat   =   new OpenLayers.Layer.Google("Google Satellite",{type: google.maps.MapTypeId.SATELLITE, numZoomLevels:20},{isBaseLayer:true});         bboxes  =   new OpenLayers.Layer.Vector("Bounding Boxes", {rederers:renderer});                    //delete feature on select     bboxes.events.on({"featureselected":function(e){bboxes.removeFeatures(e.feature); erasefeature(e);}});     //refresh on visibility change     bboxes.events.on({"visibilitychanged":refreshMap});     //do stuff on modification of vector     bboxes.events.on({"afterfeaturemodified":modbox});               //create new set of option for the map     var options = {                 projection: new OpenLayers.Projection("EPSG:3857"),                 units: "m",                 maxResolution: "auto",                 maxExtent: new OpenLayers.Bounds(miny, minx, maxy, maxx),                 theme: null,                 layers: [                     //define the base layers for the map                     g_phys,                     g_sat,                     bboxes                              ]                  };               //create new map object based on "map" div     map=new OpenLayers.Map('map',options);


title: Model Builder - unique name segment creation

text: Is it possible for us to give each output of a process a unique name segment (inline name). So for example I have a process that dissolves, then eliminates and buffers a single input file. I need to save the itermediate steps into a gdb. So we have [TownNameX_slope]--Dissolve[\townDissolves.gdb\TownNameX_dissolve]--eliminate[townEliminates.gdb\TownNameX_eliminate] etc. then TownNameY etc. So in a batch process I have to keep entering TownNameX for each output type. I know in iterators you can use %n but this is for just the output of one tool. Can we extend this to all the output where %n=townnamex and then for dissolve the output is %n_Dissolve etc I guess the other option is to run it in python in a "for" loop but I would prefer to stay in model builder as the reporting is easier. This is similar to the question in ArcGIS 10 Model Builder In-line variable substitution for input data path Thanks


title: Mapfile + OpenLayers: accessing the database by column names (hover, onclick functionality) without using MapServer templates?

text: What I'm trying to accomplish: take the functionality below (see picture) and do the same without using a MapServer template (if this is even possible). I currently have a mapfile, and an openlayers file. All the data is being displayed correctly. In my mapfile, I reference a database. I call out that database by a column name so what is what's used for the label on the map. Quite simple. Now, the thing here is that I'm accessing just a standard php/html file. Not a mapserver template, like:               http://192.168.2.164/cgi-bin/mapserv.exemap=%2Ftmp%2Fnavhome%2Fnavcore%2Ftemplates3%2FAllAirports.map&mode=nquery&imgext=-180.000000+-107.943706+180.000000+107.943706&imgxy=640+384&zoom=0&qlayer=Cities&layers=world&searchmap=true&img.x=705&img.y=347      Example: 123.456.7.890/mymap.php in the browser is what I want. Not the above. Now if I wanted to accomplish this using a MapServer template (correct me if I'm wrong on my terminology), this has proven to be quite easy, because you can reference               map_body.html     map_header.html     map_footer.html      in your Mapfile. In those files you can do some javascript and whatnot, reference the columns in the DB and it works. See below for an example of this (top left is the onClick functionality when you click on a point on the map, where the mouse is is just the mouseover). **See below: this is what I want to do without the cumbersome MapServer functionality, ie sluggish navigation buttons** ![this is what I want to do JUST using a MapFile and OpenLayers, no MapServer template](http://i.imgur.com/tonZR81.png) So what I'm trying to wrap my head around is: how do you set up the database references in your .html/php file? Or, how do you reference the templates (map_body.html, map_header.html, map_footer.html) from your MapFile in your OpenLayers file?


title: URL Hyperlinking

text: I would like to have the URL from an attribute (in the attribute table) as a 'hotlink' wherein you can click on the URL and it opens the browser with the data. As it is right now, I'm having to "toggle edit" the attribute table, R-click on the feature, open form, and copy/paste the URL into a browser. I don't write python code so does anyone have a quick workaround for me? (I can paste python script). Thanks!


title: Delete feature from geographic database

text: How do I delete a selected feature from database when it is clicked on in the ArcGIS API for javascript?


title: Wireframe display of terrain in ArcGlobe?

text: I have successfully imported a terrain dataset into ArcGlobe which provides elevation values to the globe surface. Is there an easy way I can now change the symbology/display of the terrain/globe surface to a 'wireframe' style display, ideally with customisable grid spacing etc? Or would I need to drape a grid image over the terrain surface?


title: How to specify field name in raster to feature conversion tool or intersect tool?

text: **What I want to do :** input: raster file, polygon feature layer process:   1. Convert raster to polygon using `IConversionOp.RasterDataToPolygonFeatureData` Method   2. Intersect resultant polygon feature class and input feature layer using Intersect tool from my c# code. I am doing these two steps perfectly. But I need some specific field from intersected layer for next steps. So it would be helpful if I can specify field names in intersected layer. But I could not find any way to specify field names. It would be helpful if anyone can give me some idea how to specify some field names. If there is any process to take fields from one feature class and not from another feature class in intersect tool, that will also solve my problem. Or when I am converting a raster to a feature class, can I specify the value field name? Or is there any rule for value field name in the newly created feature class from the raster layer. By default the field of this value field in feature class is "GRIDCODE".


title: How do algorithms load pieces of maps?

text: If the user moves the mouse pointer to the top left corner of the screen, I want my software to auto-load the three tiles that are closest to the mouse pointer. My assumptions are that the smaller are the map tiles to be loaded, the lesser tiles I will have to load. The question is,   1. whether there is any existing logic that I can use for this loading (something people have used already, since I need a flexible way to ) and ...   2. ...a flexible way to store the loaded tiles in C++?. My app is in C++, and I'm planning to store each tile as either a class or a struct. Is there an existing way/a better way to do it? **EDIT:** Responding to Dan: I'm still in fact-finding stage, and currently only have a map with roads stored as lines, where each line is composed of a list of line segments specified by a starting lat-long and an ending lat-long. I'm open to using any tiling scheme, but in a simple way, I'm visualizing the map to be divided into a grid of small squares, so the squares/tiles close to the mouse pointer and containing road segments should be loaded into memory. My algorithm would then process the loaded tiles to find the angle of the road and to check if the road forms a junction with other roads etc. My algorithm can't afford to wait for too many loading operations from the hard-disk. The priority is knowing how to store information so that I'll have them in the RAM as quick as possible. Whether to pre-calculate the tile area and store the tiles in a database or whether to extract square portions from the map and cache them in the RAM was a decision I was trying to take when I posted this question.


title: how to export the entire bounding box from overpass-turbo?

text: I'm new to the overpass-turbo world, so sorry if the question sounds silly. What I'm trying to do is to export all the elements (peak, tracks, routes, streets.....) of a bounding box I have selected. I'm trying to do this with the wizard, but within it I have to choose specific tags and so export one element per each time, right? Thanks!


title: Error G_Calloc v.to.rast

text: Using the grass GUI I get this error message: ERROR: G_calloc: unable to allocate 2004389888 * 8 bytes of memory at raster.c:71 What can i do to remedy this?


title: Bearing Distance to Line tool from dBASE table only draws one line in ArcGIS 10

text: I am using the **Bearing Distance To Line** tool to draw geodesic lines to points defined by a bearing and distance. This is an aside, but I will then use the **Feature Vertices To Points** tool to convert the lines to points. I have a dBase file (because my Excel and CSV tables were failing for some unknown reason which may become the subject of another question) with the bearings in decimal degrees and the distances in meters. All of the X and Y points from which the lines are drawn are the same because the bearings and distances were all taken with a total station from a single datum point. So the X and Y fields of my table are the same value for each entry. The problem arises when I run the tool -- I only get a single output line drawn when there should be 113 individual lines with different bearings and distances from the same start X,Y coordinate. Has anyone encountered this issue before or have an idea how to resolve it? Does it have something to do with the start X and start Y being the same value for each table entry? It shouldn't be the case, but am I missing something? I am not running a custom Python script, just using ArcMap 10.0. I added the dBASE file from ArcCatalog and ran the Bearing Distance to Line tool with the dBASE file as my input file. I am expecting to see something similar to the top graphic in this **post**. To work around the problem of the Bearing distance to line tool, I found and used this **script** from Susan Jones which did work with a feature class created from the data in the dBASE file. However, I'd still like to be able to figure out the Bearing Distance to Line problem. Thanks! (As an aside, see my other question for a related, though more general question)


title: Move the cursor backward - opposite of row = rows.next()

text: I'd appreciate if anyone could please tell me if there is a way to go up one row in python cursor. opposite of row = rows.next() Thank you!


title: How to determine the date of the Bing imagery for an area?

text: Is it possible to determine how current the Bing imagery for a certain area is? I've checked the community blog and searched through the monthly updates, but I could not find my study area. I found posts where they mentioned updating the imagery for the whole country, but with no specifics (is it really the entire country? Or just most of it?). In the credits in `Layout View` in ArcMap 10.1, I noticed that while Microsoft Corporation is `2012`, the other credits are `2010`. Does this indicate that the imagery is from that year?


title: Query OSM ids by name in PostGIS

text: I installed PostGIS on a server and now I want to query the OSM ids referring to a name, in order to get a list of all geometries according to this name. Is there an efficient way to query this ids in PostGIS or do I have to use a tool like nominatim? I am new to OSM and PostGIS.


title: How to change the zoom on an embedded google mymap?

text: I embedded a google mymap in my website using the embed code they provide. However it is zoomed much too far out. How can I zoom it in? I added the second character (undrscore) so it would display. <_iframe width="425" height="350" frameborder="0" scrolling="no" marginheight="0" marginwidth="0" src="http://maps.google.com/maps/ms?ie=UTF8&hl=en&source=embed&msa=0&msid=212526618623367333770.0004a2dd56ade1cf3f7eb&ll=34.931542,-81.894879&spn=0.043299,0.140592&output=embed">   View Map Name in a larger map


title: How to export one shapefile per shape in shapefile?

text: I have road shapefile consists of Road No. Suppose i have a road network consisting of 125 roads. I need to get 125 individual shapefiles from one road network shapefile without select one by one. Is it possible to do?


title: Separating a shapefile into separate shapefile based upon their features using a python script

text: I have a shapefile of a US. I need to make an individual shapefile for all the states in the US using a python in ArcGIS. How can I do it?


title: Converting Categories to shapefiles

text: I have created 2 Categories in a feature using classification.Now how do I convert them into new shapefiles?


title: Satellite Imagery Providers

text: I'm building a web application - my requirements are:   * Must have _global satellite imagery_.   * Must have easy to use API's that are capable of geoJSON data imports, mouse and keyboard gestures, etc. To my knowledge Bing, Google and Yahoo are the only companies that supply both of these. Some support the former to an extent; perhaps not global (NASA, third party mapping solutions), and some support the latter (OSM for example)...but few support both. Are my only options Bing, Google, and Yahoo? **If so, is anyone aware of any ongoing projects (satellites) that are under construction that will eventually be launched and provide global satellite coverage under open source?**


title: How to learn ArcObjects (.NET, C#) for desktop applications but having only Arcmap desktop license?

text: How can I learn .NET and C# for desktop applications having only a desktop license? I want to upgrade myself in the development side so I'm trying to learn this stuff and I'm searching all the documentation on the Esri sites; everything is there to learn but I can't test in my machine. Is there any way to learn without licenses as a beginner to get the concepts of programming in GIS? I know that I should have Visual Studio in my computer but I couldn't find 2008 Express edition on the Microsoft website. Even if I have VS 2008 Express edition will it work for me in this case? I'm using Arcmap 10. Please treat me as a complete beginner in this topic and guide in that way.


title: can arcpy detect whether it's running in Desktop or Server?

text: I'm developing some code that I want to run both as an ArcGIS Server web service and as a standalone script. I want to be able to modify the execution slightly depending on whether it's running as an AGS web service or as a standalone python/within Desktop. Is there a way to detect this? I've done some searching but haven't found anything very satisfactory. I have seen arcpy.mapping.CreateGISServerConnectionFile, but this sets connections to a server as opposed to checking the current environment. I've also looked for information on how services are run w/in AGS (e.g., this), but I'm not really seeing anything along the lines of an arcpy.isWebService() type of function. Do I just need to check the environment (like to see if "%scratchFolder% has been set or examining something related to sys.argv)? Seems kind of an unstable solution.


title: Exchange Boundaries and Wire Centers in France and Germany

text: I am attempting to find a shapefile of wire center service areas / telephone exchanges in France and Germany. Ideally, a point file showing the location of the actual equipment, coupled with a polygon that shows each center's service area. To see an example of the what I'm talking about, this is the U.S. equivalent: http://www.navmart.com/wire_center_boundaries.php Obviously, free is best, but I open to sources who might have this available for purchase. Any and all help is appreciated!


title: Python select by attribute using a variable counter

text: I am brand new to python and trying to learn how to do some simple loops but I have hit a snag. Below is what I am working with, I am trying to select the polygon row by row, but I keep getting EOL errors and various other syntax issues as I try to troubleshoot. If I remove `""" findMe """` and replace with `""" "pgNum" = 1 """` it works fine. How can I change 1 to be x? I am using the python window within ArcGIS Thanks.               import arcpy     workSpc = (r'C:\Users\adams\Downloads\tmp')     smooth = 'MNR_89_Smooth'     grid = 'test_grid'     x = 1          if x == 1:         findMe = '"pgNum" = ' + str(x)         print findMe         arcpy.SelectLayerByAttribute_management("test_grid","NEW_SELECTION",""" findMe """)         x = x + 1         print x


title: Webmapping option that allows rasters?

text: I'm trying to create an introductory GIS course that will involve an ideally simple web mapping component. I was underwhelmed that ArcGIS online doesn't support rasters, and surprised that MangoMap and GeoCommons don't appear to either at the moment. Is anyone familiar with a webmapping site with a fairly simple UI that allows the uploading of raster data in addition to vector? Many thanks, David


title: Embed an ArcGIS Map in a Crystal Report

text: Is it possible to embed a map image inside a Crystal Report using ArcGIS 10? I have seen articles about creating reports based off of map fields, but I would like to include an image of a selected area (i.e. a specific parcel on a map). Is there a way to do this?


title: How to read this GRIB data

text: on this open FTP server (ftp://ftpdatos.aemet.es/modelos_numericos/hirlam005_peninsula/) you can find HIRLAM models from the Spanish meteorological service. The documentation in the LEEME files (in spanish) says the data files are in GRIB-1 format. Unfortunally I'm unable to read it. I use differet programs, like QGIS, and commaind line tools: wgrid, grid_ls, etc without success. Anyone can take a look and let me know how read it. Thanks in advance. === Updated === I'm new to grib and I dont understand about the codes and keys. The documentation at ftp://ftpdatos.aemet.es/modelos_numericos/hirlam005_peninsula/LEEME_DATOS.txt talks about IPAR, LTYP and LEV ids. Also it contains some description about the parameters in the files:               Variable                   || Unidad  || IPAR || LTYP || LEV      Temperatura a 2m                 K        011     105     2      Componente x del viento a 10m   m/s       033     105     10      Componente y del viento a 10m   m/s       034     105     10      Precipitación total            kg/m2      061     105     0     Humedad relativa a 2m          kg/kg      052     105     2       Can someone help me reading that data ??? === UPDATED === Thanks to all for your replies. I finally found the problem which is a bit embarrassing :) A was downloading files in ASCII mode instead BINARY. I can confirm the grib_api tool, once compiled and installed works perfectly. Example:               $ grib_ls 20120220000000_sfc_fc15     20120220000000_sfc_fc15     edition      centre       typeOfLevel  level        dataDate     stepRange    shortName    packingType  gridType          1            lemm         heightAboveGround  2            20120220     15           2t           grid_second_order_SPD1  rotated_ll       1            lemm         heightAboveGround  10           20120220     15           10u          grid_second_order_SPD1  rotated_ll       1            lemm         heightAboveGround  10           20120220     15           10v          grid_second_order_SPD1  rotated_ll       1            lemm         heightAboveGround  2            20120220     15           r            grid_second_order_SPD1  rotated_ll       1            lemm         heightAboveGround  0            20120220     0-15         unknown      grid_second_order_SPD1  rotated_ll       5 of 5 grib messages in 20120220000000_sfc_fc15          5 of 5 total grib messages in 1 files      On the other hand, I'm not able to make work gdalinfo with this GRIB files, but for the moment I don't require it. Thanks !!!


title: Pyqgis 2.2 : get signal on feature added : error message

text: For my pyqgis plugin 2.2, I try to: * * *               def logFeatureAdded(self,layerId,addfeatures):          message = str(layerId) + " has features added: "          for feature in addedFeatures:              message += str( feature.id() ) + ", "          QgsMessageLog.logMessage(message)          QApplication.beep()               layer.featureAdded.connect(self.logFeatureAdded)      * * * When I try to add a feature, this following error message appears: > TypeError : logFeaturesAdded() takes exactly 3 arguments (2 given).


title: Creating a Bing and Google Tile Layer for Leaflet

text: I want to use the Leaflet Javascript Library to display both a Google and Bing Tile Layer, though I'm unsure how to construct the url. The documentation shows an example using OSM:               var cloudmadeUrl = 'http://{s}.tile.cloudmade.com/YOUR-API-KEY/997/256/{z}/{x}/{y}.png',     cloudmade = new L.TileLayer(cloudmadeUrl, {maxZoom: 18});      But for some reason I can't duplicate this for either Bing or Google. Any ideas? Thanks.


title: towgs Definition for European Datum 1950

text: I am new to all things GIS so apologies for a basic questions. I am trying to specify the CRS of a climatic raster layer in R. It is an equal area projection based on Europe Equidistant Conic and the datum is D European 1950. Can anyone help me specifiy the argument in R? projection(myraster) = CRS("+proj=eqdc +ellps=???+ datum=???") Many thanks if you can help!


title: Updated: Openlayers Vector Layers with Google Map Base Layer (via Grid): Vectors Stick but can't get correct resolution/projection

text: Revamping this question instead of opening a new one. Previous issue: vector layers were not "sticking" to the map, so when you would zoom in or out the vector layers would "come out" at you. Note their lat and long didn't change, but it looked pretty tacky when zooming-as they would not stick to the map. It seemed like the only way to do this (the stick factor) was to use a WMS or Grid Layer: correct me if I am wrong. I chose Grid in this instance. Note: the "stick factor" is not an issue when the base layer is OSM or OpenLayers, only Google base layers. Here is an example demonstrating this at bottom of page: The code I currently have for this Grid is currently this:               var options = {                 resolutions: [1.40625,0.703125,0.3515625,0.17578125,0.087890625,0.0439453125],                 maxExtent: new OpenLayers.Bounds(-180, -90, 180, 90),                 //restrictedExtent: new OpenLayers.Bounds(-360, -90, 360, 90),                 numZoomLevels: 8,                 singleTile: true,                 isBaseLayer: true,                 wrapDateLine: true,                 controls: [                 new OpenLayers.Control.OverviewMap(),                 new OpenLayers.Control.Navigation(),                 new OpenLayers.Control.ScaleLine( { maxWidth: 150 })                 ],                 getURL: function() {                         var center = this.map.getCenter(),                         size = this.map.getSize();                         return [                         this.url, "&center=", center.lat, ",", center.lon,                         "&zoom=", this.map.getZoom(), "&size=", size.w, "x", size.h                         ].join("");                 }         };      Now the first thing that comes to mind is that I feel I should have an `options` variable for the map object and an `options` variable for the Grid layer (what are your thoughts on this?), as shown below:                var grid = new OpenLayers.Layer.Grid(             "Google Streets",             "http://maps.googleapis.com/maps/api/staticmap?sensor=false&maptype=roadmap",             null, options         );         map.addLayers([grid]);      The left picture below is how the resolution was before: just using a standard OpenLayers Google layer. The picture to the right is using the Grid layer and point to the Google default layer: mentioned in the code above. ![OpenLayers Google Layer, no grid layer](http://i.stack.imgur.com/CommT.png) ![bad resolution using Grid Layer and accessing Google Maps](http://i.stack.imgur.com/zprfO.png) I just posted a snippet of Africa to show the resolution issue, and as for the points (all my points are in the US, they are completely out of whack also (the left most point on the map should be on the west coast of the USA, while the right most point should be on the east coast of the USA.: ![enter image description here](http://i.stack.imgur.com/Sbm51.png) 1) **Should I create a separate options variable for the grid layer?** 2) **Is adjusting the resolution within the options the correct way to fix this?** 3) As you can see in `getURL` there is no mention of projects or transformations: **how do I account for the points showing up correctly in the USA?** As always, any and all input is appreciated.


title: Clustering polygons by sum of values

text: I would like to group polygons that (1) are contiguous and (2) when the combined values of field X >= Y. Specifically I am grouping zipcodes that are connected. a group is when 1 or more zipcodes 'population' field equals 100,000. Can someone suggest an existing ArcGIS function or python method for doing this? Many thanks in advance! ---update---- I wrote a python algorithm to group the polygons using their respective neighbors as grouping members until threshold is met. it works nicely.


title: QGIS (getting 2 images into kml popup for each point)

text: I've got a csv file with some points with attributes. 2 of the fields reference a .jpg file which are in the same folder as the csv file. Using QGIS, how do I export a kml so that the popup in Google Earth will show the jpg's associated with each point?


title: Data loss after compressing geodatabase in ArcSDE

text: I suffer sometimes from data loss after compressing ArcSDE geodatabase, this geodatabase contains versioned feature classes, so the compressing script reconciles and compress the geodatabase:               # Execute the ReconcileVersions tool.     arcpy.ReconcileVersions_management(db_path, "ALL_VERSIONS", "sde.DEFAULT", versionList, "LOCK_ACQUIRED", "NO_ABORT", "BY_OBJECT", "FAVOR_TARGET_VERSION", "POST", "DELETE_VERSION", "c:/reconcilelog.txt")          # Run the compress tool.     arcpy.Compress_management(db_path)      it seems that the main reason is a bad use of `ReconcileVersions_management`, what are the parameters recommanded ? and how i can test if data loss would happen again ? the geodatabase is based on ArcSDE 10.1 and Postgre SQL 9.0.5.


title: How can i create single feature fwith Combining Selected feature qgis?

text: I am trying to create a single Qgsfeature from the Combination of selected features in Pyqgis. Currently i create the feature using a loop over selected features and combine them 1 by 1. but this process froze qgis UI. does anyone can suggest me another idea ?


title: Pyqgis 2.2 : zoom to layer extent?

text: Regarding my plugin pyqgis 2.2, i can load a vector layer succesfully. But when i try to zoom to the layer extent it does not work. How to do in order to zoom to the layer extent? * * * Thank you very much for your replies. I have just tried your proposition.I pasted your piece of code into my plugin. But unfortunatly, when i load a vector layer. I try to zoom to the layer extent (thanks to right click -> zoom to layer extent). I have no error message but there 's no change zoom to the extent of vector layer. You can read the python.py (github plugin link) here : https://github.com/picardie-nature/clicnat-qgis/blob/master/plugin.py Any thoughts? THank you very much.


title: Open shape file in ArcGIS Runtime WPF SDK as Graphics/Feature layer

text: I googled and followed this tutorial http://www.arcgis.com/home/item.html?id=953a530ada30477d95cc69e788268bc9 and figured out how to import Shape files and Rasters into Runtime. However, it seems most of the functionality uses a FeatureLayer/GraphicLayer instead an ArcGISLocalDynamicMapServiceLayer that the sample uses. My question is if there is a way I can import the shape file as a graphics layer or if i can somehow cast my ArcGISLocalDynamicMapServiceLayer as a Graphics layer because shape file is essentially the same as a shape file. An example is i'm trying to use the FeatureDataGrid() to show the attribute table of my shapep file, but I can't since it requires a graphics layer, not an ArcGISLocalDynamicMapServiceLayer. Thanks and Regards, Kev


title: OpenStreetMap data at different scales

text: Could anyone just clear up a question regarding openstreetmap data. Are there different OSM data sets for the various set scales, or do the data get generalized on-the-fly when the map is zoomed out to a smaller scale? I understand that the many objects are hidden in smaller scales, but are there also different generalized datasets for the smaller scales as well? Thanks, Robert


title: How do I add a Corel Draw map as .dxf file to QGIS and how do I georeference it?

text: I have a vector data map of Poland, which I exported from Corel Draw. I can display the map in QGIS. At first the boundaries and subdivisions are only displayed as lines, after saving the data as shapefile, I was able to use the "Line to Polygon"-tool to convert the lines. Now I have a shapefile with multiple features. Now I want to relocate and distort it to make it fit to my basemap. As you can see in the picture, the Poland shapefile is presented close to Africa because the location is not yet defined correctly. ![poland](http://i.stack.imgur.com/1JPSO.jpg) There are now several problems:   * I tried to move the Poland shapefile in an edit session but the shapefile contains of multiple polygons and I was not able to move them together. How can I move multiple features at the same time?   * I also tried to make a single polygon by merging the multiple features in an edit session, but when doing it with the "Merge Selected Features"-tool, an error occurs and QGIS crashes.   * Therefore, I tried it with a third way by creating a new shapefile, drawing a polygon, that covers the whole Poland polygons and then cliping both shapefiles. After Cliping, the new shapefile is created and added to the table of contents, but there is nothing displayed on my map. As you can see, there may be a mistake in the previous steps but I dont know where, I already repeated every step but there is still no solution found. I suspect the problem to be in the step of converting the Corel Draw .dxf file to a shapefile or the transformation from lines to polygons. Also if you know a more effective way to get my map to its planned location and how to reference it correctly, please let me know your thoughts.