title: ST-Links convert Feature Class to DB Table error

text: I'm using ArcMap with the ST-Links SpatialKit Conversion tool to convert a Feature Class to a DB Table. But I'm getting an error: ERROR: 42601: type modifier is not allowed for type "geometry". I can't seem to find anything on that error. Feature Class (shapefile)               Geometry Type : Polygon     Projected Coordinate System : WGS_1984_UTM_Zone_21S     Projection : Traverse_Mercator     Geographic Coordinate System : GCS_WGS_1984      DB               PostgreSQL 8.3     Postgis 1.5      Any ideas? I have other tables in the DB that were generated with Quantum GIS. I also can't open the other spacial tables on ArcMap correctly. In Quantum GIS I have no problems. But I need to work with ArcGIS. Thanks,


title: How to start from scanned maps in qgis

text: What I have: A large TIFF covering parts of several contiguous sheets of the British OS 6 inch first edition. It has no georeferencing. It's just a series of scans which someone has stitched together (actually gdal_transform extracted it from a larger .ecw) A number of scans of original sheets. These are also plain scans with no georeferencing but they do have latitudes & longtitudes along the borders so I can read off the corner coordinates in DMS. I assume that the projection of all this data will be close to if not identical to the current OS projection. What I want to achieve: The large TIFF properly georeferenced so it can be used as a basis for georeferencing various tithe maps. How I expected to go about it: Set the existing CRS of one of the scanned original sheets to OSGB70 which, as far as I can see is the OS projection using lat/long, fire up the georeferencer, set GCPs at each corner, set the output to OSGB70 & run the transform. Repeat for other sheets. Using these georeferenced sheets pick off GCPs for features on the big map to use the georeferencer to georeference this. What happened: A grossly distorted map was produced. I think I can see what the problem is; the original is in the OS projection but its CRS is purely pixels & the georeferencer treated it as degrees. But what do I set as the initial CRS to set up the correct transform? TIA


title: Zoning and land lot map in QGIS

text: Our office is using QGIS to design zoning maps. Our maps display the zones in colors (area shapefile with "Zone" text attributes) and the parcels with a thin outline and the number as label (area shapefile with "Parcel #" integer attribute). ![enter image description here](http://i.stack.imgur.com/3a433.png) What would be the best way to be able to output a list of all parcels with their area and their land use ? (note that some big land lots may be divided in different zones, so some parcel would appear twice, with let's say 35% of their area for residential zone and 65% for commercial zone) I thought of "splitting" the zone polygons with the land lots polygons (think of a puzzle-shaped cookie cutter). But I can't seem to find such a function. I thought of going the other way around, by performing a spatial join of the zone attributes into the parcel polygons. But this doesn't seem to work neither because it joins the first polygon it touches, and there seem to be no way to tell it to join the polygon which it overlaps the most. And idea ?! Thanks !! (and sorry for english mistakes ;) ) Olivier * * * Edit : Ok, UNION seems to be the function I need ! Thanks Rayer ! However, I tried it on my data, but it doesn't work at all. I get no error, but the progress bar doesn't move at all and I have to force quit QGIS. I tried on a very small dataset, and then it works well. My datasets consists of about 1000 zone polygons and 3000 land lots polygons. Is there a chance it would work with this amount of data ? Thanks ! * * * Still didn't find a way to do what I want... In fact, I didn't see there were some notifications in the console as soon as I run the UNION command : Message: "Feature geometry not imported (OGR error: )" Level : 0 The messages appears many many times (quickly fills the console), about 10 every second (when looking at the timestamp). At the end I get a shapefile with no content. Any idea about what's going wrong ? (by the way, is editing the initial message the right way to continue a discussion in this forum ? or should I rather answer to my question ?)


title: gdal_merge produces incomplete mosaic

text: I am using GDAL 1.91 on a Windows 7 machine. I have been mosaicing large numbers of rasters together and most of the time this process works exactly as expected. However a very few of my mosaics look incomplete. It is as though GDAL_merge simply ignored some of my input files (see image - notice the top lefthand corner): ![incomplete mosaic](http://i.stack.imgur.com/5QqTL.png) To try and eliminate this problem I have:   * checked that the missing rasters are valid and exist in my file system (correct on both counts - they display individually without error)   * checked my input string and that definitely includes the missing rasters.   * made all the input rasters the same type (e.g. Float32)   * explicitly specified the no_data value to ignore in the inputs   * **EDIT:** all the inputs are also the same resolution and running with verbose output gives no clues of any errors either The mosaic process involves joining a maximum of 395 rasters each of which is 980KB and the resultant files should be about 390MB. The process works very well most of the time, so this suggests that it is not a memory issue (I have 6GB of RAM and am using the 64 bit version of GDAL). The error is not isolated to only one mosaic. I have seen this elsewhere before and with GDAL 1.8 (which is why I upgraded to 1.91). In this particular case, the result using GDAL 1.8 was almost the exact opposite (i.e. the top left hand corner of this area did have data and the bottom left was empty - thus proving my data are valid but mystifying me further as to what is going wrong). Changing the order of my input files does sometimes make a difference but usually does not entirely solve the issue. Is this a bug (this is my gut feeling at the moment)? Otherwise, does anybody have any suggestions for what might be going wrong and/or work arounds? Many thanks! **UPDATE:** I have been studying my rasters in QGIS looking for differences and the ONLY difference I can see is that, although all the rasters were created in exactly the same way and are all the same resolution, same projection and all geotiffs, QGIS quotes the driver as VRT/Virtual Raster for the ones which gdal_merge is failing. QGIS correctly quotes the driver as GeoTiff for the rest!?! Gdalinfo on the other hand, quotes the driver as GTiff/GeoTiff for all the files. So - there is SOMETHING different about the tiles but what and why is not obvious... more digging to follow...


title: About zonal statistic function in ARCGIS

text: I calculated zonal statistic in Arcgis to get max/min value of each LC types. Broadly, LC image include forest, grassland, waterbody etc. But output max value of zonal statistic function give me 0.972 on the waterbody, 0.9908 on grassland, 0.9976 on forest. I think it is not reliable. What is the reason on this? Actually NDVI image has been masked on waterbody. Original value of the NDVI image range from -0.03 to 0.9976. Do you guys have another suggestion for this? Thanks a lot


title: How can models be integrated into Workflow Manager?

text: I have a ModelBuilder model that determines whether a field in a table contains null or duplicate values. The outputs of the model are two integer parameters, one with the count of nulls and one with the count of duplicates. I have a workflow in Workflow Manager that I would like to integrate the model into in order to implement logic branching. Is there a way in Workflow Manager to:   1. Run the model as a step in a workflow   2. Examine the output of the model and branch into 1 of 2 logical workflow branches. For 1) I assume that I must run the `JTXSteps.ExecuteGPTool` custom step, is that correct? For 2) (How) can I set the return code in the model? Must I refactor the model so that it only has one output parameter, or should I actually cause an error to occur for the return code to be non-zero? How can I use the return code to specify the logical branch to step into?


title: How to find the row that has a shape file containing a lat long point in my PostGIS geometry column

