title: Modeling data in PostGIS when Several Join queries with Geometry are Needed

text: Suppose you have a database where you have several tables with geometric information and you have several cases in your system where you need to query the tables with a lot of joins (ie: `LEFT JOIN table2 ON st_intersects(table1.the_geom,table2.the_geom)`). I believe this type of query is slow compared to a "normal" int field join, even if you use spatial indexes on the fields. In a system where this situation arrives a lot, how do you model the data? For example, I have many situations where you would have to query this, joining with 10 other tables. Each table holds a different type of information, but I need to see if a row in table1 intersects with the geometry of all the other tables. I thought about creating a many-to-many relationship, relating the ids to make the queries run faster. But then I would need a trigger to each time a new row is inserted, check which ones it intersects on each possible relationship (10, in the example) and save it. This would make the inserts a lot slower too. Is there a common pattern to solve this type of problem on PostGIS?


title: Get latitude and longitude of each point forming a MultiLineString

text: A MultiLineString is formed by a sequence of points and each point has latitude and longitude. Use the sql below to get the latitudes and longitudes of each point so that my MultiLineString:               SELECT         St_AsText (ST_Transform("public"."Shapes".geom,4326)),         "public"."Shapes"."id"     FROM         "public"."Shapes"      ![Result](http://i.stack.imgur.com/xXoK4.png) I want to get the sequence of points of each line in separate rows and lat / long in separate columns example:                  ID            pt_lat            pt_lon        1001             lat1             lon1        1001             lat2             lon2        1001             lat3             lon3        1001             lat4             lon4        1001             lat5             lon5        1001             lat6             lon6        1002             lat1             lon1        1002             lat2             lon2        1002             lat3             lon3        1002             lat4             lon4        1002             lat5             lon5        1002             lat6             lon6      I know I can do that consuming one GeoJSON for example, but want to know if there is a native function of PostGIS for this.


Vectorized few lines using QGIS with Google Maps as a base. If I do not Trabachar with maximum zoom does not respect the topology of the base. It is possible to make snap using Google Maps based on QGIS?


title: OpenLayers 3: context specific styling

text: I am having trouble with OL3 extracting context specific info from a nested JSON object. My feature's JSON goes like this: feature.attributes.geometry.name = "Point Alpha" How do I display the name as label when styling the feature? I tried ol.expr but it only parses the first level properties. Please advice!


title: Map Layout Elements In New Layer Automated Using Python

text: I was wondering if there was a way to put a title and integer data from an attribute table in a new layer created in Python. I was able to create a new layer using copy features, zoom to the selected features, and export the zoomed features in a map as a PNG all using arcpy, but I also wanted a title at the top and data from the attribute table of that specific layer in the map layout before it was exported. Is it possible to get all this automated using Python?


title: Can ArcGIS only merge two shapefiles each time?

text: I need to merge multiple point shapefiles (in the same directory) into one point shapefile. After inputting these shapefiles and merging it, I found the final output only contained the data from two shapefiles. How to merge multiple shapefile in ArcGIS at once instead of merging only two at each time?


title: How to make a certain attribute Toggelable?

text: I have imported data from an excel spreadsheet containing lat,long coordinates associated with other attributes such as current operator of each point (well). I imported as a tab delimited file, then saved this as a shapefile. I want to be able to organize the files by operator so that they can be toggled on and off (by operator). I have been doing this by opening the attribute table, making a selection and then saving the selection as a new shapefile for each operator, which works, but I have reached the limit of shapefiles for a single project, and am wondering if there is another way to accomplish this.


title: How to set up the layer order with Google Maps API?

text: I'm developing a web mapping site (using google maps api 3.0 and javascript) that combines wms layers and ground overlays (uploaded rasters) that are displayed on top of google maps. The software is working well except that I'm having problems controlling the draw order of the layers. I would like to have the wms layer (the NRCS soils layer) displayed on top of the custom raster images, with the google maps as the base layer. Currently, the wms layer displays as expected on top of the google maps layer, but is covered by the raster layer. The question is: does the google maps api allow control of the order that layers are displayed (in my case the vector layer on top of raster layer on top of google maps)? I have tried setting the zindex of the display order, but that has not worked (but I could easily be missing something).


title: MODIS Level 1B batch process

text: I am a new MODIS Level 1B user. I have been using the MODIS Swath tool to process my images. However, I have over 1000 images to process. I was wondering if there a script available in IDL or R (very little coding experience) where I can batch process all the images?


title: How can I find and fix broken raster metadata using QGIS?

text: I have a raster (TIFF) hillshade image that I would like to use in a mapping project. However, I'm having a hard time lining it up properly in QGIS. Through a series of investigations, I have determined that this problem might be caused by broken metadata. The raster was apparently created in ERDAS Imagine. I opened this raster in a text editor, to see if the metadata tags were visible. Sure enough, the following appears at the very end of the file:               #IMAGINE GeoTIFF Support     Copyright 1991 - 2005 by Leica Geosystems Geospatial Imaging, LLC. All Rights Reserved     @(#)$RCSfile: egtf.c $ IMAGINE 9.0 $Revision: 10.0 $ $Date: 2005/07/26 15:10:00 EST $     Projection Name = NAD_1983_HARN_StatePlane_Hawaii_3_FIPS_5103_Feet     Units = us_survey_feet     GeoTIFF Units = feet|IMAGINE GeoTIFF Support     Copyright 1991 - 2005 by Leica Geosystems Geospatial Imaging, LLC. All Rights Reserved     @(#)$RCSfile: egtf.c $ IMAGINE 9.0 $Revision: 10.0 $ $Date: 2005/07/26 15:10:00 EST $     State Plane Zone -5103     NAD = HARN      This data indicates the right projection, with the right units (feet) (this projection is functionally similar to EPSG:3760). ArcGIS is apparently able to read these tags, and readily assigns it the correct projection. However, when I view metadata on this file using QGIS or the QGIS browser I see that the projection of the raster is different. Most of the parameters are correct, but the units are set to meters. See screenshot from QGIS browser below. ![enter image description here](http://i.stack.imgur.com/UevgS.png) Using ArcGIS, I removed the spatial reference of [a copy of] this raster, and then reset it using the Define Projection tool. Unfortunately, the metadata shown in QGIS remained the same. This leads me to two questions.   1. Where is this information coming from? Is there some other set of tags that QGIS is reading?   2. How can I remove or repair this broken metadata? I'd prefer an open source solution, but I'll accept anything that works for now.


title: ArcGIS 10.0 Service Pack 5

text: Another service pack for ArcGIS 10.0 has recently been released. However, previous service pack releases have introduced problems and instability along with fixes. Have there been any reported bugs or known stability issues with the latest service pack release? Also, what are the pros/cons of upgrading to ArcGIS 10.1 over ArcGIS 10.0 SP5? Info Service Pack 5


title: Looking for pedestrian-specific routing engines

text: Most of the applications providing pedestrian routing I see are sufficient for many needs but not optimal: you cannot cross open spaces like town squares, routing is in the middle of the street etc. Do you have examples of websites or applications with routing really specific to pedestrians? Thanks


title: Python Program for "QC of PDFs Created from Data Driven Pages"

text: I need some help. I am fairly new to Python. I found a sample code and I need to change a couple things in it for it to work. I need it to make sure that all PDFs were created from the Data Driven Pages program I have now. Reason is I sometimes have about 150 PDFs created from the DDP program. For Example : Houston has 152 pages in my mxd and DDP is supposed to make 152 pages. I've made some changes already to the code I was needing someone to double check my work.               import arcpy, os, sys     from os import path as p          #Make parameters for people to choose mxd and folder of PDFs to compare     mxdList = arcpy.GetParameterAsText(0).split(';')     pdf_path = arcpy.GetParameterAsText(1)     fc_path = arcpy.GetParameterAsText(2)          report = {}     page_list = []          #Use Search Cursor to go through Attribuite Table to get Sheet number info     for mapDoc in mxdList:         mxd = arcpy.mapping.MapDocument(mapDoc)     for fcClass in fc_path:         fc = arcpy.ListFeatureClasses()                 max_list = []                 rows = arcpy.SearchCursor(mxdList, "Sheet_ID")                 for row in rows:                     max_list.append(row.Sheet_ID)                     sheet_count = max(max_list)                     arcpy.AddMessage('Sheet count: {0}'.format(sheet_count))                 else:                     arcpy.AddError('No Layers in %s match data source'%mapDoc)     pdfDoc = arcpy.mapping.PDFDocumentOpen(path)     print pdfDoc.pageCount


title: Using select by location to update a field in a feature class

text: I am using Arc 10.1 and I am working on a script that takes a value from a field within one shape file and updates a field in another shape file based on its location. Basically I have smaller parcels that reside within larger areas and I want to assign the name of the larger area to a field within the smaller parcels. The issue I am having is that when I run the script, every parcel ends up having the same "area". I believe its because its essentially writing each area to the parcel and then rewriting it again and again until there are no more areas to loop through. Can someone provide some advice?               #Make a feature layer for the neighourhoods     arcpy.MakeFeatureLayer_management(neighbourhood, "neighbourhoodLayer")     arcpy.MakeFeatureLayer_management(parcel, "parcelLayerTest1")     #Setup cursor to update new NBHD field     with arcpy.da.UpdateCursor(parcel, (parcelIDField,nbhdField)) as parRows:         for par in parRows:             # Create a query string for the current patrol zone              #parcelIDstring = par[0]             #queryString = '"' + str(parcelIDField) + '" =' + "'" + str(parcelIDstring) + "'"                  #Select parcels based on location contained within neighbourhood             arcpy.SelectLayerByLocation_management("parcelLayerTest1", "HAVE_THEIR_CENTER_IN", neighbourhood)                  nbhdTable = arcpy.SearchCursor(neighbourhood)             for nbhd in nbhdTable:                 nbhdValue = nbhd.getValue("Neighbourh")                 if nbhdValue:                     par[1] = nbhdValue                     parRows.updateRow(par)


title: Why does onDemand setting remove my buttons from my menu?

text: When I set the onDemand="false" attribute, my buttons disappear from my menu! I'm setting the attribute to false so that I control its behaviour using the OnUpdate override in code. Am I doing anything wrong?               <ArcMap>         <Commands>           <Button id="SayHi_Btn" onDemand="false" class="SayHiButton" caption="Say Hi!" image="Images\SayHi_Btn.png" />           <Button id="SayHello_Btn" onDemand="false" class="SayHelloButton" caption="Say Hello!" image="Images\SayHello_Btn.png" />           <Button id="SayAloha_Btn" onDemand="false" class="SayAlohaButton" caption="Say Aloha!" image="Images\SayAloha_Btn.png" />         </Commands>         <Menus>           <Menu id="ChooseGreeting_RootMenu" caption="Choose Greeting" isRootMenu="true">             <Items>               <Button refID="SayHello_Btn" separator="true" />               <Button refID="SayAloha_Btn"  />               <Button refID="SayHi_Btn" />             </Items>           </Menu>         </Menus>       </ArcMap>


title: Environmental styles library

text: i'm looking for environmental symbols, related to berm diagnostic (collapses, subsidences, anthropogenics deposits... Is there any style library available? Thanks in advance, Max.


title: How to create lines between pairs of two coordinates?

text: I have a table with pairs of points. Out of that I want to create a line layer (that keeps all attributes of the table). One line for each raw from point to point. I want to calculate the length of these lines later. A simplified example table would look like that: ![table](http://i.stack.imgur.com/VoGjU.png) Can I make this in QGIS? I am using QGIS 1.8.0 on a Mac 10.6.8


title: TatukGIS SQLite Spatial vs. Spatialite

text: I must apologize in advance, I'm just getting into GIS and trying to figure things out. I have an SQLite spatial database that I'm pretty sure was created in TatukGIS. I've tried to open it in QGIS but the database format is different than that of a Spatialite database. I would like to find a way to convert the SQLite database created in TatukGIS to a Spatialite Database that can be opened in QGIS. I know it has something to do with the OpenGIS simple features specification and how Spatialite may not fully comply with those specifications, but I don't have enough experience to know how to go about converting it (if it's even possible). Here is a link to the Databases to give you an idea how they are organized. https://drive.google.com/folderview?id=0B2TMeZ6iEFvHa3BBOVJpVVpBY1E&usp=sharing TatukGIS.sqlite is in epsg:3569, Spatialite is in epsg:3857


title: Custom Model should be used with ArcMap or ArcCatalog?

text: I have one basic question , when you create custom model using model builder in ArcGIS. It should be used with ArcMap or ArcCatalog? Reason I am asking this question is I have created a model and I am setting it "Store relative path names" and I have defined a directory structure for the same. Now the problem is in ArcMap when you right click on "ToolBox" it will allow you to add toolbox and you can select your tool from the specific location whereas in ArcCatalog it creates Tool (.tbx) at C:\Users\\AppData\Roaming\ESRI\Desktop10.0\ArcToolbox\My Toolboxes folder. Where I am not able to access the relative scripts because I will not have any absolute path and it should not be absolute. Also can someone share the difference between ArcMap and ArcCatalog? Any idea or suggestions on this?


title: batch attribute table export

text: How can I export multiple raster attribute tables at once? Using the batch grid control tool mentioned here? Batch Export feature attribute to ascii I can't seem to find how to access the tool in ArcMap 10.0. Thanks.


title: What are valid values for "data defined" labeling settings?

text: I'm using QGIS 1.8. I'd like to customize label features on an existing shapefile using the new label engine and it's "data defined settings" option. I'd like to know what the attribute table field types should be for each of the settings: All 7 font options, all 2 buffer options, and all 6 position options. Should they all be numeric fields, or some text, or both? I want to create these fields in the attribute table correctly so that I can map them to the label field settings. I can't find information on the proper field types anywhere. Thanks for your thoughts.


title: Is segment of polyline contained in polygon

text: I am pretty new to arcpy and am slowly learning all its functions. I have a scenario where I am traversing a polyline. I get a start point and an end point on the polyline and form a segment with those points. That segment is used to create geometries to construct a polygon around that area of the polyline. My problem is that on areas where the polyline curves, a part of the polyline is outside of the polygon. Ultimately I will need to shift my polygon downward/upward depending (change the heights of the polygon). I was looking to construct a function using the [contains] method but am lost on how to write a function to see if every point in the polyline, between that segment, fits inside of the polygon. Any help/ideas would be greatly appreciated. Thanks in advance


title: Problem in VersionEdit4.Reconcile4 method

text: I had an standalone application in which i am getting all the available versions in the database. Based on user filter he selects multiple versions and reconcile it.                if (pVersionEdit4.Reconcile4(pVesrion.VersionInfo.Parent.VersionName, true, false, false, false) == true)      I am iterating user selected versions and generating conflict report. But the previous one is **Row level** conflict check. But i want to check **column level**. For that my code changed to               if (pVersionEdit4.Reconcile4(pVesrion.VersionInfo.Parent.VersionName, true, false, false, true) == true)      But i am getting exception message **Failed Security**. Is any where i am missing.


title: Map Basic. Where-clause in Select Statement

text: what is wrong with this statement? I recon the where-clause is the problem.               Select * from "BEDIRECT_POT_" + vo_num, PLZ5_2012_bereinigt      where "BEDIRECT_POT_" + vo_num.PLZ like PLZ5_2012_bereinigt.PLZ      order by PLZ      into "BEDIRECT_POT_" + vo_num + "_PLZ"      This section seems to be wrong: "BEDIRECT_POT_" + vo_num.PLZ Thanks for your help. a.


title: How to crop an image into more than 80 small sections automatically in R

text: I'd like to crop an image into 15 sections in R. The image projection is Geog lat/long, extending 87 to 120 degree in East and 45 to 55 degree in North. I want to separate it by every 2 degree. So I have 17 degree in X coordinate and 5 degree in Y coordinate. 17 x 5= 85 sections. And then every section has own names such as crop1, 2, 3 etc. Can I do this automatically in R? If so, how to do this? Thanks a lot


title: Without Proxy page can't use IE8 and IE9 for Editing and Buffer

text: The server where my editing application is hosted does not allow a proxy page because of outgoing calls restriction. This makes it difficult for IE8 and IE9 since using CORS is not an option as there is limited support for these browsers. The limitation is that they do not support POST when using XDomainRequest and also there is a URL length limit. What is the best practice to POST Cross Origin messages from these browsers. Even if I write custom code I would like to leverage the heavy lifting done by esri's tasks like the ApplyEdits event handlers.


title: Installing Qgis2 on openSUSE 13.1

text: I need to install Qgis2 on my openSUSE 13.1 box but its standard repo (Application:Geo) on software.opensuse.org reports a build error. Has anyone already found a way to get it working? martin


title: data driven pages format

text: I have an issue, but I don't even know the terms needed to search for an answer. I'm exporting the pages to png files. The pages are divided by lots. Each page has a label with the lot number. In the case of the sample included, 18-385 is the lot being exported. Pieces of several other neighboring lots are also printed. What I would like to do if possible is format the label of the lot being exported different from the neighboring lots. The first pic is what I have so far, and the second pic is what I would like to do. Please take this into consideration: I am very new to ArcGIS.   I'm doing the exports in the Geoprocessing->Python box. I don't know yet how to get to a full IDE where I can save code.   I'm using ArcGIS 10.1 If it helps, here's the code I'm using.               mxd = arcpy.mapping.MapDocument("CURRENT")     ddp = mxd.dataDrivenPages     basePath = "C:\\Downloads\\Dev\\ArcGIS\\Fulton\\"      for pageNum in ddp.selectedPages:        mxd.dataDrivenPages.currentPageID = pageNum        distNum = mxd.dataDrivenPages.pageRow.DistNum        savePath = basePath + distNum        fileName = savePath + "\\" + mxd.dataDrivenPages.pageRow.DistPage + ".png"        print fileName        if not arcpy.os.path.exists(savePath):           arcpy.CreateFolder_management(baseDir,dirName)        arcpy.mapping.ExportToPNG(mxd, fileName,resolution=300)     del mxd      Thanks for any help. ![What I have](http://i.stack.imgur.com/YhLh2.png) ![What I want](http://i.stack.imgur.com/iJzHN.png)


title: GetFeatureInfo for gxp_olsource in OGS GeoExplorer

text: We have an OpenGeo Suite GeoExplorer 3.0.1 and want to make GetFeatureInfo- requests, via the identify button, for WMS-layers of type gxp_olsource but it doesn't work, I get no response at all and I don't get any errors in the browser debug console. It works great to make a GetFeatureInfo-request for WMS-layers with sources of type gxp_wmscsource or gxp_wmssource. Is it supposed to be working in the original version of GeoExplorer or do I have to make adjustments to the code? Layer configuration: ... sources: { ... ol: { ptype: "gx_olsource" } }, map: { ... { group: "groupfisketsgeografier", visibility: false, queryable: true, source: "ol", type: "OpenLayers.Layer.WMS", args: ["Fredningsomraden", "http://geodata.havochvatten.se/geoservices/hav-fisketsgeografier/ows", {layers: "hav-fisketsgeografier:havet-fredningsomraden", transparent: true, isBaseLayer: false, tiled: true} ] }, } ...


title: How to exclude the "More" tab where your operational layers are stored in the Flexviewer 2.2

text: I am using Robert's awesome TOC widget. I want to exclude the "more" tab where your operational layers show up, because they are listed in the TOC widget. Can I do this in the config.xml?


title: Finding DEM/HIllshade sources California

text: Creating watershed maps for schools. I want to color ramp the elevation to show low and high points. What kind of files do I need and where to find them?


title: ArcView 3.x Avenue Bitmaps (Tabs?) vs. ArcView 10 Python Cursors

text: _Note: While this question has an answer, any further tips for optimizing a cursor process would be greatly appreciated. I will be monitoring for any updates._ Currently, my boss (who works in Avenue) and I (working in Python) are both attempting to solve the same problem. Rather, we have both solved it, but the speed at which our solutions operate are... disjointed, to say the least. What his script processes in 2 hours can take mine up to 6. The only real difference in syntax and implementation in logic comes from 3.x's Bitmaps and 10.x's Cursors. We both: 1) Store values from Table 1.   2) Use those values to query a row in Table 2.   3) Store values from Table 2 for insertion into Table 3 as a new row. In both scripts, these processes are completed in two nested loops. Before I start digging into the wonderful world of code optimization, is this an expected occurrence when comparing Avenue script performance to Python? This is not the first time his scripts have greatly outperformed mine in terms of operation time, so I would like to know if there is something I should be aware of before I crucify myself for horrid scripting. Here's my script sans extraneous bits:               import arcpy     import time     import sys     import os          def recordfindcopy(inFile,query,outFile):         findRecord = arcpy.SearchCursor(inFile,query)         for record in findRecord:             copyRecord = arcpy.InsertCursor(outData) # <--- D'oh! (See answer)             field = record.FIELD             copy = copyRecord.newRow()             copy.FIELD = field             copyRecord.insertRow(copy)          StreetsFileList = [r"Path",                      r"Path"]          for sfile in StreetsFileList:         inStreets = sfile         inTable = r"Path"         outData = r"Path"         fsaEntry = arcpy.SearchCursor(inTable)         for row in fsaEntry:             id = row.ID             sQuery = "ID = %s " % (str(id))             recordfindcopy(inStreets,sQuery,outData)      **EDIT** : Given some of the comments so far, I wonder if there might be a better way to do this via joins, although I am dubious given the brobdingnagian (word of the day!) size of the tables. The heart of the processing is to append information from one table to any matching records in a second table and create a third table only containing the important fields. I wanted to try this using SDE, but that appears to not be an available option. Thoughts? I apologize if my questions are always so _involved_ , but I'm trying to get to the bottom of a long-standing annoyance. **Answered** : Jakub's simple suggestion alone decreased the processing time from **30** seconds per 500 records to **3** seconds per 500 records. Re- initiating the insert cursor on every insert slowed things down considerably (obviously). While this may not be the most optimization one can do for this process when put up against ArcView 3.x's speed, it is enough for my purposes at this time. Further suggestions are very welcome!


