title: Sum values in a field

text: Using QGIS 1.8, I've got a shapefile with an integer field called AREA and I'd like to sum the values which are in this shapefile. Does anybody know how to ? Thank's a lot. Robert


title: Leaflet, wrong API key when using CloudMade tiles!

text: I've registered on CloudMade and got some API keys but non of them work with Leaflet API. When I checked for image URLs and opened each directly using my browser it prints out the message "Wrong apikey". Anyone with the same issue?


title: How does Mobile GPS work in Windows Tablets?

text: Can anyone tell me how GPS chips generally works in Windows8 Tablets? As far as I can tell Tablets do not have a separate GPS chip like those in handheld GPS receivers - It seems to be combined into the Network chips. With my current tablet, gps doesn´t start without first logging into the network with a mobile SIM card. They seem to first get a fix according to the location of the last successful IP Address login and then refine the position using triangulation from satelites. I would be grateful if anyone with experience with such devices could clear this up for me.


title: Is there a source for free CSV files with geo data and other desired fields?

text: I'm writing some software components that process geospatial data. The input data is a CSV with general columns (e.g. `name`) and geo columns (e.g. `city`, `lat` or `zip`). I would like to test these components with some thousands of points. **Is there an online source that generates free csv files with geo data, according to a list of desired fields?** For example, I would like to have 1,000 points in North America, each with `lat`, `lon`, `city`, `zip` and `state`. Or 5,000 points in Europe, each with `lat` and `lon`. I will probably build it myself from free data sources if it doesn't exist, but using an existing source would save lots of time.


title: Distance-direction computation of lines with exactly/more than two nodes?

text: A distance and direction tutorial provided as an Answer at **How to add Direction and Distance to attribute table?** uses:   1. _xat (0 ), xat (-1 ) and yat (0 ), yat (-1)_ for **start and end points of lines** ; in addition to    2. (*atan((xat(-1)-xat(0))/(yat(-1)-yat(0)))) * 180/3.14159 + (180 *(((yat(-1)-yat(0)) < 0) + (((xat(-1)-xat(0)) < 0 AND (yat(-1) - yat(0)) >0) _2)))_ to **calculate direction** I would like to know any difference in application if lines with more than 1 part and/or nodes are used. Would there be any difference between using lines with 2 nodes vs lines with more than 2 nodes and segments? Which assumptions would QGIS automatically use/make in the following scenarios below: a - regarding **direction** : the overall direction is derived from the difference in start and end node **OR** would it be the average of each individual computed segment. b - regarding **distance** : is it the distance between start node and end node **OR** is it the sum total of each segment that is represented in the distance


title: QGIS "transform error" messages - how to stop them?

text: In QGIS, I am getting an error message every time I edit a node. Also, there is a lot of lag in performance. The message is: "Transform error caught: forward transform of (10990.308086, 100347.255186) PROJ.4: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +to +proj=utm +zone=54 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs Error: latitude or longitude exceeded limits" I hae been very careful to Set Layer CRS on all layers, to GDA94, or, in the case of the two raster layers, GDA94/MGA zone 54 (which is also the project CRS). This should all work. Shouldn't it? But it is odd, that when I was using my legacy MapInfo TAB layers, they would show up in QGIS with weird projections that were not the projection in the TAB files, like some old Dutch projection and, more usually, NAD projections. I saved those layers as SHP files (resetting the projection to GDA94) to see if QGIS works better with SHP. And it did, for a few hours. Now the error messages are back. What is wrong? Is this a QGIS bug? Or is there a backoffice setting that I've missed? * * * ## EDIT: ADDITIONAL INFORMATION. Below I provide the sequence of steps to create the transform error from a clean start.   1. Open QGIS   2. Set these options for Project CRS:      * "Automatically enable OTF if layers have different CRS"     * Do not "enable OTF by default"     * Start new projects with GDA94     * "When new layer is created, use project CRS."   3. Add raster layer 76949_zone54_mga.ecw (www.ga.gov.au/cedda/maps/1058)   4. Look at Project properties, witness that "OTF by default" is disabled, and CRS has changed to greyed-out MGA94 zone 54 (the CRS of the raster file)   5. Add shapefile layer (graticule from "Hamilton Special" GA 250K topo), downloaded from GA website.      * Note that its .prj file contents are: GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]     * That is, it has a CRS of GDA94. The file loads and displays correctly. Look at project properties, and witness that it has automatically enabled OTF projection.   6. Create new shapefile point layer, with CRS selected as GDA94   7. Edit new shp layer, by making two points.    8. Exit edit mode and save.   9. Look at message log: two error messages logged with following info: Transform error caught: forward transform of (0.000000, 0.000000) PROJ.4: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +to +proj=utm +zone=54 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs Error: latitude or longitude exceeded limits That is a replicable error, and all files used are pristine Geoscience Australia product, unquestionably ok. The fault lies somewhere in QGIS or in the settings used. Yes I could create the new file in a CRS of MGA_zone54 (to follow the raster projection in this case), and that will work ok. But I would prefer that my vector files that aren't tied to a particular UTM zone. Anyway, GDA and MGA zone 54 are based on the same datum and ellipsoid, they are two sides of the same coin. MGA is just the projection of GDA94 longlat data. So what is the problem?


title: Write raster with its attribute table and symbology

text: I am classifying lots of rasters and I'd like to automate the legend creation and symbology displaying for easy manipulation. Here is the kind of raster I'm dealing with:               library(raster)     r <- raster(nrow=10, ncol=10)     r[] = 1     r[51:100] = 3     r[3:6, 1:5] = 5     r <- ratify(r)          rat <- levels(r)[[1]]     rat$landcover <- c('Pine', 'Oak', 'Meadow')     rat$class <- c('A1', 'B2', 'C3')     levels(r) <- rat     r          levelplot(r, col.regions=c('palegreen', 'midnightblue', 'indianred1'))


title: Split an attribute Python 9.3 in modelbuilder expression?

text: **Updated question** Goal: I need to upgrade an existing Model-builder were the expressions uses Python 9.3. Workproces & Goal: My *.mxd Have 1 attribute field containing the following build: X1-X2-X3-x4_x5_x6-x7.   Example: G-OI-KI-VUUR_NAAMGEVING_0000101111_VECT_12X001111-X. I tryed to split is on the ‘-’ and ‘_’ signs because the names have a different length. My goal is putting the fields X1 - X7 in separate text-fields by using a Model-builder which is coded with Python 9.3. **Issue** : one of those 'transformations' (X4) isn’t working properly. I used _| **!LAYER!.split("-")[3]** |_ A new field is created, but I got as a result: x4_x5_x6-x7. I only wanted x4. Matthieu Brbt came with a solution to use | ***(!LAYER!.split("-")[-1]).split("_")[0]*** | instead. Unfortunately this ended up with an error: 000539: Indentation Error: unexpected indent (, line 1) Failed to execute (Calculate Field) Matthieu suggested to check for empty spaces in the expression but those were not here. This is what my expression looks like: ![enter image description here][1] **So how can I make a split based on the ‘-’ and ‘_’ for field x4 ?** Thank you for your help!


title: How to better visualize features of clusters in Openlayers

