title: POSTGIS workflow for cleaning up a noisy dataset

text: I have a conservation dataset that im trying to simplify for the purpose of using it as a GML layer in openlayers, ie reducing the total shape and node count to quite small. The data has issues: eg-   * about 15K orig polygons of which i only need the 500 largest post dissolve.   * overnoded, way too detailed   * lots of long thin slithers either inside or outside the main polygons which im not that interested in    * hundreds of invalid geometrys, rings with less than three points, self intersections etc Using this as an exercise whats a good way to clean and simplify a dataset like this? Ive been trying variations on this so far:   * st_simplifypreservetopology,10, to reduce time taken below   * hand fix invalids   * buffering +20m to swallow the internal slithers,    * st_union, trys ot make single feature, mem exhausts usually    * unbuffering 40m, to swallow external slithers.   * st_union,   * buffer +20 to return to right area.   * st_simplify,10   * st_dump   * bin any with area less than 1000ha Its slow and simplify always creates yet more invalid shapes. A subquestion is how to make postgis do the equivalent of qgis 'buffer and dissolve', ie merge shapes that adjoin or overlap, but leave other shapes as seperate features.


title: How can i upgrade geoserver on a geonetwork install?

text: I have a default geonetwork install with Jetty and geoserver but i would like to upgrade to the latest (unstable) build of geoserver instead of the stable version it comes with. Any ideas how i can do this? I don't mind reinstalling everything from scratch if i have to. I'm running Windows Server 2008r2 sp1.


title: Problem in Backing up a database in PostgreSQL

text: I need to create a backup from my database (in PostgreSQL 9.3). I have used the pg_dump tool, gives me an error. I have also used pgAdmin III for backing up which also gives me the following error: > pg_dump: [archiver (db)] query failed: ERROR: could not load library > "C:/Program Files/Boundless/OpenGeo/bin/../pgsql/lib/postgis-2.1.dll": %1 is > not a valid Win32 application. I tried checking the path but I couldn't find it. The nearest path that I reach is: > "C:/Program Files/Boundless/OpenGeo/bin" Any help would be highly appreciated.


title: KmlMouseEvent not working on Google Maps API

text:
I have this application where in I have loaded a kml file over the google map layer.I have these point layer markers from the KML. Now how do i obtain lat/long on the click of markers?. Please guide. **Please find the below code which i have used for kmMouseEvent which is not working:**               <%@ page language="java" contentType="text/html; charset=ISO-8859-1" pageEncoding="ISO-8859-1"%>     <!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">                    <meta name="viewport" content="initial-scale=1.0, user-scalable=no">     <meta charset="utf-8">     <title>Transit layer</title>          <link href="/maps/documentation/javascript/examples/default.css" rel="stylesheet"      type="text/css" />

                I want the kmlevent to work on the local kml and not on the kml which is on server.Please guide.


title: Query data without opening up database?

text: I have a bunch of polygon data that I want users to be able to sort/filter based on certain criteria. I have the data in a CartoDB table, and can pretty easily use some SQL calls to achieve what I want. The only problem is that my CartoDB table has to be set to public for this to work, meaning someone could dump the data from my table pretty easily. So I'm wondering if there's some way I can send the polygon data from CartoDB into some other service, which could filter the query/filter data without opening up my table? My apologies if this question isn't clear..


title: How to implement Routing in openlayers

text: I need to add routing functionality in my application how to add the routing in openlayers.


title: Problems with intersection in FME and ArcMap

