title: Colour matching satellite imagery from different time periods

text: Is it possible to colour match two separate satellite images taken at two different times? I am working with two high resolution, multispectral images that were taken approximately one year apart and I want to do some change analysis with them. Unfortunately one of the images is extremely washed out in urban areas and I'm trying to correct it (see attached samples). **Image From Time 1** ![Sample Time 1](http://i.stack.imgur.com/cENsF.jpg) **Image From Time 2** ![Sample Time 2](http://i.stack.imgur.com/RuCq4.jpg) Ideally, I would like to have both of the images to appear as similar as possible, preferably like the first one, so I can identify changes between the two time periods. I am using ArcGIS 10.2 Professional, but if there are other options to deal with this I would like to hear about those as well.


title: Using reader Feature Type Name as Writer Feature Type Name in FME

text: I've search and searched for a way to do this, but no method seems to work. Basically I want to use the input feature class name (feature type in FME) as the output name. This normally works, however the tool I'm using (ChangeDetector) adopts the "revised" input as the output Feature name instead of the "original" input name. Even linking the "original" input name to the output name attribute does not work. It keeps the "revised" input name. I've tried using various transformers like the SubstringExtractor to pull the name and set the writer to fan out by this attribute, but even that fails. While I would prefer to get this working within the workbench, I do plan on using this workbench in a python script, so it would be handy to get it working through that as well. I setup the workbench to run using               subprocess.Popen(["pathToFme", "pathToWorkbench"])      This works just fine. Subsequently I tried adding parameters at the end, to input the name based on code samples I found while searching this...               subprocess.Popen(["pathToFme", "pathToWorkbench"] --fileName: "nameOfOutput")      with "fileName" being the name of the parameter within the workbench, but I can't seem to get it to run, let along pick up the parameter input from the script. I couldn't get os.system to properly execute the workbench, so I want to stay away from that. Basically what I am trying to achieve is the ability to take a large number of inputs, run them through the changeDectector, and output with the same name as the "original" input. Any help or suggestions would be appreciated. Thanks.


title: What are resources to learn ArcGIS for someone who already knew Python

text: I know that there are several resources recommendation questions, however, they are all from GIS people trying to learn Python (for example, here). My situation is the opposite: I have worked extensively in Python, but know nothing about ArcGIS. What is the best way to learn? To be more specific, I want to learn the following:   * Create maps from shapefiles   * Incorporate other social science-y data on map. For example, conflict site, social network, economic activity, etc.   * Most importantly, statistical modeling with a spatial component (network, hotspot analysis, etc.)


title: Graphically depict statistical difference between thematic rasters

