title: Feature Access ArcGIS 10.1 Registered Database Publishing Error

text: I created a feature class from an XY table in ArcCatalog, and I am trying to publish a feature access enabled map to my server. When I try, I keep getting a 00090 error which is "feature service requires a registered database". I registered my database using server manager by creating a .sde file with a python script, and the table and shapefile are in my enterprise geodatabase (created and enabled with the tool in ArcGIS). The database validates on the server manager page. I have an ArcGIS server (10.1 SP1) running on a Microsoft 2012 Datacenter Edition Server. I am using MSSQL 2012 enterprise editionand trying to publish with ArcMap 10.1. I am very new to ArcGIS, so I may be forgetting/overlooking something silly, but any help would be greatly appreciated. I tried ArcGIS Server Error: "Feature service requires a registered database".


title: Zip + 4 Centroids

text: I have a data set with ZIP +4 codes, and ideally, I'd like to map them (their centroids) as points. I'm not interested in their boundaries (since they're not really geographies anyway). Does anyone know of a resource that has point data for ZIP + 4s? Or might there be a better way to go about mapping this type of data? Thanks in advance!


title: Using GDAL to connect to a WFS

text: I'm trying to connect to a wfs using GDAL, from reading some of the (fairly sparse) documentation on this i have found that this is possible if GDAL is compiled with the correct drivers (Curl, Xerces). I have compiled GDAL against both of these and have a working WFS Service to connect to. From reading the documentation it says this is how the connection should be opened.               GDALDataset  *poDataset;          GDALAllRegister();          poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );      I have done this and replaced pszFilename with a URL to my WFS which didnt work and gave the error (and i have no idea what file.dat is referring to as all im doing here is entering a URL that should return some XML):               ERROR 4: `/vsimem/http_1/file.dat' not recognised as a supported file format.     ERROR 4: Failed to create temporary file:/tmp/file.dat        alternatively it says that you could use an XML File that contains the URL like this:               <OGRWFSDataSource>         <URL>http://path/to/WFS/service[?OPTIONAL_PARAMETER1=VALUE[&OPTIONNAL_PARAMETER2=VALUE]]</URL>     </OGRWFSDataSource>      I have made this XML file and put in my URL for the WFS im trying to connect to and i get this error even though the documentation says that reading the url from XML is supported:               ERROR 4: `:/Resources/WFS.xml' does not exist in the file system,     and is not recognised as a supported dataset name.      Im wondering if anyone knows how to connect to a WFS using GDAL because as far as documentation goes this is about all its says in the way of how to connect and goes into no detail about supported file formats or how to connect to any of these. Any help would be appreciated (oh and btw my WFS is working because if i type the URL into my web browser i get the XML (or is it GML?) in the web browser for the request i send (a GetCapabilities one))


title: Using Regression Trees in ArcGIS

text: Hi I have a follow up question to this post: Statistical analyses for survey against raster? I have a very similar problem as the initial poster David. I am creating a multivariate model from raster data in order to predict a response variable. In my case I am using soil and climatic data to predict crop yields. I independently came to a very similar conclusion as the answer posted by Aaron, to use zonal stats, export to R, and create a Random Forest model. However, I am now stuck because I want to create a raster in ArcGIS of my predicted yields based on the input rasters but I can't figure out how to translate the output from R to Arc. In R I used this code:               library(randomForest)     fit <- randomForest(Yield ~ .,   data=data)     print(fit) # view results      importance(fit) # importance of each predictor     plot(predict(fit))      The output a plot that reasonably predicts yield. What I want to do is take the Random Forest model and use rasters of my predictor variables (e.g. slope) from a new location, and predict the yields for the new location. I understand how this works with simple linear regression but I can't wrap my head around how to do this with Random Forests. What equation does RF use to make the prediction plot? How can I get that equation and its coefficients from R to make a new map in Arc? Thanks for any help in advance.


title: How to edit (particularly flip direction of) polyline obtained through UpdateCursor?

text: I want to flip the direction of a polyline feature obtained as a result of an Update query; e.g.               lines= arcpy.UpdateCursor(myFeatureClass,myQuery)     dsc= arcpy.Describe(myFeatureClass)     try:         for ln in lines:             # ? Check number of parts             # ? Flip direction of each, or some of the parts             lines.updateRow(ln)     except:         print "Error:", sys.exc_info()[0]     finally:         # releasing locks         if lines: del lines         if ln: del ln      BTW, any better way to know how many rows are returned than counting them?


title: How many polygons under a point (ArcGIS solution preferred)

text: Opposite of Point in Poly. How many polygons underneath a point. I have a point layer with many points. I have a polygon layer with many polygons. I need to add an attribute to the points that states how many of the polygons it is in That is something like Point 1 18 polygons Point 2 10 polygons.


title: Continuous Polyline Network between two cities

text: I'm a relatively new user using ArcGIS 10.1 (ArcMap). I wish to construct all the rail lines between New York and San Francisco. I have all the rail lines as poly-lines that span the entire country, and the co-ordinates of the two cities. And I wish to only select those polylines that will give me a continuous network between NY and SFO. For example, if the line goes from NY to Detroit, and another from Detroit to Chicago and another from Chicago to SFO, then I wish to select that route. But if it goes from New Orleans to Austin, then I don't want it. I can't really use the "Select by Attributes" option since the attributes table doesn't have any information other than the ID and Shape. Thanks for your help!


title: Feature Tolerances and Loading Data

text: Keeping data consistent and spatially accurate is obviously important. My task on hand would be converting my current data into the Local Government Information Model (LGIM). However, upper management and myself have the concern that feature tolerances may be an issue when loading data from one feature class within a geodatabase to a new feature class within that file geodatabase. My workflow would look something like this:   1. Import empty LGIM Feature classes into current file geodatabase.   2. Load into corresponding empty LGIM feature classes data from current feature classes. Pretty simple? However, I am working with utilities data. A current feature class contains every point feature (manholes, points of connection, pump stations etc). LGIM separates these different types of point features. So the new workflow would look like this:   1. Separate current feature classes within the geodatabase.   2. Import empty LGIM Feature classes into current file geodatabase.   3. Load into corresponding empty LGIM feature classes data from current feature classes. How will my data be affected in terms of spatial accuracy and how does feature tolerance play a role in any of the workflow steps? If I end up continuing with this workflow, will my topology get messed up?


title: Does the NED use geocentric or geodetic latitude?

text: What kind of latitude does the National Elevation Dataset (NED) use? Geocentric or geodetic or something else? I can't find the information anywhere on their site (http://ned.usgs.gov/index.html) or in the documentation that comes with the data I downloaded. I'm new to GIS. If latitude type is unspecified, is there a standard default assumption?


title: Where should I create my network dataset in ArcGIS?

text: I want to create a new network dataset from a street layer. Should I create it in a normal file or in a personal geodatabase?


title: Field is not nullable but Field is nullable Box is checked

text: I had this problem already a couple of times. I add a Field (float) to an existing feature class, the `Field is Nullable` is arleady checked. But performing a calculation with the `Field Calculator` on it where the `value` would return `0` (zero, the number, because I divide 27 / 0), I will get the error message `Field is not Nullable`. How comes? It was checked that the `field` is Nullable. Any ideas where that comes from?


title: PostGIS Transform distance or create a view

text: I have a postgis table that is updated on a regular basis so it is not practical to cast my geometry column to geography So can I modify the following query to convert a distance say 0.5k into decimal Degrees on the fly at the latitude of the given point?               Select colA,colB,lat_dms,long_dms      FROM     my_geo_table     WHERE     ST_DWithin(SRID=4326;POINT(LONG LAT)’,geom,0.5);      Id be replacing 0.5 which in above equals 30 minutes with ST_Transform(500) but am unsure of the syntax Or would I be better of creating and saving a VIEW               CREATE VIEW geog_from_geom_view      Select colA,colB,lat_dms,long_dms      FROM     my_geo_table;          ALTER TABLE geog_from_geom_view  ADD column geog;     Update geog_from_geom_view set geog = ST__GeographyFromText(ST_as_Text(my_geo_table.geom));      I tried above but cauld not make it work I kept getting an error saying geog_from_geom_view is not a table also tried ALTER VIEW with Same error If I created a View what would the Syntax of the new query be? Thanx in advanced BTW I’m on windows 7 using pgAdmin posgis 2.0 with no CMD line access **UPDATE** To answer the questions martin f 1) > decimal DMS Should be decimal Degrees 2) > the 500 in ST_Transform() an SRID? no it is 0.5k in meters so should it be ST_Transform(POINT, 500, 4326); but that doesn't make sense when I'm trying to find all points within 0.5k of a given LAT Long **Update 2** in PgAdmin I created a view my view definition now looks like this               select colA,colB,lat_dms,long_dms ,geom     FROM     my_geo_table,     Geog GEOGRAPHY(POINT,4326),     geom::geography AS geog;      I now get an syntax error at or near “4326” LINE 4 Geog GEOGRAPHY(POINT,4326) any suggestions?


title: Plotting a certain range in azimuthal projection using ggplot2