text: I have one very interesting issue which I do not know how to solve/understand. In FME I used the `LineOnAreaOverlayer` transformer to get new intersected lines, these some are inside and outside polygons. When I visualize it with FMEViewer and visualize them with ArcMAP I get different results. In ArcMAP it is visible like line was not intersected according to the polygon, but at the same time in FMEViewer everything is fine. help please! ![This is situation within FMEViever](http://i.stack.imgur.com/9NhSs.jpg) ![And this is within ArcMAP](http://i.stack.imgur.com/lSU7j.jpg) Both representing same line and area. And I tried to do intersection within ArcMap and intersections occures at the same place as with FME: ![intersected with ArcMap](http://i.stack.imgur.com/7x7cg.jpg)


title: Is there a silent installation command without the desktop shortcuts?

text: I installed QGIS with Microsoft SCCM 2012 and I want to make a silent install without desktop shortcuts. I managed to do a silent install with / S but I found no s command to delete desktop shortcuts. I have Windows 7 64 client PCs Thank you


title: How can I publish MrSID files in geoserver?

text: I need to publish MrSID files using geoserver 2.1.4. I am unable to even get the plugin to add the store in geoserver. Does anyone have any suggestions on how to achieve this?


title: Given a set of coordinates, how can I calculate the minimum bounds?

text: I have a set of coordinates in lon/lat. There are always at least 3 ordered coordinates that form a poly (if they were projected flat). How can I calculate the minimum bounds for these coordinates as a set of valid longitude and latitude ranges? (by 'valid' I mean ranges that explicitly account for the antimeridian). Its kind of hard to explain what I'm looking for so I have a picture. ![enter image description here](http://i.stack.imgur.com/ih0iT.png) The case in the picture is trivial; you just find the absolute minimum and maximum for all the points. This doesn't work for all cases though. Is there a generic solution? Edit: To clarify what I mean by 'valid', say I had three longitude values in my data set: -76, -135 and 164. The values cross the antimeridian and I'd like the resulting ranges to be split: -76 to -180 AND 164 to 180. Some more clarification. The points form a polygon, so in certain cases, the required range could be from -180 to +180 (ie, the full 360 degrees): ![enter image description here](http://i.stack.imgur.com/01LE9.png) The image on the left shows the longitude of four coordinates that occur on one 'half' of the Earth. Imagine it as if you were looking down onto the north pole (black dot). The pink shows the minimum longitudinal range that encompasses the polygon (the polygon is shown between the four points in purple). The case on the left would have two longitudinal ranges: [-180 to -120] and [135 to 180] (just visually estimating it) The image on the right shows another case where the points go all the way around the Earth. This range would be [-180 to 180].


title: How to get elevation with LibGeoTiff?

text: I'm trying to get an elevation profile from a geotif file ( downloaded from http://srtm.csi.cgiar.org/). I downloaded the libgeotiff library but i can's see any function that can return an elevation for given pixel/coordinate or any useful documentation/example to how to do that. Can someone help me out how to get the elevation? i want to do it in c\c++


title: How to 'burn' a stream network into a DEM layer in ArcGIS 10?

text: I am working on a DEM that in an area where there is little relief. I would like to 'burn' the river network into the DEM so I can calculate flow accumulation and flow length accurately. I am using ArcGIS 10 but if someone knows the theory or has a good reference that would be appreciated. Thanks.


title: Parameterizing SQL Views Geoserver

text: I must do a view in geoserver. I seen some example here: http://docs.geoserver.org/stable/en/user/data/database/sqlview.html. The example show the operator like for parameterizing view. Suppose, now that in WFS request there is a value for a field (named filed1) example `&viewparams=value1:a`. If i use the operator like               Select filed1, filed2, field3            From  myTable      where  field1 like %value1%       the view return tuples where filed1 is a, aa, aaa, aaaa, aaaaa, etc, but i want tuples where filed1='a'. Can i use the operator =, if yes how write my query?               Select filed1, filed2, field3            From  myTable      where  field1 = .......


title: Iterate Workspaces in Modelbuilder not executing all tools

text: I have a folder which contains a number of file gdbs. These gdbs all have the same structure i.e. a feature dataset `FD_1` and several feature classes within the dataset (the two I am interested in are `FC1` and `FC2`). The feature classes all have the same names within the gdbs. I want to create tables from the feature classes, so my model does the following: `Iterate Workspaces` produces `Workspace` (the path of the current gdb) and `Name` (the name of the current gdb). The `Copy Rows` tool is used to produce a table for `FC1`. The input table is specified using inline variable substitution               %Workspace%\FD_1\FC1      As is the output table               %Output%\table1_%Name%.dbf      The problem is, the `Copy Rows` process is not executing. The model simply iterates over the gdbs and then ends. When I check the `Copy Rows` window, it says that the input table does not exist. I tried setting `Workspace` as the `Current Environment` for the tool, but it didn't change anything.


title: 3d arcgis raster .tiff image to CFD software

text: i am trying to export 3d raster image prepared in arcgis to cfd package such as Fluent. I was wondering if any one has experience to convert raster image to autocad 3d format, stl, vrml or any other format of its likes which could be imported to cfd package. I could import format likes IGES, .stl, .DWG, .DXF or even text files. any help provided will be appreciated.


title: osm2pgsql import planet.osm.bz2 slowness

text: I am trying to import whole planet data to my i7 3770k + postgres on SSD + 16gb RAM in order to render my custom tiles. But the process is painfully slow and multi core is not helping (though you have specified)               osm2pgsql -d planet -U osm --slim -C 8000 --number-process 4 planet-120704.osm.bz2      From `htop` what I can see is osm2pgsql take 100% of 1 core and postgres takes only 10% of cpu. and current prompt is:               Processing: Node(208320k 147.0k/s) Way(0k 0.00k/s) Relation(0 0.00/s)      Q: Is my computer too slow? Any method to speed up (without complexing too much hardware)?


title: How to create a TAB with custom CRS using QGIS?

text: I meet problem about reproject in Qgis. I created coordinate user defined. i used version Qgis 2.2. I defined coordinate: **Name:** VN 2000 / 107.75 3degees **Parameters:** +proj=tmerc +lat_0=0 +lon_0=107.75 +k=0.9999 +x_0=500000 +y_0=0 +ellps=WGS84 +towgs84=-191.90441429,-39.30318279,-111.45032835,-0.00928836,0.01975479,-0.00427372,0.252906278 +units=m +no_defs   1. when i save as (.shp format) new file to coordinate i create in here. It gets full parameters. ![enter image description here](http://i.stack.imgur.com/s5sUN.png)   2. But when i save Mapinfo (.tab) format. It doesn't get full parameters (it's missing towgs84): ![enter image description here](http://i.stack.imgur.com/OtSuI.png) And if i open file in mapinfo software and overlay in Google earth, coordinate is not correct because it doesn't get 7 parameters (towgs84). Can you help me? why it get full parameters? Thank you very much


title: Publishing shapefile on localhost

text: I need to add stores (shapefiles) to geoserver running on localhost:8080. In data_dir/data I have directories containing shapefiles. I add them as workspaces with Namespace URI set to `localhost:8080/geoserver/data/*name_of_the_dir*/`. When I try to add new store, the URL field in Connection Parameters doesn't open though and I can't choose which shapefile I want to add. I suppose the only error might be in `localhost:8080/geoserver/data/*name_of_the_dir*/`. I tried many different variations of this path, none of them worked for me. So the question is rather simple (hope the answer too): How do I add a store to a workspace when working in localhost? Thanks for any help. PS: The geoserver I'm working on is meant to run offline, thus I can't enter web URLs in Namespace URI.


title: How to get multiple SLD files of the same shapefile?

text: I need to get multiple SLD files for the same shapefile. Geometry is polygons, and they need to be colored by one attribute column. What is the best solution? I'm working with ArcMap2SLD or ArcMap to SLD Converter.


title: Radius based search with Leaflet maps

text: Using leaflet I have created a simple map, played around with markers (adding/removing) and have added some GeoJSON data polygons etc. I now need help with implementing a radius search. For example I want to restrict the markers placed on the map by radius/distance in metres/cm etc. I want to achieve this in two ways. example:   1. The user will enter a distance in metres and click submit. This should only bring the markers back within that distance/radius. My markers are populated using Lat/Long coords. Is this possible to achieve what i want or do i need additional data along with Lat/Long.    2. The user should be able to drag a marker on the map and click submit which should return all the markers within 300 metres of the marker.  Is the above possible with leaflet and can anyone point me to any plugins which can achieve the above.


title: Drawing order based on date field

text: I need a polygon layer to draw in order based on a date field... What I'm trying to do: I need to display an easement feature class with two different symbols: one for easements and one for vacations. In addition, I need this layer to draw in an order that represents the date in which they were recorded (stored in a date field). Using symbol levels is not an option because I am working with thousands of different dates (Also, ArcMap symbology options do not recognize date fields as dates...). The sort tool, which reorders the features in the database, is also not an option because this feature class is constantly being updated and it is not feasible for this tool to be ran frequently. Any ideas?


title: ArcGIS WebADF Distance tool

text: Is there a tool in the WebADF that will let me draw a straight line and calculate the distance between the two points.


title: Gdal2Tiles getting tif with dtype = uint8 but want int16

text: I have a gdal2tiles code for getting tiff file. In my code I am getting true tiff format but it makes my dype uint8. I want int16 for dtype. how can I overcome it? I know I have to add gdal.GDT_Int16 somewhere but I don't know where... the relevant parts of the my code:               default_tileformat = 'png'          tempdriver = gdal.GetDriverByName( 'MEM' )     tiledriver = gdal.GetDriverByName( default_tileformat )          ...          def writetile( filename, data, dxsize, dysize, bands):     """     Write raster 'data' (of the size 'dataxsize' x 'dataysize') read from     'dataset' into the tile 'filename' with size 'tilesize' pixels.     Later this should be replaced by new <TMS Tile Raster Driver> from GDAL.     """          # Create needed directories for output     dirs, file = os.path.split(filename)     if not os.path.isdir(dirs):         os.makedirs(dirs)          # GRR, PNG DRIVER DOESN'T HAVE CREATE() !!!     # so we have to create temporary file in memmory...          #TODO: Add transparency to files with one band only too (grayscale).     if bands == 3 and tileformat == 'png':         tmp = tempdriver.Create('', tilesize, tilesize, bands=4)         alpha = tmp.GetRasterBand(4)         alphaarray = (zeros((dysize, dxsize)) + 255).astype('l')         alpha.WriteArray( alphaarray, 0, tilesize-dysize )     else:         tmp = tempdriver.Create('', tilesize, tilesize, bands=bands)          # (write data from the bottom left corner)     tmp.WriteRaster( 0, tilesize-dysize, dxsize, dysize, data, band_list=range(1, bands+1))          # ... and then copy it into the final tile with given filename     tiledriver.CreateCopy(filename, tmp, strict=0, options=gdal_driver_options)          return 0      and another part:                # Load raster from read window.      data = dataset.ReadRaster(rx, ry, rxsize, rysize, dxsize, dysize)      # Write that raster to the tile.      writetile( filename, data, dxsize, dysize, bands)


title: Looking for a Tool to Calculate Distance Between Points (Between Layers)

text: I am a college student using ArcEditor 10.0 for a semester project measuring distance to hospitals from census tract centroids. This is my first time using ArcGIS, and I have no formal training. So this question may be very simple, but is it possible to calculate distance between two different types of points between two layers? In one layer I have a shapefile for census tracts and I've calculated the centroids. In the other layer I have geocoded hospital addresses. I would like to calculate distance from each census tract centroid to the nearest hospital. Is there a tool in Arc Toolbox that will do this/could you give a brief overview of the tool? Thank you!


title: What is the best way to manage permissions as root installing postgresql and postgis?

text: I've encountered an issue numerous times where I lock myself out of my postgres installation as a result of modifying my pg_hba.conf file and therefore screw up my postgres installation because I don't have root access. So I want to manage everything as root and create all the right users with permissions and then use my linux terminal only installation as a non-root user and still be able to manage everything necessary while limiting my user access to pg_hba.conf. Ergo, if I want to have my postgresql installation open-ended and accessible what settings do I need to make in pg_hba. And by that I mean I want to connect remotely using tcp/ip.


title: Error using SLD with MapServer based WMS (GeoBC)

text: I'm trying to customize rendering of map tiles from a government WMS server. When I use the GetStyles WMS query to retrieve the SLD for a layer, then feed that SLD directly back in a GetMap query, the server returns an error. I'm a newbie at GIS stuff so I'm probably missing something basic. Here are the gory details. Here's a request to display a tile from this server using two layers: http://openmaps.gov.bc.ca/mapserver/base3?service=wms&version=1.1.1&request=GetMap&styles=&srs=EPSG%3A4326&width=256&height=256&format=image%2Fpng&layers=BC_WATER_LINES_500M,BC_MAJOR_CITIES_POINTS_500M&BBOX=-122.343750,54.977614,-120.937500,55.776573 The capabilities for the WMS do include user styles: `<UserDefinedSymbolization SupportSLD="1" UserLayer="0" UserStyle="1" RemoteWFS="0"/>` I used this request to get the SLD XML document for just the cities layer: `http://openmaps.gov.bc.ca/mapserver/base3?service=wms&version=1.1.1&request=GetStyles&layers=BC_MAJOR_CITIES_POINTS_500M` It returns good looking XML:                     I made this file available, unmodified, on the web at `http://junkpile.net/misc/sld_unmodified.xml` If I then make a request using this SLD in the URL: http://openmaps.gov.bc.ca/mapserver/base3?service=wms&version=1.1.1&request=GetMap&styles=&srs=EPSG%3A4326&width=256&height=256&format=image%2Fpng&layers=BC_WATER_LINES_500M,BC_MAJOR_CITIES_POINTS_500M&BBOX=-122.343750,54.977614,-120.937500,55.776573&sld=http%3A%2F%2Fjunkpile.net%2Fmisc%2Fsld_unmodified.xml The server returns this error:               <?xml version='1.0' encoding="ISO-8859-1" standalone="no" ?>     <!DOCTYPE ServiceExceptionReport SYSTEM "http://schemas.opengis.net/wms/1.1.1/exception_1_1_1.dtd">     <ServiceExceptionReport version="1.1.1">     <ServiceException>     msDrawMap(): Image handling error. Failed to draw layer named 'BC_MAJOR_CITIES_POINTS_500M'.     msSDELayerNextShape(): SDE error. SE_stream_fetch(): Underlying DBMS error. (-51)     agg2RenderTruetypeSymbol(): TrueType Font error. AGG error loading font ()     agg2RenderTruetypeSymbol(): TrueType Font error. AGG error loading font ()     msDrawMap(): Image handling error. Unable to initialize image.     msPrepareImage(): GD library error. Unable to initialize image.     msDrawMap(): Image handling error. Unable to initialize image.     msPrepareImage(): GD library error. Unable to initialize image.     </ServiceException>     </ServiceExceptionReport>      What am I doing wrong? Is there some misconfiguration on the server side, or do I have a fundamental misunderstanding about how user SLD files are applied, or both? FWIW, I also tried including the SLD in SLD_BODY but ran into URL length limits... tried including the SLD in post data but couldn't find a DTD or XSD for the GetMap XML post format that MapServer expects. The lack of documentation/examples is a challenge. **EDIT:** I still don't understand why this isn't working, since the SLD came directly from the server in question. But I am able to successfully apply styles if I use ogc:PropertyIsLike instead of ogc:PropertyIsEqualTo. I guess this is a bug in MapServer?


title: Help Setting up GeoServer Proxy with OpenLayers

text: I have a GeoServer WMS in Apache-Tomcat at http://localhost:8080 and a mod_wsgi app in standard Apache at http://localhost:80/myapp. My OpenLayers map is being hosted statically by the Apache-Tomcat Geoserver app at http://localhost:8080/geoserver/www/OL.html. The JavaScript in this map has to call the mod_wsgi app for non geographic data. I've tried setting up a proxy in GeoServer the following ways               http://localhost:80/mypp          http://localhost          localhost      None of these work. After making a request some of my response object attributes look like this which means it's not really working. GeoServer logs show no errors               priority "NORMAL"     readyState 4     responseText ""     responseXM Lnull     status 0     statusText ""      Can't say I'm too much of an expert at Proxies. But it seems to me that I wouldn't need to set up a OpenLayers proxy (as mentioned here) if GeoServer can proxy for me. Ideas, hints, advice?


title: How do I get the x- and y-coordinate values of points in UTM using QGIS?

text: To get to know the size of the study area and how far objects are separated, I need to know x- and y-coordinate values. Although the $x and $y of Field Calculator returned latitudes and longitudes, x- and y-coordinate values were not obtained when it was in Universal Transverse Mercator (UTM). I think it should be possible to get x- and y-coordinate values since more complicated calculations, i.e. latitudes and longitudes, are possible.


title: How to return shapes from geo service?

text: I have a geo service that seems to finally be working. I have the service setup with a list of Polylines as the returned object. This is what I'm trying to implement but it doesn't work:                IPointCollection points = new MultipointClass();      points.AddPoint(station);      points.AddPoint(targets[index]);      messages.AddMessage(string.Format("Computing distance from origin to target {0}", index));      Polyline result = distanceComputer.CostPathAsPolyline(points,          hydrographyRasterDataset,           hydrographyRasterDataset) as Polyline;      IGeometry feature = result as IGeometry;      messages.AddMessage("Getting new row for output.");      IFeatureBuffer buffer = outputClass.CreateFeatureBuffer();      IFeature polyline = (IFeature) buffer;      messages.AddMessage("adding shape");      polyline.Shape = feature;      messages.AddMessage("adding object");      polyline.set_Value(polyline.Fields.FindField("SHAPE"), result);      messages.AddMessage("adding length");      polyline.set_Value(polyline.Fields.FindField("SHAPE_Length"), ((IPolyline)result).Length);      messages.AddMessage("Commiting to output");      outputCursor.InsertFeature(buffer);      When I try setting the value for "SHAPE" I get an exception and the message is:               This spatial reference object cannot be defined from the available information      So how exactly to I set the polyline in the return object (basically adding rows of Polylines to the OutputFeatureSet) so that it does not blow up?


title: Select Layer by location does not work

text: I am having trouble using the "Select Layer by Location" tool. I want to select everything from shapefile A (green) that overlays with Shapefile B (red). Here is the raw data: ![raw](http://s7.postimg.org/tm4vzn2jf/raw.png) And one example of the "Select Layer by location" (all three of them have the same adjustment, just different input): ![adjusments](http://s28.postimg.org/cc1wbpf25/select.png) Running this, it produces the following output. ![output](http://s28.postimg.org/vu4b9kewyu/raw.png) In a closeup one can see a few selections. But based on the input data (red - as you can see, it basically covers all of the green area) it should be exactly the same. Or am I totally mistaken? ![closeup](http://s30.postimg.org/3qudjj1gx/closeup.png) EDIT 7/14/14 4:58pm This is how it worked out. Instead of making one feature layer I made three feature layers and three different selection which need to be saved seperately. Another step was copying the features before merging them. Dont be confused, I added another step to the model (calculate field) but this has nothing to do with the trouble shooting of the original problem. ![solution](http://s14.postimg.org/6hltytbf9d/Untitled.png)


title: Whats better for selecting? WMS or WFS

text: I need to select from a layer(s) and I am wondering what way I should serve the layer. The requirements are that I need to be able to select multiple features at a time by drawing a box or polygon. They should become highlighted and I should display the data for all the features. I don't need to do any editing of them. I'm using geoserver and openlayers


title: Import Points with WKB format to QGIS

text: I have a .csv file with population data according to specific locations (Points). My problem is that all locations are given in WKB format like that: 0101000020E6100000B953EFEEEEB253C0F15BD4CCCC4C05C0 How can I visualize my data with QGIS? I don't know why there is no simple solution. Because to visualize WKB from a database is no problem. I hope you can help me.


title: Incremental network loader Error - Geometric netwotk QA

text: Trying to merge 2 geometric networks with Incremental network loader I got this error -![error picture](http://i.stack.imgur.com/uQIBQ.png) How can I find this Closed polylines at my geometric network ?


title: adding helicopter video files to Google Earth or QGIS

text: Is it possible to add a helicopter's video files to Google Earth or QGIS ? We have a flight path points derived from the xml file associated with the .mpg file which I envision being used to pilot the video. The .mpg files range in size 100mb to 2gb.


title: Re-projecting WGS84 to OSGB in SAGAGIS

text: I took some points with a GPS, and I'm trying to convert from GCS to OSGB 1936 British National Grid in SAGAGIS. I'm working with a shapefile. To reproject the vector layer to the OSGB projection I have used Proj. 4 (Command Line Arguments, Shapes). Once the vector layer is reprojected, I loaded a .tif image to check that the projected points lie in their correct places. I set its projection using Proj.4 --> Set Coordinate Reference System. However, when I try to overlay layer with the image, and the shapefile with the re-projected data, I can't overlay properly, and only appear the layer with the gps points but not the image. So my question is: How can I overlay those two properly in order to check that the points are in the adequate place? Sorry for the question, it could be a beginner one, but I appreciate any help!


title: How to obtain stream segments from complemented image(raster data)?

text: I have raster data in tiff format, i want to extract stream segments from the complement of original raster data. I am able to extract stream segments from original raster data, but not able to extract from complement data. Here complement of data means just reverse the each pixel value, i.e. for each pixel do (max_value_of_pixel - current_pixel_value)


title: Change ArcGIS Server mapservice properties programmatically

text: I am weekly caching my ArcGIS Server mapservice tiles with help of a python script using the arcpy.ManageMapServerCacheTiles_server() command (ArcGIS 10.0, SP5).   For whatever reason, this script fails from time to time. This would not be too bad since I can just restart the script in that case.   BUT after failing, the map service properties are not the same than before: The mapservice property "Maximum time a client can use a service " is now 30,000,000 seconds rather than 60s. I guess ManageMapServerCacheTiles is setting this value at the beginning of the caching and sets it back when finished? Because caching fails, the value is not reset at the end.      Is it possible to have a python script checking this value (or all mapservice properties) and if the value is not correct to stop the mapservice, reset the value and then to restart the mapservice?      Thank you!


title: symbology changes colour when harvested from WMS to online map

text: I'm harvesting a web map service and from my legend the symbology is pink but the same feature on the online map is green? I've had a look at the GetCapabilities file and opened the legend png and it's pink. Does anyone have an idea as to why the colour is changing once visualised on my map? Thanks


title: ArcMap toolbox script with draw geometries parameter

text: I have an arcpy script in toolbox, that need one additional parameter to draw geometry(ies) which will work as a filter. I have done something similar using arcobjects using GPFeatureRecordSetLayerClass as an input parameter and creating temporary layer. If feature set exist, then ArcMap allow to draw geometries in it. Something similar is done if default ArcMap tools "Server Tools>Data Extraction>Extract Data" where feature set "schema" parameter have a dynamic value [toolbox folder path][shapefile] that look something like this               C:\Users\... ...\Desktop10.1\ArcToolbox\My Toolboxes\aoi.shp      When I tried to insert scheme path manually, it tells me that files does not exist. Any ideas to how to replicate feature set behavior from the default "Extract Data" tool?


title: Export Torque Animation

text: Is it possible to export a video of a CartoDB Torque animation to a file, like in arcgis? If not, what are my options for displaying the animation in a presentation. Thanks for the help!


title: Custom UTF Grid Renderer performance improvement help

text: Back in November I asked a question (Create interactive map with no server) and had lots of great responses all of which I tried but because my dataset was so large I wrote my own custom UTF Grid renderer. It worked brilliantly and rendered a large pyramid of UTF Grid files in the projection I needed, however to render almost 6million UTF Grid json files took some time as each one is processed individually. I wondered if anyone could guide me on how best to improve the performance using multiprocessing or threading. I was thinking of calculating all the bounds into a dictionary or a list and then splitting these into a smaller number of lists and then running a process on each one? But not sure if that would work? This is my code which is built from the great code written by Dane Springmeyer, Klokan Petr Přidal and Matthew Perry.               import os     import math     import ogr     import mapnik     try:         import simplejson as json     except ImportError:     import json                    def TileBounds(tx, ty, zoom):         "Returns bounds of the given tile in EPSG:900913 coordinates"              tileSize = 250              minx, miny = PixelsToMeters( tx*tileSize, ty*tileSize, zoom )         maxx, maxy = PixelsToMeters( (tx+1)*tileSize, (ty+1)*tileSize, zoom )              return ( minx, miny, maxx, maxy )          def Resolution(zoom ):         "Resolution (meters/pixel) for given zoom level"         res = 896 / float(2**zoom)              return res          def PixelsToMeters(px, py, zoom):         "Converts pixel coordinates in given zoom level of pyramid to EPSG:27700"         res = Resolution( zoom )         mx = px * res         my = py * res              return mx, my          def MetersToTile(mx, my, zoom):         "Returns tile for given coordinates"              px, py = MetersToPixels( mx, my, zoom)         return PixelsToTile( px, py)          def MetersToPixels(mx, my, zoom):         "Converts EPSG:27700 to pyramid pixel coordinates in given zoom level"              res = Resolution( zoom )         px = mx / res         py = my / res         return px, py          def PixelsToTile(px, py):         "Returns a tile covering region in given pixel coordinates"              tileSize = 250              tx = int( math.ceil( px / float(tileSize) ) )         ty = int( math.ceil( py / float(tileSize) ) )              return tx, ty               def main():         shppath = 'aa_ward_final.shp'         layernum=0         fields=None              ds = ogr.Open(shppath)              layer = ds.GetLayer(layernum)              m = mapnik.Map(250,250)              # Since grids are `rendered` they need a style         s = mapnik.Style()         r = mapnik.Rule()         polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#cccccc'))         r.symbols.append(polygon_symbolizer)         s.rules.append(r)         m.append_style('My Style',s)              ds = mapnik.Shapefile(file=shppath)         mlayer = mapnik.Layer('poly')         mlayer.datasource = ds         mlayer.styles.append('My Style')         m.layers.append(mlayer)              if fields is None:             fields = mlayer.datasource.fields()              bbox1 = 0         bbox2 = 700000         bbox3 = 0         bbox4 = 1344000              minzoom = 0         maxzoom = 2              for tz in range(minzoom, maxzoom):             print " * Processing Zoom Level %s" % tz                  tminx, tminy = MetersToTile( bbox1, bbox3, tz)             tmaxx, tmaxy = MetersToTile( bbox2, bbox4, tz)                  list1 = list(range(tminy, tmaxy))             list2 = list(reversed(range(tminy, tmaxy)))                  counter = 0                  for ty, ty1 in zip(list1, list2):                      for tx in range(tminx, tmaxx):                     counter += 1                     output = os.path.join('utfgrid', str(tz), str(tx))                     if not os.path.exists(output):                         os.makedirs(output)                          tilefilename = os.path.join(output, "%s.json" % ty)                          tilebounds = TileBounds( tx, ty1, tz)                          box = mapnik.Box2d(*tilebounds)                     m.zoom_to_box(box)                     grid = mapnik.Grid(m.width,m.height)                     mapnik.render_layer(m,grid,layer=0,fields=fields)                     utfgrid = grid.encode('utf',resolution=2)                     with open(tilefilename, 'w') as jsonfile:                         jsonfile.write(json.dumps(utfgrid))                  print "Zoom Level %s Tile Count: %s"%(tz,counter)                    if __name__ == '__main__':         main()


title: What are some examples of interesting zoom slider customisations?

text: The default Claro zoom slider control is pretty boring. It's also hard to see on top of many basemaps: ![enter image description here](http://i.stack.imgur.com/6cJpJ.png) Are there any examples of some nice customisations I could use for inspiration?


title: Zoom to selection

text: I am trying to create an application, where in I want to zoom to selected features. Basically, I want to zoom to the extent of a selected point or polygon.


title: Difference between video captured with Google Earth and Google Earth Pro

text: Is there any difference between video which is captured with Google Earth and Google Earth Pro? Maybe some kind of tag showing 'pro'?


title: How to view or export WMS data in tabular format?

text: Is there any way to export or view data from a WMS in tabular format? The info tool brings up a table showing the attributes for the selected feature, but I want to view attributes for all features.


title: How to create an attractive choropleth map in R?

text: I've been experimenting with choropleth mapping techniques in R, having reaped the benefits of ggplot2 for creating beautiful graphs within a powerful data analysis package. When it comes to mapping I've failed to produce results comparable to those I routinely get from QGIS. Using data of flow in Sheffield (fully replicable data and .qgs file here), QGIS easily produced the following: ![QGIS image](http://i.stack.imgur.com/vgAMx.png) The best I could produce using R (using code described here) were these two maps, the first using GISTools, the second using ggplot2: ![GISTools choropleth](http://i.stack.imgur.com/bH3L2.png) ![ggplot2 choro](http://i.stack.imgur.com/jUey4.png) To me, the ggplot2 option seems far more attractive, if only I could solve the problem of the faulty lines (probably a problem with the fortify() command, or not reading in shapefiles using readOGR() described here.) So the question is 2-fold: is the ggplot2 option the best choropleth mapping solution in R and, if so, how can I solve the problem of the faulty white lines? Replicable code to find out what I've done is here. **Edit - Since made choropleth() output option more attractive:** ![New Rplot](http://i.stack.imgur.com/EfthI.png) The ugliness of this option can be reduced by exporting at higher resolution and removing the legend (other R legends are available). Still not managed to change line thickness but it's getting better! Red lines represent flows to zones which employ > 5000 people.


title: Moving cloned maps

text: Is it possible to always have exactly the same area of a map shown on the clone as on the original? As I move around the original map, I want the cloned map to move automatically to show the same area.


title: How to create/add features without geometry?

text: I have a simple SpatiaLite geodb loaded in QGIS with one table. There is a POINT geometry field and some alphanumeric fields. In some cases I can record data about a specific site, but the coordinates are not available at this time and will be retrieved and recorded later. How do I add features to this layer without adding their geometries?


title: graduated labels in QGIS

text: I've searched on Google but haven't found an answer. The QGIS Dufour GIS provide an enhanced vector-styling interface but is it also possible to graduate the labels like it is possible for the markers? Let me say: I have cities with population data: I can have the marker size according to the population value but is it also possible to adjust the lettersize of the label according the population? Best regards, Riccardo


title: Printing a Form includes data grid and axmapcontrol

text: i have created a windows form with some controls in it like data grid and axmapcontrol, it seems that printform cannot print the content of the axmapcontrol. i need to create a table of data and a map referencing the table data, like printing map and attribute table in one page. i'm looking for any suggestions, thanks


title: arcpy - Calculate field in a joined table

text: I've to calculate a field based on another field in a joined table, but an error occurs: **"Parameters are not valid. ERROR 000728: Field !Table_A.field_A! does not exist within table"** My code:               import arcpy     from arcpy import env     env.workspace = "C:/TMP/BASELINE.gdb"     Table_A = "dist"     Table_B= "matrix"     arcpy.AddJoin_management(Table_A, "id_grid", Table_B, "id_matrix")     arcpy.CalculateField_management(Table_A, "!" + Table_A + ".field_A!", "!" + Table_B + ".field_B!","PYTHON")      Both tables are in the BASELINE.gdb. The Join occurs without any problem. I forgot to say that, of course, "field_A" exists in "Table_A" and "field_B" exists in "Table_B". Last but not least, if I run the same calculation in Arcmap, using Field Calculator, no error occurs. Please, any help appreciated. Thanks in advance


title: Adding a layer to ArcMap

text: I have create "Draw polylines tool", but not know how to create spacial layer to store polylines only in this layer. It's possible locate polylines in own layer, I mean that this layer will be creating automatically when I click "Draw polylines tool" ?


title: ArcGIS 10.0 SQL logs?

text: Is it possible to have ArcGIS 10.0 dekstop log the SQL it sends to the database (direct connect with Oracle 11gr2)? I am trying to diagnose some performance issues and would like create some scripts to duplicate what this application is sending to the database. However, you cannot obtain those with a SQL trace (complex datatype). If not - I will capture the network traffic instead, although that less desirable as it is a rather manual process.


title: How to create a multimodal network?

text: How to add waterways and railroads without turn restrictions into the previously created street network dataset with turn restrictions, in order to create a multimodal network? For a university project, we use TeleAtlas Multinet data. We used the following Street Data Processing Tool: http://www.arcgis.com/home/item.html?id=755f96fcde454ece8f790fecb3e031c7 Street Network is created with all kinds of turn restrictions, then now we need to add waterways and railroads into the routable network with global turns. However we couldn't create the multimodal network. Could you please give a suggestion if you have an idea about this topic. Regards


title: How to reproject MODIS Swath data to WGS84

text: Im new to GIS and now working with MODIS04L2 Aerosol Optical Thickness image. I need to reproject the image to other coordinate system such as WGS84 for later work but now Im stuck with it. I have try two ways:   1. Write my own code using python with gdal lib: I follow the tutorial which using GCPs point as paramater for gdalwarp command. The problem is I can't get the GCPs point. The dataset.GetGCPs() function return a zero-length array of GCPs.   2. Using programs like HEG tool, MRT tool but both of it have issues. When I open MODIS file (with extension .hdf) with MRT tool I get the error "Error in Module: ReadHdrFile". I download MODIS from ladsweb.nascom.nasa.gov/data/search.html with parameter: Satellite/Instrument:Terra, Group: Terra Atmosphere Level 2 Product, Products: MOD04L2_Leve 2 Aerrosol and Collection: MODIS Collection 5.1. I want to extract subdataset "mod04:Image_Optical_Depth_Land_And_Ocean" from hdf MODIS file to geotiff which reprojected to WGS84. The sample hdf file http://sdrv.ms/JpVQNh The log from gdalinfo:http://sdrv.ms/JpTQok Anyone here has experiment with this process? Thanks UPDATED: It turns out that the gdal package on my computer is out of date so I cannt get the GCPs via gdalinfo. I updated the gdal to 1.10.0 and now everythings work fine.


title: vb.net arcobjects determine 'create features' on which layer

text: When I create new feature in editing session, I need to determine the new created feature is in which layer because after that have further more action on the layer. I wonder can this be done in vb.net arcobjects? I am using vs2012 and arcmap 10.2.


title: How to link pgrouting edges with street names to route user based in street names?

text: I have installed postgreSql + postGis + pgrouting and I have loaded OSM data into my database using osm2po and shp2pgsql. How can I search for points in the map using the street names and then route the user from point 1 to point 2. pgrouting divides the streets in edges and any given street will have multiple edges and nodes, how could I link those edges and nodes to the actual street names? Any ideas? Thanks


title: Can I use QGIS plugins in PyQGIS as a stand alone script without iface?

text: Can I use a qgis plugin especially (openlayers plugin) in stand alone pyqgis script? I don't find examples, just code using iface but that cannot used in stand alone script.


title: How to find boundaries with Fuzzy logic and GIS

text: How can I use Fuzzy logic in order to find the boundaries of an area with points? Is there a GIS software that does it? I would prefer Arcgis. example from scholar pedia ![enter image description here](http://i.stack.imgur.com/5YkXS.gif) Figure 1: Bird's-eye view on a forest: Where is the boundary of the forest? Which location is in the forest and which is out of it? (image from: http://www.scholarpedia.org/article/Fuzzy_sets#example)


title: Opening geoprocessing tools causes ArcGIS to hang

text: I am getting the same problem as http://forums.arcgis.com/threads/19468-Opening-tools-causes-ArcMap-to-freeze- crash That is, when I open any geoprocessing tool in ArcCatalog (might happen in ArcMap as well, though it hasn't to me yet), sometimes the application hangs. The timing of the hang occurs after the tool dialog opens but before it finishes drawing. I don't know exactly what conditions cause this to happen, but it feels random, maybe in 1 in 4 attempts it hangs. For me, terminating the hung process and trying again sometimes works, sometimes it doesn't. But it's a huge headache. Has anyone come across the same problem and found any resolution? I have already tried reinstalling ArcGIS. I think this started happening some time after installing SP3 for ArcGIS 10.


title: How to run spatial query get fields value of 2 layers?

text: Last month I ask a question "vb.net arcobjects spatial query get fields value of desired layer only", url vb.net arcobjects spatial query get fields value of desired layer only. Now there is some change in the requirement, I need to get value of 2 layers after use "Select Features" tool in arcmap, 1 polyline and 1 polygon. I am not sure may I use the same method? How can I do this? One thing to highlight is the editable variables in attribute table of the 2 layers are the same. My code is here:                   pMxDoc = My.ArcMap.Application.Document         pMap = pMxDoc.FocusMap         pActiveView = pMxDoc.FocusMap              Dim layerNum = GetIndexNumberFromLayerName(pActiveView, "test_db.DBO.TestTable")         Dim layerNum2 = GetIndexNumberFromLayerName(pActiveView, "test_db.DBO.TestTable2")              Dim pFeatureLayer As IFeatureLayer         Dim pFeatureLayer2 As IFeatureLayer              pFeatureLayer = pMap.Layer(layerNum)         pFeatureLayer2 = pMap.Layer(layerNum2)              Dim pFeatureSelection As IFeatureSelection         pFeatureSelection = pFeatureLayer         logger.Info("pFeatureLayer = " + pFeatureLayer.Name)              Dim pSelectionSet As ISelectionSet         pSelectionSet = pFeatureSelection.SelectionSet         logger.Info("pSelection number = " & pSelectionSet.Count)              Dim pFeatureCursor As IFeatureCursor = Nothing         pSelectionSet.Search(Nothing, True, pFeatureCursor)              Dim tempArr As New List(Of String)         Dim pFeature As IFeature         pFeature = pFeatureCursor.NextFeature              Do While (Not pFeature Is Nothing)             tempArr.Add(pFeature.Value(pFeature.Fields.FindField("Item_ID")).ToString)             pFeature = pFeatureCursor.NextFeature         Loop      EDIT: Sorry of the unclear requirement of the whole process. After performing the spatial query, I will retrieve data related to features of these 2 layers (from database) and allow user to select and zoom to any results from the spatial query.


title: location and mapset at Grass-QGIS

text: how can I indicate my drawings (unknown) location(coordinate) to GRASS at QGIS. because some time after importing the other vector or raster layers I hit an error regaring to coordinate issues. Thanks Hamid


title: Refresh leaflet marker data in place ? (on insert to sql db)

text: I have a form that sends data via ajax to be geocoded and inserted into a mysql db. Once that data is there, can I force the marker layer to refresh so that it displays the new data? The SQL API for Cartodb includes a layer method .invalidate() which refreshes the data if there has been a change. I don't care about detecting a change, but I do hope to find a similar refreshing capability in pure leaflet. The map will not be resized or zoomed when the new data is pulled in, so something like map._onResize() won't have any effect. There will be reference URLs etc that need to be pulled from the db to add to a popup, so just putting on the marker on the front end while pushing to the db via ajax won't work.


title: how do i get a GDAL(using gdalwarp) clipped layer projected in pseudo mercator ESPG:3857 to not be shifted?

text: how do i get a GDAL rasterized layer projected in pseudo mercator ESPG:3857 to not be shifted ? i have a vector layer projected in ESPG:3857,after rasterizing using this command:               gdal_rasterize -a_nodata -999 -te xMin yMin xMax yMax -a theField -tr 10 10 -l LayerName -co TFW=YES -co COMPRESS=LZW  VectorLayerPath outPutFilePath      the rasterization is perfect, but when i clip using this command:               gdalwarp -dstnodata -999 -q -cutline theMaskVectorLayerPath -crop_to_cutline -of GTiff inputrasterPath outPutRasterPath      the raster is shifted then clipped and the result is wrong. is this happening to anybody else using pseudo mercator?


title: Min/Max values from a feature class

text: I am very rusty/inexperienced in python and I was hoping someone might be able to give me some pointers on this script. I would like to create a more compact script that prints min and max values of each numeric field from a feature class. Right now I have the code below repeated for each field that contains numeric values. How would I code this to identify if a field has numeric values and then print the min and max values to the python shell or even better to a .csv/text file. Thank you for your advise!               import arcpy     data = "C:/Workfolder/Input.gdb/FeatureClass"          listname = []     cursor = arcpy.da.SearchCursor(data, "field1")          for row in cursor:         listname.append(int(row[0]))          print 'Field1_min, {0}'.format(min(listname))     print 'Field2_max, {0}'.format(max(listname))


title: Use python to execute GRASS module v.distance

text: I have successfully use python script to import a line shapefile and a point shapefile into GRASS, and the following are their attribute table in QGIS,   * line shapefile ![enter image description here](http://i.stack.imgur.com/dkn8F.png)   * point shapefile ![enter image description here](http://i.stack.imgur.com/QsOPl.png) Now I want to find the nearest distance from those points to the line feature. Here's my python code,               import os     import sys     sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))     import grass.script as g     import grass.script.setup as gsetup     gisbase = os.environ['GISBASE']     gisdb = 'C:\Users\Heinz\Documents\grassdata'     location = 'newLocation'     mapset = 'TC'     gsetup.init(gisbase, gisdb, location, mapset)          g.run_command('v.distance', fr = 'tc_station_cor', to = 'tc_bf2500', output = 'connections', upload = 'dist', column = 'dist')     print g.read_command('g.list', _type='vect')      The script ran OK, but the output of vector list has contained no map named connections. I don't know where I made mistakes. I am working with python 2.7.4, GRASS 6.4.3 under Windows 8.1 x64.


title: Why does the Clipper tool lack the mask layer option?

text: Fairly new at QGIS, but everything's gone well so far. I'm running QGIS version 1.8 on a Mac (OSX 10.8) and GDAL 1.10, and other GDAL tools (like Merge and Georeferencer) work just fine. When I try to use the Clipper tool to clip a raster image using a vector shp file, I don't see the dialog options that are shown on various forums (for example, here). There is no "No data value" control and under Clipping Mode I just have the option to clip to a rectangular extent. What am I doing wrong?


title: How to overcome the layer limit on OSX?

text: I am trying to add many layers to a project in QGIS. After adding approximately 80 layers, QGIS returns this error: "not a valid or recognized data source". There is nothing wrong with the shapefiles themselves, they will all open if done individually. It seems to just be after a certain number that it won't allow any more. The error message seems strange to me. Is there a limit to the number of layers a single project can have? Thank you for your help


title: Can't convert geometry from Dgn file to Shape file?

text: I am using below codes to try convert geometry from a Dgn file to Shape file. But, I can't exit the Feature reading loop. Are there anyone know why and how to fix it ? Please download the dgn file from the link: DgnFile               bool Convert()     {              OGRRegisterAll();              OGRDataSource       *poDS,  *poSHPDS;         const char *pszDriverName = "ESRI Shapefile";         OGRSFDriver *poDriver;              poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(                       pszDriverName );              // Open file DGN              poDS = OGRSFDriverRegistrar::Open( sourceFile, FALSE );              if( poDS == NULL )              {              cout << "Error opeing DGN file." << endl;              return false;         }          // Create datasource for file shp       poSHPDS = poDriver->CreateDataSource(destinationPath,NULL);            if( poSHPDS == NULL )            {                  cout << "Error creating shape file." << endl;                  return false;       }                 OGRLayer  *poLayer, *poSHPLayer;             // Get Layer from DGN file            poLayer = poDS->GetLayer(0);            // Create new Layer of SHP file            poSHPLayer = poSHPDS->CreateLayer (shpFileName, NULL, wkbPolygon, NULL);            if(poSHPLayer == NULL)            {                    cout << "Error creating Layer of Shape File." << endl;                    return false;       }            OGRFeature *poFeature;            poLayer->ResetReading();             **while( (poFeature = poLayer->GetNextFeature()) != NULL ) // Can't exit this loop**             {                  OGRGeometry *poGeometry;             poGeometry = poFeature->GetGeometryRef();                  // if feature is a Polygon             if( poGeometry != NULL && (wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon ) )             {                     OGRFeature *poSHPFeature;                     //Create new Feature shape file                     poSHPFeature = OGRFeature::CreateFeature (poSHPLayer->GetLayerDefn());                     poSHPFeature->SetGeometry (poGeometry);                     OGRPolygon *poPolygon = (OGRPolygon*) poSHPFeature->GetGeometryRef();                          poLayer->SetSpatialFilter(poGeometry);                     poLayer->ResetReading ();                   }              OGRFeature::DestroyFeature( poFeature );              }         OGRDataSource::DestroyDataSource( poDS );         OGRDataSource::DestroyDataSource( poSHPDS );         return true;     }       Thanks and regards, Tai


title: kriging Error dimensions do not match: locations xxxx data yyyy in R

text: I've been working (for hours) on a kriging problem. When I run an ordinary kriging model:               y.ok <- krige(y~x, locations=data1,newdata=data2),      I keep getting the message: > Error dimensions do not match in locations xxxx data yyyy I removed any duplicated from the coordinates and that doesn't help either. I checked the various websites and didn't find anything useful. Thanks for any help.


title: DTED/GeoTiff layer: Get elevation under mouse pointer

text: I have a simple desktop app with a map control. The map data that will be available to show will be either DTED files (Military Analyst's catalog / MADTEDLayer), or GeoTiffs (created by gdal). I want the user to be able to continuously see the elevation under the position of the mouse cursor. Is it possible to extract this information from the map layer somehow? PS. I suppose that GeoTiffs have some kind of elevation data embedded since I opened them with GlobalMapper and saw height information under cursor.


title: how to select line segments not overlapped by other layer's lines

text: I have two shapefiles of forest roads (mostly dirt) from different sources, covering the same area. The goal is to build one shapefile of roads that is the composite of the other two. Layer A is less complete, and Layer B is more complete, i.e. if a dirt road on A is a mile long, the same road may exist on B but is a mile and a half long, and the first mile mostly overlaps with the A line but not exactly. The trivial answer here is "just take the B road and ignore the A road" but there are cases that are reversed from the example, i.e. there are some cases where A has more detail / more roads than B.   * how would you do overlap detection? Buffer each layer and look for overlaps?   * how would you select only the portions of lines that are not in the overlapped region, for example, the last half mile of the B road from the example? Thanks


title: How many kinds of "address" are there?

text: I'm defining "address" as "a physical location, usually visibly numbered". **Numbered** was used in the definition to exclude purely geographical features not usually considered to have or need an address. **Visibly** was used in the definition to exclude entities that might be numerically coded (like census features) but the number is not typically apparent nor used for directional or delivery purposes. **Usually** was used in the definition because exceptions abound (such residences in Carmel). Enumerating the possible kinds of addresses resulted in:   * buildings.   * entire floors of buildings.   * apartments/offices/rooms in buildings.   * cubicles/mailstops in buildings.   * loading docks.   * slips on boat docks.   * slots in trailer parks.   * post office boxes and private mail drops.   * military APO/FPO addresses.   * jail/prison cells.   * campground lots.   * outbuildings (not usually visibly numbered, but often desireable in delivery databases).   * undeveloped parcels (not usually visibly numbered, but often required in address databases). What other kinds of "physical locations, usually visibly numbered" have I forgotten to include in the above list? Two questionables are:   * hospital beds.   * cemetary plots. * * * Update from comments: The enumeration (ie, determining how many codes are necessary) is desired to allow coding address kinds to differentiate records in an address database. Eg, SELECT * FROM ADR WHERE KIND IN (1,3,5). "Usually visually numbered" (alphanumerics are OK) is preferred to just "coded" as described in the original question. If the location isn't visually numbered, it's unlikely to be the target of directions/deliveries, and would not be kept in the address database. Just think of the question as "how many different kinds of places do we physically stick a number on to make directions/deliveries easier?"


title: What are some good resources related to Geodjango?

text: I am trying to learn Geodjango. Has anyone used it in their own project? The documentation on the official site is sparse. Can someone point out the best resources helpful in learning it other than official documentation ? Please also tell me how is your experience with it and other web application frameworks that you have used for mapping.


title: ArcMap: extracting overlapped buffers

text: I have a set of points in ArcMap which I create a buffer for each point. Then I use intersect to extract the overlapping buffers. The result is something like this: ![enter image description here](http://i.stack.imgur.com/r4Hk0.png) The lighter color are the overlapped layer. Weird thing is that on the left, the overlapped buffer is rather incorrect, is it? Any idea what is going on?


title: How to set code for designer legend in mapbasic

text: last couple of days i am doing work on thematic maps using mapbasic. I am getting error while set the legend code. it giving me unrecognized command or code kind of error ... my code is               Create Designer Legend From Window nMID      Portrait Default Frame Title "# Legend" Font ("Arial",0,10,0)     Default Frame Style "%" Font ("Arial",0,8,0)     Frame From Layer 2 Using column object label default Frame From Layer 3 Using column object label default Frame From Layer 1      it giving error for Create Designer legend..... does anybody knows this.....??? BTW...happy new year to all


title: Converting Ordnance Survey raster maps to WGS84 webmap tiles

text: I would like to generate Google/Bing-like map tiles from Ordnance Survey's 1:25000 and 1:50000 raster maps, which are provided as 4000x4000 UK national grid tiles (GeoTIFF). A sample tile can be found at http://www.ordnancesurvey.co.uk/docs/sample- data/25k-raster-sample-data.zip . Thanks to this site and others (this is my first time using proper GIS tools), I have succeeded in using gdalwrap to reproject the tiles, and gdal_retile.py to cut the reprojected tile into smaller pieces. I still have the following questions:   1. the source files use a color map, so based on my single-tile sample, I converted that to 24bit color in order to be able to use a better resampling algorithm than nearest neighbor, but I don't know how scalable that will be to the whole of the UK. Is there a way to get a quality resampling without first converting to 24bit color?   2. how do I get gdalwrap or gdal_retile.py (or another utility) to align the generated tiles to the webmap grid (the world divided into 2^zoom tiles of equal lat-lon delta)? Right now my tiles are the right pixel size, but most likely not aligned to the lat-lon grid nor the zoom level   3. the tiles generated by gdal_retile.py are named according to each one's location in the source file, but I need to translate that to the global tile number for my zoom level. Do the tools support doing that on the fly, or do I need to post-process?   4. using gdal_retile.py with the -levels option generates sub-directories but the tiles therein are all black; I was expecting the zoomed out pyramid (and fewer tiles, not more, than the default level), should I be using a different option?   5. even for this smallish sample, each command seems to take an inordinately long time, and this is on an octo-core MacPro and SSD. Is there a performance guide for GDAL? Here's what I have been doing, starting with the sx99.tif sample file:               pct2rgb.py sx99.tif sx99_24b.tif          gdalwarp -t_srs EPSG:3857 -r lanczos -wm 150 -multi -of GTiff sx99_24b.tif sx99_24b_warped.tif          gdal_retile.py -of png -levels 1 -targetDir tiles sx99_24b_warped.tif       * * * Instead of gdal_retile.py, I've tried using gdal2tile.py on the pre-warped raster maps. However, while the X coordinate for each of the generated tiles is correct, the Y coordinate seems to be at least double what I am expecting. I've tried explicitly setting the source SRS and target profile, to no effect.               gdal2tiles.py sx99_24b_warped.tif      I can't understand what's going wrong. The warped intermediate image looks OK:               $ gdalinfo sx99_24b_warped.tif      Driver: GTiff/GeoTIFF     Files: sx99_24b_warped.tif     Size is 4074, 4085     Coordinate System is:     PROJCS["WGS 84 / Pseudo-Mercator",         GEOGCS["WGS 84",             DATUM["WGS_1984",                 SPHEROID["WGS 84",6378137,298.257223563,                     AUTHORITY["EPSG","7030"]],                 AUTHORITY["EPSG","6326"]],             PRIMEM["Greenwich",0],             UNIT["degree",0.0174532925199433],             AUTHORITY["EPSG","4326"]],         PROJECTION["Mercator_1SP"],         PARAMETER["central_meridian",0],         PARAMETER["scale_factor",1],         PARAMETER["false_easting",0],         PARAMETER["false_northing",0],         UNIT["metre",1,             AUTHORITY["EPSG","9001"]],         EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],         AUTHORITY["EPSG","3857"]]     Origin = (-396511.819462091429159,6584391.706861065700650)     Pixel Size = (3.948971656093984,-3.948971656093984)     Metadata:       AREA_OR_POINT=Area     Image Structure Metadata:       INTERLEAVE=PIXEL     Corner Coordinates:     Upper Left  ( -396511.819, 6584391.707) (  3d33'42.93"W, 50d47'27.29"N)     Lower Left  ( -396511.819, 6568260.158) (  3d33'42.93"W, 50d41'57.18"N)     Upper Right ( -380423.709, 6584391.707) (  3d25' 2.66"W, 50d47'27.29"N)     Lower Right ( -380423.709, 6568260.158) (  3d25' 2.66"W, 50d41'57.18"N)     Center      ( -388467.764, 6576325.932) (  3d29'22.80"W, 50d44'42.32"N)      But the resulting files are in the range of (for zoom level 11): 1003 <= x <= 1004 and 1359 <= y <= 1360, where I was expecting y around 780...


title: How to choose spatial reference for GPS data?

text: We received GPs data from the vendor. I want using ArcGIS software to create XY event table and create a point feature class. What the spatial reference would be for the GPs data? Is it something I should ask the vendor about? Is it from the GPs device?


title: Uploading maps to GPS (forestry maps: roads, trails, hydro, polygons)

text: I'm using a QGIS 2.2.0 Valmiera on a Windows 8.1 system Desktop. I'm a forester and I would like to have a cheap solution to load maps into the handheld Garmin 60, 62 or 64 series GPS. Free would be even better, if possible. There is a software that costs around $350 that creates an image (.img) that is easily loaded into GPS, but it's expensive AND I don't think it works with QGIS. The best solution I've found so far would be Mapwel, which costs $45 or $90. You just take the shapefiles you need and merge them, then load them to your GPS. Does anyone know of a cheaper solution?


title: Updating SDE point feature classes with SQL SERVER 2008 R2

text: I am contacting you in hope that you can help me resolve an issue using SQL Server 2008 R2 to push updates to a SDE point feature class. We are using a C# app that uses ArcGIS Runtime to edit the attribute table and create new address records with their cooresponding points. The problem is we can see the updates in ArcMap, the new records are showing up, and when you click on the record the place where the point should be is highlighted, however it doesn't graphically display the point in ArcMap. I have been working with ESRI and they can't seem to figure out the problem. So I thought I would reach out and see if you have any suggestions.


title: What is the best way to represent bridges using a point vector layer in QGIS?

text: I'm working on a hiking/biking maps of trails in the woods which sometimes have bridges that need to be shown on the map. I'm trying to figure out the best way to capture and represent them using QGIS and a point shapefile. I'm capturing them in a general "Features/Markers" layer where I also store the location of things like information booths and notable trail features/locations. I have a attribute named "type" that I use to style each type of feature ("bridge", "information", "feature", etc.) and a "name" field I can use to optionally label the feature. Since bridges point in a certain direction I added a "rotation" field which I'm referencing in "Style, Advanced, Rotation Field". Finally I'm using "Map Unit" as my unit so the size of the bridge scales with the size of the map. I'm considering also adding a "size" field to the "feature" layer as well in case I determine different bridges need to be different sizes. I'm curious how others capture/style bridges and if there's anything I'm not considering? Here's a screenshot of how it looks: ![enter image description here](http://i.stack.imgur.com/hYwLU.png)


title: How to zoom to point boundary using OpenLayers?

text: I'm wondering what the best strategy to position a map that has two points, using OpenLayers is. At the moment, my approach is to calculate the mid-point of the two points, and to zoom the map to that mid-point using the following approach:               var centerPoint = new OpenLayers.LonLat(lon, lat);           map.setCenter(centerPoint, zoomLevel);       However, this is not ideal as I don't know what the `zoomLevel` should be, so sometimes although the map is centered, the points are not visible. Is there a better approach to dynamically position the map?


title: How do you visualize change in rasters over time in CartoDB?

text: I have a series of satellite image rasters. I need to show the change over time of these rasters. Would the best way be to change the basemap with the click of a button for each? Has anyone got an example of doing this in cartoDB that includes how it was done?


title: Best way to 'join' polygons

text: I'm looking for the best ways to 'group' a budge of block polygons in a new polygon that should be the town. My block layer is something like this: ![enter image description here](http://i.stack.imgur.com/zuvZn.png) Or zoomed: ![enter image description here](http://i.stack.imgur.com/w8c0J.png) And I'm looking to get something like this polygon (the left one): ![enter image description here](http://i.stack.imgur.com/ctUT7.png) Could somebody tell me some indications of how to get something like this? I'm not looking for the best accuracy, but something close to town or city polygons. I thought in a convez hull, but the problem is that the block layer are a lot of polygons, instead of an irregular one. I'm using QGIS and PostGIS for this project, and I'm sure that postGIS has some process that approximates to my goal, but I can't find it. Thanks!


title: What prior art can be cited to refute Apple's claim to have invented Interactive Map?

text: In U.S. Patent Application 20130339891 Apple Inc. claims "An interactive capability enables a user to dynamically adjust the content of an electronic map" as well as "Geospatial applications can be integrated with the map to provide information pertaining to the area displayed on the map, and to refine search results to those that are relevant to the area." Apple seeks to patent these ideas based on a number of claims. The claims the patent application makes seem overly broad and applicable to many GIS capabilities; starting with "A method for displaying a map on a computing device, Comprising: storing information to be displayed on the map in a memory of the computing device, the stored information comprising a plurality of different layers of information, wherein each layer contains a respective type of information." Surely there is prior art in many GIS systems that can be cited to refute this claim. Please do so.


title: How to extend polylines to the nearest point

text: Hi could anyone tell me a way of getting my polylines to search for a point within a certain number of meters of their start/end points and extend to it? edit: There was some broblem uploading a picture so I'll explain further. I have contours with gaps like here and am able to close them, but have found that other gaps, caused by edges to the area and faults are also being filled, the lines also travel completely straight when this is not quite the case. I also have a point feature class of where annotations used to be. The gives one more point when the gap should be slightly curved, and identifies the gaps that I actually want to fill. So I want to do something similar to closeing rings, but instead of extending the lines to meet eachother I want to extend them to meet the point.


title: Error while retrieving a WCS coverage with GDAL

text: I'm currently using the OSGeo-Live 7.0 (Link) image to write a PyWPS process. Part of this process is to retrieve a WCS layer with GDAL (from a GeoServer instance). Since GDAL supports WCS out of the box I only need to write the XML description file `wcs.xml`:               <WCS_GDAL>         <ServiceURL>http://localhost:8082/geoserver/sf/wcs?</ServiceURL>         <CoverageName>sf:sfdem</CoverageName>         <Version>1.1.0</Version>         <PreferredFormat>image/tiff</PreferredFormat>       </WCS_GDAL>      The `sf:sfdem` coverage is a demo layer (available in the GeoServer of the OSGeo-Live 7.0 image) which is used as example here.   But if I try to read the WCS coverage with `gdal_translate` do following error occurs:               user@osgeolive:~$ gdal_translate wcs.xml sfdem.tif       Input file size is 635, 635       0ERROR 1: Returned tile does not match expected configuration.       Got 634x477 instead of 635x635.        This is the case with every coverage I've tried. So from my point of view this should be some kind of bug in either GDAL or GeoServer (probably with the bounding box calculations). Does someone else encountered this strange error too or has an ideo how to fix it? The software versions in my OSGeo-Live 7.0 image are:   * GDAL: 1.7.3-6ubuntu3   * GeoServer: 2.3.4 (running on OpenJDK 1.7.0_25) Thanks


title: openlayers find feature by attribute

text: //ilce kontrol                 ilcecontrol = new OpenLayers.Control.GetFeature({ protocol: OpenLayers.Protocol.WFS.fromWMSLayer(ilce) });                 ilcecontrol.events.register("featureselected", this, function (e) {                     alert(e.feature.fid);                     select.addFeatures([e.feature]);                 });                 ilcecontrol.events.register("featureunselected", this, function(e) { select.removeFeatures([e.feature]); });                           map.addControl(ilcecontrol);                 ilcecontrol.activate();                 for(var f=0;f<ilcecontrol.features.length;f++) {                   if(ilcecontrol.features[f].attributes.ADI == 'TALAS') {                     selectFeatureControl.select(ilcecontrol.features[f]);                     break;                   }                 }      hi everyone, my code is like that. i want to find feature.ADI == 'TALAS' and select but its not working and there is no error. i dont understant why


title: Most efficient way of passing values from OGRGeometry to ModelField in geodjango

text: I am writing a custom parser that maps shapefile data to a geodjango model in a Django 1.6 project. To read the shapefile I use the `DataSource` object from `django.contrib.gis.gdal`, which provides access to the shapefile data through `OGRGeometry` objects. In my loop I want to map the geometries into my model as follows               dt = DataSource('path/to/my/shapefile')     lyr = dt[0]     for feat in lyr:         MyModel.objects.create(             geom = feat.geom.FOO         )      My question is, what is the most efficient conversion for `FOO` in that loop? The feature geometry can be extracted in several ways, for instance wkt, wkb, hex, etc.


title: How to model money cost as an accumulator in network analysis?

text: This is not necessarily a technical coding question, although any links that include coding are welcome if they are answering my question. I am looking for research or analysis that has been done in network analysis using the **money cost** of travel as a cost accumulator per segment of road/rail/bus route/etc. Does anybody know of a good resource such as articles/publications/blogs/journal articles showing someone who either does it, suggests how to do it, or explains the methodology behind it? I am looking into the possibility of using money cost as a cost accumulator in the **ArcGIS 10.0/10.1** network analyst least-cost path problem. My thoughts are, if a commuter rail were to be built between two metropolitan centers, which work commuters would be likely to choose the commuter rail over driving a car the whole way when they are factoring in how much it costs them to make the trip?


title: Parsing list object file paths for file name

text: I am trying to automate a naming process that takes the file name from a list of file paths to name other folders, geodatabases, zipfiles, etc. The problem I'm running into is that, I assume, I am trying to run string functions on list object and I can't figure out how to get the script to do what I need it to do. For example, I have the file path C:\TEMP\Durham.shp. I want to name a geodatabase, folder, and zipfile "Durham". Instead I get a folder and geodatabase named "C". Can someone take a look at my script and tel me what I'm doing wrong? The following script is that last incarnation of several attempts.               import arcpy, os     import zipfile     from arcpy import env     from os.path import basename          BASELAYERS = arcpy.GetParameterAsText(0)     CLIPLAYERS = arcpy.GetParameterAsText(1)     OUTFOLDERMASTER = arcpy.GetParameterAsText(2)     FILETYPES = ["*.shp", "*.dbf", "*.shx", "*.prj", "*.sbn", "*.sbx", "*.xml"]          env.workspace = OUTFOLDERMASTER     for LAYERS in CLIPLAYERS:         LAYERSSTR = ''.join(LAYERS)         BASENAME = LAYERSSTR.split("\\")[-1]         FOLDER = os.path.join(OUTFOLDERMASTER, BASENAME)         os.makedirs(FOLDER)         PGDBNAME = "NCFlood_Effective_" + BASENAME + "_PGDB.mdb"         arcpy.CreatePersonalGDB_management(FOLDER, PGDBNAME, "9.3")         PGDB = os.path.join(FOLDER,PGDBNAME)         for FILES in BASELAYERS:             FEATURENAMESTR = ''.join(FILES)             FEATURENAME = os.path.basename(FEATURENAMESTR)             SHAPEFILENAME = FEATURENAME + ".shp"             SHAPEFILE = os.path.join(FOLDER, SHAPEFILENAME)             arcpy.Clip_analysis(FILES, LAYERS, FEATURENAME)             arcpy.Clip_analysis(FILES, LAYERS, SHAPEFILE)         ZIPFILELIST = []                 for TYPE in FILETYPES:             ZIPFILELIST.extend(glob.glob(FOLDER))         ZIPFILENAME = os.path.join(FOLDER, "NCFlood_Effective_" + FOLDERNAME + "_SHP.zip")         ARCHIVE = zipfile.ZipFile(ZIPFILENAME, "w")         for ZIP in ZIPFILELIST:             ARCHIVE.write(ZIP, os.path.basename(ZIP), zipfile.ZIP_DEFLATED)         ARCHIVE.close()      I have also tried                   BASENAME = os.path.basename(LAYERS)         FOLDERNAME = os.path.splitext(BASENAME)[0]


title: Send WMS Request as JSON

text: I'm using openLayers and geoserver. I'm making a layer with openLayers.Layer.WMS to geoserver and am adding cql_filter to filter down id's. Here is my code:                       var kent = new OpenLayers.Layer.WMS(                 "Kent", "http://localhost:8080/geoserver/TK/wms",                 {                     layers: 'TK:phoenix',                     transparent: true,                     format: 'image/png',                     queryVisible: true,                     cql_filter: ("id IN (85356,85357,85358,85359,85360,85361,85362,85363)")                  }             );      I need to be able to send a large list of id's to the request but it makes the url too long. Is there a way I can make this request with JSON?


title: Shortest distance from point to line (big dataset)

text: What is the best method to calculate the distance from the nearest road for a series of address points in QGIS, having a very big dataset? (QGIS 2.0.1, Win XP) I read the answers to previous similar questions: Measuring the distance between lines and points in QGIS How to draw perpendicular lines in QGIS? Calculating the minimum distance between points and polygons in QGIS http://stackoverflow.com/questions/16448402/distance-of-point-feature-to- nearest-polygon-in-r Nevertheless: 1) MMQGIS plugin find the disance to the nearest node of a line network, not the shortest (perpendicular) distance to a line. 2) R and Python proposed scripts correctly works for small datasets and can be implemented in the modeler, however for my job I have 530 points and a regional street layer with about 135'000 features, and cycling through all the features then taking the shortest distance for all points will take an enormous amount of time 3) I have problems in making GRASS work in QGIS 2.0.1 (GRASS in QGIS not working (windows XP)).


title: Where can I get an ocean shaded relief for Central America?

text: Does anyone know where I can get an ocean shaded relief for the Central American coast?


title: GPS coordinates to neighboorhood?

text: I have tried both Google and Yahoo's Geocoding services to get the Montreal neighborhood/district from GPS coordinates but while they are both able retrieve the correct (address, city, state, country) they won't retrieve the neighborhood. Which Geocoding service would offer this service? Otherwise, where can I find the neighborhood data to do it myself?


title: How to work ArcSDEQuerier?

text: I try the ArcSDEQuerier transformer but it doesn't work when i specified the attribut : name ='TEST' for the search (Where). How does this transformer work ? An idea ? Thanks Julien


title: query raster cell values by point

text: I have a binary raster (0s and 1s) that I need to query given a set of point coordinates. This is easy enough to do with gdal (`gdallocationinfo`), but where it get's complicated is that I also need to return the sum of the 8-cell neighborhood. Is there an easy way to do this without processing the entire map? In GRASS the only way I can think to do this is to use `r.mapcalc` to return a sum of the cell I'm interested in and then change the value of all other cells to 0. This allows me to compute an aggregate sum for the map with `r.sum`. But this requires a great deal of processing and if I need to do this for hundreds of points, the task becomes enormous. Ideas?


title: Why isn't map service displaying in ArcCatalog or rest services?

text: Update: I have an issue where my map services, which pull ortho-imagery from a cache, intermittently do not display in my services directory, ArcCatalog or my web portal. Currently I go into ArcCatalog and click Service Properties > Caching > and check the box for Create tiles on demand, the map service will display in my rest directory for one of my services. If I check that Create tiles on demand box for another map services, I cannot see the services in the services directory inside the firewall, but I can see it outside of firewall. This behavior is pretty strange. If anyone has dealt with a similar issue and has some idea on resolution, I would appreciate hearing from you. Yes, I am also working with Esri on this issue. Original: I have several map services that are not displaying the image of the service in ArcCatalog, nor are the Rest Services displaying. These map services contain ortho-imagery, which are cached at several scales, and I am certain I have been able to preview them in ArcCatalog and Rest before. In fact, I can still see the correctly cached images on my Production rest services directory. I can see coordinates when I move my mouse about in the preview window of ArcCatalog of the troublesome map services. Does anyone have a suggestion of how I can preview my map services in ArcCatalog and Rest? I have tried:   * republishing service   * restarting server   * changing map service to display imagery dynamically rather than cached   * altering the output directories I am using ArcGIS 10.0 and my servers are Windows Server 2008 R2 64-bit. Thanks & regards.


title: How do you calculate lat/lons from a PostGIS table which has a projected coordinate system?

text: I have a PostGIS table in a projected coordinate system (EPSG: 28355 or MGA94 Zone55). In QGIS (with the PostGIS manager plugin) I calculate a coordinate by performing the following SQL "UPDATE gisprod.buildings SET easting = x(centroid(the_geom))". Now I want to add a "lat" and "lon" field to this table and calculate the values using the appropriate coordinate system, in this case that would be (EPSG:4283 or GDA94). Is there a way I can specify in the SQL statement to calculate the x/y coord using a different coordinate system? Thanks Andrew


title: php_mapscript.so not present after Mapserver 6.2.0 compilation

text: I use this to compile mapserver               cd /opt/     wget http://download.osgeo.org/mapserver/mapserver-6.2.0.tar.gz     tar xvf mapserver-6.2.0.tar.gz     cd mapserver-6.2.0     ./configure --with-ogr=/usr/local/bin/gdal-config --with-gdal=/usr/local/bin/gdal-config --with-httpd=/usr/sbin/httpd --with-wfs --with-wfsclient --with-wmsclient --enable-debug --with-curl-config=/usr/bin/curl-config --with-proj=/usr/local --with-tiff --with-gd=/usr/bin/gdlib-config --with-jpeg --with-freetype=/usr/include --with-threads --with-wcs --with-xml2-config=/usr/bin/xml2-config --with-postgis=/usr/local/pgsql/bin/pg_config --with-geos=/usr/local/bin/geos-config --with-sos --with-agg --with-agg=/usr/local --with-php=/usr/bin/php-config     make     make install     ldconfig      I work on Centos 6.3 x64 After the compile i do not find the php_mpscript.so in the right location (./mapscript/php). This affect my GeoMose install, giving me 500 error for apache in accesing the geomose .php files. I think this is because MapServer is not compiled. Any clue? Thanks!


title: How do I use my own basemap with Javascript API?

text: I'm trying to make a map using the Javascript API. I want to use my own basemap, but I can't find the most simple code to do this. I always see something like this:                 function init() {         map = new esri.Map("mapDiv", {           basemap: "satellite",           center: [-97.395, 37.537],           zoom: 11         });      But I don't want to use the ArcGIS online basemaps. What is the code to make a map service I am hosting my basemap? Thanks.


title: Convert image in KMZ to polygon

text: I have a KMZ file that contains several image overlays (RF coverage plots). I would like to convert the image overlays to polygons to decrease loading time in Google Maps. Currently I'm tracing a rough outline of the image but that is labor intensive and not very accurate. Is there a way to automatically convert a image to a polygon? I have Google Earth and Manifold GIS. ArcScan looks promising but it is outside of my budget.


title: Error: Exception Raised on the "for row in rows" command line?

text: I am comparing fields in layers. I have the AADT layer and AADTAnno layer and I want it to search the fields FLAG and TFLAG. For Example: If there is a 29C9 in AADT layer in FLAG fields then there needs to be the same thing in the ANNo layer in the TFLAG field. My code stops running on my for row in FRows command line. I think its because I'm not letting it know that FLAG is the field. Please help!! New Code!!               import arcpy, traceback          mxd = arcpy.mapping.MapDocument(r"K:\TASS\2 - GEO-DATA PROCESSING SUPPORT\MICHELLE'S WORK_ENTER NOT!!\Work Folder\Python Programming\Wayne's Tools\Abilene_Base_Map.mxd")     lstLayers = arcpy.mapping.ListLayers(mxd)          flayer = arcpy.mapping.ListLayers(mxd, "AADT")[0]     alayer = arcpy.mapping.ListLayers(mxd, "AADTAnnoLabel")[0]          FRows = arcpy.SearchCursor(flayer)     ARows = arcpy.SearchCursor(alayer)          ffields = arcpy.ListFields(flayer, "", "FLAG")     afields = arcpy.ListFields(alayer, "", "TFLAG")          FList = []     AList = []          for row in FRows:         Fvalue = row.getValue("FLAG")         FList.append(str(Fvalue))          for row in ARows:         Avalue = row.getValue("TFLAG")         AList.append(str(Avalue))          matched = set(FList) & set(AList)          for x in matched:         exp = "ID = " + x         arcpy.SelectLayerByAttribute_management(flayer, "ADD_TO_SELECTION", exp)         arcpy.SelectLayerByAttribute_management(flayer, "SWTCH_SELECTION")


title: Python code for locking elements in Layout view

text: Is it possible to lock elements in layout view in Python? I know it is possible in the production cartography extension from ESRI but all I need is the funcionality to lock elements. I am wondering if there is a way to do this in Python.


title: From GCS_WGS_1984 to UTM WGS 84 zone 33N

text: I was doing my project and i had to transform shapefile from bessel to gcs_wgs_1984. that went fine. Now I have to put it again to projected coordinate system. Please could someone help me transform it from WGS_1984 to UTM wgs_84 zone 33N? Thank you in advance


title: Multiple input tables for Add Field in ModelBuilder

text: My model requires `field1` to be added to `fc1` and `fc2`, and then calculated according to a Python expression `calc1`. `field2` is then added to `fc2`, and calculated using `calc1`. How do I add the multiple feature classes as the input to Add Field and Calculate Field? I have tried creating an input table variable as a list of values, and while it accepts `fc1` and `fc2` when I type it in, the model does not run ("dataset does not exist"). **edit** I'm adding these fields so that I can join the tables on a common pre-existing field and populate `fc2.field1` with the values from `fc1.field1`. ![model screenshot](http://i.stack.imgur.com/gWAzS.png)


title: How to start QGIS after building & install on ubuntu 12.04 64bit?

text: I am a newbie here. I am running ubuntu 12.04 64bit, and build after the instruction, but i don't know why can not launch qgis. Seem i successfully build $ install it. below is the linux shell:               jack@jack:~/apps/bin$ sudo su     root@jack:/home/jack/apps/bin# qgis     The program 'qgis' can be found in the following packages:      * qgis      * qgis-plugin-grass     Try: apt-get install <selected package>     root@jack:/home/jack/apps/bin# ls     qbrowser  qgis  qgis_bench     root@jack:/home/jack/apps/bin#


title: Best Landsat band for forest classification?

text: Using Qgis and the Semi-Automatic Plugin for Classification, which are the best Landsat bands to use to classify trees species (in central Argentina)? In Landsat 5, 7 and 8?


title: Basic MapServer query of PostgreSQL/PostGIS database

text: I've just installed MapServer on a linux server and I am trying to understand how to display a basic map using a PostGIS database. I have displayed a shapefile in the browser, so I know that the mapserver setup works. My mapfile is as follows:               MAP                  NAME "sample"             EXTENT -180 -90 180 90 # Geographic             SIZE 800 400                  LAYER               CONNECTIONTYPE postgis               NAME "roads"               CONNECTION "user=username password=**** dbname=databasename host=hostdb.com port=5432"               DATA "geom from public."MLSOA""                STATUS ON               TYPE LINE               CLASS                     STYLE                       COLOR 0 0 0                     END               END             END               END # end of map file      I then call this by using a url string:               http://mysite.com/mapserv.cgi?map=/home/user/webapps/htdocs/mapfile_pg.map&layer=roads&mode=map      However, I get the error:               loadLayer(): Unknown identifier. Parsing error near (MLSOA_London):(line 11)       I have also tried calling my database just `MLSOA` and `'MLSOA'`. If I use the line `DATA "geom from MLSOA"` I get the following error:               msDrawMap(): Image handling error. Failed to draw layer named 'roads'. msPostGISLayerWhichShapes(): Query error. Error (ERROR: relation "mlsoa" does not exist ) executing query: select encode(ST_AsBinary(ST_Force_2D("geom"),'NDR'),'hex') as geom,"gid" from MLSOA where geom && GeomFromText('POLYGON((-180.225563909774 -90,-180.225563909774 90,180.225563909774 90,180.225563909774 -90,-180.225563909774 -90))',find_srid('','MLSOA','geom'))       I think that the problem is with the table name, but I'm not sure how it should be named. Using pgAdminIII, if I want to run a query on this table I use the following syntax:               select getsrid(the_geom) from public."MLSOA" limit 10      What is the correct syntax for my mapfile?


title: How to print ActiveView to printer in ArcEngine 10

text: I am following a printing tutorial, but it's not working. The tutorial is here: http://edndoc.esri.com/arcobjects/9.2/NET/ViewCodePages/025b5da7-9b10-4623-a51e-c038c348ce29.PrintActiveViewCS_NetPrintActiveViewCS_Net.cs.htm Maybe printing in ArcEngine 10 is different. The initial problem is that the PageLayout is null.               //I tried this like in the sample     docPaper.Orientation = m_hookHelper.PageLayout.Page.Orientation;          //Then I tried passing in my own ActiveView obj and casting to PageLayout      IPageLayout pageLayout = docActiveView as IPageLayout;      But in both cases they are null. Does anyone have a good link for printing in ArcEngine 10?


title: ST_Intersection slow query

text: I am trying to perform an intersection between two layers:   1. Polyline layer representing some roads (~5500 rows)   2. Polygon layer representing irregularly shaped buffers around various points of interest (~47,000 rows) Ultimately, what I'm trying to do is to clip the polylines to these many (sometimes overlapping) buffers, and then sum up the total length of roadway contained within each buffer. The problem is that things are running SLOW. I'm not sure how long this should take, but I just aborted my query after > 34 hours. I'm hoping that someone can either point out where I've made some mistake with my SQL query, or can point me to a better way of doing this.               CREATE TABLE clip_roads AS          SELECT        ST_Intersection(b.the_geom, z.the_geom) AS clip_geom,       b.*          FROM        public."roads" b,        public."buffer1KM" z          WHERE ST_Intersects(b.the_geom, z.the_geom);               CREATE INDEX "clip_roads_clip_geom_gist"       ON "clip_roads"       USING gist       (clip_geom);                    CREATE TABLE buffer1km_join AS          SELECT       z.name, z.the_geom,       sum(ST_Length(b.clip_geom)) AS sum_length_m          FROM       public."clip_roads" b,       public."buffer1KM" z          WHERE       ST_Contains(z.the_geom, b.the_geom)          GROUP BY z.name, z.the_geom;      I do have a GiST index created for the original roads table, and (just to be safe?) create an index before doing the second table creation. The query plan from PGAdmin III looks like this, though I'm afraid I don't have much skill in interpreting it:               "Nested Loop  (cost=0.00..29169.98 rows=35129 width=49364)"     "  Output: st_intersection(b.the_geom, z.the_geom), b.gid, b.geo_id, b.address_l, b.address_r, b.lf_name, b.lfn_id, b.lfn_name, b.lfn_type_c, b.lfn_type_d, b.lfn_dir_co, b.lfn_dir_de, b.lfn_desc, b.oe_flag_l, b.oe_flag_r, b.fcode_desc, b.fcode, b.fnode, b.tnode, b.metrd_num, b.lo_num_l, b.lo_n_suf_l, b.hi_num_l, b.hi_n_suf_l, b.lo_num_r, b.lo_n_suf_r, b.hi_num_r, b.hi_n_suf_r, b.juris_code, b.dir_code, b.dir_code_d, b.cp_type, b.length, b.the_geom"     "  Join Filter: _st_intersects(b.the_geom, z.the_geom)"     "  ->  Seq Scan on public."roads" b  (cost=0.00..306.72 rows=5472 width=918)"     "        Output: b.gid, b.geo_id, b.address_l, b.address_r, b.lf_name, b.lfn_id, b.lfn_name, b.lfn_type_c, b.lfn_type_d, b.lfn_dir_co, b.lfn_dir_de, b.lfn_desc, b.oe_flag_l, b.oe_flag_r, b.fcode_desc, b.fcode, b.fnode, b.tnode, b.metrd_num, b.lo_num_l, b.lo_n_suf_l, b.hi_num_l, b.hi_n_suf_l, b.lo_num_r, b.lo_n_suf_r, b.hi_num_r, b.hi_n_suf_r, b.juris_code, b.dir_code, b.dir_code_d, b.cp_type, b.length, b.the_geom"     "  ->  Index Scan using "buffer1KM_index_the_geom" on public."buffer1KM" z  (cost=0.00..3.41 rows=1 width=48446)"     "        Output: z.gid, z.objectid, z.facilityid, z.name, z.frombreak, z.tobreak, z.postal_cod, z.pc_area, z.ct_id, z.da_id, z.taz_id, z.edge_poly, z.cchs_0708, z.tts_06, z.the_geom"     "        Index Cond: (b.the_geom && z.the_geom)"      Is this operation just doomed to run for several days? I'm currently running this on PostGIS for Windows, but I could in theory throw more hardware at the problem by putting it up on Amazon EC2. However, I see that the query is only using one core at a time (is there a way to make it use more?). Many thanks for any help that you can provide!


title: Does anyone know how to install the Orfeo Toolbox in QGIS?

text: I already tried this http://hg.orfeo-toolbox.org/OTB-Qgis-plugins/ repository and couldn't get it... This guy seems to have it: http://www.elance.com/samples/orfeo-toolbox-qgis/24111664/ Help, please. :=)


title: Calculate projected area of polygons with field calculator in ArcMap 10

text: I've got a polygon featurelayer in arcmap. The coordinate system of the featurelayer's featureclass is web mercator. The dataframe of the map uses an equal-area projected coordinate system. I've got a field on the featureclass called Acreage. I'd like to calculate Acreage based on the polygon's re-projected area (into that of the dataframe's projection). I don't want to create a new featureclass to do this. ~~Is there a python script that would let me do this in the ArcMap 10 Field Calculator?~~ Update: Is there an easy way to do this in ArcMap 10?


title: ArcSDE recently modified geometries

text: ArcSDE 10 with a MS SQL 2005 backend. I am trying to see the most recently modified geometries using SQL so that I can use it to update another GIS database running PostGIS. Has anyone tried this before? I am using jTDS and oracle SQL developer to connect to the MS SQL database, if that helps.