text: The graph shows a regression between two raster datasets with classified tree canopy cover: class 1 = 0-10%, class 2 = 10 - 20% ...class 6 = >60% (Figure 1, Figure 2). The spatial resolution for raster 1 is 1m while the spatial resolution for raster 2 is 30m. Note that raster 2 was resampled to 1m for analysis purposes in the script (Appendix A). I am looking for an elegant, graphical approach to statistically describe the difference between two thematic rasters at different spatial scales. My first attempt was to 1) create a raster stack, 2) randomly sample the stack, 3) run a regression and descriptive statistics, 4) graph the results (Appendix A). What other graphical methods exist for statistically comparing thematic raster data? I usually use R, Python and MATLAB, although I would be happy with generic solutions too. * * * **Figure 1** ![enter image description here](http://i.stack.imgur.com/m7Ceq.jpg) * * * **Figure 2** ![enter image description here](http://i.stack.imgur.com/ITbCt.png) * * * **Appendix A**               require(raster)     require(spatstat)     require(ggplot2)     require(gridExtra)     require(hexbin)     require(rgdal)          # Read the TIF data     file = 'C:/path/raster1.tif'       file2 = 'C:/path/raster2.tif'                    # Create raster objects from the single band canopy products     r = raster(file)     r2 = raster(file2)          # Resample r2 to 1m for analysis purposes     r2 = resample(r2, r, method = "ngb")          # Stack the two raster layers for analysis     raster <- stack(r, r2)          # sample random locations in raster stack and report values in a dataframe     df = data.frame(sampleRandom(raster, size=1000, cells=TRUE, sp=TRUE))          ## Do the regression and plot the results     new = df$X3711203_ne     old = df$X3711203_ne_30m          # Calculate RMSE and other values     fit <- lm(old ~ new)     rmse <- round(sqrt(mean(resid(fit)^2)), 2)     coefs <- coef(fit)     b0 <- round(coefs[1], 2)     b1 <- round(coefs[2],2)     r <- round(sqrt(summary(fit)$r.squared), 2)          # Build equation see ?plotmath     eqn <- bquote(italic(y) == .(b0) + .(b1)*                     italic(x) * "," ~~ r == .(r) * "," ~~ RMSE == .(rmse))          eqn_text = as.character(as.expression(eqn))          # Get max range to automatically update plot size     li = c(range(new)[2], range(old)[2])     range = tail(sort(li),1) + 0.05          # Plot results     p1 = ggplot(df, aes(new, old))  +        geom_point() +       geom_smooth(method = lm) +       stat_sum(aes(size = ..n..)) +       xlab("Cover Class (1m)") +       ylab("Cover Class (30m)") +       annotate("text", x= 0, y= range, label=eqn_text, hjust=0, size=4, face="italic", parse=TRUE) +       xlim(0, range) +        ylim(0, range) +       ggtitle("Tree Canopy Cover Class")          p1


title: How to export query result on raster table in PostGIS

text: I did success query on a raster table in postGIS(raster data clipped by a vector polygon), now I want to view the query result in QGIS. But how to export the query result to another table. In vector data, it's simple, just "CREATE TABLE schema.someTable AS (...sql in here...); How about raster data? Any suggestion? Thanks.


title: I want to get rid of unnecessary map and data from shp file, how?

text: I have loaded the base map from esri's defualt provided maps but its for the whole world, and i have gotten the roads data too, but its for whole continent. I only need for one city i am going to be working on to save my RAM. I am using arcgis 10.2


title: Unable to load PostGIS vector layer into QGIS using python cookbook recipe

text: I'm having trouble loading a PostGIS layer using Dufour (from osgeo4w), using the guide at http://www.qgis.org/en/docs/pyqgis_developer_cookbook/loadlayer.html. I just get the 'invalid layer' message. My parameters appear to be perfect.               from PyQt4 import QtGui     from PyQt4 import QtCore     from qgis.core import *     from qgis.gui import *     from ui_townplanning import Ui_TownPlanning     # create the dialog for zoom to point               class TownPlanningDialog(QtGui.QDialog):         def __init__(self, iface):             QtGui.QDialog.__init__(self)             # Set up the user interface from Designer.             self.ui = Ui_TownPlanning()             self.ui.setupUi(self)             self.iface = iface             self.canvas = self.iface.mapCanvas()             self.addPostGISLayer('myhost', '5432', 'glis', 'user', 'pass',                                  'tpadmin', 'planning_application_boundary_validation',                                  'wkb_geometry')              def addPostGISLayer(self, host, port, dbname, username, password,                             schema, table, geom_col):             '''Adds a PostGIS table to the map'''             uri = QgsDataSourceURI()             uri.setConnection(str(host), str(port), str(dbname), str(username), str(password))             uri.setDataSource(str(schema), str(table), str(geom_col), '', str(key_col))             vlayer = QgsVectorLayer(uri.uri(), str('test'), 'postgres')             # QgsMapLayerRegistry.instance().addMapLayer(vlayer)             if not vlayer.isValid():                 msgbox('Layer not loaded', uri.host(),                        uri.database(), uri.port(), uri.username(), uri.password(),                        uri.schema(), uri.table(), uri.geometryColumn())          def msgbox(*args):         a = [str(i) for i in args]         text = '\n'.join(a)         qmsg = QtGui.QMessageBox()         qmsg.setText(str(text))         qmsg.exec_()      **Code now works after correction**


title: How to return geo-tagged social media posts?

text: I have seen quite a few ESRI web apps that utilize Geo-Tagged social media posts as a point data layer. I'm wondering how one goes about querying things like Twitter, Flickr and YouTube to return posts related to a specific location. I'm also hoping to be able to filter these results based on hash- tags. Twitter is my main interest, but I know when I get that working I'll have people clamoring for the others. My hope is to add this to our public web viewer where citizens view pictures, videos and tweets about our public lands. Also, are these Geo-Tagged posts limited by your own personal followers / friends? Here is an ESRI example of what I am looking for: ESRI Australia Wildfires Map


title: Can FME Desktop (node-locked) be installed on an Amazon EC2 or VPC instance?

text: I did do some searching in the FME help and the legal notices, but I'm not getting any hits (probably due to vocabulary). Are there any legal restrictions for FME Desktop that would prevent me from installing onto an EC2 or VPC Amazon instance? I'm not sure if (or how much) that would compete against the FME Cloud offering. I'm mostly interested for test and evaluation, but I am also curious if anyone is using an EC2/VPC as a permanent home for FME Desktop. This another question that is inspired by this closely related question:   Is there an optimum RAM configuration to use for FME Desktop?


title: ArcMap 10.0 won't load

text: ArcMap 10.0 on a Windows 7, 32-bit computer stopped loading two days ago. It worked perfectly before then. ArcCatalog, etc., all continue to open fine. Task Manager indicates ArcMap is trying to load, but very, very slowly. I deleted Normal.mxt, and tried to restart ArcMap, but it has not yet successfully opened up and has not re-created Normal.mxt.


title: Error opening feature class from gdb

text: I have a file geodatabase that a person gave to me for work. I try to open the feature classes in Arcmap 10.2 and some of it opens perfectly but with the others appears a message that says: Error opening feature class. Sometimes when I try to open the properties of the layers I get "Failed to edit the selected object(s)." What can I do in this case? ![enter image description here](http://i.stack.imgur.com/uOdzG.png)


title: How to start a local GIS Users Group

text: My immediate geographic area is lacking an organized local GIS Users Group. I belong to a group that meets a couple times a year, but its about 1.5 hour drive. I'm thinking about starting up a Users Group for our County, we have about 50 users just in my organization that I think could benefit, not to mention the other users from different organizations. Problem is, I've never done anything like this and don't exactly know where to begin. Does anybody have tips or pointers for how to start a group or how to make a GIS group successful?


title: How to export Symbology from ArcGIS to Microstation or QGIS?

text: I am using ArcGIS 10.1 and trying to figure out how to export my symbols of elements in ArcGIS to CAD or QGis readable format. Problem is, that I can only export symbology to Style, or export map to DGN, but there will be export just only graphics without symobology. Please, do you know how to solve this? Thanks a lot!


title: Where is Google Maps extension for ArcGIS?

text: I need to import Google Maps into ArcGIS for Desktop and apparently I need to activate a Google Maps extension. I'm working on ArcGIS 10.0 and there are no such extension. Do I have to install a new extension?


title: Utility GIS - Type of data cleanup needed & Effort

text: I am trying to get into Utility GIS and I have been asked by my boss to size up the effort of what it takes to clean up Pipeline, Utility asset GIS data and figure out if there is a way think of a factor like a per linear mile man days of effort. The data is in the form of ESRI File G Geodatabase. Typical clean up related tasks I have are : 1\. Data consistency by using domains (Manual Effort) - Per Record Multiplier 2\. Topology checks (ArcToolbox) Per Featureclass Multiplier 3\. CAD to GIS conversion (Per Feature Class Multiplier) 4\. Visual QCs Can the hive mind suggest any other tasks typically for cleaning pipeline data and any assessment of what you would think of a man-day multiplier. Appreciate the help!


title: what raster format for r.los in qgis

text: What raster format should I use to create viewsheds in QGIS using the grass r.los plugin? I am using the Ordnance Survey Opendata DTM (tried both ascii and shapefile formats) and srtm geotiffs but when I run r.los I can't select an input raster - does that mean the raster is the wrong format?


title: Clean way to install gdaltools plugin for QGIS on OSX

text: I have qgis 1.8.0 on osx mountain lion and am installing the gdaltools 1.2.29. I have run into a problem that a number of people are seeing. The plugin cannot find my gdal python (see below). Anyway, I can fix the problem hackily (if that is a word) by symlinking                 cd /Library/Python/2.7       sudo mv site-packages/ site-packages-old       sudo ln -s /Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site-packages /Library/Python/2.7      So that fixes the prob because the python path list that the plugin knows about includes /Library/Python/2.7/site-packages But the real solution is to add a path to the list of paths that the python plugin[s] know about. How do you do that? \----start of error message---- Couldn't load plugin GdalTools due an error when calling its classFactory() method Traceback (most recent call last): File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 164, in startPlugin plugins[packageName] = package.classFactory(iface) File "/Users/aberezin/.qgis//python/plugins/GdalTools/ **init**.py", line 32, in classFactory from GdalTools import GdalTools File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 309, in _import mod = _builtin_import(name, globals, locals, fromlist, level) File "/Users/aberezin/.qgis//python/plugins/GdalTools/GdalTools.py", line 39, in raise ImportError( error_str ) ImportError: No module named osgeo [python-gdal] Python version: 2.7.2 (default, Jun 20 2012, 16:23:33) [GCC 4.2.1 Compatible Apple Clang 4.0 (tags/Apple/clang-418.0.60)] QGIS version: 1.8.0-Lisboa Lisboa, Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python', '/Users/aberezin/.qgis//python', '/Users/aberezin/.qgis//python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat- darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat- mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat- mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib- tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib- old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib- dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC']


title: How to edit Attributes on WFS Layer through Popup?

text: I am working on a project for university. So far I have displayed WfsLayers which can be select from a panel, also I can click on my points whereupon a Popup will show up with some information which is derived from my Geoserver which is connected to my Post GIS table. So far so good. I now would love to be able to edit my attributes in the popup box. I found this example (http://gis.ibbeck.de/ginfo/apps/OLExamples/OL26/examples/styles_unique_with_group_wfs.html), which looks just great and is exactly what I am looking for but being a coding newbie in general I have terrible difficulties understanding the code. I tried to use parts but just can't get it to work at all. Does someone have an easier example I could try to work with?               <!DOCTYPE HTML>               <title>OpenLayers Demo</title>


title: Draw middle ordinate perpendicular to chord, out to curve, with arcpy

text: I've been assigned to find the radius of road curves. My boss found this article describing a pretty straightforward method of doing so:   1. find the start/end points of the curve;   2. create a chord;   3. create the middle ordinate (from chord midpoint, perpendicular line that goes up to the curve)   4. do a bit of math with the chord and ordinate lengths to find the radius. Manually, this works very well. The downside is that I'm working with thousands of miles of road, so we want to build an automated tool to be able to quickly analyze our entire existing network, and also be able to check roads digitized in future. I have adapted this script to find transition points where the curve begins/ends -- comparing line segments to see if the road direction is changing (e.g. is it still turning right, or has it started curving left instead). I therefore have chords for all the curves, but now I'm struggling with how to make middle ordinates. It needs to:   1. start at the midpoint of the chord   2. be perpendicular to the chord (I'm not positive how to do this yet, but found a script that will probably work and/or be adaptable)   3. be long enough to reach the road, but no longer (???) So, long question short: how can I find the (approximate) distance between my chord midpoint and the road? (Another StackExchange question describes working with the Near_analysis function which looks like it could work, but I currently have all of those chord midpoints just stored as a python list, not a feature class. I'd like to avoid generating a bunch of new features if possible.)


title: Merge point features within a given radius

text: I have a file in QGIS with 5000 points. Some of the points are really close to each other (5m). How can I cluster the point features to one point feature if they are within 5m? The new points should be centroids from the given points.


title: remove whitespace Python code results in odd character - Field Calculator

text: I need to remove whitespace from the end of a string in a field in ArcMap. I tried first the VBA function RTrim, but whitespace still remained. I thought maybe the whitespaces were tabs, so I tried 'Replace([field], vbTab, ""), but no luck. Then I tried the Python function rstrip, as so: !FIELDNAME!.rstrip() ..but this resulted in the whitespace being replaced by this character: Â There are several of these characters in each of the fields cells trailing the string. I have looked online for some meaning of this symbol, but I haven't found any. Doe anyone know what this means? I have still not successfully removed the trailing whitespace. here is a screenshot of the attribute column after I run the Python code given above: ![enter image description here](http://i.stack.imgur.com/kPe3V.png)


title: Compiling QGIS 2.0.1 on Centos (SIP problems)

text: Trying to compile QGIS 2.0.1 on CentOS 6.4 (64 bit):               /usr/share/sip/PyQt4/QtGui/qsyntaxhighlighter.sip: In function ‘PyObject* meth_QgsExpressionHighlighter_setCurrentBlockUserData(PyObject*, PyObject*, PyObject*)’:     /usr/share/sip/PyQt4/QtGui/qsyntaxhighlighter.sip:78: error: ‘qtgui_wrap_ancestors’ was not declared in this scope     make[2]: *** [python/CMakeFiles/python_module_qgis_gui.dir/gui/sipguipart3.cpp.o] Error 1     make[1]: *** [python/CMakeFiles/python_module_qgis_gui.dir/all] Error 2     make: *** [all] Error 2      This is with SIP 4.15.2, Python 2.6.6, and PyQt-x11-gpl-4.10.3 Any help please?


title: Can I create a table showing all the Srid values for all the tables in my DB?

text: Hi I can identify the current SRID of one table (at a time) in my database, but I was wondering if there was a way to bring back a single table that would tell me all the datasets that do not presently have a SRID against the geometry? This is in SQL Server 2008. Any help would be greatly appreciated.


title: Create new geodatabase in modelbuilder and place queried data in new database

text: ![Query existing database and place queried data into newly created geodatabase](http://i.stack.imgur.com/eMPs8.jpg) I have built a model that will take an existing database and query all feature classes in that dataset using a SQL Expression and I then want all the queried data from each feature class to be exported by feature class to a newly created database. I know I can manually create a new geodatabase and then set the output to that the geodatabase in the "Select" tool but I want the model to create the database automatically. The model here shows how I have tried to do it but by doing it the way illustrated here it creates a new GDB for each iteration, exports the feature to the new GDB and then moves on the next iteration, overwriting the previous feature class. How can I create the geodatabase first (once and only once), then run the iteration and place each each iterated feature class into the newly created database? Any advice?


title: Python Field Calc for Document Hyperlink

text: I am trying to create a script that will hyperlink to a .tiff or .pdf. These files are stored in one folder on a shared server. I want the script to look in a field: 'PLAT_IMAGE', which has the plat number (ex: 12048-001-01) and find it's match in the folder on the shared server and then write for example: \\\server\engineering\Asbuilts\Plats\12048-001-01.tif Making the document hyperlink active in ArcMap when using Identify My problem is that the calculation is going through and not giving me any errors but yet it's not calculating the field, so I don't know what I'm doing wrong!               import os, arcpy, fnmatch          def  x(PLAT_IMAGE):         path = "J:\\engineering\\Asbuilts\\Plats"         plat = PLAT_IMAGE     for file in os.listdir(path):         if os.path.isdir(os.path.join(path, file)):             for file in os.listdir(os.path.join(path, file)):                 if fnmatch.fnmatch(file,plat + '.*'):                     return "\\" + "\\mercury\gisdata$\engineering\Asbuilts\Plats" + "\\" + file          x( !PLAT_IMAGE! )      Any help is much appreciated


title: Georeferencer plugin will not load in QGIS

text: I have downloaded and installed QGIS. I would like to georeference an image using the core georeference plugin. When loading it no toolbar appears and the 'update', 'reinstall' etc boxes in plugin menu are also greyed out. Any clues?


title: How to get neighborhood values count using a geoprocessing tool in ArcGIS?

text: I have a raster (black=0, red=1, and white=NoData). I am wondering if there is any tool in arcgis which can help me only count number of 0s and 1s separately in a 3*3 neighborhood. I need an attribute table with that information. Any hint of relevant tool would be appreciative.


title: How to force ArcGIS server to create KML?

text: Attempting to serve arcgis server (AGS) 10 sde data to Google Earth (GE).   I don't have any problem opening the service or the network link in GE on my desktop.   When I try to open in ipad device nothing happens.   Finally today I got a message "kmz not supported on iphone". My network links generated by AGS are all kmz.   **Is there a way to force AGS to use kml instead?**   ...or has anyone been successful with AGS service on an ipad? **EDIT:** I do have specific reasons for using GE. I hope to negate those with some time. However everything I have built to this point has been flex. (I was informed that we now have ipad devices).   So either rebuild all my widgets in java, or find a solution that allows kml on the ipad.   Update: I get an error if I try to open kmz,   if I unzip it and put the kml in my network link folder it just doesn't work.   if I try to publish a network link with the kml ags says no kmz found.   I guess I am on to rebuild all my widgets in java.


title: Alternative to Copy Features in ModelBuilder to display a result

text: I have created a model that takes in a raster, converts it to a point shapefile, then extracts a set of points based on location, then further selects a set of points based on a SQL expression. I tried to run it once but it was taking absolutely forever. I would have thought this step to be akin to 'Export Data' in the regular window. Is this the right method? Current model creation below. The goal of this procedure is to find a region of high density crops in a selected country. I will eventually need to determine the coordinate footprint of the area so I can export that information to fine tune a forecast model to focus on that particular area. ![Model](http://i.stack.imgur.com/dH0Ou.jpg)


title: What are some lightweight python modules to convert LAS data formats to text formats?

text: I searched on the net and found pyLAS, but it seems to require quite a few other modules and I couldn't get the install right. Any pointers will be much appreciated. Thanks!


In Opengeo 4.0.2 suite, enabling labels causes map to show pink screen. The map shows up properly if labels are switched off. I am using Opengeo 4.0.2 suite on Ubuntu 13.10 server. How do I resolve this problem?


title: It is possible put highcharts chart in openlayers popup?

text:
As title. I build a openlayers map on my website.   and there are many feature vector on the map. now I want to put a highcharts charts in popup.   to show the temperature in every city   something Like this   Is it possible ? How? Thank you! I used contentHTML but it is failed highcharts code               $('#chart').highcharts({         chart: {             type: 'area'         },         title: {             text: 'Area chart with negative values'         },         xAxis: {             categories: ['Apples', 'Oranges', 'Pears', 'Grapes', 'Bananas']         },         credits: {             enabled: false         },         series: [{             name: 'John',             data: [5, 3, 4, 7, 2]         }]     });      part of openlayers code               var content = "

" ;             var popup = new OpenLayers.Popup.FramedCloud("popup",                 OpenLayers.LonLat.fromString(feature.geometry.toShortString()),                 null,                        content,                 null,                 true             );


title: Error 'GDAL extension not installed' while installing cartaro

text: I am trying to install Cartaro on my Ubuntu 12.04 LTS. I have followed every step from this guide: http://cartaro.org/documentation/installation But on the 'Verify Requirements' page I get this warning which says ''GDAL extension not installed'. Also when I click on proceed with installation, on Setup database page I see mysql and SQLite as options type not PostgreSQL as shown here Here is what I did: After using               sudo unzip geoserver-2.3.5-war.zip geoserver.war -d /var/lib/tomcat7/webapps/      I opened /etc/default/tomcat7 and added the following code in the end               JAVA_OPTS="${JAVA_OPTS} -DENABLE_JSONP=true"      Then I edited geoserver/WEB-INF/web.xml and added the following code just after web-app starting tag.                 <session-config>         <session-timeout>33333</session-timeout>       </session-config>      Please tell me what did I do wrong and what is the correct way to install the GDAL extension. Thanks


title: Declination Diagram

text: 1. What point on a UTM map does the declination diagram refer to? Is it the centre? Is it the edges? Is there a standard?   2. The central meridian on a UTM map points to true north. Are grid north and true north arrows identical only at this central meridian? If they are identical then what point on the UTM map does the grid north arrow on the declination diagram refer to?   3. If all grid lines on a UTM map point to true north (causing increased distortion with distance from the central meridian) then why is there a convergence angle (grid north to true north) on the declination diagram?


title: UTM convertor from one zone to next zone

text: I want to convert the utm cordinate systenm WGS 84 "E,N " from 38 N ZONE to 39 N or to 37 N. pleas send me the formula and programs that I can do it. thank you


title: How to modify existing grass gis tools and run them in QGIS?

text: I would like to take an existing tool of the grass tool box and modify the code to the specific needs of my task. As I´m struggling to start some kick- off would be great. How to access code? What do I have to consider? How to access/run the code in QGIS?


title: Can I upload an LPA file (Long Profile) to QGIS

text: I have a .lpa file that I'd like to upload to QGIS - Can this be done and if so, how? Thanks!


title: How to export qgis as a kml file?

text: I have a .qgis contour plot.   I would like to convert this to a .kml file.   How specifically would I go about doing this?


title: gdal_proximity distance to meters

text: I followed the answer given in this topic for searching the distance to the nearest coastline given a lat/lon coordinate. The distance output obtained is in degree (I think) but I need the distance expressed in meters. By citing the answer given: > There will be a units problem if you create the distance raster in > geographic coordinates. In order to get kilimetres you will need to > rasterize the sea in a suitable projected coordinate system and then project > the distance raster to geograpic coordinates to do the lookup. This will > give you the fastest response time as you're not projecting the query > coordinate at the time of lookup. At the moment my python code is the following:               from osgeo import gdal     ds = gdal.Open("/home/user/distance.tif")          gt = ds.GetGeoTransform()          topLeftX = gt[0]     dimX = gt[1]     topLeftY = gt[3]     dimY = gt[5]          band = ds.GetRasterBand(1)     #Numpy array of the distances in the whole area     data = band.ReadAsArray()     #Indexing by longitude/latitude position            def getDistance(coord):         x = int((coord[0] - topLeftX)/dimX)         y = int((coord[1] - topLeftY)/dimY)         distance = data[y,x]         return distance      The coordinate system is WGS84. I think I have to rasterize the shape file in a Projected Coordinate System but then I don't know what to do. Any suggestions?


title: Surface distance on DEM between two coordinates

text: Given two coordinates on a digital elevation model (DEM), how do I compute the actual distance traveled between the two locations, assuming a straight line route? How will changing the resolution of my DEM affect the results?


title: openlayers image layer stretch with no blur effect

text: I use an Image Layer to draw some data at specific locations on my web application using gwt-openlayers. Assume the original image is 400x400 pixels. Whem I zoom in the map, the image is stretched. The problem is that the stretch effect is 'blurry'. How can I avoid this blurry effect and instruct openlayers to apply a normal 'raw' stretch on the image? Is there some options somewhere? When I say blurry, I mean this: ![enter image description here](http://i.stack.imgur.com/YVNNg.png) The image above demonstrate the concept. I have a 400x400 pixels image. The image is formed by 4 sub squares 200x200 pixels each. Top left square is blue, top right is red, bottom left is yellow and bottom right is red. When I draw this image as an image layer in openlayers and I zoom the map, the center of the image gets blurry and the four colors are blended. Is there a way to avoid this blur effect and have a 'raw' stretch? Thanks. Here's the code (in gwt-openlayers) I use to create the image layer:               org.gwtopenmaps.openlayers.client.layer.ImageOptions imgOptions = new org.gwtopenmaps.openlayers.client.layer.ImageOptions();     imgOptions.setIsBaseLayer(false);     imgOptions.setTransitionEffect(TransitionEffect.RESIZE);     imgOptions.setAlwaysInRange(true);     imgOptions.setLayerOpacity(0.7);     imgOptions.setProjection("EPSG:900913");          imgLayer = new org.gwtopenmaps.openlayers.client.layer.Image("my_name", null, bounds, new Size(1, 1), imgOptions);     map.addLayer(imgLayer);


title: How can I relate the osm2po osm_source_id to data in OpenStreetMap?

text: I created and imported OpenStreetMap (OSM) data into Postgres using osm2po: e.g. `java -Xmx8g -jar osm2po-4.6.9/osm2po-core-4.6.9-signed.jar config=osm2po.config.eu_xc prefix=eu tileSize=x,c europe-latest.osm.bz2` and `psql -d osm2po_eu_xc -q -f "eu_2po_4pgr.sql"` which creates the table _eu_2po_4pgr_ and has datasets with columns _"osm_source_id_ ", _"osm_target_id"_. I created and imported the exact same osm data (using the exact same bz2 file) using osm2pgsql (0.69) e.g. `time osm2pgsql -smuc europe-latest.osm.bz2 -d osm_eu` which creates table like _planet_osm_point_ which has a _"column osm_id"_. I can't find the _"osm_source_id"_ in say _planet_osm_point_. What I would like is to join the osm2po table with the original OSM Data, in order to have the complete information. Is that possible or am I missing something fundamentally?


title: How to Resize the Qgis Main Window from python

text: I would like to write a plugin that simply resizes the qgis main window to make it small. How can I access the qgis main window through the qgis api? Use Case: We are having problems with qgis taking over the screen, so that the minimization toolbar is inaccessible, and we can't use anything else on the desktop until we have exited qgis. This seems to be a result of incapatability between Qt and our very old windowing system. I know how to fix this problem using the QGIS.conf file or the registry, but we have many users of varying ability, and many are switching between windows and linux (tunneled through Solaris) servers, which can cause the problem to recur. Fixing the problem with repeated edits of the config file (or the windows registry) is really too much to ask of them. So, I want to make a button that they can just press, and fix it.


title: GeoWebCache: misalignment of tiles

text: I'm trying to test simple web map based on PostGIS/Geoserver/OpenLayers stack on Windows Server 2008. When I load my WMS layer by:               var gwctest = new OpenLayers.Layer.WMS(         "gwctest",         "http://localhost:1979/geoserver/wms",         {         'layers': 'Index2000:index_vect_s100',         'format': 'image/png',         'srs': 'EPSG:3857',          'units': 'm',         'transparent': true         },         {'opacity': 0.8, 'isBaseLayer': false, 'visibility': false}          );      Everything seems to be fine. ![alt text](http://i.stack.imgur.com/pdHAh.jpg) However, when I try to add WMS layer to my map using GeoWebCache, replacing               "http://localhost:1979/geoserver/wms"      with:               "http://localhost:1979/geoserver/gwc/service/wms"      My data suddenly move from correct position: ![alt text](http://i.stack.imgur.com/p6H7y.jpg) Any ideas on what may be the cause? EDIT: some additional info on GWC mailing list post.


title: Error in performing transect in ambur R?

text: I'm trying to filter a transect, window size is 5 and filter individual baselines is set to 0. I select the transect shapefile, then the inner baseline shapefile and it gives me this error.               Error in SpatialLinesDataFrame(shape.prep2, trandata) : row.names of data and      Lines IDs do not match.      I did everything correct up to this point, add fields to shapefiles, construct transects, stuck on filter transects. Anyone know how to fix this error or knows what it means? The syntax for SpatialLinesDataFrame is               function (sl, data, match.ID = TRUE)     {     if (is.character(match.ID)) {         row.names(data) = data[, match.ID[1]]         match.ID = TRUE     }     if (match.ID) {         Sl_IDs <- sapply(slot(sl, "lines"), function(x) slot(x,              "ID"))         data_IDs <- row.names(data)         mtch <- match(Sl_IDs, data_IDs)         if (any(is.na(mtch)))              stop("row.names of data and Lines IDs do not match")         if (length(unique(mtch)) != length(Sl_IDs))              stop("row.names of data and Lines IDs do not match")         data <- data[mtch, , drop = FALSE]     }     if (nrow(data) != length(sl@lines))          stop("length of data.frame does not match number of Lines elements")     new("SpatialLinesDataFrame", sl, data = data)     }     <environment: namespace:sp>      not too familiar with code sorry. The error is happening because one of those if statements. Looking at the documentation on the R website I see that shape.prep2 is the object in SpatialLines class and trandata is an object in the data.frame class. "object of class data.frame; the number of rows in data should equal the number of Lines elements in sl." It must be something wrong with one of the shapefiles.


title: How can to merge DEMs in ArcGIS for Desktop 10.2?

text: I'm trying to merge some dems downloaded from USGS. When I used arcgis 9.2 there was an expression to be used in raster calculator to do that. I used it again but unfortunately it doesn't work in arcgis 10.2. What can I do?? I do not have the toolbar Production Contouring


title: Generating new shp files according to attribute in QGIS

text: I know that it is possible to color a vector map by attribute. I separated a map into 4 colors based on an attribute and now I want to generate a new map for each color. Is there a way to do this in QGIS?


title: Looking for web and desktop tools with GML support

text: I am a graduate student and a novice in GIS field. I have few GML files that contain vector polygon data and I want to see them in a map format. The problem is that I do not have access to ArcGIS/Grass GIS software installed in my machine. I am wondering if there are online GIS website where I can upload my file and visualize it. Also is there any simple and small size free/open source software that I can install on my computer?


title: How to import excel point and line data into ArcGIS 10.1 as shapefile?

text: I have two separate excel documents that I'd like to import into ArcMap, but I don't exactly know how to create a shapefile. The two files are as follows: 1) Point data. This spreadsheet has a long list of points (road intersections) with unique IDs and XY coordinates. This one shouldn't be too hard to import. 2) Line data. This spreadsheet has "From" and "To" columns with the unique IDs from the point data spreadsheet. It also has the XY coordinates of the points. Is there any way to import the line data spreadsheet so that ArcMap draws a line between each of the points? Another thing to consider is that many of these roads are not simply straight lines, so I don't know how to account for that. I'm using ArcGIS 10.1. Thanks!


title: Agent Analyst - How to reduce memory use during many, many loops?

text: I am currently using the Agent Analyst adaptation of Repast within the ArcGIS framework. The model I've built is capable of doing what I'd like (ie. randomly generating biologically-sensible migration pathways in order to find NDVI- maximized routes). However, in the loop I'm doing where I make an ArrayList for each randomization ('Positions'), and then iterate this into another, sorted ArrayList ('Movements', I find that after about 25-35k loops (ie. TickCount) the Java program crashes. Upon inspection it looks like the RAM used by Java during this process rapidly grows to take over all available memory. I'm looking for any input into what might be going on that is overloading the RAM (e.g., retained ArrayLists in the cache, 'new' rasters loaded rather than reloaded during each Week iteration). This Noob is very grateful for any guidance. Here is the relevant part of the NqPy code created within the Vector Agent part of the model (where I suspect the problem is occurring):               def step():       #self.model.loadTracker()       tc = self.model.getTickCount()       self.PABUid = int(tc)       self.totalNDVI = 1000000.0       point = self.the_geom.coordinate       point.x = self.model.startX       point.y = self.model.startY       self.setPositions(ArrayList())       positions = self.getPositions()       positions.add("PABU"+self.PABUid)       self.model.iterationWeek = self.model.startWeek        wk = self.model.iterationWeek       eW = self.model.endWeek       while wk <= eW:         self.rotate()         self.distance()         #movement output from the 'distance' action is evaluated to see if the agent would         #move into an inhospitable enviroment (eg. ocean)         self.propMove()         #print "CheckingNDVI"         nX = self.model.newX         nY = self.model.newY         self.model.rasterName = "NDVI"+wk+".tif"         rasterName = self.model.rasterName         raster = self.model.getRaster(rasterName)         val = raster.getPixelValueAtMap(nX,nY,0)         #print "NDVI value: ",val         self.PixelValue = val         if (val > self.model.minNDVI):           self.model.isSuitable = true         else:           self.model.isSuitable = false         self.calcMove()              positions.add("Week: "+wk)         point = self.the_geom.coordinate         positions.add(Double(point.x))         positions.add(Double(point.y))         positions.add(Double(self.PixelValue))         self.totalNDVI = self.totalNDVI+self.PixelValue         #self.model.trackPABU()         wk = wk+1         self.model.iterationWeek = wk         #self.model.writeAgents()         #self.model.updateDisplay()            positions.add(0, Double(self.totalNDVI))       #print positions            movements = self.getMovements()       iterator = movements.listIterator()       iterator.add(str(positions))            #print movements       Collections.sort(movements)       if tc>self.model.SimsToKeep:         movements.remove(0)                 if tc % 5000 == 0:         for strAgent in movements:           file = File(self.model.outputDir+tc+"output.txt")           fw = FileWriter(file.getAbsoluteFile(),true)           bw = BufferedWriter(fw)           bw.newLine()           bw.write(str(strAgent))           bw.close()         print tc


title: What is most suitable DEM to delineate watershed (from spot height, ASTER and SRTM)?

text: I have three DEMs of our country which are spot height DEM, SRTM DEM and ASTER DEM. SRTM has 90m resolution and ASTER has 30m resolution. My supervisor told me to delineate watershed using these three DEMs separately, and also using merged DEMs. What is the main differences between the three types of DEM that I mentioned above? How do I merge these three DEMs? What is the most suitable DEM to delineate watershed?


title: Alternative to ST_ClosestPoint, one which returns a vertex of the geometry

text: ST_ClosestPoint returns the closest point on a geometry, as in the example at http://postgis.net/docs/ST_ClosestPoint.html What I need is to get the closest point which is _already_ a part of the geometry. Example:               geometry = [1,2,3,4]     Point = any position      `ST_ClosestPoint(geometry,point)` returns something like `[2.5]` which is _between_ existing vertexes. How can I get a result without interpolation (so it should be 1,2,3 or 4 and nothing else possible)? Thanks!


title: Lost relation when importing from OSM to PostgreSQL with ogr2ogr

text: I am currently trying to import all administrative boundaries (countries, states, cities, etc) from OSM. I downloaded a small file from Geofabrik (Corsica), in order to run my tests. I understand that I could run the same operation worldwide using planet- latest.osm. First, I extracted all boundaries using the following command like               osmosis --read-xml file="corse-latest.osm" --tf accept-ways boundary=administrative --used-node --write-xml out1.osm      Then, I import everything into PostgreSQL using               ogr2ogr -f PostgreSQL PG:"[connection details]" out1.osm -lco COLUMN_TYPES=other_tags=hstore      Which fills up the multipolygon table with the list of cities, and that is exactly what I wanted. However, if I look for "Ajaccio", there is a lot of useful tags associated with the node 51341115, that do not show up in the multipolygon. Instead, they are associated to a point linked to that node ID. In the XML file, I have been able to link that point to the relation 73283, inside of which are standing all the ways of the city, and the tags that I actually find in the multipolygon table. The question is then the following: how can I link the multipolygon to the point in order to extract the tags, using solely the database? I'm not specifically tied to those tools, I'll use anything that works and that spares me writing code. Many thanks :)


title: How to display images as cross sections in 3d

text: I want to display a base map and several vertical cross sections in 3D. The base map image is easy to rectify as a geotiff. How can I also display the cross sections? The cross sections are at present just rectangular images. I can work out the depth/pixel conversion and the surface xy locations for each column of pixles. I have no DTM so the base map is flat and the cross sections simply reside directly underneath. My preferred options are (1) anything written in Python (2) open source GIS software. eg QGIS, uDig etc. (3) perhaps point cloud software like meshlab? Regards


title: Workflow in ArcGIS 10

text: Our agency has collabrated with six townships and one city for preparing an inventory of stormwater network. We are using ArcGIS 10 (2 ArcView licenses) and ArcPad 7.1 (stuck with it since our Mobilemapper CX is running on Windows CE). We are using ESRI File based Geodatabase (chose it because it allows multiple readers and one editor per feature datasets and it has one feature dataset for every political subdivision). We have finished 25% of the project. My primary concerns are: a) We cannot edit same feature dataset at a time. Mostly, my intern and I collect and edit data for same area/region at a time. b) Collecting data in the field is not efficient enough since we cannot see each others edits. (using ArcPad Data Manager) c) Sharing data with client is not a smooth task (converting the data into KML format so that it can be viewed on ArcGIS Explorer). I understand that the solution to my problem is ArcSDE and ArcGIS server but we don't have funding. We want to finish the project ASAP. Please recommend a workflow and open source software where we can work and share the data with clients more efficiently. Thank you.


title: How to use the output file of a script as an input file for the rest of the script in arcpy?

text: If my question does not make sense or needs clarification please let me know. I am struggling to solve it so I was struggling to phrase it into a question. I have this code here:               import arcpy,os     from arcpy import env          #Name .gdb and set location. This will create the .gdb     out_folder_path = ("C:/Users..../Documents/test")     out_name= "Test.gdb"     arcpy.CreateFileGDB_management(out_folder_path, out_name)          # This will load the .xml into the .gdb created above     target_gdb = "C:/Users/..../Documents/test/Test.gdb"     in_file = "C:/Users..../MASTER.xml"     import_type = "DATA"     config_keyword = "DEFAULTS"          # Execute ImportXMLWorkspaceDocument     arcpy.ImportXMLWorkspaceDocument_management(target_gdb, in_file, import_type)      I am struggling to convert this into a tool. The problem I am encountering is I want the user to just set the .env settings, then input the database name and then location. Then the tool will take over and create the file gdb and then load the .xml workspace into it. The .xml workspace has been hardcoded to a certain location. I am just struggling with the target_gdb for the .xml import. It does work as is with hardcoded locations and for one time use. But for those who do are not python saavy having to change code everytime is not feasible. Again I want this to become a tool. I did find this formum: Automating script to set geodatabase? but this just confused me even more. Updated code:02/4/14               import arcpy,os     from arcpy import env      # Name .gdb and set location. This will create the .gdb               out_folder_path = arcpy.GetParameterAsText(0)     out_name= arcpy.GetParameterAsText(1)     arcpy.CreateFileGDB_management(out_folder_path,out_name)      # This will load the .xml into the .gdb created above               target_gdb = out_folder_path + os.path.basename(out_name)     in_file = "C:/Users/.../MASTER.xml"     import_type = "DATA"     config_keyword = "DEFAULTS"      # Execute ImportXMLWorkspaceDocument               arcpy.ImportXMLWorkspaceDocument_management(target_gdb, in_file, import_type)


title: Is there a Deployment build of GeoExt 1.1?

text: The latest version of GeoExt that is available for download at http://geoext.org/downloads.html (version 1.1) appears to be a development version. They are individual files with whitespace and comments. Usualy there is a minified, compress version of javascript libraries available, to be used in production environment. Nowhere could I find such a version for GeoExt 1.1. Is there a deployment build available?


title: Removing duplicate rows from a shapefile

text: I would like to write some c# code to delete duplicate rows from a shapefile stored in a folder. Can anyone provide me with an example or pointers on the best way to attempt this, please?


title: Identifying Features With a Relationship Class

text: I have a feature class which represents points where vegetation clearing has occurred. I also have a table which has a list of vegetation species which were cleared at each location and the amount of vegetation which was in the sampled area. I have a relationship class in the geodatabase which links the two based on a `Clearing_ID` field. When the user clicks on a point with the Identify tool, they can drill down to see the list of plant species present. However, the list appears as numbers which is the `OBJECTID`. Ideally, I'd like the user to see the list of species present as the default field instead of the object ID. ![enter image description here](http://i.stack.imgur.com/v9YZt.jpg) Is there a way to do this?


title: Google Maps Markers inconsistent as they're drawn?

text: I've playing around with the google maps api and am puzzled at the following behavior. If I use `mPoint` as the LatLng for my marker, the marker is rendered on a different point on the map as opposed to putting the same value directly into the properties of the marker. Code chunk is as follows:                       var mPoint = [new google.maps.LatLng(38.991300,-76.936165)];     //        var GPoint = new google.maps.Marker({map:map,position:AVPoint});       var marker = new google.maps.Marker({         map: map,         position: new google.maps.LatLng(38.991300,-76.936165),         //position: new google.maps.LatLng(mPoint),         draggable: true       });      Why would this happen?


title: How can I develop a vehicle tracking application on a tight budget?

text: I need to develop and GIS application in .net with following Requirement, Can you please suggest the options how can i start with. I am new to GIS World, but I have a good knowledge. I neither developed any GIS application before! 1) We need to track vehicles on the map, the data recieved from Multiple Controllers connected Over Ethernet. 2) Display of Movements of Multiple Vehicles to be Provided for simultaneous operation. 3) Displaying Latitude, Longtitude, Altiitude, Roll, yaw, Pitch, Velocity with Every Vechicle. 4) Unique Icons for Multiple vehichles (icon should be configurable) 5) Provision for the modes to facilitate selection of vehicles for simultaneous display on Map. 6) Automatically move the Vehicle icon representing the respective vehicle, to the correct geographic location on map when vehicle change the position. 7) Automatic PAN movement of the area concerned on the digital map during Online Vehicle Tracking. 8) Calculate Actual disptance traveled by vehicle along with Aerial Distance, 9) Calculate relative distance between two vehicles. 10) Digital Map retrieval and display as per the area of operation. 11) Provision for additional icons creation on the map. (These icons should be displayed for futher use) 12) Provision of Information Layer (Peliminary/Detailed information about Terrain etc), 3D ANALYSIS 13) Geo-Referencing Facility. It should be able to digitize it by georeferencing. 14) Import digital masp irrespective of their original/native format. 15) OS Platform independency.   1. It should support Military grid coodrinate Data. 17) Provision for Defininf Waypoints with positioning Information 18) Multiple Access to various Digital map,image map like geotiff,.jpeg,.dted image 19) Integration of scanned maps, satellite images . 20) Zoom in Zoom out, Zoom in Max Scal1 1:1 and zoom out 1:50000 21) User Specified scale and Dynamic display of scale in the Mapview 22) Exchange of GIS Data for generation of Reports. 23) Import/Export Maps using .dgn (Microstation DGN) format.


title: How do I get the path of a layer file (.lyr)?

text: In C#, how do I get the path of a layer file? I can get the datasource and name of the layer but I need to get the path to the .lyr file. I want to record the path so my extension can copy and load the .lyr file later. This gives the path to the datasource (.mdb geodatabase in my case) that the .lyr file references:               var dataset = (IDataset) map.Layer[0];     var name = dataset.Workspace.PathName + "\\" + dataset.Name      Using that info plus the FeatureClass name I can open the layer and get the layer's polygons on the map but it is not the same as opening the .lyr file.


title: Why does field calculator expression sometimes return a binary value (0 or 1)?

text: I am running simple field calculator expressions and experiencing some of the silliest results I have ever seen. I am working on census data (living as a shape file) and trying to perform a simple [fieldName]/[fieldName] calculation on a newly created field. ArcGIS happens to enjoy playing around with me... and returns values 0 and 1 in no particular order/reason that I can really justify. Both input fields from the field calculator expression are of double data types. I have created a new field of double data type as well. However, the problem persists with other data types, including long/short integer. Is this a bug? Has anyone else experienced this? Shortest route to create + update a new field other than moving between SPSS/Excel and ArcGIS all day long?


title: How can I use different config.yaml files for application using mapfish print in geoserver?

text: System: Ubuntu 10.04 server / Tomcat6 / geoserver2.1.2 with integrated GeoWebCache I have various webmapping applications running which utilize the mapfish print module.Some are in tomcat6/webapps some are in my apache2 DocRoot. I need to have various mapfish print layouts to show different logos, legends, scalebars etc depending on which webmapping application is using it....but I only have one config.yaml file to style the pdf layout. Qustion: How can Have I various versions of config.yaml and link them to the different webapps? I´d be grateful for any replies, yours, Rob


title: How to create a maplayer from mapfile?

text: I am new in Mapserver, now I want to create maplayer , so its required shapefile for this, and if yes how can i used it. Map file is given below               eg2.map     MAP      NAME ENG1      STATUS ON      SIZE 800 500      IMAGETYPE GIF      IMAGECOLOR 240 240 240      SHAPEPATH "/www/geog5870/geog5870test/data"      EXTENT 87000 8000 660000 660000      UNITS METERS      WEB      IMAGEPATH "/www/tmp/"      IMAGEURL "/tmp/"      END      LAYER      NAME "EnglishGOR"      TYPE POLYGON      STATUS DEFAULT      DATA "England_gor_2001_area"      CLASS      NAME "English Government Office Regions"      OUTLINECOLOR 60 60 60      COLOR 255 255 0      SYMBOL 0      END      END      END


title: Transform an Irish grid reference (29903) into UK grid reference (27700) on insert

text: I have an insert script to input several hundred points, however some of the points use the Irish Grid Reference system. I have separated these out however was unsure how to convert or transform? these into 27700 on insert. I have the line below that grabs the UK reference, and the POST_GIS column is set to 27700. `st_geomfromtext('POINT(" . addslashes($data[6]) . ")',27700)`


title: Customise the splash image of QGIS

text: Where can I find the splash image which is shown while starting QGIS, and how to customise it?


title: What projection appeared in the NYTimes this morning?

text: ![enter image description here](http://i.stack.imgur.com/GmSG6.png) Perhaps not as GIS related as some, but definitely about maps. Can anyone identify the map projection that appeared in the New York Times this morning? It's not one I've seen before, and a quick scan through _Map Projections: A Working Manual_ did not yield the result I was looking for.


title: Programmatically start and stop map services with Python

text: Is it possible to programmatically stop and start ArcGIS Server map services using Python?


title: Reusing Layer Selection on Toolbar

text: I'm working on a custom toolbar for an extension which will have operations that work upon a specific, user-selectable layer. Is it possible to directly use the layer selector used by ESRI's toolbars in your own toolbar? For example, the 3D Effects toolbar contains a "Layer:.." combo box for selecting the layer on which to operate. It includes a full tree of the scene, filtered to certain types of layers. Is this control a reusable object?


title: How to generate a grid of horizontal equal-area rectangles in a Mercator projection?

text: I'm new to GIS so pardon my ignorance and bear with me. I'm generating a grid relative to a LonLat point. To make the cells all equal area, I'm creating it in 2163 projection. The problem is that when I reproject to 900913 and stick it on a Google Maps baselayer in OpenLayers, the whole grid is rotated clockwise. I understand this is a purely cosmetic issue, but it's a project requirement. Is it possible to change the map's projection to display it vertically? If not, is there another map source I can use that uses 2163 or can be projected to it? If all of this is out of a question, is there another way to generate equal area grids for Mercator 900913 which will line up vertically?


title: Does PostGIS not work with psql?

text: I have a sql script that creates several tables, a trigger, and a trigger function. All the commands work in pgAdmin when I execute them in the "Query" section. My goal is to put all these commands in a script so that I can set up the database on any machine relatively easily and quickly. When I attempt to run the script through psql, however, I get the following errors on all of the lines where I attempt to use `Geography` types:               ERROR:  type "geography" does not exist      It doesn't make sense because I can copy and paste the exact same queries into the query editor of pgAdmin and it works just fine. Does psql not support PostGIS? Or maybe just not the most recent version of PostGIS, 1.5.3, which includes the new `Geography` types? Is there something I can do to fix this? An example of one of my queries:               CREATE TABLE source_imagery (         id SERIAL PRIMARY KEY,         image_type VARCHAR(1000),         image_path VARCHAR(1000),         boundary GEOGRAPHY(POLYGON, 4326),         image_time TIMESTAMP,         catalog_time TIMESTAMP DEFAULT CURRENT_TIMESTAMP     );      Error message:               ERROR:  type "geography" does not exist     LINE 5:     boundary GEOGRAPHY(POLYGON, 4326),                          ^


title: ArcPad 10.0.2 - Removing shortcuts from ArcPad today

text: I have created a shortcut that goes to ArcPad today menu, but now I want to remove it. How does one do this?


title: Intersection problem: too many elements?

text: I have this code where an IGeometry shape (4-5 circles not filled) should be intersected with another layer's features                               // Convert from IGeometry to ADF Geometry Polygon                     ESRI.ArcGIS.Geometry.IPointCollection com_pointcollection =                             (ESRI.ArcGIS.Geometry.IPointCollection)multipointPoints2;                     ESRI.ArcGIS.ADF.Web.Geometry.Point[] new_adf_points =                                             ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.Local.Converter.FromIPointCollection                         (com_pointcollection);                     ESRI.ArcGIS.ADF.Web.Geometry.PointCollection new_adf_pointcollection = new                         ESRI.ArcGIS.ADF.Web.Geometry.PointCollection();                     for (int i = 0; i < new_adf_points.Length - 1; i++)                     {                         new_adf_pointcollection.Add(new_adf_points[i]);                     }                          ESRI.ArcGIS.ADF.Web.Geometry.Ring new_adf_path = new                         ESRI.ArcGIS.ADF.Web.Geometry.Ring();                     new_adf_path.Points = new_adf_pointcollection;                     ESRI.ArcGIS.ADF.Web.Geometry.RingCollection new_adf_pathcollection = new                         ESRI.ArcGIS.ADF.Web.Geometry.RingCollection();                     new_adf_pathcollection.Add(new_adf_path);                     ESRI.ArcGIS.ADF.Web.Geometry.Polygon new_adf_polygon = new                         ESRI.ArcGIS.ADF.Web.Geometry.Polygon();                     new_adf_polygon.Rings = new_adf_pathcollection;                                    // Now get the intersect query                          ESRI.ArcGIS.ADF.Web.SpatialFilter spatialFilter =                         new ESRI.ArcGIS.ADF.Web.SpatialFilter();                     spatialFilter.Geometry = new_adf_polygon;                     spatialFilter.ReturnADFGeometries = false;                     spatialFilter.MaxRecords = 10000;                          System.Data.DataSet queryResultsDataSet = new System.Data.DataSet();                     for (int index = 0; index < layerIDs.Length; index++)                     {                         System.Data.DataTable queryResultsDataTable =                       queryFunctionality.Query(queryFunctionalityq.Name, ConvertFromLayerNameToLayerId(m_Map, m_ResourceMapItem, m_firstLayerName), spatialFilter);                         queryResultsDataTable.TableName = gisResource.Name + "_" + layerNames[index];                         queryResultsDataSet.Tables.Add(queryResultsDataTable);                     }      I'm converting the IGeometry (which I've drawn on the screen and I know it's correct) to a web ADF Polygon and then using this polygon for an intersection query However some features which are clearly not intersecting the circles' borders are being included in the result To better explain the situation I did this screenshot: ![enter image description here](http://i.stack.imgur.com/tYr3u.jpg) Why is that? The red squared features are all different so it's not a problem of "same feature" Update: with less features it's the same, if they intersect the circle line they get colored, but the others get colored too.. simply adding them seems like telling the system "I'm intersected too"


title: How to select all files in a directory for merging in Processing?

text: I'm writing a generic script which involves writing shapefiles to a directory which are then merged together. After writing the files out to the _Output_ folder, I'm trying to get the **saga:mergeshapeslayers** algorithm to merge all the files in the _Output_ folder. I have used the Model Builder and although it is helpful to an extent, I find that it is used for specific purposes whereas I am attempting to make a script for generic purposes. Any advice would be more than welcome. (Code edited)               ##Test=name     ##Select_folder=folder     ##Result=output vector          import os     import glob          path_1 = Select_folder     path = path_1     os.chdir(path)          def function():              output = glob.glob(path_1 + './*.shp')         x = 0              while output[x]:             for fname in glob.glob("*.shp"):                 outputs_1 = processing.runandload("qgis:fieldcalculator", output[x], 'Number', 1, 10, 0, True, 1 , "C:\Users\Me\Desktop\Output\\"  + fname)                 multiple_0 = glob.glob("*.shp")                 x = x + 1                  if x + 1 > len(output):                 processing.runalg("saga:mergeshapeslayers", output[0], ";".join(multiple_0) , Result)                 break             else:                 continue          if path_1:         function()     else:        pass


title: popup.contentDiv OpenLayers Popup window

text: Could somebody help with popup.contentDiv()? I don't understand how can I define the HTML div class here. In general I would like to define some customs characteristics via HTML. The source http://dev.openlayers.org/docs/files/OpenLayers/Popup-js.html doesn't help much without an example. I found some link where it is used, but I didn't get how it works either. Thanks in advance.                                            popup=new myPopupClass("infoPopup",olMap.getLonLatFromPixel(event.xy),new OpenLayers.Size(10,10),popupInfoText, null, true);                                     popup.setBackgroundColor("#CCD8E7");                                     popup.setOpacity(0.9);                                     popup.setBorder("1 px solid #8DB2E3");                                     popup.contentDiv();                                     olMap.addPopup(popup);


title: Change points to polygon

text: I have a map which has one data point per zip code. I would like to change those points to zip code polygons. Any ideas on how to go about it?


title: How can I locate features along multiple routes in a consistent direction?

text: I have a large network of routes that I've created based on an NHD plus stream layer and a large number of points that are located throughout this network. I've got the routes segregated by stream as perfect as I can imagine getting them and the points layer is snapped to the routes. I've QAQC'd and the routes match the streams without duplicates..etc. Everything looks good. My goal is to determine how far upstream each point is on it's respective stream route. Ultimately this will lead me to stream distances from several fixed points in the watershed to every sample point. I've been using the Locate Features Along Routes Tool in ArcGIS 10 to do this. At first it seemed like it worked perfectly. But when I looked at the M values for each point I realized that for some routes, measurement started on the downstream end and progressed upstream, giving me what I needed, while on others, measurement started at the upstream end and worked downstream, effectively giving me the opposite of what I need. I think the definition of the "start" and "end" may be inconsistent among routes (I've checked on this via the Construct Points tool and am pretty sure this is the issue) Is there a way that I can ensure distance is being calculated in a certain direction on each route? Any input would be appreciated, Thanks for your time, Vin


title: Spike finder in PostGIS

text: I tried to modify spikeRemover from Schmidt & Krüger. My goal it's to obtain, for each polygon, the position of the spikes. I don't need to remove it. For some reason, that I can't figure out, my script fails and I don't obtain any position. Probably there is a problem at the union of different points but I don't know how to manage it in another way.                   CREATE OR REPLACE FUNCTION ST_SpikeFinder (geometry, angle double precision)          returns geometry as           $body$          DECLARE          ingeom alias for $1;          angle  alias for $2;          lineusp geometry;          newgeom geometry;          numpoints integer;          point_id integer;          result_pnt geometry;                   begin                 -- input geometry or rather set as default for the output              newgeom := ingeom;                  -- check polygon             if (select st_geometrytype(ingeom)) = 'ST_Polygon' then                 if (select st_numinteriorrings(ingeom)) = 0 then                      lineusp := st_boundary(ingeom) as line;                         numpoints := st_numpoints(lineusp);                 point_id := 0;                               -- the geometry passes pointwisely                 while (point_id <= numpoints) loop                             result_pnt := ST_Collect(result_pnt, st_pointn(lineusp, point_id));                         -- the check of the angle at the current point of a spike including the special case, that it is the first point.                         if (select abs(pi() - abs(st_azimuth(st_pointn(lineusp, case when point_id= 1 then st_numpoints(lineusp) - 1 else point_id - 1 end),                                      st_pointn(lineusp, point_id)) - st_azimuth(st_pointn(lineusp, point_id), st_pointn(lineusp, point_id + 1))))) <= angle then                  (probably the problem is here)      --->        result_pnt := ST_Union(result_pnt, st_pointn(lineusp, point_id));          end if;                 point_id = point_id + 1;                 end loop;             end if;             end if;                 return result_pnt;          end;          $body$            language 'plpgsql' volatile;      Thank you in advance. C.


title: Is there a Time Slider feature in MapInfo?

text: Is there an Time Slider Feature in MapInfo? The only conclusion I have found in my research is that you have to buy the MapInfo Profiler arm?


title: Sending xml post request to geoserver

text: This may be simple question but I can't get started. This tutorial from pen state university has xml request               <wfs:GetFeature service="WFS"  version="1.0.0"       outputFormat="GML2"  – we would like gml2 as the return format       xmlns:topp="http://www.openplans.org/topp"       xmlns:wfs="http://www.opengis.net/wfs"       xmlns:ogc="http://www.opengis.net/ogc"       xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"       xsi:schemaLocation="http://www.opengis.net/wfs       http://schemas.opengis.net/wfs/1.0.0/WFS-basic.xsd">       <wfs:Query typeName="topp:bc_roads">  we want to query                      bc_roads features         <ogc:Filter>           <ogc:FeatureId fid="bc_roads.1"/> only return features                                              with an ID of bc_roads.1         </ogc:Filter>         </wfs:Query>     </wfs:GetFeature>      My question here is how do I send this request to geoserver using geoext or javascript?


title: ArcGIS 10 - disappearing rasters?

text: From time to time it happens to me that rasters disappear at certain zoom level in ArcMap, or only some parts of the raster (strips) show up. I'm working on WinXP, ArcGIS / ArcView 10, Sp2, rasters are jpg files, both with and without pyramids. Anyone facing a similar behaviour? Any tips what is the reason and solution for this? It's really annoying.


title: Intermittent 'Could Not Find Layer' and 'Bad Request' Errors on GeoServer

text: Morning StackExchangers, We are having some intermittent faults with GeoServer at the moment. We have a mapping application - http://explorer.geowessex.com \- that makes WMS requests to our Geoserver to retrieve tiles. The layers are properly loaded and work via GeoServers built in 'Layer Preview' function, but sometimes in the application the tiles don't load. In the Firebug Net Panel you will see all the WMS Get Requests being fired when a layer is loaded and they say they return with no errors. In reality if you open one of these response in a new tab you get the following XML file.               <?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://wms.geowessex.com:80/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" ><ServiceException code="LayerNotDefined">       Could not find layer 'workspacename':'layername'      Which is stupid because it definitely exists, I can see it load perfectly well in GeoServer layer preview. What's even more annoying is that sometimes it doesn't even return this, it just gives the error '400:Problem communicating with GeoServer'. On top of that I cannot reproduce it regularly, I can have the map up, turn on a layer and it wont work. Turn on another layer and it will, go back and turn the other layer back on and it will start working again.   I'm really banging my head against a brick wall here. GeoServers own documentation is not that good for troubleshooting and I must've done enough Google searches to bring down their servers in the last two days and am still no closer as to why this happens.   Anyone experienced this themselves? Suggestions, comments anything is hugely welcome, and if you need any more info let me know. Cheers, Rob


title: How to set QGIS Python path so that I can use system packages on Windows?

text: Several Python packages are best installed with a Windows installer. But those installers install to the system-wide Python directory. How do I set the Python path in QGIS? I tried setting PYTHONPATH as a Windows environment variable. It can import netCDF4 (which is installed globally), but the import doesn't succeed because it cannot load a DLL. So perhaps also some DLL load path needs to be adjusted?


title: Incorporating known GPS/INS data to automatically produce geo-registered 3D models

text: I am trying to reconstruct a 3D building model from 2D video frames using an open-source Structure-From-Motion(SFM) software (Bundler). The final reconstructed 3D model will be in Cartesian coordinates but ultimately it needs to be geo-registered so that I can overlay it on a virtual globe or in 3D analyst. As far as I know, the typical SFM approach solves for _relative_ camera poses and the 3D points. However, since I already have GPS/INS metadata for each video frame, is it possible to use those as an initial estimate for the extrinsic camera parameters (camera pose) to input into bundler, and get my 3D point cloud directly in the same world coordinate system directly? I know a common approach (correct me if I am wrong) after the SFM process is to manually try to fit the model or find a rigid transformation to the world coordinate system, by using some marked ground control points on a reference imagery, but I just want to know whether this geo-registration can be done automatically instead by specifying known rough initial camera poses as input to Bundler (or any other similar sFm software) for subsequent refinement through bundle adjustment. Thanks!


title: Select Layer By Attribute using loop

text: I have 273 features in the feature class and i want to select them one by one and use the statement in the following code. At the line:               *arcpy.SelectLayerByAttribute_management (fc1,"NEW_SELECTION","OBJECTID=79")*,       The feature whose object id is 79 one of the features. When i use this code, I do it for only that feature. I want to apply the statement to all features one by one with a loop. I try to do it, but i am new for programming with python.               import arcpy     from arcpy import env     env.workspace = "D:\pomme de terre 2\pomme de terre2.gdb"     env.overwriteOutput = 1          fc1 = "Flowlin24K_sl_Uns"     fc2 = "New_Poyline"          arcpy.SelectLayerByAttribute_management (fc1,"NEW_SELECTION","OBJECTID=79")     arcpy.SelectLayerByLocation_management(fc1,"INTERSECT",fc1)     arcpy.SelectLayerByAttribute_management (fc1,"REMOVE_FROM_SELECTION","StreamLevel>6")     arcpy.CopyFeatures_management(fc1,"hidden")     arcpy.SelectLayerByAttribute_management (fc1,"CLEAR_SELECTION")          coordList = []          arcpy.FeatureVerticesToPoints_management("hidden", "vertices", "MID")     arcpy.AddXY_management("vertices")     rows = arcpy.SearchCursor("vertices")          for row in rows:         X = row.getValue("POINT_X")         Y = row.getValue("POINT_Y")         coordList.append([X, Y])          del row, rows          coordList.sort()          point = arcpy.Point()     array = arcpy.Array()          for feature in coordList:         point.X = feature[0]         point.Y = feature[1]         array.add(point)          polyline = arcpy.Polyline(array)          array.removeAll()     arcpy.Delete_management("vertices")     arcpy.Delete_management("hidden")     arcpy.CopyFeatures_management(polyline, fc2)


title: How do I restore the docking guides in ArcMap 10?

text: I lost my docking guides quite some time ago and I can't figure out how to get them back. Initially, I thought maybe it was an OS (Windows 7) configuration but my visual studio docking guides work correctly. I can still dock, but it's the old style of docking (no docking guides). ArcGIS 10 SP2 in Windows 7, x64 ![Viewing Multiple Tables](http://i.stack.imgur.com/PPYr4.png)


title: Is there a simple way to compute the UTM Zone from a lat/long point?

text: I am trying to convert a latlong point to UTM. To define UTM projection, I need to calculate zone for the point. I am trying to figure out the best way to do this. One way to do this would be use longitude values to find the proper zone. This would require lot of coding. I am curious if there is a better way of doing this. Thanks


title: Dynamic ArcServer Layer

text: I have the following setup: ArcGIS Server + ArcSDE on top of Oracle RAC. I created a view merging some data that is in Oracle and ArcSDE. In ArcMap it all works just fine. Whenever a new location is added to the locations table, I can see the joined features, etc. However I need this to happen in ArcServer. I mean, I have a MSD published and I can consume it via REST and OpenLayers. It seems to me that AGS is caching the generated tiles, although it's set to NOT cache anything. Any ideas why this is happening?


title: How to display IGeometry?

text: I've done an intersection with the ITopologicalOperator6 interface, and I now have an IGeometry with the result. How to draw this in a graphics layer on the map? I want to see what was intersected (possibly colored). According to this topic: CustomLayer Draw Com Exception it seems that I could create a new layer from the IMap interface.. is this possible? Is there any other more convenient way to do this? How would you draw an IGeometry item to your map?


title: Arc ASCII Grid output problem in QGIS 2.0

text: I used to convert tif files to ascii in Qgis 1.8 using gdal translate without having any problems with the output i.e. I have the normal heading and values. When I upgraded to 2.0 I started to get the output below (second image) which do not have the normal headings and values. The ascii grid output looks okay when loaded in Qgis but it seems broken when you check it on Textedit/Notepad/Terminal. I tried to convert the tif file using a newly installed Qgis 2.0 on windows, mac and ubuntu but they have the same results. Any ideas how to fix it? I used the data from worldclim http://www.worldclim.org/ as sample. Ascii output from Qgis 1.8 ![](http://i.stack.imgur.com/GGw8O.png) Ascii output from Qgis 2.0 ![](http://i.stack.imgur.com/DYgMJ.png)


title: Developing standalone gis applications using pyqgis

text: I am trying to develop a standalone application using pyqgis but I have no clue on where to start and where to write the code to make it happen so if anyone can provide me a tutorial to getting started, it'll be a huge benefit for me. Thanks.


title: Edit SQL Server table (not registered with geodatabase) with arcpy

text: I am trying to update some rows in a SQL Server database table (not registered with geodatabase) with arcpy.da.UpdateCursor. However, when running the `deleteRow()`, I get the `error return without exception set` error. Updating geodatabase tables works fine. I guess it might be a limitation of arcpy which can't edit a non-registered with geodatabase RDBMS table (in this case, SQL Server), but I want to make sure I don't miss anything. Can ArcSDESQLExecute class help me to update SQL Server table that are not registered in geodatabase (have found no info on that in the Help)? Any other workarounds? **EDIT:** I have ended up using ArcSDESQLExecute which works fine for executing SQL commands (including deleting rows) for tables that are not registered with geodatabase.               sde_conn = arcpy.ArcSDESQLExecute(r"C:\GIS\Temp\nongeo.sde")     sql_statement = """delete from added_table WHERE AddID = 2 """


title: Description CRS EPSG:6372

text: does someone know where could I find a description of CRS EPSG:6372. I searched in the EPSG page but seem not to find something that quick.. Any help would be really appreciated. Regards!


title: Complexity of A-star algorithm using pgRouting in a graph

text: I'm using the A-star algorithm in pgRouting. It uses the Manhattan distance as a heuristics ((|dx|+|dy|)/2). What would be the time complexity in a standard route network? Based on Why is the complexity of A* exponential in memory?, the complexity should be O(n) with n the maximum shortest path in the network. Is this correct?


title: Feature Cursor Count on mutiple layers

text: I would like to get a total count of multiple layer features (defined) within the ArcMap view extent. The code works fine with one layer but when I add a second or third layer with an "OR" operator I only get the count of 0 because one layer may be visible within the extent but the others are not. e.g. If pLayer.Name = "Layer1" Or pLayer.Name = "Layer2" Here is the code. Any suggestion would be greatly appreciated. thx!               Dim pMxDoc As IMxDocument     Set pMxDoc = ThisDocument          Dim pAv As IActiveView     Set pAv = pMxDoc.FocusMap     Dim pMap As IMap          Dim pEnv As IEnvelope     Set pEnv = pAv.ScreenDisplay.displayTransformation.VisibleBounds          Dim pSF As ISpatialFilter     Set pSF = New SpatialFilter     Set pSF.Geometry = pEnv     pSF.GeometryField = "Shape"     pSF.SpatialRel = esriSpatialRelIntersects          Dim pEnumLayer As IEnumLayer     Dim pLayer As ILayer      Set pMap = pMxDoc.ActiveView.FocusMap Set pEnumLayer = pMap.Layers               Dim pFLayer As IFeatureLayer     Dim pFClass As IFeatureClass     Dim pFCursor As IFeatureCursor     Dim pF As IFeature     Dim intCount As Integer      Set pLayer = pEnumLayer.Next               Do Until pLayer Is Nothing              If pLayer.Name = "Layer1" Or pLayer.Name = "Layer2" Then '?MULITI CUR                        Set pFLayer = pLayer          Set pFClass = pFLayer.FeatureClass          Set pFCursor = pFClass.Search(pSF, False)              End If         Set pLayer = pEnumLayer.Next     Loop             ' Loop through the cursor and count number of selected features.     Set pF = pFCursor.NextFeature     Do Until pF Is Nothing         intCount = intCount + 1         Set pF = pFCursor.NextFeature     Loop          If intCount = 0 Then     MsgBox "Count = 0"     Else     MsgBox ("Count = " & intCount)     End If          pAv.PartialRefresh esriViewGeoSelection, Nothing, Nothing


title: Selecting by bounding box with PostGIS

text: I have a large dataset in a PostGIS database. I am viewing it via a web interface that loads the data as GeoJSON. There is too much data to give the client all at once (10s of MBs worth...) so I would like to get an area of data within the bounding box formed by the edges of the window (the dotted line in the image below). Getting the coordinates for points f and g is easy. **Question 1:** Is this a good way of doing things? Should I be thinking of caching them as tiles or is this method likely to be efficient enough? **Question 2:** How do I retrieve just the data within this bounding box? **Question 3:** If a shape overlaps the edge of the bounding box (e.g. shape A below) is there a simple way to crop it as it is queried like in the second image? ![enter image description here](http://i.stack.imgur.com/y2o4Z.png)


title: Why is the OGC Symbology Encoding standard not being implemented and is SLD an alternative?

text: I searched exhaustively over the internet as well as at this How to export Symbology from ArcGIS to Microstation or QGIS? but didnt find a satisfactory answer. Hope to get the answers from the experts over here. Suppose I have styled a polygon shape file (eg: Border width = 0.26, border color = green, red colored simple fill, label displayed, etc ). When I save the style in QGIS, it is saved in xml format which is specific only to QGIS. Same is the case with other softwares which they store it in their own format. Which means, I cannot exchange styling information between different softwares. Though the OGC SE standard is in place, I was just confused why this is not being implemented by any software. Qgis 1.8 supports SLD but then what is the purpose of storing the style as qgis specific .qml file ? ..... Also, since SLD standard is used for web-services, is it okey to use this standard in this case(for exchanging data between different softwares) ?


title: usa electric transmission powerlines shape files sources

text: Does anyone have a good source for these? I have some old (1993) data from the US National Renewable Energy Laboratory (NREL) and Cubewerx' WMS (Cubewerx are a spatial data provider); neither of which are very accurate. Pipeline shape files other than the 10 year old NPMS data would also be helpful. Thanks


title: Force SRID in views in Postgis 2.0

text: I have a very tricky problem: I have a Postgis 2.0 database, used as datastore for Geoserver 2.1.3 I have a table with a geometry column (named shape), with SRID 4326 I have created a view like this:               CREATE VIEW pippo AS SELECT st_geometryn(shape,1) as geom FROM events      And then I published the view using geoserver. Everything works, except for GetFeatureInfo requests, that fails with the error 'Operation on mixed SRID geometries A simple               select * from geometry_columns      returns that the SRID of pippo.geom is 0 (but should be 4326) I tried also with this definition for the view               CREATE VIEW pippo AS SELECT st_setsrid(st_geometryn(shape,1), 4326) as geom FROM events      but I still have SRID=0 in the geometry_columns....any help? Here is how the layer is defined in Geoserver: ![enter image description here](http://i.stack.imgur.com/x7aYg.png) # Workaround I found a workaround to fix not the problem itself, but to let Geoserver to perform spatial search on GetFeatureInfo requests: instead of creating a view and publish it, I published an SQL-layer defined directly in geoserver. In the layer-definition form it allowed me to select the right gemoetry type and SRID and everything is working now!


title: Error using Grass v.net

text: I'm trying to connect my network to a points layer with a threshold of 4500 meters. I have about 11000 or so entries in my polyline and about 24 in my points layer. I have alreday used v.distance, hence the threshold. However I recieve an error (image attached) any time I try the following:               v.net input=network points=cdw1c output=network_v1 operation=connect threshold=4500      Here is the error: ![enter image description here](http://i.stack.imgur.com/JT0bb.png) Error content is: > Assertion failed: !p->taken[i], file split_q.e, line 69 Can anyone help me to understand what it means? * * * **Update:** After a long consulation with a generous user from GRASS community, I still havn't came to a solution. I've tried on two computers with WIN 7 as operation system, the following GRASS versions: 6.4.2 (through QGIS and without the plugin), 6.4.3 RCs, 7.0.0 For the first two I have received the same error, for version 7 I didn't recieve the error but the process have never finished. On the other hand, my helper managed to complete it in seconds using my files. Can the problem have something to do with windows 7, some background or system process, or so? Should I install any extention or plugin in order to run network analysis programs on win 7? Thanks a lot in advance. Dor


title: What are the FOSS equivalents to these ArcGIS products?

text: I'm a long-term user of ESRI software, and I now have a need to use free and open source software. I've been reading about Open Layers, PostgreSQL, PostGIS, GeoServer and MapServer but I can't find a good overview stating what each product does, why it's needed, and how they all fit together. For example, http://www.osgeo.org/ and http://freegis.org/ both list a bunch of products but don't give enough information for me to decide which ones are relevant. In ArcGIS, I would use the following:   * **file geodatabase** or **ArcSDE geodatabase** for data storage   * **ArcMap desktop** to edit the data and compile the map document   * **ArcGIS Server** to create web services   * **ArcGIS Server JavaScript API** to create the end-user maps from the web services What are the equivalent FOSS products?


title: How to deal with Roundabouts in PgRouting

text: Can anybody guide me to to deal with roundabouts in pgrouting?? I have some roundabouts in my road network and when i try to make route from one point to another, it gives me wrong direction where it finds the shortest path because length on road is my cost. In the picture below, blue cross is my starting point and other side is mt destination. ![enter image description here](http://i.stack.imgur.com/nbGkW.png) I am using simple pgr_dijkstra function and i am not considering any direction. boolean parameter for direction is set false. And there is no information about direction in my data. I just want to know what should I do if I want correct routing on my roundabouts. thanks in advance


title: How to serve a geotiff file using WMS?

text: I want to display a GeoTIFF file which is a air pollution map by use of WMS. This GeoTIFF file is an output of WPS. The following is my .map file and the URL request and the output. I use MapServer for windows.               MAP     NAME PM10     IMAGECOLOR 255 255 255      SIZE 600 800      IMAGETYPE PNG24 ## use AGG to for anti-aliassing      OUTPUTFORMAT       NAME 'AGG'       DRIVER AGG/PNG       MIMETYPE "image/png"       IMAGEMODE RGB       EXTENSION "png"      END # outputformat       PROJECTION      "init=epsg:3035" #latlon on etrs 1989 laea       END       EXTENT 3487500 2297500 4402500 3202500 # meters extents of region2       WEB       IMAGEPATH "c:/tmp/ms_tmp/"      IMAGEURL "/ms_tmp/"       METADATA      "ows_enable_request"   "*"       "map" "C:/ms4w/apps/airpollution/config.map"         "ows_schemas_location" "http://schemas.opengeospatial.net"      "ows_title" "Sample WMS"      "ows_enable_request"   "*"      "ows_onlineresource" "http://localhost:7070/cgi-bin/mapserv.exe?   map=C:/ms4w/apps/airpollution/config.map&"                     "ows_srs" "EPSG:3035 " #latlon            "wms_feature_info_mime_type" "text/plain"      "wms_feature_info_mime_type" "text/html"      " wms_server_version" "1.1.1"     "wms_formatlist" "image/png,image/gif,image/jpeg, image/geotiff"     "wms_format" "image/png"      END #metadata      END #web      LAYER     NAME "pm10"     DATA "pm10.tif"     TYPE RASTER     STATUS ON     METADATA      "ows_title" "pollution"       END #metadata       PROJECTION      "init=epsg:3035"      END #projection      END #layer pm10      END #map      I want the output format to be PNG. This is the URL request I am trying:               http://localhost:7070/cgi-bin/mapserv.exe?map=C:/ms4w/apps/airpollution/config.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=pm10&STYLE=&SRS=EPSG:3035&BBOX=97,5,106,21&WIDTH=600&HEIGHT=800&FORMAT=image/png      what I get as a response is an empty white image. What I want is an map like image which different values pixels displayed with different colors. I did the same thing by opening the same GeoTIFF file in arcmap. I know that I have to classify the raster file by use of SLD. but I get error there that "Could not open SLD sld/sld.xml and save it in temporary file " I want to be sure that the GeoTIFF file is accepted as an input data. The reason that I think it is not accepted is that when I change the name of GeoTIFF file to something else in .map file I still get the same output (white colored image).


title: Difference between Geostatistical Analyst and Spatial Analyst toolboxes for interpolation?

text: I did interpolation on point data using "Geostatistical Analyst" tool. It gives me surface layer. Then I exported this layer into raster image (by right click on this layer=>export data into raster) for convenience of raster calculation. But few minutes ago I found another approach on internet to create raster image from point data. That is to use "Spatial Analyst tool". So I'd like to ask there is any difference with output of further calculation between the two raster images created from 2 difference interpolation approach in ArcGIS for Desktop.


title: ArcObject 10.1 Java : Raster Analysis

text: I have been working with ArcObjects Java for a short time and have not been able to find a solution to the following problems: **How to get the Raster object from a IRaster reference?**               IRaster raster = rasterDataset.createDefaultRaster()          Raster rs = (Raster)raster      This gives classcastexception: _com.esri.arcgis.geodatabase.IRasterProxy cannot be cast to com.esri.arcgis.datasourcesraster.Raster_ I have been unable to find the documentation of IRasterProxy. Apparently it is deprecated!! **How to over write pixel values of an existing raster and store it permanently?** I am following this example but to no avail. I can see the pixel array is altered but when I open the raster in ArcMap I see no change. **How to get Raster Info?**               RasterDataset streamOrderRas = (RasterDataset)streamOrderRaster          Raster raster = new Raster()          raster.appendBands(streamOrderRas)          IRasterInfo rasterInfo = raster.getRasterInfo()      I get AutomationException: _No such interface supported, stacktrace shows it goes inside com.esri.arcgis.datasourcesraster.Raster.getRasterInfo()_ Please forgive if the questions appear stupid. If some one can show some working examples of the above problems, i'd be very grateful.


title: Can I create a 3D model in QGIS or is there some model for 3D city model?

text: Can I create a 3D model in QGIS or is there some model for 3D city models? I need these information for my project.


title: tools for sorting/filtering based on map state

text: Im about to embark on a project that list restaurants on a map. Id be grateful to find out what the best tools for the job are? The restaurants will be in a list beneath a map. The list will be refreshed based on the map zooming. The map will also have selectable/hoverable bounding boxes, so the user can select a location (london) and then drill down on the location to get neighbourhoods/boroughs in that area (westminster, clapham etc) Im looking at the following tools for the job - Gisgraphy (for data), Openstreemap (for maps and data) and leafletjs (for the boxes). Ive looked into geoserver but this seems to be something I cant host locally. Is there anything Im missing, or are there any better tools to achieve this? Also are there any examples I could look at where this has been done already? Cheers Ke


title: Esri Query Task service does not display

text:
I am trying to make my way through some ArcGIS server Javascript API tutorials. My current issue is that when I try to display an Esri.Map, all I see is the "powered by esri" logo in the bottom right. My init is below as is my mapDiv.               function init () {              map = new esri.Map("mapDiv",{           basemap: "hybrid",           center: [-80.94, 33.646],           zoom: 8         });              //Create your query         var queryTask = new esri.tasks.QueryTask("http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Petroleum/KSFields/MapServer/0");         //set the onComplete event handler, in this case, when the query is complete, call initLineID,         //production code would handle handle the error callback as well         dojo.connect(queryTask, "onComplete", initLineID);              //build and execute your query         var query = new esri.tasks.Query();         query.outFields = ["field_name"];         query.text = "all";         query.returnGeometry = true;         queryTask.execute(query);          }

      Any reason why this would not be working?


title: Why does Field Calculator say failure during processing, check geoprocessing results window for details?

text: I am using the Field Calculator and keep receiving the following error: 'There was a failure during processing, check the geoprocessing results window for details'. I am not very tech-savvy but I have tried a few simple fixes and still have had no luck. It's only a simple division equation so i'm confused. My equation is [Sum_NO_OF_] / [KM2AREA]


title: Calculate field based on most frequently occuring text/value in other fields in attribute table (ArcGIS)

text: I was wondering if anyone knows if it is possible to calculate a field based on the most frequently occurring value in other fields. What I have is 5 fields with crop information over 5 consecutive years. What I want is to create a new field with the most frequently occurring crop over that 5 year period. So for example if there were fields for Crop 2008, Crop 2009, Crop 2010, Crop 2011 and Crop 2012 with respective crop types for the first polygon as Winter Wheat, Spring Barley, Winter Barley, Winter Wheat and Grass the function would return Winter Wheat because it appears twice. I am aware there is a function called Counter in the Collections package that can do this for a list but I'm unsure how this would work within the arcpy field calculation tool? Likewise there also needs to be a consideration for values that tie i.e. appear the same amount of times across the 5 fields. I can't seem to find anything else from my searching about. I'm using ArcGIS 10 and hoping to calculate fields in the python window. Apologies for not providing any attempts at coding this so far, I really am stumped on this one!