text: I would like to limit a map to a particular range of values (i.e. northern hemisphere) under an azimuthal or orthographic projection in ggplot2. When I set ylim, I receive an error message that relates to `geom_path()` and certain portions of the map are missing (note the triangle cut out of the United States). Any ideas on how to plot a particular region using a non-cartesian projection without losing portions of the shapefile?               world<-map_data('world')           p<-ggplot(data=world,aes(x=long,y=lat,group=group)) +     ylim(30,90) +     geom_polygon(fill='skyblue') +     geom_path(color='gray60',aes(x=long,y=lat),size=0.1) +     coord_map("orthographic")     p          Warning message:     Removed 276 rows containing missing values (geom_path).       ![example](http://i.stack.imgur.com/WLAcR.png)


title: Rasters added to mosaic dataset only display at large acale

text: I have a set of 357 rasters in .jp2 that I've added to a mosaic dataset. All but 7 of them display properly, but the 7 only display at larger scales, around 1:6500 and under. I've tried adding them several times, deleting duplicates, overwriting, building overviews, etc., but can't figure out how to get the 7 to display normally like all the rest. Any clues? Thanks, much appreciated. Using ArcGIS 10.0, sp5, Windows 7 Enterprise.


title: How do I Fix Projection Issues in QGIS

text: I have some exiting shape files that supposedly have a CRS of OSGEO:41001. If I import new data, and set the CRS to OSGEO:41001 they do not align. If I change the CRS of the new layers to EPSG:4326, everything appears to line up. While I can create pretty pictures, I can't do any distance or co-location calculations between layers. Apparently, I incorrectly transformed some data once in the past. Starting for scratch, import old files (a base map of US) loads and displays. Screen coordinates for St George UT show as -12460233, 4682824. Load up new data consisting of Lat Long delimited text, the screen coordinates for St George show as -113.55, 37.02. Both layers have a CRS of EPSG:4326 Question: How can I fix my old files to be correct? Question: What did I do wrong?


title: Locating the centroid (center of mass) of spherical polygons

text: I'm trying to work out how best to locate the centroid of an arbitrary shape draped over a unit sphere, with the input being ordered (clockwise or anti-cw) vertices for the shape boundary. The density of vertices is irregular along the boundary, so the arc-lengths between them are not generally equal. Because the shapes may be very large (half a hemisphere) it is generally not possible to simply project the vertices to a plane and use planar methods, as detailed on Wikipedia. A slightly better approach involves the use of planar geometry manipulated in spherical coordinates, but again, with large polygons this method fails, as nicely illustrated here. On that same page, @Cffk highlighted a paper which describes a method for calculating the centroid of spherical triangles. I've tried to implement this method, but haven't yet succeeded in calculating centroids that seem to fall in the proper place. My attempt is implemented in Python and posted on StackOverflow here. I can copy the code over if wanted, but basically what I want to know is if the method described in the paper is applicable to a more complicated polygon boundary. Specifically by (1) stepping the base of the triangles along paired vertices (c.w. or c.c.w. ordered) of the polygon boundary, (2) calculating sub-centroids for each of these triangles and then (3) adding the area-scaled elements of the sub-centroids together to get a centroid value for the polygon as a whole. In short, I am treating the equations in the paper as a 3d extension of the commonly-used method of determining the area and centroid of a planar polygon, but something doesn't appear to be holding up, and at the moment it seems more likely a flaw in the math or my logic than in the construction of the python code.


title: Go to the URL when clicked on Leaflet marker

text:
I want to implement marker click event so when user click on marker it goes to url. Url is written inside an array (a[3]). However, this is my code which doesn't work, it always goes to the link of the last marker only. Can you help me? I need a function which will open a new window with url in an array of that specific marker.                for (var i = 0; i < markers.length; i++) {                 var a = markers[i];                 var title = a[2];                 var link = a[3];                 var image = a[4];                        var list =                  "
"                     + "
" + title + "
"                     + image                          +                "
";                           var marker = new L.Marker(new L.LatLng(a[0], a[1]), { riseOnHover: true }).bindLabel(list).addTo(map);                           marker.on("click", function() {                          window.open(link, "_blank");                          }                 );                    }


title: Modelling Zombie Outbreak

text: I know that there are many scripts (mostly in matlab) out there that address the spread of disease/zombie infection, but does anyone know of ways this has been developed/implemented in ArcGIS? Essentially, how would you deal with taking into account issues of:   1. Terrain   2. Population density (by polygon or raster) Would you write it in something like python or are there tools/models/scripts that already exist that I don't know about? From what I can find, a majority of the modelling is done in math, like here. What steps would you take personally to answer the question: Given point X, how fast would a zombie outbreak spread given Y speed, Z terrain, or other variables? The reason I ask, I teach GIS to a lot of young people and would love to use this as an opportunity to inspire them. Thanks!


title: Dxf to Shp - Coordinates wrong

text: I've converted a dxf file to a shape file using the folowing ogr2ogr command:               ogr2ogr -s_srs EPSG:31464 -where "LAYER='0'" -f "ESRI Shapefile" dgm.shp DGM1.dxf      Or               ogr2ogr -where "LAYER='0'" -f "ESRI Shapefile" dgm.shp DGM1.dxf      In both cases the coordinate values seem to be confused. In the original dxf file the x-coordinates range around 450xxxx and the y-coordinates around 536xxxx in the resulting shape file x-coordinates range around 901xxxx and the y-coordinates 1073xxxx. That means coordinate values are somehow doubled(?). Can anybody explain what I have done wrong? Thanks and regards


title: Plotting spatial data when two spatial objects have different CRS

text: I have a spatial polygon object and a spatial points object. The latter was created from xy latlong data vectors (called latitude and longitude, respectively) in a dataframe, while the former was simply read into R directly using rgdal. My code is as follows:               boros <- readOGR(dsn = ".", "nybb")     rats <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv", header = TRUE)     coordinates(rats) <- ~longitude + latitude      At this point neither spatial object is projected. If I project these objects as follows:               proj4string(boros) <- CRS("+proj=lcc")     proj4string(rats) <- CRS("+proj=lcc")      Both objects are now projected, and both will successfully map with the plot() function as follows:               plot(boros)     plot(rats)      However when I try to plot them together:               plot(boros)     plot(rats, add = TRUE)      I get the first plot only, without the rats object superimposed over boros. However, and this is the big problem, I get NO error message, so I have been unable to determine where the disconnection is between these two spatial objects being able to speak to each other. Both commands run smoothly without error or warning, yet I am left with just the single plot. And when I check the projections of each object with proj4string() I get the same projection returned for each object. I have now spent many, many hours over several days trying various ways of creating two spatial objects whose CRS and projections match such that they can be mapped together using plot(). Incidentally, one approach I took was to create a shapefile in ArcGIS for the rats object, which worked fine to create the file. But I am still left with the same inability of the two spatial objects to work together in R. I have tried many different projections for both objects, spTransform on both objects, nothing seems to work. Any help would be most appreciated. I have also included a dropbox link with the 2 data files I have described above: https://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384


title: Why do Shapefiles and GeoJSON behave differently in GDAL Python?

text: I have a `test.shp`. I loop over the file, buffer the features, and write the buffers to a new file. I tested writing the buffers to a Shapefile and to a GeoJSON. If I try to use the new buffer layer afterwards, for example to get a feature count, I get a number for the Shapefile `90 features`, but not for the GeoJSON `0 features`. What is the reason for that? The code looks like this               import ogr     test = ogr.Open('test.shp', 0)     lyrTest = test.GetLayer()      For the Shapefile               shpdriver = ogr.GetDriverByName('ESRI Shapefile')     ds = shpdriver.CreateDataSource('Buffer.shp')     lyrBuffer = ds.CreateLayer('buffer', geom_type=ogr.wkbPolygon)     featureDefn = lyrBuffer.GetLayerDefn()          featureTest = lyrTest.GetNextFeature()     while featureTest:       geomTest = featureTest.GetGeometryRef()       geomBuffer = geomTest.Buffer(250)       outFeature = ogr.Feature(featureDefn)       outFeature.SetGeometry(geomBuffer)       lyrBuffer.CreateFeature(outFeature)       outFeature.Destroy()       featureTest = lyrTest.GetNextFeature()          print lyrBuffer.GetFeatureCount()      For the GeoJSON               GeoJSONdriver = ogr.GetDriverByName('GeoJSON')     ds = GeoJSONdriver.CreateDataSource('Buffer.geojson')     lyrBuffer = ds.CreateLayer('Buffer.geojson', geom_type=ogr.wkbPolygon)     featureDefn = lyrBuffer.GetLayerDefn()          featureTest = lyrTest.GetNextFeature()     while featureTest:       geomTest = featureTest.GetGeometryRef()       geomBuffer = geomTest.Buffer(250)       outFeature = ogr.Feature(featureDefn)       outFeature.SetGeometry(geomBuffer)       lyrBuffer.CreateFeature(outFeature)       outFeature.Destroy()       featureTest = lyrTest.GetNextFeature()          print lyrBuffer.GetFeatureCount()


title: Time attribute problems with Network Analyst

text: My network dataset is in a geographic coordinate system, and when I project into UTM to calculate Miles for a length field to use with my Speed Limits to create a time attribute it does not seem to work. It ends up being that my numbers are so ridiculously small. I have a massive network dataset and it ends up saying that my total miles is like 1 when I check the statistics in the data. Network Analyst obviously has a built in length attribute that calculates miles fine and allows me to run tools just using distance, but how can I incorporate time with these measurement problems? Is there any way I can create an attribute using the network's self-created length attribute?


title: Can I make a clickabe image map in geoserver?

text: I mean when a region of a map (for example, a country) is clicked to redirect me to another URL, not only to give me some information like `fid`,`name`, etc. You know when a region of a shapefile is clicked below the map is retrieved some information and my question is : would it be possible when the same region is clicked the user to be redirected to another page ? I am talking about SHAPEFILES not IMAGES(PNG,JPEG,etc) !! Can this be done?


I have been struggling with a csv file I need to import into QGis 2.2.0 (on Windows 7). The coordinates are in this format:               LATITUDE  LONGITUDE     424605    861204     395222    861206     303051    861219     32155     861220     395541    861229     425119    861242      Latitude and Longitude are in units of degrees, minutes and seconds according to the person who provided the data. I tried different formats (placing commas, dots, ',", etc at different places) to import this dataset into QGis, checking the "DMS coordinates" box in the "Create a layer form Delimited Text File". I always had a Delimited text file error message. What is the format expected by QGis for Latitude/Longitude in DMS? Thank you


title: Loading data and symbol them to the map dynamically

text: I am using ArcGIS SDK for Android 10.2.2, and we want to add some geographic data(graphic) to the map. While our data is deployed by a common web service(no geographic web service like wms or wfs) which can be queried by a bbox. Now we want to add the queried point to the map after the map zoom,pan or something else which make the map update. But it seems that the only map event we can use is `OnZoomListener` and `OnPanListener`, but we found that we the view(no matter center or the resolution) of the map is changed by program like this:               mapView.setResolution(...); // or     mapView.centerAt(...);      These listeners will not be triggered. So I tried to make the points a kind of `Layer`, because I think the layer must be sensitive to the map update event, but I have see almost all the api about the layer, I can not found a class to extend. In fact, the `ArcgisFeatureLayer` is the most close to our requirement, it seems that it can automatically load new data according to the map view, but we do not deploy our data by arcgis feature service and we do not plan to. So I wonder if there is an alternative way to make it work? * * * Update the data model: for example, we get data from this url:               http://server/ws?bbox=30,30,40,40&xx=xx      Then we will get:               [{id:1,x:33,y:44,name:'xx'},{id:2.....}]


title: Geoserver WFS Fails to Load in OpenLayers

text: I'm working with Postgis->Geoserver->openlayers on localhost. I followed a Geoserver tutorial to import some nyc sample data into my postgres/postgis database and link it to Geoserver. The layer seems to display correctly in Geoserver using the openlayer layer preview. I then attempted to use OpenLayers to display a WFS-T from Geoserver on localhost, routing the LH:8081 (where my geoserver is located) through a proxy.cgi script to avoid Cross- Scripting. When I visit: proxy.cgi I do get the openlayers home page (without images). But when I try to display the wfs vector layer in my openlayers html/js page in either apache LH or jetty LH:8081/geoserver/www/ the POST call to LH:8081/geoserver/wfs doesn't complete. I am dumbfounded as to why this is happening. https://www.dropbox.com/s/7j3v3riymzytyq5/openlayers%20POST.png   https://www.dropbox.com/s/n5y44bg0rz2qm8c/geoserver%20getfeature.png


title: Determining average rainfall for cumulative watershed, where only sub-watersheds defined?

text: I have a polygon layer with sub-catchments (sub-watersheds) and a grid of rainfall. I have used the Zonal Stats as a table to calculate the average rainfall for each sub-catchment, but I now need to get the cumulative average rainfall for each sub-catchment as I move down through the catchment. I'm not sure how to prepare a layer with the "accumulated" sub-catchments, and then run the Zonal statistics (which would then be on overlapping polygons which I don't think it likes!)


title: Convert GEDCOM to KML

text: Does anyone know of a way to convert GEDCOM files to KML? I am building a family tree and I would like to be able to upload the infromation I have gathered to Google Earth. I already tried to use GedGenie but it didn't work. Many thanks!


title: Calculating "travelled distance" of a GPS track without including "stops"

text: I have a set of GPS tracks from a cycling event and I want to correlate their distance along the route. Problem is: when a rider stops, generally by entering a covered place, GPS keeps logging (low precision) points around its location, and by summing those noisy segments I get a false increase in travelled distance. If I look at the plotted tracks, it's easy to "see" that actually the rider didn't move, but how could I create some cumulative distance algorithm that "doesn't count" some points because the rider is "obviously not moving" there? One interesting thing I noticed: if I consider only one coordinate, for example, only latitude for a road heading straight north, the fuzzy points at a stop keep wandering back and forth around a definite latitude, the same applying for longitude (although in this case that would be the same longitude of every point of that road segment, because it is "vertical" on the map). But if I use euclidean distance, then this random movement is converted always into "forward" movement and keeps increasing the distance. If the stop were in the middle of a straight road, I could get the average direction and "normalize" and average the positions along the road direction, but oftentimes the stop point is at a sharp turn (a cross-road) or it is the turning back point of a route. Is there any known method for this, or does anyone have a suggestion?


title: OpenLayers 3 mousePosition projection on the fly

text: I'm displaying coordinates in epsg:3857 using `ol.control.MousePosition()` with OSM. I'd like to convert these on the fly in another epsg:2154. I call in the header :               <script> src="http://cdnjs.cloudflare.com/ajax/libs/proj4js/2.2.2/proj4-src.js" type="text/javascript"></script>      Then my html :               var map = new OpenLayers.Map({             div: "map",             projection: "EPSG:3857",             controls:[                     new OpenLayers.Control.LayerSwitcher({'ascending':false}),                     new OpenLayers.Control.ScaleLine({                         'geodesic': "true"                     }),                      new OpenLayers.Control.Navigation(),                      new OpenLayers.Control.KeyboardDefaults(),                      new OpenLayers.Control.MousePosition(                          {displayProjection: new OpenLayers.Projection("EPSG:2154")}                      )                     ]             });      but it doesn't seem to work : coordinates are false, even negative where they shouldn't. Any idea ? thx


title: How to exactly georeference in QGIS

text: When I georeference a raster image in QGIS, it's always a bit displaced as in the attached image (the cyan polygon is the template, the red border is the georeferenced image), although I try to set the georeference points as exact as possible. Is there a trick to set the georeference points more exact or maybe move the raster layer a bit? Maybe there is a different software, you can georeference by moving and then import it to QGIS?![enter image description here](http://i.stack.imgur.com/kllR1.png)


title: OpenLayers 2.13.1, how could I make a click, get coordinates, send it to geoserver function?

text: I am trying to make this very simple function, at the moment I am quite stuck when I have to place the markers as start/end points. I would like to add a marker as start and another as end and then do not allow other markers. I don't know if start and end point have to stay on different layers or the same and I don't know how to make the code reckon which is the first or second click. Here is what I have done so far.                function add_start (e){                             var position = map.getLonLatFromPixel(e.xy);                             var size = new OpenLayers.Size(21,25);                             var offset = new OpenLayers.Pixel (-(size.w/2), -size.h);                             var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker-green.png',size,offset);                             //var icon_end = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png',size,offset);                             var two_points = new OpenLayers.Geometry.Collection(Array(icon));                              var start_point = two_points[0];                             //var end_point = two_points[1];                             var markerslayer = map.getLayer('Markers');                             markerslayer.addMarker(new OpenLayers.Marker(position,icon));                             var start_point = markerslayer.markers[0];                             //var end_point = markerslayer.markers[1];                                  //arkerslayer.addMarker(new OpenLayers.Marker(position,end_point));                             start_coords = start_point.lonlat;                             alert("You clicked near " + start_coords.lat + " N, " +                                                   + start_coords.lon + " E");                                              //var two_points = new OpenLayers.Geometry.Collection({});                              //var start_point = two_points[0];                             //var end_point = two_points[1];                             }


title: How to detect openlayers feature vector intersections in any zoom level?

text: I have to change vector feature when vectors are intersect each other, I wrote the following code:               // check all vectors features (arrows) in the current layer for intersections       allArrowFeatures = pointsLayer.features;       var hideArr = new Array();       var l = 0;       for(var i = 0 ;i < allArrowFeatures.length;i++){              arrowObj = allArrowFeatures[i];         for(var j = 1 ;j < allArrowFeatures.length;j++){           if(j===i){ // same vector             continue;           }           if(allArrowFeatures[j].geometry.bounds.intersectsBounds(arrowObj.geometry.bounds)){             hideArr[l] = pointsLayer.features[i].style = {display:'none'};             l++;           }              }       }      When you zooming in to the max zoom you wont see any intersections so I assume this code is OK but the problem is when I zooming out the vectors still intersects each other. ![zoom out example](http://i.stack.imgur.com/Lkk7h.png) ![zoom in example](http://i.stack.imgur.com/HDSBe.png) Its look like I don't use the correct function for my case, and I tried several ways this was the closest to my case. Any idea what openlayers method I have to use to detect if there was intersection also in zoom out probably not base on the bounds?


title: Maximize clusters of points from coordinate data

text: So I am definitely new to the GIS world, so bear with me if I don't know the best way to articulate my problem. That said, I am ready to learn, I would definitely welcome any direction that you all could point me in. **The problem:** I have a set of coordinate points (just three attributes: name, lat, long). I hope to identify the minimum number of regions (circles with a predefined radius) that would be required to encompass XX percent of the sites (Where XX is configurable). **Initial thoughts:** My thoughts were to first create clusters (of set size) of all data points, then sort the clusters by the number of points they enclose. That way I can just count the number of clusters needed (in order of priority) until they represent the desired XX percent of the points. Thus far I have tried to do some research, but I haven't found any resources where the clustering is done within a set radius. (I have found plenty where the number of clusters is set, and the size is variable - I need the reverse of that.) I have tried using QGIS, but I am open to any suggestions. I am new to GIS, but I have a decent knowledge of programming if that could be of any help.


title: ArcGIS Javascript API Query & Add features without map using REST API

text:
I need to query arcGIS map service and add this result to another feature service using ArcGIS JavaScript API. I should not use MAP object. I did query data without map from ArcGIS sample: Now i need to post this data to feature service:               <!DOCTYPE html>                   <meta http-equiv="Content-Type" content="text/html; charset=utf-8">                  <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">         <title>Query State Info without Map</title>                                               US state name :

                Any Idea?


title: Polygon Shared Boundaries in GML

text: Last week I was told on a phone conference that GML supports shared boundaries in Polygon's. I have been unable to find any documenatation to substantiate the claim. Here is the situation: I have a series of Polygons which represent areas in a forest. As a requirement for our operating license we have to submit these administrative boundaries to the government in a custom xml schema with the spatial data included as GML objects. I have a shape file which contains our polygons which was exported from our SDE. I then ran ogr2ogr on this shapefile to convert it into a GML file. I have a custom program that reads that GML file and then inserts the other information that the government requires. The problem is the file is too large for the government site to accept, the original shapefile is 8.5 MB and the GML file is roughly 12.5 MB. Generalizing the shapefile isn't an option unfortunately because that alters the boundaries of the polygons and could lead to compliance issues. Can GML Polygon, or MultiPolygons share boundaries? If so can someone point me to some documentation to that effect.


title: Iterating field values using ModelBuilder?

text: I'm creating a model that takes field values from points and groups them together, creates a middle point for each group as well as creates a middle point between each group and then creates a spanning tree to all the points. Right now I'm stuck on only being able to perform this model on two different groups and would like to ultimately run this model until all the groups in the input feature field are used. I would like to know if its possible to do a field value iteration process that can complete this task or if anyone would have any suggestions on how I can create my final output. ![enter image description here](http://i.stack.imgur.com/5xjxB.jpg) ![enter image description here](http://i.stack.imgur.com/L9U9F.jpg)


title: Geoserver LegendGraphic support for Rules in SLDs

text: I want to serv custom legend via GetLegendGraphic WMS request. I find two answers:   1. http://jira.codehaus.org/browse/GEOS-3747?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel   2. https://github.com/geoserver/geoserver/pull/397 On first site there is patch file to geoserver, but I don't know how to use it (where to paste it). On second site is github so I clone repo and build geoserver by myself, but what now? I try add `<LegendGraphic> -> <Graphic> -> <ExternalGraphic> -> <OnlineResource>` to SLD file but geoserver don't validate it. So can anyone explane how to get custom legend?


title: Merge (a lot of) features using QGIS or other open source

text: I have one shape-file with a lot of features (polygons) and if I try to merge these into one big polygon in QGIS it takes days to complete the merge. So I wonder if there is another way? I saw this question: How to "Merge Selected Features" with python? The problem is that I don't know python and I have no idea of what they are talking about... So, could someone either give me some advice to perhaps do this operation from command line (I'm using Mac) or maybe if someone has a ready python script that does the trick?


title: Local Map Web API Without internet connection

text: I have a single-page HTML5 app that consumes a bunch of different map providers. One of my requirements is that this work with no internet connection and some sort of map provider running on the box itself. What , if any, providers -- free or commercial can give me something along the lines of an offline google maps or google earth that I can use in a web application? Thanks!


title: Register and display real time tracking of my smartphone (gps) on a map

text: I have 2 smartphones (android galaxy and iphone) and i would like to display the position, according to the gps coordinates, on a map. Is there any app to download or any kind of service where i could register? Every time I start the app, the coordinates of the smartphone would be sent to some web service and few seconds later that actual position would be displayed on a map. The map should be accessed/viewed by a user from a desktop pc through a normal Internet Browser. edit: any known open source projects where you could setup your own "tracking system" and expand/modify it?


title: How to find the nearest line from a point?

text: I have a table of Roads in PostGIS which geometry is " **geometry(MultiLineString,4326)** ". Its attribute are > gid osm_id name ref type oneway bridge maxspeed geom Now i want to know the closest Road from latitude and longitude. Also I want to set some distance from that **Line/Road** like **20m** both side and if the point is within that distance then it should give the name of that line else not. Thank you


title: How to correlate raster map and elevation data?

text: I have been provided with a set of raster maps. Thanks to another helpful user here I have been able to acquire some elevation data for the same region from gis.utah.gov in asc format. Now I need to correlate the two data sets. Can this be done in Matlab? If not, I am a proficient programmer and would appreciate instructions on how to do this in any language.


title: How to handle 500 Internal Server Error with ESRI Silverlight API

text: I have a simple AGS map service. I accidentally specified an invalid SRID explicitly in my Silverlight XAML:               <esri:Map x:Name="Map" IsLogoVisible="False">         <esri:Map.Extent>             <esri:Envelope XMin="526801.81" YMin="489775.0655" XMax="593105.41" YMax="523904.2045">                 <esri:Envelope.SpatialReference>                     <esri:SpatialReference WKID="3807" />                 </esri:Envelope.SpatialReference>             </esri:Envelope>         </esri:Map.Extent>         <esri:Map.Layers>             <esri:ArcGISDynamicMapServiceLayer x:Name="MyLayer" ID="MyService" Url="http://localhost/ArcGIS/rest/services/MyService/MapServer" Initialized="MyLayer_Initialized" InitializationFailed="Layer_InitializationFailed" />         </esri:Map.Layers>     </esri:Map>      This results in a RESTful request that looks like this:               http://localhost/ArcGIS/rest/services/MyService/MapServer/export?bbox=526801.8125,479254.134305427,593105.4375,534425.146944573&size=1078,897&format=png24&transparent=true&layers=show:0,1&f=image&imageSR=3807&bboxSR=3807      This request fails with the following IIS error:               Error Summary          HTTP Error 500.0 - Internal Server Error     The page cannot be displayed because an internal server error has occurred.          Detailed Error Information          Module  ESRI.RestModule-Integrated     Notification    ResolveRequestCache     Handler StaticFile     Error Code  0x00000000      If I replace the invalid SRID, everything works fine. Is there any way to get the Silverlight API to report exceptions like this? As it is I have to use Fiddler to find it. ![enter image description here](http://i.stack.imgur.com/cu65k.png)


title: protecting kml access to google only

text: I have a large amount of data (thousands of polygons) that I want to share with people through my site. I want them to be able to view the polygons, and due to the large amount of data, I am not worried about them stealing this data through the displayed google map. However, since the kml file must be accessible to google, a proficient coder could discover the url of the kml files and simply download the entire file. Is there any way to restrict access to these kml files only to google so no outside people can get to them. I have tried .htaccess, but I can't figure out how to restrict access to everyone but google. Any ideas?


title: Best routing algorithm for this scenario

text: Trying to work out the most appropriate routing algorithm to use for this situation. Even better if it's available as part of pgrouting or similar.. Here's the requirements:   * Specified date/time for start point   * Specified date/time for end point   * Constant speed of travelling   * Possibility for a few additional minutes added per hour of travelling   * Take into account impedences/weights of certain routes   * The number of stages / legs doesn't matter   * Preferable not to retraverse stages / legs but can be done if required   * Start / end point could be the same   * The algorithm should return the top 3 matching routes (if possible) The idea is that the routing algorithm would work out a near optimal route to take that will take from the start time until the end time. So for example, if the start time is 9am and the finish time is 12pm then a 3hr3min route would take preference to a 2hr route. Any thoughts?


title: ArcGIS javascript. Trying to use a dojo table of contents in CakePHP

text: Sorry for asking so many questions, but I am trying to learn arcgis, dojo, and cake all at the same time :/ Basically what I want to do is get this to work in the latest CakePHP. Where do I place the files and how do I call it?


title: ArcObjects - Export map to PDF - graphics not appearing

text: I'm having an issue with our ArcObjects code for exporting a map to PDF. There is a section of code that adds graphics on to the map (polygons, lines, text etc). This seems to work fine for the first map generated after the map service is started, but on subsequent prints the markup does not get displayed on the map. This code worked previously but no we've converted it to run as an SOE on ArcGIS Server 10.0. The code in question:                       // IMap map, IGraphicElements graphicElements                  GraphicUtil.GetGraphicsContainer(map).DeleteAllElements(); // clear any previous markup                  _logger.Debug("Add graphic elements to container");             IGraphicsContainer graphicsContainer = GetGraphicsContainer(map);             graphicsContainer.Reset();             for (int i = 0; i < graphicElements.Count; i++)             {                 IElement element = graphicElements.get_Element(i) as IElement;                 _logger.Debug("Adding geometry type: " + element.Geometry.GeometryType.ToString());                 graphicsContainer.AddElement(element, 0);             }                  ((IActiveView)map).PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);              ...              public static IGraphicsContainer GetGraphicsContainer(IMap map)         {             return (IGraphicsContainer)map;         }      If I take out the DeleteAllElements the markup stays on the map between requests but any new markup doesn't get added. If I restart the map service, then it will add the markup on the next request, then not add anything after that. In the ArcGIS Server logs I can see the following items which indicates nothing is being drawn, even though the graphic elements have been correctly added to the graphics container in the code (as far as I can tell)               INFO3   End of the graphics drawing.     INFO3   Graphic count: 0.     INFO3   Beginning of the graphics drawing.      Any suggestions would be much appreciated!


title: Where can I download ArcGIS Engine 10?

text: Can anyone help me find the link to download ArcGIS Engine 10? I have tried to find it, but unfortunately I have not found it.


title: List of available online WMS services (Weather, Land Data, Place Names...)

text: I'm Looking for WMS services freely available to the public to overlay on my application, specifically weather and AIS information. Perhaps this is a good place to list some available online WMS services you have stumbled across?


title: Get points along a line at defined distance

text: I was reading a lot of similar questions here but no one ask about postgis (or I didn't find it). I have a line in postgis and I want to get the point along the line at defined distance. It is possible? I found st_locate_point but it doesn't work for I want to do. Any idea?


title: arcpy.mapping.ListLayoutElements() takes longer in a tool than in ArcMap

text: Has anyone else noticed that accessing elements from arcpy.mapping.ListLayoutElements() through a tool script takes longer than within ArcMap? Whether using a specific name, wildcard or nothing, the call takes on the order of 90 seconds or longer to return:               arcpy.AddMessage("  Start: " + str(datetime.now()))     mapDoc = arcpy.mapping.MapDocument(mapDocumentPath)     elements = arcpy.mapping.ListLayoutElements(mapDoc, "GRAPHIC_ELEMENT")     arcpy.AddMessage("  Finished: " + str(datetime.now()))          Start: 2011-05-24 12:29:28.479000     Finished: 2011-05-24 12:31:10.510000      When the call is made within ArcMap the return is instantaneous.               mapDoc = arcpy.mapping.MapDocument("CURRENT")     elements = arcpy.mapping.ListLayoutElements(mapDoc, "GRAPHIC_ELEMENT")      In both instances, I'm referencing the same map document. Using ArcGIS 10 SP2. **Update:**   1. It may not be obvious but the results are instantaneous when the script is run against "CURRENT". Specify any map document other than "CURRENT" and the script will slow down.   2. I executed a loop as Jason suggested and found little difference between the first call to ListLayoutElements() and the last. The first selection took 28 seconds while the other nine averaged 24 seconds.   3. The map document has nine data frames and nine corresponding legend elements that were converted to graphics. When these "legends" are present in the document then ListLayoutElements() is slow; when they're deleted then ListLayoutElements() returns immediately.


title: How to insert 3d line into postgis database?

text: I have 3D line which I exported from GRASS in .shp format. When I try to insert this line to Postgis database with QuantumGIS, it gives an error: > ERROR: Invalid type name - valid ones are: POINT, MULTIPOINT, LINESTRING, > MULTILINESTRING, POLYGON, MULTIPOLYGON, CIRCULARSTRING, COMPOUNDCURVE, > MULTICURVE, CURVEPOLYGON, MULTISURFACE, GEOMETRY, GEOMETRYCOLLECTION, > POINTM, MULTIPOINTM, LINESTRINGM, MULTILINESTRINGM, POLYGONM, MULTIPOLYGONM, > CIRCULARSTRINGM, COMPOUNDCURVEM, MULTICURVEM CURVEPOLYGONM, MULTISURFACEM, > or GEOMETRYCOLLECTIONM CONTEXT: SQL statement "SELECT AddGeometryColumn('', > $1 , $2 , $3 , $4 , $5 , $6 )" PL/pgSQL function "addgeometrycolumn" line 4 > at SQL statement Can someone help me? How can I insert this 3d line into postgis database?


title: How to import qgis.core in a Osgeo4W installation?

text: when I try to import qgis core module from python:   `from qgis.core import *`   I get this error:   `ImportError: DLL load failed: The operating system cannot run %1.` Searching the web I found that I need to point my PATH to the qgis installation, so I did that:   `set PATH=C:\OSGeo4Windows\apps\qgis-dev\;%PATH%` But I still get the error. obs, i'm using: Osgeo4w Qgis 1.7, Windows 7, Aptana Studio 3 with Pydev.


title: Definition Query to Python script on Oracle 9i and ArcGIS 9.3.1

text: SUBTYPECD=1 AND SERIESNAME NOT LIKE '%-%' AND OPERATIONALSTATUS IS NULL OR SUBTYPECD=1 AND SERIESNAME NOT LIKE '%-%' AND OPERATIONALSTATUS = 976 The above is my def query in ArcMap against ArcSDE 9.3.1 on oracle 9i with python 2.5. I want to be using the gp.makefeaturelayer with the sql parameter filled in with a rendition of the above. ps I started to ask this question because it seemed daunting and I am in a crunch..then figured it out so I thought I'd finish the post for those to come.


title: Access is Denied error on ProjectEx call

text: I need to be able to convert a lat/long coordinate into a point coordinate and am having trouble getting my code to work. I created a test app isloating the code that is giving me problems. I'm running ArcSDE 9.3.1 SP2 on my local machine and also have it running on another box on the network, Jester. If I run my test code with passing in "localhost" to be the host, it works fine. If I pass in "Jester", I get the error: System.UnauthorizedAccessException: Access is Denied (Exception from HRESULT: 0x80070005 (E_ACCESSDENIED)) and it occurs on the line where I do geom.ProjectEx. The S4_SPATIAL geodataService is created from the same database on both boxes. In the code this snippet came from, I use the same server context created from the same geodata service on Jester and it works fine ( I can get to all my tables, feature classes and locators, I can query them, etc.), except when I try to do geom.ProjectEx. I've looked online for this error and haven't been able to find anything helpful. Apparently its a pretty generic message. Does anyone have any idea why it would work on my local server but not on the other? What am I missing?               static void Main(string[] args)     {              double latitude = 25.73338;         double longitude = -80.3207233333334;              IGISServerConnection conn = new GISServerConnection();              string host = args[0];;              conn.Connect(host);         IServerObjectManager som = conn.ServerObjectManager;         IServerContext gdbContext = som.CreateServerContext("S4_Spatial", "GeodataServer");                   IPoint pt = (IPoint)gdbContext.CreateObject("esriGeometry.Point");         pt.Y = latitude;         pt.X = longitude;              // Create spatial reference for the latitude and longitude point         SpatialReferenceEnvironment spatialReferenceEnv = (SpatialReferenceEnvironment)gdbContext.CreateObject("esriGeometry.SpatialReferenceEnvironment");         ISpatialReference latlongRef = spatialReferenceEnv.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);         pt.SpatialReference = latlongRef;              // Get spatial reference of the map          ISpatialReference outgoingCoordSystem = spatialReferenceEnv.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_NAD1983HARN);              NADCONTransformation datumConversion = new NADCONTransformation(spatialReferenceEnv.CreateGeoTransformation((int)esriSRGeoTransformation2Type.esriSRGeoTransformation_NAD_1983_TO_HARN_FL));              IGeometry2 geom = (IGeometry2)pt;         geom.SpatialReference = latlongRef;              // Project the latitude and longitude to map coordinate system         geom.ProjectEx(outgoingCoordSystem, esriTransformDirection.esriTransformForward, (IGeoTransformation)datumConversion, false, 0.0, 0.0);              IPoint newPt = (IPoint)geom;         Console.WriteLine(newPt.X + " , " + newPt.Y);     }


title: Batch Georeferencing in ArcMap 10.2

text: I have about 50 pdf maps that I would like to batch georeference and/or georectify. I've saved one as a jpg and georeferenced it manually and would like to automate the rest which are in the same location and have the same extent. ![enter image description here](http://i.stack.imgur.com/dOhmU.jpg) The density gradient will change in each case, but the map (saved as jpg now) will have the same extent for all 50 cases. Does anyone know of a link to a **step-by-step** tutorial that deals with batch processing in Arcmap 10.2 only? I'm thinking of using the Warp from File code or possibly the Autoregistration option, but am not sure where to start? The sample code is:               import arcpy     arcpy.WarpFromFile_management(     "\\cpu\data\raster.img", "\\cpu\data\warp_out.tif",     "\\cpu\data\gcpfile.txt", "POLYORDER2", "BILINEAR")


title: How to disable overlay check box after a certain zoom level?

text: How do I disable a particular overlay checkbox in the Openlayers Layerswitcher?


title: How to compile Degrib on Ubuntu 12.04 LTS?

text: I'm trying to compile degrib on a typical Ubuntu 12.04 LTS Server. However, I get an error "undefined reference to gzopen64" which appears to be common when compiling degrib. Do you have any suggestions? I see Salda got it working at How to convert GRIB-1 to GRIB-2?. Any suggestions will be much appreciated.


title: Is adding a new table for every layer a good design idea?

text: Is it a good design idea to maintain separate tables for each layer in a PostGIS database?. Every layer is defined by users. A user creates his/her layer say ATM location. Thats meant by defined by users(I dont have control over the table definition of the layers much) This makes it impossible to predict fields(attribute fields and type are not known initially) used in the table/layers. School vs roads have different geometry and set of attributes. The problem is traditional databases wont advice for having dynamic tables very much. Does the same apply to spatial databases?


title: SLD filter with embedded strTrim function in Geoserver

text: Using Geoserver 2.2.2 (GeoTools Version 8.4) on Windows 7, I'm trying to create a SLD that trims the leading and trailing spaces of a string property (using the strTrim function) and compares it to a given string pattern (using PropertyIsLike) for filtering. Here is the code:               <?xml version="1.0" encoding="ISO-8859-1"?>     <StyledLayerDescriptor version="1.0.0"      xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"      xmlns="http://www.opengis.net/sld"      xmlns:ogc="http://www.opengis.net/ogc"      xmlns:xlink="http://www.w3.org/1999/xlink"      xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">            <NamedLayer>         <Name>style1</Name>         <UserStyle>           <Title>style1</Title>           <Abstract>style1</Abstract>           <FeatureTypeStyle>             <Rule>               <Name>Fauche</Name>               <Title>Fauche</Title>               <ogc:Filter>                 <ogc:PropertyIsLike wildCard="*" singleChar="#" escape="!">                   <ogc:Function name="strTrim">                   <ogc:PropertyName>type</ogc:PropertyName>                   </ogc:Function>                   <ogc:Literal>Testword</ogc:Literal>                 </ogc:PropertyIsLike>               </ogc:Filter>               <PolygonSymbolizer>                       <Fill>                   <CssParameter name="fill">#FF8C00</CssParameter>                   <CssParameter name="fill-opacity">0.6</CssParameter>                 </Fill>               </PolygonSymbolizer>             </Rule>            </FeatureTypeStyle>        </UserStyle>      When validating the SLD, I get the following error: > cvc-complex-type.2.4.a: Invalid content was found starting with element > 'ogc:Function'. One of '{"http://www.opengis.net/ogc":PropertyName}' is > expected. Where is the mistake?


title: openstreetmaps offline tile packs

text: Is there a place where I could download openstreetmaps tile packs for the united states? I would like to have offline access for the whole country. Thanks,


title: How to Display the World Map in Projected Coordinate Systems

text: I am trying to create a simple presentation map and I download the global level country boundaries shapefile from http://www.diva-gis.org/Data. The layer coordinate system is in Geographic Coordinate System: GCS_WGS_1984. Can you please let me know how I can display the map on UTM or any projected coordinate System using ArcMap? Thanks


title: OpenStreetMap Tiles Not Loading in Firefox 17.0.1 - Pink Background

text: I recently ran into an issue where only a pink background is displayed in the openstreetmap on my webpage (tiles are not displaying). I found this topic OpenStreetMap tiles not loaded despite being received which would seem to address the issue but the solution presented no longer fixes it. I tried adding the suggested code to my layer but the issue persists.               map.addLayer(new OpenLayers.Layer.GeoRSS("My Points", "georss", {         tileOptions: {             crossOriginKeyword: null         }      }));      Any ideas on how to fix this in the newer version of Firefox?


title: How to get max value from 3 vector files?

text: I'd like to know how to the get max value from 3 vector files. Originally, the vector file are interpolated images. (First, point data are interpolated into surface data by IDW method and then exported into vector file). Now I need a vector file with max value obtained from the 3 vector files. When I work on raster image to get max value, firstly, I always do layer stack operation on all images from which max value will extract. After that, I get max value image from stacked image. So what I want to ask is do I need to do layer stack on 3 vector images? If so, what is the operation name/tool in ArcGIS and how can I find max value from stacked vector files?


title: How to convert overlapping polygons in single layer to raster by count in QGIS?

text: I'm trying to do this in QGIS 2.0 (Mac Mt Lion). But cannot do step 2, as I can't dissolve _original_ polygons by 'count' in QGIS (in that tool at least). Advice?


How can I use osm2po to calculate distance matrix in Km? There is an example that calculates distance matrix in time but I cannot make it turn into kilometers. I dont understand how should I use the int[] path parameter into the calcPathLength function.


title: How to create a raster with rotated cell orientation from lattice points - ArcMap 10.1

text: I am trying to create a raster layer from gridded points where the cells have a rotated orientation. Each point is the center of a 1.6 km by 1.6 km cell, rotated to the left. Is there a way to do this in ArcMap? I know I do not want Rotate(data management) as the cells themselves are rotated, not the entire raster. I attached a screenshot of the point layer I'm working with. ![Points example](http://i.stack.imgur.com/QPFdb.png)


title: How to specify units when using SelectLayerbyLocation_management WITHIN_A_DISTANCE?

text: Can anyone help by suggesting how to specify the units when using the SelectLayerbyLocation_management tool. I am selecting by "WITHIN_A_DISTANCE" and the ArcGIS resources do not suggest how to do this. I am getting some funny results at the moment.


title: What are the differences between WMS, WFS, WCS, WPS?

text: By looking at ![enter image description here](http://i.stack.imgur.com/IO2oR.png) I do not really understand the difference of WMS, WFS, WCS, WPS. May you please explain the differences between them?


title: Using OpenLayers with OSM as map

text: Got this code from a teacher and I don't know what I am doing wrong! I have searched around but I can't find the error. Visual Basic shows no error.               <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"      "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">     <html xmlns="http://www.w3.org/1999/xhtml">     <head>     <title>OpenLayers map</title>     <script src="http://www.openlayers.org/api/OpenLayers.js" type="text/javascript"></script>     <style type ="text/css">         /*Stilsättning*/         mapContainer         {             position:absolute;          /*Elementets position i förhållande till andra element*/             padding:10px;               /*Elementets avstånd från andra element (10px åt alla håll)*/             width:600px;                /*Elementets bredd */             height:500px;               /*Elementets höjd*/             border:1px solid black;     /*Ram runt elementet */         }     </style>     <script type="text/javascript">         //Variabel som kommer hantera OpenLayers.Map objektet.         var map;              //Funktion som anropas från <body> taggen, kör koden som finns i funktionen när sidan laddas (onload).         function initMap() {             //Initiera OpenLayers.Map objektet och tilldela till variabeln map.             //"mapContainer" är en referens till den div (container) där du vill ha din karta.             map = new OpenLayers.Map('mapContainer');                  //Variabeln osm tilldelas ett objekt med en OpenLayers-Lagertyp OSM(openstreetmap).             //Ex på olika lagertyper:  OpenLayers.Layer.WMS  (vilken WMS-tjänst som helst kan kopplas in)             //                         OpenLayers.Layer.Vector  (ett vektor lager där du kan lägga in tex punktdata)             //                         OpenLayers.Layer.Google  (om du vill använda Googles bakgrundskartor)             var osm = new OpenLayers.Layer.OSM();                  //Variabeln osm som innehåller Openstreetmap objektet läggs till kartobjektet genom map.addLayer()             map.addLayer(osm);                  //När kartan initieras måste OpenLayers veta var du vill börja någonstanns, i det här fallet zoomas kartan till              //full utzoomning i förhållande till Openstreetmaps utbredning (extent).             //Ett alternativ är att använda map.setCenter() då du anger exakta koordinater och zoomnivå men mer om det senare.               map.zoomToMaxExtent();         }     </script>


title: Grid Analysis from samples data

text: I'm trying to analyze set of samples into vector grids in order to generate average data in corresponding grid.I'm working with 50m separated grids but the size of grid may change over time.That's why we need to generate grids on the fly.There is a snapshot of data distribution below. ![vector data](http://i.stack.imgur.com/psgDk.png) There are two important things for us,one of which creating grid,as you may know we can create grid from minimum bounding rectangle (MBR) of our data.The trick is let say we have 1000 rows of distributed samples (points) for a large area,then the processing of this MBR plus dividing into grids will be too long.Since,what we want is intersecting grids alone,we must generate them from samples which means **we must** count as the length of our sample data.That's what I'm asking,is there an easy way to do it without processing all grids? The other is bound to first where our samples contain various data that we should aggregate into the single grid.How can we do it with or without processing grid at the same time ? NOTE: This question may require programming or querying which we are wiling to migrate to solution in any ways.Our samples are in text format that we can transform them to any other solutions you can suggest. Best Regards


title: how to align vector and raster images?

text: How to align vector and raster images? For example vector images will be static/rigid and raster image will be morphable. For example I manually do initial guess and algorithm must converge to optimal solution by morphing raster image. I have raster and vector data represented as point clouds in 2d. Raster data extracted as contours so it contain some noise. Raster and vector data slightly differs in scale and position,raster data can also have some small perspective distortion and barrel distortion, but I can provide initian alignment. I want to raster data "converged" to unchangeble vector data, so I think I want some kind of non-rigid point cloud matching. By alignment of 2 point clouds then I want to determine some parameters(maybe TPS spline?) and unwarp raster image and then draw vector image above it. this video using ceres solver shows something like I want to achive other maybe related imgs 1 2 3 I tried to use KCReg it use some specific metric for point clouds, but it give me unsatisfactory results because my distortion can be only roughly expressed by perspective transform and I need some more non-rigid alignment.I also tried gmmreg. I seems it uses some Broyden Fletcher Goldfarb Shannon minimization from VNL, but I managed to achive good results only for simple tasks for example here some my attempts using gmmreg for example affine http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/%D0%B0%D1%84%D0%B8%D0%BD%D0%BD%D0%BE%D0%B5.png http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/%D1%82%D1%80%D0%B0%D0%BF%D0%B5%D1%86%D0%B8%D1%8F.png barrel-barrel backwards-barrel rigid http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/barrel_test.png http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/barrel_test_back.png http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/barrel_test_rigid.png rigid scale rotate shift http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/rigid%3Dscale_rotate_shift_noise.png but it fails on real problem, maybe I don't set constraints/bounds or initial guess well. rigid shift http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/rigid_no_guess.png TPS_L2 no guess http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/TPS_L2_no_guess.png TPS_L2 initial guess http://dl.dropbox.com/u/8841028/ICP/gmmreg_tests/TPS_L2_initial_guess.png


title: Feature Class Extents Changed After Import

text: I have been tasked with editing the attributes and fixing geometry of a feature class in a shared geodatabase. The feature class is not in a feature dataset. I decided to export the feature class into a separate geodatabase with a feature dataset in order to set up topology rules on the data and clean it up. Once I finished cleaning the data, I deleted all rows in the original dataset and imported my cleaned up data. Everything worked perfectly. However, now there is something strange happening with the data. When I "zoom to layer" on my copy of the data, it does what is expected and shows the full extent of the data: ![Correct Zoom](http://i.stack.imgur.com/y4cf3.jpg) However, when I zoom to layer of the corrected data after I imported it back into the shared geodatabase, it shows an extent much larger than the actual data: ![Incorrect Zoom](http://i.stack.imgur.com/ZfmSy.jpg) Each feature class contains identical data. But when I compare the extents of each feature class, they are different. The extents of the data in my copy: ![Correct Extents](http://i.stack.imgur.com/qjb9D.jpg) The extents of the data after it was imported back into the shared geodatabase: ![Incorrect Extents](http://i.stack.imgur.com/j0T2v.jpg) Is there a way to reset the extent of the feature class? I'm reminded of a similar issue that has occurred in MapInfo that was fixed by packing the table's geographic data. Thanks.


title: How to upload TIGER shapefiles into MongoDB?

text: I have the school district shapefiles from TIGER for Michigan, from: ftp://ftp2.census.gov/geo/tiger/TIGER2012/UNSD/ I'd like to use MongoDBs geoquery features to interrogate this data. How do I upload the .sh, .sh.xml, .dbf, .prj files into MongoDB so that I can execute these queries?


title: How to find number of fields in shapefile using Python?

text: I've been looking at the ArcGIS resources for python but I can't seem to figure this one out.   I'm looking for a method to gain access to the number of fields that are within a shapefile / feature class. Is this information available via Python?


title: Can QGIS use a non-spatial Oracle Database?

text: I'm building my first plugin in qgis, that needs to read some coordinates from an oracle database (not oracle spatial database) and create a shapefile with this points. I'm searching but only find to in oracle spatial! Is there a way to do this in a simple oracle database? Thanks


title: Creating a BATCH with multiple INPUTS (variables) to run a workbench that intersects LINES AND POLYGONS

text: I'm trying to create a batch file to run my workbench.fmw in which I have two input files: 8 pipeline shapefiles and 96 communal areas shapefiles. The workbench allows cutting pipeline relative to the areas by using "Spatial Relator" Transformer. I need to do batch processing in order to run this workbench several times with different inputs automatically (by double clicking on the .bat). I tried to create my DOS script in a very complicated way so I come to you for help and suggestions and if anyone has ever created a similar batch?


title: How to save colour (color) ramps in QGIS

text: This strikes me that it should be simple, but efforts in QGIS and research here and elsewhere are not providing a solution. I have a (vector) map which I have to mark out in colours according to the values of the individual sectors. To aid the maps use for those with colour blindness issues I am only using one colour - blue. The colours in the standard colour ramp are too close so I have made my own thus: - New Colour Ramp Gradient - OK Colour 1 - white Colour 2 - blue Name the Colour Ramp - Test - OK I have now got 5 colours in front of me. Only the White is acceptable, so I double click on each colour in turn and change their values to make them more widely spaced - OK - OK. When I print out the map all is well, but when I save that project, then open another and try to apply the colour ramp the colours have gone back to the default values in that colour ramp. I have messed about in the style manager with no success. Will I have to manually set all of the values for this colour ramp in each project that I wish to apply this colour ramp to? Surely not! There must be a simpler way.


title: "Item Description" from Catalog window in ArcMap won't open

text: I can't access metadata information from the catalog window in ArcMap. ![enter image description here](http://i.stack.imgur.com/ziGm9.jpg) I can click on it but nothing happens. I haven't used this option before, so I'm not sure how long I've had this problem. Already checked on ESRI's site to see if there are any patches out for this issue and on GeoNet to see if anyone else is having this issue. Any ideas?


title: How to control Image Output Settings in QGIS?

text: I am exporting a QGIS map to an image format. I would like to have control over the output settings of the image (dpi), as currently it is only exporting the images at 96 dpi. Does anyone know if this is possible?


title: Styling a polygon layer using external graphics

text: I have a layer which is comprised of a number of different polygons symbolised according to their type. I have set up a style which successfully filters by type then applies a corresponding fill (in some cases overlays this with a hatch style) and this works great. Now I want to apply the same filters but in some instances apply an external graphic with a transparent background instead of a hatch style. I have used the sld code found in the Geoserver SLD Cookbook and this works fine on its own but when I try and use it in conjunction with a filter it fails to work. Has anyone come across this same issue and if so can you help me please? Thanks


title: Global grid, which reference system to use?

text: I have seen a couple global studies using WGS84 as geographic reference, with the problem that the area of raster cells is not constant. Which reference system would you use for a global analysis in order to have a constant cell size?


title: Is there a way to create raster by interpolating point values using GeoTools 8.2?

text: im new in using GeoTools (8.2). I was able to create shp from asc file which contains my data. Now my question is is there a simple way to create (interpolate) a simple shapefile (only points) to raster image? Thx in advance for all your answers. Best regards


title: UK cities boundries in Ordnance survey

text: Are there any city boundaries of UK cities in the Ordnance Survey?


title: Convert Landmask file into WKT Polygon

text: I have a land mask file in netCDF (also in ascii, if needed) which basically consists of a lat/long pair and a boolean (1 or 0) determining if that value is over land or not. The thing is, I'm not really sure how to use this in my Java applications. If I can create a WKT polygon from the boundary of land/water, that would be a huge help, but other wise I'm looking for answers on how to use the net CDF files in JAVA (I looked at the netCDF java stuff but wasn't sue how to use it for my purposes). I grab points from another source and must determine if they are on land or not.


title: ArcGIS 10 Search Services

text: I am trying to learn more about the new ArcGIS Server Search that comes with ArcGIS 10. Does anyone know where I can find detailed instructions on how to implement and use this feature? Specifically, how to build the seach index? Thanks, Seth


title: Migration from postgresql 8.4 to postgresql 9.0, postgis 1.5 involved

text: Im migrating from 8.4 to 9.0, in pg8.4 i have postgis 1.5.1 and in 9.0 im just installed 1.5.2. I decided to dump the postgis databases before the pg_upgrade command, and delete them, planning to create my new postgis_template and then create the databases:               pg_dump -C -h host -U postgres -f databasegis.dump databasegis     psql -U postgres < databasegis.dump      But I have errors :               ERROR:  could not access file "$libdir/liblwgeom": No such file or directory     ERROR:  function public.st_text(boolean) does not exist      What can i do? Thanks a lot.


title: Postgis: Is there a way to do an intersection query based on lines while taking into account the dateline?

text: _This is an exact copy of a question I asked on Stackoverflow a few weeks back and still haven't found a good answer for it. The original question is from here:original question._ I've been searching about this for a while and can't seem to find anything useful. In my database, I have an outline of the world coastline saved as lines. I want to take another line and determine where on the coastline that line intersects. My current query is like this:               SELECT st_intersection(makeline(p1.latlong, p2.latlong)::geography, the_geom::geography)     FROM worldcoastline, points p1, points p2      WHERE p1.id = 1 AND p2.id = 2      All is well and good, and everything seems to line up when I test it in KML and such, except when the dateline is involved. Instead of going across the dateline, it will go all the way around the world, and show me all these collisions it couldn't have had (across Russia, Europe, America). The most interesting thing I find though is if I throw in a st_intersects on the `makeline` as a geography, it returns false. I'm considering adding it to the where clause, but the problem is lines that cross the dateline and intersect land will return really whacky results anyways! Is there anything I can do in PostGIS to fix this problem? EDIT: I just found this in the reference: > ST_Intersection — (T) Returns a geometry that represents the shared portion > of geomA and geomB. **The geography implementation does a transform to > geometry to do the intersection and then transform back to WGS84**. So I guess that's why it's having dateline problems. Any alternatives or suggestions? EDIT: I just wanted to point out that I did add the st_instersects WHERE clause, but of course if a line does intersect the dateline as well as a coastline I will end up with a messed up ST_Intersection because it will still try to find the intersections going the long way around the world instead of across the dateline


title: Parameterized insert query to PostGIS fails

text: I have a PostGIS database table that I need to insert data into. I'm using Python 3 and the pg8000 interface to Postgres. Somehow, parameterized inserts do not work. If I do               iq="INSERT INTO points(location) VALUES(ST_GeographyFromText('SRID=4326;POINT(:lat :lon)'))"     cursor.execute(iq, {"lat":lat,"lon":lon})      or               iq="INSERT INTO points(location) VALUES(ST_GeographyFromText('SRID=4326;POINT(? ?)'))"     cursor.execute(iq, (lat, lon))      Then the query fails with a ProgrammingError:               pg8000.errors.ProgrammingError: (b'ERROR', b'XX000', b'parse error - invalid geometry')      This, however, works:               iq="INSERT INTO points(location) VALUES(ST_GeographyFromText('SRID=4326;POINT({0} {1})'))".format(lat,lon)     cursor.execute(iq)      What am I doing wrong in the parameterized query? Edit: Forgot to mention that, the database query works, if done with PgAdmin against the PostGIS base with some random point coordinates.


title: QGIS attribute action, open email client from Attribute table

text: I'm adding an `QgsAttributeAction` programmatically with the code below.               actions = iface.activeLayer.actions()     actions.addAction(QgsAction.OpenUrl, 'Email', 'mailto:[% "address" %]')      Once added, I inspect the `QgsAction` with               al = actions.listActions()     if len(al) > 0:       r = al[0].runable()       a = al[0].action()       n = al[0].name()       t = al[0].type()      Everything seems to be in order, but once I activate the action the "To" field in the mail client (outlook) is just the attribute field name (address) i.e., it's not substituted with the value of the attribute field. Another attempt is by looking at how the style file saves the action               <attributeactions>         <actionsetting action="mailto:[% "address" %]" capture="0" type="5" name="Email"/>     </attributeactions>      I've tried to mimic this by use the `QgsAttributeAction.readXML()`. I've manage to create a `QDomNode` but it fails to register as an action. Any thoughts or suggestions?


title: Data driven pages to line up?

text: When I create a map book, using data driven pages, when I place teh actual pages side by side there is some overlap, so the pages do not line up exactly. Is there a way to clip the overlap so that the pages are fairly seamless in the map book?


title: How do I add multiple OpenLayers Vector Layer markers to a map using AJAX, JSON, and jQuery

text: I am attempting to parse a JSON file and place markers on an OpenLayers based map application. I am using jQuery to perform the AJAX and iteration on the JSON file. I am having a problem with the code below. When inspecting the element, I can see in the console log that the latitude and longitude are being iterated through correctly. However, only one marker appears on the map at the 0,0 location. Can someone point me in the right direction to make this happen? A question similar to this was answered before (How to add multiple OpenLayers Vector Layer markers with data parsed from JSON file), however, the jQuery based answer was replaced with pure Javascript.                // Setup marker layer for all transmitters from DATA.json file     var markers = new OpenLayers.Layer.Markers( "Markers" );     map.addLayer(markers);          var size = new OpenLayers.Size(21,25);     var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);     var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png', size, offset);     $.getJSON('json/DATA.json', function(data) {         $.each(data.POI, function() {             markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(this.longitude,this.latitude),icon));             console.log(this.longitude, this.latitude);             markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(this.longitude, this.latitude),icon.clone()))         });     });      Here is a sample of the JSON               "LOCATION_A": {               "latitude": 10.2070,               "longitude": 25.3215,               "altitude": 1.2,               "status": "CLOSED",               "serverAddress": "192.168.0.1"      }      Any assistance with this matter is greatly appreciated! Thanks in advance!


title: return when asynchronous task has completed

text: I'd need a function to wait until the asynchronous task has completed and then return. In my example, I'm getting undefined returned from the "execute" method:               var result = myClass.execute(obj);  //undefined      How can I get the actual result returned?               return declare(null, {             constructor: function (obj) {                 this.mapServiceUrl = obj.mapServiceUrl;                 this.geometry_json = obj.geometry_json || undefined;                 this.where = obj.where || "1=1";             },                  execute: function () {                 var queryTask = new QueryTask(this.mapServiceUrl),                     query = new Query()                 ;                      query.returnGeometry = true;                 query.outFields = ["RDUWI"];                      if (this.geometry_json) {                     query.spatialRelationship = Query.SPATIAL_REL_INTERSECTS;                     query.geometry = Polygon(this.geometry_json);                 };                      query.where = this.where;                 var deferred = queryTask.execute(query);                 deferred.then(function (featureSet) {                     return featureSet;                 });             }         });


title: MODIS products reprojection problem

text: I am working with Modis Data (MOD02, MOD11, MOD03, MOD35). In particular I am looking for a solution to reproject the subdatasets below, into albers equal. Because I am using NLCD dataset which uses albers.(PROJ4: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs) area in 1000 m resolution. I will run a model for whole USA so I thats why I chose Albers. Can somebody help me? By the way I am aware of the MRT Swath but I want to do this in Python environment. So far I tried this but there was shifts in the pixels. (Alignment problems). gdalwarp -tr 1000 1000 -of GTIFF -tps -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs” -r cubic 'HDF4_EOS:EOS_SWATH: "MOD021KM.A2007181.1750.005.2010199133141.hdf”:MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB’ x.tif This one works for the MOD021KM but I couldn't work it out for the rest. I am confused. These are the subdatasets that I want to extract as geotiffs. HDF4_EOS:EOS_SWATH:"MOD021KM.A2007157.1700.005.2010197153031.hdf":MODIS_SWATH_Type_L1B:EV_250_Aggr1km_RefSB HDF4_EOS:EOS_SWATH:"MOD021KM.A2007157.1700.005.2010197153031.hdf":MODIS_SWATH_Type_L1B:EV_500_Aggr1km_RefSB HDF4_EOS:EOS_SWATH:"MOD03.A2007157.1700.005.2010196183307.hdf":MODIS_Swath_Type_GEO:SensorZenith HDF4_EOS:EOS_SWATH:"MOD03.A2007157.1700.005.2010196183307.hdf":MODIS_Swath_Type_GEO:SolarZenith HDF4_EOS:EOS_SWATH:"MOD11_L2.A2007157.1700.005.2007158222249.hdf":MOD_Swath_LST:LST HDF4_EOS:EOS_SWATH:"MOD11_L2.A2007157.1700.005.2007158222249.hdf":MOD_Swath_LST:View_angle HDF4_EOS:EOS_SWATH:"MOD35_L2.A2007157.1700.005.2010197153939.hdf":mod35:Cloud_Mask Sorry for the simple question. I am a newbie in the gdal and python.


title: Token Error with ArcGIS Server 10.2.2

text: I am using the ArPy example from esri sample. Script was working fine but then I updated my ArcGIS Server to 10.2.2 and I am getting this error: **Token Expired, Code: 498**. I have proper credentials as I am able to login to ArcGIS via Manager. http://resources.arcgis.com/en/help/main/10.1/index.html#/Example_Publish_service_definitions_listed_in_a_text_file/01540000061r000000/


title: Difference between OGRGeometry and Geometry in C++ and python

text: According to the documentation, there is a python version of OGRGeometry's C++ class (http://www.gdal.org/ogr/classOGRGeometry.html) which is http://gdal.org/python/osgeo.ogr.Geometry-class.html However, I don't understand why there is a lot of differences between these two class. For instance, in Python there aren't any import method (importFromWkt). Is it normal or maybe I'm on on the wrong way. Thanks


title: PostGIS Intermittent INDEX Performance

text: I have a table containing approximately 55 million data points ( point is a geometry with SRID 4326 ) and for my query I need to join this to an area table ( currently have ~1800 areas ) which contains a variety of different ranging from large polygons ( 2000 km square ) to fairly small ( small being about 100 km square ). The initial query selected by the user narrows down the initial 55 million points to around ~300,000 points dependent on date range etc they select. Then the join is done and dependent on what area set they have selected to use once the query is complete this normally narrows it down to a ~150,000. The problem I'm having is that some times the query just grinds to a halt and instead of taking the expected ~25 seconds it can take up to ~18 minutes. At this point have to usually do a VACUUM ANALYSE and then run a few queries before it starts behaving itself again. No data has been added, updated or removed from the data or areas tables at this point. I've played around with everything I can think of and this still seems to keep happening with no constancy too it. Both the data.point column and the area.polygon column have GIST INDEXES on them. I found removing the INDEX from the data.point column seemed to make things a bit more stable but is slower ~35 seconds normally. However removing an INDEX seems to be a very bad choice as should it not be helping not hindering? I'm using PostgreSQL 9.1.4 with PostGIS 1.5 Here's the query I'm running                   select * FROM data, area WHERE st_intersects (data.point, area.polygon) AND      (readingdatetime BETWEEN '1948-01-01' AND '2012-11-19') AND datasetid IN(3) AND      "polysetID" = 1 AND area.id IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)      EXPLAIN               Nested Loop  (cost=312.28..336.59 rows=5 width=2246) (actual time=1445.973..11557.824 rows=12723 loops=1)       Join Filter: _st_intersects(data.point, area.polygon)       ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.017..0.229 rows=35 loops=1)             Index Cond: ("polysetID" = 1)             Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))       ->  Bitmap Heap Scan on data  (cost=312.28..316.29 rows=1 width=297) (actual time=328.771..329.136 rows=641 loops=35)             Recheck Cond: ((point && area.polygon) AND (datasetid = 3))"             Filter: ((readingdatetime >= '1948-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-11-19 00:00:00'::timestamp without time zone))             ->  BitmapAnd  (cost=312.28..312.28 rows=1 width=0) (actual time=328.472..328.472 rows=0 loops=35)                   ->  Bitmap Index Scan on data_point_index  (cost=0.00..24.47 rows=276 width=0) (actual time=307.115..307.115 rows=1365770 loops=35)                         Index Cond: (point && area.polygon)                   ->  Bitmap Index Scan on data_datasetid_index  (cost=0.00..284.37 rows=12856 width=0) (actual time=1.522..1.522 rows=19486 loops=35)                         Index Cond: (datasetid = 3)     Total runtime: 11560.879 ms      My create tables               CREATE TABLE data     (       id bigserial NOT NULL,       datasetid integer NOT NULL,       readingdatetime timestamp without time zone NOT NULL,       value double precision NOT NULL,       description character varying(255),       point geometry,       CONSTRAINT "DATAPRIMARYKEY" PRIMARY KEY (id ),       CONSTRAINT enforce_dims_point CHECK (st_ndims(point) = 2),       CONSTRAINT enforce_geotype_point CHECK (geometrytype(point) = 'POINT'::text OR point IS NULL),       CONSTRAINT enforce_srid_point CHECK (st_srid(point) = 4326)     );          CREATE INDEX data_datasetid_index ON data USING btree (datasetid);     ALTER TABLE data CLUSTER ON data_datasetid_index;          CREATE INDEX "data_datasetid_readingDatetime_index" ON data USING btree (datasetid , readingdatetime );     CREATE INDEX data_point_index ON data USING gist (point);          CREATE INDEX "data_readingDatetime_index" ON data USING btree (readingdatetime );          CREATE TABLE area     (       id serial NOT NULL,       polygon geometry,       "polysetID" integer NOT NULL,       CONSTRAINT area_primary_key PRIMARY KEY (id )     )          CREATE INDEX area_polygon_index ON area USING gist (polygon);     CREATE INDEX "area_polysetID_index" ON area USING btree ("polysetID");     ALTER TABLE area CLUSTER ON "area_polysetID_index";      Hope that all makes a degree of sense if need to know anything else please ask. Short summary is really that the INDEXES seem to work some times but not others. Could anyone suggest anything I could try to work out what's happening? Thanks in advance. EDIT: Another example               select * FROM data, area WHERE st_intersects ( data.point, area.polygon) AND      (readingdatetime BETWEEN '2009-01-01' AND '2012-01-19') AND datasetid IN(1,3) AND      "polysetID" = 1 AND area.id IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)       Run on copy of table with point index               Nested Loop  (cost=0.00..1153.60 rows=35 width=2246) (actual time=86835.883..803363.979 rows=767 loops=1)       Join Filter: _st_intersects(data.point, area.polygon)       ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.021..16.287 rows=35 loops=1)             Index Cond: ("polysetID" = 1)             Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))       ->  Index Scan using data_point_index on data  (cost=0.00..1133.30 rows=1 width=297) (actual time=17202.126..22952.706 rows=33 loops=35)             Index Cond: (point && area.polygon)             Filter: ((readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone) AND (datasetid = ANY ('{1,3}'::integer[])))     Total runtime: 803364.120 ms      Run on copy of table without point index               Nested Loop  (cost=2576.91..284972.54 rows=34 width=2246) (actual time=181.478..235.608 rows=767 loops=1)       Join Filter: ((data_new2.point && area.polygon) AND _st_intersects(data_new2.point, area.polygon))       ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.149..0.196 rows=35 loops=1)             Index Cond: ("polysetID" = 1)             Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))       ->  Bitmap Heap Scan on data_new2  (cost=2576.91..261072.36 rows=90972 width=297) (actual time=4.808..5.599 rows=2247 loops=35)             Recheck Cond: ((datasetid = ANY ('{1,3}'::integer[])) AND (readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone))             ->  Bitmap Index Scan on "data_new2_datasetid_readingDatetime_index"  (cost=0.00..2554.16 rows=90972 width=0) (actual time=4.605..4.605 rows=2247 loops=35)                   Index Cond: ((datasetid = ANY ('{1,3}'::integer[])) AND (readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone))     Total runtime: 235.723 ms      As you can see the query is significantly slower when the point index is being used. EDIT 2 ( Pauls Suggested Query ):               WITH polys AS (       SELECT * FROM area       WHERE "polysetID" = 1 AND area.id IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)     )     SELECT *      FROM polys JOIN data ON ST_Intersects(data.point, polys.polygon)     WHERE datasetid IN(1,3)      AND (readingdatetime BETWEEN '2009-01-01' AND '2012-01-19');      EXPLAIN               Nested Loop  (cost=20.04..1155.43 rows=1 width=899) (actual time=16691.374..279065.402 rows=767 loops=1)       Join Filter: _st_intersects(data.point, polys.polygon)       CTE polys         ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.016..0.182 rows=35 loops=1)               Index Cond: ("polysetID" = 1)               Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))       ->  CTE Scan on polys  (cost=0.00..0.02 rows=1 width=602) (actual time=0.020..0.358 rows=35 loops=1)       ->  Index Scan using data_point_index on data  (cost=0.00..1135.11 rows=1 width=297) (actual time=6369.327..7973.201 rows=33 loops=35)             Index Cond: (point && polys.polygon)             Filter: ((datasetid = ANY ('{1,3}'::integer[])) AND (readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone))     Total runtime: 279065.540 ms


