title: AGS10 Soap Geometry/Project function vanishes

text: So I am doing a pretty simple process that looks a lot like right out of the SOAP API/SDK for AGS10. My webservices takes a Lat/Lon coord from a mobile app, I reproject it to Web- Mercator to identify against my mapservice data and return a result. Pretty simple, but for some reason that I am not seeing any place, the project method is returning blank, no exception raised, no errors on the server/service side.                           SpatialReference inputSpatialReference = new GeographicCoordinateSystem();                 inputSpatialReference.WKID = 4326; //WGS 84                 inputSpatialReference.WKIDSpecified = true;                      PointN llPnt = new PointN { X = lon, Y = lat, Z = 0 };                 llPnt.SpatialReference = inputSpatialReference;                      SpatialReference outputSpatialReference = new ProjectedCoordinateSystem();                 outputSpatialReference.WKID = 3857; //Web Mercator (Auxiliary Sphere)                 outputSpatialReference.WKIDSpecified = true;                      Geometry[] inputGeometry = new Geometry[] { llPnt };                      bool transformForward = false;                      GeoTransformation transformation = new GeoTransformation();                 // NAD1983_To_WGS1984_1                 transformation.WKID = 108100;                 transformation.WKIDSpecified = true;                      EnvelopeN extent = null;                      Geometry[] outputGeometry = geometryService.Project(inputSpatialReference, outputSpatialReference,                                                                     transformForward, transformation, extent,                                                                     inputGeometry);                      PointN wmPnt = outputGeometry[0] as PointN;            Now; when my defining my points works great; I am able to see the point, my X, Y and Z values as well as my existing coordinate system look great. So I define where my coordinates are going, to, define my transformation since I need to get from Lat/Lon to WebMercator and then after that, fire the GeometryService.Project call.      Now when I look at the output in the debugger of the outputGeometr the point is blank; not even the existing coordinates the new point is null.      So I look to the server, to see what it says; no errors listed, so I am very puzzled. Any thoughts or suggestions?


title: Arcobjects: Zoom to first queried feature

text: I'm attempting to query a feature class for all features that meet a certain criterion (in this case, the `Status` field is `NULL`). However, I'm not getting very far into the attempted method before ArcMap crashes completely. I have no reasonable idea on how to trace errors that cause a program crash, for starters, so I'm more confused than necessary, I expect. I'm building this in VS 2008 Express, for what it's worth.                public void PerformAttributeQuery(ESRI.ArcGIS.Geodatabase.IFeatureClass fc)         {             ESRI.ArcGIS.Geodatabase.IQueryFilter queryFilter = new ESRI.ArcGIS.Geodatabase.QueryFilterClass();             queryFilter.WhereClause = "[Status] IS NULL"; // create the where clause statement       Everything is fine for that first couple of statements, but this query causes a crash:                       // query the table passed into the function and use a cursor to hold the results             ESRI.ArcGIS.Geodatabase.IFeatureCursor featurecursor = fc.Search(queryFilter, false);      It crashes before reaching the next line (I put in debug message boxes for testing), but I'll put it here for propriety:                       ESRI.ArcGIS.Geodatabase.IFeature feature = (ESRI.ArcGIS.Geodatabase.IFeature)featurecursor.NextFeature();             ESRI.ArcGIS.ArcMapUI.IContentsView currentContentsView = ArcMap.Document.CurrentContentsView as IContentsView;             ESRI.ArcGIS.Carto.IFeatureLayer featureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)currentContentsView.SelectedItem;             ESRI.ArcGIS.Carto.IFeatureSelection fSelection = featureLayer as ESRI.ArcGIS.Carto.IFeatureSelection;             fSelection.Add(feature);             FindCommandAndExecute(ArcMap.Application as ESRI.ArcGIS.Framework.IApplication, "{AB073B49-DE5E-11D1-AA80-00C04FA37860}");         }      I cannibalized the `PerformAttributeQuery` snippet to work with FeatureClasses instead of Tables (or so I had hoped), if you'd like a proper reference. My questions are: 1) How can I get better reporting of crash errors? 2) Where did I go wrong? **Answer** : I must shamefacedly admit that I used the wrong characters to denote the location of the feature table I was querying. Changing the `[Status]` in my `queryFilter.WhereClause` to `\"Status\"` quite literally solved all obvious problems. Thank you Kirk!


title: Configuring ArcMap applications to display data from postgis/geoserver

text: I'm afraid im quite new to geoserver and i have been set the task of configuring ArcMap applications to display data from postgis/geoserver. I realize that I will need to set up a WMS service in geoserver but I'm not really sure how one goes about doing this. If any one can give me some pointers or step by step guide it would be much appreciated. I have tried googling it but haven't had much luck. Thanks in advance


title: Style editor in uDig

text: I added layers from Geoserver, I have line objects. But when I open style editor in uDig, I have only WMS Named Style. I thought I could change for example colors by attributes, but I don't how to do it when it offers me only WMS Named Style. I'm working with uDig and Geoserver for the first time so please help me. Thank you.


title: ArcObjects DataStat Error

text: I have some code that adds a bunch of group layers into ArcMap based on a table within ArcMap. The code works fine in VBA, however in VS 2010 I get an error on the **"pEnumVar = pDataStat.UniqueValues"** line (see code below). ERROR: InvalidCastException was unhandled by user code Unable to cast object of type 'System.Runtime.InteropServices.CustomMarshalers.EnumeratorViewOfEnumVariant' to type 'ESRI.ArcGIS.esriSystem.IEnumVariantSimple'. I'm guessing .NET does not like to use the IEmumVariantSimple object? Does anyone have a suggestion of how I can work around this or update this line with a different enumeration object? Thank You!                       Dim pEnumVar As IEnumVariantSimple, value As Object                  Dim pODGSLyr As New ODGSLayer                  'Start             'Find the Layer Files metadata table             pStdTableColl = m_map             For i = 0 To pStdTableColl.StandaloneTableCount - 1                 pStdTable = pStdTableColl.StandaloneTable(i)                 If pStdTable.Name = "ODGSLAYERS" Then                     pTable = pStdTable                     lngFldLayerName = pTable.FindField("LAYERNAME")                     lngFldPath = pTable.FindField("PATH")                     lngFldGroupOrder = pTable.FindField("GROUPORDER")                     lngFldGroupName = pTable.FindField("GROUPNAME")                     lngFldGroupTOCOrder = pTable.FindField("GROUPTOCORDER")                     lngFldGroupVis = pTable.FindField("GROUPVISABLE")                 End If             Next i                  If pStdTable Is Nothing Then                 Exit Sub             End If                  'Sort the Table             pTableSort = New TableSort             With pTableSort                 .Fields = "GROUPTOCORDER, GROUPNAME"                 .Ascending("GROUPTOCORDER") = True                 .Ascending("GROUPNAME") = True                 .QueryFilter = Nothing                 .Table = pTable             End With             pTableSort.Sort(Nothing)                  pCursor = pTableSort.Rows                  'Find Unique Values in the Table             pDataStat = New DataStatistics             pDataStat.Field = "GROUPNAME"             pDataStat.Cursor = pCursor                  pEnumVar = pDataStat.UniqueValues '<<<<<<<<<< ERROR HERE             value = pEnumVar.Next      **Update 1** Update 1 is an update from @Kirk Kuykendall comment. The new error is on the **"lngGroupTOCOrder = pRowLyrs.Value(lngFldGroupTOCOrder)"** line (see code below). I'm not sure what I need to do next. ERROR: NullReferenceException ws unhandled by user code Object reference not set to an instance of an object.                pEnumVar = pDataStat.UniqueValues             value = pEnumVar.MoveNext             Do Until IsNothing(value)                 'Now resort the table based upon the layer order in the group                 pQf = New QueryFilter                 pQf.WhereClause = "[GROUPNAME] = '" & value & "'"                 pLyrCursor = pTable.Search(pQf, False)                 pTblSortLyrs = New TableSort                 With pTblSortLyrs                     .Fields = "GROUPORDER"                     .Ascending("GROUPORDER") = True                     .QueryFilter = pQf                     .Table = pTable                 End With                 pTblSortLyrs.Sort(Nothing)                      'Get the newly sorted rows and create the new Group and Layers inside the Group                 pLyrCursor = pTblSortLyrs.Rows                 pRowLyrs = pLyrCursor.NextRow                      'Create the new Group                 lngGroupTOCOrder = pRowLyrs.Value(lngFldGroupTOCOrder) '<<<<<<< NEW ERROR                 blnGroupVis = pRowLyrs.Value(lngFldGroupVis)


title: Python Efficiency -- Need suggestions about how to use OGR and Shapely in more efficient way

text: I'm looking for some suggestions about how to make my python code more efficient. Normally efficiency doesn't matter for me but I am now working with a text file of US locations with over 1.5 million points. With the given setup it is taking about 5 seconds to run operations on one point; I need to get this figure way down. I'm using three different python GIS packages to do a few different operations on the points and output a new delimited text file.   1. I use OGR to read a county boundary shapefile and get access to the boundary geometry.    2. Shapely checks to see if a point is within any of these counties.    3. If it is within one, I use the Python Shapefile Library to pull attribute information from the boundary .dbf.    4. I then write some information from both sources to a text file. I suspect that the inefficiency lies in having a 2-3 tiered loop... not quite sure what to do about it. I'm particularly looking for help with someone experienced in using any of these 3 packages, as it is my first time to use any of them.               import os, csv     from shapely.geometry import Point     from shapely.geometry import Polygon     from shapely.wkb import loads     from osgeo import ogr     import shapefile          pointFile = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\NationalFile_20110404.txt"     shapeFolder = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New"     #historicBounds = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\US_Counties_1860s_NAD"     historicBounds = "US_Counties_1860s_NAD"     writeFile = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\NewNational_Gazet.txt"          #opens the point file, reads it as a delimited file, skips the first line     openPoints = open(pointFile, "r")     reader = csv.reader(openPoints, delimiter="|")     reader.next()          #opens the write file     openWriteFile = open(writeFile, "w")          #uses Python Shapefile Library to read attributes from .dbf     sf = shapefile.Reader("C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\US_Counties_1860s_NAD.dbf")     records = sf.records()     print "Starting loop..."          #This will loop through the points in pointFile         for row in reader:         print row         shpIndex = 0         pointX = row[10]         pointY = row[9]         thePoint = Point(float(pointX), float(pointY))         #This section uses OGR to read the geometry of the shapefile         openShape = ogr.Open((str(historicBounds) + ".shp"))         layers = openShape.GetLayerByName(historicBounds)         #This section loops through the geometries, determines if the point is in a polygon         for element in layers:             geom = loads(element.GetGeometryRef().ExportToWkb())             if geom.geom_type == "Polygon":                 if thePoint.within(geom) == True:                     print "!!!!!!!!!!!!! Found a Point Within Historic !!!!!!!!!!!!"                     print str(row[1]) + ", " + str(row[2]) + ", " + str(row[5]) + " County, " + str(row[3])                     print records[shpIndex]                     openWriteFile.write((str(row[0]) + "|" + str(row[1]) + "|" + str(row[2]) + "|" + str(row[5]) + "|" + str(row[3]) + "|" + str(row[9]) + "|" + str(row[10]) + "|" + str(records[shpIndex][3]) + "|" + str(records[shpIndex][9]) + "|\n"))             if geom.geom_type == "MultiPolygon":                 for pol in geom:                     if thePoint.within(pol) == True:                         print "!!!!!!!!!!!!!!!!! Found a Point Within MultiPolygon !!!!!!!!!!!!!!"                         print str(row[1]) + ", " + str(row[2]) + ", " + str(row[5]) + " County, " + str(row[3])                         print records[shpIndex]                         openWriteFile.write((str(row[0]) + "|" + str(row[1]) + "|" + str(row[2]) + "|" + str(row[5]) + "|" + str(row[3]) + "|" + str(row[9]) + "|" + str(row[10]) + "|" + str(records[shpIndex][3]) + "|" + str(records[shpIndex][9]) + "|\n"))             shpIndex = shpIndex + 1         print "finished checking point"         openShape = None         layers = None               pointFile.close()     writeFile.close()     print "Done"