title: Standard SQL with ArcGIS

text: Is it possible to use Python to run actual standard SQL on an enterprise geo- database? Think about batch deleting, joints, etc moving beyond just the cursor objects in Arcpy. Note that one can shortcut ArcGIS and directly work with the back end database, but that is not the purpose of the questions here (and it I'd not recommended if you have n advanced geodatabase, like a versioned one)


title: NHD Flowline legend colors?

text: Does Anyone have legend colors for the NHD Flowline that USGS sets colors for theirs ? I'd like to get a copy of their color legend. They can be from their website, brouchure, lyr files or anything that I can set it up easily to match their colors. Thanks.


title: Create vector grid from canvas using script

text: I've modified "Hex grid from layer bounds" script to my need. Just added `cellsizex` and `cellsizey`. Changed the hexagonal - 3 to "1" to make it as a square grid. Is it possible to include an option to select extent from canvas? Instead of getting it from layer extent, I want to manually select the extents from the canvas by drawing a square or rectangle. Edit: I am new to scripting and though I can understand what a line of code does, I find it difficult to implement it by myself. For eg, with the help of the above said code, I was trying to get the extents from the canvas through a square or rectangle drawing. What I did was, Create new script -> just started to type (after the initial declaration or whatever) `## get extent`, `##select extent`...things like that. The problem is, i don't know which will be the right method, what is for what? Is script editor's code different from python's? (Don't know python either). What do I have to do to build small scripting tools by myself? If anyone can shed light as to where to start, it'll be immensely helpful.


title: Where can I find shapefiles for New Zealand Cities and Suburbs/Neighbourhoods?

text: I'm trying to find some shapefiles (administrative boundaries) of New Zealand CITIES and SUBURBS/NEIGHBOURHOODS. I've found some boundaries for the Country, Regions, Districts (not useful for me) .. but not cities and suburbs. Can anyone lead me in the right direction, please?


title: how to get ALL FIELD NAMES from table? (using arcpy.ListFields)

text: I have a table: ![enter image description here](http://i.stack.imgur.com/wLrUD.jpg) , Script at ModelBuilder's "calculate values": Expression:               a("%path%")      Code BLOCK:               def a(path):       fields = arcpy.ListFields(path,"FID_*")       for field in fields:         return field.name      Gives me only ONE field name!


title: Displaying text from dxf in Openlayers

text: I managed to import .dxf into Postgis table, and now have, among others, rows with Point geometries and values in text column (which corespodent with Label values from ACAD .dxf). What is the best way to display these geometries in Openlayers so they look the same as in any other dxf viewer (format, size, color)? Fonts, sizes and colors are imported into postgis using sql "select *,ogr_style from entities". The rest of the geometries (lines, polygons) are correctly displayed using HTTP protocol and GeoJSON format, but with labels missing. Also, is there a way to import Text from .dxf into Postgis as set of lines/polygons? That would also solve the problem.


title: Looking for qml specification

text: Where I can find qml specification which is used for store symbology information in QGIS? Try to googe it - no result.


title: ECW/MrSid plugin for GDAL from Kyngchaos not working on Mountain Lion?

text: Subject says it all. I have installed all of Kyngchaos packages, but               gdalinfo --formats      does not list ECW or MrSid formats. Could someone else with Mountain Lion try if they can get _any_ JP2 driver installed into GDAL?


title: Find Duplicate Polygons Algorithm

text: I have a group of polygons. My objective is to find the duplicate polygons from these group of polygons. Using NetTopologySuite(A c# port of JTS), it is possible to compare two geometires and check if they are equal. The brute-force method(checking each polygon against every other polygon) is the only idea that comes to my mind but it is not usable if there are large number of polygons. Are there any algorithm's which would better the brute-force in this case? I found this script which seems to contain ideas to what i am looking for(the comments indicate a divide and conquer approach..but not much details in the comments). But I have to admit i cannot get anything out of it :) .. It was done for ArcView, which I am not familiar with. **Note:** I am not looking for a PostGIS/database solution as in this question. I am looking for something that can be integrated to AutoCAD, Quantum GIS or such desktop GIS products through customization(c#,c++, python etc..)


title: Generating tiles for osm/esri shapefiles

text: I've generated tiles for OpenStreetMap .OSM files using Maperitive and quickly consumed them with leaflet (with no tile cache yet). I'm now looking for the right and simple way to do that also with the shapefiles given by OSM. (Tried TileMill but the generated mbtiles file needs further processing)


title: Drawing a rectangle that will always be on certain area of map

text: Using OpenLayers is it possible to draw a rectange that will always be on one portion of map without regards to zooming and panning of map?


title: How To Highlight Roads?

text: I wish to change the colors of various roads on a map and save the changes to either a "personal map" or a "map layer." The personal map could look something like this, except it would need to have the "snap to roads" feature, which is currently bugged on google mymaps. The map layer option seems a lot more complicated, but would be ideal if combined with a wiki function as the entire point of this exercise is to color roads according to a cyclist's perspective: according to speed and volume of traffic. It would be ideal if the road itself could be selected and highlighted so that any changes to the road would automatically be incorporated into the color scheme--instead of having colors going where the road no longer goes. I am not necessarily asking for a detailed step-by-step walkthrough of a specific possibility (though that would be nice), but rather a general overview of what my options are and hopefully some idea as to which options are most promising. I have been researching various mapping systems for several hours and feel as if I'm re-inventing the wheel. Please suggest a method for highlighting roads.


title: How do i popup a context menu with custom command items in a Select tool OnMouseDown (Right Click) with VB .Net?

text: This should be simple... I am converting VBA customization to VB .NET in VS 2008 Express Edition. In VBA i have a custom selection tool. Aside from just selecting features, when OnMouseDown receives Button = 2 parameter (right click) I popup a context menu with some macro items that call various other procedures. I need some help with making this work in .Net... I used the ESRI .NET sample project to get the feel for this and i am sure i can re-create the select tool but i am having problems with the context menu. I would also like to know how to call procedures when user clicks the context menu item.               Private Sub UIToolControl1_MouseDown(ByVal button As Long, _         ByVal shift As Long, ByVal x As Long, ByVal y As Long)          Dim pLayer As ILayer     Dim pContextMenu As ICommandBar     Dim pAppPosition As IWindowPosition     Dim pPoint As POINTAPI          If pMxDoc.SelectedLayer Is Nothing Then     MsgBox "Select a layer in the TOC!": GoTo ep     End If          Select Case button     Case 1     Set pLayer = pMxDoc.SelectedLayer     If pLayer Is Nothing Then     MsgBox "You must select 1 feature layer in the TOC!": GoTo ep     End If          If Not TypeOf pLayer Is IGeoFeatureLayer Then     MsgBox "You must select 1 feature layer in the TOC!": GoTo ep     End If          Set m_pPoint = _         pActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(x, y)     m_bIsMouseDown = True     Case 2     'on right click = popup     Set pContextMenu = CommandBars.Create("Go2Grid", esriCmdBarTypeShortcutMenu)          Set pAppPosition = Application          ' Add 3 built in commands to the new context menu using the     ' built-in ArcID module.          pContextMenu.CreateMacroItem "New Feature Class", 8, _         "Normal.Module1.CreateNewFeatureClass"     pContextMenu.CreateMacroItem "Delete Layer", 13, _         "Normal.Module1.DeleteFeatureClass"     pContextMenu.CreateMacroItem "Copy To...", 20, _         "Normal.ThisDocument.CopyFeatures"     pContextMenu.CreateMacroItem "Draws Around...", 26, _         "Normal.ThisDocument.DrawsAroundFeatures"     pContextMenu.CreateMacroItem "Clear Tool", 10, _         "Normal.ThisDocument.ClearCurrentTool"     pContextMenu.CreateMacroItem "Delete Features", 13, _         "Normal.Module1.DeleteFeatures"          ' Popup the menu.          pPoint.x = x     pPoint.y = y     ClientToScreen pActiveView.ScreenDisplay.hwnd, pPoint          pContextMenu.Popup pPoint.x, pPoint.y     Case Else     GoTo ep     End Select          End Sub


title: QGis: scripting the map rendering (multiple similar maps from one shapefile)

text: I would like to ask for advice on creating a **Python script** that will provide the opportunity to **render a plenty of similar maps** (with the same extent, projection, layout, color scheme and composition, but with differing titles and data used for rendering). Say we have a single shapefile or Spatialite DB with a huge attribute table inherited from csv file like this: ![enter image description here](http://i.stack.imgur.com/9fx5I.png) I want to turn them into a set of choropleth maps like this: ![enter image description here](http://i.stack.imgur.com/VeexB.png) So I guess that I need a script that will just:   1. Iterate through specified range of columns in the attribute table;   2. Call some (yet not written?) rendering function that will produce the choropleth map with a legend and title;   3. Dump the output to the vector formatted-file (PDF, EPS, etc.).  I'm going to use QGis for that purpose and would be glad if anyone could outline the main ideas and key QGis API functions that I can rely on at first. Thank you.


title: Export table structure information to Excel QGIS or MapInfo

text: I want to create a data dictionary of the tables of our primary layers. So I want a way to export the table structure to Excel. I don't need the actual data, I just want to know what columns are in each table and the data types for each column. Our layers are in MapInfo, and I also have access to QGIS. I can use the MapInfo Table Manager to export to a .met file, but when I copen that with Excel, the information is imported into one column - then requiring manual processing.


title: Modifying the Score text in the Locate Widget of Flexviewer 2.4

text: I want to add text where it says "score" to make it more clear that this is an address matching score. ![enter image description here](http://i.stack.imgur.com/nRgu7.jpg) I am using the LocateWidget_US.xml and do not think I can do it here. I think the place to do it may be in the LocateWidget.mxml in this line????                locateResult.content = getDefaultString("resultScoreText") + " " +       (addrCandidate.score % 1 == 0 ? addrCandidate.score : addrCandidate.score.toFixed(1));


title: Merging US counties in QGIS?

text: I'm doing research on US cotton price fluctuations from 1877-1930. Here's the problem: my data is county-level, but county borders changed from 1877 to 1930! To solve this, I plan to "aggregate up" to the smallest consistent unit (eg. if a county border splits, use the un-split original county as the unit of observation). How can I do this on QGIS? I am using data from the Historical County Project here: http://publications.newberry.org/ahcbp/ Any help would be greatly appreciated.


title: Is there a way to get a good human readable error from a COM error in arcEngine?

text: I'm trying to perform the following operation:                IPropertySet addressProperties = reverseGeocoding.ReverseGeocode(point, false);      And I get this error: Exception from HRESULT: 0x80040217 Is there some function to get a human readable form of this error?


title: Add toolbar - error "no object added"

text: I'm trying to add Hawths Tools to by Toolbar selection. I'm running ArcMap/ArcInfo v.9.3.1. Hawths Tools is not supported in v.10. I've dowloaded the zip file, extracted the files, adn ran teh .exe file. When I go to "customize" Tools, add from file, and select under Program files "HawthsTools3.dll", an "Added Objects..." menu pops up that reads: No new objects added. If you then open ArcMap, select View, then Toolbars, Hawths Tools does not show up. Can anyone provide assistance?


title: Want to convert Lat/Long to UTM

text: I want to convert Lat/Long to UTM. When I search on net I am getting completely different formulas. e.g. Lat/Long, 21.0909,79.067635 How to convert it to UTM? Thanks


title: How to install Apache Web Server from OSGeo4W?

text: Good morning, I'm trying to create a QGIS web client on my Windows 8 But with OSGEO4W not appear the option to make Apache installs Web Server, as of this writing in several tutorials and documentations. A print screen of my Windows 8 ![enter image description here](http://i.stack.imgur.com/CVUea.jpg) Should this well, this print belongs to a tutorial found on the web ![enter image description here](http://i.stack.imgur.com/bnolO.jpg)


title: Is there a way to have a choropleth map in Excel reading data from a sheet?

text: Is it possible to have a solution similar to ESRI maps for office in which I want basically a geojson file or any file for that matters to render a choropleth map based on data in a sheet? Tried already this tutorial http://www.clearlyandsimply.com/clearly_and_simply/2012/09/faster-choropleth- maps-with-microsoft-excel.html but I have a to do this in a way that can be dynamic. I was thinking exporting to geojson and loading on LeafletJS and be able to display it within Excel, is this even possible?


title: Mapserver and SLD Files

text: I have a set of data layers that define a map from small to large scales. These layers are available through mapserver and are grouped in the mapfile so:-               LAYER       GROUP large_scale       NAME roads       ..     END     LAYER        GROUP large_scale       NAME water      ...     END     LAYER        GROUP medium_scale       NAME roads      ....     END     ....      I have been given SLD files for these layers but they break down into the NAME level not the group level. So I **large_scale_roads.sld** **large_scale_water.sld** **medium_scale_roads.sld** etc.. As I understand it mapserver cannot use an SLD file from the mapfile itself, this has to come from the request. I am hoping to publish this all through a WMS with the layer name as either **large_scale** or **medium_scale** and using the SLD files rather than the mapserver file. Can I do this given multiple SLD files, or do I have to combine them all into a single file which is going to be huge to cover all of these layers? Or is there a better approach? I do not want to re-do all of the styling in the SLD files and am using mapserver behind a mapcache tile cache.


title: How do I stop null values from being assigned extreme color values?

text: I have a raster surface that I kriged in Surfer and imported to Grass GIS but when I view it in QGIS, the null values outside of the interpolated surface take an extreme value and wash out the resolution in the rest of my surface. Is there an easy way to stop the null values from being assigned off the scale values?


title: Autoincrement field when New Feature added

text: I know that this is probably a very basic programming but I am brand new to Python. I have a layer in a geodatabase that I need a unique ID field (for joining). I can't use the ObjectID because that can change. I can't use a GlobalID because the data will be used by databases that can't handle the "{" brackets or dashes. I have already populated the field. When I add a new feature I want to have the value of the max existing value + 1 for the new feature. Thanks to Jason here is what I have so far-               import arcpy     gdb_full_path = r'D:\GIS\ZA.gdb\Current'     feat = 'ZA_Merge'     arcpy.env.workspace = gdb_full_path     with arcpy.da.SearchCursor(feat,["ANID"]) as cursor:          for row in cursor:              unique_id_counter = row.count     with arcpy.da.InsertCursor(feat,["ANID"]) as cursor:          cursor.insertRow((unique_id_counter + 1))      Here is the error that I get - Runtime error Traceback (most recent call last): File "", line 9, in TypeError: unsupported operand type(s) for +: 'builtin_function_or_method' and 'int'


title: How to represent trend over time?

text: I have a time series of rasters where the cells of each raster represents some value at a certain time. I'd like to generate a map that highlights the magnitude of positive or negative trends in this value over time for each cell. My somewhat naive approach is to fit a simple linear regression (X=time and Y=value) to each cell and output the array of slopes to a raster (as per example images below). This can be filtered by only exporting significant results. ![time series from four individual cells](http://i.stack.imgur.com/XEPTu.jpg) ![Slopes from linear regression](http://i.stack.imgur.com/KIbyS.jpg) How else might I represent trend over time in a raster timeseries? NB: I'm interested in general techniques not software specific instructions.


title: Python GDAL: Save array as raster with projection from other file

text: I have an array of data, and for each datapoint I know the latitude and longitude. I'd like to save it as a GTiff with the same projection as other rasters I have. This is what I've tried so far, but no luck.               import numpy as np     import gdal     from gdalconst import *     from osgeo import osr          def GetGeoInfo(FileName):         SourceDS = gdal.Open(FileName, GA_ReadOnly)         GeoT = SourceDS.GetGeoTransform()         Projection = osr.SpatialReference()         Projection.ImportFromWkt(SourceDS.GetProjectionRef())         return GeoT, Projection          def CreateGeoTiff(Name, Array, driver,                        xsize, ysize, GeoT, Projection):         DataType = gdal.GDT_Float32         NewFileName = Name+'.tif'         # Set up the dataset         DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType )                 # the '1' is for band 1.         DataSet.SetGeoTransform(GeoT)         DataSet.SetProjection( Projection.ExportToWkt() )         # Write the array         DataSet.GetRasterBand(1).WriteArray( Array )         return NewFileName          def ReprojectCoords(x, y,src_srs,tgt_srs):         trans_coords=[]         transform = osr.CoordinateTransformation( src_srs, tgt_srs)         x,y,z = transform.TransformPoint(x, y)         return x, y          # Some Data     Data = np.random.rand(5,6)     Lats = np.array([-5.5, -5.0, -4.5, -4.0, -3.5])     Lons = np.array([135.0, 135.5, 136.0, 136.5, 137.0, 137.5])          # A raster file that exists in the same approximate aregion.     RASTER_FN = 'some_raster.tif'          # Open the raster file and get the projection, that's the     # projection I'd like my new raster to have, it's 'projected',     # i.e. x, y values are numbers of pixels.     GeoT, TargetProjection, DataType = GetGeoInfo(RASTER_FN)     # Meanwhile my raster is currently in geographic coordinates.     SourceProjection = TargetProjection.CloneGeogCS()          # Get the corner coordinates of my array     LatSize, LonSize = len(Lats), len(Lons)     LatLow, LatHigh = Lats[0], Lats[-1]     LonLow, LonHigh = Lons[0], Lons[-1]     # Reproject the corner coordinates from geographic     # to projected...     TopLeft = ReprojectCoords(LonLow, LatHigh, SourceProjection, TargetProjection)     BottomLeft = ReprojectCoords(LonLow, LatLow, SourceProjection, TargetProjection)     TopRight = ReprojectCoords(LonHigh, LatHigh, SourceProjection, TargetProjection)     # And define my Geotransform     GeoTNew = [TopLeft[0],  (TopLeft[0]-TopRight[0])/(LonSize-1), 0,                TopLeft[1], 0, (TopLeft[1]-BottomLeft[1])/(LatSize-1)]          # I want a GTiff     driver = gdal.GetDriverByName('GTiff')     # Create the new file...     NewFileName = CreateGeoTiff('Output', Data, driver, LatSize, LonSize, GeoTNew, TargetProjection)      But this results in the following error message:               File "TES2GtifBBB.py", line 25, in CreateGeoTiff       DataSet.GetRasterBand(1).WriteArray( Array )     File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal.py", line 1082, in WriteArray       return gdalnumeric.BandWriteArray( self, array, xoff, yoff )     File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal_array.py", line 256, in BandWriteArray       raise ValueError("array larger than output file, or offset off edge")     ValueError: array larger than output file, or offset off edge


title: Will Google Earth see my kml - on different versions?

text: Is there any difference in the kml script to allow different versions of Google Earth to see it ? If I write               <?xml version="1.0" encoding="UTF-8"?>     <kml xmlns="http://earth.google.com/kml/2.1">     <Document>          </Document></kml>      How would the kml file content change if using a different Google Earth version ? (I would like to be sure that the kml files I create are visible on currently downloadable Google Earth versions, e.g. 5 to 7.03) I plan to use simplekml (pyPI) to generate the files. They will use GroundOverlay, NetworkLink, Region... (as created by simplekml, which I trust to do a better job than my own brain) I tried to find any information on version changes (for Google Earth) - but have been unsuccessful. Thank you for help.


title: QGIS 1.8.0 on Mac - how to cut/paste/delete?

text: I have a line vector layer selected, editing is on, I select one vector (a dozen points) with the node select tool so the vertices are red squares.   1. the cut icon is greyed out, command/X does nothing, delete part isn't greyed out but does nothing; so I can't cut/copy/delete it   2. I select a vertex so its square turns to blue, then hit the delete key - nothing happens. Remember, a PC delete key is forward delete, back delete is the backspace key; a Mac has no forward delete, the delete key is a back delete; is this the problem or is it related to 1.?   3. I can't select more than one vector by holding down the control key (or shift or command - I've tried everything I can think of) - when I click on the second vector the first is deselected. Again, is this a Mac-only problem?


title: How to fix erros loading result layers from SAGA functions in Sextante?

text: Using QGIS 1.8.0 and SAGA When trying to use the IDW function from Sextante SAGA 2.1.0 64bit, I get the following error: `Could not load layer: C:/test.tif`. This error seems mysterious to me, any ideas? ![enter image description here](http://i.stack.imgur.com/JzH9v.png) **Sextante error logs** With vector layer loaded as a PostGIS layer: `Could not load model watersheds.model ALGORITHM:ftools:export/addgeometrycolumns` `sextante.runalg("saga:inversedistanceweighted","dbname='mydb' host=www.mywebsite.com port=5432 user='postgres' password='mypassword' sslmode=disable key='id' estimatedmetadata=true srid=2163 type=POINT table="analysis"."scores" (geom) sql=","score",0,0,2,1,0,100,0,0,10,"2297980.75,2325375.5,108272.492188,134982.28125",0,None)` With vector layer loaded as a shp: `sextante.runalg("saga:inversedistanceweighted","C:/test.shp","score",0,0,2,1,0,100,0,0,10,"2297980.80156,2325375.34865,108272.496186,135921.444199",0,None)` **System** I am using QGIS 1.8.0 with SAGA 2.1 b3 64bit on Windows 8 64bit. **No problems running it directly in SAGA without SEXTANTE** I cant figure out how to import the raster `.sgrd` file back into QGIS... _Some how the points loaded into SAGA are skewed at an angle.. thats for another question_ ![enter image description here](http://i.stack.imgur.com/TSGVp.png) * * * ## Update Updated to QGIS v.1.9.90. When activating SAGA under SEXTANTE, I get the following error. I get the same error when I choose any of SAGA's function under SEXTANTE's toolbox.               An error has occured while executing Python code:          Traceback (most recent call last):       File "C:\Users\0_o/.qgis2/python/plugins\sextante\gui\SextanteToolbox.py", line 138, in executeAlgorithm         dlg = ParametersDialog(alg)       File "C:\Users\0_o/.qgis2/python/plugins\sextante\gui\ParametersDialog.py", line 43, in __init__         self.paramTable = ParametersPanel(self, alg)       File "C:\Users\0_o/.qgis2/python/plugins\sextante\gui\ParametersPanel.py", line 74, in __init__         self.initGUI()       File "C:\Users\0_o/.qgis2/python/plugins\sextante\gui\ParametersPanel.py", line 116, in initGUI         widget = self.getWidgetFromParameter(param)       File "C:\Users\0_o/.qgis2/python/plugins\sextante\gui\ParametersPanel.py", line 257, in getWidgetFromParameter         item.addItems(self.getFields(layers[0], param.datatype))       File "C:\Users\0_o/.qgis2/python/plugins\sextante\gui\ParametersPanel.py", line 325, in getFields         for idx, field in fieldMap.iteritems():     AttributeError: 'QgsFields' object has no attribute 'iteritems'          Python version:     2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]               QGIS version:     1.9.0-Master Master, 404d89f          Python path: ['C:\\Users\\0_o/.qgis2/python/plugins\\sextante', 'C:/OSGeo4W/apps/qgis-dev/./python', 'C:\\Users\\0_o/.qgis2/python', 'C:\\Users\\0_o/.qgis2/python/plugins', 'C:/OSGeo4W/apps/qgis-dev/./python/plugins', 'C:\\OSGeo4W\\bin\\python27.zip', 'C:\\OSGeo4W\\apps\\Python27\\DLLs', 'C:\\OSGeo4W\\apps\\Python27\\lib', 'C:\\OSGeo4W\\apps\\Python27\\lib\\plat-win', 'C:\\OSGeo4W\\apps\\Python27\\lib\\lib-tk', 'C:\\OSGeo4W\\bin', 'C:\\OSGeo4W\\apps\\Python27', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\OSGeo4W\\apps\\qgis-dev\\python\\plugins\\fTools\\tools', 'C:/Users/0_o/Desktop/test/Scores']      * * * ## Update #2 Tried an older version of QGIS 1.9.0 298f0c0 **SAGA execution console output**               SAGA execution console output          C:\>set SAGA=C:\Program Files\saga          C:\>set SAGA_MLB=C:\Program Files\saga\modules          C:\>PATH=PATH;C:\Program Files\saga;C:\Program Files\saga\modules          C:\>saga_cmd grid_gridding "Inverse Distance Weighted" -SHAPES "C:/Users/0_o/Desktop/mytest/Scores/MA Full test/MA Scores.shp" -FIELD "score" -TARGET 0 -WEIGHTING 0 -POWER 2 -BANDWIDTH 1 -RANGE 0 -RADIUS 100 -MODE 0 -POINTS 0 -NPOINTS 10 -USER_XMIN 2135113.26163 -USER_XMAX 2341329.19555 -USER_YMIN 88475.4178879 -USER_YMAX 185492.363448 -USER_SIZE 80 -USER_GRID "C:\Users\0_o\AppData\Local\Temp\sextante\MA IDW.tif.sgrd"     Unknown option 'POWER'     Usage: saga_cmd -SHAPES <str> [-FIELD <str>] [-TARGET <str>] [-SEARCH_RANGE <str>] [-SEARCH_RADIUS <str>] [-SEARCH_POINTS_ALL <str>] [-SEARCH_POINTS_MAX <num>] [-SEARCH_DIRECTION <str>] [-WEIGHTING <str>] [-WEIGHT_POWER <str>] [-WEIGHT_BANDWIDTH <str>] [-USER_XMIN <str>] [-USER_XMAX <str>] [-USER_YMIN <str>] [-USER_YMAX <str>] [-USER_SIZE <str>] [-USER_GRID <str>] [-GRID_GRID <str>]     -SHAPES:<str>               Points     Shapes (input)     -FIELD:<str>                Attribute     Table field     -TARGET:<str>               Target Grid     Choice     Available Choices:     [0] user defined     [1] grid     -SEARCH_RANGE:<str>         Search Range     Choice     Available Choices:     [0] local     [1] global     -SEARCH_RADIUS:<str>        Maximum Search Distance     Floating point     Minimum: 0.000000     -SEARCH_POINTS_ALL:<str>    Number of Points     Choice     Available Choices:     [0] maximum number of nearest points     [1] all points within search distance     -SEARCH_POINTS_MAX:<num>    Maximum Number of Points     Integer     Minimum: 1.000000     -SEARCH_DIRECTION:<str>     Search Direction     Choice     Available Choices:     [0] all directions     [1] quadrants     -WEIGHTING:<str>            Distance Weighting     Choice     Available Choices:     [0] inverse distance to a power     [1] linearly decreasing within search radius     [2] exponential weighting scheme     [3] gaussian weighting scheme     -WEIGHT_POWER:<str>         Power     Floating point     -WEIGHT_BANDWIDTH:<str>     Bandwidth     Floating point     Minimum: 0.000000     -USER_XMIN:<str>            Left     Floating point     -USER_XMAX:<str>            Right     Floating point     -USER_YMIN:<str>            Bottom     Floating point     -USER_YMAX:<str>            Top     Floating point     -USER_SIZE:<str>            Cellsize     Floating point     Minimum: 0.000000     -USER_GRID:<str>            Grid     Data Object (optional output)     -GRID_GRID:<str>            Grid     Grid (output)          C:\>saga_cmd io_gdal 1 -GRIDS "C:\Users\0_o\AppData\Local\Temp\sextante\MA IDW.tif.sgrd" -FORMAT 1 -TYPE 0 -FILE "C:/Users/0_o/Desktop/mytest/Scores/MA Full test/MA IDW.tif"          error: Grid file could not be opened.          error: input file [C:\Users\0_o\AppData\Local\Temp\sextante\MA IDW.tif.sgrd]          error: empty input list [GRIDS]     Usage: saga_cmd -GRIDS <str> [-FILE <str>] [-FORMAT <str>] [-TYPE <str>] [-OPTIONS <str>]     -GRIDS:<str>    Grid(s)     Grid list (input)     -FILE:<str>     File     File path     -FORMAT:<str>   Format     Choice     Available Choices:     [0] ERDAS JPEG2000 (SDK 3.x)     [1] Flexible Image Transport System     [2] Network Common Data Format     [3] Virtual Raster     [4] GeoTIFF     [5] National Imagery Transmission Format     [6] Erdas Imagine Images (.img)     [7] ELAS     [8] In Memory Raster     [9] MS Windows Device Independent Bitmap     [10] PCIDSK Database File     [11] ILWIS Raster Map     [12] SGI Image File Format 1.0     [13] Leveller heightfield     [14] Terragen heightfield     [15] USGS Astrogeology ISIS cube (Version 2)     [16] ERMapper .ers Labelled     [17] Raster Matrix Format     [18] Idrisi Raster A.1     [19] Intergraph Raster     [20] Golden Software Binary Grid (.grd)     [21] Portable Pixmap Format (netpbm)     [22] ENVI .hdr Labelled     [23] ESRI .hdr Labelled     [24] PCI .aux Labelled     [25] Vexcel MFF Raster     [26] Vexcel MFF2 (HKV) Raster     [27] VTP .bt (Binary Terrain) 1.3 Format     [28] Erdas .LAN/.GIS     [29] Image Data and Analysis     [30] NOAA Vertical Datum .GTX     [31] NTv2 Datum Grid Shift     [32] ARC Digitized Raster Graphics     [33] SAGA GIS Binary Grid (.sdat)     -TYPE:<str>     Data Type     Choice     Available Choices:     [0] match input data     [1] 8 bit unsigned integer     [2] 16 bit unsigned integer     [3] 16 bit signed integer     [4] 32 bit unsigned integer     [5] 32 bit signed integer     [6] 32 bit floating point     [7] 64 bit floating point     -OPTIONS:<str>  Creation Options     Text          C:\>exit      **SAGA execution commands**               SAGA execution commands     grid_gridding "Inverse Distance Weighted" -SHAPES "C:/Users/0_o/Desktop/mytest/Scores/MA Full test/MA Scores.shp" -FIELD "score" -TARGET 0 -WEIGHTING 0 -POWER 2 -BANDWIDTH 1 -RANGE 0 -RADIUS 100 -MODE 0 -POINTS 0 -NPOINTS 10 -USER_XMIN 2135113.26163 -USER_XMAX 2341329.19555 -USER_YMIN 88475.4178879 -USER_YMAX 185492.363448 -USER_SIZE 80 -USER_GRID "C:\Users\0_o\AppData\Local\Temp\sextante\MA IDW.tif.sgrd"     io_gdal 1 -GRIDS "C:\Users\0_o\AppData\Local\Temp\sextante\MA IDW.tif.sgrd" -FORMAT 1 -TYPE 0 -FILE "C:/Users/0_o/Desktop/mytest/Scores/MA Full test/MA IDW.tif"


title: Creating polygons for a choropleth map using PostGIS

text: I have a grid of points, each with a value to it as shown in the image below. The data was retrieved from PostGIS and rendered in Quantum GIS. The final result should be a polygons of a choropleth map that will be used to render on the browser. **Question:** How can PostGIS use the grid of weighted points to create the choropleth polygons? **Initial** ![enter image description here](http://i.stack.imgur.com/uKwcr.png) **Desired result** Polygons to be created in PostGIS, rendered in Leaflet. Classes group according to some rule like:               d > 1000 ? '#800026' :     d > 500  ? '#BD0026' :     d > 200  ? '#E31A1C' :     d > 100  ? '#FC4E2A' :     d > 50   ? '#FD8D3C' :     d > 20   ? '#FEB24C' :     d > 10   ? '#FED976' :      ![enter image description here](http://i.stack.imgur.com/Fh71J.png)


title: How to find the width and height of the html content dynamically to pass to esri map infowindow.resize?

text: There will be a html content stored in a variable with div id to pass into `infowindow`               map.infoWindow.resize(width,height);     map.infoWindow.setTitle("Callout");     map.infoWindow.setContent(calloutData );     map.infoWindow.show(event.screenPoint,map.getInfoWindowAnchor(event.screenPoint));      How can i pass the width and height to `infowindow.resize` dynamically? Is there any way the width and height of the html content can be calculated and pass without passing the value by us?


title: How to calculate Line Density in QGIS?

text: I have a map of fracture distributions over an area and I would like to calculate the density distribution of fractures that intersect a moving window with a specified radius. Is there a way to calculate this in QGIS, similar to the line density calculation in ArcGIS Spatial Analysis? Thanks


title: OSM To Postgis & WGS84 problems

text: I’m trying to import OpenStreetMap (OSM) data into postGIS it doesn’t appear to be working in that features do not appear to be in the correct place.. Using the QGIS OSM add in and a coastal shape file imported into PostGIS setting the SRS to 4326 on import (all my data uses WGS84 Longitude & Latitude) the two data sets are in the same space. Importing the OSM data set into postgis with the following code & adding it as a QGIS layer I would expect lines points etc of the OSM file & PostGIS layers to be adjacent if not on top of each other but they are not, I can’t actually see the coastal shapefile & the OSM PostGIS layer together without turning one layer off and selecting zoom to layer on the other. BTW I have attempted to use the -E 4326 switch with out -m but it failed So what am I doing wrong?               createDB osm     createlang plpgsql osm –U username      psql –U username –d osm  -f /path/to/postgis1.5/postgis.sql     psql –U username –d osm  -f /path/to/postgis1.5/spatial_ref_sys.sql     psql –U username –d osm  -f /path/to/osm2pgsql/900913.sql          osm2pgsql –S /path/to/osm2pgsql/default.style –U username –d osm  -l –c –m /path/to/ cambridgeshire.osm.bz2      I then tried following the instructions on http://spatialmounty.blogspot.co.uk/2010/12/transform-osm-data-in-postgis-to- wgs84.html to re-project the data to WGS84 Long Lat it failed at The first line with invalid SRID               SELECT AddGeometryColumn( 'planet_osm_line', 'geom_wgs84', 4326, 'LINESTRING', 2);      My database has the following tables               Geography_columns  view     Geometry_columns table     Planet_osm_line table     Planet_osm_point table     Planet_osm_polygon table     Planet_osm_roads table     Spatial_ref_sys table      Can somebody give me the steps / syntax i need to follow to accurately import OSM data into PostGIS so that it uses ESPG:4326? Any one know of any ground truth data I can use to test the process?


title: Multiple instance of leaflet on my page

text:
I need a few pointers really; I've got a series of tabs running down the left hand of my page and on each page I have a leaflet map. Each map is based on the choropleth example but I'm running into problems. Only one of the three maps enables me to do the custom info call to put the attributes into the info panel. I've quickly worked out that each map has to have its own 'map' ID e.g

     ...var map = L.map('map1').setView([37.8, -96], 4);...

     ...var map = L.map('map2').setView([37.8, -96], 4);...

     ...var map = L.map('map3').setView([37.8, -96], 4);...      so I'm guessing that I would need to make each instance of the info.onAdd function unique? I've played around do this but with no luck.               var info = L.control();     info.onAdd = function (map) {         this._div = L.DomUtil.create('div', 'info'); // create a div with a class "info"         this.update();         return this._div;     };     // method that we will use to update the control based on feature properties passed     info.update = function (props) {         this._div.innerHTML = '
US Population Density
' +  (props ?             '' + props.name + '
' + props.density + ' people / mi2'             : 'Hover over a state');     };     info.addTo(map);      Any ideas or pointers welcomed Cheers Si


title: How to query aggregate max, min, median latitude of polygons in PostGIS?

text: I am just starting to learn how to use PostGIS (1.5.4) and my SQL skills are essentially those of a competent beginner. Having perused the PostGIS documentation and searched these forums, I'm not finding the answers to what are surely some very simple questions. Given a single table loaded with polygons from shapefiles, how can I query the latitudinal maximum and minimum and the centroids of single and aggregate rows? I'll note that I have found the ST_Centroid() function, but haven't figured out how to interpret its output (CRS=WGS84).


title: Deploy new WPS at deegree without UI

text: How can I deploy a custom WPS on deegree without using the UI?


title: Best method to convert CSV file of lat/lon and depth data to bathymetry polygon?

text: I'm trying to create bathymetry polygons based on lat/long and depth data. I have tried several methods and the resulting polygons are clunky. Make XY Event Layer-->Copy Features--> Project--> Interpolate with IDW-->Make values Integeter with (Int tool) --> Raster to Polygon Dissolve Polygons by depth, so that all areas of equal depth are in same feature in the bathymetry feature class. This method results as you might imagine is a clunky, blocky polygon. I have tried several routes, including Contour, Create Filled Contour - I'm sure there must be a simpler way to create a bathymetry polygon?


title: QGIS I would like to get a single point from a polygon Centroid if within the polygon any point if not

text: I am trying to generate label points from my polygons. The Polygon->Centroid works great except that if the centroid happens to land in a hole it is still returned. I need to get a single point (Centroid if that point is within the polygon) a point within the polygon if the centroid happens to be outside of the polygon. Thank you,


title: ArcGIS 10 Add Field to Feature Class

text: In ArcGIS 10, is there a way to add a field to a feature class even if it currently contains data?


title: PostGIS: calculate width/length of rectangular surface

text: I have rectangular shaped surfaces (multipolylines) and I need to calculate the width and/or length.. they have a rectangular shape but they're not real rectangles. Is there a way to do this ?


title: How can I generalize "field name" variable?

text: Here is my code:               import arcpy     import math          arcpy.env.workspace = "C:\RG\Kings County\Test.gdb"     fcs = arcpy.ListFeatureClasses("NVF_*")          try:         for fc in fcs:             rows = arcpy.UpdateCursor(fc)             for row in rows:                 if row.Field_ID == 1:                     rows.deleteRow(row)             del rows, row     # This code adds up several fields like Ratio_NDVI, Ratio_Acres etc.             arcpy.AddField_management (fc, "Ratio_NDVI", "float")             arcpy.AddField_management (fc, "Ratio_Area", "float")             arcpy.AddField_management (fc, "Final", "integer")-      # This code calculates area in square meter using existing datum coordinates of the vector file             arcpy.AddField_management(fc, "NewArea", "DOUBLE")             geometryField = arcpy.Describe(fc).shapeFieldName             cursor = arcpy.UpdateCursor(fc)             for row in cursor:                 AreaValue = row.getValue(geometryField).area                 row.setValue("NewArea", AreaValue)                 cursor.updateRow(row)             del row, cursor       # This code converts square meter into acres             rows = arcpy.UpdateCursor(fc)             for row in rows:                 row.NewArea = (row.NewArea*0.00024711)                 rows.updateRow(row)             del rows, row        # This code calculates ratio of NDVI of child_fields to the average NDVI of parent field             rows = arcpy.UpdateCursor(fc)             for row in rows:                 if row.NDVI090606 == 0:                      row.Ratio_NDVI = (row.M090606/row.M090606)                  else:                     row.Ratio_NDVI = (row.M090606/row.NDVI090606)                  rows.updateRow(row)             del rows, row        # This code calculates the ratio of acres of child fields to their corresponding parent field             rows = arcpy.UpdateCursor(fc)             for row in rows:                 row.Ratio_Area = (row.NewArea/row.Acres)                 rows.updateRow(row)             del rows, row             print ("Completed calculating ratio of acres")        # This code finds that field that has either Ratio_NDVI higher than 1.5 to its parent #field or less than 0.5 to its parent field and Ratio_Acres higher than 0.05. Then #comments this field as "Change"             rows = arcpy.UpdateCursor(fc)             for row in rows:                 if (row.Ratio_Area > 0.05 and row.Ratio_NDVI > 1.5):                     row.Final = 1                     rows.updateRow(row)                 if (row.Ratio_Area > 0.05 and row.Ratio_NDVI < 0.5):                     row.Final = 1                     rows.updateRow(row)             del rows, row     except Exception as e:         print e      In this code, how can I generalize the field name variable (row.NDVI090906). In above, NDVI090906 is a field name, where this field name follows the pattern of the name of the file. For example, if the file name is 'NVF_090906' the corresponding field name is NDVI090906. How can I make this field name as a variable so that when it loops through, it automatically performs computation of field names of other feature classes.


title: What causes QGIS 2 Plugin Manager Error from showPluginManagerWhenReady?

text: there. when i open the plugin manager in qgis, i get this: Ocurrió un error durante la ejecución del siguiente código: pyplugin_installer.instance().showPluginManagerWhenReady() > A Google translation of the Spanish is: An error occurred during the > execution of the following code Traceback (most recent call last): File "", line 1, in NameError: name 'pyplugin_installer' is not defined any help?


title: How does one build a Python script that exports to several output files?

text: I'm trying to build a Python Script that will receive one feature layer that will be clipped by several adjacent grid squares (each in a different shapefile) and output several shapefiles. Now, I've built in the past a loop script that receives a list of files and does a single action on each file, like this for file in os.listdir(ws): if file.endswith(".shp") or file.endswith(".dxf") or file.startswith("str"): print file But I'm not succeeding in making a loop that will receive the input shapefile and clip features (shapefiles) and EXPORT to different shapefiles (I mean how do I change the output shapefile's name for every grid square ?). Can someone show how it is done? Thanks in advance. (Oh, and I'm using python 2.6 and Arcgis 9.3.1)


title: producing a bathymetric map in qgis with limited depths measurements

text: i am new to your community and hope you to find tips on how to work with the following project. i want to produce a bathemetric map (depths profile) of a water area. I do have aerial images and 3 cross sections with depth measurements. at first i would like to "connect" the depth lines to receive contour lines as used in mapping elevation and form there on estimate the depth at places where no measurements are taken. i do really have an idea on how to start as i only have done some remote sensing before. Please ask if i have not described it accurately enough. Does anyone have suggestions on how to tackle this? much thanks in advance!


title: Install plugins for qgis in centos

text: I followed the instructions here: How to install QGIS on Centos 6? If I run:               QgsVectorLayer * vpoly = new QgsVectorLayer("Polygon", "pointbuffer", "memory");     QgsVectorDataProvider * provider = vpoly->dataProvider();      Then provider is a NULL pointer Also, if I do:               QString myPluginsDir = "/usr/lib64/qgis";     QgsProviderRegistry * preg = QgsProviderRegistry::instance(myPluginsDir);     QString pluglist=preg->pluginList();     printf("plugins: %s\n",pluglist.toStdString().c_str());      Then it prints:               plugins: No data provider plugins are available. No vector layers can be loaded      Also, when I try to load a QgsVectorLayer with the following:               QgsVectorLayer * mypLayer = new QgsVectorLayer(myLayerPath, myLayerBaseName, myProviderName);     if (mypLayer->isValid()){         qDebug("Layer is valid");       } else {         qDebug("Layer is NOT valid");         return;       }      Then it says that the layer is NOT valid. This all works just fine in Ubuntu, but I can't get it to work in centos. I think maybe it's missing all the plugins. What am I missing? Please help.


title: Converting GPS heights to AHD heights

text: I was wondering if there is a way to convert heights that are logged from a normal bluetooth GPS to AHD heights.


In ArcGIS, I have only a point layer relating to a hydrography network's pour points. Now I want to write code in VBA that selects a point from this layer, then highlight the closest points to this point. One of this points must be downstream of the desired point and another point must be upstream. In this work ,I must measure the distance between two points on the network, not a direct distance between two points.Also both points are related to a layer. Also The program must be able to determine one of points in downstream and another points in uostream. I don't know How do I determine the closest points on the network? please help me again.


title: Which API should we use for printing a PageLayout?

text: We are building an standalone ArcGIS Engine application and we are going to add support for printing. There are several APIs for printing a PageLayout in ArcObjects, in ArcGIS 10 there are at least three different ways to do this:   * The new PrintAndExport class, in the documentation they seem to promote this API Conceptual help - Printing maps.   * It's also possible to print a layout using the PrintPageLayout method PageLayoutControl Sample: Printing with the PageLayoutControl   * But also Output method on IActiveView can be used. Can someone give recommendations on which API we should choose? Is there any known pitfalls in these APIs? Update: Found one pitfall during my research:   * When working with ArcGIS Server layers, only PrintAndExport draw the patches/swatches on the legend when printing the layout.


title: Upgrading Postgres database to Postgis 2.1.1

text: I executed the postgres 9.1.6 with Postgis 1.5.8 dump from a remote server on Linux pg_dump -h localhost -p 5432 -U postgres -Fc -b -v -f "/somepath/tumidb9.1geom1.5.8.backup" tumidb9 So I have a valid dump Can I just run the executable postgresql-9.3.0-rc1-linux-x64.run, select the Postgis 2.1.0 option with the check box that specifies prep for restore, create and intialize a new spatial database newdb in the 9.3 instance and then execute perl ../postgis_restore.pl "/somepath/tumidb9.1geom1.5.8.backup" | psql -h localhost -p 5432 -U postgres newdb 2> errors.txt on the 9.3 instance. or do I have to install postgresql-9.1.6-1-linux-x64.run, with the Postgis 2.1.1 option, run the postgis_restore.pl on the 9.1 instance, and then use pg_upgrade to get us to 9.3? thanks


title: How to get a list of features from Oracle via FME?

text: ![Added Image shows the workbench used: note in the lower left corner the concatenated result the FeatureReader is supposed to query](http://i.stack.imgur.com/TsE3P.jpg) I have a table in Oracle containing only the names of other tables of varying feature types, which I'd like to process (i.e.export to csv). this task is to be scheduled via FME Server. the feature types to work on are only known at runtime. Therefore a test like Mapperz suggested would not work in this situation. There are several examples how to do this based on names in text files or using cx_oracle (see http://fmepedia.safe.com/articles/How_To/Pass-a-List-of- Tables-or-Layers-to-an-FME-Reader-using-Python-Scripted-Parameters). My situation is that cx_oracle isn't availabe and a text file won't do it as entries will change frequently (but the latter approach via scripted parameter is actually working fine). The link above states in the beginning that: (...) For many situations, using an FME reader to read the table list and then the FeatureReader transformer to read the data may be a better approach than the one described below. (..) I tried to follow this approach, which is also described by Mark Irland FME Features to read parameter, set up the reader for the table with the names of the other tables and connected it to a FeatureReader. It is configured to query the feature types specified in the attribute below (which is set to a concatenated attribute as described at "Feature Type List" in http://blog.safe.com/2011/05/fmeevangelist78/ or in another test case to the relevant column of the table ). All ports are then connected to a shapefile- writer. The result is that the entries of the table with the names are output via the initiator port and the features (with geometries) themselves to the 'other' port. As result the geometries are shown in the inspector but aren't contained in the shapefile. I also tried exposing fme_feature_type on the reader then taking it on the FeatureReader, but this makes only appear the feature names in the .dbf file with no geometries at all. the hole thing is quite related to How to get a list of features from FME? Python? and Mark Ireland answer, and might be only causing troubles because Oracle is involved. I also checked: http://fmepedia.safe.com/articles/How_To/How-to-Read-All- Feature-Classes-from-Multiple-ESRI-Geodatabases / http://screencast.com/t/5OjFTrkiZL: this approach won't work with oracle as you have to specify credentials. so how could I get the FeatureReader to query the table names and put each in different shapefiles? Unfortunatly the workbench code I used cannot be pasted here because of size limitation. Any thoughts are greatly appreciated! Added Image shows the workbench used: note in the lower left corner the concatenated result based on which the FeatureReader is supposed to query the database Version information: #! LAST_SAVE_BUILD="FME(R) 2013 SP2 (20130605 - Build 13499 - WIN64)"


title: LiDAR Basics - How classes are determined at LiDAR data?

text: LiDAR data have classes as Medium vegetation, Buildings, Ground, Model Key and more. What determines these classes ? is it combination of intensity and scan angle rank ? how to reclassify these classes using ArcGIS ?


title: Calculating longest flow path and avg. slopes for catchments with qgis?

text: I am a hydrologist and their are 2 important pieces of information I need for modeling. I need to know the longest flow path for each sub catchment in my area of interest and I need to know the elevations of points specific distances along those paths. As an example, I may have a shapefile with specific boundaries that represent my catchment and an elevation raster. I would like to be able to determine the longest path that water flows in each of these defined areas. Once I have a shape file containing the longest flow paths for each catchment, I would like to find the elevation at 10% and 90% the length of each flow path and the distance between each point. I have not found either of these functions in anything except the ArcHydro tool bar for ArcGIS. I know that the individual tools for creating these things exist in Qgis ( through SAGA and GRASS ) but I am new to python and not a programer in any respect. I am hoping these tools are in an existing package and I am just too new to qgis/grass/saga to find them. Any help would be greatly appreciated. edit: When I mentioned that I was new to Qgis, I meant I have only been using it for about 18 months so I have some idea what I am doing.


title: Why use " from arcpy.sa import * "

text: Previously on the following post: How to get an extent for Raster using python A user edited the original Python script by adding the following code:               from arcpy.sa import *      Why did he do this and what exactly is it doing? I also received an error when I tried running the original script w/o the edit.


title: How to create non-overlapping street polygons from centerlines?

text: I want to create non-overlapping polygon layer from center line by applying the buffer to line layer in postgis. I am able to create polygon layer but how to avoid or resolve overlapped polygons.               select gid,st_multi(ST_union(ST_BUFFER(geom,10,'endcap=flat join=round'))) from roads group by gid      Using this query I am able to create polygon which is overlap at intersection ![overlap polygon at intersection](http://i.stack.imgur.com/UFxwH.png)   Due to 10 meter buffer distance polygon overlaps how to avoid or resolve this ![expected output](http://i.stack.imgur.com/Pq4K1.png)   I need expected output as non overlapping polygon with "gid" attribute. Any suggestion will always welcome. Thanks,Amit


title: ArcGIS Server 10.1 login error

text: Trying to login into my AGS server 10.1 via server manager throws an error see atatched pic ![enter image description here](http://i.stack.imgur.com/w32Pq.jpg)


title: How to motivate co-workers to "think outside the box"

text: Maybe this isn't really the right place to ask this, but I figure worst case scenario my question is closed as not relevant or too abstract. I am a GIS administrator, I manage our SDE and web map services, do ad hoc- projects for teams without their own GIS staff and assist numerous other GIS users across our organisation when they don't know how to do something. This last task is becoming more and more frequent as 90% of the other GIS staff seem incapable of thinking inventively. Unless they have every single button click fully documented, with screen-shots, they won't do it. I would like to know how other people encourage their co-workers to think through problems to come up with a method that works, rather than panicing or saying "I haven't done it before/for a long time/with this software" and passing it on to me. How can I get them to practice and play around with tools they haven't used before? At this point they won't even join two tables together unless I hold their hand. Thanks for any advice/suggestions


title: MODIS LDOPE mask_sds Error

text: I am using the LDOPE mask_sds tool to mask a MODIS image based on the quality flags for MOD13A3 NDVI. The full command that I am using is:               mask_sds -of=mask.hdf -fill=-3000 -sds='1km monthly NDVI' -mask="MOD13A3.A2000336.h25v02.005.2007112105435.hdf,'1km monthly VI Quality',0-1==00,AND,*,*,11-13==001" -meta MOD13A3.A2000336.h25v02.005.2007112105435.hdf      I am getting the following error:               Cannot find the SDS '1km monthly VI Quality'      Ignoring input masking option MOD13A3.A2000336.h25v02.005.2007112105435.hdf,'1km monthly VI Quality',0-1==00,AND     Cannot open the input HDF file *     Cannot open the HDF file *     Ignoring input masking option *,*,11-13==001     usage:          mask_sds -help [filename]               mask_sds -of=<output filename> -sds=<SDSname1>[,<SDSname2>[,...]]>          [-fill=<mask fill value>] -mask=<mask1>[,AND|OR,<mask2>[,...]] [-meta]         filename             where maskn=< filename>,<SDSname>,<bit_numbers operator bit_values>      It looks like it can not find the SDS (science data set) for the quality flag. I have tried several variations of '1km monthly VI Quality', including without quotes. This website lists SDS's for each MODIS product (I am using MOD13A3). I would appreciate any suggestions.


title: How do I reduce the number of drawn polygons on an Android Map?

text: We're using spatialite-android and Nutiteq to render Archaeological datasets. One of our sample data sets has 7000 polygons that turn into 400,000 objects in our renderer. Is there a common GIS method that deals with this problem? Unfortunately, our team hasn't done much GIS development and we don't know the right question to ask.


title: How to Set up Data Driven Pages using arcpy?

text: I've read a lot about the topic data driven pages. I'm writing a python script for ArcGIS and I've defined my stripmap using the `arcpy.StripMapIndexFeatures_cartography` command and now I have to set up Data Driven Pages with all the parameters it needs. My only problem is that I cannot understand how to do this. Is there a specific command in the arcpy library or a specific workaround to do this? Thanks in advance, Annalisa


title: Finding the Geometry within which a Point falls

text: I have the following query, which takes on average about 4 seconds. How might one improve the query such that the same response comes back faster?               SELECT * FROM geometries g          ORDER BY            ST_Distance(g.wkb_geometry,               ST_SetSRID(ST_Point(-77.8252923, 39.4328016), ST_SRID(g.wkb_geometry)))          LIMIT 1      (I understand that this actually finds the _nearest_ geometry to a point, which is the same as the geometry in which the point falls unless the point doesn't fall in any geometry. I'm open to seeing other approaches as well.)


title: Probability distribution from list of rasters

text: I have 100s (exactly 224) of wind distribution rasters. I want to create a probability distribution raster of wind speed using all these rasters, so we can predict the possible wind speed in any location. Is there any sufficient way in ArcGIS, Python, or R? if not is there any other solutions? Thank you very much is advance Kalus


title: QGIS cannot view raster loaded to PostGIS through raster2pgsql

text: I tried to load a raster GeoTIFF image to a PostGIS database using "raster2pgsql" following the OpenGeo Suite Manual. The command I used is:               C:\>raster2pgsql -s 4283 -I -C -M "C:\Program Files\GeoServer 2.2.5\data_dir\data\gcur\interp_clipped.tif" -F public.myraster | psql -d MyPostGISDB -h localhost -p 5432 -U postgres      The table was successfully created and found in "pgAdmin III" but not accessible through QGis. When I tried to load PostGIS layers in QGIS, it only showed "raster_columns" not the "myraster" table. In the PostGIS Manager window in QGIS, this table is marked as plain table and with "No Geometry". Can anyone give me some hints about it? Thanks in advance!


a simple question. Did anyone manage to use time manager with dates before 1900? Other than the solution suggested by Anita in Is it possible to change the year range of time manager in QGIS?


title: How to calculate midpoint given north/south latitude and east/west longitude?

text: First, sorry if this questions sound very simple. I am very new to GIS. Given four values of north/south latitude, east/west longitude how to calculate the center point that is correct for all cases on earth?                                 ----------(north latitude)---------------                      |                                         |                      |                                         |      (west longitude)|                ? (calculate this point) |(east longitude)                      |                                         |                      |                                         |                       ----------(south latitude) --------------      The dateline is an issue if you just use the following formula (N+S)/2 and (E+W)/2 to calculate the center.


title: Is it possible to perform 3D vector editing in QGIS?

text: I've searched the manuals and plugins and done some tests.. Is it in any way possible to perform 3D drawing and editing of vector objects for 3D mapping e.g. in a 3D shapefile or PostGIS 3D table?


title: QGIS 1.8 will not start on Windows 7

text: I have downloaded QGIS 1.8 and installed it on my windows 7 system. When I click on the desktop icon, nothing happens. I am currently running QGIS 1.7.4 with no problems and the installs of previous versions have always worked with no problems. I'm not sure what's going on here. Has anyone else experienced this and have a solution?


title: Running ArcPy script from ArcObjects?

text: the standard method requires to run ArcPy script in command line, i've struggled with it to run but without a result, i can `import arcpy`, set the workspace to the geodatabase path and describe the env variable using `arcpy.Describe(arcpy.env.workspace).name` it gives me the geodatabase name, but when i try to get featureclasses using `arcpy.ListFeatureClasses()` i get `[]`, i have made sure that all the environement variables in PATH and PYTHONPATH are ok . i would like to know if there's an alternative way to run ArcPy script, like it could be run in ArcMap console, is there a way to send scripts to ArcMap console using ArcObjects, what's the difference between ArcMap console and standard command line console? so any way to run an ArcPy script using ArcObjects 10 or Arcgis Engine 10 will be very welcome.


title: Google earth elevation data

text: I have a polyline in google earth. I want the elevation data (say in excel). I am able to see the elevation profile in google earth when I move along the polyline but Is there a way to extract the lat,long and elevation data in an external software? (say ArcGIS or excel) Any help is appreciated.


title: WCS row based security ?

text: I have an application where each user has individually generated geotiff files. Since the users pay for the generated geotiff, I have to make sure that the users can only see their own geotiff and not the geotiff of any other user. I'm using geoserver and PostgGIS and I plan to user WMS / WCS to deliver the data. The geotiff files are stored in the database, together with some additional data like an id. How can I setup such an security? Is it really necessary to create an own layer and role for each user? Or is there any way to create a security based on the id of the files in the database? best regards, DDD


title: 2D/3D Area parameter differences in Surface Volume calculation dialog box?

text: In dialog box of volume surface calculation and also in ASCII file there are 2 out come with names 2D Area and 3D Area. What is difference between these two numbers? What is included in 3D that is not in 2D? For example :                Plane_Height, Reference, Z_Factor, Area_2D, Area_3D, Volume      45,00, ABOVE, 1,000000, 644929,78213426, 716708,47208592, 17433798,859008      Note: Decimal point is a comma.


title: Points disappear when zooming in a certain scale

text: My points are disappearing when I zoom in past 300,000 in ArcMap. I tried to repair geometry, but it doesn't fix anything? I also checked the scale extents and those are turned off. I can export the points to a new shapefile/feature class, but I'd like to not have to do that. Any ideas what it could be? I'm using 10.1 on Windows 7 x64.


title: Spatial Analyst

text: does anyone know of some reasources to learn more about Spatial Analyst and raster analysis? I'd like to sharpen up for some job stuff, but I'd like to find some resources to practice on. Thanks.


title: WFS Layer Not displaying in Openlayers

text: I've been pulling my hair out for a week trying to figure this problem out. perhaps you can help. I'm trying to display a shapefile in the browser. It works great if I use WMS, but recently I've wanted to switch to WFS but have been unable to get it to work. What I get now is a map that is centered where I want it and an OSM background, but no WFS data In geoserver, I'm pretty sure that everything is set up OK. Here is the code I've used to request the WFS               var wfs = new OpenLayers.Layer.Vector("WFS",          {             strategies: [new OpenLayers.Strategy.Fixed()],             protocol: new OpenLayers.Protocol.WFS(             {                 url: "http://localhost:8080/geoserver/wfs",                 featurePrefix:"Van263kTweet",                 featureType: "vantest",                 featureNS: "file:///data/2636k_shp",                 maxFeatures: 100,                 geometryName: "the_geom",                 //featurePrefix: "",                 version:"1.0.0",                 srsName: new OpenLayers.Projection("EPSG:900913"),             }),stylemap: stylemap         });      in the browser it generates this request:               <wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.0.0"          maxFeatures="100"      xsi:schemaLocation="http://www.opengis.net/wfs      http://schemas.opengis.net/wfs/1.0.0/WFS-transaction.xsd"      xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">     <wfs:Query typeName="Van263kTweet:vantest" xmlns:Van263kTweet="file:///data/2636k_shp"/>     </wfs:GetFeature>      and the response I receive is this:               <?xml version="1.0" encoding="UTF-8"?><wfs:FeatureCollection      xmlns="http://www.opengis.net/wfs" xmlns:wfs="http://www.opengis.net/wfs"      xmlns:gml="http://www.opengis.net/gml" xmlns:Van263kTweet="file:///data/2636k_shp"      xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"      xsi:schemaLocation="file:///data/2636k_shp http://localhost:8080/geoserver/wfs?     service=WFS&version=1.0.0&request=DescribeFeatureType&     typeName=Van263kTweet%3Avantest http://www.opengis.net/wfs      http://localhost:8080/geoserver/schemas/wfs/1.0.0/WFS-basic.xsd">     <gml:boundedBy>         <gml:null>unknown</gml:null>     </gml:boundedBy>     <gml:featureMember>         <Van263kTweet:vantest fid="vantest.1">             <Van263kTweet:the_geom>                 <gml:Point srsName="http://www.opengis.net/gml/srs/epsg.xml#4326">                     <gml:coordinates xmlns:gml="http://www.opengis.net/gml"                          decimal="." cs="," ts=" ">                             -123.12660112,49.28351467                     </gml:coordinates>                 </gml:Point>             </Van263kTweet:the_geom>             <Van263kTweet:id>                 1             </Van263kTweet:id>         </Van263kTweet:vantest>     </gml:featureMember>     </wfs:FeatureCollection>      So as I'm generating the javascript is generating the request and the response returns GML, why isn't openlayers displaying the data? I've tried just about everything I can think of, and I'm out of ideas (changing the SRS, checking the namespace, reading geoserver log files, using different versions of WFS...). the complete html file that the page is based on is as follows:               var map;              // prepare to style the data         stylemap = new OpenLayers.StyleMap({           'strokeWidth': 11,           'strokeColor': '#ff0000'         });              map = new OpenLayers.Map('map-id', {projection: "EPSG:900913"});         osm = new OpenLayers.Layer.OSM( "Simple OSM Map");              var wfs = new OpenLayers.Layer.Vector("WFS",          {             strategies: [new OpenLayers.Strategy.Fixed()],             protocol: new OpenLayers.Protocol.WFS(             {                 url: "http://localhost:8080/geoserver/wfs",                 featurePrefix:"Van263kTweet",                 featureType: "vantest",                 featureNS: "file:///data/2636k_shp",                 maxFeatures: 100,                 geometryName: "the_geom",                 //featurePrefix: "",                 version:"1.0.0",                 srsName: new OpenLayers.Projection("EPSG:900913"),             }),stylemap: stylemap         });              map.addLayers([osm, wfs]);                   //set the map extent to vancouver on load         var proj = new OpenLayers.Projection("EPSG:4326");         var point1 = new OpenLayers.LonLat(-123.312607,49.065769);         var point2 = new OpenLayers.LonLat(-122.424088,49.390418);         point1.transform(proj, map.getProjectionObject());         point2.transform(proj, map.getProjectionObject());              //zoom/pan controls         map.addControl(new OpenLayers.Control.PanZoomBar({                     position: new OpenLayers.Pixel(2, 15)                 }));                 map.addControl(new OpenLayers.Control.Navigation());              //set default extent         map.zoomToExtent([point1.lon,point1.lat,point2.lon,point2.lat]);             //If provided as an array, the array should consist of four values              (left, bottom, right, top).      If anyone can help me out, I would really appreciate that. I've been through loads of other questions on StackExchange, but as of yet, none of the advice here has helped. Cheers.


title: QGIS 2.0.1 connecting to ESRI WMTS jpeg vs jpg

text: I'm currently trying to troubleshoot an issue I have connecting to ESRI WMTS services. Esri ArcGIS Server seems to provide the format as JPEG via the rest service, but in the get capabilites document they specify JPG - which it seems QGIS and Global Mapper don't want to recognise. Would the norm amongst mapping clients be JPG, JPEG or Both? I've included the URL to our server if anyone is keen to try it out, I've also attached an image of what I get in QGIS http://www.topofthesouthmaps.co.nz/arcgis/rest/services/CacheAerial/MapServer/WMTS ![QGIS to ESRI WMTS](http://i.stack.imgur.com/e8TZv.jpg) Thanks in Advance


title: How to simplify osm2po's routable graph?

text: I'm building a pgrouting routable network thanks to osm2po, the problem is that my network is too big. Indeed I want to create isochron and 1 node = 1 canvas circle, so I don't want too many nodes. I managed to reduce my network with grass (v.clean) but when I do that I'm losing all my costs (and reverse_costs). Is it another solution?


title: Overlay polygon on raster and save as new raster with same resolution

text: New to GIS and trying to learn. I'm trying to take a blue marble raster, overlay a countries .shp file, and save an image with the original blue marble image and the country outlines "burned into" the image. I've tried both ArcMap and QGIS and I know I can place both layers on the map and then save the map, but the resolution is modified to that of the screen. Ex. original raster was 5400x2700, but output raster is 1564x978. I don't see a way to modify the output files resolution either?? I searched the forum and didn't see anything like this, so I hope this is not something that has been covered ad nauseam. Thanks.


title: Can MapInfo Professional use data direct from a File Geodatabase?

text: Esri released its **File Geodatabase API** in mid 2011, and it seems that QGIS via GDAL's OGR FileGDB Driver added support for the File Geodatabase within about a year. By support, I mean simply being able to display, identify and query features using a **direct read** , with editing them not necessarily being required. Has such support of File Geodatabases also been implemented for MapInfo Professional at any release (or via an extension)? Please note that this question is only about `File Geodatabase (*.gdb)` rather than `Personal Geodatabase (*.mdb)` support because there is a big difference between the two formats.