title: Open-source software for managing terabytes of different formats of geospatial data?

text: I have a problem that I am sure many GIS Teams or Departments have which is how to find, keep track of and manage multiple different datasets from multiple sources. I have read the responses from a similar question asked two years ago Best practices for managing geospatial data I was wondering if there are any other software solutions, I know about VoyagerGIS, which looks good, but it is linked with ESRI and wondered if there was a similar Open Source version. To give you an overview of the current situation:   1. Terabytes of Imagery which changes every 6 weeks, I would like to keep historic   2. 1 Terabyte of MapInfo TABs   3. ESRI Shapefiles again changes every 6 weeks   4. 100s Gbs of CSV, GeoJSON, KML, SQLite, Excel spreadsheets, JPEGs, TIF files (some georeferenced, some with world files), these come in from multiple different sources. I guess my perfect solution would be a GUI interface to allow me to update metadata, search based on filename or even attribute, a visual image (like seen in VoyagerGIS), it show metadata, a list of attributes. However, I would like to keep the original format of the data rather than loading it all into PostGIS.


title: How to upload polygon shapefiles to my GPS?

text: I used QGIS to create a table (*.shp) containing 3 polygons. These polygons represent "forbidden zones" that I would like to export to my Garmin GPSmap62st so once on the countryside, I'll can see my position directly besides these "forbidden zones". Can you tell me how please?


title: Openlayers display a markers popup using a button out of the map

text:
I have the following code which addes 3 markers to the map along with there popup boxes what I want to do is have a list of location at bottom of page and using the id of the marker when click a place in the list it just make that places popup appear on the map. Code:                                               <meta http-equiv="content-type" content="text/html; charset=utf-8"/>                  <title>Open Street Map</title>

                                 The above code adds three markers to the page and when each is clicked displays there popup. What I want to do is add a list of locations along with a button for each and when the user clicks the location button that locations popup box appears. These button will sit of the map with in there own section of the site page. Can anyone help with how I can adapt my code. UPDATES Added:                 function clickmarker()                                                       var lon = -0.12519;         var lat = 51.51112;             popup = new OpenLayers.Popup("1",          new OpenLayers.LonLat(lon,lat).transform(new OpenLayers.Projection("EPSG:4326"),                                                     map.getProjectionObject()),                                 new OpenLayers.Size(200,100),                                 "example popup",                                 true);           map.addPopup(popup);          }


title: PgRouting: Driving_distance() or shootingstar_sp()?

text: I use OSM data imported into PgRouting database. I developped a web application with a map to use OSM data. I trace routes using shooting_star_sp() function. For example, to have time cost between 2969 and 70890 nodes, I use : `SELECT sum(time_cost_car) FROM ways, (SELECT * FROM shootingstar_sp('ways', 2969, 70890, 0.1, 'time_cost_car', true, true))as rt WHERE ways.gid=rt.gid;` And the result is 861 (seconds). I think it's correct. Now, I want to create polygon on my map to create isochron area. I followed this tutorial : http://underdark.wordpress.com/2011/09/25/a-closer-look-at- alpha-shapes-in-pgrouting/ It uses driving_distance() function. But the produced polygon seems to be wrong. I have performed the previous time cost calculation using the tutorial approach :               SELECT cost FROM driving_distance('SELECT gid AS id, source::int4 AS source, target::int4 AS target, time_cost_car::float8 AS cost FROM ways', 2969, 1000, false, false) WHERE vertex_id=70890;      And the result is 600 (seconds). That's wrong. So what is the difference between shootingstar_sp() and driving_distance() ? How can I create a correct isochron polygon ? Best regards, RudyWI * * * I think a part of the problem is resolved. The following request returns 1337 (seconds), it seems to be correct : `SELECT cost FROM driving_distance('SELECT gid AS id, source::int4 AS source, target::int4 AS target, time_cost_car::float8 AS cost FROM ways', 2969, 100000, true, false) WHERE vertex_id=70890;` However, the polygon produced by the following request seems to be wrong : `SELECT * FROM alphashape('SELECT id, ST_X(the_geom) AS x, ST_Y(the_geom) AS y FROM nodes JOIN (SELECT * FROM driving_distance(''SELECT gid AS id, source::int4 AS source, target::int4 AS target, time_cost_car::float8 AS cost FROM ways'', 2969, 30, true, false)) AS route ON nodes.id = route.vertex_id');` Do you kwow why ? I'm not sure about JOIN with "nodes.vertex_id". What is "vertex_id" ? Why not "edge_id" ? how to be sure "nodes.id" semantic corresponds to "vertex_id" ? The source parameter used into "driving_distance()" ("2969" in my example) is vertex_id ? edge_id ? or another id ? * * * I have performed other tests to illustrate the problem : `SELECT sum(cost) FROM (SELECT * FROM shortest_path('SELECT gid AS id, source::int4 AS source, target::int4 AS target, length::float8 AS cost FROM ways', 25069, 13530, true, false)) AS foo; Result => 16.7790377846424 SELECT sum(length) FROM ways, (SELECT gid, the_geom AS the_geom FROM shootingstar_sp('ways', 25069, 13530, 0.1, 'length', true, false)) as rt WHERE ways.gid=rt.gid; Result => 63.5863605213523` The result retuned by shootingstar_sp() seems to be OK. But to generate the polygon, I use driving_distance() which returns the same result as shortest_path(). This result seems to be wrong. So, what is the difference between shortest_path() (or driving_distance()) and shootingstar_sp() ?


title: Problems while installing Postgres

text: I was trying to install Postgres several times, but it didn't work out. I am using Mac and I tried to install Postgress.app, I got problems with database connection. I decided to install Postgress using HomeBrew. At the beginning, I uninstalled all the versions: Using `brew remove --force postgresql` to remove all versions I used brew cleanup posgtres as well Then, launched the installation using home brew               > ==> Downloading https://downloads.sf.net/project/machomebrew/Bottles/postgresql-     > Already downloaded:     > /Library/Caches/Homebrew/postgresql-9.3.2.mountain_lion.bottle.tar.gz     > ==> Pouring postgresql-9.3.2.mountain_lion.bottle.tar.gz     > ==> Caveats If builds of PostgreSQL 9 are failing and you have version 8.x installed, you may need to remove the previous version first. See:   https://github.com/Homebrew/homebrew/issues/issue/2510     >      > To migrate existing data from a previous major version (pre-9.3) of     > PostgreSQL, see:       > http://www.postgresql.org/docs/9.3/static/upgrading.html     >      > When installing the postgres gem, including ARCHFLAGS is recommended:      > ARCHFLAGS="-arch x86_64" gem install pg     >      > To install gems without sudo, see the Homebrew wiki.     >      > To reload postgresql after an upgrade:     >     launchctl unload ~/Library/LaunchAgents/homebrew.mxcl.postgresql.plist     >     launchctl load ~/Library/LaunchAgents/homebrew.mxcl.postgresql.plist      Any recommendation to install Postgres right way on my machine.


title: How to define own icon for qgsvertexmarker?

text: I want to use a small image instead of predefined icons of vertex marker in qgis. how to do this?               QgsVertexMarker * marker = new QgsVertexMarker(mPlugin->iface()->mapCanvas());     marker->setColor(Qt::black);     marker->setIconType(2);   // i want to use my icon here     marker->setCenter(A);       Any help is appreciated. Thank you.


title: Qgis Dufour crashes when using QgsMarkerLineSymbolLayer from a python plugin

text: I'm writing a python plugin to create symbology for Line Layers. I want to use Svg file for line symbology. I tried with default symbol, and also get a crash. The layer get rendered, but when I zoom or pan the map, QGis crashes. Here's my code:               # -*- coding: utf-8 -*-     from PyQt4.QtCore import *     from PyQt4.QtGui import *     from qgis.core import *     import os.path          class TestPlugin:              name = u'TestPlugin'              def __init__(self,iface):             self.iface = iface             self.pluginPath = os.path.dirname(__file__)              def initGui(self):             self.test1Action = QAction(u'Test1: line symbology',self.iface.mainWindow())             self.test1Action.triggered.connect(self.testLineSymbology)             self.iface.addPluginToMenu(self.name,self.test1Action)              def unload(self):             self.iface.removePluginMenu(self.name,self.test1Action)              def run(self):             pass              def testLineSymbology(self):             # setup properties of the symbol layer             name = os.path.join(self.pluginPath,'svg','symbols','SYMB1.svg')             prop = { u'name': name, u'size': '8.26', u'size_unit': 'MM' }             lineSymbolLayer = QgsSvgMarkerSymbolLayerV2.create(prop)             # Get a default marker symbol             markerSymbol = QgsMarkerSymbolV2()             # delete the default symbol layer             markerSymbol.deleteSymbolLayer(0)             # append the svg symbol layer             markerSymbol.appendSymbolLayer(lineSymbolLayer)             #             lineLayer = QgsMarkerLineSymbolLayerV2()             # this return True             lineLayer.setSubSymbol(markerSymbol)                  # Get a line symbol             symbol = QgsLineSymbolV2()             # Delete default symbol             symbol.deleteSymbolLayer(0)             # set my symbol layer             symbol.appendSymbolLayer(lineLayer)             # Setup renderer             renderer = QgsSingleSymbolRendererV2(symbol)                  # Setup a Memory Layer             layer = QgsVectorLayer(u'LineString?crs=epsg:4326',u'LineTestLayer',u'memory')             layer.dataProvider().addAttributes([                 QgsField(u'test_field',QVariant.Int),                 QgsField(u'angle_field',QVariant.Double)             ])             layer.updateFields()             # Setup a Feature             fet = QgsFeature()             fet.setGeometry( QgsGeometry.fromPolyline([                 QgsPoint(10,10),                 QgsPoint(12,12),                 QgsPoint(12,10)]))             fet.setAttributes([1,30.0])             fet.setValid(True)             layer.dataProvider().addFeatures([fet])             # Prepare layer for rendering             layer.setRendererV2(renderer)             layer.dataProvider().updateExtents()             layer.updateExtents()             #             QgsMapLayerRegistry.instance().addMapLayers([layer])      Is this wrong? Where do I mistake? Thanks for your answer. Al


title: How to limit decimal digit of GeoJSON

text: I am using postGIS. I am using ST_AsGeoJSON(geom) function to get GeoJSON data of Geom field. Its returning a long lat value that have 13 digits after decimal point. Is there any way to reduce or limit the decimal digits?


title: Openlayers How would i get a hover label over my polygon

