title: Why can't I open postgis layers in QGIS?

text: I've set up a table in PostGIS using OSM data (converted using osm2pgsql). I know this database is good, because I'm able to import it and work with it in TileMill, but whenever I try to interact with any of the layers in QGIS (1.8.0, OSX), the Add button is greyed out (even if I put in a query). I've tested the connection and when I test my queries it shows that lots of rows were returned, so (again), the data's clearly there—I just can't open it up. CLARIFICATION: It's also worth noting that QGIS is able to see the database just fine. I can build queries and read out data in the Add PostGIS layer interface and everything comes back correctly. The only thing I can't do in the Add PostGIS Layer interface is add the layer. CLARIFICATION 2 (from the comments): `postgis_version` returns 2.0 USE_GEOS=1 USE_PROJ=1 USE_STATS=1. The extensions active on the database are plpgsql (1.0), postgis (2.0.3), and postgis_topology (2.0.3). Postgres and postgis were installed using Homebrew. The table does (now) have a primary key, but QGIS doesn't appear to find it. Screenshots: ![](http://ge.tt/api/1/files/4chLoGn/0/blob/x675?noinc=1) ![](http://ge.tt/api/1/files/4chLoGn/1/blob/x675?noinc=1)


title: Using pgr_drivingDistance with large ids

text: I have little experience with postgres and none with pgrouting. I'm trying to use the `pgr_drivingDistance` function with Open Street Map data (generated through `osm2po`), but whenever I run my query:               SELECT pgr_drivingDistance('SELECT id, osm_source_id::int4 AS source, osm_target_id::int4 AS target, cost::float8 AS cost FROM nyc_2po_4pgr', 85, 2, false, false);      I get `ERROR: integer out of range`—I assume because some of the osm_source_ids I'm getting in don't fit in an `int4`. Is that what's going on? I'm using a small subsection of the OSM data, so I assume I could squash the source_ids/target_ids into an `int4` range. How would I do that? Alternatively, is there a way to make the `pgr_drivingDistance` work with the larger source/target ids?


title: Setting coordinate system in PostGIS

text: The geographical information data that INEGI (part of the Mexico's goverment) provides are in the following coordinate system (or is projection?)               Lambert Conformal Conic, double standard parallel     Units: Meter     Coordinate System Type:     Datum: World Geodetic System of 1984.     False Origin     Northing: 0.0000     Easting: 2500000.000     Projection Parameters     Origin latitude: 12° 0’ 0.0’’     Origin longitude: -102° 0’ 0.0’’     Northern standard parallel: 29° 30’ 0.0’’     Southern standard parallel: 17° 30’ 0.0’’      I already put the information in the PostGIS/PostgreSQL database, now, How I can change the projection/coordinate system to that projection? Could you provide the code for doing it? And How I can change it to the coordinate system that Google use?


title: How to use ArcObjects to limit classes visible in legend?

text: Is there a method use ArcObjects for .NET to manipulate the "Only show classes that are visible in the current map extent" option for legend items? The GUI for this functionality is accessible via the legend properties: ![enter image description here](http://i.stack.imgur.com/pXX9Y.png) I have read through the ESRI documentation for the interfaces supported by Legend, and ILegendItem3 but have not found anything relevant.


title: Installing the development files libqgis-dev in Ubuntu 12.04

text: I am trying to install development files which include the headers and the libraries for qgis in Ubuntu 12.04. I used Ubuntu software center, synaptic and the terminal and am still getting the same result. The following packages have unmet dependencies: libqgis-dev: -Depends: libgeos-dev (>= 3.0.0) but 3.3.3-2~precise2 is to be installed -Depends: libgsl0-dev but it is not going to be installed -Depends: libqgis1.9.0 (= 1.9.0+git20120902+5a0d2f5~precise-ubuntugis1) but 1.9.0+git20120902+5a0d2f5~precise-ubuntugis1 is to be installed -Depends: libqt4-dev (>= 4.4.0) but 4:4.8.1-0ubuntu4.2 is to be installed -Depends: python-qt4 (>= 4.1.0) but 4.9.1-2ubuntu1 is to be installed -Depends: python-qt4-dev (>= 4.1.0) but it is not going to be installed -Depends: qt4-designer (>= 4.4.0) but 4:4.8.1-0ubuntu4.2 is to be installed Can anyone help me with this installation? Thank You.


title: How to create a customizabe table in QGIS Print Composer?

text: I would like to add a custom table to a map layout created in QGIS. I would like a fair bit of control over the look of the table (borders, background colours of the cells, etc.), and it appears that the "Add Attribute Table" tool doesn't let you control those things. I tried created what I wanted in Excel, and saving as an image, but of course it's pixelated. I would like a vector-based solution if at all possible. Any thoughts? I considered creating an .html table, and inserting that; is this the best method?


title: Arcpy Error 999999: Error executing points to line. Coordinates or Measures are out of bounds

text: I am selecting two points from a table and then doing a points to line function. I was originally getting this error because they were multipoints with multiple locations. I ran Multipart to Single Part on that layer and now am getting the same error on only a select few points but there are no multiples. What else can cause this error? Thanks! temporarytodraw holds 2 points and that's the result I get. I'm sorry I don't have much more to post than this. It selects points that match certain criteria and places them in the temporarytodraw but there is an if statement check to kick out any records that match less than 2 or more than 2.               arcpy.PointsToLine_management("temporarytodraw", "offset_line_thing", "", "", "NO_CLOSE")


title: gdal_translate KMLSUPEROVERLAY creates blank black image

text: I am trying to convert a geotiff to kmz file to view in GE. The geotiff is 878,720 KB. Using a 'regular kml' - the resolution is very low in GE. Having read around - it seems as conversion to a superoverlay is the best method? I have chosen to use gdal_translate above other methods (e.g. MapTiler or Arc's tool that costs) because I have been using the GDAL python library already and I am trying to automate the entire process. To test, I type into UNIX:               gdal_translate -of KMLSUPEROVERLAY filepath/2010_full_mosaic.tif filepath/2010_full_mosaic.kmz      A kml file is created....but on loading it into GE it is blank (black) - although extent/projection appear correct. Would anyone be able to advise what I may have done wrong? When using the gdalinfo command, I am provided the following details:               Driver: GTiff/GeoTIFF     Files: Filepath/2010_full_mosaic.tif            Filepath/2010_full_mosaic.tif.ovr            Filepath/2010_full_mosaic.tif.aux.xml     Size is 10572, 22630     Coordinate System is:     PROJCS["WGS 84 / UTM zone 36S",         GEOGCS["WGS 84",              DATUM["WGS_1984",                 SPHEROID["WGS 84",6378137,298.257223563,                     AUTHORITY["EPSG","7030"]],                 AUTHORITY["EPSG","6326"]],             PRIMEM["Greenwich",0],             UNIT["degree",0.0174532925199433],             AUTHORITY["EPSG","4326"]],         PROJECTION["Transverse_Mercator"],         PARAMETER["latitude_of_origin",0],         PARAMETER["central_meridian",33],         PARAMETER["scale_factor",0.9996],         PARAMETER["false_easting",500000],         PARAMETER["false_northing",10000000],         UNIT["metre",1,             AUTHORITY["EPSG","9001"]],         AUTHORITY["EPSG","32736"]]     Origin = (390313.896406000014395,8164584.107470000162721)     Pixel Size = (25.000000000000000,-25.000624509551976)     Metadata:       AREA_OR_POINT=Area     Image Structure Metadata:       INTERLEAVE=BAND     Corner Coordinates:     Upper Left  (  390313.896, 8164584.107) ( 31d58'18.28"E, 16d35'54.43"S)     Lower Left  (  390313.896, 7598819.975) ( 31d56'22.37"E, 21d42'38.47"S)     Upper Right (  654613.896, 8164584.107) ( 34d26'57.65"E, 16d35'45.40"S)     Lower Right (  654613.896, 7598819.975) ( 34d29'41.00"E, 21d42'26.42"S)     Center      (  522463.896, 7881702.041) ( 33d12'49.09"E, 19d 9'27.96"S)     Band 1 Block=10572x1 Type=Float32, ColorInterp=Gray      Overviews: 5286x11315, 2643x5658, 1322x2829, 661x1415, 331x708, 166x354      Becky


title: Is there a no point symbology option in ArcGIS?

text: I'm sure it exists, but is there a way to set a point layer to have no symbology. This works best for cities as you just want to display the name, and not necesserily show the location. The way I do it now, is by picking a symbol and making the fill and outline to have "no color"


title: What is the smallest contour interval possible when derived from a DTM?

text: I am converting a 2m DTM to contours, and would like to have the smallest interval possible. Is there a theoretical limit, or even a rule of thumb for defining contour intervals based on the resolution of the source data?


title: buildings shadow areas calculation

text: is it possible to calculate shadow areas of buildings or simulate shadows of buildings in a city using e-cognition software, with the DSM of the buildings and the sun angle and azimuth known?....if not, could you please tell me the software that can do this work. thanks in advance...


title: Copy a polygon geometry in Arcpy

text: This to me sounded very simple so I am probably missing something basic here. What I need to do is to create a copy of a polygon shapefile to a new shapefile with slightly different attributes. Where I got stuck is in copying the geometry - Geometry doesn't have to change at all in the new shapefile. The actual 'copying' of the polygon geometry is not the problem. I loop through the input shapefile:               > for row in rows:     >     partnum = 0     >     feat=row.getValue("SHAPE")      And then I enter this shape in the output feature class:               outFeat.Shape=feat      Problem is, I also want to calculate the Area, Perimeter and centroid of the polygon and I can't seem to 'cast' the **_feat_** variable to a polygon so I can then get the polygon properties (area, centroid, etc). Can you do that somehow? One of the things I tried was to create the polygon from scratch, i.e. loop though the vertices of the input polygon and create a new polygon/array to use for the new polygon. This method apart from being IMO over the top, it also had the problems with the multi-part polygons. My polygons ended up corrupted - although I followed the steps outlined here I didn't dwell too much with it since what I did instead is create a dummy Polygon which I then passed to it the current geometry value in the cursor.               # Create a dummy polygon                 dummyArray=arcpy.Array()                 dummyPoint = arcpy.Point()                 coordList = [[[1,2], [2,4], [3,7]],[[6,8], [5,7], [7,2], [9,5]]]                 for feature in coordList:                     for coordPair in feature:                         dummyPoint.X = coordPair[0]                         dummyPoint.Y = coordPair[1]                         dummyArray.add(dummyPoint)                     dummyArray.add(dummyArray.getObject(0))                      parcel=arcpy.Polygon(dummyArray)                 parcel=feat ### --> This is where I replace the dummy geometry value with the current value in the loop                 parcel_c=arcpy.Point()                 parcel_c=parcel.centroid # I can now get to the polygon properties                 c_x=parcel_c.X                 c_y=parcel_c.Y      Although this works fine, I am not really happy with it as its a bit of a hacking and -as mentioned in the beginning- I think I am missing something vital! Any insights?


title: ArcGIS JS API Pop-Up w/TabContainer

text: I am attemtping to call a dojo tab container within an ESRI pop-up. The code that I have (see below) is able to create a pop-up with a tab container and the contents, but also creates three 'dead' dropdown menus. Is there a way to create the pop-up without the dropdown menus?                 var map;       var resizeTimer;       var identifyTask,identifyParams;            function init() {                 //setup the popup window          var popup = new esri.dijit.Popup({           fillSymbol: new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID, new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new dojo.Color([255,0,0]), 2), new dojo.Color([255,255,0,0.25]))         }, dojo.create("div"));              map = new esri.Map("map",{           infoWindow:popup,         });                   dojo.connect(map,"onLoad",mapReady);              //Add the imagery layer to the map. View the ArcGIS Online site for services http://arcgisonline/home/search.html?t=content&f=typekeywords:service             var basemap = new esri.layers.ArcGISTiledMapServiceLayer("http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer");         map.addLayer(basemap);              var landBaseLayer = new esri.layers.ArcGISDynamicMapServiceLayer("http://<someServer>/ArcGIS/rest/services/test/MapServer");         map.addLayer(landBaseLayer)       }                function mapReady(map){                 dojo.connect(map,"onClick",executeIdentifyTask);                 //create identify tasks and setup parameters             identifyTask = new esri.tasks.IdentifyTask("http://<someServer>/ArcGIS/rest/services/test/MapServer");                 identifyParams = new esri.tasks.IdentifyParameters();            identifyParams.tolerance = 3;            identifyParams.returnGeometry = true;            identifyParams.layerIds = [0];            identifyParams.layerOption = esri.tasks.IdentifyParameters.LAYER_OPTION_ALL;            identifyParams.width  = map.width;            identifyParams.height = map.height;                      //resize the map when the browser resizes            dojo.connect(dijit.byId('map'), 'resize', map,map.resize);                }                 function executeIdentifyTask(evt) {                 identifyParams.geometry = evt.mapPoint;                 identifyParams.mapExtent = map.extent;         var deferred = identifyTask.execute(identifyParams);              deferred.addCallback(function(response) {                // response is an array of identify result objects               // Let's return an array of features.           return dojo.map(response, function(result) {             var feature = result.feature;             feature.attributes.layerName = result.layerName;             console.log(feature.attributes.OBJECTID);             var template = new esri.InfoTemplate();             template.setTitle("Hello World");             template.setContent(getWindowContent);             feature.setInfoTemplate(template);             return feature;           });         });                   // InfoWindow expects an array of features from each deferred         // object that you pass. If the response from the task execution          // above is not an array of features, then you need to add a callback         // like the one above to post-process the response and return an         // array of features.         map.infoWindow.setFeatures([ deferred ]);         map.infoWindow.show(evt.mapPoint);       }            function getWindowContent(graphic) {         //make a tab container          var tc = new dijit.layout.TabContainer({         });              return tc.domNode;       }            dojo.addOnLoad(init);


title: Problems with labels halo in OpenLayers

text: Draw label as follows:               var styleMarkersLabels = new OpenLayers.Style(                     {                         strokeWidth: 2,                         labelYOffset: 14,                         label: "${label}",                         fontColor: 'red',                         fontSize: 16,                         fontWeight: "bold",                         labelOutlineColor: "black",                         labelOutlineWidth: 1                     });                 var vectorPoint = new OpenLayers.Layer.Vector("Узлы (надписи)",                 {                     styleMap: new OpenLayers.StyleMap(                     { "default": styleMarkersLabels,                          "select": { pointRadius: 20}                     })                 });                 map.addLayer(vectorPoint);                      function addPoint(lon, lat, title, ident, layr){                     var ttt = new OpenLayers.LonLat(parseFloat(lon), parseFloat(lat));                     ttt.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));                     for (var k = 0; k < layr.features.length; k++) {                         if(layr.features[k].attributes.PointId==ident) {                             layr.features[k].move(ttt);                             layr.features[k].attributes.label=title;                             return false;                         }                     }                     var point0 = new OpenLayers.Geometry.Point(parseFloat(lon), parseFloat(lat));                     point0.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));                     layr.addFeatures(new OpenLayers.Feature.Vector(point0, { label: title, name: title, PointId: ident }));                 }                      var lat2, lon2, title, ident;                 for (l = 0; l < nodesLabels_Count; l++) {                     lat2 = nodesLabels_arr[l]["points"][0]["lat"];                     lon2 = nodesLabels_arr[l]["points"][0]["lon"];                     title = nodesLabels_arr[l]["title"];                     ident = nodesLabels_arr[l]["ident"];                     addPoint(lon2, lat2, title, ident, map.layers[4]); // Here adds the label                 }      ![enter image description here](http://i.stack.imgur.com/aI8Bg.png) But the halo is shifted slightly. Screen: http://softez.pp.ua/Image.PNG How to fix it? Thanks!


title: Geographical metadata for traditional and non-traditional media outlets (print, tv, radio, blogs, etc)

text: A data analytics company I am consulting for has a media metadatabase of about 4,000 print publications with associated names and another 56,000 other media entities (blog, forum, news, radio, tv) with associated URLs. They have about 200 data entry staff to update metadata for these media. The specific focus of this question is **what is the best method of finding and capturing specific geo metadata, given the unique identifiers I have?** **for print media the unique identifier is only the company's proprietary publicationid number. you can also get a unique identifier by combining the publicationtitle with the country**               publicationid   publicationtitle        country          6541449         Auto Bild               Germany     15728025        Auto Bild               Slovakia     15703996        Auto Bild               Hungary     6542021         Auto Bild               Switzerland     6542549         Auto Bild               Mexico     6541450         Auto Bild allrad        Germany      The best solution I've come up with so far to semi-automate finding more specific geo-metadata for these is to utilize the google or facebook api for each of these rows, using those two fields (publicationtitle and country) as a query. **The query would hopefully return the offical facebook page of the queried media outlet, and my script could save the metadata of all these facebook pages, hopefully containing the specific geo metadata I seek (i.e. coordinates, address, city, state, etc)** -- **For non-print media, we have the same option to combine mediaOutletName with country.** **We also have a unique value for the domain property, which likely is more useful.** **What associated metadata can one find with readily available APIs using a domain as the input? I am specifically interested in metadata of other properties currently in our system (though incomplete). Other properties listed below** ...... extra info for the curious... If you can provide information for open source or proprietary APIs, json, xml, web scraping business logic for how to help automate the initial capturing and future updating of any of the following fields, that would be wonderful. Or if you can think of any other properties we could associate with our media outlet objects, then I'd be happy to consider. It could be useful to open source a project like this, and I may discuss this as an approach with the company. The fields of metadata for the print (many incomplete/out of date):               publicationid     publicationtitle     country     language     general focus (ex. consumer)     specific focus (ex. Fashion & Lifestyle)     image path     city     supplements     print format     subscription source     freqdays     freqinterval     time     description     [...many fields for accounting of subscriptions to physical print media (for in-house scanning of images)]      And the metadata for non-print media types:               mediaOutletName (unique string)     country     language     logo     modified by (employee who last edited)     modified at (time of last edit)     domain     lastcheck      media type (ex. tv, online news, etc)     specific focus (ex. Fashion & Lifestyle)     dataProvider1ID     dataProvider2ID     hasRSSFeed     pageImpressions     uniqueVisitors     metricSource     metricUpdatedTimestamp


title: Identify Polygons on raster image

text: I have OS street view mapping as raster images. Is it possible (without the obvious of tracing each dwelling) to identify all polygons on each tile. To clarify, My ultimate aim is to produce a shapefile layer of polygons for each component....particularly dwellings in order to create a 260m buffer around each which is easy done. Using QGIS


title: How to clip all shapefiles, layers, etc associated with global map of Mars

text: I'm working on a project where I've found an open source NASA .mxd of shapefiles, layers, metadata, geology, a 30x30 graticle that splits the lat and long into 30 degree square polygons, etc. I only want to use data that exists between the equator and 30 degrees latitude in both northern and southern hemisphere. How would you clip this data without having to select a region to clip by hand for all the different raster and vector layers (there are a lot of them)? Thanks


title: Ways to change ArcMap 10 editing tool back to ArcMap 9.3 version

text: I'm trying to get used to using ArcGIS 10, but in the meantime I still need to actually get work done. I stumbled across a way to change the editor tool from its 10 version back to the 9.3 version (with the drop-down task bar etc). Now I can't find out ow to do it anywhere, -1 internet for me. Can anyone can give me the steps to change the editor toolbar back to the 9.3 version?


title: merging ecw images

text: Is anyone able to recommend a program/process/method for merging ecw images together. I have approximately 2000 ecw images (tile) that 1.5x1.5km each which together and put to combine them all together to make one mosaics there is no overlap but each tile abuts exactly to the next. Collectively they equate to approximately 100gb. Any ideas?


title: Upgrade strategy for geoserver 2.1.4 to 2.5.x

text: I have a 2.1.4. geoserver installation and I wish to upgrade to the latest version which is currently 2.5.2. Can I do this at once by installing the latest version or is it wise to do an intermediate step via geoserver 2.2. and upgrade further from there?   I noticed for instance that since geoserver 2.2.0 some security changes have been applied. Is it wise to get past this upgrade first and then upgrade further?


title: Is it possible to rotate the geographic projection in QGIS?

text: I'm trying to create a Qgis plugin that can rotate a raster layer in the canvas so that the north is pointing at another direction. Unlike the thread below, I'm open to defining a custom projection. > Is it possible to rotate the *canvas* in QGIS The closest feel to what I wanted was using: > qgis.utils.iface.mapCanvas().rotate(x) But it literally rotates the map canvas, causing tools such as the map panning and map extents go haywire. This suggests that a custom map projection might be the best solution out there. However, I'm quite new to GIS so I'm not too sure if it is actually possible to derive a custom projection by rotating the existing crs. Any guidance will be deeply appreciated. Also print composer is not really an answer.


title: Version problems with Geotools CSV2SHP example

text: I'm trying to use the CSV2SHP tutorial. I spent some time trying to get it to run but using v10.x of Geotools causes the shapefile to be output with no geometries. For instance if I open it in QGIS, the layer is present and the attributes can be viewed in the table but there are no points on the map. QGIS reports the extent as xMin,yMin -1,-1 : xMax,yMax 0,0. I've just tried using version 9.0 and it works as expected. Can anyone explain what changed between 9 and 10 to cause this to break? It's quite puzzling as the tutorial seems to have been written for v10. I've pasted the full code from the tutorial (including some minor edits) below.               import java.io.BufferedReader;     import java.io.File;     import java.io.FileReader;     import java.io.Serializable;     import java.util.HashMap;     import java.util.Map;          import org.geotools.data.DefaultTransaction;     import org.geotools.data.Transaction;     import org.geotools.data.shapefile.ShapefileDataStore;     import org.geotools.data.shapefile.ShapefileDataStoreFactory;     import org.geotools.data.simple.SimpleFeatureSource;     import org.geotools.data.simple.SimpleFeatureStore;     import org.geotools.feature.DefaultFeatureCollection;     import org.geotools.feature.simple.SimpleFeatureBuilder;     import org.geotools.feature.simple.SimpleFeatureTypeBuilder;     import org.geotools.geometry.jts.JTSFactoryFinder;     import org.geotools.referencing.crs.DefaultGeographicCRS;     import org.geotools.swing.data.JFileDataStoreChooser;     import org.opengis.feature.simple.SimpleFeature;     import org.opengis.feature.simple.SimpleFeatureType;          import com.vividsolutions.jts.geom.Coordinate;     import com.vividsolutions.jts.geom.GeometryFactory;     import com.vividsolutions.jts.geom.Point;          /**      * This example reads data for point locations and associated attributes from a comma separated text     * (CSV) file and exports them as a new shapefile. It illustrates how to build a feature type.      * <p>      * Note: to keep things simple in the code below the input file should not have additional spaces or      * tabs between fields.      */     public class Csv2Shape {              public static void main(String[] args) throws Exception {                  File file = JFileDataStoreChooser.showOpenFile("csv", null);             if (file == null) {                 return;             }                  /*              * We use the DataUtilities class to create a FeatureType that will describe the data in our              * shapefile.              *               * See also the createFeatureType method below for another, more flexible approach.              */             final SimpleFeatureType TYPE = createFeatureType();                   /*          * We create a FeatureCollection into which we will put each Feature created from a record          * in the input csv data file          */         DefaultFeatureCollection collection = new DefaultFeatureCollection();              /*          * GeometryFactory will be used to create the geometry attribute of each feature (a Point          * object for the location)          */         GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);              SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);              BufferedReader reader = new BufferedReader(new FileReader(file));         try {             /* First line of the data file is the header */             String line = reader.readLine();             System.out.println("Header: " + line);                  for (line = reader.readLine(); line != null; line = reader.readLine()) {                 if (line.trim().length() > 0) { // skip blank lines                     String tokens[] = line.split("\\,");                          double latitude = Double.parseDouble(tokens[0]);                     double longitude = Double.parseDouble(tokens[1]);                          System.out.println("lat : "+latitude);                     System.out.println("lon : "+longitude);                          String name = tokens[2].trim();                     int number = Integer.parseInt(tokens[3].trim());                          /* Longitude (= x coord) first ! */                     Point point = geometryFactory.createPoint(new Coordinate(longitude, latitude));                          featureBuilder.add(point);                     featureBuilder.add(name);                     SimpleFeature feature = featureBuilder.buildFeature(null);                     collection.add(feature);                 }             }         } finally {             reader.close();         }              /*          * Get an output file name and create the new shapefile          */         File newFile = getNewShapeFile(file);              ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();              Map<String, Serializable> params = new HashMap<String, Serializable>();         params.put("url", newFile.toURI().toURL());         params.put("create spatial index", Boolean.TRUE);              ShapefileDataStore newDataStore = (ShapefileDataStore) dataStoreFactory.createNewDataStore(params);         newDataStore.createSchema(TYPE);              /*          * You can comment out this line if you are using the createFeatureType method (at end of          * class file) rather than DataUtilities.createType          */         newDataStore.forceSchemaCRS(DefaultGeographicCRS.WGS84);              /*          * Write the features to the shapefile          */         Transaction transaction = new DefaultTransaction("create");              String typeName = newDataStore.getTypeNames()[0];         SimpleFeatureSource featureSource = newDataStore.getFeatureSource(typeName);              if (featureSource instanceof SimpleFeatureStore) {             SimpleFeatureStore featureStore = (SimpleFeatureStore) featureSource;                  featureStore.setTransaction(transaction);             try {                 featureStore.addFeatures(collection);                 transaction.commit();                  } catch (Exception problem) {                 problem.printStackTrace();                 transaction.rollback();                  } finally {                 transaction.close();             }             System.exit(0); // success!         } else {             System.out.println(typeName + " does not support read/write access");             System.exit(1);         }     }          /**      * Prompt the user for the name and path to use for the output shapefile      *       * @param csvFile      *            the input csv file used to create a default shapefile name      *       * @return name and path for the shapefile as a new File object      */     private static File getNewShapeFile(File csvFile) {         String path = csvFile.getAbsolutePath();         String newPath = path.substring(0, path.length() - 4) + ".shp";              JFileDataStoreChooser chooser = new JFileDataStoreChooser("shp");         chooser.setDialogTitle("Save shapefile");         chooser.setSelectedFile(new File(newPath));              int returnVal = chooser.showSaveDialog(null);              if (returnVal != JFileDataStoreChooser.APPROVE_OPTION) {             // the user cancelled the dialog             System.exit(0);         }              File newFile = chooser.getSelectedFile();         if (newFile.equals(csvFile)) {             System.out.println("Error: cannot replace " + csvFile);             System.exit(0);         }              return newFile;     }          /**      * Here is how you can use a SimpleFeatureType builder to create the schema for your shapefile      * dynamically.      * <p>      * This method is an improvement on the code used in the main method above (where we used      * DataUtilities.createFeatureType) because we can set a Coordinate Reference System for the      * FeatureType and a a maximum field length for the 'name' field dddd      */     private static SimpleFeatureType createFeatureType() {              SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();         builder.setName("Location");         builder.setCRS(DefaultGeographicCRS.WGS84); // <- Coordinate reference system                  // add attributes in order             builder.add("Location", Point.class);             builder.length(15).add("Name", String.class); // <- 15 chars width for name field              // build the type         final SimpleFeatureType LOCATION = builder.buildFeatureType();              return LOCATION;        }          }


title: Geotools writing shapefiles without geometries

text: When I am writing a shapefile according to CSV2SHP tutorial, it contains features without geometries. I use version 10-SNAPSHOT. For example the result of ogrinfo looks like below:               D:\>ogrinfo locations1.shp locations1     INFO: Open of `locations1.shp'           using driver `ESRI Shapefile' successful.          Layer name: locations1     Geometry: Point     Feature Count: 9     Extent: (0.000000, 0.000000) - (-1.000000, -1.000000)     Layer SRS WKT:     GEOGCS["WGS84(DD)",         DATUM["WGS84",             SPHEROID["WGS84",6378137.0,298.257223563]],         PRIMEM["Greenwich",0.0],         UNIT["degree",0.017453292519943295],         AXIS["Geodetic longitude",EAST],         AXIS["Geodetic latitude",NORTH]]     Name: String (15.0)     Number: Real (15.13)     OGRFeature(locations1):0       Name (String) = Trento       Number (Real) = 140.0000000000000          OGRFeature(locations1):1       Name (String) = St Paul       Number (Real) = 125.0000000000000          OGRFeature(locations1):2       Name (String) = Bangkok       Number (Real) = 150.0000000000000          OGRFeature(locations1):3       Name (String) = Ottawa       Number (Real) = 200.0000000000000          OGRFeature(locations1):4       Name (String) = Minneapolis       Number (Real) = 350.0000000000000          OGRFeature(locations1):5       Name (String) = Lausanne       Number (Real) = 560.0000000000000          OGRFeature(locations1):6       Name (String) = Victoria       Number (Real) = 721.0000000000000          OGRFeature(locations1):7       Name (String) = Cape Town       Number (Real) = 550.0000000000000          OGRFeature(locations1):8       Name (String) = Sydney       Number (Real) = 436.0000000000000      Where could be the problem?


title: Adding special characters (e.g. $) to double value in Identify window?

text: I don't know what I want is possible in ArcMap? I have a field whose type is "Double". The value in it is accurate to 2 decimal places. When I use Identify to see the detail information of a feature, I wonder whether there is a way to add a $ in front of the double value? What I see in Identify is as follows: ![Identify window](http://i.stack.imgur.com/j8gWA.jpg) What I intend to see in the CONTRACT_SUM field is $24222568, is this possible?


title: Draw linestring programmatically using draw control

text: I want to draw linestring feature programmatically which can be done from draw control of openlayers which can inherit default behavior of draw control that is modification, selection, undo, redo etc I found one method in openlayers but seems not working for me               drawControls.line.drawFeature(geometry)


title: How to edit a layer in Geoserver, when we change table schema in a database (SQL Server)?

text: So, is there a way to change/edit SQL Server layer definition in Geoserver, when a table schema has been changed in a database ? Or you have to delete that layer and then build a new one ? Thanks!


title: How tow use an ArcGISCache layer along with an ArcGIS93Rest layer?

text: I want to have two layers in my OpenLayers map. The BaseLayer, as well as the overlay layer are being served by ArcGIS Server. I have a problem with ArcGIS93Rest when using ArcGISCache as basemap. The problem is `ArcGIS93Rest` doesn't show up and there is no HTTP request to the service. Is it possible to use `ArcGISCache` and `ArcGIS93Rest` together? Here is the code :               var baseLayer = new OpenLayers.Layer.ArcGISCache("AGSCache", "http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer", {         isBaseLayer: true,                              resolutions: resolutions,         tileSize: new OpenLayers.Size(layerInfo.tileInfo.cols, layerInfo.tileInfo.rows),         tileOrigin: new OpenLayers.LonLat(layerInfo.tileInfo.origin.x, layerInfo.tileInfo.origin.y),         maxExtent: layerMaxExtent,         projection: 'EPSG:' + layerInfo.spatialReference.wkid     });     var forestLayer = new OpenLayers.Layer.ArcGIS93Rest("ForestLayer", "http://192.168.1.36:6080/arcgis/rest/services/forest_editing/MapServer/export", {         layers: "0,1,2",         TRANSPARENT: true     });


title: Error in displaying wfs layer from mapserver/postgis

text: I have read and tested all the relative posts but I still don't have a solution for my problem. I have a layer suscribers (points) stored in postgis database with srid=32232. I would like to display it in openlayers in geographic projection 4326. I can see the response with firebug but no display. my code : mapfile * * *               MAP      IMAGETYPE PNG      EXTENT -180 -90 180 90      UNITS DD      SIZE 1500 650      IMAGECOLOR 255 255 255      FONTSET "fonts/fonts.list"      SYMBOLSET "symbols/symbols35.sym"          WEB IMAGEPATH "tmp/tmp_img/"      IMAGEURL "**/guide/tmp/tmp_img/"      HEADER "query_header.html"      FOOTER "query_footer.html"      METADATA          "wfs_title" "WFS Guide" ## REQUIRED     "wfs_onlineresource" "http://*********/cgi-bin/mapserv.exe?" ## Recommended     "wfs_srs" "EPSG:4326" ## Recommended     "wfs_abstract" "This is the wfs service of guide" ## Recommended     "wfs_enable_request" "*"  # necessary     END     END     PROJECTION      "init=epsg:4326"      END          LAYER # suscribers begin here      PROJECTION       "init=epsg:32232"      END      METADATA      "wfs_title" "suscribers" ##REQUIRED      "wfs_srs" "EPSG:4326" ## REQUIRED      "gml_include_items" "all" ## Optional (serves all attributes for layer)      "gml_featureid" "gid" ## REQUIRED      "wfs_enable_request" "*"       END      CONNECTIONTYPE postgis      NAME suscribers      CONNECTION "user=postgres dbname=postgis20 password=pass" DATA "geom FROM                  public.suscribers USING UNIQUE gid USING SRID=32232"     STATUS ON      TYPE POINT      DUMP TRUE          END # suscribers ends here          END # end of map      openlayers * * *                     var map = new OpenLayers.Map('map',{ controls: [], allOverlays: false, tileSize:       new OpenLayers.Size(100, 100), maxExtent: OpenLayers.Bounds.fromString("-180,-90,180,90"), projection: new OpenLayers.Projection("EPSG:4326"), displayProjection: new OpenLayers.Projection("EPSG:4326") //mouse position });                var basic = new OpenLayers.Layer.WMS("Basic", "http://vmap0.tiles.osgeo.org/wms/vmap0", { layers: 'basic' },{singleTile: true, ratio: 1, transitionEffect: 'resize', opacity: 0.5 });               var wfsSuscribers = new OpenLayers.Layer.Vector("States", {protocol: new OpenLayers.Protocol.WFS({     url: "http://*********/cgi-bin/mapserv.exe?map=C:/ms4w/Apache/htdocs/routing/index.map&",     featureType: "suscribers",     featureNS: "http://****/routing"     }),     projection: new OpenLayers.Projection("EPSG:4326"), strategies: [new OpenLayers.Strategy.BBOX()]     });                map.addLayers([basic, wfsSuscribers]); map.zoomToMaxExtent(); map.setCenter(new OpenLayers.LonLat(11.52, 3.83), 3);      * * *               <gml:boundedBy>         <gml:Box srsName="EPSG:4326">             <gml:coordinates>11.416136,3.714272 11.554812,3.969178</gml:coordinates>         </gml:Box>       </gml:boundedBy>     <gml:featureMember>       <ms:suscribers fid="suscribers.2">         <gml:boundedBy>             <gml:Box srsName="EPSG:4326">                 <gml:coordinates>11.521489,3.869029 11.521489,3.869029</gml:coordinates>             </gml:Box>         </gml:boundedBy>         <ms:msGeometry>         <gml:Point srsName="EPSG:4326">           <gml:coordinates>11.521489,3.869029</gml:coordinates>         </gml:Point>         </ms:msGeometry>         <ms:gid>2</ms:gid>         <ms:osm_id>331136682</ms:osm_id>         <ms:addr_house></ms:addr_house>         <ms:amenity></ms:amenity>         <ms:name>Yaoundé</ms:name>         <ms:barrier></ms:barrier>      Thank you.


title: Does ArcObjects ConstructUnion operation sometimes fail to merge touching parcels?

text: I have some code which identifies touching parcels and uses ConstructUnion to produce a single polygon from them. Just occasionally the touching parcels remain unmerged even though they have been selected with esriSpatialRelTouches. Is this a known bug? Or, the reason might be that I am using a geometryBag to store the parcels to be merged. The constructUnion documentation states that it does not work on a geometrybag but I cannot find a single example on the web of anything else being used instead - can anyone advise what I should be using. Val.


title: Why does compressed & zipped/unzipped file geodatabase give Open Failed Error opening feature class Cannot access secured data?

text: In Arc 10.2, a colleague compressed a file geodatabase and then zipped it and sent it to my via email. I unzipped it and was unable to add any of the feature classes into my map document. I received this error: "Open Failed. Error opening feature class. Cannot access secured data." I uncompressed it in ArcCatalog and got the same error when trying to add the data to ArcMap. This also happened when the colleague zipped the uncompressed geodatabase. Thanks in advance for any help on this.


title: Multiple input shapefiles for ogr2ogr sqlite spatial operations

text: I'm running up against an issue doing a spatial operation requiring multiple input shapefiles using ogr2ogr and sqlite. I have:               ogr2ogr -f "ESRI Shapefile" output.shp input.shp -dialect sqlite -sql "CREATE VIRTUAL TABLE 'overlap' USING VirtualOGR('overlap.shp'); SELECT ST_SymDifference('input','overlap')"      Where 'output' is my output shape, 'input' is one input shape and 'overlap' is the other input shape. I can't seem to get any output, or if I do, it's just a DBF file? Any assistance? I suspect it's to do with the Virtual Table not having any geometry but I can't seem to do a "RecoveryGeometryColumn" that actually works. Both 'input' and 'overlap' have the same projection and overlap nicely in QGIS. I'm on version "GDAL 1.10.0, released 2013/04/24" Thanks! Edit: Fixed typos in my sqlite command Edit2:               ogr2ogr -f "ESRI Shapefile" output.shp input.shp -dialect sqlite -sql "SELECT ST_SymDifference(GEOMETRY,'overlap.sqlite.GEOMETRY'),* FROM input"      Just returns the input shape. I think I've completely missed something obvious...


title: Coordinates misplaced in openlayers

text: I am using a timeline from timeline-github. I used json for showing temporal data in timeline but the features are clustered into (0,0) i.e features are misplaced. I think it's because of projection problem but don't know how to debug. Any help and suggestion is appreciated. My code is               var map, timeline;     OpenLayers.ImgPath = "green_theme/green_img/";          function init() {         OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";         map = new OpenLayers.Map('map', {             projection: new OpenLayers.Projection("EPSG:4326"),             displayProjection: new OpenLayers.Projection("EPSG:4326"),             numZoomLevels: 20,             controls: [             new OpenLayers.Control.Navigation(), new OpenLayers.Control.PanZoomBar(), new OpenLayers.Control.LayerSwitcher(), new OpenLayers.Control.MousePosition()]         });         mapnik = new OpenLayers.Layer.OSM.Mapnik("Mapnik");         map.addLayer(mapnik);         map.setCenter(new         OpenLayers.LonLat(84.39226, 27.57728) // Center of the map         .transform(         new OpenLayers.Projection("EPSG:4326"), // transform from WGS 1984         new OpenLayers.Projection("EPSG:900913") // to Spherical Mercator Projection         ), 5 // Zoom level         );         timeline = new OpenLayers.Timeline({             map: map,             timeline: "#timeline",             date_key: "when",             //  date_funct:current_date,             format: OpenLayers.Format.GeoJSONTimeline         });         //   var url = "http://toolserver.org/~sonet/api_geojson.php?article=OpenLayers&callback=?";         var url = "example.json";         //  url: "http://localhost:8080/geoserver/world/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=world:cities&maxFeatures=100&outputFormat=json&format_options=callback:processMarkers",         $.getJSON(url, function(data) {             $("#status").text("Data loaded");             timeline.initTimeline(data);         });     }          function toggleCumulative() {         if ($("#cumulative").attr("checked")) {             timeline.cumulative = true;         } else {             timeline.cumulative = false;         }         timeline.update();     }          $(document).ready(function() {         init();         timeline.cumulative = true;         $("#cumulative").attr("checked", "checked");         //  timeline.defaultStyle = new OpenLayers.Style({});     });


title: How to configure spawn-fastcgi Mapserver with Cherokee webserver in ubuntu?

text: Has anyone successfully configured fastcgi mapserver with Cherokee webserver. I'm able to run the mapserver as fastcgi daemon (thanks to https://github.com/dzizes/cgi-mapserver-spawn-fcgi-start-stop- daemon/blob/master/mapserver-fcgi-daemon) but couldn't link it to Cherokee.


title: Convert GeoJSON to TopoJSON

text: # Question How do you convert GeoJSON data to TopoJSON? ## Why Topojson files are often smaller than geojson and the conversion allows you to incrementally simplify a dataset, if you wish.


title: Problem with projection of layers in QGIS

text: I'm having a problem with projections in QGIS that i hope someone here can help me with. I have been searching and found other people having the same problem, but have not been able to solve mine based on those topics. I have a .shp, converted from a .SOSI file, and a .csv file made by me. Both layers are displayed, but they are two different places like shown in this picture: ![Picture](http://i.stack.imgur.com/hW5v4.png) I have also two pictures of the some general info of the .shp file and .csv file. Hope they can give some info. ![layer properties - test](http://i.stack.imgur.com/QaUyJ.png) ![layer properties - 02_admin...](http://i.stack.imgur.com/kS8qn.png) Hope you can help me out here, cant figure this out. On the fly is enabled, i have tryed to choose both same CRS and different CRS, nothing is working. If you need any more info to solve this, just ask me and i will provide it.


title: How to quickly calculate bounding box coordinates of a selected polygon with qgis

text: I just manually (and somewhat inaccurately I'm sure) retrieved the bounding box coordinates for a county polygon, which is a selection from a larger layer of county shapes. I'd like to be able to use a qgis plugin or other function to calculate this quickly and not-by-hand. This is almost certainly possible but I couldn't find anything with a google search thanks


title: Extract raster values at points

text: How can I extract values from a raster by points. I prefer not in Arcgis. I prefer in Qgis or Mapwindow or other open source gis.


title: Load Raster Values into Shapefile

text: I have a shapefile containing many points. Some of these points fall within an area for which I have raster data. I would like to load the values of the raster which underlie each point into the shapefile such that they appear as a column in its DBF. Is there a way to do this with, say, GDAL?


title: What is an efficient Algorithm to draw a straight line from points?

text: I have some cross section survey data, which is like this:               X Y Z pt_dt     2970969.635 359725.0088 83.4242 1-x     2970968.278 359722.2182 83.2591 1-x     2970941.771 359670.127  83.0655 1-x     2970961.369 359708.6424 83.4785 1-x     2970950.888 359688.048  85.2994 1-x     2970955.084 359696.2922 83.6479 1-x     2970953.47  359693.1165 84.9628 1-x     2970948.813 359683.977  83.1451 1-x     2970931.904 359650.7413 83.1042 1-x      When I use the X & Y to make points, I can see the points on the map like this: ![enter image description here](http://i.stack.imgur.com/Opy9N.png) What I need to do now, is to get the points in linear order (i.e. from left to right, or right to left), so that I can create a 3D line, which can serve as the Cross sectional profile. My current algorithm, is to calculate the distance between every pair of points(For one given cross section, which is indicated by pt_dt). The pair which has the maximum distance, are the outermost points. I then arrange the points based on distance from one of the outermost points. I have implemented this algorithm in python, and it works. This algorithm is O(n2); My entire data consists of about 5000 such surveys, so any complexity that I can reduce will be beneficial. So Is there a more efficient way, or algorithm to do this?


title: How we can get output for os.system() in python?

text: I wrote the following command in Geoserver rest-api using curl:               cmd=curl -v -u admin:geoserver -XGET -H "Accept: text/xml" http://localhost:8080/geoserver/rest/workspaces/myworkspace/datastores/mydatastore/featuretypes?list=available         I want to execute this command through python like os.system(cmd).   Now How can I get the list in python itself.   I tried as follow:               import subprocess     import os     from subprocess import Popen,PIPE,STDOUT,call     cmd=r"""curl -v -u admin:geoserver -XGET -H "Accept: text/xml" http://localhost:8080/geoserver/rest/workspaces/myworkspace/datastores/mydatastore/featuretypes?list=available"""       p=sybprocess.Popen(cmd,stdout=subprocess.PIPE,stderr=None,shell=True)        output=p.communicate()       print output


title: QGIS 2.2: prevent labels being cut off / clipped at edge of map

text: I've read a few other issues about trying to prevent labels from rendering if they run off the edge of the map canvas. This is particularly a problem in Composer, where I get lots of country labels which are only half on the page. https://hub.qgis.org/issues/10314 In the above post in the QGIS issue queue, Nyall Dawson says "There's an option in the global label settings to prevent clipped labels at the edge of maps". I can't find this setting. Does anyone know where it is? Or any other solution for this problem?


title: Export shapefile to filegdb with specific projection

text: I tried to export a simple shapefile with polygons via "save as" > Format "ESRI fileGDB" and choose "DHDN / Gauss-Kruger zone 3". So there is the error message : "Error: Failed at creating table for \project (General function failure.)" With choosing the projection "WGS84" it works proper, but the data are needed in "DHDN / Gauss-Kruger zone 3". Any suggestions about that? Or solutions? Thanks!


title: How to edit the position of a banner text's pointer?

text: I have build an ITool where one can add text (`ITextElement`) to the map. The user can choose the symbol from the ESRI symbology. There one can choose a symbol called "Banner Text" which looks like this in the preview: ![Banner Text Symbol Preview](http://i.stack.imgur.com/SbTpk.jpg) When one adds the text with this symbol, only the rectangle is shown without the 'pointer'. How can I edit this pointer? I didn't find any information about it. Maybe I'm simply missing the right name or search phrase for it.


title: write rasters from estudm2spixdf with significant digits

text: I have estimated my kernels using kernelUD from the adehabitatHS package which produces information from multiple animals as a class estUDm. I am trying to write out from R to get into ArcGIS the raster from each of the animals.               myrkde <- kernelUD(mytrial[,1], h="href", same4all=T, grid=200, extent = 0.2)      So, I turn the kernelUD into a spixdf               b <- estUDm2spixdf(myrkde)      Only, I want a raster of each animal so I thought to stack the SPIXDF and then write each out.               rastyrday <- stack(b)     writeRaster(rastyrday, filename=paste("MooseYrHR", names(b)) , format="raster", bylayer=T, overwrite=T)      The problem is when I write out the raster, ArcGIS says Invalid raster dataset. Failed to create raster layer. So I change the format to "img" file...OR I change format as "ascii". Both of these produce something that ArcGIS can read/visualize - **however** , the next problem is that the values are all sooooo small that ArcGIS only registers each value as 0. When I use ArcGIS to look at the pixel values it tells me the value is 0.000000, though the max value is often x*10^-7. Basically, its cutting off after 6 significant digits and my values are 7 to 9. If I writeRaster and format to ascii, the output indicates min and max values are all NA               writeRaster(rastyrday, filename=paste("atesttest", names(b)) , format="ascii", type = "Float32", bylayer=T, overwrite=T)     class       : RasterStack      dimensions  : 200, 138, 27600, 29  (nrow, ncol, ncell, nlayers)     resolution  : 159.6563, 159.6563  (x, y)     extent      : 1470049, 1492082, 6602809, 6634741  (xmin, xmax, ymin, ymax)     coord. ref. : NA      names       : atesttest_Aug2010_F07002_1, atesttest_Aug2010_F07002_2, atesttest_Aug2010_F07002_3, atesttest_Aug2010_F07002_4, atesttest_Aug2010_F07002_5, atesttest_Aug2010_F07002_6, atesttest_Aug2010_F07002_7, atesttest_Aug2010_F07002_8, atesttest_Aug2010_F07002_9, atesttest//_F07002_10, atesttest//_F07002_11, atesttest//_F07002_12, atesttest//_F07002_13, atesttest//_F07002_14, atesttest//_F07002_15, ...      min values  :                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                    NA,                    NA,                    NA,                    NA,                    NA,                    NA, ...      max values  :                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                         NA,                    NA,                    NA,                    NA,                    NA,                    NA,                    NA,      If I writeRaster as "raster" or "HFA" at least the min/max values are indicated - but still I get the failed to create raster when using the raster output type or 0's when converting to the image.               writeRaster(rastyrday, filename=paste("atesttest", names(b)) , format="HFA", type = "Float32", bylayer=T, overwrite=T)     class       : RasterStack      dimensions  : 200, 138, 27600, 29  (nrow, ncol, ncell, nlayers)     resolution  : 159.6563, 159.6563  (x, y)     extent      : 1470049, 1492082, 6602809, 6634741  (xmin, xmax, ymin, ymax)     coord. ref. : NA      names       : atesttest_Aug2010_F07002_1, atesttest_Aug2010_F07002_2, atesttest_Aug2010_F07002_3, atesttest_Aug2010_F07002_4, atesttest_Aug2010_F07002_5, atesttest_Aug2010_F07002_6, atesttest_Aug2010_F07002_7, atesttest_Aug2010_F07002_8, atesttest_Aug2010_F07002_9, atesttest//_F07002_10, atesttest//_F07002_11, atesttest//_F07002_12, atesttest//_F07002_13, atesttest//_F07002_14, atesttest//_F07002_15, ...      min values  :                          0,                          0,                          0,                          0,                          0,                          0,                          0,                          0,                          0,                     0,                     0,                     0,                     0,                     0,                     0, ...      max values  :               2.351000e-07,               4.006146e-07,               6.310832e-07,               3.072779e-07,               2.306967e-07,               2.216463e-07,               2.233004e-07,               2.483713e-07,               4.006146e-07,          2.593051e-07,          1.439714e-07,          1.587329e-07,          2.105937e-07,          2.340155e-07,          2.380073e-07, ...      if I convert the estUDm2spixdf to a raster before stacking -               b <-raster(estUDm2spixdf(myrkde))       It only writes the very first animal - all other animals "layers" are lost. And I still end up with the same problem of the outputs either being invalid or 0's for pixel values. How do I get each animal's UD out of adehabitatHS? Or how do I get the ascii/img files to have significant digits small enough to recognize that there are values in the raster layer? or how do I get the raster file output to be recognized in ArcGIS? I have posted this question to the ade4 listserve, but because writewRaster is not part of their "thing" I thought I would try here too. As some of these commands come from the sp package and the raster package. Im using R 3.0 and ArcGIS 10.2


Using OpenLayers with Bing maps, I'm trying to load JSON data from a service URL and assign an array of clustering strategies to the resulting layer. There are plenty of examples of using cluster strategies with **OpenLayers.Format.GeoJSON** loaded via **OpenLayers.Protocol.HTTP** , but my data is Plain Ol' JSON, not GeoJSON (and I have no control over the format of it). So using this code, the points load on the map correctly. But I get an error: **"Uncaught TypeError: Cannot call method 'read' of null** ". When I open up the OpenLayers source code to the error-causing line, it's this:                layer.protocol.read(OpenLayers.Util.applyDefaults({      and indeed, **protocol** is null. So my question is: How to use OpenLayers cluster strategies without using a Protocol? Here's the code I'm using:               map = new OpenLayers.Map("project-map");     osm = new OpenLayers.Layer.Bing({ name: 'street', type: 'Road', key : bingApiKey});     map.addLayer(osm);     var cp = new OpenLayers.LonLat(0,20).transform(fromProj, toProj);          var pointLayer = new OpenLayers.Layer.Vector("Markers",{                 strategies: [new OpenLayers.Strategy.Fixed(),                  new OpenLayers.Strategy.Cluster({distance: 25})]     });          OpenLayers.Request.GET({         url: "ServiceHandler.json",          success: function(response) {           var jsonReader = new OpenLayers.Format.JSON();           var sparkle = jsonReader.read(response.responseText);             $.each(sparkle, function(i,f){               var point =                   new OpenLayers.Geometry.Point(f.Long, f.Lat).transform(fromProj, toProj);               var pointFeature = new OpenLayers.Feature.Vector(point, null,              {               externalGraphic: "http://localhost:8081/Icons/175.png",               pointRadius: 8             });                 pointFeature.attributes = f;                 pointLayer.addFeatures(pointFeature);             });             }         });              map.addLayer(pointLayer);


title: QGIS 2.2 problem, loading layer into road graph

text: I have the road graph plugin installed, but when I try to configure it with the Settings dialogue, no layer names appear in the Layers field. I have tried various shapefiles, but the Settings dialogue doesn't see any of them. Anyone else seen this ? Thanks. * * * Addition: I am having the same problem. A few more details regarding my experience of the issue, since not too many are provided above: -Road Graph plugin version 0.1 with QGIS 2.2.0 32-bit via OSGeo4W install on a 64-bit Windows 7 machine. Tried same software versions on another Win7-64 machine (thinking maybe some setting was tweaked wrong). No difference. I have tried shapefile and spatialite line layers of topologically sound street networks. Tried having just one layer at a time loaded as well as multiple, with no change. Made sure CRS for layer and project were matched and that 'on the fly projection' was disabled. Also de/re activated plugin. Under Vector>Road Graph>Settings all 3 dropdown menus under the Transportation Layer tab are grayed out/inaccessible with no options shown. All other settings in this dialogue are accessible. Clicking on the canvas to choose Start and Stop points registers coordinates and shows selected points on the canvas. When 'Calculate' is attempted, "Plugin isn't configured" error (not surprising since layer couldn't be chosen).


title: Can someone clear up a lat/long issue? with google maps/arcmap

text: I know a bit about lat/long, projected and geographic coordinate systems etc, and a bit about geodesy and datums, but I would like someone well informed to help me out. I have a lat-long in deg/min/sec provided to me by someone else (point location = northern Italy). I have a map in ArcGIS in the same region and im constructing a grid over it for some analysis, centred on the lat-long point above. In ArcGIS, im working with projection Roma 1940 Gauss Boaga Ovest - the data frame and map layer are both in it - and when adding the XY point also in this projection, it ends up pretty much where it should be (but not exactly) The thing is, when putting these coordinates into Google Earth, the point is around 25 or 30 metres south east of where it is in ArcGIS. The provider is saying the Google Earth location is where it should be. So, where does the difference derive from? I understand Google Earth is a WGS84 geographic coordindate system so i tried a point file in WGS84 and added the XY and it still ends up in the slightly incorrect position in ArcGIS as outlined above. Also have tried various combinations of WGS84 and Roma 1940 Gauss Boaga Ovest and cant get it right. Would the Roma 1940 datum skew it slightly or is it nothing to do with that? Any expert knowledge greatly appreciated. EDIT: I see that google maps actually uses WGS 84 Web Mercator projected coordinate system, which has its own problems. Would ArcGIS actually be more accurate then? * * * i asked the provider and they said the coords provided have come from google maps. I've tried to reproject etc and nothing will change the position from the red dot (see left of image). I've since discovered google maps is a WGS_1984_Web_Mercator projection, based on a sphere. I've attempted to create a custom transformation of datum in ArcMap from 'GCS_WGS_1984_Major_Auxiliary_Sphere' to 'GCS_Roma_1940', no change to dot position when entering the coordinates. As you can see, the google maps on the right is in a slightly different position, by roughly 35m, to the arcmap image on the left. Im pretty much out of ideas. ![enter image description here](http://i.stack.imgur.com/wjX5U.png)


title: Acquiring ArcGIS-like speed in Postgis

text: I have been using Postgis 2.0 for 3/4 of a year now and while I really enjoy using it, excessive query processing time has rendered it basically unusable for my use case. I tend to do heavy geoprocessing on municipal datasets which often have hundreds of thousands of multipolygons. These multipolygons are sometimes shaped very irregularly and can vary from 4 points to 78,000 points per multipolygon. For example, when I intersect a parcel dataset with 329,152 multipolygons with a jurisdiction dataset containing 525 multipolygons, I get the following stats for total time consumed:               ArcGIS 10.0 (on same host with windows 7 OS): 3 minutes     Postgis:56 minutes (not including geometry pre-processing queries)      In other words, **it requires 1500% more time to do this intersection in Postgis than in ArcGIS** --and this is one of my more simple queries! One of the reasons ArcGIS supposedly runs faster is due to better indexes. Some programmers recently figured out how these indexes work and I am wondering if anyone knows how to build these indexes in Postgis (or build tables that would mimic the indexes). Perhaps this would solve most of the speed issues in Postgis. I can only hope there must be some way, especially since ArcGIS can only use 4 GB of RAM while I could use up to 4 times that for my postgis server! Of course there are many reasons postgis can run slowly, so I will provide a detailed version of my system specs:               Machine: Dell XPS 8300      Processor: i7-2600 CPU @ 3.40 GHz 3.40 GHz      Memory: Total Memory 16.0 GB (10.0 GB on virtual machine)          Platform: Ubuntu Server 12.04 Virtual Box VM          Potgres Version: 9.1.4     Postgis Version: POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER      I also detail the entire installation process that I used to set up postgis including the creation of the VM itself. I also increased shared memory from the default 24MB to 6 GB in the conf file and ran the following commands to allow postgres to run:               sudo sysctl -w kernel.shmmax=7516192768 (I know this setting is deleted every time you restart the OS)     sudo /etc/init.d/postgresql restart      As far as I can tell this does absolutely nothing noticeable in terms of performance. Here are links to the data I used for this test:   1. Parcels: tcad_parcels_06142012.shp.zip from City of Austin, TX   2. Jurisdictions: Jurisdictional Boundaries from City of Austin, TX Here are the steps I took to process the data: ArcGIS   1. Add datasets to ArcMap   2. Set coordinate system to central texas feet (srid 2277)   3. Use intersection tool from dropdown menu Postgis Import parcels using:               shp2pgsql -c -s 2277 -D -i -I -W UTF-8 "tcad_parcels_06142012.shp" "public"."tcad_parcels_06142012" |psql -d postgis_testing -U postgres -h local_ip -p 5432      Import jurisdictions using:               shp2pgsql -c -s 2277 -D -i -I -W UTF-8 "jurisdictions.shp" "public"."jurisdictions" |psql -d postgis_testing -U postgres -h local_ip -p 5432      Clean invalid geometry in parcels:               DROP TABLE IF EXISTS valid_parcels;     CREATE TABLE valid_parcels(       gid serial PRIMARY KEY,       orig_gid integer,       geom geometry(multipolygon,2277)     );     CREATE INDEX ON valid_parcels USING gist (geom);     INSERT INTO valid_parcels(orig_gid,geom)       SELECT          gid          orig_gid,         st_multi(st_makevalid(geom))        FROM          tcad_parcels_06142012;     CLUSTER valid_parcels USING valid_parcels_geom_idx;      Clean invalid geometry in jurisdictions:               DROP TABLE IF EXISTS valid_jurisdictions;     CREATE TABLE valid_jurisdictions(       gid serial PRIMARY KEY,       orig_gid integer,       geom geometry(multipolygon,2277)     );     CREATE INDEX ON valid_jurisdictions USING gist (geom);     INSERT INTO valid_jurisdictions(orig_gid,geom)       SELECT          gid          orig_gid,         st_multi(st_makevalid(geom))        FROM          jurisdictions;     CLUSTER valid_jurisdictions USING valid_jurisdictions_geom_idx;      Run cluster:               cluster;      Run vacuum analyze:               vacuum analyze;      Perform intersection on cleaned tables:               CREATE TABLE parcel_jurisdictions(       gid serial primary key,       parcel_gid integer,       jurisdiction_gid integer,       isect_geom geometry(multipolygon,2277)     );     CREATE INDEX ON parcel_jurisdictions using gist (isect_geom);          INSERT INTO parcel_jurisdictions(parcel_gid,jurisdiction_gid,isect_geom)       SELECT         a.orig_gid parcel_gid,         b.orig_gid jurisdiction_gid,         st_multi(st_intersection(a.geom,b.geom))       FROM         valid_parcels a, valid_jurisdictions b       WHERE         st_intersects(a.geom,b.geom);      Explain Analyze intersection query:               Total runtime: 3446860.731 ms             Index Cond: (geom && b.geom)       ->  Index Scan using valid_parcels_geom_idx on valid_parcels a  (cost=0.00..11.66 rows=2 width=1592) (actual time=0.030..4.596 rows=1366 loops=525)       ->  Seq Scan on valid_jurisdictions b  (cost=0.00..113.25 rows=525 width=22621) (actual time=0.009..0.755 rows=525 loops=1)     Nested Loop  (cost=0.00..61428.74 rows=217501 width=24213) (actual time=2.625..3445946.889 rows=329152 loops=1)       Join Filter: _st_intersects(a.geom, b.geom)      From everything I have read, my intersection query is efficient and I have absolutely no idea what I am doing wrong for the query to take 56 minutes on clean geometry!


title: How can I upload more than 10 images in a category map?

text: I did a category map with 30 different categories. For each, I would like to upload an image but the "wizards" permits me to upload only ten images and then one image for "others". Is it possible to upload more? Thank you very much, Simon


title: How to store autocad data in ArcSDE database

text: Hello everybody, Please does anyone have an idea how i can create a database instance on my ArcSDE geodatabase to store my .dwg files. I have over 1 million of these files and would be happy to automatically extract and migrate them into a specified schema in the Geodatabase. I am using ArcSDE 10.1 on SQL Server2008 R2 Thanks


title: ST_Transform not converting from utm to lat/lon

text: This line:               ST_AsText(ST_Transform(the_geom, 4326))      Returns the "coordinates" as utm coordinates. Although, I thought ST_Transform was supposed to convert the coordinates from utm to lat/lon. I found a tutorial where the behavior I need is shown in an example: http://workshops.boundlessgeo.com/postgis-intro/projection.html Is there any reason why it would not be doing this?


title: Why does the file size increase after I clipped the polygons to be smaller?

text: I have a polygon layer that shows properties along a road. I want to create a kml-file of the polygons, and to downsize the file, I tried clipping the properties with a buffer. But, the filesize of the shapefile increased from 370 kB to 1.97 MB, and when I saved the clipped file to kml it ended up being almost 5 MB? Why does the filesize increase and what can I do differently to get a smaller kml-file?


title: Outputting results message from spatial statistics tool to text file using python

text: I'm using the following tool in ArcGIS 10 - Calculate Distance Band from Neighbor Count (Spatial Statistics). I'm using a python script to execute this tool several times on different shapefiles:                   import arcpy         arcpy.env.workspace = "c:/data"         mindist, avgdist, maxdist = arcpy.CalculateDistanceBand_stats("Blocks1", 1, "EUCLIDEAN_DISTANCE")         mindist, avgdist, maxdist = arcpy.CalculateDistanceBand_stats("Blocks2", 1, "EUCLIDEAN_DISTANCE")         mindist, avgdist, maxdist = arcpy.CalculateDistanceBand_stats("Blocks3", 1, "EUCLIDEAN_DISTANCE")      At the moment, this outputs the 3 result values to the Results-Messages window in ArcMap 10, overwriting the past results each time. Is there any way I can write these values to a text file after each execution of the tool using python code. If possible I would like to run this on many shapefiles and write the output values to the same text file. If a text file isn't feasible, then anything that can store this information will suffice. Any help much appreciated!


title: Grouping map points into fixed cluster sizes

text: I have a dataset of 655 lat/long pairs which I'd like to split into about 100 groups. A group should have 5-10 pairs which are geographically close to each other. Dense groups should have more points, sparse groups should have less. For example urban groupings should be larger, rural ones smaller. Is there an established algorithm for doing this kind of grouping, or am I going to have to design one from scratch? I'm using google maps v3 api to display this data, but as it's a fixed dataset I'm prepared to do some offline number crunching.


title: GPS drift during the receiver is stationed

text: I have developed a GPS Tracking device it is working fine when the receiver is in motion, when it stationed in a position for couple of hours i get some drift for around 10-15meters. The receiver send different lat/lon every 30seconds, which makes it impossible for me to calculate the total IDEAL time of the tracker device. Is there any GPS accuracy filter which can be implemented in my software to overcome this issue?


title: Spatial Join in PostGIS according to Line Attribute

text: I am in the process of building a data warehouse in PostGIS, that will store the length of the road network of a country (Lines) per Administrative Area (i.e. Province or State, referred to as AA hereafter and is a Polygon). This is my first time of building a large scale Relational data warehouse so I am not really familiar with this kind of operations. The "id" field of the spatial join resulting table is used to relate to the unique "id" field of the AA polygon layer. The length of the road segments is already calculated and is derived from the "len" attribute of the Road Network table.   So far the operation yields the desired results by performed the following query:               SELECT         m.id,         m.order08,         sum(r.len) as length,       FROM       public.and_lines_2012_03 AS r,        public.and_a8_2010_12 AS m     WHERE       ST_Intersects(m.geom,r.geom)      GROUP BY        m.id, m.order08     ORDER BY        id;      And the resulting table looks something like this:               gid |   id   | length        1  | 000002 | 74118      2  | 000012 | 128131      3  | 000013 | 174296      4  | 000014 | 82240      ...|  ...   |  ...      The Road Network is divided in 3 different line types that are attributed in the Road Network table by the name "line_type" and can have values 0,1,2 My intention is to create a table that will have a column for the length of each line type for each Administrative Area. That is the table would ideally look something like this:               gid |   id   | length_type_00 | length_type_01 | length_type_02      1  | 000002 |     30028      |     11564      |     32526      2  | 000012 |     64270      |     15602      |     48259      3  | 000013 |     104492     |     51981      |     17823      4  | 000014 |     38153      |     9524       |     34563      ...|   ...  |      ...       |      ...       |      ...      I have experimented with different approaches in SQL statements, but so far have not managed to grasp the approach of this kind of queries. The most relative resource I have found on the Internet is this one from OpenGeo: Any help would be highly appreciated Kind regards   Sergios **Update 2012/05/23** @Nicklas I have managed to produce what seems to be very accurate results with the way you have suggested. The query used is the following:               SELECT         b.id AS id,       b.order08 AS order08,       round(sum(b.line_type_00)) AS line_type_00,       round(sum(b.line_type_01)) AS line_type_01,       round(sum(b.line_type_02)) AS line_type_02     FROM     (SELECT         m.id,       m.order08,       (line_type=0)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_00,         (line_type=1)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_01,         (line_type=2)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_02     FROM       public.and_lines_2012_03 AS r,        public.and_a8_2010_12 AS m     WHERE       ST_Intersects(m.geom,r.geom)     GROUP BY        m.id, m.order08,line_type     ORDER BY       m.id     )      AS       b     GROUP BY       b.id,b.order08     ;        The query is quite slow though (4846 tuples take around 41 seconds), since from what I understand it re-projects on the fly to produce the the results in meters. Have I understood/am I doing something wrong? Thank you in advance   Sergios


title: How can I save raster dataset using arcpy

text: I'm trying to save filled raster object to geo database by my specified filename.               import arcpy     from arcpy import env     from arcpy.sa import *          env.workspace = r"C:\GIS Data\Cameron_Run\Topo_NED\Projected"          # Set local variables     surfaceRaster = "raw_10m_dem"          # Check out the ArcGIS Spatial Analyst extension license     arcpy.CheckOutExtension("Spatial")          # Execute FlowDirection     fillRaster = Fill(surfaceRaster)     arcpy.RasterToGeodatabase_conversion(fillRaster, r"C:\GIS Data\Cameron_Run\Topo_NED     \Projected\fillRasterGDB.gdb")      But this code is saving raster file by its own specified name, for this case `Fill_raw_10m1` is the default name of saved raster file. Is there any way so that I can specify the name of saved raster file?


title: gdalwarp cutline along with shapefile

text: I am trying to clip my Raster file according to a shape file geometry. I am using the following code               gdalwarp -cutline INPUT.shp INPUT.tif OUTPUT.tif      But it is resulting a black color outside the shape file geometry extent. I am giving the example here. In the first Image that I want to clip. The second image is the resultant raster but I want the 3rd Raster where the outside of the geometry will be null. ![1st Image: INPUT Raster](http://i.stack.imgur.com/osaj8.jpg) 1st Image: INPUT Raster ![2nd Image: OUTPUT Raster](http://i.stack.imgur.com/cKs1v.jpg) 2nd Image: OUTPUT Raster ![enter image description here](http://i.stack.imgur.com/34Jh6.png) 3rd Image: I want to get output like this