text: I am fairly new to PostGis and I'm trying to return a row in my table that has the shape file containing a lat long point I pass in via st_contains. Below is a query I wrote myself that brings back 0 rows with no errors as seen below.               =# SELECT * FROM es_zones      WHERE ST_Contains(geom, ST_GeomFromText('POINT(-73.952545 40.774576)'))=true;      Having found a query similar of that to which I am trying to accomplish here, below is my last query ran.               =# SELECT * FROM es_zones      WHERE ST_Contains(es_zones.geom,ST_transform(ST_GeometryFromText('POINT(-73.952545 40.774576)',26918),26918))=true      AND ST_Distance(es_zones.geom,ST_transform(ST_GeometryFromText('POINT(-73.952545 40.774576)',21918),21918))=false;      which also brings back 0 rows as seen below > > ----+-------------+-------+--------+--------------+-------------+--------------+-------------+---------------+---------------+ > gid | schooldist | boro | initials | creat_date | edit_date | zoned_dist | > boro_num | shape_leng | shape_area > > -----------+-----+----------+-------+-------- > > remarks | dbn | esid_no | label | geom > > (0 rows) Really not sure where I am going wrong. Any help is appreciated. I was able to get QGIS setup and connected to my PostGIS data. I don't understand what the 997487, 221261 coordinates are. How do they correspond to the lat lon points I got back from my GeoCoder Gem? I passed them in as my point but got the following error; ERROR: Operation on mixed SRID geometries ![QGIS view of PostGIS data](http://i.stack.imgur.com/YS5jq.png) Verifying the SRID shows               =# SELECT ST_SRID(geom) FROM es_zones LIMIT 1;      st_srid      ---------       926918      (1 row)      So it sounds like something may have happened to my spacial_ref_sys table. Or perhaps I didn't get it configured correctly to begin with. A Google search for my last error `ERROR: GetProj4StringSPI: Cannot find SRID (4326) in spatial_ref_sys` returned this StackOverflow question. Which points to running `spatial_ref_sys.sql` to regenerate the rows in that table. I found some elaboration on how to go about executing `spatial_ref_sys.sql` in the PostGIS documentation here. But I'm confused why it says it can't find SRID 4326 when looking in PGAdmin I can see it as shown below.![enter image description here](http://i.stack.imgur.com/WA7tG.png)


title: FME area on area overlayer dropping components of some features

text: As it says in the title, is this expected behaviour? What happens is that the dataset before the area on area overlayer has some features (in the image below, the large green parts are the same feature, there is a point where they touch). The lower part of the green feature is lost when the transformer is run. This shouldn't happen! I have an example dataset and workspace, here: https://www.dropbox.com/s/n50h39yftp7nhtx/fme_bug.zip. There are a number of other cases where this happens in my dataset. Any thoughts? ![This image shows, circled in blue, the component of the green area that is lost after running the aeraonareaoverlayer transformer.](http://i.stack.imgur.com/4tH8G.jpg)


title: Is it possible to save frequently used expressions in QGIS?

text: Is there any way to save frequently used expressions in the selection or field calculator dialogs? Obviously I can cut and paste them into a file where I might save them, but I mean within the context of QGIS like you can do in ArcMap.


title: Do pad functions of QGIS Field Calculator have bug?

text: I have this expression in the field calculator: `regexp_replace('1', '^(\\\d+)', rpad('\\\1', 3, '0') )` gives `01` rather than `001` as expected. the *pad functions work fine when called outside regexp_replace Weird! This has to be a bug? or am I missing something? If it is a bug I'll log something in the issue tracker. BTW I was delighted when I found I could use complex expressions including functions in the third parameter of `regexp_replace`! kudos to the QGIS developers!


title: How come "Zoom to Polygon" always puts me in the middle of Nebraska?

text: This is an issue related to this question I'd asked a few days ago. I've got that code working, to an extent. Every time, though, that I double click on one of the Polygons in the attribute table that I'd recently populated, ArcMap takes me to a point in the middle of Nebraska. ![attribute zoom](http://i.stack.imgur.com/Mdw7h.png) ![1:2,500 scale](http://i.stack.imgur.com/4OOOc.png) That's at 1:2500. If I zoom out, it just appears as an X in the middle of Nebraska. It's nearly the same for any of the Polygons I choose, except there are a different configuration of rectangles within that picture. But it's always in the middle of Nebraska. I can go into the attribute table using python to pull out the polygon information with the Polygon getPart() method, and it looks like this:               Out[132]: <Array [<Array [<Point (-123.070983887, 43.0001220703, #, #)>, <Point     (-123.070983887, 45.5001220703, #, #)>, <Point (-121.890991211, 45.5001220703, #     , #)>, <Point (-121.890991211, 43.0001220703, #, #)>, <Point (-123.070983887, 43     .0001220703, #, #)>]>]>      ...which, based on the points, should be in Oregon. But if I double click on it, it takes me to Nebraska. The other problem is that the polygons should be spread across the country (near national parks) but I think they're all getting bunched into that little square in Nebraska. Does anyone know why this might be happening? Is the array-within-array throwing it off? I'm adding the code that I used to create and fill the feature class. I created it with:               fcname = r'C:\Users\user\Documents\ArcGIS\parkpolys.shp'     output_fc = arcpy.CreateFeatureclass_management(os.path.dirname(fcname), os.path.basename(fcname), "Polygon")      And then filled it with:               def push_coords(coords):         """Push a set of corner coords to featureclass record"""         cursor = arcpy.InsertCursor(fcname)         feat = cursor.newRow()         print 'new row created'         for corner in corners:                 point.X = corner[0]             point.Y = corner[1]             array.add(point)         array.add(array.getObject(0))         polygon = arcpy.Polygon(array)         feat.shape = polygon         cursor.insertRow(feat)         array.removeAll()         del feat         del cursor              # Corners list format: ["lowerLeft", "lowerRight", "upperRight", "upperLeft"]        #new format: ["lowerLeft", "upperLeft", "upperRight", "lowerRight"]     cursor = arcpy.SearchCursor("gmaps")          for row in cursor:         corners = [[row.getValue('xmax'), row.getValue('ymin')],                 [row.getValue('xmax'), row.getValue('ymax')],                 [row.getValue('xmin'), row.getValue('ymax')],                 [row.getValue('xmin'), row.getValue('ymin')]]         push_coords(corners)


title: How do I drag item/file from ArcCatalog to custom UI form created with Toolbar Add-in?

text: I have created a button add-in in ArcMap that pops up custom UI/Form. I want to drag and drop file from ArcCatalog into the textbox on the Form. I am able to do so with windows explorer but not from ArcCatalog. I'm writing code using C#. When I drag the file from catalog, I am sure DataFormats type is not FileDrop. Any suggestions?               private void openShapeFileTextBox_DragEnter(object sender, DragEventArgs e)         {             if (e.Data.GetDataPresent(DataFormats.FileDrop))             {                 e.Effect = DragDropEffects.Copy; // Okay             }             else             {                 e.Effect = DragDropEffects.None; // Unknown data, ignore it             }             nameTextbox.Text = temp;         }


title: OpenLayers 3 and Interactions - Select and Draw

text: I'm trying to make small app with OpenLayers3. I want to use two interactions - select and draw. Unfortunately when selet adding it's style to draw and remove other (selected earlier) features from selection (more on screens below). How to prevent this? I want to have two independent tools.   1. Only select interaction ![Only select interaction](http://i.stack.imgur.com/Dr2xT.png)   2. Select and Draw ![Select and Draw](http://i.stack.imgur.com/oNwkW.png)   3. End drawing with select interaction ![On end drawing](http://i.stack.imgur.com/SnWKW.png)   4. Only draw ![Only draw](http://i.stack.imgur.com/tsp3w.png)   5. End drawing without select interaction ![End drawing without select interaction](http://i.stack.imgur.com/2L7Zf.png)


title: Converting rasters to polygons through cell centers

text: With arcgis, it is possible to convert a raster image to a polygon representing the boundary of the raster using the Raster To Polygon tool. The output polygon will follow the cell edges of the raster. I'd like to create a polygon which follows the centres of the raster's edge cells. Any idea how to do this? I'm agnostic as to whether it uses arcgis, postgis, python, C, whatever... The input can be an Esri 'Raster Dataset' or any gdal-supported format...


title: Bulk SHP/KML Processing Tools

text: I need to accomplish two things in **bulk** :   1. Convert SHP to KML   2. Modify the fill and outline color and opacity either of SHP or KML files **What tools can you recommend for these two tasks?** I've tried Google Earth Pro, which allows me to import in bulk and set the fill/outline color/opacity in bulk. But, I'm unable to export to KML in BULK. For the collection of data I have right now I have 193 *.shp files. In Google Earth Pro, I would have to export each one manually. In the future I will have 1000's of *.shp files and one-at-a-time won't cut it then. **Can you bulk export to KML in Google Earth Pro?**


title: arcmap export layer as a image with different scale

text: In our web application, we load too many layers cause the map loading very slow. I am thinking of export some static layers as images with different scale. There are totally 10 scales in our application, hence I want to export 10 images for a static layer (each scale 1 image). When the application load map data, I only load the image file for the static layer. When the map zoom, I will change different image according to the zooming scale. I am using arcmap 10.2. I feel it can be performed, but I haven't get the way to do it. Please advise me, thanks.


title: Create Addons for GRASS GIS

text: Is there a complete guide - example describing how to create an add-on/plugin for the GRASS GIS? I would like to extend the GUI with extra menus and toolbars creating extra functionality.


title: How To Show Map In bigg window after click on it by api for javascript

text: i want user when click on map open in full screen view as popup to show map bigger how to do this


title: Problem with variable for arcpy.select_analysis whereclause

text: I have a problem with the `whereclause` on python for `arcpy.Select_analysis`. I try to use a variable for condition:               for row_stat in rows_stat:               % for calculated the variable"maxi"         maxi= row_stat.MAX_GRID_CODE         % for calculated the variable"maxi"     del rows_stat                            % for calculated the variable"maxi"     where = '"GRID_CODE" = maxi'             % whereclause     arcpy.Select_analysis("point_acc.shp", "point_max.shp", where)      Ofcourse, I get an error message telling me that 'maxi' can't be found in "point_acc.shp" but I don't want 'maxi' but the number in the variable maxi! Who can help me write this?


title: Can a routing service take polygon geometry (for stops)?

text: I have built a small road featureclass, created a network off of it using network analyst, and created routes successfully in ArcMap. I'm now trying to create these routes in a silverlight web application and in order to do that I need to pass the geometry of two stops (polygons) to the routing service. I'm not having any luck and I think it's because the centroid being calculated by the web application isn't "exactly" on the routing service. So I'm wondering instead of passing a point centroid, I could instead pass the polygons geometry of my two stops to the routing service? Does anyone have an example of how to do this? Or is this possible?


title: python tool parameter: what datatype to list unique values of field?

text: I'm attempting to convert a 9.3 arcobjects tool to a 10.2 python tool. The form needs to allow user input of a shapefile, a field, and then the user needs to enter numbers corresponding to unique values of the selected field. Does anyone have any advice for prompting the display of the field's unique values? I imagine it would look something like the Reclassify tool form (but the unique values would be string and the list limited to ~10 values). (The input_shapefile and input_field are pretty straightforward with parameter datatypes 'Shapefile' and 'Field.') Many thanks!


title: PostGIS - ST_Difference query returns a dissolved result

text: I am using the ST_Difference to erase multipolygon feature by a mask. However, the query returns a dssolved multipolygon. The query is the following:               '''SELECT ST_AsGeoJSON(ST_Collect(geom))::json        FROM (SELECT (          ST_Dump (              ST_Difference(ST_GeomFromGeoJSON('%s'),ST_GeomFromGeoJSON('%s')))).geom)                  As f''' % (feature,mask)      Please take a look at those images to understand the problem: Input: ![enter image description here](http://i.stack.imgur.com/QO4MN.png) Dissolved output: ![enter image description here](http://i.stack.imgur.com/Slw32.png)


title: Python, GDAL, VRT: help wth code efficiency and processing time

text: I'm going to attempt to make this question as succinct as possible while still including relevant information, as I have trouble with being concise and am trying to avoid writing unnecessary stuff. Anyways, I need help with increasing the efficiency of my image processing python code. I'm relatively new to python/gdal/coding in general, and don't have very much knowledge about performance, processing time, efficiency, etc. -- so I need your help. However, I do have two specific ideas that could probably speed things up, if I knew how to implement them. In a nutshell, my code does the following:   1. Takes a set of classified Landsat geotiffs (say number of images in set = N), whose sizes all vary by 50-100 pixels, and extracts one element, water, from the classified image. It does this by utilizing GDAL to open and read the classified tiff as an array, and using the where function to produce an array of 1's and 0's (1's where pixel = water and 0's where pixel = not water). Currently, my code writes these arrays as tiffs using the same geospatial information as the original classified tiff, and then uses the set, as a whole, for input of the next step. However, it this step of writing the 'extracted' arrays as tiffs that takes the majority of processing time (8-9 of the 11-12 total minutes per individual N), and so I would like to try a different approach (see bottom) since these 'extracted tiffs' are only needed temporarily for step 2.   2. Next, I use a list pointing to these individual 'extracted' tiffs as the input for gdalbuildvrt, so that a VRT "stack" is built, encompassing the extent of all the individual files. This step is important for implementing step 3 properly, because of the individual file's varying sizes. It is also the main reason why I take the time to write out the 'extracted' arrays as tiffs-- to preserve their spatial information so that gdalbuildvrt will accurately stack each image on top of each other according to their spatial reference. [NOTE: When calling gdalbuildvrt, I use a mask so that the "empty" edges (as in, the empty space around the smaller tiffs) of the individual tiffs are filled with 0's]   3. After creating the N-layer virtual stack of the individual images, and thus 'normalizing' their sizes accurately according to their spatial info, I use a for loop (done since GDAL can only read one band at a time) to read out each resized 'extracted tiff' as a layer of a 3-d numpy array [this array would be of shape (N files, Nlines, Nsamples) if the overall size of the stack is Nlines by Nsamples]. Then, using numpy.sum, I get a 2-d array of size (Nlines, Nsamples), that represents the per-pixel total number of water observations seen in the N files over the course of the set. It then writes this total array as a tiff, which is the final output that I am interested in. Now knowing my overall goal, hopefully you understand my code explanation. I will attach it if it will be helpful, but it's fairly long and has a lot of comments so I'm leaving it out for now. There are two ways I have thought about that could (or could not) help improve the efficiency of my code:   1. Skip the step of writing out the 'extracted' arrays as tiffs and either feed them directly into gdalbuildvrt after extracting (while still retaining the spatial information, which it absolutely has to do) or write each file as some other virtual/in-memory format that won't take as long as writing tiffs, which can then be fed as inputs into gdalbuildvrt. I've read that you can't write an array as a VRT, so if that's true, are there any suggestions for alternative methods?   2. I was also thinking of rearranging the code to do this: Feed the classified tiffs directly into gdalbuildvrt which will resize them like it already does, read each classified band one by one into the large 3-d numpy array (in the loop i described above), THEN using the where function to extract water and not water. After that I would total the numpy array, and then write that result to tiff, and be done. The only issue I'm afraid of with this method is running into memory problems when trying to do the where function over a large 3-d numpy array. I imagine the volume of this 3-d numpy array might be the same volume as the numpy arrays in the original method, but I don't know how having values ranging 1-7 (each class) instead of 0's and 1's will affect performance. I also don't know how numpy's where function will act in terms of performance over the large 3-d array. Can anyone shed some light on what might be the best improvement, in terms of performance and/or efficiency? I realize this is a long post, but I figure being more explicit is better than not enough detail. Notes: I'm using python and gdal for the majority of my code, and cannot resort to arcpy or any other module that requires a license. Also, in case you are wondering, I am going to be processing thousands of 'sets' (~4050) with N ranging anywhere from 10 to 30, so I really need to find ways to improve efficiency. Right now, 11-12 minutes per n puts processing time at a few years.... so, I need your help! Thanks so much!


title: arcpy geoprocessing service with sde database produces no output url

text: I'm trying to use a geoprocessing script to clip some data and return to the user (via url) the path to a zip file with the results of the clip. I also need to connect to an sde database. I can get a url to show up if I use shapefiles and but it seems to loose the plot when i use the sde connection (probably because i have to set the workspace to the sde connection and then the server doesn't have a place to write temp files to?!). It runs okay but doesn't produce any output on the server folders and i do not get a url. I have pasted the code below (hopefully it shows up ok). Thanks. ps: I copied a lot of the code from esri's ship and clip example. Also i registered "D:\arcgisserver\directories\service_output" on the server as a folder so it is able to see it.               import arcpy     import sys     import os     import traceback     import zipfile          arcpy.env.overwriteOutput = True          def zipUpFolder(folder, outZipFile):         # zip the data         try:             zip = zipfile.ZipFile(outZipFile, 'w', zipfile.ZIP_DEFLATED)             zipws(unicode(folder), zip, "CONTENTS_ONLY")             zip.close()         except RuntimeError:             # Delete zip file if exists             if os.path.exists(outZipFile):                 os.unlink(outZipFile)             zip = zipfile.ZipFile(outZipFile, 'w', zipfile.ZIP_STORED)             zipws(unicode(folder), zip, "CONTENTS_ONLY")             zip.close()             #Message"  Unable to compress zip file contents."             gp.AddWarning(get_ID_message(86133))     def zipws(path, zip, keep):         path = os.path.normpath(path)         # os.walk visits every subdirectory, returning a 3-tuple         #  of directory name, subdirectories in it, and filenames         #  in it.         for (dirpath, dirnames, filenames) in os.walk(path):             # Iterate over every filename             for file in filenames:                 # Ignore .lock files                 if not file.endswith('.lock'):                     #gp.AddMessage("Adding %s..." % os.path.join(path, dirpath, file))                     try:                         if keep:                             zip.write(os.path.join(dirpath, file),                             os.path.join(os.path.basename(path), os.path.join(dirpath, file)[len(path)+len(os.sep):]))                         else:                             zip.write(os.path.join(dirpath, file),                             os.path.join(dirpath[len(path):], file))                          except Exception as e:                         #Message "    Error adding %s: %s"                         arcpy.AddWarning(get_ID_message(86134) % (file, e[0]))         return None               if __name__ == '__main__':         try:             arcpy.env.workspace = os.path.join(r"D:\arcgisserver\directories\service_scripts", u'db_Conn.sde')             ws = "D:\\arcgisserver\\directories\\service_output"              arcpy.env.scratchWorkspace = "D:\\arcgisserver\\directories\\service_output"             outputZipFile = "D:\\arcgisserver\\directories\\service_output\\out.zip"             zipFolderPath = "D:\\arcgisserver\\directories\\service_output\\output"                  aoi = "455719.332353, 6300257.134719, 499534.419983,6269598.547920"             serviceDataset="ADMINISTRATIVE_AREA_HISTORICAL"                  dsList = arcpy.ListDatasets("*","Feature")             for ds in dsList:                 if serviceDataset in ds:                     fcList = arcpy.ListFeatureClasses("*","",ds)                     for fc in fcList:                         featureclass = fc.split(".")[2]                          arcpy.AddMessage(fc)                              exts = []                         exts = aoi.split(",")                         extent = arcpy.Extent(exts[0],exts[1],exts[2],exts[3])                         pts = [arcpy.Point(extent.XMin, extent.YMin),                         arcpy.Point(extent.XMax, extent.YMin),                         arcpy.Point(extent.XMax, extent.YMax),                         arcpy.Point(extent.XMin, extent.YMax)]                         array = arcpy.Array(items=pts)                         poly = arcpy.Polygon(array)                              sr = arcpy.SpatialReference()                         sr.factoryCode = 26912                         sr.create()                         poly = arcpy.Polygon(array,sr)                                tmp_poly = os.path.join(ws, "polygon.shp")                         tmp_orig = os.path.join(ws, "to_clip.shp")                         outputpath = os.path.join(zipFolderPath, fc.split(".")[2]+".shp")                              arcpy.CopyFeatures_management(poly, "D:\\arcgisserver\\directories\\service_output\\tmp_poly.shp")#tmp_poly)                         arcpy.Clip_analysis(fc, tmp_poly, outputpath)                  zipUpFolder(zipFolderPath, outputZipFile)                  arcpy.SetParameterAsText(0, outputZipFile)             arcpy.SetParameterAsText(1, arcpy.env.scratchWorkspace)               except:             tb = sys.exc_info()[2]             tbinfo = traceback.format_tb(tb)[0]             pymsg = "ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n    " + \                     unicode(sys.exc_type)+ ": " + unicode(sys.exc_value) + "\n"


title: save IGeometry to disk as a Layer file

text: How can I save an IGeometry to disk as a Layer file? ArcEngine 10, C#, VS2010 Edit 1: With Jason Scheirer 's suggestion I managed to convert a point into a shapefile - and I'll post the code here for the community.               string strFolder = @"C:\temp";     string strName = "NewShapeFile"; //don't include .shp extension     string strShapeFieldName = "Shape";     ShapefileWorkspaceFactory pWorkspaceFactory = new ShapefileWorkspaceFactoryClass();     IFeatureWorkspace pFWS = pWorkspaceFactory.OpenFromFile(strFolder, 0) as IFeatureWorkspace;      IFields pFields = new FieldsClass();     IFieldsEdit pFieldsEdit = pFields as IFieldsEdit;     IField pField = new FieldClass();     IFieldEdit pFieldEdit = pField as IFieldEdit;     pFieldEdit.Name_2 = strShapeFieldName;  //Name is read-only Name_2 is writeable -- don't ask me why     pFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;      IGeometryDef pGeomDef = new GeometryDefClass();     IGeometryDefEdit pGeomDefEdit = pGeomDef as IGeometryDefEdit;     pGeomDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPoint;     pGeomDefEdit.SpatialReference_2 = CrowWingCoordinates;     pFieldEdit.GeometryDef_2 = pGeomDef;     pFieldsEdit.AddField(pField);     pField = new FieldClass();     pFieldEdit = pField as IFieldEdit;     pFieldEdit.Length_2 = 25;     pFieldEdit.Name_2 = "TextField";     pFieldEdit.Type_2 = esriFieldType.esriFieldTypeString;     pFieldsEdit.AddField(pField);     IFeatureClass pFeatClass = pFWS.CreateFeatureClass(strName, pFields, null, null, esriFeatureType.esriFTSimple, strShapeFieldName, "");     IFeatureBuffer fb = pFeatClass.CreateFeatureBuffer();     IFeature feature = pFeatClass.CreateFeature();      IFeatureCursor featureCursor = pFeatClass.Insert(true);     double x = -94.105797;      double y = 46.559214;      IPoint p = new PointClass();      p.SpatialReference = gcs;     p.PutCoords(x, y);      p.Project(myCoordinates);      IGeometry geometry = p as IGeometry;     fb.Shape = p;     featureCursor.InsertFeature(fb);     featureCursor.Flush();


title: Create Shapefile with auto increment primary key in QGIS

text: Is there a way to create an auto increment primary key column when creating a new Shapefile in QGIS?


title: How to edit symbol diameter related to the catchment area?

text: I have a table with column Catchment area (in km2), and I need to set symbol diameter based on that column. I need to get .sld file so I'm working with uDig but if you know better software for this, I will accept your advice. EDIT: I found out how to do it, but result is not acceptable because symbols are too big. There is no any normalization option in uDig or QGIS. Thank you.


title: Trouble with PostGIS and multiple selection sources

text: I'm attempting to take geometry columns from two different tables, perform basic PostGIS operations on them, and have them returned in a single result set. My SQL statement looks like the following:               SELECT (ST_Dump(ST_Polygonize(exterior_rings.the_geom))).geom AS exterior,            (ST_Dump(ST_Polygonize(interior_rings.the_geom))).geom AS interior          FROM (SELECT the_geom FROM geodata."_unioned_exterior_rings") AS exterior_rings,           (SELECT the_geom FROM geodata."_unioned_interior_rings") AS interior_rings;      When I execute this command in pgAdmin III's arbitrary SQL query tool, after roughly 8 minutes of execution, I receive:               ********** Error **********      ... and the query tool complains it is no longer connected to the database. This leads me to believe the operation is timing out. The respective sizes of the tables are fairly humble. `_unioned_exterior_rings` has ~4000 rows, and `_unioned_interior_rings` has ~100. I've had other queries succeed after twice as much run time so I don't feel as if I'm hitting some preset limit. When I use pgAdmin's console tool to execute the query, after some time it becomes apparent that the entire tool has been disconnected without any kind of error message. On some reconnect attempts I've received notices that the reconnection fails because the database is in a fatal recover mode. If I try again shortly there after it usually works. When I execute the statements separately (doing the polygonize/dump on an individual table) I get the expected results after only a few seconds. The version of PostgreSQL is 8.4, PostGIS is 1.4.2, and pgAdmin III is at 1.12.2. I can't claim any real SQL proficiency, so I'm willing accept that my mistake may be simple or fundamental. Any help would be appreciated.


title: IApplication reference with ArcObjects Class Extension

text: I have an ArcObjects (9.3) Class Extension which has an event handler for OnCreate so that feature properties can be manipulated / assigned when a feature is first created. Some of the functionality needs enclosing in an edit operation, which requires an IWorkspaceEdit object. As far as I know the only way I can get the IWorkspaceEdit is by getting the IEditor (extension) from IApplication, then grabbing the IEditor's EditWorkspace. In an ArcObjects command (extending BaseCommand) I get the IApplication 'hook' passed on create (that is, the command's creation), so I can retain a reference to it. However, in the class extension's OnCreate method I am only passed the IObject which has been created. I can't see a way to get the IObject's workspace to work from that direction. Can anybody please tell me   1. How to access the calling IApplication from a class extension or,   2. How to get an IObject's EditWorkspace or,   3. Any other way to create an edit operation within a class extension's OnCreate method? Any help very much appreciated.


title: Use Dojo to load modularized widget config files

text: I've taken a basic webmap 'ESRI HomeButton' and I'm trying to create a modularized webmap setup from it. I want to break down this simple webmap to have it utilizing a module system with config files and more. Then I can take other single-use webmaps and breakout the widgets into separate files much more easily and approach everything very much like the FlexViewer environment. The problem I am having is that in my widgetHomeButton.js file, I took the HomeButton.js code and added in three lines... In my require, I added "widgets/homebutton/config" to point to my config.js file. I should mention I have a package.js file that sets up the folders so that 'widgets' points to package_path + "path/to/widgets" just like in the example. In my function, I added **widgetConfig** as the variable to store the config.js file contents, which happen to be JSON style data. Then in the HomeButton widget code I added a single line of code that reads...               console.log( widgetConfig );      And all I get are errors. Any ideas what I am doing wrong?


title: ESRI : 54019 or World winkel II in QGIS

text: Currently I didn't find World winkel II or ESRI : 54019 coordinate system support in qgis. One of our client using this.. Can any one tell me a solution to load shape files,which are using his coordinate system. I didn't find .proj4 file for this.


title: What are the coordinates of the mysterious hole in Siberia?

text: I am interested in taking a Google Earth aerial view of the mysterious hole recently discovered in Siberia and perhaps doing some hydrological analysis of the area in a GIS. However, there appears to be very little information on the actual coordinates of the hole. What are the coordinates of this interesting geological feature? More details on the Siberian hole:   * Siberian Times   * NBC ![enter image description here](http://i.stack.imgur.com/Etspy.jpg) Source: Associated Press


title: ArcGIS for Javascript Template available?

text: ESRI has a stack of templates in the ArcGIS Online section. But when you download them it turns out that in order to make them run, you have to first create a webmap via their WebMap Viewer application. And even then the map that you then host only references the webmap application you create on the ESRI Portal. Coming from a FlexViewer background where I could use the Viewer app to create an application and then all files and folders were generated locally and I could edit or tweak as needed I was hoping for something similar with the WebMap Viewer. I am looking for a blank webmap using the ArcGIS API for Javascript complete with widgets, place to insert basemaps, operational maps, etc. Does such a map template exist?


title: Handling Popups in Leaflet with Overlapping Features

text: I am new to Leaflet as well as Javascript, and I'm looking for some advice on a basic topic. Currently, I am trying to create a spatial index of library holdings for old topographic maps that people can access and download online. The source of the features is a GeoJSON file. Sample Map The problem is that the school library has multiple types of the same map over the years, but when I click on a map, only one popup shows up. I want the option to cycle through multiple maps, but I am unsure on how to approach the problem. Is there a special term in the local jargon for cycling through overlapping polygons, or is there a stronger approach to this problem? Any advice on the direction I should take is appreciated. Thanks!


title: Projection mismatch when plotting precipitation with ggplot

text: I'm trying to plot accumulated precipitation in Brazil for a given period, see code below: **Load packages**               library(raincpc)     library(SDMTools)     library(raster)     library(ggplot2)     library(rgdal)     library(maptools)     library(rgeos)     library(mapproj)      **Use _cpcrain package_ to obtain data**               cpc_get_rawdata(2014, 1, 01, 2014, 1, 08)          rain1 <- cpc_read_rawdata(2014, 1, 01)     rain2 <- cpc_read_rawdata(2014, 1, 02)     rain3 <- cpc_read_rawdata(2014, 1, 03)     rain4 <- cpc_read_rawdata(2014, 1, 04)     rain5 <- cpc_read_rawdata(2014, 1, 05)          rain_tot <- rain1 + rain2 + rain3 + rain4 + rain5      **Zoom in Brazil**               lon_vals <- seq(275, 330     , 0.5)     lat_vals <- seq(-60, 20      , 0.5)     reg_box <- expand.grid(lon = lon_vals, lat = lat_vals,     stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)     reg_box$rain <- extract.data(reg_box, rain_tot)     reg_box$rain_chunks <- cut(reg_box$rain, breaks = c(0, 80, 160, 240, 320), include.lowest = TRUE)      **Use ggplot for visualization**               gfx_gg <-        ggplot(data = reg_box) +                   geom_raster(aes(lon, lat, fill = rain_chunks)) +                   scale_fill_manual(values = c("lightgrey", "red", "green", "blue")) +                   theme(axis.text = element_blank()) + # , axis.ticks = element_blank()                    labs( x = NULL, y = NULL, fill = " Precipitation (mm)") +                   ggtitle("Rainfall over South America for selected period")      which, as expected, yields: ![enter image description here](http://i.stack.imgur.com/y9sHU.png) **Add state boundaries**               brazil1 <- getData("GADM", country = "BRA", level = 1)          gfx_gg +       geom_polygon(data = brazil1, aes(x = long, y = lat, group = group),                    color = "black", show_guide = FALSE, fill=NA) +       labs(x = "x", y = "y") +       theme_bw()      But I'm getting an unexpected error: ![enter image description here](http://i.stack.imgur.com/QfSNx.png) How can I fix this (projection?) error?


title: How to create polygon feature class from UTM coordinates?

text: I am working to create a polygon from UTM coordinates. The script I am using fails with a "must be a sequence, not float" error:               Traceback (most recent call last):       File "F:/sample/clipTiles.py", line 104, in <module>         arcpy.Array([arcpy.Point(*coords) for coords in feature])))     TypeError: type object argument after * must be a sequence, not float      I have been working off the ESRI Polygon help section. My UTM coordinates are generated earlier in a script and take the following form:               [[734855.5142000001, 4744379.424500001], [740607.6685000001, 4744379.424500001], [740607.6685000001, 4736862.1753], [734855.5142000001, 4736862.1753]]      These are the min/max point pairs of a custom bounding box I need to convert to polygon. * * *                               # feature_info eventually takes the following form: [[734855.5142000001, 4744379.424500001], [740607.6685000001, 4744379.424500001], [740607.6685000001, 4736862.1753], [734855.5142000001, 4736862.1753]]                     feature_info = [[xmin, ymax], [xmax,ymax], [xmax, ymin], [xmin, ymin]]                          # A list that will hold each of the Polygon objects                     features = []                          for feature in feature_info:                         # Create a Polygon object based on the array of points                         # Append to the list of Polygon objects                         features.append(                             arcpy.Polygon(                                 arcpy.Array([arcpy.Point(*coords) for coords in feature])))                          # Persist a copy of the Polyline objects using CopyFeatures                     outshp = os.path.join(outws, row[2] + "_shp")                     arcpy.CopyFeatures_management(features, outshp)      How can I write these UTM coordinates to a polygon featureclass?


title: How to customize inputs in Iterate Feature Class - ModelBuilder 10.1

text: I am using the Iterate Feature Class tool in ModelBuilder to iterate through feature classes within a Geodatabase. I would like to add the capability to custom select which feature classes within a single GDB to iterate through. This desired capability is very similar to the "Input Datasets" step of Merge (Data Managment). ![enter image description here](http://i.stack.imgur.com/KMZUg.jpg) The Iterate Feature Class tool only provides you with 3 ways to filter: -Feature Type (line, polygon, etc)   -Wildcard   -Recursive  The ability to custom select which feature classes to iterate does not exist, can someone please help me figure out how to customize the input? Thanks!


title: OpenLayers WFS - GetFeature Draw Circle

text: I am developing GIS Portal using OpenLayers, I use WFS GetFeature to retrive some datas like bellow:               var layer = map.getLayersByName("pf_trecho_at");          var polygonLayer = new OpenLayers.Layer.Vector("Polygon Layer");          var protocol = new OpenLayers.Protocol.WFS.fromWMSLayer(layer[0], {                         url: "/geoserver/wfs",                         featurePrefix: 'light',                         readFormat: new OpenLayers.Format.GeoJSON(),                         outputFormat: "JSON",                         extractAttributes: true                     });     control = new OpenLayers.Control.GetFeature({                                    protocol: protocol,                 box: true,                 hover: true,                 multipleKey: "shiftKey",                 toggleKey: "ctrlKey"             });          control.events.register("featureselected", this, function(e) {                                   polygonLayer.addFeatures([e.feature]);                      this.displayPopup(e, 'pf_trecho_at');             });             control.events.register("featureunselected", this, function(e) {                 polygonLayer.removeFeatures([e.feature]);             });               map.addLayer(polygonLayer);     map.addControl(control);     control.activate();      My doubt is, the selection control draw a rectangle to select the features, is it possible do draw a circle to select then?


title: How to render a multi-layered (raster + vector) map to a file (PDF, SVG or JPEG, PNG)?

text: When using OpenLayers I can compose a map with several layers which is displayed nicely by my browser, but I need to do the same server side (no browser composing the layers), and I'm waiting a image file as output (ie PDF, SVG, JPEG or PNG). I need to do composition with at least OSM and Vector layer types. What are my solutions ?


title: Algorithm for splitting polygons/polylines into multiple sections at 180° E/W

text: I have code (in C++) that draws various polygons on a map that uses a cylindrical projection (typically Mercator, sometimes Miller Equidistant if you want to be precise). For simple polygons, curves are drawn between vertices through a set of tiepoints such that the line roughly forms either the great circle segment or rhumb line between the two vertices. The more complex polygons also include arcs which generally don't form part of a great circle or rhumb line. At the time of drawing, the tiepoints that make up the polygons are translated into screen positions and drawn using built-in Windows C++ drawing commands. For the most part, all my polygons and polylines are drawing fine, but I'm encountering problems when at least one line crosses the 180° meridian. My display is effectively limited to 180° W to 180° E. One solution I've considered is to unwrap the points that make up the polygon, perform some mathematical wizardry to calculate the screen position for points outside of [-180, 180], and then draw the object twice, once for 180° W and then again for 180° E. Unfortunately this method falls apart if there are 2 consecutive crossings in the same direction. As unlikely as that is to happen in my case, the mere possibility has me worried. It also doesn't work for objects that enclose one of the poles. Is there an algorithm for splitting a polygon along the 180° meridian while also properly closing it? Bonus points if it also works for polygons that enclose one of the poles. Extra bonus points if it'll also work for polar and/or conic projections (not currently supported, but may be required in a later project).


title: Programmatically load composer from template and generate atlas using pyQgis

text: I Am trying to build a plugin to load a print composer from file, generate an atlas and export to image. So far I have been successful in loading the template and exporting it to image. I have been unable to add any of the layers in the legend (which are also in the toc) to the exported map, which results in a blank map and none of the field expressions working.                   # Get layers in the legend and append, must be a cleaner way to do this?         layers = self.iface.legendInterface().layers()         layerStringList = []         for layer in layers:             layerID = layer.id()             layerStringList.append(layerID)              # Add layer to map render         myMapRenderer = QgsMapRenderer()         myMapRenderer.setLayerSet(layerStringList)         myMapRenderer.setProjectionsEnabled(False)              # Load template         myComposition = QgsComposition(myMapRenderer)         myFile = os.path.join(os.path.dirname(__file__), 'MMR_Template.qpt')         myTemplateFile = file(myFile, 'rt')         myTemplateContent = myTemplateFile.read()         myTemplateFile.close()         myDocument = QDomDocument()         myDocument.setContent(myTemplateContent)         myComposition.loadFromTemplate(myDocument)              # Save image         myImagePath = os.path.join(os.path.dirname(__file__), 'come_on.png')         myImage = myComposition.printPageAsRaster(0)         myImage.save(myImagePath)      Here is a snippet from the loaded template which should setup the atlas:                <Atlas hideCoverage="false" featureFilter="reference = '61922'"    coverageLayer="desktop_search20130615160118593" fixedScale="true" composerMap="0" singleFile="false" filenamePattern=""reference"" enabled="true" filterFeatures="true" sortFeatures="true" sortKey="0" sortAscending="true" margin="1"/>      I am also unsure of the best way to add all the layer in the toc to the instance of QgsMapRenderer(). Thanks


title: Where is the pdf of the US Census data stored that contains all the codes?

text: I'm working with US census data and I can't locate the documentation that explains all of the codes (e.g. P01001 etc) that is used to describe the attribute tables in the shapefiles. I haven't been able to find it on the www.census.gov site and it's tricky to search for if you forget the exact name (which I have!). Thanks.


title: Create python script to update field based on change in another field ArcGIS

text: I am trying to update a field "TPN" based on another field "PASS_NUM". In my script "TPN" is a newly added number field and thus blank, "PASS_NUM" is also a number field but it contains the wrong numbers. I want to assign values to "TPN" when "PASS_NUM" is no longer equal the the value in the row above it. The feature class table is sorted by "PASS_NUM". See example of attribute table I want to achieve here: ![enter image description here](http://i.stack.imgur.com/WrwmX.jpg) Is there any way to check if the value in the row above is equal to the next value (some sort of onchange event)? The PASS_NUM field does not increment by 1 each time (i.e. it may skip 7 as in the example) but I want to TPN to only go up one. Is this possible using the UpdateCursor in arcpy? Or is there another way to achieve this through scripting? Thanks Much!


title: How to Calculate fields for attribute table using Model Builder

text: I am working with the Model Builder tool in ArcMap. I have created a simple model that takes the Strip Map Index Features as an input. I have added fields to the model to produce more fields in the attribute table. When I run my model, it works as the original Strip Map Index Features tool should work, however, in the attributes table, I see my fields that I added, but they are not being populated / calculated. I am not sure how to calculate / populate these fields. I have searched around the net for solutions but am not getting too may concrete answers as I am still pretty new to using model builder. Any help on this issue would be greatly appreciated. Thanks in advance. =========================== My apologies as I have not updated / seen this post in a while. I have in fact found a solution to my problem. To reiterate the problem,... I was using model builder to create a tool for arcmap. I was basing my tool off of the Strip Map Index Features tool in arcmap. I wanted to add extra fields to the attribute table. When my tool ran, I would look in the attribute table and saw that the fields were there, however, they weren't being populated with any values. In my model, I had added the 'Add Field' tools but I hadn't added any 'Calculate Field' tools. Once I attached the 'Calculate Field' tools to my 'Add Field' tools, I was then able to open (double-click) the 'Calculate Field' tool to enter in a formula/expression/function. This helped add values to the fields in my attribute table. Again, for those concerned with this post, sorry I hadn't responded to it sooner.


title: Loading not advertised layers into opengeo SDK app

text: I`ve made several apps using OpenGeo suite SDK but so far never needed to keep layers from being exposed in the capabilities. Now I´ll like to "hide" some of them but still being able to load them. I taught that the solution was to set the layers as not advertised but if I do that the app will not load those layers!! As far as I understand I should be able to call a layer that is not advertised!! Using geoexplorer and saving the map I can do the same with no problems…. Any idea what may be wrong?


title: Remove Layer preview in Geoserver

text: I'm using Geoserver 2.1.1, Is there way to remove the Layer preview?, actually I'm only using the gmaps web services, I want to remove the KML preview or at least add some sort of authentication.


title: Controlflow for Geoserver

text: How to make use of controlflow.properties file for geoserver. I tried pasting it in \geoserver\WEB-INF\lib... But I couldn't feel any difference whether it is working or not ?


title: How to use navteq information for turn restrictions in pgrouting?

text: I'm developing a pgrouting application using navteq information but I don't know how to use the navteq information to add turn restriction in pgrouting. Even, I don't know where to find this kind of information in the navteq tables. I guess this data would be in condition/Driving Manoeuvres. Also, I'm tried to figure out how to use the divider information to add some turn restrictions in pgrouting but I guess that must be another source of information in this tables (navteq) that can be more easy and useful for turn restrictions.


title: WFS-T Edit using OpenLayers, Geoserver, PostgreSQL  ("error performing insert")

text: I tried following the tutorials: http://www.gistutor.com/openlayers/22-advanced-openlayers- tutorials/47-openlayers-wfs-t-using-a-geoserver-hosted-postgis-layer.html and: http://apps.who.int/tools/geoserver/www/wfs-t.html With the following results:   * the delete-function works, ergo I am able to alter my data in the postgres file.   * when inserting or editing a polygon I always get the answer: " Error performing insert: Error inserting features" or " Update error: Error occured updating features" I made a simple file, using only one kind of projection (EPSG:32647)- still, same problem. Where do I go wrong?! My Postgres Data-Table:   * gid - Integer, not null   * fid - double precision   * geom - geometry(MultiPolygon,32647) Here is my code: function init(){               //Save-strategy     var saveStrategy = new OpenLayers.Strategy.Save();          //empty map, bounds are test-layer bounds (EPSG:32647)       map = new OpenLayers.Map({         div: "map",         allOverlays: true,         maxExtent: new OpenLayers.Bounds(         653237.69439077,1519879.063165,655229.57939001,1520825.6733868         )         });          //WFS-Layer Test= editable data     var test = new OpenLayers.Layer.Vector("Editable Features", {         strategies: [new OpenLayers.Strategy.Fixed(), saveStrategy],             protocol: new OpenLayers.Protocol.WFS({         url: "http://..../wfs",         featurePrefix: 'testkf',         featureNS: "http://.../testkf",         featureType: "test",         geometryName: "geom",         })     });          map.addLayer(test);          //Toolbar:         var panel = new OpenLayers.Control.Panel(         {'displayClass': 'customEditingToolbar'}         );              var navigate = new OpenLayers.Control.Navigation({         title: "Pan Map"         });              var draw = new OpenLayers.Control.DrawFeature(         test, OpenLayers.Handler.Polygon,         {         title: "Draw Feature",         displayClass: "olControlDrawFeaturePolygon",         multi: true         }         );              var edit = new OpenLayers.Control.ModifyFeature(test, {         title: "Modify Feature",         displayClass: "olControlModifyFeature"         });              var save = new OpenLayers.Control.Button({         title: "Save Changes",         trigger: function() {         if(edit.feature) {         edit.selectControl.unselectAll();         }         saveStrategy.save();         },         displayClass: "olControlSaveFeatures"         });               panel.addControls([save, edit, draw, navigate]);     panel.defaultControl = navigate;     map.addControl(panel);       map.addControl(new OpenLayers.Control.LayerSwitcher());     map.zoomToMaxExtent();      } Thanks for your help!! Edit: I use Geoserver 2.1.4 and PostGIS 2.0 for PostgreSQL 9.1


title: Opening a featureclass not being the owner of the data in SDE

text: How can you open a `featureclass` if you are not the schema owner but have been given privs to select or view the data? So given a name of a `featureclass`, say "STREETS" - how can I open this if I am using an SDE connection with username "VIEWER" (given privs to see "EDITOR" owned data) but the owner of "STREETS" is "EDITOR"? - using the following code:               IFeatureWorkspace.OpenFeatureClass(string name)      This will work if you are the owner, you just put:               IFeatureWorkspace.OpenFeatureClass("STREETS")       But if you are not the owner then you have to fully qualify it, a la:               IFeatureWorkspace.OpenFeatureClass("EDITOR.STREETS")      Is there anyway round this? Bearing in mind my users will not have access to an sde connection for the data owner.


title: GeoTools 10 - Convert Raster To Vector

text: I am fairly green with GeoTools, and having a hard time converting a raster (GridCoverage2D) into a vector format. In the hours I've spent searching this site and Google, there seems to be multiple ways this could be done: PolygonExtractionProcess, RasterToVectorFactory, RasterToVectorProcess to name a few. However, regardless of what JAR's I reference in my project, I can't seem to find any of these classes. So, I am hoping someone can tell me the easiest way to perform raster to vector with GeoTools 10.


title: Moving point on OpenStreetMap route

text: Basically I generate random numbers in R language which are indeed times between bus arrivals. I need to simulate a moving point(a bus) in OpenStreetMap by using this generated numbers. I have read that R language communicates very good with OSM but not sure how to realize this goal. Any idea?


title: ESRI geodatabase, Oracle spatial performance

text: We are storing data in a Geodatabase on Oracle 11g using Oracle Spatial with SDO_Geometry, and then consuming that using ArcGIS server 10.2.1. The performance is very slow, because we are talking about between one to two million records. When I loaded the data, I didn't pay attention to creating or configuring indexes, as I just used FME to load the data from a file geodatabase. What are the best practices when you use Geodatabase on Oracle with SDO geometry? Which indexes are being used, the oracle ones, or sde ones? I am very new to database performance issues because I am mostly a developer. Any help, best practices, guidance, will be appreciated. Thank you.


title: Versions of Python & QGIS for new user to both to install?

text: I have used GIS (very basic functions), but now want to be able to customize to my needs. My initial research leads me to believe that Python + QGIS is the way to go. I understand that there are several plugins that will be useful along the way, but do not have enough knowledge to discern which ones and if they are all compatible with Python V3. The install will be on Windows 7 machine. Since I am learning Python and QGIS for the first time I want to use the newest versions possible-provided they (and their respective plug-ins) are compatible. What versions of Python and QGIS should I, as a new user to both, install?


title: Copying raster from Raster Catalog

text: I have a Raster Catalog containing 4 File GeoDatabase Raster Datasets. Each raster is about 150-250 GB. I want to convert these to rasters with .TIFF extension so I can publish it in a map service. The .gdb however, that contains the rasters, are strangely enough around 65 MB. How can this be? I have tried using Copy Raster tool, and also exported directly from TOC. Both of these are taking ages to even get started.


title: Problem with mosaic orthophoto

text: ![This should be a better mosaic :\)](http://i.stack.imgur.com/afnZ1.jpg) I am trying mosaic my orthophotos with ArcMap and Qgis as well. But every times the intersected area will be white color. I know its something about this frame around my pics but normally it is ignored. Have somebody idea or same problem before?


title: Routing result not showed up

text: I'm following the http://workshop.pgrouting.org/ but my routing result not showed up at all, like the picture below: ![enter image description here](http://i.stack.imgur.com/Wk8EG.jpg) my firebug also didn't give me an error, I'm using geoserver 2.4.2, postgresql 9.2 x64, postgis 2.1.1, and pgrouting 2.0. Here are the code that originally came from the pgrouting latest workshop:               <!DOCTYPE html>      ol3 pgRouting client #ol-map { width: 100%; height: 400px; } start final               <button id="clear">clear</button>          <script type="text/javascript">            var map = new ol.Map({         target: 'ol-map',         renderer: ol.RendererHint.CANVAS,         layers: [           new ol.layer.TileLayer({             source: new ol.source.OSM()           })         ],         view: new ol.View2D({           center: [-127998, 6974591],           zoom: 10         })       });            var params = {         LAYERS: 'pgrouting:pgrouting',         FORMAT: 'image/png'       };            var startPoint = new ol.Overlay({         map: map,         element: document.getElementById('start-point')       });       var finalPoint = new ol.Overlay({         map: map,         element: document.getElementById('final-point')       });            var transform = ol.proj.getTransform('EPSG:3857', 'EPSG:4326');            map.on('click', function(event) {         var coordinate = event.getCoordinate();         if (startPoint.getPosition() == undefined) {           // first click           startPoint.setPosition(coordinate);         } else if (finalPoint.getPosition() == undefined) {           // second click           finalPoint.setPosition(coordinate);                // transform the coordinates from the map projection (EPSG:3857)           // into the server projection (EPSG:4326)           var startCoord = transform(startPoint.getPosition());           var finalCoord = transform(finalPoint.getPosition());           var viewparams = [             'x1:' + startCoord[0], 'y1:' + startCoord[1],             'x2:' + finalCoord[0], 'y2:' + finalCoord[1]           ];           params.viewparams = viewparams.join(';');                // we now have the two points, create the result layer and add it to the map           result = new ol.layer.ImageLayer({             source: new ol.source.SingleImageWMS({               url: 'http://localhost:9090/geoserver/pgrouting/wms',               params: params             })           });           map.addLayer(result);         }       });            document.getElementById('clear').addEventListener('click', function(event) {         // hide the overlays         startPoint.setPosition(undefined);         finalPoint.setPosition(undefined);              // remove the result layer         map.removeLayer(result);       });          </script>      am i missing something? I've been followed the workshop and check it twice. Please help me I've been through so much error whether it using the older tutorial sample or this is the latest tutorial. best regards


title: I need to extract a specific layer from a Geospatial PDF using GDAL/OGR

text: I have GDAL/OGR 1.10 installed on a 64bit Windows 7 machine. I have some Geospatial PDFs that I want to crop to the outline of the map and convert to Rasters. However I cannot clip to the neatline. They have been created using ESRI Production Mapping, and so the neatline is offset from the actual edge of the map, and so clipping by the neatline includes the Axis labels and a lot of whitespace. I can interrogate the Layers of the PDF by using gdalinfo -mdd LAYERS and I know which layer represents the actual neatline/the edge of the map. What I need to know is how to use GDAL or OGR to convert this layer to a shapefile. Thanks!


title: How to pass expression to where clause in ArcPy?

text: I want to pass a user defined variable into a where clause > Syntax > > SelectLayerByAttribute_management (in_layer_or_view, {selection_type}, > {where_clause})               #user input     MTRsel = arcpy.GetParameterAsText(0)          #put user input into where clause     arcpy.SelectLayerByAttribute_management("Townships","NEW_SELECTION", "MTR" = + MTRsel)      I have tried many variations of this but it will not accept it, what am I missing


title: Geographic Transformation Qgis

text: How to Define projection and enter the Geographic Transformation parameters from SGL-GEO to WGS84 in QGIS? Thanks


title: Dynamically set field information in popup template discussion (Automate the process)

text: I want to create a certain function to streamline some of my customization work. I work with a web app that's built using ArcGIS JavaScript API which contains a few layers. Each layer is configured to an identity task that uses a popup template to display attribute information from the layer. Now, there are lots of attributes PER layer that I manually configure up in the popup template declaration. Ex.               var popupContent = {     title: "Business ID#: {UNIQUEID_1}",     fieldInfos: [     {fieldName: "BUSNAME_1", visible: true, label:"Business Name: "},     {fieldName: "Standard_Name", visible: true, label:"Standard Name: "},     {fieldName: "LATITUDE_1", visible: true, label:"Latitude: "},     {fieldName: "LONGITUDE_1", visible: true, label:"Longitude: "},     {fieldName: "INFO_EMP_1", visible: true, label:"Employment Info: "},     {fieldName: "Est_Employment", visible: true, label:"Estimated Employment: "},     {fieldName: "INFO_SAL_1", visible: true, label:"Salary Info: "},     {fieldName: "Est_Revenue", visible: true, label:"Estimated Revenue: "},     {fieldName: "PRMSIC4_1", visible: true, label:"Primary SIC: "},     {fieldName: "NAICS_1", visible: true, label:"NAICS: "},     {fieldName: "NAICS4_1", visible: true, label:"NAICS 4 Code: "},     {fieldName: "NAICS4DESC_1", visible: true, label:"NAICS 4 Desc.: "},     {fieldName: "NAICS2_1", visible: true, label:"NAICS 2 Code: "},     {fieldName: "NAICS2DESC_1", visible: true, label:"NAICS 2 Desc.: "},     {fieldName: "Use_Type", visible: true, label:"Use Type: "},     {fieldName: "UseSubcategory", visible: true, label:"Use Subcategory: "},     {fieldName: "UseCategory", visible: true, label:"Use Category: "}     ]     };      This is one of the shorter templates that i use, but as you can see, it gets pretty tedious to set these manually. Once in a while is okay, but the layers are constantly getting fields added/removed. I was wondering if any of you had any clever ideas as to how I could automate this. One possible way I was thinking is to have a table in my ArcGIS server that acts as a lookup table for the fields. Table structure:   1. The column field names in this table will be the "Layername_fieldname" and "LayerName_label" of the layers in my web app (from the layers' corresponding attribute tables). Ex.               || Layer1_fieldname || Layer1_label || Layer2_fieldname || Layer2_label ||     --------------------------------------------------------------------------     || BusinessLoc      || Biz Location:|| CountyBoundary   || County Name: ||     ||                  ||              ||                  ||              ||        1. Rows can be added or deleted as needed to each column as my web app layers grow or shrink in attribute size. Basically, I want to know if this is a viable option or if there are other methods out there to perform this task. If this is a good solution, is it possible to query the field lookup table and use those values in my popup template constructor. Thanks!


title: How to copy closed loop lines into polygon feature class

text: I'm trying to utilize line contour data to find the volumes of local depressions. To do this, I want to copy the closed loops formed by depression contours into a polygon feature class to find the area of the loops. I've done this several times in ArcGIS 10.0 successfully by selecting the contours I want to copy, then using "Ctrl+C" and "Ctrl+V". My organization just upgraded to 10.1 and the ability to copy from line data to polygon data doesn't seem to work. Is this ability no longer supported, or am I going about this the wrong way?


title: new field names after joining 2 tables in QGIS

text: When joining two files (A.csv and B.shp with a common unique ID) in QGIS the new fields in B.shp (from A.txt) include the name of the field + the name of the file where it comes from (A). For example if we import the field "gender" it will appear in b.shp as A_gender not gender only. Now, the problem comes when the field name es very long because when you save the new joined file (B + A) as a .shp, because it only allows field names with a maximum of 10 characters, the field name is truncated. ¿Any solutions apart for ? Many thanks in advance, Eduardo


title: Printing/exporting WMTS tile layer with desired resolution / zoom level in QGIS

text: I have a project which uses WMTS and some vector layers in QGIS. I can freely zoom in and out, and tiles are refreshed with appropriate zoom level, just as expected. Now I would like to export part of this map (a given bounding box) with the tiles at a fixed zoom level, so that I can better match the printer resolution when map is printed in the desired physical size. How should I tell QGIS to use that given zoom level, and which file format is best suited to do what I want?


title: ArcPy loops/cursors faster without background processing?

text: I am running loops via cursors and raw Python in Arcpy. Within the loop several processes take place, selecting shapefiles and converting to raster, clipping, running euclidean distances, extracting by mask etc. When I would run with background processing enabled things would often crash and background server would throw an error. Now that I have disabled background processing, things seem to be faster and not throwing errors. Is it fundamentally better to run long, complex loops without background processing enabled?


title: Data Feed for Space Weather / Aurora

text: I am building a map application for a northerly area and would love to bring in a feed of northern lights prediction/probability (like below). http://www.gi.alaska.edu/AuroraForecast Does anyone know of a service for space weather that I could digest into leaflet or something?


title: How to register dll extension with RegAsm.exe on Windows 7 64bit ArcGIS 10 SP1?

text: I am having some issues with registering a dll:   1. run DOS as Administrator   2. CD "C:\Program Files (x86)\Common Files\ArcGIS\bin\" (where the ESRIRegAsm.exe is)   3. I attempt to register the dll as follows: Syntax: ESRIRegAsm.exe assembly_name /p:product_name ESRIRegAsm.exe "C:\Program Files (x86)\MyExtension\myextension.dll" /p:desktop I get a warning raised by the registration utility that says: Specified path not found. However, the path is positively correct. What am I doing wrong?


title: Why is my XY Event Source layer only displaying at full extent?

text: My XY Event Source layer is only displaying at full extent and when I zoom in, the features disappear. If I select a feature in the attribute table, the feature is highlighted. When I see this for a normal featureclass, it usually is a problem with spatial index. I first experienced this issue in early Oct. and I did open a ticket with ESRI Tech Support. While on the phone with the analyst, I fired up my MXD and no longer had the problem. I've since created several of these XY layers and all have worked perfectly. Some details: the X,Y theme is created from a query layer (new to ArcGIS 10), which is a SQL Server query to an external SQL Server 2005 database (not SDE and a DB that I'm only a read-only user on). The X,Y values are character values casted to numeric ( _not efficient, but this will only be used occasionally by a few key users_ ). Is anyone else having this issue? I had this problem once and I did unknowingly fix it, but I have no clue as to how. My current workaround is to export to a local file/personal geodatabase--the exported featureclass displays just fine. But I would rather have live access to the data. * * * **Update:** I found a layer file from back in Oct which is working (I must have saved one out). When I add the layer to the map, it works. I'm thinking that there's something in the query that ArcMap isn't liking. However, I'm not sure how to access that query. When I add the layer They query layer, from which the x,y event theme was originally created, isn't in the new MXD. I'll have to dig into this further, but I'm in meetings this afternoon. * * * **Update #2:** On a hunch, I thought the issue might have something to do with the number of rows (71847). So in the query layer, I added `TOP 5` my select statement. I do a `Display X, Y Data...` and I can zoom in or out with no trouble--everything works as expected. So I kept increasing the `TOP N` and each time it worked. Then I put in a `TOP 71847` (the total number of records) into my Select statement. And that works!! That still doesn't answer my question, but it does me get where I want to be (I have to put a `SELECT TOP (SELECT COUNT(*) FROM TABLE) field1,field2,... FROM TABLE` which is messy.


title: Blank extensions box in ArcMap 9.3

text: After a reinstall of ArcInfo 9.3 (which works fine for mapping etc), the box for all of my extensions is blank. I've checked, and they're all licensed in Desktop Administrator, but not showing at all in ArcMap. Does anyone know how I can fix this please?


title: Error in polygonising raster data with QGIS on Mac

text: I'm using QGis Chugiak on a Mac and I'm trying to converting a raster layer in a vector data by: Raster/Conversion/Polygonize. So, I select the layer needed from the drop-down menu, I give a name to the output, layer but every time an error heppens. Probably there's something that doesn't work in GDAL. I leave there the error informations: ![Error in polygonising a raster](http://i.stack.imgur.com/Xk78i.png) Thank you!


title: problem with viewing ASCII data files in ArcGIS

text: I have some problem viewing an ASCII file in a structured way in arcGIS. I have tried running this code but it seems problematic at least with me. I am unable to run this. I think the code needs modification. Can you please tell me the requirements to run this code... and some other help. the link contains the file as well as the IDL source file https://drive.google.com/?authuser=0#folders/0B8iAKVklcaMUVU5qY2NNVFBRbDQ


title: How to share *.sql files between ArcGIS and MapInfo?

text: In Denmark many municipalities have the same provider of public data. These data are pushed to a local database – often a SQL Server. The DBMS layout is the same for all municipalities. When I create a linked table from the SQL Server, I can save (or load) a .sql file in the MapInfo dialog. This .sql-file is portable and I can share it with municipalities with the same setup – but only with MapInfo-municipalities. Recently I had a post at a conference on this item and I urged my colleagues to use my .sql-files and/or my source code. I was however unable to guide my colleagues from the ArcGis-municipalities. Can anyone give me some heads up on this matter? Is it possible to do In ArcGis? The .sql file from MapInfo contains the plain SQL query which is sent to the SQL Server: I.e. :               Select      "JY72000V"."VEJ_NAVN"+ ' '+ "JY72000V"."HUS_NR" As Adresse,      "JY70100V"."KOMMUNENUMMER" AS Kommunenummer,      "JY70100V"."EJENDOMSNR" AS Ejendomsnummer,      Cast("JY72000V"."KOMMUNENUMMER" As varchar(3)) + Cast("JY72000V"."VEJKODE" as varchar(4)) + "JY72000V"."HUS_NR" As ADRESSE_ID     From "ql5020pd"."dbo"."JY70100V", "ql5020pd"."dbo"."JY72000V"     Where "JY70100V"."ADG_ID"="JY72000V"."ADG_ID"      Thanks for any help Søren Breddam, Stevns Kommune


title: It is possible to use different fields for symbology and legend labeling in QGIS?

text: [Edit : working with Quantum GIS version 2.0.1-Dufour] I wonder if it is possible to use some field to set the symbology and some other to label the map. [ **Edit 1 :** I mean labeling the legend in the map designer tool to output as JPG or PDF, not labeling objects in the canvas for working purpose.] [ **Edit 2 :** The purpose of my question is to get an automatical process so I could save time and not manually change all the labels of the legend.] My problem is : I have to make a map of natural habitats, based on european "CORINE biotopes" typology which uses entities with number and names, e.g. : "37 hygrophilous meadows", "38 mesophilous meadows", "41 deciduous forests", "42 coniferous forests", "51 active bogs", etc. So I would like to use shemes of the same color for habitats which are similar : all meadows in different tons of yellow, all forest in different tons of green, mires and bogs in some brown scale etc. To do this I need to base my styling work on the NUMBERS, since following numbers identify similar habitats (meadows are 3X, forests are 4X, etc). But for reading and understanding purpose, the printed map should be labelled with the NAMES. And the labelling should not get "noised" by the numbers, so I separated NAMES and NUMBER of the same entities in 2 different fields. Is there some way to say QGIS to let me do the symbology on the basis of the "NUMBER"-field, but then to label the map on the basis of the "NAME"-field ? I would like to set the "label" (like here from the Layer > Properties > Symbology) ![Layer - Properties - Symbology](http://i.stack.imgur.com/hAKCO.png) To set the label of the legend in the Map Designer Tool : Here I would like to get the names of the habitats, and not the numbers. ![Map designer - Legend](http://i.stack.imgur.com/FzYF2.png) Thanks for all kind of help.


title: QGIS radius buffer question, can it be more circular?

text: Conundrum for the day, When I request a radius buffer in QGIS (Dufor). my buffer comes out as a 20-sided polygon rather than a true "circle". This has lead to some missed blocks when clipping underlying vector layers for population density analyses (e.g., they would be within a circle, but this 20-sided polygon's straight sides preclude fringe blocks from being included in the radius). is there a way I can go into the radius buffer command lines and either increase the number if sides to make the area more circular, or make it truly a "radius" buffer? Clear as mud? I'm trying not to ramble. Thanks. AY


title: How to get Drainage Divide and Streams from DTM

text: Does anyone know how you can get a map like this attached, having the digital terrain model?      I'm sorry, but I do not know what is the exact technical term in English for them.   **I hope you understand from the images.**      (maybe drainage divide, water divide, divide, ridgeline, watershed, water parting) My professor called it _compluvium and impluvium_. But another professor colled it _articulation of the morphological type of the soil_. _(in the Italian language: Carta dei compluvi e displuvi / Articolazione morfologica terreno )_   Thanks in advance for the answers ![enter image description here](http://i.stack.imgur.com/TZ4cR.jpg) ![carta compluvi e displuvi](http://i.stack.imgur.com/hw4tX.jpg) ![carta morfologica del terreno](http://i.stack.imgur.com/65Osh.jpg)


title: MemoryLeak - http://help.arcgis.com/en/webapi/javascript/arcgis/samples/fl_any_projection/index.html

text: Attention ESRI - Loading http://help.arcgis.com/en/webapi/javascript/arcgis/samples/fl_any_projection/index.html inside of sIEve shows a significant amount of memory loss each time the page loads (+10MB). Any workarounds?


title: How do I create a network from connected links

text: I have a bundle of links in 'some' networks. I need to find out which ones are connected to each other i.e. I need to find out which links are not connected to the 'main' network. I have been working on a flood fill recursive method, which effectively goes from link to link, finding out which ones are connected, and recursively doing the same. But the numbers are quite high, and my machine just bailed, so I wondered if there were other, better, or more efficient, ways of doing this. Using arcpy Thanks in advance


title: Wind Farm Viewshed/ZTV

text: I've created an ArcGIS tool using ModelBuilder and **Viewshed** to create a wind farm ZTV and whilst each raster cell tells me whether there's visibility or not it doesn't tell me what the vertical angle of visibiliy is. Does anyone have any ideas how to do this ?


title: Optimize openlayers / Geoserver 2.5 with GeoWebCache - reprojection necessary?

text: I have an openlayers map setup with several base layer options (some served by my own GeoServer) and many optional layers served by my own GeoServer. I've selected Google's EPSG:900913 as the map's projection. And I've specified specific resolutions to optimize use of the GeoWebCache (see below for explanation). Here are the relevant portions of the map's options when created in OpenLayers:                    options = {             projection: new OpenLayers.Projection("EPSG:900913"),             displayProjection: new OpenLayers.Projection("EPSG:4326"),             numZoomLevels: 12,             resolutions: [152.8740566, 76.43702829, 38.21851414, 19.10925707, 9.554628536, 4.777314268, 2.388657134, 1.194328567, 0.597164283, 0.29858214, 0.14929107, 0.07564554],      ...      These resolutions come from looking at the Tile Matrix Set for GridSet EPSG:900913 in GeoServer. These correspond to zoom levels 10 through 22. I did this to optimize the tile requests to match the seeded cached tiles in GeoWebCache. Turns out that many of my layers, served by my GeoServer, have native projection of EPSG:2234 (NAD83 / Connecticut (ftUS)). I see that I can, if I wish, create a new GridSet for caching tiles in this projection. If I just let it create default Tile Matrix resolutions, they will be different from the resolutions for EPSG:900913 GridSet. At least I believe so - these will be in feet rather than in meters, so I'm comparing apples to oranges anyway. But now I'm thinking this whole step is completely uneccessary... Can I still simply seed my native EPSG:2234 layers with the EPSG:900913 Tile GridSet, and then expect tiles requested by OpenLayers to the GWC service using srs:'EPSG:900913' will get created and cached properly? Or should I force reprojection at the layer level, so it is both published and cached at the same projections? I don't want to go through this step if it is completely unnecessary... Thanks for any advice you can provide.


title: Help with ArcGIS10 Python and SearchCursor reading Null values in field

text: I'm using a Phython script to count the number of unique string values in a feature class. My script works if there are values in the field but if there are any NUll values in the field it crashes. Any Ideas?               import arcpy     import operator          inTable = arcpy.GetParameterAsText(0)     stringField = arcpy.GetParameterAsText(1)               stringList = [row.getValue(stringField) for row in arcpy.SearchCursor                                            (inTable, "", "",stringField)]


title: Interaction with two variables as parameters, saving the raster properties output into a table (.dbf)

text: The draft of the model below **TASK** : I need to calculate the percentage of river channel that is covered with shadow in each timestep **RESOURCS:** I have values for azimuth and elevation of the sun (hourly timestep for daylight hours within a year – over 3000 values) **INPUTS:** 1. A table (.dbf) with azimuth and elevation values 2. A river chanel mask (single polygon) 3. A DSM **PROBLEMS** **How to make an interaction with two lists as parameters** The elevation value and azimuth value are assingned to a certain time and date, .I don’t want to calculate hillshade for e.g elevation row 2 and azimuth row 6 –I want only elevation row 1 + azimuth row 1, elevation row 2 + azimuth row 2, elevation row 3 + azimuth row 3… etc… **How to save an output – all the means - back to the input table** – I had an idea to use a calculate field tool, but I don’t know how to “import” the values to the expression field. * * * I've figured out the solution to my first problem. It was easier than I thought. I Do not need to iterate. The simple for of input variable as a list and specifying the table to derive the values. But there is still the other problem - how to write the output back to the input table as a next column![model](http://i.stack.imgur.com/zKRia.jpg)


title: XY data in sketch properties

text: Does anyone know if there is a way to export the xy data from sketch properties in the editor to a .txt file?


title: Extract individual land cover classes from NLCD dataset in ENVI

text: I have a National Land Cover Dataset with 16 different classes. Is there a way to extract the pixels from each individual class and save it as a GeoTiff using ENVI (either Classic or 5.0) so that in the end i'd have 16 GeoTiffs, one for each class?


title: What are the areas, relating to database or server, that should be viewed while making map based webapplication?

text: This question is relating to GIS web application that are based on relational database. I have database that is not indexed with spatial column and I use geoserver to display these data as layers. But for a single layer that is displayed zooming in and out takes time(Intranet or localhost). Are there anything that can be done so that the application can respond faster? I have used:               1.Geoserver     2.SQl-SERVER 2008     3.Openlayers and geoExt     4.No of rows nearly 3000 only for each layers or less (200)


title: Passing parameters to PostGIS SQL queries in Python: weird behaviour

text: The following code always goes to the except section and I do not understand why..Any help would be much appreciated.               l = None     c = conn.cursor()     try:       s = "SELECT ST_NPoints(ST_GeomFromText(%s))"       d = (line,)       c.execute(s,d)       l = c.fetchall()     except:       print "could not get the line vertices"          print line           MULTILINESTRING((-3.16420835153456 55.9269166007097,-3.164222 55.926918),     (-3.1642070167833 55.9269296196706,-3.16421351659546 55.9268662214904),     (-3.16421351659546 55.9268662214904,-3.16421636372824 55.9268384509897),     (-3.16421636372824 55.9268384509897,-3.16422182573761 55.9267851753802),     (-3.16422182573761 55.9267851753802,-3.16422870102352 55.926718114886),     (-3.16422870102352 55.926718114886,-3.16423309121073 55.926675293667),     (-3.16423309121073 55.926675293667,-3.16423565148822 55.9266503211093))      If I execute this in psql, it works fine. I also have the impression that I am following documentation perfectly regardging passign parameters to sql queries etc. When I call the function, I pass the argument as a string like this `'MULTILINESTRING((-3.164....565148822 55.9266503211093))'`


title: How to convert excel csv file to filled contour Google Earth Map?

text: I have a csv file http://www.geocities.ws/radar/data2.csv.   The first column is latitude.   The second column is lonitude.   The third column is precipitation type.   5 represents freezing rain. 15 represents sleet. 25 represents rain. 35 represents snow. 45 represents hail.   How would I translate that into a filled contour plot like the following?   (http://www.mathworks.com/help/examples/graphics/CreateAFilledContourPlotExample_01.png)   I would like this on a Google **Earth** kml file. How would I edit the spacing between the contours?


title: C# library for converting geom field of SQLServer into set of google latlng

text: Is there any C# library which can convert the geometry field into set of google latlngs? I have stored a shapefile into SQLserver table. From Sqlserver i want to read the geometry field data and want to convert this geom data into set of google latlngs in my C# application. Finally i want to save these google latlngs into another table in SQLServer.


title: Python script fails on stand alone but works in ArcMap as a toolbox

text: I have a fresh install of ArcGIS Server and ArcGIS for Desktop both 10.2, latest patches applied. The data is served in this application from a PostgresSQL enterprise geodatabase. The problem is that I can connect to the database using ArcMap, on the ArcGIS Server (the maps that have their data source set in the database) and I can run python scripts that interact with the database in the python console in ArcMap. Yet, the same scripts, or even the smallest interactions (e.g. listing the feature classes with `arcpy.ListFeatureClasses()`) does not work stand alone in python command line. I have tried both the 32 bit and 64 bit Python installations that came with the server and desktop and none run any smallest interactions with the geodatabase. Any thoughts? For reference, a small snippet of code that fails:               import arcpy     arcpy.env.workspace = 'path/connection.sde'     arcpy.ListFeatureClasses()      returns:               []


title: Is the QGIS API suitable for my application

text: I am working on an Linux/QT application that includes a real-time map based display of multiple robotic vehicle positions, where the positions are derived from a variety of sources--some are GPS, others are custom positioning systems that do not produce NMEA type strings. I need to be able to display both custom icons and vehicle track histories. The rest of my application is very domain specific, but it seems that the map display portion ought to be doable using a relatively generic mapping API. I need to be able to interactively zoom and pan, and refresh whenever I get new data. I also need to support geographically based raster underlays. My number of vehicles and history points is relatively small--say several thousand I do this now using a custom C++/Open GL QT application, but will soon be embarking upon a rewrite. Is the QGIS API suitable for this application? Is it fast enough to be interactive?


title: Split lines into non-overlapping subsets based on points

text: Given a table with line geometry, and one or more points that are snapped to this line in a separate table, I would like to split each line with one or more intersecting points at each of the locations where the line intersects a point. For example, there is a line, L, with three intersecting points, A, B, and C in order along the line geometry. I would like to return L as four distinct geometries: from the beginning of L to A, from A to B along L, from B to C along L, and from C to the end of L. In the past I have used shapely for this task, which is a linear referencing problem (http://sgillies.net/blog/1040/shapely-recipes/). However, this would not be practicable in this case, which has many millions of lines and points. Instead, I'm looking for a solution using PostgreSQL/PostGIS. Note that points are constrained to be on a line. Further, a point can validly be on the start or end of a line, in which case the line need not be split (unless there are other points which are not coincident with the same line's start or end points). The subset lines need to retain their direction and their attributes, but the attributes of the point features do not matter.


title: Can a custom add-in execute code upon install of the .esriaddin

text: I have made a custom add-in using C#.NET(a toolbar that contains buttons) and I was wondering if it is possible to have some code execute upon the install of the .esriaddin file. Basically I want to create a directory with some map templates on the users computer. I can just add the               DirectoryInfo di = Directory.CreateDirectory(@"Path"); etc....      to the on click event of the buttons and have it only execute if the folder doesn't exist but I was wondering if there was a way to execute it on install so I don't have to include it within the buttons. Edit: I have looked at this and it is similar though I'm not using a custom installer.


title: directional line symbol ArcGIS vs ArcView 3

text: Is this possible in ArcGIS to have variable number of arrows on the line to show direction? This is an ArcGIS example, where I can have fixed (1, **2,3** ,4,5...) number of arrows and it stays the same doesn't matter what the line length is. First and last always stick to the end, doesn't help with network mapping, where there is usually the node to show. Having just one in the middle is never enough, check line at the right... ![enter image description here](http://i.stack.imgur.com/BzV92.png) No such issues with ArcView 3: ![enter image description here](http://i.stack.imgur.com/g6wu9.png)