title: Create a Python script to merge broken polylines from a microstation CAD file in ArcGis

text: I have tried using the dissolve tool with unsplit_lines option but it does not merge the lines... maybe there is a step I am missing?


title: Why latitude/longitude and not the other way round?

text: When you google "latitude/longitude" you get 10 times more results than "longitude/latitude". This I find very confusing as "x/y" seems 20 times more common than "y/x". But on a map Latitude is on the Y axis and longitude on the X. Maybe I'm just venting, but maybe I'm missing a clue that makes all this madness logical. Any explanations?


title: IFeatureClassWrite and ITableWrite doubt

text: I'm looking at these specific interfaces ITableWrite and IFeatureClassWrite. I have a need to store some features (mostly a custom log) and I can't having them triggering my editor events (OnCreate, OnChange and OnDelete - but more specifically OnCreate). The document is clear that IFeatureClassWrite will NOT trigger these events (IRow.Store). ITableWrite just says: > " .WriteRow | Write a row to the database without polymorphic object update > behavior. " Do they mean the same thing? e.g. using WriteRow or WriteRows will not trigger the OnStore event (which BTW is not listed on IRowEvents) and using IFeatureClassWrite will not trigger OnStore events? Another question: do I need to bypass these events? I'm only listening to OnCreate, OnChange and OnDelete events of the editor. Also, is that the best tatic to capture these specific events? I'm afraid that programmatically, someone else might use IWorkspaceEdit, and my editor event handler won't listen to those (events). Sorry for the long question, but the setup to test this is kinda of big, and I'm open to advice. Thanks!


title: Rescale Raster Band With GDAL Python Bindings

text: I am trying to figure out how to rescale a raster band's pixel value using the GDAL python bindings. I am looking to emulate the -scale switch in gdal_translate.exe; rescaling the pixel values in my raster to 0-255. EDIT: Some goofing around today has produced this, seems to work thus far.               from osgeo import gdal     import numpy     import math          def singleBandToRGB(source,output):         src_ds = gdal.Open(source)         band = src_ds.GetRasterBand(1)         min = band.GetMinimum()         max = band.GetMaximum()         data = src_ds.ReadAsArray(0,0,src_ds.RasterXSize,src_ds.RasterYSize)         data = ((255-0)*((data-min)/(max-min)))+0         data = data.astype(numpy.uint8)         driver = gdal.GetDriverByName("GTiff")         dst_ds = driver.Create(output,src_ds.RasterXSize,src_ds.RasterYSize,3,gdal.GDT_Byte)         for i in range(1,4):             dst_ds.GetRasterBand(i).WriteArray(data)             dst_ds.GetRasterBand(i).ComputeStatistics(False)         dst_ds.SetProjection(src_ds.GetProjection())         dst_ds.SetGeoTransform(src_ds.GetGeoTransform())              # housekeeping         driver = None         dst_ds = None         src_ds = None


title: arcpy taking longer and longer in a loop

text: I am looping through an amount of point arrays, creating data on the fly. The features/arrays, are largely similar in size. One thing I have noticed, is that when I start the process, the time it takes to perform, rises from about 8 seconds, to around 1 minute. I can't see why this is happening, but you cna see the times slowly increase if you put a timer around the loop. Is there a way to increase this performance? Is there a flush method? I believe I am properly handling my objects and deleting them after use, but I cannot for the life of me see why it is happening? Any ideas? The code is refactored, so there may be some typos, etc.               def createGeom(geom,header,band,scratchDB,buffer, dist):          filetime = (str(time.time())).split(".")     outfile = "fc" + filetime[0]+filetime[1]     outpath =  scratchDB + "tmp.gdb/Polygon/"     print outpath     sCon = ["1", "16", "2.0"]      ts = header.split(" ")     ti = ts[6]     Current = ts[9]     LL = ts[2][2:5]     LU = ts[2][8:11]          outFeatureAggClass = outpath+outfile+"_agg"     outFeatureClass = outpath+outfile          arcpy.AggregatePoints_cartography(geom, outFeatureAggClass, dist)     arcpy.Buffer_analysis(outFeatureAggClass, outFeatureClass, buffer)                    arcpy.AddField_management(outFeatureClass, "Name", "Text")     arcpy.AddField_management(outFeatureClass, "ID", "SHORT")     arcpy.AddField_management(outFeatureClass, "LL", "SHORT")     arcpy.AddField_management(outFeatureClass, "LU", "SHORT")             arcpy.AddField_management(outFeatureClass, "Tim", "SHORT")      arcpy.AddField_management(outFeatureClass, "Band", "SHORT")     arcpy.AddField_management(outFeatureClass, "Lvl", "TEXT")      arcpy.AddField_management(outFeatureClass, "Cur", "SHORT")               rows = arcpy.UpdateCursor(outFeatureClass)     for row in rows:         row.name = outfile         row.ID = 1000         row.LL = int(LL)         row.LU = int(LU)         row.Tim = int(ti)            row.Band = int(band)         row.Lvl = sCon[int(band)]         row.Cur = Current         rows.updateRow(row)          del row     del rows          **UPDATE** I am now running this in ArcMap, and it's taking significantly less time, to a factor or about 3 or 4; what took 30 mins, now takes 7. What's that all about? **EDIT** ESRI have agreed to look at the issue. I'll park this for now, ready to update it with any progress.


title: Using ESRI.ArcGIS.Display.ITemplate to draw cutom lines?