text: I am using the cluster strategy for mapping a large number of features in Openlayers, now I would like to find a good solution for providing info about the features contained in a cluster selected by the user. In other words whenerver the user selects one cluster I would like to explode the cluster with as many markers as the features contained in it, of course selectable markers. Anybody could help me? I attached an image about how the markers around the selected cluster would look like. thanks in advance. Defelix ![enter image description here](http://i.stack.imgur.com/sFnoc.png)


title: Requirements for BaseCustomLayer in ArcScene?

text: What must be implemented/included to add a layer deriving from BaseCustomLayer to ArcScene? I have had success adding a BaseCustomLayer to ArcMap by implementing Draw, AreaOfInterest, Cached, MaximumScale, MinimumScale, SpatialReference, Visible, and Name. However, when I attempt to add this layer to ArcScene it gives me the error "ArcScene cannot add this layer".


title: How do I add a KML layer using the ArcGIS Silverlight API?

text: I am trying to create a map which will include radar rainfall data in a KML format. How can I easily do this using the Silverlight API (the version doesn't matter)?


title: What is the maximum number of features a SDE featureclass can hold?

text: I have a SDE 10.0 Geodatabase based on MS SQL 2008. I'm running a continuous process where about 6400 point features are written to a feature-class every hour. I need it to be a featureclass, since I'm doing some krigging on the feature-class, after applying an attribute query based on a time field in the data. I was wondering what is the limitation on the number of features a SDE Feature-class can hold?


title: Convert point to grid: Numbers of cells differ

text: In R I am trying to aggregate points data to a grid. The data for the grid (G-Econ) comes as a csv-file. In this file each row represent a grid-cell at 1x1 degree resolution and the given coordinates, in longitude and latitude, label the south-west corner of the grid-cell. In trying to convert these points to a grid I run into the issue that the number of grid-cells in the output (so the eventual grid-raster) does not correspond with the number of cells in the input. Basically the short version of my question is: How do I solve this? The simplified G-Econ csv-file can be found here and below I describe step-by- step what I've been doing. First I load the csv-file, limit the number of observations to only include cells in African countries, and create a spatial object out of the data frame.               # Load data     gecon<-read.csv("gecon.csv",header=TRUE,sep=",",row.names=NULL)     gecon$Cell<-1:nrow(gecon) # Create a grid cell identifier           # Limit data to African countries     require(countrycode)     gecon$cont<-countrycode(gecon$NEWCOUNTRYID,"iso3n","continent",warn=TRUE)     g.africa<-na.omit(gecon[gecon$cont=="Africa",])           # Create spatial object     require(raster)     coordinates(g.africa)<-~LONGITUDE+LAT # Creates spatial object     proj4string(g.africa)=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")          summary(g.africa)      From this we see that Africa is covered by 3550 grid cells in the dataset   Next step it to transform the points into a grid:               # Transform points to grid     extent<-extent(g.africa)     extent # Determine number of cols. and rows for grid     length(-26:63) # Number of cols: 90     length(-38:37) # Number of rows: 76     rast<-raster(extent,ncol=90,nrow=76,              crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")     grid<-rasterize(g.africa,rast,g.africa$Cell,fun="last") # Create grid     grid      During this step the number of cells has increased -to 6840- this is due to the inclusion of cells that cover the oceans. I included a cell identifier so I can check how many of the new cells do not belong to a country.               summary(grid)       There are 4053 NAs. Subtract these from the total number of grid-cells and the result is 2787 grid-cells, which is way less than the original number of cells for Africa. I transform the grid into polygons and this number is confirmed.               africa<-rasterToPolygons(grid) # Create spatial polygons     africa # See nfeatures       So I guess something goes wrong when I rasterize the data, that somehow the interpolation of points to grid-cells results in this reduction. Due to my inexperience with GIS in R -and in general- I have been unable to solve this issue. Does anyone have any suggestions on how I could solve this? Or maybe recommendations for a better method?


title: How to install python-mapscript in a virtualenv?

text: I am on a django project where I want to use the python implementation of the mapserver mapscript library. I was able to install mapserver and python- mapscript from source. However, due to the project setup I am forced to run django in a virtualenv that does not allow site-packages in my virtual environment. Since mapscript is installed as a global site-package my django app can't see the mapscript installation. When I try to install mapscript through `pip install mapscript` I get an error `'mapserver.h' file not found`. So how can I install python-mapscript in my virtual environment? I thought of copying the installed egg into my virtual environment, but not sure if that is possible.


title: How is it possible to publish data with ArcGIS Server without creating ArcSDE service or using ArcSDE?

text: I am still relatively new to the ESRI server techologies. I have server1 which runs ArcGIS server 10.0 and when I look at which processes running on this server I see ArcSOM.exe, ArcSOCMON.exe and many ArcSOC.exe (probably around 50 instances) I have another server (server2) which runs 6 instances of giomgr.exe 11 instances of gsrvr.exe and 2 Imgrd.exe instances. There are also 2 DBMS Servers which hold sde-geodatabase data. My question is how have I been able to publish data from a file geodatabase stored on a partition on server2 although there is no SOC, SOM, DBMS or AGS running from it?


title: Content for Introduction to GIS course using QGIS?

text: For teaching an Intro to GIS class using QGIS, what could be its content? The class will be for an Masters in Public Administration program, 20 hours over the course of a weekend, with the students often looking for jobs as government analysts and managers or nonprofit executive directors. I want them to be able to make a simple map that can go in a document and be able to talk to technical GIS staff intelligently. It will probably be the only course in GIS they ever take and they probably won't be looking for work as GIS techs. I will do a couple of hours demonstration on QGIS to ESRI parallels during the weekend, but teach the course with QGIS. QGIS pluses: it is a fine platform for teaching the concepts, and it is free. Granted ESRI will be more recognizable on a resume, but if they get a job at a place with no GIS, it will be easier to get a budget for QGIS. I have also decided using QGIS somewhat selfishly -- I run FreeBSD at home, and I want to be able to prep on my machine. My basic text will actually be "Mapping It Out" (Monmonnier, 1996 ??), which they will read before class and write a short page on, since it is most important they learn how to make an argument from a map as part of a bigger decision making process. I will assign the QGIS tutorial to be skimmed, and have as optional Glaser's book on QGIS and old editions of Dent's Cartography. I think I will do a demonstration with questions (~2 hours), a little bit of lecture on basic concepts (~3 hours, distributed over the weekend), a lot of practice time (~12 hours), and then student presentations and discussion (~3 hours). The demo and student maps will start with some addresses, make a locator map of the points, edit them a little for practice, do something with the field calculator, figure out the local demographic characteristics of the points by which tract they are in, calculate some point-in-polygon rates based on local population and point/ event counts, and finish by embedding the PNG in a document and writing a paragraph that refers to the map.


title: Why does AddGeometryColumn error using the libspatialite-2 DLL with pysqlite2?

text: What I'm trying to achieve is to create a geometry column of a point based on the latitude and longitude of that particular row in the already created table. Given the following code:               #!/usr/local/bin/python     from pysqlite2 import dbapi2 as sqlite3          conn = sqlite3.connect('adde.sqlite')     conn.enable_load_extension(True)     conn.execute("SELECT load_extension('libspatialite-2.dll');")     conn.execute("SELECT InitSpatialMetadata();")     conn.enable_load_extension(False)     conn.execute("SELECT AddGeometryColumn('coords','geom_column', 4326, 'Point', 2);")     conn.execute("UPDATE coords SET geom_column = GEOMETRYFROMTEXT('POINT('||longitude||' '||latitude||')', 4326);")     conn.close()      I receive the error:               AddGeometryColumn() error: "insert on table 'geometry_columns' violates constraint: 'spatial_ref_sys.srid'"      I've been able to create the column in the Spatialite GUI no problem, but I'd like to do all of this programmatically using Python. I've searched for a while now and I'm coming up a bit short. Any guidance as to why I'm receiving this particular error would be greatly appreciated. I'm using Python 2.6 and pysqlite2 on Windows 7 64-bit


title: How can I create a pollutant dispersion raster in QGIS?

text: As my previous questions have indicated, I am working on an air quality problem where I have a dataset of industrial emissions that consists of a location (long, lat), and an annual volume of emissions in pounds or tonnes. The end result that I am looking for is a 'dispersion raster' that will provide an aggregate value of emission dispersion (see below) based on overlapping 11KM buffers, where the dispersion value of each buffer is aggregated wherever the buffers (2, up to 57!) overlap. **Details** I am trying to understand how I can model the dispersion of emissions, from the emissions source, through a simple linear decay model where the dispersion at a point (away from the emissions source) is ((1-d)/(11000m)) (from Hystad 2011 - "Environmental Health Perspectives") where d = the distance to the emissions source in meters. I envision creating a buffer around each emissions point (and have done so, see below) but do not know how to apply the above formula to the buffer. Once the buffer has that formula, thus providing a range of values that degrade to zero as you move from the emissions source, I would like to aggregate all the values (through the vertical or Z) of overlapping buffers. For instance, the attached image provides a view of the emissions sources (red dots) and the buffers (the grey buffers are 'low values' and the red 'high'). How would I create buffers that decay from the emissions point by ((1-d)/(11000m))? Once these are created, how could I calculate the aggregate (i.e. the exposure from all the buffers that overlap a point) exposure value at a point, as represented by the blue dot. This is hard to explain, so please ask questions! Thanks! ![Red= Emissions Sources. Buffers = Range of dispersion. Blue = How do I calculate an aggregate dispersion volume for here?](http://i.stack.imgur.com/DsqXq.png)


title: Using polygons to apply unique value to all raster cells within those polygons?

text: In ArcGIS 9.3, I have a raster file that I would like to "mask" with polygons, giving a single value to all the raster cells that overlap the polygons. How can I do this? Here's a screen shot, where the black polygons are my polygon layer, and the rest is raster-based. ![enter image description here](http://i.stack.imgur.com/in1z9.jpg)


title: Cluster processing - what resources to ask for 50 Gb DEM processing

text: For "renting purposes", I have to specify the computing resources I need from the computing center at my University. I have very little idea of what to ask for. Considering the information below may you suggest something? This is the computational cluster: "A ten blade server with each blade consisting of 2 Quad Core Intel Xeon 2.33 Ghz CPUs with 16 GB of memory for a total of 80 CPUs and 160 GB of memory. Accounts on a 2 TB disk system. Running 64-bit Linux" With lossless compression, the original data is one 50gb tiff file. At some point in time i shall be working with several (up to ten) files 50gb in size. I would use GDAL, Python scripting and maybe C++ scripting. If allowed, I would use softwares like Grass GIS and Saga GIS. Also, R language with spatial libraries. I will be deriving the usual terrain parameters, trying to apply object-oriented algorithms for the extraction of specific features (landforms) and using statistical analysis for descriptive and modeling purposes. Thank you in advance, Marco


title: is it possible to obtain layerInfo from the Layer's URL

text: I am trying to fetch the names of the layers in the MapService, so that I can populate them in the dijit.form.select. I am able to get the `layerInfo` when the ArcGISDynamicMapServiceLayer is added to the map by using,               map.getLayer(map.layerIds[0]).layerInfos[0].name      But this is only after the layer has been added to the map. Is there some way by which I can get the layerInfo names directly from the URL?


title: Stitching thousand of 32 bit Tiff without geographic information

text: I'm working on a project which goes beyond my knowledges, I need your help. I'm working with FLIR movies (thermal imagery recorded from plane), and I have to find a way to make a Tiff containing the temperature information and geographic information. I can extract images from the movies in different ways (png, jpeg, csv, Tiff, 32 bit tiff, FITS). My plan was to extract one jpeg/second and then stitch those thousand of pictures with hugin, then georeference the stitched result on ArcGis. My problem is that I suppose the JPEG will be in RGB format so the temperature won't be displayed... So the only way might be to export 1 32 bit tiff/second (which is the only format with the real temperature information) and then try to stitch them together... But I don't see how. The study area is really huge so I have 5000 thousand pictures each day, so I can't afford any manual method. Does someone has a genuis idea?


title: Mapping and Pedestrian Routing on Android (guidance needed)

text: I 'm new to this stuff (GIS). I was assigned a project in which i have to implement for the Android platform:   1. showing up maps   2. showing up routes   3. showing up nearby POIs. Everything has to be offline!!(loaded on the Android device!) with NO connection to Wifi or GSM network. Please correct me if i'm wrong (as i stated above this is brand new stuff for me - i never did a project in this field, though 20 years in business): Can the above points be solved with the following tools (respectively)?   1. using OSMDroid (with ZIPped Openstreet maps)   2. using Spatialite (just quickly read through some articles but not sure if i can use it!) - or should i use pgRouting (i understand i have to port it to SQLite - but is it possible?))   3. using "simple maths" (great-circle distance,...) The above points have to be implemented (as far as the amount of data concerned separately) for maps like Italy, Spain, France, Germany, and Greece. I need some guidance!! Is there anyone who has done something similar? Thanks in advance!


title: Shape file polygons to STL format conversion

text: I have Building polygons like Ground, wall and Roof polygons from this need to generate STL format( 3d Printer format). Please can anyone suggest how to do this?


Hy, I am running the following code which I found on: GDAL - Perform Simple Least Cost Path Analysis which I adapted to my case, but I get the following error when running the code in the python console of QGIS. I am rather new to the programming & QGIS world so if somebody would be able to help me that would be wonderful!               import gdal, osr     import osgeo.gdal     osgeo.gdal.GetDriverByName     from skimage.graph import route_through_array     import numpy as np          def raster2array(rasterfn):         raster = gdal.Open(rasterfn)         band = raster.GetRasterBand(1)         array = band.ReadAsArray()         return array            def coord2pixelOffset(rasterfn,x,y):         raster = gdal.Open(rasterfn)         geotransform = raster.GetGeoTransform()         originX = geotransform[0] # East/West location of Upper Left corner         originY = geotransform[3] # North/South location of Upper Left corner         pixelWidth = geotransform[1] # X pixel size         pixelHeight = geotransform[5] # Y pixel size         xOffset = int((x - originX)/pixelWidth)         yOffset = int((y - originY)/pixelHeight)         return xOffset,yOffset          def createPath(CostSurfacefn,costSurfaceArray,startCoord,stopCoord):                 # coordinates to array index         startCoordX = startCoord[0]         startCoordY = startCoord[1]         startIndexX,startIndexY = coord2pixelOffset(CostSurfacefn,startCoordX,startCoordY)              stopCoordX = stopCoord[0]         stopCoordY = stopCoord[1]         stopIndexX,stopIndexY = coord2pixelOffset(CostSurfacefn,stopCoordX,stopCoordY)              # create path         indices, weight = route_through_array(costSurfaceArray, (startIndexY,startIndexX), (stopIndexY,stopIndexX),geometric=True,fully_connected=True)         indices = np.array(indices).T         path = np.zeros_like(costSurfaceArray)         path[indices[0], indices[1]] = 1         return path          def array2raster(newRasterfn,rasterfn,array):         raster = gdal.Open(rasterfn)         geotransform = raster.GetGeoTransform()         originX = geotransform[0] # East/West location of Upper Left corner         originY = geotransform[3] # North/South location of Upper Left corner         pixelWidth = geotransform[1] # X pixel size         pixelHeight = geotransform[5] # Y pixel size         cols = array.shape[1]         rows = array.shape[0]              driver = gdal.GetDriverByName('GTiff')         outRaster = driver.Create(newRasterfn, cols, rows, gdal.GDT_Byte)         outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))         outband = outRaster.GetRasterBand(1)         outband.WriteArray(array)         outRasterSRS = osr.SpatialReference()         outRasterSRS.ImportFromWkt(raster.GetProjectionRef())         outRaster.SetProjection(outRasterSRS.ExportToWkt())         outband.FlushCache()              def main(CostSurfacefn,outputPathfn,startCoord,stopCoord):              costSurfaceArray = raster2array(CostSurfacefn) # creates array from cost surface raster              pathArray = createPath(CostSurfacefn,costSurfaceArray,startCoord,stopCoord) # creates path array              array2raster(outputPathfn,CostSurfacefn,pathArray) # converts path array to raster               if __name__ != "__main__":         CostSurfacefn = "Administratief_perseel_1polygoon_PXLKOST_6000_6000.tif"         startCoord = (174084.004,176969.786)         stopCoord = (173697.916,175206.172)         outputPathfn = 'Path.tif'         main(CostSurfacefn,outputPathfn,startCoord,stopCoord)      Error Code:               Traceback (most recent call last):       File "<input>", line 1, in <module>       File "/Users/PeterVanvoorden/Documents/GroepT/Thesis/c)QGIS/CLIP_AGIV_GRBgis_e40179/Projects/python_console_scripts/least_cost_path0.1.py", line 77, in <module>         main(CostSurfacefn,outputPathfn,startCoord,stopCoord)       File "/Users/PeterVanvoorden/Documents/GroepT/Thesis/c)QGIS/CLIP_AGIV_GRBgis_e40179/Projects/python_console_scripts/least_cost_path0.1.py", line 65, in main         costSurfaceArray = raster2array(CostSurfacefn) # creates array from cost surface raster       File "/Users/PeterVanvoorden/Documents/GroepT/Thesis/c)QGIS/CLIP_AGIV_GRBgis_e40179/Projects/python_console_scripts/least_cost_path0.1.py", line 9, in raster2array         raster = gdal.Open(rasterfn)     AttributeError: 'module' object has no attribute 'Open'      I assume there is a problem with opening the file that I want to analyze "Administratief_perseel_1polygoon_PXLKOST_6000_6000.tif" but I just don't know how to resolve this.. Thanks in advance!


title: How to remove <VALUE> text from legend in layout view?

text: I'm producing a map in the layout view which I will eventually export as a jpeg. I've added a legend which includes 15 classes but when I insert a legend in layout view its including the word VALUE at the top of the legend column - I had a look at the properties and can't seem to figure out how to remove this - I'm using ArcMap 10.


title: Is it possible to use RPy/RPy2 in a Geoprocessing Script that is exposed as a GP Task on ArcGIS Server?

text: Is it possible to publish a Script tool that uses the RPy/RPy2 module (or more generally any Python-to-R link) to ArcGIS Server 10 SP2? If so, what specifically do I need to do, i.e. Install R on the server machine(s). Thanks.


title: Use Ushahidi api in another map using openlayers

text: I have a ushahidi deployment and i am using Openlayers + ushahidi api to show the reports in a separate map along with other features. i have queried the data from ushahidi api and made a vector layer. the layer has data but it not showing up. This part gets the data from ushahidi               var roadreportrequest = OpenLayers.Request.GET({     url: '../roadreport/api/?task=incidents',     async: false,     success: function(Request) {         // alert("success");         roadreportjson = JSON.parse(Request.responseText);     },     failure: function() {         alert("failed to load data from roadrepoter");         roadreportjson = {}     }      }); and this part converts it to vector layer               var roadreport_vector = new OpenLayers.Layer.Vector('RoadReporter Reports');          for (incident_index in roadreportjson.payload.incidents) {     // console.log(incident_index, roadreportjson.payload.incidents[incident_index]);     incident = roadreportjson.payload.incidents[incident_index];     roadreport_vector.addFeatures([         new OpenLayers.Feature.Vector(             new OpenLayers.Geometry.Point(incident.incident.locationlongitude, incident.incident.locationlatitude),             incident)     ])     }     map.addLayer(roadreport_vector);      There is data in the layer, i have checked the console and this code to pan the map also works.               map.panTo(new OpenLayers.LonLat(roadreport_vector.features[2].geometry.x, roadreport_vector.features[2].geometry.y))      My issue is that i am expecting the data to be shown up as on the map(the default orange circles), but it is not showing up.


title: PostGIS: Insert a Point using Python

text: What is the right way to insert a Point to a PostGIS database using Python?


title: Efficiency of python script using SAGA and arcpy

text: === SUMMARY === I have created a process to import a LAS point cloud, clip it, select the ground points and save the output as an XYZ.txt for clients to use. This uses a python based script (will become a tool) that utilizes some of the better features in SAGA GIS as well as arcpy. SAGA seems to first create/save the spc (saga point cloud) and then re-import it before creating a shape file. Can I get it to work only in the memory? Also once the spc is cut by the AOI, I want it to split the file into 200,000 record segments (max allowable by some clients sofware) -can this be done? === UPDATE 2 === based on @johanvdw posting of Wiki with answers == I'm still having issues with what the sys arg values are and how to set the input/output. I now have the following test code...I have added questions to it.               import saga_api, sys, os, arcpy          ##########################################     def las2spc(fLas):     # fLas: python string containing the filename     # returns a pointer to a pointcloud         print saga_api.SAGA_API_Get_Version()              Library = saga_api.CSG_Module_Library()         if not Library.Create(saga_api.CSG_String('C:\SAGA\modules\io_shapes_las.dll')):             print 'unable to load SAGA module library'             return 0              Module = Library.Get_Module('Import LAS Files')         print "Module Loaded"         pSpc  = saga_api.SG_Create_PointCloud()         Parms = Module.Get_Parameters() # default parameter list         Parms('FILE').Set_Value(saga_api.CSG_String(fLas)) #filename in saga string format         Parms('POINTS').Set_Value(pSpc) #pointcloud data object as required by the module (see above)         Parms('T').Set_Value(0) #This are booleans, either on or off         Parms('i').Set_Value(0) #T eg is the time, see above         Parms('a').Set_Value(0)         Parms('r').Set_Value(0)         Parms('c').Set_Value(0)         Parms('u').Set_Value(0)         Parms('n').Set_Value(0)         Parms('R').Set_Value(0)         Parms('G').Set_Value(0)         Parms('B').Set_Value(0)         Parms('e').Set_Value(0)         Parms('d').Set_Value(0)         Parms('p').Set_Value(0)         Parms('C').Set_Value(0)         Parms('VALID').Set_Value(0)         Parms('RGB_RANGE').Set_Value(0)              if not Module.Execute():             print 'module execution failed'             return 0         print 'module successfully executed'              return pSpc                if __name__ == '__main__':      # I have commented the following out as it wasn't clear what it does. Are the sys.arg values from running this as a tool?               ##    if len( sys.argv ) != 2:     ##        print 'Usage: las2shape.py <in: shape> <out: pointcloud>'     ##    else:             fLas   = '335000_8084000_1K_class.las'       # I put in the name of a las file in the same directory as the py file.**                       if os.path.split(fLas)[0] == '':                 fLas   = './' + fLas             pSpc = las2spc(fLas)             fSpc   = 'file.spc'       # this is the name I want to give the output                       if os.path.split(fLas)[0] == '':                 fSpc   = './' + fSpc                 print "Saving"             pSpc.Save(saga_api.CSG_String(fSpc))             print "Saved"     #or instead of saving: use pSpc as on of the inputs to your next module; the wiki gives an example.      which gives the following result. > SAGA Application Programming Interface - Version: 2.0.7 Module Loaded module > execution failed Traceback (most recent call last): File > "P:\2011\Job_193_LAS_2_SHP_Procedure\Working\NewScript\SAGA.py", line 57, in > pSpc.Save(saga_api.CSG_String(fSpc)) AttributeError: 'int' object has no > attribute 'Save' === UPDATE 1 === based on @johanvdw posting of Wiki with answers == Based on the input that we can't save in_memory using the cmd line, I have switched to python. I am still having issues with figuring out what some of the code means.   1. What do all the sys.arg[1-4] values mean? Is this designed to run as an arcgis tool?    2. "OUTPUT .Save(saga_api.CSG_String(sys.argv2))" is the last line -  When I try the code below it gives me an error -               # Python script template for SAGA module execution (automatically created, experimental)          import saga_api, sys, os, csv          def SAGA_importLas(SHAPES, ZFIELD, OUTPUT):         print saga_api.SAGA_API_Get_Version()              Library = saga_api.CSG_Module_Library()         if Library.Create(saga_api.CSG_String('C:\SAGA_win32\modules\pointcloud_tools.dll')) == 0:             print 'unable to load SAGA module library'             return 0              Module = Library.Get_Module('Point Cloud from Shapes')         #print "Opened module"         Parms = Module.Get_Parameters() # default parameter list         Parms('SHAPES').Set_Value(0) # input NOT optional shapes         Parms('ZFIELD').Set_Value(2)         Parms('OUTPUT').Set_Value(1)              print SHAPES, ZFIELD, OUTPUT         if Module.Execute() == 0:             print 'module execution failed'             print sys.stdout             return 0         print 'module successfully executed'         return 1          import os, arcpy, sys, datetime          top = os.getcwd() # change to a specific path if required.     # This otherwise starts with the directory the script is in.     RootOutput = top # change if you want output somewhere else     FileTypes=['las'] # add filetypes or * as required.     SearchStrings=['Temp_Pol', 'Spatial_Ex'] # add strings as required.     RootDirectory = r'L:\Vector_Data\LAS'          filecount=0     #successcount=0     #errorcount=0          print "Working in: "+os.getcwd()          List =[]     f = open(RootOutput+'\\Success_LOG.txt', 'a')     f.write("Log of files Succesfully processed. RESULT of process run @:"+\             str(datetime.datetime.now())+"\n")     f.close()          FileList = csv.reader(open('Request.csv'))     SearchStrings=[]     for File in FileList:         #SearchStrings.append(File[0]+","+File[1])         SearchStrings.append(File[0])          for root, dirs, files in os.walk(RootDirectory, topdown=False):         for fl in files:           currentFile=os.path.join(root, fl)           for FileType in FileTypes:               status= str.endswith(currentFile,FileType)               if str(status) == 'True':                   for SearchString in SearchStrings:                       if str(SearchString in currentFile) == 'True':                         #print str(currentFile)+str(status)                                List.append(currentFile)           filecount=filecount+1          #print 'File List: ' + str(List) + 'containing ' + str(len(List)) + ' files.'          #for item in List:     Item=List[45]     print "Working on : " + Item     if __name__ == '__main__':     ##    if len( sys.argv ) != 4:     ##        print 'Usage: this_script.py <in: gridfile> <out: gridfile> <in: shapefile> <out: shapefile>'     ##     ##    else:             SHAPES = Item             ZFIELD = ''             PositionDot = Item.find('_1K_class')             PositionDir = Item.rfind('\\')             OUTPUT = Item[(PositionDir+1):PositionDot]+'.spc'                       if SAGA_importLas(SHAPES, ZFIELD, OUTPUT) != 0:                 OUTPUT  .Save(OUTPUT)      Gives result > Working in: P:\2011\Job_193_LAS_2_SHP_Procedure\Working\scripts Working on : > L:\Vector_Data\LAS\Rangeview_Tolga\338000_8095000_1K_class.las SAGA > Application Programming Interface - Version: 2.0.7 > L:\Vector_Data\LAS\Rangeview_Tolga\338000_8095000_1K_class.las > 338000_8095000.spc module execution failed What's happening? I also don't understand what if "len( sys.argv ) != 4:" - does? I know it's checking if 4 system arguments have been called but why? Is this in/out grid code that's automatically created wrong? the python of the two process mentioned below also has the same info repeated. Once it completed it should run Point Cloud Reclassifier / Subset Extractor where the input is the in_memory from the above process which is then used to run Point Cloud from Shape. === ORIGINAL POST === new process is only in python. The py, and bat files are in a zip in google docs. Thanks for your advise and help, Code snippets from bat file (I know there is the python saga_api but I can't seem to download it from http://www.saga-project.org/download/python as it gives a 404 error.               @ECHO OFF          SET SAGA_MLB = C:\SAGA\Modules     SET PATH = %PATH%;C:\SAGA          ECHO %1 %2 %3          PAUSE            C:\SAGA\saga_cmd io_shapes_las "Import LAS Files" -POINTS= -FILE=%3\%1_1K_class.las -T -i -a -r -c -u -n -R -G -B -e -d -p -C -VALID -RGB_RANGE=0          PAUSE     c:\SAGA\saga_cmd pointcloud_tools "Point Cloud to Shapes" -POINTS= -SHAPES=%2\%1_fc.shp           PAUSE      This is called from a python script that also uses arcpy               [CODE}     import os, datetime     from subprocess import Popen          LASList = ['329000_8119000','330000_8118000','330000_8119000','331000_8118000','331000_8119000']     #LASList = ['329000_8119000']     workingDIR = r'P:\2011\Job_183_Contours_DaneCraze\Working'     outDIR = workingDIR+"\\test3"     AOI=workingDIR+'\\Data_Request_Extents_region.shp'     SourceDIR=r'L:\Data_Admin\temp_dataset_storage\LiDAR\FinalSupply\Atherton_Tolga_Walkamin_Mareeba_Bilwon\LAS'     #SourceDIR=outDIR     spatial_reference = r'C:\Program Files\ArcGIS\Desktop10.0\Coordinate Systems\Projected Coordinate Systems\National Grids\Australia\GDA 1994 MGA Zone 55.prj'                    print "Starting process on "+ str(len(LASList)) +" tiles @ "+str(datetime.datetime.now())+ "\n"          for LASTile in LASList:          # SAGA Processes          # Import LAS              p = Popen("importLAS_exportSHP.bat "+LASTile+" "+outDIR+" "+SourceDIR, cwd=r"P:\2011\Job_183_Contours_DaneCraze\Working")         stdout, stderr = p.communicate()           print "Imported LAS @ "+str(datetime.datetime.now())+ "\n"               ### Convert LAS to shp     ##     ##    p = Popen("PointCloud_shape.bat "+LASTile+" "+outDIR, cwd=r"P:\2011\Job_183_Contours_DaneCraze\Working")     ##    stdout, stderr = p.communicate()     ##    print "Converted LAS to SHP @ "+str(datetime.datetime.now())+ "\n"     ##              # ARCGIS Processes          ### Set local variables     ##    in_features = outDIR+"\\"+LASTile+"_fc.shp"     ##    clip_features = AOI     ##    out_feature_class = outDIR+"\\"+LASTile+"_clip.shp"     ##    xy_tolerance = ""     ##     ### Execute Clip     ##    arcpy.Clip_analysis(in_features, clip_features, out_feature_class, xy_tolerance)     ##    arcpy.DefineProjection_management(out_feature_class, spatial_reference)     ##    print "Executed Clip @ "+str(datetime.datetime.now())+ "\n"     ##     ### Execute Subset     ##         ##    arcpy.MakeFeatureLayer_management(out_feature_class, "lyr")     ##    arcpy.SelectLayerByAttribute_management("lyr", "NEW_SELECTION", ' "CLASSIFICA" = 2 ')     ##    print "Selected Subset @ "+str(datetime.datetime.now())+ "\n"     ##         ##    arcpy.CopyFeatures_management("lyr", LASTile+"_class2")     ##    arcpy.DefineProjection_management("lyr", spatial_reference)     ##    print "Executed Subset @ "+str(datetime.datetime.now())+ "\n"     ##     ### Create XYZ - in SAGA     ##     ##    p = Popen("export_xyz.bat "+LASTile+" "+outDIR, cwd=r"P:\2011\Job_183_Contours_\Working")     ##    stdout, stderr = p.communicate()     ##    print "Created XYZ txt and exiting @ "+str(datetime.datetime.now())+ "\n"     ##          [/CODE}      One reason for me to post this is that I can't convert spc to shp on large clouds as I get the attached error. It runs fine on smaller tiles. ![enter image description here](http://i.stack.imgur.com/TEkWH.gif)


title: Where can I get (USA) county neighbour relationship data?

text: I need a graph where each node is a US county and each edge represents a border shared between counties. I do not particularly care about the absolute position or shape of each county (though that would be a plus). Where can I find that information (for free)? Army Corp of Engineers? US Geo Survey? Ideally, I could just get a csv list something like   * FL-Polk, FL-Lake   * FL-Polk, FL-Orange   * FL-Polk, FL-Osceola   * etc for every county, including "kitty-corner" borders, borders across water, or across state borders I'm a programmer so expect I can handle any standard exchange format. (That could just be hubris though). (I don't have any cool GIS apps. Just humble scripting in perl and python.)


title: What does ArcGIS Server error (The required port '4000 ' is already in use) indicate?

text: I tried to publish my map in ArcGIS for Server 10.1, but when doing so the following error appeared: The required port '4000 'is already in use.' For your information I use ArcGIS for Desktop 10.1.


title: QGIS - Unable to get OSM data on re-opening project

text: I have now successfully got OSM maps to import as I would like, but if I save the project then later re-open, I always get a bad layer message and I have to reload the OSM once the project is open. How can I save the datasource reference with the project so I don't get this every time?


title: Polygon to point conversion in ArcMap

text: I tried to convert polygons to point using `Feature to Point` under data management toolbox in arcGIS 10.2 (got it from uni. licensed). I got the point layer but it is not showing anything. I added X and Y fields and calculated geometry but still it is the same. I imported the coordinate system of the other shapefiles but no luck. I am just wondering where is the problem? Why is it not showing the points?


title: Importing GDAL from Python window of ArcGIS for Desktop?

text: I am trying to run some Python code using the Python window in ArcGIS 10.1 that uses both the `arcpy` and `gdal` modules. However, when I try and import the `gdal` module I get an error:               ImportError: No module named osgeo      Obviously it can't find the module, so I have added my main python site- packages directory to the `sys.path` list:               sys.path.append(r"C:\Python27\lib\site-packages")      However, when I try and run `import osgeo` or `from osgeo import gdal` I get the following error:               Runtime error      Traceback (most recent call last):       File "<string>", line 1, in <module>       File "C:\Python27\lib\site-packages\osgeo\__init__.py", line 21, in <module>         _gdal = swig_import_helper()       File "C:\Python27\lib\site-packages\osgeo\__init__.py", line 17, in swig_import_helper         _mod = imp.load_module('_gdal', fp, pathname, description)     ImportError: DLL load failed: %1 is not a valid Win32 application.      I've found various resources on the internet that seem to talk about similar problems (for example this question and this forum post, but they seem to be out-of-date (ie. not using ArcGIS 10.1) or not asking quite the same question. Of course, importing `osgeo` from a normal non-Arc Python window works fine, and I have compared `sys.path` and `os.environ['PATH']` between the scripts and updated them so that they are the same, and it doesn't seem to fix the problem. Does anyone have any idea how I can get this to work?


title: Is there a way to convert 3D polylines to 3D points in ArcGIS?

text: I need to convert a huge number of very dense (densified to 1m) 3D polylines to 3D points. I need each vertex along the line and no duplicates (two of the same XYZ coordinates) if possible. Is there a tool for that? Each feature class to be converted yields about 1 to 2 million points. I did write my own add-on for this and although I am using the insert cursor (as oppose to iFeature:Store) it's quite slow so I am hoping there is a better way.


title: How to combine polygons in a vector with GRASS?

text: I need to combine polygons of similar geological classes which are stored in a vector file called geology. I was looking for somethiing like "union selected features" or similar. In the GRASS Book 3ed I found that v.overlay could be useful, however I did not manage to merge/union my classes. To illustrate it, here an example: I have a polygon with class (column MAIN_LITHO) 'biotite granite' but also 'biotite granite (foliated) and others. These should be all combined into one single multipolygon 'biotite granite'. Later, I want to turn that vector into a raster map but that should be easy. here is also an image of that example ![enter image description here](http://i.stack.imgur.com/XF3v6.jpg) thank you in advance!


title: How to create reports in ArcGIS using Python?

text: Does anyone have experience in creating reports from your ArcGIS data (geodatabase, map exports) using a python reporting engine? In my current project, we are using Jasper Reports in the following way: A custom geoprocessing tool creates a intermediate XML file containing the data from the Oracle SDE and an exported map image as encoded byte array. This file and a template from iReport is used to generate a PDF file by Jasper Reports. I am looking for a different solution using a report engine based on python (Jasper is based on Java), I can imagine that a python approach can be better integrated in ArcGIS. Did anyone tried this before?


title: Is it possible to create an updateable spatial view?

text: 1- I have encountered a problem relating with sde commands for creating sde_ view and also sde_table . I used appropriate codes but an error massage occured during execution. the related codes:               sdetable -o create -t table22 -d"name string(20), value integer(5)" -i service -u username -p password -s server_name -q      the error is: > opration failed(-1) > > cannot get configuration information 2- I would like to know is it possible to create a spatial sde view ,register it by `-o register`,from an existing gdb, and replace it with the original attribute table of a feature class? I am looking for a possible solution to retain relationships belong to the related feature class.I need an updatable view being produced for esri products. I have done this with arcinfo toolbox extentions(with out relationships) .If yes how it is possible? The current enviroment is sde9.2, sql server 2008, arcgis9.3. thanks for anyone who reply


title: OpenLayers: Error message 404 (Not Found)

text: I've created a map using openlayers contains osm, google layers and some vector features, On loading it shows two errors like: **GEThttp://tile.openstreetmap.org/0/0/1.png 404 (Not Found)**. How to prevent this , I've googled it, and i can't understand what is the exact problem and how to solve this please help me.


title: Open source GIS tool for mapping health data

text: Im looking for an open source desktop gis tool that can be used for visualization and exploratory analysis. I want it to be interactive and perhaps display two or more maps at the same time to allow for comparative analysis. I know CommonGIS allows for this, but are there any others? thanks!


title: Creating color table file in PCI Geomatica

text: What is the simplest way to create a stand-alone color-table file in PCI Geomatica (PCT or STR)? I have some classified imagery with values ranging between 0-255 and I would like a way to automatically have them consistently colored the same way...


title: Extent in ArcGlobe

text: Does anyone know if it's possible to get an Envelope/Extent in ArcGlobe based upon what's visible on the globe. For instance, I'm looking to be able to zoom in over the state of Pennsylvania and get an Envelope back with that part of the globe? That way I could have a boundary-box/envelope of what's currently visible.... I've looked at GlobeDisplay.CurrentViewer.Camera I've looked at Camera.Target.QueryCoords(out x, out y) I've looked at GlobeDisplay.ActiveViewera.Camera.QueryEnvelope None of these seem to work. If anyone could point me in the right direction that would be awesome.


title: How can I have ArcGIS not automatically add raster clip output from Python window to the current map

text: I have a collection of rasters that I'm writing code to automate clipping based on a user defined shapefile. I am currently doing this in the Python window using the following as the main processing function:               arcpy.Clip_management( inputRasterFilename, "#", outputRasterFilename, "shapefile" )      Basically I want to read the raster, clip it, and have it saved to outputRasterFilename for later use. This is exactly what happens except that ArcGIS also adds the result of clip operation to the current map that I have open. Is there a way to disable this?


title: Select points within 2 buffers

text: I have 2 huge files of points (hundreds of thousands of points). I need to determine a way to loop through all the points in one file and find all the points in the second file that are only a certain distance away. So like within a distance outside the buffer. Or on the outer edge of the buffer. I am thinking of 2 buffers, one larger buffer and then one smaller buffer and then somehow subtracting the values. No idea how this could be accomplished though. The files are huge so I am also concerned in doing this the quickest way possible. Any ideas? Thanks Just Me


title: Is Orbit EOS Lite Free GIS Map Server an OpenSource software?

text: I'm writing an essay about web mapping software and I can't find out is Orbit EOS Lite Free GIS Map Server an OpenSource software or not? Also, I would be thankful if someone gives me some links with some list of web mapping programs. Thank you.


title: dwg from ArcGIS to Autocad 2012

text: I haven't used CAD in years and am quite rusty about its features. But a client who uses CAD has been experiencing difficulty lining up *dwg files exported from ArcGIS in their CAD workspace. How can I export tha data in GIS so that it is referenced to fit their workspace? Or how can they make adjustments to fit the data themselves? * * * Sorry I got caught up working on other projects. Thank you for the response guys. The client is not using Map 3D so projections are not important as I understand it from my CAD days. The problem is when the shapefile is exported in ArcGIS to *.dwg the client is asking for a reference point. According to Brad's instructionsI have another question. Does that mean we should agree on a point in the shapefile as a reference point that they have as well in their CAD data, and I then send them the coordinates for it to use to make adjustments in CAD? Does CAD even use geodetic coordinates?


title: Creating styles from scratch

text: I have been tasked to create a list of point symbols for the counties Watershed Improvement Program from scratch. They don't want modified icons from the programs database, but custom drawn. I already have the hand sketched symbols approved and I'm using Inkscape to create the vectors. Now to the actual question: Is there a step-by-step process from converting the .bmp or .png files to Esri symbols? Ie. what size and format should the vector be created in, how do you upload into Arcmap (10.1), and how do you access your folder in the Style Manager? I know this may seem simple to a lot of you, but the Googled help I've seen is about as simple as _Greek Calculus_ to me.


title: Adding field to table without path to table (python)

text: I want to add a field to my table in the FieldCalculator.               import arcpy     arcpy.AddField_management("veg8", "test", "TEXT")      An error message occured doing this and I figured out that I need to set a path to my table Like this for instance               C:\\user\\...      But, since I work in the Modelbuilder and the Model will be running on other computers, the path to the table will be a different one and I dont want to update all my paths everytime I run the model on a different computer. Is there a way to tell python to use the table to which the FieldCalculator is "connected"? Thanks Chris


title: XYZ Clustering in Python and numpy arrays

text: I am looking for a method of clustering XYZ data in Python that works on numpy arrays. I know of a few sources, such as clusterpy and Pysal but have had little success with them as they seem to still be in the development stage and do not support numpy arrays. EDIT: My numpy arrays are of XY points with a Z column of a value at the point such as soil moisture or plant height. I want to spatially cluster the Z values into similar regions.


title: Best practices to minimize error when buffering unprojected data

text: If I buffer hundreds of unprojected WGS84 points across the US by say, 3 miles, will I get accurate areas returned (3 x pi x r) for each of the circles created? I suspect not. If not, what would be some best practice to minimize error for each buffer. Basically, I'd like to minimize error across the US (perhaps including Hawaii) and processing time (in ArcGIS, for instance). For context, my goal is to extract demographics (also in unprojected geography) for each point lying within a constant distance of each point for the entire nation. My fear is that using unprojected data for this task means that if I'm comparing the demographic info captured by each buffer, I'm not comparing apples to apples. Any help is appreciated. Thank you.


title: Openlayers Map Marker / Popup Behaviour Mobile Devices

text:
Quick overview. i am creating a site with a map using Openlayers and Open Street Maps. User clicks on map and fills out a popup form and submits form. Markers appear at lonlat points where map was clicked and form filled out. Things were generally working as expected - when page loads it is populated with Markers. click on markers and view FramedCloud popup to see information entered by user. The only issue was after closing the popup i was not able to click on the same marker again (to create a popup) until i had clicked on either another popup or the map itself. i finally fixed this my adding the onPopupClose() method. After this everything was working as expected in the page on desktop browsers but now behaviour on mobile devices is not as exptected ( clicking the close popup cross closes the popup but makes the user form appear / multiple clicks on markers yeild no results until i click or drag on map and more odd behaviour). Does anyone have any ideas what changes i might need to make to get the behaviour as expected on Mobile devices. As im not certain what parts of the code are the most relevent to fix this issue ill post the whole thing. thank you i hope you can help or give me some guidance as to where to look. i also tried including openlayersmobile.js in the head of the index page but this didnt seem to make any difference               $(document).ready(function() {            map = new OpenLayers.Map("mapdiv");       map.addLayer(new OpenLayers.Layer.OSM());       navigator.geolocation.getCurrentPosition(function(position) {                var lonLat = new OpenLayers.LonLat(position.coords.longitude,position.coords.latitude)             .transform(                 new OpenLayers.Projection("EPSG:4326"), //transform from WGS 1984                 map.getProjectionObject() //to Spherical Mercator Projection                 );                        map.setCenter(lonLat, 17 // Zoom level         );            });       OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {     defaultHandlerOptions: {         'single': true,         'double': false,         'pixelTolerance': 10,         'stopSingle': false,         'stopDouble': false     },     initialize: function(options) {         this.handlerOptions = OpenLayers.Util.extend({}, this.defaultHandlerOptions);         OpenLayers.Control.prototype.initialize.apply(this, arguments);         this.handler = new OpenLayers.Handler.Click(this, {         'click': this.trigger         }, this.handlerOptions);     }            });          // Create click control       var clickControl = new OpenLayers.Control.Click( {      trigger: function(e) {       var x = e.pageX;       var y = e.pageY;         $('#popupform').css({         "position": "absolute",         "left": x,         "top": y,         "visibility": "visible"         });        var toProjection = new OpenLayers.Projection("EPSG:4326");       lonLat = map.getLonLatFromPixel(e.xy).transform(map.getProjectionObject(), toProjection);       longi = lonLat.lon;       lati = lonLat.lat;     }       });       map.addControl(clickControl);       clickControl.activate();                      $('#when_spotted').datepicker({ dateFormat: 'yy-mm-dd' }).val();       $('#time_spotted').timepicker();       var frm = $('#sub');                      frm.click(function() {            });             $("#fileuploader").uploadFile({         url:"/Streetfind/streetfind.php",         fileName:"file",         dynamicFormData: function () {           var data = $('#popform').serializeArray();           data.push({name: 'longitude', value: longi},{ name: 'latitude', value : lati});             return data;         }        });                 $('#popupform').submit(function(){     $('#popupform').css('visibility','hidden');     return false;       });            $( "#close_box" ).click(function() {     $( "#popupform" ).css('visibility','hidden');       });                      var customIcons = 'http://labs.google.com/ridefinder/images/mm_20_blue.png';            var vector = new OpenLayers.Layer.Vector('vector layer', {     styleMap: new OpenLayers.StyleMap({       externalGraphic: 'http://openlayers.org/dev/img/marker.png',       graphicWidth: 30, graphicHeight: 34, graphicYOffset: -10,       title: '${tooltip}'     })       });       map.addLayer(vector);                 var kontrol = new OpenLayers.Control.SelectFeature(vector,     {onSelect: onPopupFeatureSelect, onUnselect: onFeatureUnselect});       map.addControl(kontrol);       kontrol.activate();                      downloadUrl("phpsqlajax_genxml3.php",  function(data) {     var xml = data.responseXML;     var markers = xml.documentElement.getElementsByTagName("marker");     for (var i = 0; i < markers.length; i++) {           var item_name = markers[i].getAttribute("item_name");       var item_description = markers[i].getAttribute("item_description");       var item_condition = markers[i].getAttribute("item_condition");       var when_spotted = markers[i].getAttribute("when_spotted");       var time_spotted = markers[i].getAttribute("time_spotted");       var point = new OpenLayers.Geometry.Point(         parseFloat(markers[i].getAttribute("longitude")),         parseFloat(markers[i].getAttribute("latitude"))).transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject());       var loon = markers[i].getAttribute("longitude");       var laat = markers[i].getAttribute("latitude");       var file = markers[i].getAttribute("file");       var road = markers[i].getAttribute("road");       var html = "" + road + "
" + item_name + "
" + item_description + "
" + item_condition + "
" + when_spotted + "
" + time_spotted + "
" + '';       var size = new OpenLayers.Size(15,22);       var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);       var icon = new OpenLayers.Icon('http://labs.google.com/ridefinder/images/mm_20_blue.png',size,offset);                 var features = new OpenLayers.Feature.Vector(point,{         description:   html        });            vector.addFeatures(features);     }       });            function onPopupClose(features) {     var featured = features;      if (featured.layer) { // The feature is not destroyed        kontrol.unselect(feature);          } else { // After "moveend" or "refresh" events on POIs layer all               //     features have been destroyed by the Strategy.BBOX               this.destroy();             }          }                      function onFeatureUnselect(features) {     feature = features.feature;     map.removePopup(features.popup);     features.popup.destroy();     features.popup = null;       }                        function onPopupFeatureSelect(features) {     feature = features.feature;     selectedFeature = feature;               popup = new OpenLayers.Popup.FramedCloud("popup",     features.geometry.getBounds().getCenterLonLat(),     null,features.attributes.description, null, true, onPopupClose);     popup.panMapIfOutOfView = true;     popup.autoSize = true;     features.popup = popup;     map.addPopup(popup);            }                           function downloadUrl(url,callback) {     var request = window.ActiveXObject ?     new ActiveXObject('Microsoft.XMLHTTP') :     new XMLHttpRequest;          request.onreadystatechange = function() {       if (request.readyState == 4) {         request.onreadystatechange = doNothing;         callback(request, request.status);       }       else {         console.log('this was the status ' + this.status);       }     };          request.open('GET', url, true);     request.send(null);       }                      function doNothing() {}           });


title: QGIS - select/aggregate by attributes and calculate

text: I need to know how I can select/aggregate equal attributes in a table. The column "Number" has attributes like 1.1 or 1.3 or 4.3. There are a lot of 1.1 or 1.3....in the column. In step1 I want to find all the 1.1 or 1.3 ...for calculating later the weighted average of a value by using the area. The problem is that the "Numbers" are different in all my tables. In EXCEL I could use a data filter to select all the equal attributes. Its not easy to describe...sorry. To say it in one sentence: If „Numbers“ are equal than calculate the weighted average of the Values and create a new table ![screenshot_explanation](http://i.stack.imgur.com/6tt13.jpg) But hopefully the screenshot describes the problem better than my words.


title: how to Find n nearest lines for the given points from a collection of lines?

text: I have two tables, points and lines. For each point in point table, I want to find out nearest 3 points in the nearest 3 cables. I almost developed a query, but it is best suitable for single point, but i need a query for all the points in the point table. My query is               select * from (     with index_query as (     select l.data,p.data, ST_Distance(c.geom,p.geom) as distance,     ST_AsText(st_line_interpolate_point(l.geom,      st_line_locate_point(l.geom, p.geom))) as nearest_point_in_line,      from lines l,points p order by l.geom <#> p.geom limit 100)     select *,RANK() OVER (PARTITION BY pid order by distance ) AS r     from index_query ) as s where r < 4      LIMIT 100 is the limiting factor, which is limiting total output rows to 100 irrespective of number of points in Point table. I need nearest 3 point_in_line for all points in my table. Total 3*no of rows in point table. How to achieve it?


title: What is exactly Geotrigger service by ESRI ?

text: I am going through the FAQ's of ESRI latest Geo-trigger service & according to FAQ it provides location-based alerts for mobile platform. Also I am looking for an example regarding this on Google and ESRI site but not able to find so any help regarding explanation with example will be great !!! Thanks in advance !!


title: TileMill: How to access geojson properties with CartoCSS

text: When trying to display a geojson layer in TileMill, is it possible to access the 'properties' of each feature? My geojson looks like this:                       "type": "Feature",              "geometry": {                 "type": "Point",                  "coordinates": [2.22614, 48.8765]                 },              "properties": {                 "type": "stop",                 "name": "SOME_NAME",                 "color": "#ff0000"                 }      and I would like to display the name in properties/name, using the color in properties/color. Is it possible? Properties are clearly accessible since one can write this and it works:               #id[type='stop'] {       /* various css rules */       }      But is it possible to use properties in the css rules themselves? This for instance does not work:               #id[type='stop'] {       marker-line-color: [color];       }      (or any variation with or without quotes inside or outside brackets for 'color'). If it's not possible, in what other format should I transform the geojson in order to access properties of each element from the css?


title: How to subset X number of features from Shape file in QGIS

text: I'm looking for a means to take a rather large Shape file and split it out into separate, smaller shapefiles based on number of features. So for example, if I have a file with 1 million features, subset it out into smaller files of 100k features each. I know how to grid a shapefile and subset into smaller files, but I don't want to grid it as I want to preserve the original feature geometry (e.g. gridding polygons will ruin the geometry). Does anyone know how this can be done with QGIS?


title: How to couple GeoTools GML bindings for encoding to GML 3.1.1

text: I've successfully created GeoTools (8.3) bindings for parsing and encoding XML bindings as described in their tutorial: http://docs.geotools.org/latest/userguide/library/xml/internal/bindings.html I can parse the XML data + GML 3.1.1 geometries successfully, including combining multiple XSD's and using the binding context for factories, etc.. but the encoding seems to result in XML data with geometries that are not linked to a particular namespace ('null:'). Since the above tutorial (or anything else I could find) doesn't mention specifics about encoding GML- geometies I have the following question(s):   * _Which bindings do I have to make myself and which can be handled by GeoTools?_   * _In the repository I see a number of util classes (e.g. GML3EncodingUtils), Geometry type bindings (implemented) and normal geometry bindings (unimplemented). Does anyone have experience using these? Again, which to implement, and which not?_ I'm only interested in a limited part specifying the geometries. Here's the fragment of the XSD that connects to GML geometries (Dutch Cadastral data):               <?xml version="1.0" encoding="UTF-8"?>     ...     <xs:import namespace="http://www.opengis.net/gml" schemaLocation="../../gml/bag-gml.xsd"/>     ...     <xs:complexType name="PuntOfVlak">         <xs:annotation>             <xs:documentation>Keuze tussen een punt (gml:Point) of een vlak (gml:Surface).</xs:documentation>         </xs:annotation>         <xs:choice>             <xs:element ref="gml:Point"/>             <xs:element ref="gml:_Surface"/>         </xs:choice>     </xs:complexType>     ...     <xs:complexType name="VlakOfMultiVlak">         <xs:choice>         <xs:element ref="gml:_Surface"/>         <xs:element ref="gml:MultiSurface"/>         </xs:choice>     </xs:complexType>     ...      So basically, its about how to implement the encode(..) or getProperty(..) method of my PuntOfVlakBinding and VlakOfMultiVlakBinding. The original XSD's can be found here:   * The GML XSD: http://pastebin.com/nMisVg2Q The VlakOfMultiVlakBinding:               ...          public class VlakOfMultiVlakBinding extends AbstractComplexBinding {              private static final String SURFACE = "_Surface";         private static final String POLYGON = "Polygon";         private static final String MULTI_POLYGON = "MultiSurface";              @Override         public Object getProperty(Object object, QName name) throws Exception {                  if (SURFACE.equals(name.getLocalPart())) {                 return object;             }                  if (POLYGON.equals(name.getLocalPart())) {                 return object;             }                  if (MULTI_POLYGON.equals(name.getLocalPart())) {                 return object;             }                  return null;         }              @Override         public QName getTarget() {             return BAGType.VlakOfMultiVlak;         }              @Override         public Class<?> getType() {             return Geometry.class;         }              @Override         public Object parse(ElementInstance instance, Node node, Object value) throws Exception {                  if (node.hasChild(SURFACE)) {                 return node.getChildValue(SURFACE);             }                  if (node.hasChild(POLYGON)) {                 return node.getChildValue(POLYGON);             }                  if (node.hasChild(MULTI_POLYGON)) {                 return node.getChildValue(MULTI_POLYGON);             }                  throw new IllegalArgumentException("Unsupported GML type detected for binding " + VlakOfMultiVlakBinding.class.toString());         }     }      I can also post attempts at the _SurfaceBinding.java if that helps **Edit** : This post already scratches the surface of the use of bindings, but doesnt elaborate yet: http://gis.stackexchange.com/a/6430/14253 **Edit 2** : I've implemented my own bindings where each element with subelements now has it's getProperty() method implemented and each 'end'-element it's encode() method implemented: What this looks like for the polygon binding tree:               <gml:Polygon srsName="...">                            --> PolygonBinding + getProperty()         <gml:exterior>                                     --> ExteriorBinding + getProperty()             <gml:LinearRing>                               --> LinearRingBinding + getProperty()                 <gml:posList srsDimension="2" count="..."> --> PosListBinding + encode()      This seems to work, encoding wise, but I'm still stuck with the 'null:' namespace. Rather than the expected 'gml:' namespace. This is strange, because debugging through the bindings shows the namespaceURI is at all times known as "http://www.opengis.net/gml". The current result for the VlakOfMultiVlak binding for a test case:               <bag_LVC:woonplaatsGeometrie>         <null:Polygon srsName="urn:ogc:def:crs:EPSG::28992">             <null:exterior>                 <null:LinearRing>                     <null:posList count="6" srsDimension="2">193273.335 389233.922 194582.516 389880.369 194456.597 389867.675 194207.272 389712.575 194201.368 389710.434 193273.335 389233.922 </null:posList>                 </null:LinearRing>             </null:exterior>         </null:Polygon>     </bag_LVC:woonplaatsGeometrie>


title: How to set ODBC connection in MapBasic

text: i am trying to set code for ODBC connection in Mapbasic. Instead of using Username and password i wanted to set code for Window authentication.my sample code is               Register table "tablename"     Type ODBC     Table "Select * From Database.tablename"     CONNECTION     "DSN=Database; UID=Username; PWD=Password"     into "C:\temp.table.tab"      so in this case instead of UID and PWD i wanted to set as window authentication to connect to database...anybody knows this.


title: Versioning with ArcSDE - Is it possible to cancel or reject edits once they've been posted?

text: I am using ArcGIS 9.3.1 and attempting to work with an SDE geodatabase (with one polygon feature class) that has already been registered as versioned. I'm new to versioning and am still trying to figure out some of it's basic functions. So far, I haven't been able to discover if it is possible to "cancel" or "reject" certain edits once they've been posted to a parent version. For example, let's say we have three versions: the original SDE.DEFAULT that was created when it was registered as versioned, a child version of the default called SDE.QA (for Quality Assurance), and a child version of the QA called SDE.Edit1 (where the edits first take place). If certain features of SDE.Edit1 were edited (e.g., to keep it simple, let's say one polygon was added and one was removed) and then SDE.Edit1 was reconciled with the SDE.QA and subsequently posted to SDE.QA, would there be any way to later undo this change? Following up on this question, would it possible to reject only _some_ changes? For example, accepting adding the first poly, but rejecting removing the second poly? As far as I can tell, once edits have been posted to the parent version, all these changes are now a "permanent" (for lack of a better word) part of the parent version. I am aware of the fact that these changes are all recorded within two tables, the "ADD" and "DELETE" tables (often referred to as the "delta" tables), and don't actually change the original FC itself. I considered looking into manually altering these delta tables, but I found enough people warning against that to know that it's likely not the right solution. Perhaps it is my understanding of versioning that needs some work, but I couldn't seem to find a way to reject a change or undo a change once it's been posted. This seems strange to me, as this would mean that there is no way to undo a post which contained an error. I also can't seem to find a way to trace the lineage of these versions (i.e., which version is the child of which parent). While I'm on the topic, if anyone might know of any particularly useful ArcSDE references (links, articles, books, etc) that might help with my understanding of ArcSDE (and perhaps answer some of these questions), it would be much appreciated! ** **UPDATE** ** Although answers so far have been helpful (thank you for the links), I still cannot find an answer to the core of my question. Again, perhaps it is my own misunderstanding of the situation. Here is what I want to know: **Can you reverse (by reverse, I mean _undo_ ) a post once it has been made from a child version to a parent version? In this scenario, the parent can be, but doesn't have to be, the SDE.DEFAULT version. Even better, I would like to know if you can reverse a _part_ of a post (say, a single edit to a polygon), after it has been posted? I would also like to know if this can be done without the need to have had any conflicts detected.** The fact that I can't find a clear answer to this question (i.e., "yes" or "no") documented anywhere makes me think I am missing something important about versioning in the ArcSDE. I would also prefer to avoid manipulating the A and D tables manually.


title: How to use Google Default Pan/Zoom controls in openlayers

text: I am using Google maps version 3 with Openlayers,My question is how do you use the default google control for Panning and Zooming and not the default controls Openlayers panning and zooming(The dark blue controls) I know that it is controlled by the CSS for Openlayers, But what is the default CSS for the Google Pan/Zoom Controls / Where do I find this CSS for the Google controls in the image below. ![Google Default Controls](http://i.stack.imgur.com/Sv4dS.png) ![Openlayers Pan/Zoom Control](http://i.stack.imgur.com/mKswg.png)                var map, layer;                  function init(){                  layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",                         "http://vmap0.tiles.osgeo.org/wms/vmap0", {layers: 'basic'} );            // Add Styles       var myStyles = new OpenLayers.StyleMap({                 "default": new OpenLayers.Style({                     fillColor: "#000",                     fillOpacity: 0.3,                     strokeColor: "#000",                     strokeOpacity: 0.3,                     strokeWidth: 1,                     graphicZIndex: 1                 }),                 "select": new OpenLayers.Style({                     fillColor: "#ccc",                     fillOpacity: 0.6,                     strokeColor: "#ccc",                     strokeWidth: 1,                     label: "Label for Mouseover",                     labelAlign: "cc",                     fontColor: "#333333",                     fontOpacity: 0.9,                     fontFamily: "Arial",                     fontSize: 14,                     graphicZIndex: 2                 })             });                  vlayer = new OpenLayers.Layer.Vector( "Editable",{                 styleMap: myStyles,                 rendererOptions: {zIndexing: true}             });                  map = new OpenLayers.Map( 'map', {                  // Add Openlayers Map Controls              controls:                            [                         new OpenLayers.Control.LayerSwitcher(),                         new OpenLayers.Control.PanZoom(),                         new OpenLayers.Control.MousePosition(),                         new OpenLayers.Control.ScaleLine(),                         new OpenLayers.Control.OverviewMap(),                         new OpenLayers.Control.EditingToolbar(vlayer)                            ]                 });          // Make Variables for Google Version 3      var gphy = new OpenLayers.Layer.Google(         "Google Physical",         {type: google.maps.MapTypeId.TERRAIN}     );     var gmap = new OpenLayers.Layer.Google(         "Google Streets", // the default         {numZoomLevels: 20}     );     var ghyb = new OpenLayers.Layer.Google(         "Google Hybrid",         {type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20}     );     var gsat = new OpenLayers.Layer.Google(         "Google Satellite",         {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22}     );              // Add layers to map         map.addLayers([gphy, gmap, ghyb, gsat]);              map.addLayers([layer, vlayer]);              // Google.v3 uses EPSG:900913 as projection, so we have to transform our coordinates         map.setCenter(new OpenLayers.LonLat(10.2, 48.9).transform(         new OpenLayers.Projection("EPSG:4326"),         map.getProjectionObject()     ), 5);             var panel = new OpenLayers.Control.Panel();         panel.addControls([             new OpenLayers.Control.Button({                 displayClass: "helpButton", trigger: function() {alert('help')}, title: 'Help'             })         ]);          map.addControl(panel);


title: Why is transparancy of GRASS rasters missing when imported into QGIS?

text: I have a raster map in QGIS imported from a GRASS mapset with the QGIS GRASS plugin. The problem is that I have null() cells in the GRASS map that appear "colored" in QGIS. How can I set the transparency? I have tried in 'Proprieties' but it seems that it doesn't work. I would like to know a simple way to solve my problem.


title: 3D modelling a snow patch

text: In my studies, I work with the internal dynamics of different snowpatches. In this study I wish to generate 3D models of snow patches. I have some problems when it comes to the actual work flow of creating this 3D model: Data I have available are GPS points (x, y, and z values) over large parts of the research area (these points are collected by using a Real Time Kinematic GPS, on the surface of the snow patches) and GPR (ground penetrating radar) data, which tells me what's under the ground. For example, how deep is the ice, how thick are the snow layers, etc. My thought is that I want to create a TIN from gps points using height values. The problem here is that I do not have GPS points for the whole research area, about 80% is covered. Question 1: Are there any interpolation tools that can estimate the points where I do not have the point data, based on the points I already have? When it comes to the actual 3D modeling, that I have a TIN of the glacier, but will now set the GPR data. Question 2: Would it be best to create eg. three shapefiles one for ice, one for snow and one for firn? Then put different base heights on these files according to the Z value in the GPR data? Or can I do this in a simpler way? Highly appreciate all the help I can get. Regards, Jorn


title: projection discrepancy with gdal (1.9) and proj4

text: I boiled the problem down to that python example:               import sys     try:     from osgeo import gdal     from osgeo import osr     except:     import gdal     import osr          if __name__ == '__main__':     version_num = int(gdal.VersionInfo('VERSION_NUM'))     print( "GDAL Version: " + str(version_num))          wkt = """     PROJCS["Popular Visualisation CRS / Mercator",         GEOGCS["Popular Visualisation CRS",             DATUM["Popular_Visualisation_Datum",                 SPHEROID["Popular Visualisation Sphere",6378137,0,                     AUTHORITY["EPSG","7059"]],                     TOWGS84[0,0,0,0,0,0,0],                 AUTHORITY["EPSG","6055"]],             PRIMEM["Greenwich",0],             UNIT["degree",0.0174532925199433],             AUTHORITY["EPSG","4055"]],         PROJECTION["Mercator_1SP"],         PARAMETER["central_meridian",0],         PARAMETER["scale_factor",1],         PARAMETER["false_easting",0],         PARAMETER["false_northing",0],         UNIT["metre",1,         AUTHORITY["EPSG","9001"]],         EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],         AUTHORITY["EPSG","3785"]]"""          X = 2662981.191207001451403     Y = -1935374.062384711578488          srs = osr.SpatialReference()     srs.ImportFromWkt(wkt)          srsLatLong = srs.CloneGeogCS()     ct = osr.CoordinateTransformation(srs, srsLatLong)     (int, lat, height) = ct.TransformPoint(X, Y)          # Report results     print('X: %s\t\tY: %s' % (X, Y))     print('Lat: %s\t\tLon: %s' % (gdal.DecToDMS(lat, 'Lat', 2), gdal.DecToDMS(int, 'Long', 2)))      my output is:               GDAL Version: 1900     X: 2662981.19121                Y: -1935374.06238     Lat:  17d 1' 1.01"S             Lon:  23d55'19.08"E      However, when I run this from proj4, I get a different answer:               invproj +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs     2662981.191207001451403 -1935374.062384711578488     23d55'19.081"E  17d7'29.797"S      The latitude is off by 6'28" or so... ?????? What could be going wrong? Many thanks, Pat.


title: How to reduce the file size of GeoJSON data?

text: For my project I am trying to convert the tiger shapefile information for us states into GeoJSON format However when I use the ogr2ogr tool the data file is about 20 MB. I have seen the geoJSON file of US States elsewhere of around 88KB. How do I reduce the output to geoJSON file?


title: How to find names of Unincorporated parts of a county in Navteq data?

text: I am currently using Navteq Data and am unable to find how to label unincorporated places within the LA County. Most of these spots are small and spread throughout the entire county however they do have well known names, most times these are same as neighboring cities. Has any one had similar problem? i.e. how to label the unincorporated parts of any county? Thanks, Shakira


title: How to generalize polygon attributes in a regular polygon grid?

text: I'm trying to generalize the attributes of a set of polygons(let's say they're census blocks with a population property) into another, overlapping set of polygons(let's say it's just a grid with no properties) _I'm certain there's a word for this problem that I've never heard of and never thought to search... If that's so, please just point me in the right direction!_ I want to make a choropleth map that shows the population density of the area without showing the shape of the census blocks themselves. I want to put their attributes into the regular grid. BUT, there's of course a lot of overlap, and irregular shapes don't fit neatly within a square. I need something that can:   1. Assume that population is evenly distributed across the census block polygons.   2. for each square in the grid, calculate..      1. what polygons are even a little bit under it, and for each one, what portion of that irregular polygon is indeed under it.     2. Then, lets say we have an irregular census block polygon that has a population of 100, and is 37% covered by the grid-block of interest... It would need to calculate the amount of the population attribute to assign to the block, in this case, exactly 37, assign it to that block, and then do the same, adding for every other shape under that block before moving on.    3. Then do that for all the other polygons in the grid layer.  I should end up with the exact same total population, just abstracted and generalized into regular shapes. How can I do this, particularly in QGIS? I really just want to count up some point data (points in polygon) and create a choropleth corrected for population density.


title: How to reproject from Michigan Georef to WGS84 in QGIS?

text: I am trying to reproject a michigan state land shapefile which is the NAD83/ Michigan Oblique Mercator EPSG:3078 to the WPG84 EPSG:4326 projection. However when i try to do this my map gets squished down. I have tried the commands in "Reprojecting from MGRm to WGS84 using open-source tools" with no luck. I need to get these coordinates into a decimal degrees format. Thanks


title: How to set up a plugin repository for QGIS?

text: I am working with python plugins for QGIS.I developed my plugin and created .zip file including my .ui,.py and resources files.In order to incooperate this plugin into QGIS i created .xml file plugins.xml as shown below:               <?xml version = '1.0' encoding = 'UTF-8'?>      Example #1 testing 0.1 1.0 ftp://softwarindia.com               <file_name>about.zip</file_name>     <author_name>aaaa</author_name>          <download_url>ftp://softwarindia.com/abt/about.zip</download_url>     <uploaded_by>me</uploaded_by>     <create_date>2011-08-31</create_date>     <update_date>2011-09-11</update_date>      I uploaded my .zip file and .xml file on ftp://softwarindia.com. In qgis i went to 'fetch python plugins'-->repositories.Added repository as name:About and url:ftp://softwarindia.com/plugins/plugins.xml.but when i go to plugins nothing is displayed in it. what can be the problem??


title: Openlayers - Loop through control collection

text: Is there a way to loop through the collection of map controls that have been added to an OpenLayers map?


title: Is it possible to try the ESRI EDN without purchasing a license?

text: Albeit using QGIS for most of my GIS tasks, I still use ArcMAP for final map production. I am curious and would like to test out EDN. According to ESRI the annual subscription is $1,500. I'm not sure I'm ready to invest that much money into it as I'm not sure I'll even use it. Do I have any options to develop on the platform without subscribing to EDN?


title: How to get lat/long from Google Analytics when country but no city given?

text: About a year ago, I created a web-based map of the visitors to a web site. I've been asked to switch it from using Urchin and web server log data, to retrieving data from Google analytics. When I was using the older data sources, I had the IP address of each visit, and I would use the IP address to look up city name / country code / latitude / longitude. Sometimes the city was unknown in the database I was using, but I would almost always get the country code, with a centroid latitude/longitude. This is important data, because I would put markers on a map showing where visitors were from, and sometimes a large proportion of visits from a given country were from unknown cities. For example, I have results like               Country: Anguilla  Visits: 2  City: (not set) Latitude: 0.0000 Longitude: 0.0000      This is the **only** data point from Anguilla. So if I don't plot it, Anguilla is completely omitted. Now with the Google Analytics Data Export API, IP address is unavailable. Thankfully, city name / country name / latitude / longitude are available most of the time, so I don't have to look those up by IP. But in the cases when the city is unknown, but country name is given, lat/long is given as 0.0000,0.0000. How can I plot these country points on a map? It's difficult to look up country locations in a database because the country **name** is given, rather than a country **code**. Country names of course are somewhat variable. It really surprises me that Google would not make country code available at all, but I don't see it listed under the allowed "dimensions" for an API query. I guess I could build up my own lookup table, fine-tuned to match Google's country names, but it would be subject to breakage if Google ever changed spelling or phrasing (e.g. "Democratic Republic of the Congo" vs. "Congo- Kinshasa"). Thanks for any suggestions.


title: How to change the the style of a vector layer programmatically?

text: I want to change the style of the polygon contained in the vector layer. Is there any API available? (I can click the properties of the shapefile, then select the style and change the symbol properties, but I want to do this programmatically.) I am using Python to customise QGIS.


title: Are there any freely available WMS tile servers?

text: I am looking for an open public endpoint to retrieve WMS map tiles from, for use with leaflet maps For clarity, I am not looking to set up a tile server or for tile server software. It looks like the NASA service is now closed


title: Why doesn't the split tool return any records in the attribute table?

text: I've been trying to split a large land use shapefile using a set of study area polygons as the split feature, but when the function finally finishes running (the land use dataset is quite large), all of the shapefiles created by the split are empty. We also tried using another computer, but found the same result, both using ArcGIS 10 with SP2. Is there any reason why this is the case? Are there any alternate methodologies which would achieve a similar result to split? Edit: To be more specific: the split field is the FID in text form and we didn't run into any errors relating to the input fields. The study area polygons completely contain the land use.


title: Using the US census to analyze data nationally

text: I'm trying to make a table to relate the STFID and LOGRECNO values for the US 2000 census at the block group level, nationally. I realize that the LOGRECNO value is a state specific unique variable, so I am joining it to the state ID code. The steps that I have done so far are: Importing all the national `geo` files into a database (e.g. akgeo.txt,.... to wygeo.txt) Extracting all the blockgroup values using the following `SQL` command:               SELECT SF1GEO.LOGRECNO, SF1GEO.REGION, SF1GEO.STATE, SF1GEO.COUNTY, SF1GEO.TRACT, SF1GEO.BLKGRP     FROM SF1GEO     WHERE SF1GEO.BLOCK IS NULL      AND     SF1GEO.BLKGRP IS NOT NULL     AND      SF1GEO.TRACT IS NOT NULL     AND      SF1GEO.COUNTY IS NOT NULL;      I want to ignore all the block level information, and only take the relevant LOGRECNO for each blockgroup (BLKGRP). Then to produce the STFID I concatenate the variables:               STATE + COUNTY + TRACT + BLKGRP       which gives me a 12 character code (STATE is 2 characters, COUNTY is 3 characters, TRACT is 6 characters, BLKGRP is 1 character). I also combine the `STATE` and `LOGRECNO` values for future analysis, based on advice from census.gov as I want to perform national analysis.                 LOGRECNO REGION STATE COUNTY  TRACT BLKGRP        STFID LOG_STATE     1  0000038      3    01    001 021100      2 010010211002 010000038     2  0000077      3    01    001 021100      3 010010211003 010000077     3  0000107      3    01    001 021100      1 010010211001 010000107     4  0000169      3    01    001 021100      2 010010211002 010000169     5  0000231      3    01    001 021100      3 010010211003 010000231     6  0000308      3    01    001 021000      2 010010210002 010000308      At this point I have 546293 values, but only 208144 unique STFID values. One problem I have is that I can't figure out why all my STFID values are not unique. Then, when I try to merge this table with a shapefile of blockgroups it only matches on 72972 polygons. I also can't figure out why the STFID does not match all. This is what the map looks like: ![enter image description here](http://i.stack.imgur.com/HKLBz.png) I'm not sure if there is an obvious pattern to the blockgroups that are included/excluded. * * * EDIT: I used Sean's suggestion of selecting by `SUMLEV` asthat makes a lot more sense than what I was doing earlier. So my SQL query is now:               SELECT SF1GEO.LOGRECNO, SF1GEO.STATE, SF1GEO.COUNTY, SF1GEO.TRACT, SF1GEO.BLKGRP     FROM SF1GEO     WHERE (((SF1GEO.SUMLEV)="150"));      However, I'm still mystified when I join this to the shapefile as I end up with the same map as above, and the same number of joins. I am manually making the code to join within the shapefile using the information within by concatenating STATE + COUNTY + TRACT + BLKGRP. What I can't understand is why 72972 values successfully join, and the rest do not.


title: Problem with runing a loop in a query in PostGIS

text: I am new to PostGIS and I have a question which can't solve. I have a table of buildings each having a value for number of people living in it. I would like to calculate the total value (sum of number of people) in all buildings which fall into a defined circular buffer. I have written this code for it which uses a loop, however it gives an error which is followed. It would be highly appreciated if you can help me regarding this problem. I run the following query in PgAdminIII. query:               BEGIN     FOR building_pop IN SELECT total_pop FROM Buildings WHERE ST_DWithin(ST_Transform(geom, 26918), (SELECT ST_Transform(ST_GeomFromText('POINT(1116858.6062789 7086290.6680572)', 900913), 26918)), 1000)         LOOP         totalpop = totalpop+building_pop;         NEXT building_pop;         END LOOP;         UPDATE public.h_grid200m SET grid_pop = totalpop WHERE gid =1       END      I receive the following error: > ERROR: syntax error at or near "FOR" LINE 2: FOR building_pop IN SELECT > total_pop FROM h_buildings WHERE ...


title: Runtime Error 424 when embedding MapInfo map in Access form

text: I'm kind of new to VB and DAO so please bear with me on this! I'm not sure if this belongs in gis.se or stackoverflow. I am trying to embed a simple MapInfo map in an Access form based on the code provided in the MapInfo user guide. It seems I've fallen at the first hurdle. So far I have this:               Dim mi As Object      Private Sub Form_Open(Cancel As Integer)          Set mi = CreateObject("MapInfo.Application")          MsgBox (mi)          mi.do "Set Application Window " & Form1.Hwnd          mi.do "Set Next Document Parent " & Form1.Hwnd & " Style 1"          mi.do "Open Table ""World"" Interactive Map From World"          mi.RunMenuCommand 1702          mi.do "Create Menu ""MapperShortcut"" ID 17 As ""(-"" "      End Sub       I am getting runtime error 424 "object required" on this line: `mi.do "Set Application Window " & Form1.Hwnd` The MsgBox outputs "MapInfo Professional" which leads me to believe that the CreateObject command has worked successfully. However if I comment out the problem line I get the same error on the next line. Therefore I assume that the mi variable isn't correctly set up. Can anybody shed any light on this? Are there any good resources on line for MapInfo integrated mapping? Thanks in advance!


title: Edit SDE table with ArcObjects

text: The link below provides a code snippet for editing an SDE table with ArcObjects. Are there any other or newer code examples out there of how to edit a SDE table by using ArcObjects? Edit SDE geodatabase table w/ ArcObjects Thank you!