title: Programmatically Downloading Landsat Images

text: There was a very similar question to this one here , but it specified Landsat 8 images; I'd just like to be able to download ANY landsat images, given latitude and longitudes as inputs. Just wanted to see if anyone knew of something that does this already, although I don't hold much hope. I'm thinking of writing an application that uses usgs's Earth Explorer tool, but I don't have much experience with writing applications that make HTTP requests and such. If someone could at least help in going down this path, it would be something. I don't necessarily have any preference where language is concerned, but for what it's worth I'm writing the rest of the related software in C++. I've looked around at various sources of satellite imagery and other spatial data, but everything I've found falls into either of two categories; 1) it's a website that allows you to search for datasets/images and download them individually, or 2) is integrated into a GIS software system and doesn't appear to have any information for developers looking to do do this integration themselves. Ideally, I'd like something like google earth's static maps API, where I just specify latitude and longitude and I get the an image corresponding to my input. However, I don't want pretty pictures, I want geographic data, including data that I can't get from google earth (anything that's not visible light imagery). I'm getting away from the question-answer format of SE (I was mostly asking whether there were any existing solutions for automatically retrieving satellite imagery) but I might as well tell you my intent with all this; I'm trying to write software that reads in a variety of geographic data (imagery, height, temperature, etc) and attempts to classify the terrain into types (such as forest, grassland, ocean, etc). I'm not looking at a specific area; it's a general method and the software should be able to pull data from arbitrary locations.


title: plypython set-up

text: I am trying to set up plypython on the following system: Windows 8.1 PostgreSQL 9.2 64-bit PostGIS 2.0 Python 2.7.05 I have added plpython2.dll to \PostgreSQL\9.2\lib When I run:               CREATE EXTENSION plpython2u;      I get               ERROR:  could not load library "C:/Program Files/PostgreSQL/9.2/lib/plpython2.dll": %1 is not a valid Win32 application.      Is this a problem with the set-up or the .dll? Cheers.


title: Setting up pre-requisites for creating viewer using suite-sdk

text: I am creating a Viewer using the `suite-sdk` from the following tutorial, which says The Client SDK should be installed onto a machine where the development will occur. This is typically bundled with the OpenGeo Suite itself (under the sdk subdirectory) I have installed the OpenGeo and so I wanted to know,   1. where is the sdk sub directory in the Opengeo suite that i have installed located in Ubuntu?   2. I have JDK version, java version "1.6.0_24" OpenJDK Runtime Environment (IcedTea6 1.11.1) (6b24-1.11.1-4ubuntu3) OpenJDK 64-Bit Server VM (build 20.0-b12, mixed mode)   1. ant -version -> 1.8.2 **Edit:** I am not sure what the doc meant when it said extract the sdk dir to /opt/opengeo/suite/sdk, but I copied the /usr/share/opengeo-suite/sdk to the /opt/ dir, by creating opengeo, suite dirs. Is this the right method?


title: Use Python to change datasource from a feature class to a view built on the feature class

text: I'm trying to replace the datasource for all of the layers in my map document using Python but I keep running into an error. Right now, all of the layers point to a feature class called GISDB.DBO.TESTFEATURECLASS and I want them to point to a view that I built using the sdetable -o createview command. It has both the object id and geometry columns from the original feature class. It shows up in ArcCatalog with the same icon as a feature class and I can add it as a layer in ArcMap. > Runtime error <type 'exceptions.ValueError'>: Layer: Unexpected error               mapdoc = arcpy.mapping.MapDocument("CURRENT")     workspace_path = r'Database Connections\TestGIS.sde'     workspace_type = "SDE_WORKSPACE"     dataset_name = r'GISDB.DBO.TESTFEATURECLASS_VIEW'     for layer in arcpy.mapping.ListLayers(mapdoc):         layer.replaceDataSource(workspace_path, workspace_type, dataset_name)      How can I determine which property is causing the error?