text: I am using ESRI.ArcGIS.Display.ITemplate to draw custom lines such as dot/dash lines on a map control. Then, I used below codes               ESRI.ArcGIS.Display.ITemplate templateCls = new ESRI.ArcGIS.Display.TemplateClass();     templateCls.Interval = 1;     templateCls.AddPatternElement(1, 3);      The below image 1 will show the result for this. ![enter image description here](http://i.stack.imgur.com/4fxfV.png) But, I am want to the line include round elements (not rectangle) as the below image 2. ![enter image description here](http://i.stack.imgur.com/dZdy6.png) Are there anyone know how to do that ? Thank you very much !


title: How do I create a custom geostamp in the TerraGo toolbar

text: The blurb for the new 6.0 tool bar says it supports custom geostamps, but I can't find any info on how to create one. Can anyone point me to the method for creating a custom geostamp that the toolbar will recognize? Edit: So far, this is what I've figured out: The geostamp files on Windows are in C:\Program Files\Common Files\TerraGo\GeoStamps   1. Each collection of stamps in the Geostamps menu corresponds to a pdf file in that directory.   2. Each file begins with a sequential letter of the alphabet. The installation has A thru E.   3. The name of the collection corresponds to the document name (which in the default files is the file name minus the first letter).   4. If you copy one of the existing collection files, rename the document (in Acrobat), rename the file to start with the next letter of the alphabet (F), you will see a new collection in the GeoStamps menu in Acrobat Reader.   5. The pdf file consists of one 8.5x11 page (apparently blank), plus a series of pages each .65" square containing a single symbol, which is a glyph from a font whose name essentially matches the collection name.   6. If you delete one (or more) of these pages, the file still works, but the corresponding stamp is missing. And here's where it breaks down - 8. If you add a page to this file, it does NOT appear as a symbol available to the Geostamp function in the toolbar. Interestingly, the added page does not break the collection, even if placed before the original pages. If you add markup to a symbol from the original document (like an arrow, for example), the markup appears in the symbol preview, but no in the symbol "stamped" onto the document. I've checked all the document properties and can't find anything that seems out of the ordinary. But there must be something magic on the page to make it work.


title: Why is polygon used as range in GIS and maps?

text: I am enroute to create an app that finds nearest venues. I already have another SO question on this regard. I don't have any experience with such features, and some of the concepts are not being very obvious to me. So, I have been thinking that I am not taking things on right direction. I have seen in number of places with both GIS and in case of mobile towers, the range is calculated in polygons. I have also been reading a bit about Haversine Formula, and it appears it does not use Polygon, it is basically two-point distance (from my understanding). Why use Polygon for reference ? Why is circle not used for this ? How does Haversine formula integrate with this polygons concept ? #Copied : Original question in SO as refered by Nate


title: Retrieve road segment from Google or Bing Maps

text: Let me see if I would be able to explain this correctly... Both Bing and Google Map geocoders can return you and a point given an address. Is there any way to retrieve an 'address range' instead i.e. to get returned a LINE rather than a point? Its not necessarily address ranges I am after- the key word here is 'line' (along a road). Another example would be to get a road line given 2 points so I would want to get returned the road segment between these 2 points. I assume for this I would need to use the routing functions only? Another scenario (i guess even more difficult) is given 1 point (along a road) to ask for the road segments for x Km 'up or down' this same road Hope I am making sense...


title: How can I convert esriMeters to Lat/Lng?

text: I am pulling a JSON list of features into a web page. The points come as esriMeters and the rest of my app uses Lat/Lng. How can I do this conversion in the browser (javascript)?


title: Geoserver Problem

text: The following error message shows during an attempt to open geoserver url... "HTTP ERROR: 404 NOT_FOUND RequestURI=/geoserver" Note: I have restared the Start/Stop .bat files several times & also restarted the PC.


title: Arcpy: Troubleshooting parameter error for Export Report

text: I'm running into what I believe is a parameter error in writing an arcpy exporttoreport script. I've written several successful scripts of this same type, but the difference is that this one has a "DEFINITION_QUERY", which requires a query. Because of the order of parameters, I need to have something between the where_clause and field map, but using "" doesn't seem to be working. Below are the script and error messages. For the second error message I tried something I saw on the help page of including a double bracket and then ...",True) at the very end, that prompted an error about 3 versus 2 arguments. Any suggestions would be much appreciated! **Current and clean version of code**               import arcpy          arcpy.env.overwriteOutput =True          mxd = arcpy.mapping.MapDocument("C:\\DM_Final_submission\\Graffiti_DM3_5.mxd")     df = arcpy.mapping.ListDataFrames(mxd, "*")[0]          New_layer= "C:\DM_Final_submission\\New_Graffiti occurrence.lyr"     zip_code='"Incident_Z" = ' + "'10002'"       where_clause='"Incident_Z" = ' + "'10002'"       field_map ={"Incident_Z":"zipcode", "Incident_A":"Address", "Precinct": "pp"}     finalpdf="final.pdf"          graffitilyr=arcpy.mapping.ListLayers(mxd, "*", df)[0]     expdir="Z:\\DM_Final_submission"     rlf= "C:\\DM_Final_submission\\2013_Graffiti.rlf"     export= expdir + "\\" + "repexpnamepdf.pdf"     newpointout="Z:\\DM_Final_submission\\newpointout1.shp"      finalpdf="Z:\\DM_Final_submission\\" + finalpdf          newpointlayer="Z:\\DM_Final_submission\\newpointlayer.lyr"     pdfddp="Z:\\DM_Final_submission\\pdfddp2.pdf"     ddp = mxd.dataDrivenPages     nycstreets="C:\\DM_Final_submission\\NYCStreets.lyr"     newnycstreets="Z:\\DM_Final_submission\\newstreetsout.shp"       # this code block selects a series of points and outputs them as a new shapefile based upon a query (zip_code)               arcpy.Select_analysis(New_layer, newpointout, zip_code )      # this code block selects a series of streets and outputs it as a new shapefile based upon a query (zip_code)               arcpy.Select_analysis(nycstreets, newnycstreets, zip_code)       # this code block takes the previously created new point .shp file and and outputs it as a new layer               arcpy.MakeFeatureLayer_management(newpointout, newpointlayer)      # this code block creates the report (graffitilyr is the input, rlf is the .rlf, export is the export path, "DEFINITION_QUERY" is the parameter for data_set_option, New Graffiti Report" is the new title, 1 is the starting page, 1-5 is the range, where_clause is a query for definition query, field map is field map)               arcpy.mapping.ExportReport(graffitilyr, rlf, export,"DEFINITION_QUERY","New Graffiti Report", 1, "1-5",where_clause,field_map=field_map)      # the following code does DDP, creates a blank pdf, and then appends a variety of pdfs, including the report to it.               ddp.exportToPDF(pdfddp, "ALL")     mapbookPDF= arcpy.mapping.PDFDocumentCreate(finalpdf)     mapbookPDF.appendPages (expdir + "\\overviewmap.pdf")     mapbookPDF.appendPages(expdir + "\\repexpnamepdf.pdf")     mapbookPDF.appendPages(expdir + "\\pdfddp2.pdf")     mapbookPDF.saveAndClose()      **Edited, and edits at bottom**                arcpy          arcpy.env.overwriteOutput =True          mxd = arcpy.mapping.MapDocument("C:\\DM_Final_submission\\Graffiti_DM3_5.mxd")     df = arcpy.mapping.ListDataFrames(mxd, "*")[0]          New_layer= "C:\DM_Final_submission\\New_Graffiti occurrence.lyr"     zip_code='"Incident_Z" = ' + "'10002'"       where_clause='"Incident_Z" = ' + "'10002'"       field_map ={"Incident_Z":"zipcode", "Incident_A":"Address", "Precinct": "pp"}     finalpdf="final.pdf"          graffitilyr=arcpy.mapping.ListLayers(mxd, "*", df)[0]     expdir="Z:\\DM_Final_submission"     rlf= "C:\\DM_Final_submission\\2013_Graffiti.rlf"     export= expdir + "\\" + "repexpnamepdf.pdf"     newpointout="Z:\\DM_Final_submission\\newpointout1.shp"      finalpdf="Z:\\DM_Final_submission\\" + finalpdf          newpointlayer="Z:\\DM_Final_submission\\newpointlayer.lyr"     pdfddp="Z:\\DM_Final_submission\\pdfddp2.pdf"     ddp = mxd.dataDrivenPages     nycstreets="C:\\DM_Final_submission\\NYCStreets.lyr"     newnycstreets="Z:\\DM_Final_submission\\newstreetsout.shp"           arcpy.Select_analysis(New_layer, newpointout, zip_code )     #arcpy.Select_analysis(nyczip, newzipcode, specifyzipzip)     arcpy.Select_analysis(nycstreets, newnycstreets, zip_code)      arcpy.MakeFeatureLayer_management(newpointout, newpointlayer)     arcpy.mapping.ExportReport (graffitilyr, rlf, export,"DEFINITION_QUERY","New Graffiti Report", 1, 1-5,where_clause, " ", field_map)     ddp.exportToPDF(pdfddp, "ALL")     mapbookPDF= arcpy.mapping.PDFDocumentCreate(finalpdf)     mapbookPDF.appendPages (expdir + "\\overviewmap.pdf")     mapbookPDF.appendPages(expdir + "\\repexpnamepdf.pdf")     mapbookPDF.appendPages(expdir + "\\pdfddp2.pdf")     mapbookPDF.saveAndClose()               Traceback (most recent call last):       File "C:\\New_file_New_query1.py", line 30, in <module>         arcpy.mapping.ExportReport (graffitilyr, rlf, export,"DEFINITION_QUERY","New Graffiti Report", 1, 1-5,where_clause, " ", field_map)       File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_         return fn(*args, **kw)       File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\mapping.py", line 514, in ExportReport         return report_source._arc_object.ExportReport(*gp_fixargs((report_layout_file, output_file, dataset_option, report_title, starting_page_number, page_range, report_definition_query, extent, field_map), True))     ValueError: Can't set these parameters with selected dataset_option     >>> ================================ RESTART ================================     >>>           Traceback (most recent call last):       File "C:\\New_file_New_query1.py", line 30, in <module>         arcpy.mapping.ExportReport ((graffitilyr, rlf, export,"DEFINITION_QUERY","New Graffiti Report", 1, 1-5,where_clause, " ", field_map), True)       File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_         return fn(*args, **kw)     TypeError: ExportReport() takes at least 3 arguments (2 given)     >>> ================================ RESTART ================================      Based upon the suggestion below of @GetSpatial, I edited the code to read               arcpy.mapping.ExportReport(graffitilyr, rlf, export,"DEFINITION_QUERY","New Graffiti Report", 1, "1-5",where_clause,field_map=field_map)      However that produced the following error. In reviewing the traceback, what is interesting is A) that it says the error is that it can't open the report layout template (note I checked in ArcMap and the .rlf opens fine, I also reran another version of the script that uses the .rlf and it ran correctly), and particularly interesting is B) something I hadn't previously seen in the traceback above is that it has a different order for the parameters than that of the help page:               (*gp_fixargs((**report_layout_file**, output_file, dataset_option, report_title, starting_page_number       Specifically, the report source file, which should be first isn't there, but it starts with the .rlf. Could this be the cause (note have also tried taking out the input and replicating that order above from the error, but that doesn't work either) Traceback (most recent call last):               File "C:\DM_Final_submission\New_file_New_query1.py", line 30, in <module>         arcpy.mapping.ExportReport(graffitilyr, rlf, export,"DEFINITION_QUERY","New Graffiti Report", 1, "1-5",where_clause,field_map=field_map)       File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_         return fn(*args, **kw)       File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\mapping.py", line 514, in ExportReport         return report_source._arc_object.ExportReport(*gp_fixargs((report_layout_file, output_file, dataset_option, report_title, starting_page_number, page_range, report_definition_query, extent, field_map), True))     IOError: Could not open report template


title: How to draw perpendicular lines in QGIS?

text: I have a situation like this: ![enter image description here](http://i.stack.imgur.com/eCCHj.png) What I need to do is to connect each point to every line that is at most, let's say 200 m, away from the point. In other words, I need to draw a perpendicular line from each point to every line that is in the buffer. Is there a way to do this in QGIS and PyQGIS? Thanks


title: Select by location of doughnut hole polygon

text: I have vectorized the Corine Land class raster and want to select certain classes that are "inside" other classes. ![Example: selecting polygons circled in blue](http://i.stack.imgur.com/TgvYH.png) How can I select the polygons circled in blue? All the select by location options do no work since there are no overlaps. Explanations per the comments: What you see is the area around Barcelona, Spain. the red outlines urban polygons, consisting of several types of CLC classes that have been polygonized from the original Corine raster. The green polygons are "green" areas in and around the city. CLC, as a raster dataset, has only one class for every 100*100 meter pixel through the whole of Europe. What I'm trying to do is select "green" polygons that are completely surrounded by "urban" polygons, and then export them as a separate layer.


title: Finding lines in point patterns?

text: I have a point pattern (see below) in which I assume to find "significant" lines of points. However, I did not find any method that is integrated in a GIS or in R that searches for such line pattern in a spatial point pattern. There is one method that is implemented in the free software called PAST (http://folk.uio.no/ohammer/past) called point alignment. the method is described in : > Hammer, Ø. (2009). New statistical methods for detecting point alignments. > Computers & Geosciences, 35(3), 659-666. beyond this method I am not aware of any tool that finds lines in point patterns. Any hints on possible tools, r-packages, etc. are most welcome. ![Spatial pattern of individual plants](http://i.stack.imgur.com/2slf2.jpg)


title: Can I use leaflet to access a WMS service published as an OpenLayers javascript API?

text: I would like to add this layer to a site I'm working on that uses leaflet:               http://116.240.195.134/vicmapapi/UserGuide.jsp      The API itself is a javascript file based on OpenLayers: http://api.maps.vic.gov.au/vicmapapi/js/vicmapapi.js (I don't know OpenLayers at all, so I'm not really sure what I'm looking at). Am I able to access the tiles using leaflet? This didn't work: var map = L.map('map').setView([-37.3, 145.3], 5);                       var wms = L.tileLayer.wms("http://api.maps.vic.gov.au/vicmapapi/map/wms", {             LAYERS: 'CARTOGRAPHICAL',             FORMAT: 'image/png',             FIRSTTITLE: false,             transparent: true,             WIDTH:512,             HEIGHT:512,             EXCEPTIONS:"application/vnd.ogc.se_inimage",             SRS:"EPSG:3111"             });             wms.addTo(map);      In particular, it didn't work because the tiles it tried to fetch had URLs like this:               http://api.maps.vic.gov.au/vicmapapi/map/wms?service=WMS&request=GetMap&version=1.1.1&layers=&styles=&format=image/jpeg&transparent=true&height=256&width=256&LAYERS=CARTOGRAPHICAL&FORMAT=image/png&FIRSTTITLE=false&WIDTH=512&HEIGHT=512&EXCEPTIONS=application/vnd.ogc.se_inimage&SRS=EPSG:3111&srs=EPSG:3857&bbox=18785164.071364917,-3757032.814272983,20037508.342789244,-2504688.542848655      instead of:               http://api.maps.vic.gov.au/vicmapapi/map/wms?LAYERS=CARTOGRAPHICAL&FORMAT=image%2Fpng&FIRSTTILE=false&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&SRS=EPSG%3A3111&BBOX=2327867.7504022,2539132.2495979,2869735.5008043,3081000&WIDTH=512&HEIGHT=512&cache_bbox=2296369.732799%2C2279293.001349%2C2931371.002802%2C2702627.181351      The server then returns error messages like "400: (x + width) is outside of Raster". If I increase WIDTH to 10000, then Apache returns 500 Internal Server error with a complaint about Java heap space. Oops. So - am I barking up the wrong tree here, or is there a way to do this?


title: Need to make a buffer around a part of a feature, keeping other features attached

text: Basic problem: I have two layers, a layer containing river basin outlines and a layer of wetland outlines and types. I need to make a buffer around one river basin, and then create 100 random points within a single type of wetlnad within the river basin + buffer. I guess the first problem is how to create a buffer around just one river basin when the layer contains many (they do not overlap). The second problem is that I need to isolate just one wetland type within the total area contained by the buffer (so that river basin + the buffer around it), so that I can make random points inside that wetland type. I've been fighting with trying to make a buffer around just that one river basin for days, I'm sure it's possible but I've never used buffers before. How do I do this? Thanks in advance for any help!   * Blue


title: How to have QGIS show the number of points that overlay the raster polygon?

text: I have a number of single-value raster layers defining a area mask and a point layer pulled from PostGIS. Is there a way to have Quantum GIS show me the number of points that overlay the raster polygon? Or is my only option to convert the raster to a vector, pull into the database, and run a query against it?


title: Map formats supported by QGIS

text: May i know whether QGIS supports the following map formats? (1) S-57 Electronic Navigation Charts (Naval Maps) (2) Army/IAF Maps in DGN (3) DTEDO0, DTEDO1 and DTEDO2 (4) GeoTIFF (5) DTED


title: Trouble storing an IMap with ArcObjects 10 and IPersistStream

text: I'm storing a reference to an IMap in a custom object (which is stored by a custom extension) so that it can listen to map events. When I save the custom object using IPersistStream and then reload the MXD, it appears that I now have a reference to a COPY of the original IMap. I tested this by making a simple custom extension that saves the first IMap in the document... when reopening the saved MXD, the saved IMap != the first IMap in the document. Is this how it should be working? Should I just be saving the location of the map in the document (rather than the actual IMap) and then resetting the reference during IPeristStream.Load? I'm using C#, ArcGIS 10.0, and VS 2008. Thanks for any help.


title: QGIS processing does not support jpg

text: I would like to work with the gdal functions provided in the processing toolbox. It is easy to call:               import processing     processing.alglist()     processing.runalg("gdalogr:warpreproject","/home/USER/.qgis2/python/plugins/qgis2leaf/test_data/src_img.tif","EPSG:4284","EPSG:4326",0,1,"","/home/USER/Desktop/export_image.img")      as an example. unfortunately I need jpg in the ouput or png. Is there any reason why I cannot save as jpg. I can use the warp or translate function in the GDAL toolbox with a jpg but it is hard for me to call it from within python. Back to the question: Is there a support for jpg in the GDAL processing functions in QGIS?


title: Python ogr2ogr language character encoding

text: In python script i'm exporting .shp files to .csv/.gpx/.kml               import ogr, ogr2ogr     ...     ogr2ogr.main(["","-f", "CSV", csvfile, shpfile, "-lco", "GEOMETRY=AS_WKT", "-lco", "SEPARATOR=SEMICOLON"])     ogr2ogr.main(["","-f", "GPX", "-dsco", "GPX_USE_EXTENSIONS=YES", gpxfile, shpfile])     ogr2ogr.main(["","-f", "KML", kmlfile, shpfile])      The files are exported, but the language(Latvian) specific string field characters are incorrect. Shape files have string field "LABEL" with a value i.e. **"marķieris"**. In case of .kml file it is turned into **"markieris"** , and in case of .csv and .gpx it is turned into **"maríieris"**. All exported files have utf-8 (without BOM) encoding. So far i had no luck finding if i have to pass some kind of additional parameter, or do something else to fix this language issue. Scrip is executed in Windows.


title: "Shapefile to WMS" with FME

text: I'm using FME Desktop and Server 2014.   At the moment my workbench looks like the following. It reads a .csv file and transforms this to a ESRi Shapefile. This is working. I can use the shapefile with QGIS or ArcGIS Desktop. ![enter image description here](http://i.stack.imgur.com/U29ij.png) The next step, which I want to perform: Create a WMS Service Link, which uses my FME or ArcGIS Server. Because FME provides the functionality to publish the workbench as an OGC Web map Service, I would assume that this works. ![enter image description here](http://i.stack.imgur.com/jhjQf.png) I've already tried to publish the workbench as shown in the first pic. The WMS Service is up and running and I can load the WMS Service to ArcMap, I can query the WMS link and the layer is shown. Unfortunately when I want to add the layer to the map, I get the error **"A selected item could not be added to the map. Could not connect to server".** So my first question: Is there any error with my workbench? Second "question": If searched for some solutions and I found that I need to have a result as PNG, JPEG, TIF.   So my question would be: What do I need to add/edit in my workbench to get my idea to work? Any help is much appreciated!


title: QGIS V. parallel command using GRASS

text: I cannot offset a line to a specific decimal place, (ie only whole numbers can be inputted) Am I missing something here??


title: How to run a Spatial Query

text: I have a shapefile of a forest inventory of the uk. I want to find out the total area of woodland within a 60km radius of a specific point. I have added a layer for the point and added a 60 km buffer. so i now want to know the total area of woodland within this areas


title: Why would an ArcGIS 10 Python custom tool fail without throwing a "failed to execute" or "completed script" message at the end?

text: I have written a custom tool in Python for ArcGIS that converts in batch SHP files to GML files. It processes one file at a time in a folder, file, or ZIP file and then moves on to the next. If the tool experiences an exception it states that an error happened and just continues to the next file. What one of the users of the tool experienced after kicking it off overnight is that it ended in a strange way. Usually when the tool ends (prematurely or not) it has a message at the end saying "failure to execute " or "completed script ". This time, however, it did not (as the image below shows). ![Image showing lack of confirmation of completion or error.](http://i.stack.imgur.com/X8LLw.jpg) Is there a specific reason why the script would just end abruptly like this, like maybe the process just quit for some reason? I asked the user if the computer restarted or if he noticed if ArcCatalog closed and he said no.


title: Z Values - KML to Layer

text: I have short question concerning the z - values of KML - Files. I am importing a KML - File of linestring s with the KML to Layer Function (ArcMap 10), but it seems that the z values of the KML (10.3035,47.0045, **1656** ) is lost during the import? I tried different approaches for environment settings of the z-values, but nothing seems to help?


title: Delete layer from ArcMap using Python

text: **Problem:**   * I'm trying to loop through all my layers looking for layer named `"CADAnnotation"`.   * If the layer exists then remove the layer from the mxd **Notes:**   * Running from Stand-Alone script (ie. NOT within Arcmap)   * `"CADAnnotation"` Data Type is a CAD Annotation Feature Class   * `"CADAnnotation"` is NOT in a geodatabase, it's created from a AutoCAD .dwg   * `"CADAnnotation"` is within a Group Layer named `"ACAD"`   * If the group layer `"ACAD"` can be deleted that also removes `"CADAnnotation"` that would be great. **Code thus far:**               for item in mxds:         print (item)         mxd = arcpy.mapping.MapDocument(item)         df=arcpy.mapping.ListDataFrames(mxd,"Project Area")[0]         for lyr in arcpy.mapping.ListLayers(mxd, "*",df):             if lyr.name == "CADAnnotation":                 print(lyr.dataSource)                 arcpy.Delete_management("CADAnnotation")                 print("Layer Deleted")             else:                 pass      **Notes on Code:**   * I can find the layer no problem   * the line `arcpy.Delete_management("CADAnnotation")` does not work throws an error. **Question:**   * How to I remove `"CADAnnotation"` and/or `"ACAD"` group layer?


title: How to overlay points from MySQL served by Geoserver on an Openlayers map?

text: This might be a very basic question as I am new to GeoServer. **What I am trying to do:**   1. I have a MySQL table that stores data attached to a spatial point and have latitude and longitude columns to store the location attribute.   2. Use GeoServer to display these data points as layer. **What has been done:**   1. Create a new table which is duplicate of the original table but also have POINT field which is generated using POINT(orig_lat, orig_lon).   2. Use the MySQL extension to create a new data store linked with this spatial_table.   3. Create a new layer here we need two fields called as native and declared SRS, in the demo I had seen EPSG:4326 as the declared SRS and hence I have added that SRS. **Problem:**   1. When I view the layer preview in openlayer I do not see any map? Not sure if this is because of SRS or do I need to specify a map layer.   2. Most of the points seem like they are overlayed on top of each other and zoomin and panning does not make much difference. This seems like SRS issue.    3. Is it possible to overlay the points on top of google maps using geoserver?


title: Low resolution raster conversions to vector yielding millions of single pixel polygons. Trying to merge all directly adjacent polygons

text: I'm trying to scan plain pdf maps containing color-coded geographic regions and convert them into undifferentiated (I don't actually need them to contain any attribute data) overlays. My workflow is to take one of my color-coded maps, select a few of the color-coded regions via Gimp's "Select by Color" tool, cut those out into a separate layer, then export just that layer with a transparent background as a tiff. I then georeference that tiff, and proceed to polygonize it. Polygonization of poor quality images is causing what should be 100-1000 polygons (map regions) into 1 million polygons, most of which are 1-4 pixels in size (area). So far a solution hasn't been forthcoming. So on the GIS front, I'm hoping I can at least merge directly adjacent polygons. Here's what they look like: https://drive.google.com/file/d/0B5inC0VAuhH1bHhuRXVTMHZkOEk/edit?usp=sharing In QGis, I've tried the "dissolve" tool, the "buffer tool" & the "Multipart(s) to singlepart(s)" tools, as well as the "Densify" tool & and "Simplify" tool. None have really worked so far. The buffer tool, which has worked in the past, turns the entire area into a single polygon, even when I try to get it to approximate the maximum (99) number of segments. Other methods either do nothing, fill the entire screen with a single polygon, or take literally days to process. So maybe there is a "Polygonize" method (from raster) that will best do the job here, or maybe its a post-polygonize job. Thanks for taking a look


title: finding administrative boundaries in Guinea using OSM

text: I am trying to map the quartiers (neighborhoods) in the city of Conakry, Guinea. As far as I can find, OMS contains information on the location of each quartier, and displays this with a point in the center of each quartier. I need the boundaries of each quartier. I have checked www.gadm.org, they have boundaries up to the point of commune (sub-prefecture), which is one level above quartier. My question is, is this information in OSM, and did I miss it (no purple lines), if not, is there an other way to find it?


title: How to Create Custom ZoomIn Button that is Windows Tools Button with Custom Cursor?

text: I want to create large zoomin Button without toolbarbutton. I created a 'Zoomin Script' on axMapContorol. It worked is fine but I can't set cursor. My code is as follows:               IMapControl3 mapControl = (IMapControl3)FMap.Instance.axMapControl1.Object;     MemoryStream ms = new MemoryStream(LargeZoominCursor);     Icon icon = new Icon(ms);     mapControl.MouseIcon = (stdole.IPictureDisp)ESRI.ArcGIS.ADF.COMSupport.OLE.GetIPictureDispFromIcon(icon);     mapControl.MousePointer = ESRI.ArcGIS.Controls.esriControlsMousePointer.esriPointerCustom;       This code doesn't call errors during the first time only. After checking SelectFeatureTool on ToolbarControl and pushing my custom zoomin botton, the above script doesn't work anymore. The cursor does not change. The 'SelectFeatureTool' and 'Zoomin' still work with 'SelectFeatureTool' cursor. How to create change cursor?


title: Google Earth alternate tile set

text: I have a Google Maps compatible tileset running on a local server and I can browse this in a browser with the Google Maps API and OpenLayers. Can I force Google Earth to use this tileset so I can navigate it in 3D or will it be incompatible? If so, where do I set the function or string that takes a zoom,row,column index and turns that into a URL? Is there some KML format somewhere?


title: How to change overlays as radio button and accordian list in Layer Switcher Control?

text: Could you please guide me how can I change the default properties of base layer and overlays to check box(Must be always be checked when I play with radio button list of overlays) to radio button(This overlay radio button list has be like accordion type of style with static data information) Please have a look into the attached imaged. so if anyone has any hints on this requirement I'd more appreciate it. ![enter image description here](http://i.stack.imgur.com/fFx4g.png)


title: Geostatistical Analyst raster export looks completely different

text: I'm having trouble trying to re-create how a Kriged layer looks when interpolated using Geostatistical Analyst. I like the level of generalization (smoothness) of the filled contours in the "temporary" GA layer prior to export, but when I "Export to Raster" I get a jagged pixelated raster. It doesn't seem to make a difference whether I use a large (50m) or small (20m) cell size. Here's an example, with both the GA and raster layers at the same contour interval. Are there other settings I should be changing for export? Geostatistical layer (smooth contours): ![GA_20m](http://i.stack.imgur.com/Tmm3G.png) Raster export (jagged contours): ![enter image description here](http://i.stack.imgur.com/IgpAA.png)


title: R-Studio - Bare Soil Line using landsat data in R

text: I am using R-studio with `landsat` package for computing and plotting the Bare Soil Line (BSL) on landsat tif image running on my HP Pavilion i7 laptop. The codes for running the `BSL()` function is as under:               b3 <- as.data.frame(raster("LT51500371993170ISP00_B3.TIF", na.rm=TRUE))     b4 <- as.data.frame(raster("LT51500371993170ISP00_B4.TIF", na.rm=TRUE))          bsl.1993 <- BSL(b3, b4, "quantile", ulimit=0.99, llimit=0.005)      I have received the error message when I write the code without `na.rm=TRUE` as by default argument is `FALSE` for `na.rm`. After using `na.rm=TRUE` I received the same message as under. But one thing I want to mention here although my system is 64 bit machine with 8GB RAM and running R-studio 32 bit version when I execute the code for whole tile of landsat RAM reaches up to approximatively 7GB, it gives the message to running the memory low message after 5 minutes of the execution then the same message appeared the second time as under: > Error in quantile.default(ratio43, llimit) : missing values and NaN's not > allowed if 'na.rm' is FALSE Please suggest me the solution to solve this matter and how do I improve the low memory. I have also running the linux machine on VMware worksation.


title: Python Field Calculator ArcGIS 10 expression for a RANGE

text: what is the python expression for defining a range of numbers? I have tried numerous combinations but cannot get a simple range of angles (100 -200) to compute correctly. Screen shot : ![Range in Field Calculator needed](http://i.stack.imgur.com/XbM1O.jpg) screen shot here


title: writing a 3D polygon (arcpy)

text: I currently working on some scripts with 3D polygons using arcpy (arcgis10). Reading and data manupulation works fine, however writing the output back into a shapefile doesn't seem to work, it only generated 2D files. for writing i use the following function:               import arcpy         def write_new_MP(data,outfile):             array = arcpy.Array()             point = arcpy.Point()             features=[]             for feature in data:                 for part in feature:                              point.X=part[0]                     point.Y=part[1]                     point.Z=part[2]                     array.add(point)                 MP=arcpy.Polygon(array,None,True,True)                 array.removeAll()                 features.append(MP)             arcpy.CopyFeatures_management(features, outfile,"", "0", "0", "0")             return      the data is a list like`data=[[[5,5,5],[10,5,5],[10,10,5],[5,10,5],[5,5,5]],[[15,15,15],[10,15,15],[10,10,15],[15,10,15],[15,15,15]]] ` Does somebody know why the Z-values aren't used when generating the output?


title: Add AutoCAD .dwg annotation file to Arcmap using python (Same as +ADD Data button)

text: I would like to add an AutoCad .dwg file to Arcmap using python the same way that the +ADD Data button does. I can add an Annotation layer with python but only as a feature layer not as an annotation layer... even thought in the properties it says it's a CAD Annotation Feature Class               targetGroupLayer = arcpy.mapping.ListLayers(mxd, "ACAD", df)[0]     addLayer = arcpy.mapping.Layer(PlanPath)     arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM")      The files share the same properties, but are different types (Feature, Annotation)... ![enter image description here](http://i.stack.imgur.com/1O8xs.jpg) ![enter image description here](http://i.stack.imgur.com/RBoDl.jpg) NOTE: I realize in my screencaps the files have different R1, R2 numbers ignore this. The documentation for layer says: > "There are essentially three categories of layers in a map document: feature > layers, group layers, and raster layers. There are a few specialized layers > and datasets that don't fall into one of these three categories: annotation > subclasses Group layers and other sublayers (e.g., annotation classes) are > treated just like ordinary layers." Is there any way to add it as an Annotation layer? The same way the +Add Data button works? EDIT: This is what I tried doing.               targetGroupLayer = arcpy.mapping.ListLayers(mxd, "ACAD", df)[0]     arcpy.MakeFeatureLayer_management(PlanPath,"CAD")          for lyr in arcpy.mapping.ListLayers(mxd,"CAD",df):         lyr2 = lyr         #Include the Add Layer to Group code as part of the For Loop         arcpy.mapping.AddLayerToGroup(df,targetGroupLayer,lyr2)      From what I could tell from your code (Get Spatial). The annotation layer is added to a temp feature layer. You then search for that temp feature layer "cadname" and then add it to the target grouplayer.I tried to recreate this but using and for loop to loop through each layer and find the one called "CAD" @GetSpatial I ran the new code and looks like it made the feature layer yet did not correctly add it to the map. In the report below is the "CAD" after "\Annotation" part of the file name it is trying to find, or is it simply the name of the Layer and has no effect on the file path? ![Report](http://i.stack.imgur.com/0mFVS.jpg)


title: Methods for optimizing multcore processing in ArcGIS

text: I am interested in learning methods to utilize the full extent of multicore processing power available on a desktop computer. Arc states that background geoprocessing allows the user to utilize multiple cores, however, tasks essentially have to wait in line for the previous task to be completed. Has anyone developed parallel or multithreaded geoprocessing methods in Arc/Python? Are there hardware bottlenecks that prevent multicore processing on individual tasks? I found an interesting example in Stackoverflow that caught my interest, although it is not a geoprocessing example:               from multiprocessing import Pool     import numpy          numToFactor = 976          def isFactor(x):         result = None         div = (numToFactor / x)         if div*x == numToFactor:             result = (x,div)         return result          if __name__ == '__main__':         pool = Pool(processes=4)         possibleFactors = range(1,int(numpy.floor(numpy.sqrt(numToFactor)))+1)         print 'Checking ', possibleFactors         result = pool.map(isFactor, possibleFactors)         cleaned = [x for x in result if not x is None]         print 'Factors are', cleaned


title: How do I use Oracle Locator's mdsys.sdo_geometry with ArcSDE?

text: I'm using the Oracle Locator MDSYS.SDO_Geometry type to store my spatial information and then accessing it through ArcSDE. I'm confused about how the precision component of the metadata is used and more specifically how it relates to GPS accuracy. What is the maximum precision that I can use with ArcSDE? This is a sample of my code: Create sequence for object id column (oid)               CREATE SEQUENCE seq_myTable_oid      MINVALUE 0      MAXVALUE 99999999      START WITH 1      INCREMENT BY 1      CACHE 20;      Create table               CREATE TABLE myTable     (      id INT NOT NULL,      geom MDSYS.SDO_GEOMETRY NOT NULL,      oid INT NOT NULL     );     ALTER TABLE myTable     ADD CONSTRAINT pk_myTable PRIMARY KEY (id);      Create insert trigger on myTable for OID               CREATE OR REPLACE TRIGGER myTable_i       BEFORE INSERT ON        FOR EACH ROW     DECLARE     BEGIN       SELECT seq_myTable_oid.NEXTVAL INTO :new.oid FROM DUAL;     END;     /      Insert geom metadata for myTable               INSERT INTO USER_SDO_GEOM_METADATA     VALUES(     'MYTABLE','GEOM',         MDSYS.SDO_DIM_ARRAY(         MDSYS.SDO_DIM_ELEMENT('X', -180, 180, 0.00005),         MDSYS.SDO_DIM_ELEMENT('Y', -90, 90, 0.00005),         MDSYS.SDO_DIM_ELEMENT('Z', -10000, 100000, 1)          ),        4326     );      Create spatial index on table               CREATE INDEX sp_idx_myTable_geom ON myTable(geom) INDEXTYPE IS MDSYS.SPATIAL_INDEX;      Register table with SDE               sdelayer -o register -l myTable,geom -C oid,USER -t SDO_GEOMETRY -P HIGH -e l3 -s serverName -i instance -u userName -p password


title: Read KML file with GDAL/OGR for Python

text: I have file : test.kml and I try to open it with GDAL/OGR (1.6) librairie in Python (2.5). But I don't succed can you help me ? This is my code :                from osgeo import ogr          adresse = 'test.kml'     driver = ogr.GetDriverByName('KML')     datasource = driver.Open(adresse)     layer = datasource.GetLayer()     layerDefn = layer.GetLayerDefn()     featDefn = layer.GetLayerDefn()     feat = ogr.Feature(featDefn)     nbFeat = layer.GetFeatureCount()          print nbFeat      > ERROR Traceback (most recent call last): File "C:\Documents and > Settings\Guilhain\Mes documents\My Dropbox\dev\test.py", line 6, in layer = > datasource.GetLayer() AttributeError: 'NoneType' object has no attribute > 'GetLayer'


title: How to make a one part polyline from a multipart polyline?

text: I have a polyline which is divided into multiparts and these parts are not well ordered so that the vertex's ID are not following the polyline's plotting. This polyline is a representing a bus line and I need to have the plotting following the driving way of the bus. **I would like to know how to merge the multiparts of this polyline?** I have already tried to Split and Unsplit, as well as Multipart to Single Part, Union, Merge, etc. but no one is giving me the result I want. Could anybody help me? Thanks! BM


title: Can I change the GDAL/OGR library that QGIS uses?

text: I have been adding a new format to OGR, all seems well and OGRInfo works. Now I would like to test my code out using QGIS. Is there an easy way to switch which GDAL lib QGIS uses on a windows box? I've read Tim's blog post but would rather not go down the recompiling from source route if possible.


title: GeoExt printPage cannot handle WFS vector features

text: I recently added the print function to my GeoExt web application, using               <script type="text/javascript"src="http://localhost:8080/geoserver/pdf/info.json?var=printCapabilities"></script>      to define the printCapabilities variable,                       var printProvider = new GeoExt.data.PrintProvider({                 method: "GET",                  url: "localhost/geoserver/pdf",                  capabilities: printCapabilities             });      for the print provider and within the mapPanel               bbar: ["->", {                     text: "Print...",                     handler: function() {                          var printDialog = new Ext.Window({                                 items: [new GeoExt.PrintMapPanel({                                     sourceMap: mapPanel,                                     printProvider: printProvider,                                         listeners: {                                         "loadcapabilities": function() {                                         var printPage = new GeoExt.data.PrintPage({                                             printProvider: printProvider,                                             customParams: {                                                 mapTitle: "LVR WebGIS",                                                 comment: "created with GeoExt"                                             }                                         });                                             printPage.fit(mapPanel, true);                                             printProvider.print(mapPanel, printPage);                                         }                                         },                                          })                                 ],                              bbar: [{                                     text: "Create PDF",                                     handler: function() {                                         printDialog.items.get(0).print();                                     }                                 }]                         });                     printDialog.show();                     }                 }],      to add the function to my mapPanel. Regarding this, I have two issues to be solved: 1.) the printPage popping up does not contain anything but the button saying "Create PDF". There should be a preview of the map to print, right? 2.) when some WFS vector features are displayed in the map, the print function does not work properly anymore. It opens a new tab in my browser but does not create the PDF. Instead, the new tab is left about:blank. What am I doing wrong here? Is the information I provided enough for you guys to help me solve the problem? Thanks a lot for your help!


title: Kriging with barriers in ArcGIS Desktop?

text: I am trying to do interpolation in a coastal area with data points in several rivers emptying into a bay and on opposite sides of peninsulas. I have created a coastline layer to use it as a barrier but cannot figure out how to do this in ArcGIS 10. I do not see a barriers object anywhere in any of the kriging tools. I tried just doing the kriging then clipping it to the land, but I get unrealistic results where data points close together but separated by land create problems. I see online documentation for doing this in geoprocessing scripts in 9.3 but nothing similar for 10.


title: How to connect points with a polyline in ArcMap 10?

text: I have a layer containing a few hundred points representing a survey track. The points are laid out in chronological order in the corresponding attribute table. How can I have ArcMap 10 automatically "connect the dots" so that I essentially turn the points into a polyline representation of my survey track? Thanks! BTW, is there a way for me to specify a few of those points for ArcMap to IGNORE when it makes the polyline?


title: How to calculate percentages using QGIS field calculator expression?

text: ![Here is an example of the table i am editing](http://i.stack.imgur.com/vZtwc.jpg)I am new to QGIS. I am wanting to calculate a field using other fields and some operators. Here is my expression:               case       When "RVM_CAT" =  'B' OR "RVM_CAT" = 'C' Then ("PC1" /100) * "RE_ha"       When "RVM_CAT" = 'R' Then ("PC1_1" /100) *  "RE_ha"       When "RVM_CAT" = 'X' AND "FPC2011"  = 'gr12' Then ("PC1_1" /100) *  "RE_ha"       ELSE 0       end        When it runs, it looks like its only calculating values when the PC1_1 value or the PC1 value is 100. If I take out the Else 0, its the same thing.   There are over 370,000 features   The PC fields are percentages that i need to use to calculate the percentage of area from the M column. The RVM_CAT column tells me if I use PC1 or PC1_1. So the answer for row 2 uses PC1_1 and is (50/100)*32.03. WHile the answer for row 3 uses PC1 and is 100/100)*254.73


title: Creating polygons from a one to many relationship of points

text: I have a set of polygons that are in a one to many relationship with points. The polygons are administrative boundaries and the points have information on companies and percentage of ownership . I need to get a summary of the company data for each polygon. The way I know how to do it is to do a relate. I did a spatial join on the points to attach the polygon information and then did a relate on the polygon table to highlight points that fall within it. From there I can do a summarize to get the values I need (sums/averages of company data within each polygon). But that process is only good for one at a time processing. I have several hundred polygons I want the information for. Any ideas? I feel like there is something obvious I am overlooking. Thanks!


title: Creating a User Defined Point Grid

text: Im trying to write some code (in C#, but any language will do) where the user defines and area, based on Min/Max Long/Lat's, and specifies a distance between points in meters, and the code returns a table of long/lats (WGS84/ESPG:4326) to create a grid. I wrote some code a while ago that I was under the impression did the job, but after a quick bit of analysis, doesnt _quite_ measure up to the job. Please see the code below:               decimal minX = -74.002747535706M;     decimal minY = 40.722282672831M;     decimal maxX = -73.98386478424M;     decimal maxY = 40.733470232685M;     decimal distance = 100;     decimal x = 0.0M;     decimal y = 0.0M;     decimal length = 0.0M;     decimal lengthIncrease = 0.0M;          decimal increasePercentage = 0.0M;     int i = 0;          decimal yIncrement = distance / (decimal)111111.111111;          DataTable dt = new DataTable();     dt.Columns.Add("Longitude" , Type.GetType("System.Decimal"));     dt.Columns.Add("Latitude" , Type.GetType("System.Decimal"));          y = minY;          while (y < maxY)     {         while (x < maxX)         {             length = (decimal)111.325 * (decimal)Math.Cos((double)y * (double)0.0174532925199433);             increasePercentage = (decimal)0.001 / length * 100;             lengthIncrease = increasePercentage / 100 * distance;             x += lengthIncrease;             i++;             DataRow dr = dt.NewRow();             dr["Longitude"] = x;             dr["Latitude"] = y;                  dt.Rows.Add(dr);         }         y += yIncrement;         x = minX;     }          return dt;      The above works fairly well. It returns long/lats seperated by 100m lat, but by 90m long. When i change the `distance` to 25m, it seperates latitude by 25m, but longitude by 22.3m. **EDIT** : C Sharp Version Solution Based on yosukesabai Answer:               int radius = 6378137;     double pi = Math.PI;     double deg2Rad = pi / 180;     double dst2Lat = 360 / (2 * pi * radius);          DataTable dt = new DataTable();     dt.Columns.Add("ID" , typeof(int));     dt.Columns.Add("Longitude", Type.GetType("System.Decimal"));     dt.Columns.Add("Latitude", Type.GetType("System.Decimal"));          x = minX;     y = minY;     while (y < maxY)     {         double rLat = radius * Math.Cos(deg2Rad * y);         double dst2Lon = 360 / (2 * pi * rLat);         while (x < maxX)         {             i++;             x += dst2Lon * distance;             DataRow dr = dt.NewRow();             dr["ID"] = i;             dr["Longitude"] = x;             dr["Latitude"] = y;                  dt.Rows.Add(dr);         }         x = minX;         y += dst2Lat * distance;      Any help here would be fantastically helpful! Regards AM


title: Python script to check node / edge coincidence

text: I have written a script with arcpy to test node and edge coincidence. Each edge (line) should have a node (point) at either end. If this is true, the script does not write anything. If there is node at one end it is marked as closed; if there is not a node at the other end it is marked as open. The script works, but I'd like to optimize it. For 100 edges it is fine, but my dataset has 89,000 edges and it takes over an hour to run:               import arcpy     import os          node_featureclass = arcpy.GetParameterAsText(0)     edge_featureclass = arcpy.GetParameterAsText(1)     folder = arcpy.GetParameterAsText(2)     filename = arcpy.GetParameterAsText(3)               openfile = open(folder + (os.sep) + filename + '.csv', 'a')     openfile.write('PIPEobjectID, ERRORtype1, ERRORtype2' + '\n')          nodeDict = {}     for row in arcpy.da.SearchCursor(node_featureclass, ["SHAPE@XY"]):         nodeXY = row[1]         nodeDict[nodeXY] = row[0]               for row in arcpy.da.SearchCursor(edge_featureclass,["OID@", "SHAPE@"]):         edgeOID = row[0]         firstPt = row[1].firstPoint         lastPt = row[1].lastPoint         firstPointXY = firstPt.X, firstPt.Y         lastPointXY = lastPt.X, lastPt.Y               for key in nodeDict:         if nodeDict.has_key(firstPointXY):             upStreamMH = "CLOSED"         else:             upStreamMH = 'OPEN'              if nodeDict.has_key(lastPointXY):             downStreamMH = "CLOSED"         else:             downStreamMH = 'OPEN'          if upStreamMH == "CLOSED" and downStreamMH == "CLOSED":        pass     else:         openfile.write('%i'% edgeOID + ',' + upStreamMH + ',' + downStreamMH + '\n')          openfile.close()      Any suggestions on how it might be re-written to make it run faster? Thanks in advance.


title: Using MakeFeatureLayer to select by location in arcpy

text: I'm trying to assign a state plane zone based on if a point intersects it. I make a feature layer from a feature class. It gets through the loop once and calculates everything correctly, but when it loops again it says that the dataset already exists. Can I put the MakeFeatureLayer command at the top and have it loop through underneath that? This is the line that it gets stuck on after it calculates the first group of points:   arcpy.MakeFeatureLayer_management(fc,"fcTemp")                   for fc in SMON:             print fc + " Started Processing"             arcpy.AddField_management(fc, 'LATDMS','text')             arcpy.AddField_management(fc, 'LONGDMS','text')             arcpy.AddField_management(fc, 'DESC_','text')             arcpy.AddField_management(fc, 'StatePlane','text')             arcpy.AddField_management(fc, 'State','text')             arcpy.AddField_management(fc, 'Unique_ID','text')             arcpy.AddField_management(fc, 'DateImport','text')             print "Fields added"             arcpy.CalculateField_management(fc, "LATDMS", "decdeg2dms ( !GlobalLati! )", "PYTHON", "def decdeg2dms(dd):\\n    negative = dd < 0\\n    dd = abs(dd)\\n    minutes,seconds =  divmod(dd*3600,60)\\n    degrees,minutes =   divmod(minutes,60)\\n    if negative:\\n        if degrees > 0:\\n            degrees = -degrees\\n        elif minutes > 0:\\n            minutes = -minutes\\n        else:\\n            seconds = -seconds\\n            seconds = round(seconds, 5) #my attempt to define it\\n   \\n    dms = '{0:.0f} {1:.0f} {2:.5f}'.format(degrees, minutes, seconds)\\n    \\n    return (dms)")             arcpy.CalculateField_management(fc, "LONGDMS", "decdeg2dms ( !GlobalLong! )", "PYTHON", "def decdeg2dms(dd):\\n    negative = dd < 0\\n    dd = abs(dd)\\n    minutes,seconds =  divmod(dd*3600,60)\\n    degrees,minutes =   divmod(minutes,60)\\n    if negative:\\n        if degrees > 0:\\n            degrees = -degrees\\n        elif minutes > 0:\\n            minutes = -minutes\\n        else:\\n            seconds = -seconds\\n            seconds = round(seconds, 5) #my attempt to define it\\n   \\n    dms = '{0:.0f} {1:.0f} {2:.5f}'.format(degrees, minutes, seconds)\\n    \\n    return (dms)")                       #Assigning State Plane and Unique ID for SMON             #fcTemp = "fcTemp"             arcpy.MakeFeatureLayer_management(fc,"fcTemp")             arcpy.SelectLayerByLocation_management("fcTemp", "INTERSECT", CON, "", "NEW_SELECTION")             urows = arcpy.UpdateCursor("fcTemp")             for urow in urows:                 urow.StatePlane = 'Colorado NorthTEST'                 urow.State = 'Colorado'                 urow.Unique_ID = str(urow.PointID) + str(Con)                 urows.updateRow(urow)                  arcpy.SelectLayerByLocation_management("fcTemp", "INTERSECT", NDN, "", "NEW_SELECTION")             urows = arcpy.UpdateCursor("fcTemp")             for urow in urows:                 urow.StatePlane = 'North Dakota North'                 urow.State = 'North Dakota'                 urow.Unique_ID = str(urow.PointID) + str(Con)                 urows.updateRow(urow)


title: Labeling in TileMill with CartoCSS - Add line break between variables

text:
I am labeling some points in TileMill using CartoCSS. The labels are a composite of two variables text-name: [NAME]+' '+[TYPE]; The rule above inserts a space between the two variables. How can I add a line break? I have tried               text-name: [NAME]+
+[TYPE];      .               text-name: [NAME]+'
'+[TYPE];      .               text-name: [NAME]+'_'+[TYPE];     text-wrap-character: '_';      All to no avail. Any suggestions?


title: How to create a polygon feature using Python?

text: I'm making a plugin for Quantum GIS in python. I have a polygon layer and want to create a new feature. So far, I have this because I need a polygon geometry to add the feature:               surface = QgsGeometry.fromPolygon (polygon)      The problem is that I do not create a polygon. I have a list of points where I want to place the polygon. Can anyone help?


title: Python/Arc10- Loop through multiple feature datasets and Intersect Files

text: I am getting an error with a python script I have created. I have 1 geodatabase with 3 feature datasets. What I would like to do, is loop through 1 feature dataset (NB) to get a list of feature classes. One by one, I would like to get the name of the feature class (ID_NB) and search for its corresponding feature class in another dataset (called CCH) so that I can perform the Intersect function on the two files and save the output in the 3rd feature dataset and call it ID_Ints. **There may not always be a matching feature dataset between the two either. The script works fine the first time through, but when it comes to the second feature class, it gives me an error saying "Dataset does not exist or is not supported". I have a feeling it is because I have switched workspaces in a for loop, but I can't think of another way around this. Anyone have any ideas? Any help would be so much appreciated. Thanks!               #Define Environments     CCHws = "C://TEST.gdb//CCH"     PerMap = "C://TEST.gdb//NB"     opws = "C://TEST.gdb//opws"          #Set workspace     arcpy.env.workspace = PerMap     arcpy.env.overwriteOutput = True     NBList = arcpy.ListFeatureClasses()          #Start forloop, find first NB Feature Class and extract the ID (the first 5 characters)     for NB in NBList:         NBname = NB[:5]         NBname2 = NBname + "*"         NBnameEx2 = "\"" + NBname2 + "\""          #Calculate Original Area of NB feature Class         arcpy.AddField_management(NB, "NB_Area", "DOUBLE")         arcpy.CalculateField_management(NB,"NB_Area","!shape.area@meters!","PYTHON")          #Set variables to use in Intersect         CCH = str(CCHws) + "//" + str(NBname) + "_CCH"          opresult = opws + "//" + str(NBname) + "_Ints"         NBfile = str(PerMap) + "//" + str(NBname) + "_NB"          #Re-set workspace to look for matching ID with an '_CCH' in the CCH dataset         arcpy.env.workspace = CCHws         arcpy.env.overwriteOutput = True         cchList = arcpy.ListFeatureClasses(NBname2)         for cch in cchList:             if len(cchList) > 0:                 arcpy.Intersect_analysis([NBfile,CCH],opresult)               else:                 print NBname + " Does not have a matching CCH"


title: TIN Creation from 3D Contour File

text: I'm currently trying to import a 3D contour file into a Quantum GIS project. My question is, how can I export the DXF file from Microstation with height information? I already have imported it to the Quantum GIS project but can not create a TIN because the height info is missing. I'm using Quantum GIS Lisboa 1.8.0


title: Raster data on PostGIS

text: I need to publish raster and vector data through GeoServer WMS service and I have already chosen PostGIS as datasource to keep vector data. Now I would like to know what is the best practice to work with raster data. I see two options:   * store raster data in PostGis and then use ImageMosaicJDBC or Postgis Raster Plug In (I read about some limitations and instability)   * store raster data somewhere on file system and setup GeoServer to read data from there. I have little experience with GIS and I wonder if someone could suggest the best choice?


title: Where are the Expression-based Labeling options in QGIS?

text: I have imported a layer into my map to show the location of towns/cities. However I had to set graduated styling so that only cities over a certain population show. I now want to label the ones that show but if I add labels all of the labels show rather than just the ones which I want to appear in my map. I think I can solve the problem using Expression-based labeling. Following this tutorial: https://underdark.wordpress.com/2011/10/27/expression-based- labeling-for-qgis/ But I can't find the settings. When I go into the labelling options I see this: ![Label options](http://i.stack.imgur.com/EYluD.jpg) What do I do next? Any help would be amazing.


How to draw a line with pop-up in OpenLayers? For markers with pop-up i do so:               var lonLat = new OpenLayers.LonLat(lon, lat).transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject());                 var markers = new OpenLayers.Layer.Markers( "Markers" );         map.addLayer(markers);              map.events.register("click", map , function(e){        var opx = map.getLonLatFromPixel(e.xy) ;        var marker = new OpenLayers.Marker(lonLat);        markers.addMarker(marker);        marker.events.register("click", marker, function(e){        popup = new OpenLayers.Popup.FramedCloud("chicken",                              marker.lonlat,                              new OpenLayers.Size(200, 200),                              "example popup",                              null, true);               map.addPopup(popup);         });      });


title: QGIS Plugin development version compatibility

text: At this time latest QGIS version is 2.4 , QT Version is 5.3.1 and PyQT5.Are these version compatible with each other ? Once QGIS 2.4 is installed , how does one know which version of QT and PyQT(and Python too..) are required to develop Python,C++ plugins and C++ standalone applications ?


title: Best way to manage large number of geoTiff rasters

text: I have around 2-3000 1-bit geoTiffs of size 11000x9400. Color 0 should be transparent, and QGis gets it wrong in the first place, so I used a .vrt file to set all the parameters in one single place. Problem is, even if I set the layer to be visible only at a small scale, the project takes a lot of time to load in QGis. Which is the best way of dealing with all these files? I tried using grass, but I don't see if it can really be faster than QGis.


title: WFS-T Namespace

text: I've faced with this problem: The WFS POST request that OpenLayers performs does not contains following namespaces: xmlns:wfs="hxxp://www.opengis.net/wfs" xmlns:ogc="hxxp://www.opengis.net/ogc" xmlns:gml="hxxp://www.opengis.net/gml" xmlns:xsi="hxxp://www.w3.org/2001/XMLSchema-instance" The weird thing is that these namespaces are mentioned in test expected value at OpenLayers\tests\Format\WFST\v1_1_0.html ... STATE_NAME STATE_FIPS STATE_ABBR --> ... I have analysed this test with FireBug and get following: expected value (exp) is: STATE_NAME STATE_FIPS STATE_ABBR And real value (got) is: STATE_NAME STATE_FIPS STATE_ABBR As you can see there is no topp, wfs, ogc, xsi namespaces but test is resulted as PASSED. The fact that mentioned namespaces are missed caused following problem with the Geoserver: org.xmlpull.v1.XmlPullParserException: could not determine namespace bound to element prefix wfs (position: START_DOCUMENT seen ...www.opengis.net/wfs hxxp://schemas.opengis.net/wfs/1.1.0/wfs.xsd">... @1:139) could not determine namespace bound to element prefix wfs (position: START_DOCUMENT seen ...www.opengis.net/wfs hxxp://schemas.opengis.net/wfs/1.1.0/wfs.xsd">... @1:139) I heve reported this issue on the Geoserver users mail list but no response yet. Please if anybody have any idea on this issue help me. I've search all the internet with no luck. Maybe I've missed something? Thanks, Sunny


title: Drag and Drop with Arcengine MapControl within WindowsFormsHost in WPF

text: I have an application that requires dragging and dropping onto Arcengine's Mapcontrol in a Windowsformhost. I have set the Allowdrop property on both the host and the control itself and it still does not allow me to drop on the control. Is there a way to enable dropping on the MapControl? thanks in advance.


title: How to collect the graphics element drawn on a graphics layer?

text: The graphics element may consist of point, circles, polygon, free hand line drawing or marker symbols. I m using arcobjects.


title: changing ArcGIS Viewer for Flex page title

text: perhaps I'm just missing the obvious. when working with the ArcGIS Viewer for Flex application is it possible to change to title page tab to something else e.g. sometitlepage.com or will the Internet police come for me if I try to alter this setting? was going to post an image but cannot because I'm 'new'. link here - http://i1100.photobucket.com/albums/g417/jdsegars/junkphotos/pagetitletoremore.png?t=1311342928 thanks! John


title: Dynamic scale dependency from attribute query using ArcGIS for Desktop?

text: If I have a feature class called roads with an attribute for type e.g. highway etc. Is there a way to setup a dynamic scale dependency on the layer based on the type of road? Like do a SQL query and apply a specific scale dependency range on the records that match?


title: PostGIS 2.0 ST_Transform with Rasters cannot find custom SRID

text: I'm trying to calculate an average value in a buffer around a point. Given a raster in EPSG 4326 and a similar point, this is fine, except that the values will be skewed as they are not projected in an equal-areas projection. So - I've added to my spatial_ref_sys the SRID '97099' from saptialreference.org: http://spatialreference.org/ref/sr-org/7099/ This works fine for transforming points [ **edit** - added example ]:               -- get all the diff areas of diff pix values in one tile in which a point lies     SELECT ST_Area((ST_DumpAsPolygons(r.rast)).geom),            (ST_DumpAsPolygons(r.rast)).val,            r.rid     FROM rasters AS r     INNER JOIN (         SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-2, 51), 4326), 97099) as pt     ) AS p     ON ST_Intersects(r.rast, p.pt) AND r.filename = 'someraster.tif';      But using it with the raster fails with an error saying it cannot find SRID 97099 in spatial_ref_sys: [ **edit** - added error message ]               ERROR: rtpg_getSR: Cannot find SRID (97099) in spatial_ref_sys     SQL state: XX000     Context: SQL function "st_transform" statement 1      A `SELECT ... WHERE SRID = 97099` confirms it is there... So I'm a bit lost. Is it querying from somewhere else, or maybe ignoring anything that isn't EPSG? I've followed the example here, for rasters: http://www.postgis.org/documentation/manual-svn/RT_ST_Transform.html Here's the full test SQL, the commented out lines show what I want to do, and the line above them will actually work (but is wrong, as it will still be 4326)               DROP TABLE IF EXISTS tmp_points;     CREATE TEMP TABLE tmp_points ("gid" int, "point" geometry);     INSERT INTO tmp_points ("gid", "point") VALUES       (100231, ST_MakePoint(172.3, -43.76667)),       (100233, ST_MakePoint(171.71667, -43.91667)),        (100246, ST_MakePoint(170.92009, -44.09834)),        (100290, ST_MakePoint(170.95009, -44.44835)),       (100497, ST_MakePoint(173.28271, -42.77409)),        (100498, ST_MakePoint(171.90983, -43.9236)),       (1, ST_MakePoint(-2.136154, 49.220288))        ;               -- get the average value arount a point for a given radius     --      SELECT       gid,       filename,       lon,       lat,       CAST(SUM(ST_Area((foo.gv).geom)*(foo.gv).val)/SUM(ST_Area((foo.gv).geom)) AS decimal(9,7)) AS avgimr,       COUNT(foo.gv) as parts     FROM (       SELECT       ST_X(p.point)::NUMERIC(9, 5) AS lon,        ST_Y(p.point)::NUMERIC(9, 5) AS lat,       p.gid,       sn.filename,       ST_Intersection(         sn.rast,         -- THIS COMMENTED LINE BELOW FAILS, THE ABOVE WORKS (is wrong)         -- ST_Transform(sn.rast, 97099, 'Bilinear'),         ST_Buffer(ST_Transform(ST_SetSRID(p.point, 4326), 97099), 25000)       ) AS gv       FROM rasters sn, tmp_points p       WHERE ST_Intersects(         ST_Buffer(ST_Transform(ST_SetSRID(p.point, 4326), 97099), 25000),         sn.rast         -- THIS COMMENTED LINE BELOW FAILS, THE ABOVE WORKS (is wrong)         -- ST_Transform(sn.rast, 97099, 'Bilinear')       )       AND sn.filename = 'someraster.tif'     ) AS foo     WHERE (foo.gv).val >= 0 AND (foo.gv).val <= 10     GROUP BY filename, gid, lon, lat     ORDER BY filename, gid, lon, lat;      Many thanks


title: Create Spatial Index in PostGIS on an entire Schema

text: I have loaded a number of shapefiles using SPIT (QGIS plugin) to my PostGIS database. These layers did not have a Spatial Index created upon loading. I am wondering if there is a way to create a Spatial Index for each layer in the schema without writing a query for each layer. I am not a good PostGIS script writer, so any help would be greatly appreciated. Thanks


title: Add record in related table ArcPad

text: i have a layer form (named form_layer) containing a combobox and allows to add records into the related table, by showing a form of the related table (named form_table) When i try to add record showing first the layer form , then showing the form_table, the row is not saved when i click "ok" on the form_table (so i can't see the record in the combobox when just added despite the combobox is reloaded ) but until i go back to the form_layer and click ok. i tried to show the form a toolbar button without calling a layer form, it works great :) i dont know if there's an inner lock when displaying a form layer any suggestion will be very helpful


title: Calculating the space a feature class will use in ArcSDE?

text: Is it possible to get a good estimate on how much space a feature class will take once loaded into ArcSDE from how big they are in a FGDB? I have some feature classes in a FGDB, where some are a few GB's. I'm using Oracle/ ArcSDE 9.2 Thanks


**The Problem:** Every time I try to solve the service area around locations, ArcGIS crashes. It takes forever to load facilities and draw just a single polygon. How do I optimize the speed and enable ArcGIS to create the Service Area around my locations? **Method**   I have been following this tutorial: http://www.youtube.com/watch?v=10yAQNiDzNU, which describes how to create aservice area around points in a shape.shp   * Open map file (all of Denmark)   * Open location file (~450 locations in Denmark)   * Verify that they have the same coordinate system   * Add and calculate the following fields to the map: Length, Speedlimit, Minutes   * Create "New Network Dataset..." from the modified map file   * Settings: Model Turns "no", Model elevation "no", specific attributes "Minutes"   * Build and Add the "New Network Dataset" and the feature classes   * Network Analyst "New Service Area"   * Polygon Generation Settings: Polygon Generation "Generalized", "Overlapping", "Disks"   * Analysis Settings: Default breaks "1;2;5;10;20", Direction "away from facility", U-turns "not allowed", "Ignore Invalid locations"   * Load Facilities   * Solve **Data**   I have a polygon shp file of the roads in Denmark (available here: http://download.kortforsyningen.dk/content/fotkort10, Im using vejmidt- brudt.shp) I also have a point.shp file of ~450 locations in Denmark. Projected Coordinate System: WGS_1984_UTM_Zone_32N   Geographic Coordinate System: GCS_WGS_1984 **Ressources employed**     * ArcGIS: Configuring the Network Analyst environment   * ArcGIS: Exercise 5: Calculating service areas and creating an OD cost matrix   * StackExchange: ArcGIS 10 Network Analyst Service Area Problem/Peculiarity with Solving   * StackExchange: How to Enable DriveTime (Minutes) in Network Analyst Service Area Creation? **Problem Outline**   * How do I speed up the Service Area Generation?   * How do I optimize the map?   * When I create a service area for one location it takes at least 60 minutes, and the boundaries of the created service area goes into the water - how do I avoid that? Can anyone sheed some light on what Im doing wrong, or help me reach my goal? Any help is very much appreciated. Thank you in advance.


title: Is it possible to split a polygon inside another one in QGIS simply?

text: I can split polygon if I put the points of start and end outside the boundaries. How can I do the same if I have a part that I want to split completely inside a polygon. In Arcgis I can do it if the splitting line is self-intersecting but it doesn't work in QGIS. I'm interesting in the existenceof a simple and quick way.


title: Drawing a road network: offsetting overlapping lines (directed graph)

text: Say I have a road network defined by a table, indicating, for each arc, its ID, and its head and tail node. In other words, my network is defined as a directed graph. For example, we have two points, A and B, linked by a two-way road ; each way is thus represented by a specific arc. It gives us the following table :               Arc_ID | Head_node | To_node     Arc_01 |     A     |   B     Arc_02 |     B     |   A      Each node (A and B) has (X,Y) coordinates, so I drew the network using "LINESTRING" in WKT. The problem is, Arc_01 and Arc_02 lines are overlapping. I am trying to find a way to draw them so that they do not overlap, and that they are drawn in respect with the way traffic flows on them (right-hand driving). For that, I thought about calculating the normal vector of each arc, and forcing an offset of a certain number of pixels based on that (so that changing the scale of the map doesn't cause the lines to overlap again), but I have no idea on how to implement it, since I am fairly new to QGIS (I have version 2.0.1 installed) and programming in general.


title: How to prevent ArcGIS SOC from locking Oracle tables?

text: I have a python geoprocessing that performs CRUD operations on feature classes and oracle tables. Sometimes, when handling Oracle tables, ArcGIS SOC gets hung and locks the Oracle table. I need then to restart ArcGIS SOM service to unlock the table. Is there a way to prevent ArcGIS SOC from locking tables and getting hung?


title: Lon/lat transformation

text: First I should clarify I don't have previous experience with the field, so I don't know the technical terminology. My question is as follows: I have two weather datasets:   * The first one has the _regular_ coordinate system (I don't know if it has an specific name), ranging from -90 to 90 and -180 to 180, and the poles are at latitudes -90 and 90.   * In the second one, although it should correspond to the same region, I noticed something different: latitude and longitude were not the same, as they have another _reference point_ (in the description is called a _rotated grid_ ). Together with the lat/lon pairs, comes the following information: southern pole lat: -35.00, southern pole lon: -15.00, angle: 0.0.  I need to transform the second pair of lon/lat to the first one. It could be as simple as add 35 to the latitudes and 15 to the longitudes, since the angle is 0 and it seems a simple shifting, but I'm not sure. Edit: The information I have about the coordinates is the following http://www.ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/grids/0/ for the first one http://www.ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/grids/10/ for the second one Apparently, the second coordinate system is defined by a general rotation of the sphere "One choice for these parameters is:   * The geographic latitude in degrees of the southern pole of the coordinate system, thetap for example;   * The geographic longitude in degrees of the southern pole of the coordinate system, lambdap for example;   * The angle of rotation in degrees about the new polar axis (measured clockwise when looking from the southern to the northern pole) of the coordinate system, assuming the new axis to have been obtained by first rotating the sphere through lambdap degrees about the geographic polar axis, and then rotating through (90 + thetap) degrees so that the southern pole moved along the (previously rotated) Greenwich meridian." but still I don't know how to convert this to the first one. Any help on the subject will be much appreciated


title: Is there any way of merging two sets of tiles with gdal2tiles?

text: I've executed gdal2tiles against 2 raster coverages:   * One big raster coverage, composed by several TIFF files. Example of tile: http://dl.dropbox.com/u/6599273/319.png   * One single TIFF file, geographically placed inside the coverage above, but that wasn't procesed in the previous step. Example of tile: http://dl.dropbox.com/u/6599273/319_bis.png So, I first executed gdal2tiles against a "puzzle" of TIFF files, with a missing piece, and then I executed gdal2tiles against that single missing piece. I want to merge both tiles set. Zoom levels were 4-14, only Google Maps output. Only levels 4-12 have the problem (tiles with the same coords). Is there any way of merging tiles or should I execute gdal2tiles against the whole "puzzle"? Thanks in advance, and best regards


title: Count the number of points of a polygon

text: Actually, this is a continuation of a previous post entitled "Count the number of corners of a polygon". Although I managed to count the number of points (i.e. saying corners) of each polygon the problem is that even a polygon may be approximately rectangular actually has more than four points because of the digitization process i.e. sometimes an approximate straight line consists of more than 2 points. How can I ignore the intermediate additional points that are not really needed for my purpose? Thanks.


title: Call tools from ArcToolBox in VBA form

text: I am now looking for a code in VBA that can call Tools from ArcToolBox(The little red box) I have find a Arcobject code(I thought it's can run under VBA) but it's don't work! I have to delete the URL here.... in the end of the page it's have Arcobject code But it's don't work! can anyone tells me how to use it? thanks for the help!! @ Mike Toews Compile error: Sub or Function not defined at OpenPointDataset @artwork21 I call it from Tools>marcos>visual basic editor It's a simple compiler in GIS not a ui control I have tried the method with the page http://edndoc.esri.com/arcobjects/9.0/Samples/Geoprocessing/RunToolWithGPDispatch.htm and I try my Kriging Tool               Sub RunGPTool()      'Create the Geoprocessor object        Dim GP As Object        Set GP = CreateObject("esriGeoprocessing.GpDispatch.1")            'Set the toolbox        GP.Toolbox = "C:\Program Files\ArcGIS\ArcToolBox\Toolboxes\Spatial Analyst Tools.tbx"             'Execute tools       GP.Kriging "C:\sample.shp", _         "F30", "C:\AA\456.shp", "250", "12"                End Sub      but it's highlight error at GP.Kriging "C:\sample.shp", _ "F30", "C:\AA\456.shp", "250", "12" and says   automation error   unspecified error   how to fix it?? update 2011/06/11 10:00PM(+8) @Casey http://edndoc.esri.com/arcobjects/9.0/Samples/Geoprocessing/RunToolWithGPDispatch.htm in this page there is no "_analysis"in the sample code... and I tried GP.Kriging_sa but it's still don't work.... thanks for you support! by the way I have tried buffer it's works!!


title: Where to find a polygon database of golf courses in the US?

text: Does anyone know where I can find an accurate polygon database of all golf courses in the us?


title: How to represent Typhoon(s) wind speed using Quantum GIS / PostGIS?

text: I am not completly new to GIS (was taught ArcGIS at Uni) but I would still consider myself (especially compared to all these experts here) a newbie in QGIS or PostGIS. I want to display the windspeed of typhoons with the help of buffers along their tracks. I used PostGIS to combine the points of the typhoons to Lines ( like in this **post** ). Anyhow, if I buffer the lines, buffer segments from the same typhoon overlap or leave gaps. And I need to change that (Here is a **link to a picture** showing my problem) My idea was to create a single line for each typhoon, buffer the line and then split the buffer with the help of my point layer. I need the information on windspeed broken down to the points rather than just an average for the complete typhoon track. Is there a way to split polygons with the help of point layers? And if not, a way to somehow interpolate the buffer segments that I have right now to close gaps and cut overlapping parts? Doing this manually is not an option at I have to work with typhoon data from 1945-2013. EDIT: So... if I buffer the lines in Quantum I kind of"solved" the problem of gaps. So I "just" need to interpolate/cut the overlapping parts at the beginning/end of each segment (see picture: http://postimg.org/image/maca8nrc7/ )