text: make a hover/mouseover label for a polygon in Openlayers? I made a map with draw, delete and drag polygon function. Now I need to make a hover mouseover for this. How to get a hover label for polygon                       var polygonLayer = new OpenLayers.Layer.Vector("Polygon Layer");                  map.addLayers([wmsLayer, polygonLayer]);                  polyOptions = {sides: 40};             polygonControl = new OpenLayers.Control.DrawFeature(polygonLayer,                                             OpenLayers.Handler.RegularPolygon,                                             {handlerOptions: polyOptions} );                 var myStyles = new OpenLayers.StyleMap({                 "default": new OpenLayers.Style({                     pointRadius: "${type}", // sized according to type attribute                     fillColor: "#ffcc66",                     strokeColor: "#ff9933",                     strokeWidth: 2,                     graphicZIndex: 1                 }),                 "select": new OpenLayers.Style({                     fillColor: "#66ccff",                     strokeColor: "#3399ff",                     graphicZIndex: 2                 })             });                   var select = new OpenLayers.Control.SelectFeature(polygonLayer, {hover: true});             map.addControl(select);             select.activate();


title: Cost Surface No sense Results. GRASS

text: I have some doubts with the generation of a cost surface with _r.walk_ of GRASS. It's assumed that the resulting layer is measured in seconds, but when I try to convert seconds in days march (6h per day), the resulting layer has no sense. This is the workflow I'm following: **1st.-** I have several vector layers as roads, rivers, bridges, and a raster layer of elevation. **2nd.-** Rasterize the vector layers with a field named "weight" (PESO in spanish) where I have weighted this variables according to if they facilitate the walking or not. Reclassify is made if it's needed. **3rd.-** Find Slope (in percent) layer and reclasify as follow:   * 0-10=1;    * 11-maxValue=100 **4th.-** r.cross with all the rasterized layer. (roads, bridge, river and reclassified slope). **5th.-** I create a new friction layer using this formula proposed by Pandolf (1977) with _r.mapcalc_ :               M=1.5W+2.0(W+L)(L/W)^2+N(W+L)(1.5V^2+0.35V*abs(G+10))      **6th.-** With this layer measured in Watts I use _r.walk_ with default parameters. **7th.-** Once obtained this layer measured in seconds I try to convert it in hours dividing 3600 and then convert these hours in days march dividing 6 (6 hours per day). Well, the results are not what I expected, the max data of this layer is about 8.000... No sense, it should be something like 3 or 4 days, maybe. **_Here arethe files I'm using to do this workflow._** I don't know what I'm missing, maybe is the 4th step wrong? Maybe the 5th? Another point of view would be nice to me. If you need more info, please ask for it and I will add to the question. Thanks a lot!


title: Geoprocessing Tools on ArcGIS 10.2 and Visual Studio 2012 -issues

text: I am migrating one of the existing project (custom Geo-processing tool- _class library_ project) source code which runs on ArcGIS 10 SP4 and Visual studio 2010. This is working fine. I am could run the project and able to add this tool in ArcTool box. Publish this as GP service and consumed in application. I am moving this code to new machine which Windows 2012 server and ArcGIS 10.2 and Visual studio 2012. The source code is fetched from TFS and and code is built successfully after updating ArcGIS libraries to 10.2 When I run the project in VS2012 - class library (starting - ArcMap.exe with specific mxd document as input) break points are not hit. After doing some research I have found that ESRI Link here I have to update ArcMap.exe.config file. After updating config file to use Framework 4.0 break points are set but not getting hit. Whenever I run the class library project MXD is opens, but the break points are not getting hit and not able to debug the code. Project properties in VS2012 is same as I used in VS2010 (checked Register COM Exception) and enabled debugging. Question: a) Why VS2012 is not hitting breaking points? b) Dll is successfully built, but I am not seeing this custom tool box in "Add tool " box in ArcTool box. I am not sure what am missing or some changes in ArcGIS 10.2. Please help.


title: How to get feature fields and values?

text: How can I get features from EID? I cannot define a Geodatabase connection since it changes. I have to get it from `IUtilityNetworkAnalysisExt` or IEditor. I figured out how to get EIDs on downstream. I would like to get fields and values by EIDs.


title: How to do LOS Analysis in QGIS?

text: I am not able to do LOS analysis in QGIS. I have point Features and Raster SRTM Data. I have installed the LOS Plugin in QGIS but not getting any results. Kindly suggest?


title: Invalid MXD filename in a python script

text: I have a script that prints and exports a batch of MXDs. It runs fine as long as the lenght of the path and mxd are shorter. I can't tell what the 'sweet spot' is, but I know it runs fine on my c:/temp and doesn't run well when it gets to over 30 characters. The file structure is setup very lengthy and cumbersome and probably not going to change. Is there a way to temporary copy these files to temporary location on the C: drive then loop into the rest of script? Something like this.                import arcpy      import shutil           temp = r"C:/temp"      for Copymxd in InputMXDs           shutil.move(Copymxd, temp)           for MXDPath in temp:                #This is the rest of the script


title: ArcGIS 10.1 Custom Toolbar issue with addin wizard

text: I am having some issues. In ArcGIS, I created a custom toolbar (with the python addin wizard) to do the following: 1. Button to connect to feature class, 2-5. A series of comboboxes to view, edit data in a field in the feature class. I want each of the combo box classes to directly connect to the field. And be able to write new data from the boxes into the fields. Example: Combobox 1 to connect to field 1. So that when I deploy this a user can edit data without doing this the traditional way. Its a way to have inexperienced users edit data. I am new to python and willing to pay for the time required to help me as I am out of ideas. And ESRI doesn't have a solution either.               import arcpy     import pythonaddins     from arcpy import env               class ButtonClass1(object):         """Implementation for QuickTableEdit_addin.button (Button)"""         def __init__(self):             self.enabled = True             self.checked = False         def onClick(self):             print("Button1 Clicked")             env.workspace = "C:/data"             arcpy.CreateTable_management("C:/data", "test_table.dbf","template.dbf")             mxd = arcpy.mapping.MapDocument("CURRENT")             df = arcpy.mapping.ListDataFrames(mxd)[0]             dbf_Table = arcpy.mapping.TableView(r"C:\data\test_table.dbf")             arcpy.mapping.AddTableView(df, dbf_Table)             arcpy.RefreshTOC()             pass          class ComboBoxClass2(object):         """Implementation for QuickTableEdit_addin.combobox_1 (ComboBox)"""         def __init__(self):             self.items = [1,2]             self.editable = True             self.enabled = True             self.dropdownWidth = 'WWWWWW'             self.width = 'WWWWWW'         def onSelChange(self, selection):             global cb1             cb1 = self.value             pass         def onEditChange(self, text):             pass         def onFocus(self, focused):             pass         def onEnter(self):             pass         def refresh(self):             pass


title: Linking demographics data to zip code, census data

text: Where do I find population data to link to DMA / MMA or zip code? I need to be more granular then state. Any free data source? ArcGIS?


title: Attribute Indexes ArcSDE 9.3.1 for Attribute Joins on GUID field types

text: I create my joins on primary/foreign keys that are guids...I don't seem to be able to create an index on this field in my target table. After joining to this table it takes 4 minutes to open the source table. Can indexes be built on guids fields..if not what work arounds do you guys see? Thanks.


title: Get geometry of a Line Callout

text: How can I extract the geometry of a line callout and the border of an annotation featureclass in ArcGIS desktop. I want to transfert those to a separate featureclass. I'm able to get the Anchor point of the callout, but it's not enough for my needs. I have try to use following code, but it doesn't return nothing. In the code, I try to get a polyline, but I've also try to get a polygon with the same results. I'm using ArcMap version 9.3 but I also tried with 10 and I have same results. Quite strange also that the method QueryBoundary on IQueryGeometry says at run-time that it's not implemented!               Dim pActiveView As IActiveView = pMap     Dim pScreenDisplay As IScreenDisplay = pActiveView.ScreenDisplay     Dim pDisplayTransform As IDisplayTransformation = pScreenDisplay.DisplayTransformation          Dim pFDOGraphicsLayer As IFDOGraphicsLayer = pMap.Layer(0)   'The Annotation featureclass     Dim pFDOGraphicsLayerRead As IFDOGraphicsLayerRead = pFDOGraphicsLayer          pFDOGraphicsLayerRead.StartGeneratingGraphics(Nothing, pScreenDisplay, True, True, False)          Dim pTextElement As ITextElement = pFDOGraphicsLayerRead.NextGraphic          Do Until pTextElement Is Nothing           Dim pFormattedTextSymbol As IFormattedTextSymbol = pTextElement.Symbol           Dim pTextBackground As ITextBackground = pFormattedTextSymbol.Background                If TypeOf pTextBackground Is ILineCallout And pTextBackground IsNot Nothing Then                 Dim pLineCallout As ILineCallout = pTextBackground                 Dim pQueryGeometry As ESRI.ArcGIS.Display.IQueryGeometry = pLineCallout                 Dim pElement As IElement = pTextElement                 Dim pPointGeometry As IGeometry = pElement.Geometry                 pScreenDisplay.StartDrawing(pScreenDisplay.WindowDC, ScreenDisplay.ActiveCache)                 Dim pPolyline As IPolyline = New Polyline                      pPolyline = pQueryGeometry.GetGeometry(pScreenDisplay.WindowDC, DisplayTransform, pPointGeometry)                      If pPolyline IsNot Nothing Then                          'At this point the polyline is always nothing***                    MsgBox(pPolyline.Length)                 End If                      pScreenDisplay.FinishDrawing()             End If         pTextElement = pFDOGraphicsLayerRead.NextGraphic     Loop      Source: http://forums.esri.com/thread.asp?c=93&f=993&t=275411 Any input would be of great help.


title: Why does Erase give ERROR 999999 (Error executing function Invalid Topology Negative point valence)?

text: I am getting the following error when I try to erase a set of data from a file within a gdb. Any idea on what to do? I have done repair geometry multiple times on the input/erase datasets. I have 16GB RAM and close to a TB of free space. The dataset is approx 1.6m records. Arcmap 10 sp3. BTW how do I see the size of a dataset within a gdb? > Executing (Erase): Erase D:\SCRATCH\Projects\106\data\7_Merge.gdb\mergeall_1 > D:\SCRATCH\Projects\106\data\5b_ReclassVector.gdb\Lidar_d2b > d:\SCRATCH\Projects\106\data\6_Erase.gdb\srtm_m10m_e "1 Meters" Start Time: > Sat Oct 15 02:53:38 2011 Reading Features... Processing Tiles... ERROR > 999999: Error executing function. Invalid Topology [Negative point valence.] > Failed to execute (Erase).


title: Editing a shapefile in ESRI ArcMap (programmatically)?

text: Is this something that can be done? For instance if I created an overlay on a map and saved it as a shape file, I want to go and edit this ( **pre-existing file.** ) overlay...Either add points/remove points change color/fillstyle/etc. Ideally I want to avoid creating a new shp file and deleting the old one. Thoughts? Edit 1: OK, now that someone has said it can be done. The obvious question is How do you edit the shape file (in ArcObjects pro-grammatically) Edit 2: Looking at comments it appears I need to do something like this, but I'm not sure. I currently have a basic shapefile that I want to edit in ArcObjects in C#. I've created a shape file by creating a FeatureClass , Workspace and then obviously passing in my shape. In this case it's just a simple Polyline with a few points on it. What I'd like to do is have the ability to update this shapefile. I'm assuming I need to re-use that FeatureClass, Workspace and pass in a new shape. However, what I'm doing doesn't seem to work. This is where I'm at on the edit - any advice would be great.               IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)inWorkspace;     workspaceEdit.StartEditing(true);     workspaceEdit.StartEditOperation();          ComReleaser comReleaser = new ComReleaser();          IFeatureCursor featureCursor = inFeatureClass.Update(null , true);          IFeature feature = null;     while((feature = featureCursor.NextFeature()) != null)     {        int id = feature.Fields.FindField("shape");        IFeatureBuffer featureBuffer = inFeatureClass.CreateFeatureBuffer();        featureCursor.DeleteFeature(); //Delete the old         featureBuffer.Shape = (IGeometry)inShape;        featureCursor.InsertFeature(featureBuffer);     }          workspaceEdit.StopEdidtOperation();     workspaceEdit.StopEditing(true);      EDIT 3: _Final Solution_               featureCursor = inFeatureClass.Update(null ,true);     //DELETE old shape     featureCursor = inFeatureClass.Insert(true);     //Insert new      workspaceEdit.StopEditOperation();     workspaceEdit.StopEditing(true);


title: Can you search WMS services for a particular projection?

text: Is there an interface or website to search public WMS services for location/projection?


title: What Danish position format is this?

text: I have been given a list of locations within Denmark. Most of the are in decimal degrees. However, some of the positions are in different formats and I can't figure them out. Does anyone recognize these or have a guess as to what they might be? I have some that look like this:               lat: 55,800,643,991,678,900 lon: 12,297,413,013,875,400     lat: 55,804,257,011,041,000 lon: 1,229,334,301,315,240     lat: 55,807,527,964,934,700 lon: 12,333,790,976,554,100      and some that look like this:               lat: 583625 lon: 6236347     lat: 583722 lon: 6236625     lat: 562843 lon: 6218293


title: PostGIS Type Modifier Error

text: I'm trying to work an example from the "PostGIS Cookbook" and I'm getting an error I don't understand. I'm running the example using Postgres.app on a MacBook. Postgres.app is running PostgreSQL 9.1.12 with PostGIS 1.5.3. Here's the PL/Python function:               -- load.sql     CREATE OR REPLACE FUNCTION Geocode(address text)         RETURNS geometry(Point, 4326)     AS $$         from geopy import geocoders         g = geocoders.GoogleV3()         place, (lat, lng) = g.geocode(address)         plpy.info('Geocoded %s for the address: %s' % (place, address))         plpy.info('Longitude is %s, Latitude is %s.' % (lng, lat))         plpy.info("SELECT ST_GeomFromText('POINT(%s %s)', 4326)" % (lng, lat))         result = plpy.execute("SELECT ST_GeomFromText('POINT(%s %s)', 4326) AS point_geocded" % (lng, lat))         geometry = result[0]["point_geocoded"]         return geometry     $$ LANGUAGE plpythonu;      When I try to load the function into the database with this command:               psql geodjango < load.sql      I get this error:               ERROR: type modifier is not allowed for type "geometry"      I suspect this error has something to do with the fact that the book example assumes you're using PostGIS 2.0 and I'm using 1.5.3. However, it still doesn't make sense to me. I'm read about the ST_GeomFromText function in the PostGIS 1.5 Manual and this code appears to be correct. What am I doing wrong? Thanks.


title: How to add polygon layer into raster data?

text: I performed a supervised land-use classification with Erdas. Landsat data was used in ArcGIS to Digitize urban areas as polygon shapefile. I want to insert these polygons into the land-use raster data. How can I do this with ArcGIS?


title: Create a PDF with map imagery and feature data from ArcGIS API for Javascript?

text: I would like to allow the user to zoom to a particular zoom, within their preferred extent, and then send this information to a PDF as well as relevant information from visible features. The end result would be a PDF with map imagery, showing the location at the users zoom level, as well as grabbing data from the visible features and filling in a report-like set of data within the PDF as well. Is something like this possible or maybe even being done already by anyone?


title: ArcObjects - Geometry type select with Vbnet

text: How can I read the geometry type of shapefiles with ArcObjects and Vbnet? I want to check whether it is a point, polygon, or a line.


title: How to count occurrences of one field grouped by values of another?

text: Is there a way of calculating the number of unique text strings in a field in ArcGIS 10.1? I have a number of repeated strings in this field, but would like to calculate the number of unique strings in a field based on groupings from another field. E.g - **Current table:**               Field 1    Field 2        Apple      Green          Apple      Red            Apple      Red     Citrus     Yellow     Citrus     Green     Citrus     Orange     Citrus     Orange      **Result:**               Field 1    Field 2     Apple      2     Citrus     3


title: MapServer rendering at the lowest zoom levels

text: I am developing a web mapping application based on OpenLayers and using a WMS layer from MapServer (the source data for this layer is stored in PostGIS database) and I am facing the following problem: MapServer eliminates a lot of features at the lowest zoom levels. Zoom Level 10: ![Zoom Level 10](http://i.stack.imgur.com/971Ck.png) Zoom Level 7: ![Zoom Level 7](http://i.stack.imgur.com/oERje.png) Zoom Level 5:   ![Zoom Level 5](http://i.stack.imgur.com/AaS7b.png) Zoom Level 4:   ![Zoom Level 4](http://i.stack.imgur.com/4XCZM.png) Zoom Level 3:   ![Zoom Level 3](http://i.stack.imgur.com/9nF4n.png) I guess it is because MapServer doesn't display features with area less than one pixel. Is it possible to render such features at lowest zoom levels (for example as one pixel features)? UPD.: I'm talking about orange features, green ones - it is part of OSM tiles. Here is my mapfile:               MAP     SIZE 800 600     MAXSIZE 4096          IMAGECOLOR 255 255 255     IMAGETYPE PNG          OUTPUTFORMAT         NAME "png"         EXTENSION "png"         MIMETYPE "image/png"         DRIVER AGG/PNG         IMAGEMODE RGBA         FORMATOPTION "INTERLACE=OFF"     END          EXTENT -180 -90 180 90     PROJECTION         "init=epsg:4326"     END          WEB         METADATA             wms_onlineresource "${request.route_url('mapserver.ows')}"             wfs_onlineresource "${request.route_url('mapserver.ows')}"             ows_title "mcx_pk"             wms_enable_request "*"             # ows_service_onlineresource  "${request.route_url('mapserver.ows')}"             wms_srs "epsg:4326 epsg:3857 epsg:900913"         END     END          LAYER         NAME "fields"         CONNECTION "user=pk password=******* dbname=pk host=***.info"         CONNECTIONTYPE postgis         PROCESSING "CLOSE_CONNECTION=DEFER"         DATA "wkb_geometry from mask_f using unique ogc_fid using srid=4326"         TYPE line         DUMP TRUE         TEMPLATE dummy.html         PROJECTION             "init=epsg:4326"         END         EXTENT -180 -90 180 90         METADATA             gml_include_items   "all"             wms_include_items   "all"         END         CLASS             STYLE                 COLOR 255 136 0             END         END     END


title: Time Manage + Twitter in 2013

text: I am a geography teacher and would love to map tweets in QGIS for my pupils. I followed all the instructions but it seems Twitter have changed the way they stream their tweets. I can cope with the GIS side but I am not a Twitter developer. Is there any advice you can give a very busy teacher on getting those tweets?


title: Overlapping map tiles with mapnik's generate_tiles.py script

text: I'm trying to **generate map tiles** using **the mapnik's generate_tiles.py script**, but **the resulting tiles are overlapped**. ## Example The tiles have the following directory structure: `zoom_level/x_coordinate/y_coordinate.png`. Tiles **3/4/1.png** , **3/4/2.png** and **3/4/3.png** have some **countries shown twice** (for example: Turkey and a part of Italy) like in the following image: ![overlapping tiles](http://i.stack.imgur.com/7Zt8H.png) * * * ## Begin Edit Here is how the tiles should look at the **same zoom level** as the previous image (zoom level: 3; these tiles were generated using **MOBAC** with the " _OpenStreetMap Mapnik_ " Map source): ![OpenStreetMap Mapnik tiles generated with MOBAC](http://i.stack.imgur.com/50K9n.png) As you can see, in this case there aren't any tiles with overlapping content (or duplicate, if you prefer). Another strange thing is that **only the tiles from the middle of the map** (3.png and 4.png) **look similar** (cover _approximately_ the same area). The other tiles (0.png , 1.png and 2.png) **look totally different and I don't understand why!!!** **Note:** this is the first time I generate my own tiles so it's very likely that I'm doing a stupid mistake. It might be related to the **projection** I use (see below for the **osm.xml** ). ## End Edit * * * I tried setting the `self.m.buffer_size` to `0 (zero)` (like this page says; anyway, that's regarding labels, which I don't need) in:               def render_tile(self, tile_uri, x, y, z):         # ...              render_size = 256         self.m.resize(render_size, render_size)         self.m.zoom_to_box(bbox)         self.m.buffer_size = 128              # Render image with default Agg renderer         im = mapnik.Image(render_size, render_size)         mapnik.render(self.m, im)         im.save(tile_uri, 'png256')      Could it be a **problem with the projection?** Here are the important bits from my **osm.xml** file:               <Map bgcolor="#f2efe9" srs="+proj=latlong +datum=WGS84">      and               <Layer name="countries" status="on" srs="+proj=latlong +datum=WGS84">      **Can someone help me? I already spent almost a week searching and documenting myself about this process and I got stuck here.**


title: Calculate View Angle

text: I am looking for the most accurate way of calculating the view angle (Azimuth, Elevation) from a reference point (my location) A to another point B. I have the Lat, Lon, and Altitude for each point. Currently I am converting to local ENU coordinates using the Wikipedia formulas for ECEF to ENU and using that to calculate the angles.               Az = atan2(E,N)     El = atan(U / sqrt(E^2 + N^2)      Is there a more accurate way of calculating these Angles? Maybe a different coordinate system I should be projecting to? Is this the correct type of Angle, should a rhumb line azimuth or something else be used instead? I have also looked at using Vincenty's formula for calculating the azimuth but it does not take into account the altitude of each point. A change in Altitude using ENU seems to slightly change the Azimuth which seems correct(?). When comparing Vincenty Azimuth with the ENU Azimuth on the surface (Altitude 0), the answers are extremely close but start to differ if the points are not on the surface.


title: How does mapscript queryByShape work?

text: I am trying to return the results of queryByShape call from a (large) raster file given a user specified polygon. ![enter image description here](http://i.stack.imgur.com/1cQne.png) For example, when a user plots the following polygon on an openlayers layer as shown in this diagram, I send the coordinates of the polygon vertices to the python mapscript backend. With python mapscript, I re-create the polygon as a mapscript shapeObj.               line = mapscript.lineObj()     line.add( mapscript.pointObj( 16.484375, 59.736328125 ) )     line.add( mapscript.pointObj( 17.4951171875, 55.1220703125 ) )     line.add( mapscript.pointObj( 24.74609375, 55.0341796875) )     line.add( mapscript.pointObj( 22.5927734375, 61.142578125) )     line.add( mapscript.pointObj( 16.484375, 59.736328125) )     poly = mapscript.shapeObj( mapscript.MS_SHAPE_POLYGON )     poly.add( line )      Then, I experiment with the queryByShape call like this:               mapfile = '/path/to/my/mapfile.map'     mymap = mapscript.mapObj( mapfile )          layer = mymap.getLayer(0)         layer.queryByShape( mymap, poly )      The returned result is an integer "1", which isn't what I was expecting. Unfortunately, there isn't much by way of documentation for queryByShape so I am not sure what I am doing wrongly (http://mapserver.org/mapscript/querying.html#by-shape) How does mapscript queryByShape work? I am actually trying to return a shape object where I can extract information about the selected shape (polygon defined by user) and determine   1. The total area of this shape,   2. The % in this shape which is of a specific RGB ("blue") value,   3. The % in this shape which is of a specific RGB ("green") value etc.  Here's the LAYER definition in my mapfile.               LAYER         NAME 'GLOBCOVER_L4_200901_200912_V2.3.color'         TYPE RASTER         DUMP true         TEMPLATE 'fooOnlyForWMSGetFeatureInfo'         EXTENT -189.001389 -115.106069 188.998611 140.108847         DATA '/media/GLOBCOVER_L4_200901_200912_V2.3.color.tif'              METADATA             'ows_title' 'GLOBCOVER_L4_200901_200912_V2.3.color'         END              STATUS OFF         TRANSPARENCY 100              PROJECTION             'proj=longlat'             'ellps=WGS84'             'towgs84=0,0,0,0,0,0,0'             'no_defs'         END     END


title: How to crop google hybrid map in QGIS?

text: I download google map from openlayers plugin in QGIS, it was like this: ![enter image description here](http://i.stack.imgur.com/hVk7c.jpg) How can I crop the view to get a certain area from the map like this: ![enter image description here](http://i.stack.imgur.com/0hmLi.jpg)


title: Count the area of rasterdata in specific shapefile polygons (residential density in different districts)

text: I'd like to count the exact area of three different residential classes (low, medium, high) by district of a shapefile which contains 225 districts. like: district 1 has: XXX sqm low, XXX sqm medium, XXX sqm high. distict 2 has: XXX.... the residential classes are a file system raster, the districs are one shapefile. I can count the total area of the residential classes of the whole study area. but I can't count the area for every unique district I have in the shapefile. (Maybe I could count the area for 1-3 specific districts I pick but for 225 it would be too time-consuming) After all I'd like to see the area of each district in the attribute table of the district-shapefile and be able to export the data to Excel. I hope you can help, thanks.


title: ArcPy not running on standard console

text: I know this is a newbie ArcPy question, but here is the deal: While I'm inside ArcMap 10, the Python window gives me all output correctly, like this:               >>> arcpy.env.workspace = r"C:\uber.gdb"     >>> arcpy.ListDatasets("*","ALL")     [u'uberlandia_urbano', u'Idw_pop2006', u'Idw_pop2000']      But, if I open cmd window (Run > cmd) and type the exact same commands, it will give me an empty list as a result. I'm pretty sure that the Python version used is the one tied up with arcgis, since I can import arcpy. Any chance that arcpy does not work on the standard python console? Thanks


title: Divide a Square to a serval equal squares using Acrgis

text: hi i'm looking for an algorithm to split a square to a SEVERAL equal squares: i have only :   * x min , y min    * x max , y max please see picture (not valid picture only to describe problem): ![enter image description here](http://i.stack.imgur.com/GW6IS.png) **any suggestions please.**


title: Problem when calculating longitude and latitude for points

text: I have a PostGIS database with SRID=900913 (Google Mercator) and I would like to know the latitude and longitude of points. I use this command:               ST_X(ST_TRANSFORM(the_geom,4283))      Following what is written here, but it doesn't work. I get the error:               ERROR: AddToPROJ4SRSCache: Cannot find SRID (4283) in spatial_ref_sys     SQL state: XX000.       I tried to add 4283 in the spatial_ref_sys:               INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 94283, 'epsg', 4283, '+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs ', 'GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]]');      But I still have the same error. What does that mean?


title: GIS applications used in Civil Engineering & sources for acquiring more GIS knowledge

text: I have finished my masters in GIS in August 2009. Since then I have been working full time. My background is in Civil engineering. I do not know many people from similar background getting into this field. I truly want to learn GIS applications used in my field. I have been listening to few podcasts (Esri, Geospatial) that barely keeps me up with rapidly changing GIS technology. I hoping one day I could start a consulting business of mine own. Anyway, the question that I want to ask is where can acquire more knowledge related to my field? Do most us acquire GIS knowledge by attending seminars, online courses, webniars, and carrying out projects at work? I would love to get some responses.


title: How to edit ArcGIS Server .json?

text: Is it possible to edit the .json format, which is created by default when publishing a service on ArcGIS for Server? For a project we have the following problem:   * Our by default created .json for a Feature Service Layer looks e.g. like **this**.    * We want to change some of these .json elements like edit the objectID name or adding another element like "api_key"? Because there is also another group, which works simultaneously on another application, both teams should use the same data coming from an ArcGIS Server. So....is there a way to edit the json schema?


title: How to develop configuration interface for custom widget in ArcGIS Flex Application builder?

text: I'm trying to know if it's possible to develop a custom configuration interface for custom widgets in ArcGIS Flex Application Builder. The standard widgets have a nice and shiny configuration interface, with buttons and text-box, but I would like to know if it's possible (and how) to have similar GUI for my custom widgets as well. I mean, editing XML is note very user friendly and will not fit in what I want to do. If it's not possible at the moment, is this something planned for this builder?


title: Is schema definition compulsory for WFS-T

text: Most tutorials on WFS transaction uses a schema and geom attributes in WFS request like this               var wfs_layer = new OpenLayers.Layer.Vector("Editable Features", {     strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],     //projection: new OpenLayers.Projection("EPSG:26910"),     protocol: new OpenLayers.Protocol.WFS({     version: "1.1.0",     url: "http://localhost/geoserver/wfs",     featureNS :  "http://www.opengeospatial.net/cite",     featureType: "wfst_test",     **geometryName: "the_geom",**     **schema: "http://localhost/geoserver/wfs/DescribeFeatureType?version=1.1.0&;typename=cite:wfst_test"**     })     });      can we omit this or is it compulsory to include them?


title: Hibernate Spatial Maven Issue

text: I'm following the tutorial at http://www.hibernatespatial.org/tutorial- hs4.html. I've created a maven project, and imported it into Eclipse. However, when I run `mvn dependency:list`, I get the error `Missing artifact org.postgis:postgis-jdbc:jar:1.5.3`. The same error is reported from within Eclipse. Any idea on what I may be doing wrong?


title: How to create KML file with OpenLayers?

text: I try to create kml file from features. My code               var feature=myVecLayer.getFeatureByFid('filedata.1');                         var format;                             function GetKMLFromFeatures(feature) {                                  format = new OpenLayers.Format.KML({                                     'maxDepth':10,                                     'extractStyles':true,                                     'internalProjection': new OpenLayers.Projection("EPSG:4326"),                                     'externalProjection': new OpenLayers.Projection("EPSG:4326")                                 });                             alert(format);                                    alert(format);                             return format.write(feature);                             }      After i use this library to save file https://github.com/dcneiner/Downloadify                                           Downloadify.create('downloadify',{                                 filename: function(){                                                             alert("okey");                                     return  'kml_test.kml';                                 },                                 data: function(){                                     return GetKMLFromFeatures(features);                                 },                                 onComplete: function(){                                     alert('Your File Has Been Saved!');                                 },                                 onCancel: function(){                                     alert('You have cancelled the saving of this file.');                                 },                                 onError: function(){                                     alert('You must put something in the File Contents or there will be nothing to save!');                                 },                                 transparent: false,                                 swf: 'media/downloadify.swf',                                 downloadImage: 'images/download.png',                                 width: 100,                                 height: 30,                                 transparent: true,                                 append: false                             });      But nothing heppens. If firebug i get error               TypeError: f.el is null           ...tainer.id};if(f.options.enabled===false)f.enabled=false;if(f.options.transparent...      But more strange that i not get arert `alert("okey");` Whats wrong? Or if someone know better way to get kml file please give me advice.


title: Summary Statistics tool for Mean Z values

text: I'm using ArcGIS 10.2 The summary statistics tool is not functioning in the way that I want. I made an elevation layer that is a DSM, and then in my buildings layer I determined the height for buildings by creating a random points sample. Now that all the random points in each building have an elevation, or z-value, I want to get the mean for each building. I used the summary statistics tool, but it calculates a summarized mean for all 482100 buildings combined. I want the z-mean summarized for each individual building, this is possible because the buildings all have an unique ID number. See Attached Photo, Thank you. ![Z mean value operation](http://i.stack.imgur.com/iUz2a.jpg)


title: How to combine use gradient symbols for layer properties excluding 0?

text: I want to visualize my data by using graduated symbols.   Now, there are some objects with value "0", they should not be printed on the map (which I solved by using rule-based symbols).   But, anyway, "0" is included when automatically generating the (equal interval) classes. This is important because my data are values of more or less 1800 to 2000.   Also three different sizes are applied, but I don't think that will be important for this question. How do I exclude "0", when assigning gradual colours? ![enter image description here](http://i.stack.imgur.com/9ucAO.jpg)


title: Removal of Qgis 2.2 after 2.4 install, while keeping WFS and WMS features

text: Can i safely uninstall Qgis 2.2 after I've installed the 2.4 version and keep my WFS and WMS features? I'd rather be certain, so I don't have to setup the whole thing again.


title: How can I trace an ANALYSIS_EXPORT error when building QGIS?

text: When building qgis1.8 using QT 4.7.4, I encountered the following error. I am unable to trace it:               src/analysis/raster/qgsrastercalculator.h:31:24: error: variable ‘ANALYSIS_EXPORT QgsRasterCalculatorEntry’ has initializer but incomplete type      Earlier there was an error for `CORE_EXPORT` too, i have solved it using `DEFINES += GUI_EXPORT= CORE_EXPORT=` Any help is greatly appreciated.


title: merge_management

text: Here is the edited code: import arcpy import os from arcpy import env # set the environment settings env.workspace = "C:\RG\Kings County\Test.gdb" # fcs = arcpy.ListFeatureClasses("Sub_NVF_*") output = "C:\RG\Kings County\Test.gdb\Merge3FC" arcpy.env.overwriteOutput = 1 # this code overwrites existing output try:                   FC1 = "Sub_NVF_050413"         FC2 = "Sub_NVF_052503"         FC3 = "Sub_NVF_052509"              fldMap_Mean = arcpy.FieldMap()         fldMap_NDVI = arcpy.FieldMap()         fldMap_ID = arcpy.FieldMap()         fieldMappings = arcpy.FieldMappings()              Mean1 = "M050413"         Mean2 = "M052503"         Mean3 = "M052509"         NDVI1 = "NDVI050413"         NDVI2 = "NDVI052503"         NDVI3 = "NDVI052509"         Fld_ID = "Field_ID"              fldMap_Mean.addInputField(FC1,Mean1)         fldMap_Mean.addInputField(FC2, Mean2)         fldMap_Mean.addInputField(FC3, Mean3)              fldMap_NDVI.addInputField(FC1, NDVI1)         fldMap_NDVI.addInputField(FC2, NDVI2)         fldMap_NDVI.addInputField(FC3, NDVI3)              fldMap_ID.addInputField(FC1, Fld_ID)         fldMap_ID.addInputField(FC2, Fld_ID)         fldMap_ID.addInputField(FC3, Fld_ID)              #set name of new output field         Mean_type = fldMap_Mean.outputField         Mean_type.name = 'Mean'         Mean_type.aliasName = 'Mean'         fldMap_Mean.outputField = Mean_type         fieldMappings.addFieldMap(fldMap_Mean)              NDVI_type = fldMap_NDVI.outputField         NDVI_type.name = 'NDVI'         NDVI_type.aliasName = 'NDVI'         fldMap_NDVI.outputField = NDVI_type         fieldMappings.addFieldMap(fldMap_NDVI)              FieldID_type = fldMap_ID.outputField         FieldID_type.name = 'Field_ID'         FieldID_type.aliasName = 'Field_ID'         fldMap_ID.outputField = FieldID_type         fieldMappings.addFieldMap(fldMap_ID)              arcpy.Merge_management([FC1, FC2, FC3], output, fieldMappings)      except Exception as e: print e Above code works fine. Although mean and NDVI field name are feature class specific, it perfectly does the job. However, when I use for loop and try to read field name as as Mean = "M"+ fname (as shown below), it gives me an error saying "Error in adding input field to field map" import arcpy import os from arcpy import env env.workspace = "C:\RG\Kings County\Test.gdb" fcs = arcpy.ListFeatureClasses("Sub_NVF_*") output = "C:\RG\Kings County\Test.gdb\MergeAllV1" try: for fc in fcs: fcname = fc.split('_', 1)[1]                   fldMap_Mean = arcpy.FieldMap()         fldMap_ID = arcpy.FieldMap()         fieldMappings = arcpy.FieldMappings()              Mean = "M"+fcname         Fld_ID = "Field_ID"              fldMap_Mean.addInputField(fc, Mean)         fldMap_ID.addInputField(fc, Fld_ID)                   Mean_type = fldMap_Mean.outputField         Mean_type.name = 'Mean'         Mean_type.aliasName = 'Mean'         fldMap_Mean.outputField = Mean_type         fieldMappings.addFieldMap(fldMap_Mean)              FieldID_type = fldMap_ID.outputField         FieldID_type.name = 'Field_ID'         FieldID_type.aliasName = 'Field_ID'         fldMap_ID.outputField = FieldID_type         fieldMappings.addFieldMap(fldMap_ID)          arcpy.Merge_management(fc, output, fieldMappings)     print ("done")      except Exception as e: print e However, if I try to merge all feature classes using fields whose name is same in all files (as shown below), it writes only the last file. Any clue? import arcpy import os from arcpy import env env.workspace = "C:\RG\Kings County\Test.gdb" fcs = arcpy.ListFeatureClasses("Sub_NVF_*") output = "C:\RG\Kings County\Test.gdb\MergeAll" try: for fc in fcs:                   fldMap_ID = arcpy.FieldMap()         fldMap_Name = arcpy.FieldMap()         fieldMappings = arcpy.FieldMappings()              Fld_ID = "Field_ID"         F_Name = "File_Name"              fldMap_ID.addInputField(fc, Fld_ID)         fldMap_Name.addInputField(fc, F_Name)              FieldID_type = fldMap_ID.outputField         FieldID_type.name = 'Field_ID'         FieldID_type.aliasName = 'Field_ID'         fldMap_ID.outputField = FieldID_type         fieldMappings.addFieldMap(fldMap_ID)                   FName_type = fldMap_Name.outputField         FName_type.name = 'File_Name'         FName_type.aliasName = 'File_Name'         fldMap_Name.outputField = FName_type         fieldMappings.addFieldMap(fldMap_Name)               arcpy.Merge_management(fc, output, fieldMappings)      except Exception as e: print e


title: GeoTools Omitting Quotes in its Queries?

text: I have a table in my PostgreSQL 8.4 with a column 'Shape' of type 'geometry'. When I start up Tomcat my GeoServer throws an exception:               24 Sep 16:20:25 WARN [geotools.jdbc] - Error occured determing srid for my_table.Shape     org.postgresql.util.PSQLException: ERROR: column "shape" does not exist      In my PostgreSQL log I can see the corresponding query and error:               ERROR:  column "shape" does not exist at character 77     STATEMENT:  SELECT ST_SRID("Shape") FROM "core"."my_table" WHERE Shape IS NOT NULL LIMIT 1      Just to explore, I tried that same query from psql, and it failed in the same way. I tried it again, except I put double quotes around the 'Shape' in the WHERE clause:               gtda=# SELECT ST_SRID("Shape") FROM "core"."my_table" WHERE "Shape" IS NOT NULL LIMIT 1;      st_srid     ---------         4326     (1 row)      What is generating that query, why is it omitting the quotes, and is there something I can do about it? I have a number of similar tables, all of which are popping this exception.


title: How can I remove the menu bar in QGIS?

text: I want to create a clean user interface for my students to use QGIS. I tried to remove the menu bar using the customization in setting but I am not able to remove everything, e.g. Database, plugin, etc. Is there a way to remove them?


title: Comparing TIGER 2013 with TIGER 2006

text: I'm doing a project when I want to compare TIGER 2013 data with the TIGER 2006 data. Essentially, the Census updated the data between 2002-2008, and a number of the counties were updated in the 2006 version and many were updated after 2006. By comparing the two versions I hope to figure out which counties are similar in the 2006 dataset and which have been updated. Trouble is, I've never worked with the TIGER data, or done comparisons of such large datasets. To make it more complicated, the 2006 data is in TIGER/Line ASCII format, while the 2013 data is a shapefile. I'm only interested in the streets layer. Any pointers on getting this done efficiently? I have access to QGIS 2 and have some experience manipulating large data files using python tools like shapely. My university has the license for ARCGIS, but I use Linux primarily, so would like to avoid that if possible.


title: Can't run node-mapnik-sample-code

text: I'm using ubuntu 12.04. I installed mapnik from packages as described here: https://github.com/mapnik/mapnik/wiki/UbuntuInstallation now I want to run node-mapnik-sample-code from here: https://github.com/mapnik/node-mapnik- sample-code The example is:               node ./examples/simple/render.js examples/stylesheet.xml map.png      but if I run this command it shows the errormessage               /home/tiletest/node-mapnik/examples/simple/render.js:27     map.loadSync(stylesheet);         ^     Error: Could not create datasource for type: 'shape'  encountered during parsing of layer 'world' in Layer at line 11 of     'examples/stylesheet.xml'         at Object.<anonymous> (/home/tiletest/node-mapnik/examples/simple/render.js:27:5)         at Module._compile (module.js:456:26)         at Object.Module._extensions..js (module.js:474:10)         at Module.load (module.js:356:32)         at Function.Module._load (module.js:312:12)         at Function.Module.runMain (module.js:497:10)         at startup (node.js:119:16)         at node.js:902:3      do I have to edit the stylesheet.xml or how can I run the example ?


Apparently if you stand at Greenwich then your GPS will not say you are at longitude zero. This is because the GPS reference meridian is about 100 metres to the east of the prime meridian. My question is: **Why is the GPS Reference Meridian 100m to the East of the Prime Meridian?**


title: Can Geodatabase files from ArcMap be imported and used in Qgis?

text: I'm looking into Qgis and am doing some tests, but one thing I need to be able to do is use geodatabase files from ArcMap and ArcPad. Is that something that can be done and with relative ease? I have been playing around with it for a few days and I haven't been able to figure it out. I've also looked at several threads, but I haven't seen anything specific this topic. Any help or advice would be greatly appreciated.


title: Using AcuGIS with QGIS

text: Has anyone use **AcuGIS** , a web hosting service for PostGIS and GeoServe? I am curious to see how it works with QGIS. I am looking at some cloud hosting services for PostGIS and I would like to be able to access the database in QGIS and pgAdmin. Has anyone had any success or issues with AcuGIS. Or can you recommend any other cloud hosting services.


title: Software for Stitching Aerial Photographs?

text: I am looking for software that is capable of stitching together large number of aerial photographs. The photos are taken with a high quality DSLR from a small helicopter, they cover about 100 meters, they have a small amount of overlap, they also have location data stored in the exif header, are about 7 mb in size, and there are usually around 30000 of them. I was just wondering if anyone had any suggestions. So far I have tried:   * Autopano Giga 3 - no good results yet just some weird stitches   * Palentier (yes I know it doesn't really stitch the photos) - Hard to use as I am lacking the direction info, and height above ground info in my exif headers. The height above ground can be calculated by subtracting the ground elevation from the elevation stored in the photo exif but all photos must have the same subtraction value and the terrain height varies greatly.   * Agisoft Photoscan - Doesn't like the large amount of images, but I am currently trying it with smaller batches.


title: Which Bands from LISS III image are to be used for NDVI?

text: I have some IRS LISS III & IRS P6 IMAGES. What will be the correct bands to make NDVI image and how can I calculate from the data the amount of vegetation present or the change detection?


title: Printing possible for openlayers GeoJson WFS with Mapfish printing service and wkhtmltopdf?

text: I searched lots of places but I couldn't find proper examples or documentation to do printing using Mapfish Printing service or Wkhtmltopdf. My application has a few WFS layers from openlayers and geoserver only (no wms ). So is printing possible from the above two methods or not? Also, I would like to add that I am using osm as baselayer. Actually, I am requesting mapfish printing (integrated with geoserver) service through url. Through my url, I am getting wms layer in print out but not osm layer and also if you can tell me how I can request for geojson wfs layer through url? Mapfish printing working properly and this is my info.json details: > > {"scales":[{"name":"1:25,000","value":"25000"},{"name":"1:50,000","value":"50000"},{"name":"1:100,000","value":"100000"},{"name":"1:200,000","value":"200000"},{"name":"1:500,000","value":"500000"},{"name":"1:1,000,000","value":"1000000"},{"name":"1:2,000,000","value":"2000000"},{"name":"1:4,000,000","value":"4000000"}],"dpis":[{"name":"75","value":"75"},{"name":"150","value":"150"},{"name":"300","value":"300"}],"layouts":[{"name":"A4 > portrait","map":{"width":440,"height":483},"rotation":true},{"name":"Legal","map":{"width":440,"height":483},"rotation":false}],"printURL":"http://xxxxxxxxxxxx.xxxxxx.xxx:8080/geoserver/pdf/print.pdf","createURL":"http://xxxxxxxxxxxx.xxxxxxxx.xxx:8080/geoserver/pdf/create.json"} Details for osm layer which I am requesting in url is: > {"baseURL":"http://tile.openstreetmap.org/",' + '"opacity":1,' + > '"singleTile":false,' + '"customParams":{},' + '"type":"OSM",' + > '"extension":"png",'+ '"tiled":"false",'+ '"tileSize" : [256, 256],'+ > '"maxExtent" : [-20037508.34,-20037508.34,20037508.34,20037508.34],'+ > '"resolutions" : [19567.87923828125],'+ '}],' + //am getting error,could not > be able to fetch tiles 404 error. One more thing I would like to ask: is it feasible to have the same screenshot (display) on pdf which the user can see at that moment in browser and click for print request? Since I am using the service and reading it from server side, how can I get the pdf of what the user is seeing in client side if this is possible?


title: Converting XYZ tile request to WMS request?

text: I am wondering if it is possible to convert XYZ layer request e.g (http://dev.openlayers.org/docs/files/OpenLayers/Layer/XYZ-js.html) to a regular WMS request? (Basically if it is possible to convert the numbers Z/X/Y to a bounding box) I have a working WMS on the server side that can return the needed tiles, if the convertion is possible then I would convert the Z/X/Y on the server side and pass it to my WMS service to get the tile. I need this because I have an Arcgis generated cache and it would be nice to also support the XYZ standart for the same cache.


title: Openlayers - custom property for layer object

text: I'm fairly new to OpenLayers, and have written a little function to switch the legend image for WMS layers on and off if the layer is visible or not. This is all much neater if I can define my own custom property for the WMS layer object, e.g legendurl: 'http//....' While this works, I'm wondering if it is ok/advisable to create custom properties on the fly in this way? Cheers Matty


title: Raster domain output polygon has is an empty feature class. Shape_Length and Shape_Area = 0

text: I am trying to convert a 32-bit floating point raster with 205289 columns and 120615 rows to a polygon using Raster Domain. It creates a feature class that does not contain any polygon. There is a record in the attribute table but it says zero in the Shape_Length and Shape_Area fields. Any ideas why I am getting this output? Or other suggestions how to convert this raster to a polygon? I just need a simple polygon to delineate the boundary of my raster.


title: Signing ArcGIS Add-In with a post-build event

text: I'm attempting to automate signing of an ArcGIS Add-In using a post build event in visual studio that triggers a batch file (which subsequently triggers ESRISignAddIn utility). Running the batch file manually through command prompt works fine, the Add-In is signed successfully. However if I attempt to run the batch file in the post-build event, it returns an error saying the Add-In file cannot be found. I've noticed that by the time the build process gets to the post-build stage, only the DLLs are created - the Add-In file is not, hence the batch file cannot find it. It appears the Add-In file is created AFTER the post-build event (illogical, I know). I have tried adding the batch commands to the post-build event command line (in project settings > compile > build events), and have also tried adding this directly to the project file XML:               <Target Name="AfterBuild">         <Exec Command=""$(ProjectDir)Post_Build.bat" "$(TargetDir)my.esriAddIn"" />     </Target>      Seems that 'post build' and 'after build' events are both executed before the esriaddin file is built, useless! Has anyone had any experience working with post-build events and ArcGIS Add- Ins? **I want to know if it's possible to force the Add-In file (my.esriaddin) to be created before the post-build event.** My last resort solution is to execute the build and then the batch file using a build runner (e.g. finalbuilder), but would ideally like to keep everything in the same project.


title: Sunlight Layer into CartoDB

text: I have a list of tweets on a torque map made with CartoDB. I want to add an additional layer that shows the sunlight around the world at certain times, similar to the Twitter #sunrise map found here http://cartodb.s3.amazonaws.com/static_vizz/sunrise.html?title=true&description=true&search=false&shareable=true&cartodb_logo=true&layer_selector=false&legends=false&scrollwheel=true&sublayer_options=1%7C1&sql=&zoom=2&center_lat=22.917922936146045&center_lon=51.328125# Does anyone know how to add a layer like this?


title: Raster' object is not iterable

text: I'm trying to calculate the percentage of cell counts in new raster and the new raster is created from an old raster based on the Value > 14. In my code, 1st for loop is for original raster and keeping the value and count fields in a dictionary. Then, creating new raster using Value > 14. The 2nd for loop is for the new raster and keeping its value and counts of in a new dictionary and finally calculating the percentage. However, I am getting an error on the 2nd for loop. The error says " 'Raster' object is not iterable". Anybody has any suggestion. Thank you for your help! Code:               import arcpy, os     from arcpy import env     from arcpy.sa import *          #To overwrite output     arcpy.env.overwriteOutput = True          #Set environment settings     env.workspace = "C:/Subhasis/Test/Neshanic_Python"          outws="C:/Subhasis/Test/Neshanic_Python/extract"          #checkout ArcGIS spatial analyst extension license     arcpy.CheckOutExtension("Spatial")          # set local variable     inraster = ["01367620-r-r"]          for i in inraster:         flds = ("VALUE", "COUNT")          dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(i, flds)}          sumcnt = sum(dct.values())         newraster = ExtractByAttributes(str(i), "VALUE>=14")         for j in newraster:             flds1 = ("VALUE", "COUNT")             dct1 = {row[0]:row[1] for row in arcpy.da.SearchCursor(j, flds)}             sumcnt1 = sum(dct1.values())             percentage=(sumcnt/sumcnt1)             print percentage


title: How to create input data for the LRS plugin from an excel table without coordinates?

text: I'm working with QGIS 2.2.0. I want to add a shapefile from an excel file that has informations I want to use in the LRS plugin. The excel file has no geographic information. Every time I try to save it as a shapefile, the program saves in a diferent format that I can not use in the LRS plugin. I import the file as a text layer, with no geographic information, than I save it as shapefile, but that doesn't happen. How can I do it? To use the LRS plugin, I have defined a road, made by a single line with a length. My Excel file has points with informations about specific kilometric points of that road, but don't have the coordenates, they have the length of the road matching that point. Do I have to add the respective coordenates? Thank you!


title: plotting a SpatialPolygons object using lines() without lines connecting separate polygons

text: This is almost an extension of a question asked on stackoverflow (http://stackoverflow.com/questions/12196440/extract-feature-coordinates-from- spatialpolygons-and-other-sp-classes), so I'll use the demo dataset from there. I'm attempting to plot these polygons with lines(), so I can overlay these with an image.plot. To do this, I use fortify, like so                   # taken from the linked stackoverflow question         Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))         Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))         Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))         Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)              Srs1 = Polygons(list(Sr1), "s1")         Srs2 = Polygons(list(Sr2), "s2")         Srs3 = Polygons(list(Sr3, Sr4), "s3/4")         SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)              image.plot(x = c(0:10), y = c(0:10), z = matrix(runif(100, 0,1), nrow = 10),             col = terrain.colors(20)) # fake data so lines() plays nice.         lines(fortify(SpP))         # there's a weird line connecting the square island in the middle to the polygon surrounding it      What I would like to get is an image like this, which I can produce with ggplot...                   require(ggplot2)         ggplot(aes(x = long, y = lat, group = group), data = fortify(SpP)) + geom_path()      I can see two possible solutions to this, either an analog to the "group =" argument in the aes() call in ggplot(), or some way to get fortify to insert 'NA' between different groups (this is what the lines() function interprets as breaks). I can produce an effect manually, but I couldn't do this with the real dataset.                   SpP3 <- rbind(SpP2[1:5,], rep(NA, 7), SpP2[6:10,], rep(NA, 7), SpP2[10:14,], rep(NA, 7), SpP2[15:19,])         image.plot(x = c(0:10), y = c(0:10), z = matrix(runif(100, 0,1), nrow = 10), col = terrain.colors(20))         lines(SpP3)      edit: I'd like to be able to produce these plots using image.plot, since it's faster for my real dataset, and I will be more compatible with what others I'm working with are using. Thanks!


title: What's the difference between vincenty and great circle distance calculations?

text: Python's geopy package features two distance measurements techniques: Great Circle and Vincenty's formulae.               >>> from geopy.distance import great_circle     >>> from geopy.distance import vincenty     >>> p1 = (31.8300167,35.0662833)     >>> p2 = (31.83,35.0708167)     >>> vincenty(p1, p2).meters     429.16765838976664     >>> great_circle(p3, p4).meters     428.4088367903001      What is the difference? Which distance measurement is preferred?


title: Export ArcMap cartographic representation to a symbol?

text: Using _Cartographic Representation_ (CR) I've created a nice feature symbol with a filled interior and sharp outline. It looks great on a map, on screen and printed. However I can't duplicate that symbol in the legend, as the regular Symbol Property Editor doesn't have the same level of control and edit options as the CR marker editor does. Here's the CR marker editor, and the resultant symbol: ![Representation Marker Editor result](http://i.stack.imgur.com/3To4V.png) And here's the closest I can get with regular symbol editor, just not good or close enough. With stacked level character symbol the two don't line up at all, and with halo effect the effect is thicker on the right side than the left. ![symbol editor, unsatisfactory results](http://i.stack.imgur.com/dOQeK.png) * * * So, how can I export the CR symbol so I can use it in the legend (and elsewhere)? We're not using live legend as it doesn't afford enough control over the result, for instance the size of the CR symbol can't be changed, even after converting to graphics.


title: Why does Raster Calculator (RasterCalculator_sa) give ERROR 000539 in ArcPy?

text: I am getting a not executed error for a raster calculator. Here is part of the output: > File "C:\projectTemp\FloodAnalysis\floodxsectionanalysis.py", line 107, in > makeSeemless arcpy.gp.RasterCalculator_sa(Con(IsNull(rasterIn), > FocalStatistics(rasterIn,NbrRectangle(3, 3), "MEAN"), rasterIn), > "{0}yrSeemless.img".format(increment)) File "C:\Program Files > (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing_base.py", line 498, in > return lambda *args: val(*gp_fixargs(args, True)) > arcgisscripting.ExecuteError: ERROR 000539: Parsing error SyntaxError: > invalid syntax (line 5) Failed to execute (RasterCalculator). This code I have setup to put into this calculator is:               def makeSeemless(outputPath, increment, rasterIn):         arcpy.gp.RasterCalculator_sa(Con(IsNull(rasterIn), FocalStatistics(rasterIn,NbrRectangle(3, 3), "MEAN"), rasterIn), "{0}yrSeemless.img".format(increment))


title: Where can I download OpenStreetMap data by geographic region?

text: I want to use OpenStreetMap data for an ArcEngine application I'm working on, and I need the data to be stored locally. I don't want to download the entire planet.osm file, but rather by state, county, region, etc. **Are there any resources where I can download OSM data by geographical region?**


title: Is there a way to set a mask(halo) for labels in a python script?

text: I currently am writing a python script and using a lot of the arcpy mapping module functionality. I have been using the LabelClass of the arcpy mapping module to automate the display of labels with my ptyhon script. Is there a way to set a mask(halo) for labels with arcpy?


title: How to buffer raster pixels by their values?

text: The pixels to the left represent tree locations and their associated crown radii (i.e. pixel values ranging from 2 - 5). I would like to buffer these raster pixels by their crown radius value. The image to the right is what I am hoping to accomplish _using only raster processing methods_. I would initially think to use a circular focal sum in ArcGIS, although the neighborhood setting is a fixed value, which would not take into account the variable sized crown radius. What is a good method to "buffer" pixels by their values? ![enter image description here](http://i.stack.imgur.com/rb33F.jpg)


title: How to Improve a Distance Model

text: I am working on a distance model that creates distance rings around a given point and assigns those rings a probability. The point is in the center of a 25 x 25 km area, and the model fails to assign a probability to the area that is outside of the rings but inside the 25 x 25 km area. How can I build a simple model in model-builder where I can merge the two layers together and assign probability to the area that is not included in the model? The probability that I want to assign to the outside area is 5%. Thank you


title: Intersection between a ray from lat, lon, alt, with pan, tilt and earth's surface

text: I have a problem that consists of finding the intersection point between a ray from a given lat, lon, alt, that is looking to a given pan and tilt angles and the earth's surface. I need the lat, lon point where this ray intersects with the earth surface. Pan is a horizontal clockwise angle from true north, and tilt is a vertical angle: 0 pointing to earth's center, 90 to horizon, 180 sky. Depending on the given altitude, and the pan and tilt values, the ray may not intersect in any point with the earth. Is there an easy way to compute this intersection?


title: How to get tiff file of project limits at specific scale using QGIS 2.0?

text: I have shapefile which has a layer extent scale of 1:551920. I want the TIF file exports of the file at a scale of 1:200000. Is there way to do it automatically? Or do I have to do it manually by moving the screen? I would like to use QGIS 2.0 for my problem.


title: How to solve Geoserver Labeling issue

text: In google map all the labels are rendering well and there are none missing. In my map some times only the point marker is showing and the label is not displaying. I am asking about "Auto align and AUTO OFFSET" OPTION in geoserver like mapserver. In this image all the labels are shown in different offset. ![enter image description here](http://i.stack.imgur.com/9G40M.jpg) What do I have to change in my SLD?


title: ArcObjects ExportMapImage modifies image extent

text: I am trying to extract a map image at a certain range. I am using IMxdServer's ExportMapImage method, documented here. Here is the relevant snippet of code:                   IMapDocumentPtr ipMapDoc(CLSID_MapDocument);         IMapDescriptorPtr ipMapDescriptor(CLSID_MapDescriptor);         IMxdServerPtr ipMapServer(CLSID_MxdServer);         IArrayPtr ipArray(CLSID_Array);         IUnknownPtr ipUnknown;         int dpi = 300;                 IEnvelopePtr ipPixelBoundsEnv(CLSID_Envelope), ipReturnExtent(CLSID_Envelope);         ipPixelBoundsEnv->PutCoords(mUpperLeftLong, mLowerRightLat, mLowerRightLong,  mUpperLeftLat);              ipMapServer->Start(MxdPath);         ipMapServer->get_MapDescriptors(&ipArray);              ipArray->get_Element(0, &ipUnknown);         ipMapDescriptor = (IMapDescriptorPtr) ipUnknown;              IColorPtr ipColor(CLSID_RgbColor);         ipMapDescriptor->get_BackgroundColor(&ipColor);         ipMapDescriptor->put_TransparentColor(ipColor);         ipMapDescriptor->put_Extent(ipPixelBoundsEnv);         ipMapDescriptor->put_FullExtent(ipPixelBoundsEnv);         ipMapDescriptor->put_KeepAspectRatio(VARIANT_FALSE);              HRESULT hr = ipMapServer->ExportMapImage(ipMapDescriptor, dpi, (256*(hortizonalTileBlockCount+2)), (256*(verticalTileBlockCount+2)), CComBSTR("PNG8"), bstrFileLocation, &ipReturnExtent);         double ux, uy, lx, ly;         ipReturnExtent->get_XMax(&ux);         ipReturnExtent->get_XMin(&lx);         ipReturnExtent->get_YMax(&uy);         ipReturnExtent->get_YMin(&ly);         std::cout<<"return extent "<<lx<<","<<ly<<","<<ux<<","<<uy<< " \n";      In my current example the following bounds are passed in:               -75.9375,39.6395,-77.3438,40.714      The return extent is shifted a touch:               -77.3438,39.4736,-75.9375,40.8799      I'm extracting quite a few images but I put one example extent for reference.


title: How to change a schema name within a MXD with ArcObject .NET?

text: I have a MXD file where all layers point to an Oracle Spatial database. Recently our DBA changed the schema name. So now my MXD file has broken links, each time I open it, all my layers have a exclamation mark beside them. In order to bulk change the new schema name in the MXD file, I've written an .NET ArcObjects application. But since all my layer data links are broken in my MXD, I do not have the capability to get a reference on IFeatureClass from a IFeatureLayer. I've found on the web that is normal if a IFeatureLayer is broken within the MXD that the IFeatureClass always returns null. Actually, all MXD's layers point on the schema "A" and now I want them to point on schema "B". How to do that programmatically? Can you help me to get the good object that should give me access to change the QueryDef of my layers? Here is the code that I've written until now:                        //Exit if no map document is selected             string sFilePath = openFileDialog1.FileName;             if (sFilePath == "") return;                  IMapDocument ipMapDoc;             ipMapDoc = new MapDocumentClass();                       ipMapDoc.Open(sFilePath);                  IMap map = ipMapDoc.get_Map(0);                  IEnumLayer layers = map.get_Layers();             ILayer layer;             while ((layer = layers.Next()) != null)             {                 if (layer is IFeatureLayer)                 {                     IFeatureLayer featureLayer = layer as IFeatureLayer;                     IFeatureClass featureClass = featureLayer as IFeatureClass;                          if (featureClass == null)                      {                         Console.WriteLine("FeatureClass where FeatureLayer name is " + layer.Name + " is null");                     }                      }             }


title: Assigning a changing value to a variable in ArcMap using Python

text: I have two lists/strings/sets (the exact format can easily be changed) of values (u and v). I would like to do a cost distance analysis between each point in u and all points in v. For this, I need to assign a value to the source and destination variables (n and m). I would like to do this using an automatically revolving value for n and m based on the values in u and v, because of the size of my set (4000 entries). I cannot, however, find the code which would allow me to do this. So far, I have u and v, which are made using this:               fc=villagesnear     field='IN_FID'     values = [row[0] for row in arcpy.da.SearchCursor(fc,(field))]     u=str(values)     print(u) ##prints a list of the ID  of every village with another within 500m          fc=villagesnear     field='NEAR_FID'     values = [row[0] for row in arcpy.da.SearchCursor(fc,(field))]     v=str(values)     print(v) ##essentially prints the same list of numbers as IN_FID, but to be used as destination data      and n and m need to be defined to then be used in the following code:               costpath=arcpy.sa.CostPath('n', 'costD', 'costB','EACH_ZONE','m')     cursor=arcpy.da.InsertCursor('output','n','m')     threshhold_cost=1000     while (n<=4909):         if n == m:             skip         if 'costpath'<='threshhold_cost':             cursor.insertRow([n],[m])             del cursor         else:             print "done"             n+=1      any help would be much appreciated!!


title: How to set a fixed scale to a layer in ArcMap?

text: I have a polygon layer in ArcMap. To full extent the scale is 1:35,000. It only fits in a 210cm x 160cm page. I want to fix the scale to 1:10,000 and make it fit to a 84.1cm x 59.4. Is there a way to do it in ArcGIS 10.1?


title: How to convert MXD to DWG File

text: I have been trying to export from MXD format to DWG format without much success. I used the Export to CAD conversion tool in GIS 9.3 but it only converts the lines and polygon.. etc to DWG format and NOT the colors and texts that are also included in the original MXD file. Any idea of how I might be able to convert the WHOLE MXD to CAD with all the features included?


title: [ArcGIS toolbars.Draw.POLYGON ]Object #<BarProp> has no method 'activate'

text: I am facing this error while activating the toolbar. This is my body code                 <button dojoType="dijit.form.Button" onClick="toolbar.activate(esri.toolbars.Draw.POLYGON); map.hideZoomSlider();">Polygon</button>      And this is the javascript i am loading:                //Polygon & FreeHand Polygon------------------------------------------------------------------------------------                var toolbar, symbol, geomTask;                  function createToolbar(map) {                 toolbar = new esri.toolbars.Draw(map);                 dojo.connect(toolbar, "onDrawEnd", addToMap);             }              dojo.require("esri.tasks.geometry");         dojo.require("esri.toolbars.draw");                function addToMap(geometry) {                      switch (geometry.type) {                          case "polygon":                         var symbol = new esri.symbol.SimpleFillSymbol                         (esri.symbol.SimpleFillSymbol.STYLE_SOLID, new esri.symbol.SimpleLineSymbol                         (esri.symbol.SimpleLineSymbol.STYLE_DASHDOT, new dojo.Color([255, 0, 0]), 2), new dojo.Color([255, 255, 0, 0.25]));                              break;                      }                 var graphic = new esri.Graphic(geometry, symbol);                 map.graphics.add(graphic);                      var theme = document.getElementById('SelectTheme').value;                 var token = "xkg8VRu6Ol+gMH+SUamkRIEB7fKzhwMvfMo/2U8UJcFhdvR4yN1GutmUIA3A6r3LDhot215OVVkZvNRzjl28TNUZgYFSswOi";                 var queryTask = new esri.tasks.QueryTask("http://www.onemap.sg/DataService/Services.svc/" + theme + "?token=" + token);                      var query = new esri.tasks.Query();                      query.geometry = graphic.geometry;                 query.outFields = ["DESCRIPTION"];                 query.geometryType = "";                 query.outSR = "";                 query.returnGeometry = true;                 query.spatialRel = "esriSpatialRelIntersects";                 query.text = "";                 query.where = "";                      //Execute QueryTask                 queryTask.execute(query);                 dojo.connect(queryTask, "onComplete", function (fset) {                     //create symbol for selected features                               //                var symbol = new esri.symbol.SimpleMarkerSymbol();                     //                symbol.style = esri.symbol.SimpleMarkerSymbol.STYLE_CIRCLE;                     if (theme == "Disability") {                              var symbol = new esri.symbol.PictureMarkerSymbol('/images/disability.jpg', 25, 25);                              symbol.setColor(new dojo.Color([255, 255, 255, 1]));                     }                     else if (theme == "VoluntaryWelfareOrgs") {                              var symbol = new esri.symbol.PictureMarkerSymbol('/images/Welfare.jpg', 25, 25);                              symbol.setColor(new dojo.Color([1, 1, 1, 1]));                     }                      else if (theme == "ElderCare") {                              var symbol = new esri.symbol.PictureMarkerSymbol('/images/elder.jpg', 25, 25);                              symbol.setColor(new dojo.Color([1, 1, 1, 1]));                     }                     else if (theme == "CommuinityClubs") {                              var symbol = new esri.symbol.PictureMarkerSymbol('/images/ramp.jpg', 25, 25);                              symbol.setColor(new dojo.Color([255, 255, 255, 255]));                     }                          var infoTemplate = new esri.InfoTemplate();                     infoTemplate.setTitle(theme);                     var resultFeatures = fset.features;                     for (var i = 0, il = resultFeatures.length; i < il; i++) {                         var graphic = resultFeatures[i];                         graphic.setSymbol(symbol);                         graphic.setInfoTemplate(infoTemplate);                         OneMap.map.graphics.add(graphic);                     }                 }                      );             }      can anyone help me out?


What is currently the highest value of OSM_ID in openstreetmapping? What is the highest possible value of OSM_ID? 2^31 (2.147.483.648)? And here I add some more text or my question is rejected because it is too short...


title: Arcpy Script tool

text: I'm creating a script tool and using the SelectLayerBy Attribute function (below), arcpy.SelectLayerByAttribute_management("cities_lyr", "NEW_SELECTION", '"CRIME_INDX" <= 0.02') However, I keep getting this error - ERROR 000358: Invalid expression. Can anyone please offer some clarity here. Thanks.


title: How to link Geoserver WMS layer to Google earth?

text: I tried with what was available in Geoserver Layer preview as kml. It gives the geom itself at higher zooms, not the image for that area. What to do to get only the png/jpeg images from Geoserver and show it in Google Earth,over normal sat. view?


title: can i merge resolution with float single data type

text: Can I merge resolution with both image is float single data type? I have done this, but the output image become white. Do the data type must be in unsigned-8 bit for the merge to be success?


title: Looking for non-Earth GIS Web Services

text: Hey I'm building a web application and want to use some GIS services for other planets/moons. In general, I want to find GIS web resources which apply to bodies apart from Earth. Anybody have url endpoints to a:   * Mars tile service    * Titan WMS service   * Jupiter's Red Spot place finder   * Milky Way routing service   * Pluto geocoder I'm more interested in services on the web, rather than data I would need to download and tile/host myself. Thanks in advance for all answers!


title: Using ArcGIS COGO (Coordinate Geometry)

text: I have used COGO to create two polygons in a proposed outgrant file for some land parcels. The polygons are correctly shaped and in the correct area but they are spatially skewed and oversized. How do I fix the projection or coordinate system, IF that is the problem? Currently the data frame is set to NAD83-UTM Zone 14N.


title: Calculating field length in certain directions (raster)

text: Currently I'm working on winderosion. For my purpose I like to calculate the length of agricultural fields in certain wind direction. I would like to do it for raster. To be more precise: I have agricultural areas as SHP and can convert them to raster. After it I'll calculate the distance of one boarder of the area to another in a certain direction (for instance 270° West). At the end I'll weight the different rasters (with field length in different directions) to get a mean feld length. So my idea was to designate those cells within the area "1", those outside "0" and use the raster calculator to start calculating with the first "1" and add 20 (m = the resolution of one cell). The next cell should take the value of the previous one plus 20 and so on until it reaches "0" again. This is just the idea but I have no idea how to do it in the raster calculator of ArcGIS/QGIS/SAGA and use a certain direction for calculating. Have you any proposals for me? I would be also very thankful if you give me a complete different way of soultion besides my ideas. I also thought about neighborhood or zonal statistics.


title: Using ArcObjects control commands without ToolbarControl

text: I want to use some ArcObjects predefined commands for my buttons actions..I used this code               `ICommand command = new ControlsSceneZoomInToolClass();                 command.OnCreate(axMapControl1.Object);                  axMapControl1.CurrentTool = (ITool)command;`      when I press button the muse chenges to zoom in but it doesn't work..Is there some thing wrong with my code?


title: Creating surface from points in ArcGIS for Desktop?

text: how to create surface from points in ArcGis? I have a points (centroids of the neighbourhood and number of people living in every neighbourhood) and i would like to create surface that would present density of population.


title: Print custom feature form using PyQGIS

text: I am looking for a way to print a custom feature form created in PyQt Designer for QGIS 2.4 using PyQGIS. My goal is to print/save image as the current extent of the map with the selected feature's attributes of interest on the same page. I'm experimenting with two different ways: First, is to use QgsComposition and QPrinter to add and print items to a pdf: including the map, legend, and attributes. My issue here is getting the only the desired attributes to display dynamically changing depending on the active layer. If I could add my feature form as a map item, I would be displaying the attributes added to the custom form. UPDATE: The following code may be able to display my map with a legend and selective table, but I still have many issues with the map items (no text in the table, positioning, unable to include only desired layer in legend). I've included all the code associated with this print composer:               layer = qgis.utils.iface.activeLayer()          qgis.utils.iface.actionZoomToSelected().trigger()     qgis.utils.iface.mapCanvas().zoomScale(1000)          mapRenderer = iface.mapCanvas().mapRenderer()     c = QgsComposition(mapRenderer)     c.setPlotStyle(QgsComposition.Print)     x, y = 0, 0     w, h = c.paperWidth(), c.paperHeight()     composerMap = QgsComposerMap(c, x,y,w*.75,h)     c.addItem(composerMap)          legend = QgsComposerLegend(c)     c.addComposerLegend(legend)     legend.model().setLayerSet(mapRenderer.layerSet())     legend.setTitle('')          legendModel = QgsLegendModel(c)     legendModel.setLayerSet([layer.id()])          legend.updateLegend()          selected_ids = [f.id() for f in layer.selectedFeatures()]     filter_id_string = ','.join([str(id) for id in selected_ids])          table = QgsComposerAttributeTable(c)     table.setVectorLayer(layer)     table.setMaximumNumberOfFeatures(20)     table.setFeatureFilter("id in (" + filter_id_string + ")")     table.setFilterFeatures(True)     table.setItemPosition(x+50, y+50)     table.applyDefaultSize()     col1 = QgsComposerTableColumn()     col1.setAttribute('Diametre')     col2 = QgsComposerTableColumn()     col2.setAttribute('Nature')     table.setColumns([col1, col2])     c.addItem(table)          printer = QPrinter()     printer.setOutputFormat(QPrinter.PdfFormat)     printer.setOutputFileName('out.pdf')     printer.setPaperSize(QSizeF(c.paperWidth(), c.paperHeight()), QPrinter.Millimeter)     printer.setFullPage(True)     printer.setColorMode(QPrinter.Color)     printer.setResolution(c.printResolution())     pdfPainter = QPainter(printer)     paperRectMM = printer.pageRect(QPrinter.Millimeter)     paperRectPixel = printer.pageRect(QPrinter.DevicePixel)     c.render(pdfPainter, paperRectPixel, paperRectMM)     pdfPainter.end()      Second, is to embed a image of the current extent (using mapCanvas().saveasImage()) using a dynamic QPixMap element on the custom form. This brings me back to the same question as to whether or not the custom feature form could be printed (to pdf, and eventually paper).


title: Python - OGR: Distance different than in ArcGIS

text: I want to determine the shortest distance between a point and a line. I use the following Python script to do so. If I compare the result of the script with the result of ArcMap's Near Analysis, I get a completely different results. The shapefiles are in the CRS GCS_North_American_1983 and projected with Albers. What could be the reason for the difference?               import ogr          driver = ogr.GetDriverByName('ESRI Shapefile')          lineshp = driver.Open('U:/My Documents/Tool/shp/line.shp', 0)     linelyr = lineshp.GetLayer()          pointshp = driver.Open('U:/My Documents/Tool/shp/point.shp', 0)     pointlyr = pointshp.GetLayer()          line = linelyr.GetNextFeature()     point = pointlyr.GetNextFeature()          line_geom = line.GetGeometryRef()     point_geom = point.GetGeometryRef()          dist = point_geom.Distance(line_geom)          print dist


title: Distance Between Two Points Same Value Different Feature Classes

text: I have two point feature classes - one is the original data and one is the data after editing. I want to know how far each of the points moved in the editing process. How do I compare all of the points in the new data to the points in the old data that have the same ID# and display the distance between them? Also I have the standard license, so I don't have access to the Near or Generate Near Table tools.


title: Qgis - unwanted offset when using raster interpolation plugin

text: I'm using QGis and I have a set of points with a depth attribute. I'm trying to make an interpolation using the Raster interpolation plugin (Raster -> interpolation -> interpolation). But when I do so I get a misplacement on the produced raster. It's positioned too far north. Do you have any ideas of what may be causing this? To make it clear what I am doing and what my problem is I've created a short video where I'm showing the issue: http://jimutt.com/interpolationProblem.swf As you can see in the video I'm using the same projection in both the viewport and for the point layer. Please tell me if you can't display the video.


title: Display labels over a graphic (ArcObjects)

text: I've created a graphic polygon that is actually a polygon with the center cut out (I created it using ITopologicalOperator.Difference to remove the smaller polygon from the larger one). Now I want the labels on my map to display in the cutout area of my polygon. I found a suggestion to set IBarrierProps.Weight on my graphics layer and tried setting the weight to 0 hoping that the labels would still display but no luck. Any suggestions for me?


title: Load .shp file in Hive(Hadoop)

text: I need to load .shp files in Hive more precisely in HortonWorks sandbox through Beeswax, can anyone provide example of Hive commands or steps for Beeswax ?


title: Setting Rotation of ExportWebMap

text: I see from the ExportWebMap specification that rotation can be set in the JSON object sent to the sever. Does anyone know if this can be accessed via the ESRI JavaScript API (v3.9)? I do not need to display the map with a rotation, just honor it in a printout Thanks, Joe


title: GeoDjango / PostGIS - Cut geometry by bounding box

text: I am looking for a solution to retrieve the geometry(ies) from a polygon/multipolygon within a bounding box in a GeoDjango application. The application I am developing has a requirement to serve out tiles of geographical data, and many of our geometries are very large and span across multiple tiles. To reduce the payload of each tile I would like to splice(cut) a geometry shape(s) from a bounding box. I am storing Planet.OSM data in a PostGIS database. Any suggestions would be greatly appreciated! Thanks


title: Label using VBscript in ArcMap - multiple IF statements?

text: I have 4 fields that hold number values. I only want a label returned if the field <>0. For instance, if col1<>0, then label = col1 AND if col2 <>0 then label = col2, etc. So, I want all four columns to show up as labels, but only if they are not equal to 0. I have tried finding info online and am having issues. Thought maybe everyone could help. The following script works, but only for one at a time. I need to combine them so that this script works with Seq2, Seq3 and Seq4. I am thinking there is a way to do multiple IF statements?               Function FindLabel ( [Seq1_TractDepthSequence], [Seq1_BegInterval] ,[Seq1_EndInterval], [Seq1_DeckCorpRITot], [Seq1_DeckCorpORRITot] , [Seq1_DeckCorpIntTot] , [Seq1_DeckCorpGWITot]   )     if ( [Seq1_TractDepthSequence] = 1) AND ([Seq1_DeckCorpRITot] <> 0) then       FindLabel = [Seq1_BegInterval] &" - "& [Seq1_EndInterval] &" : "& [Seq1_DeckCorpRITot] &" RI"     end if     End Function      The following won't work:               Function FindLabel ( [_Seq1_TractDepthSequence], [_Seq1_BegInterval] ,[_Seq1_EndInterval], [_Seq1_DeckCorpRITot], [_Seq1_DeckCorpORRITot] , [_Seq1_DeckCorpIntTot] , [_Seq1_DeckCorpGWITot]   )     if ( [_Seq1_TractDepthSequence] = 1) AND ([_Seq1_DeckCorpRITot] <> 0) then       FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpRITot] &" RI"     if ( [_Seq1_TractDepthSequence] = 1) AND ( [_Seq1_DeckCorpORRITot]  <> 0) then       FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpORRITot] &" ORRI     if ( [_Seq1_TractDepthSequence] = 1) AND ( [_Seq1_DeckCorpIntTot]  <> 0) then       FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpIntTot] &" CorpInt"     if ( [_Seq1_TractDepthSequence] = 1) AND ( [_Seq1_DeckCorpGWITot]  <> 0) then       FindLabel = [_Seq1_BegInterval] &" - "& [_Seq1_EndInterval] &" : "& [_Seq1_DeckCorpGWITot] &" WI"     else     end if     End Function