title: Initializing Network dataset workspace extension VB.NET problem

text: I have an add in tool that is basically compiling routes of people travelling to work. It is a large tool ( a lot of code) and it worked a month ago. When I came back to use it yesterday it suddenly doesn't work any more. I think I've narrowed it down to the RouteClass which was designed from the Route Layer sample here http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//00010000041z000000 My code seems to crash when I initialize the routeclass. I'll post the section where I think it is going wrong.               MsgBox("6 " + startpoint.X.ToString)                          routeClass.Init()                     route1 = routeClass.SolveRoute()                     If Not route1 Is Nothing Then                              startOID = OutputDatabase1.WriteToRouteTableRowField(person_number, route1)                         route1length = route1.Length                         MsgBox("99" + route1length.ToString)                     Else                         route1length = 0                     End If                          crowflies = crowFliesDistance(startpoint, endpoint)                     MsgBox("22" + crowflies.ToString)                 End If             End If      This is only a small section of the code, the message boxes are for me to see where the code stops. The last message box that worked is MsgBox("6 " + startpoint.X.ToString) . Here is my RouteClass               Public Class RouteClass     Public Const SHAPE_WORKSPACE As String = "C:\Users\Paul\Downloads\Geo\"     Private Const INPUT_STOPS_FC As String = "StopsInput"     Private Const SHAPE_INPUT_NAME_FIELD As String = "StopsInput"     Private Const NETWORK_DATASET As String = "DTO_Network_ND"     'Public naLayerStore As INALayer          Dim workspaceFactory As IWorkspaceFactory = New ShapefileWorkspaceFactoryClass()          Dim featureWorkspace As IFeatureWorkspace = workspaceFactory.OpenFromFile(SHAPE_WORKSPACE, 0)          Dim pUID As UID = New UIDClass()     Dim networkDataset As INetworkDataset     Dim inputStopsFClass As IFeatureClass     Dim NALayer As INALayer     Dim naContext As INAContext     Dim stopsNAClass As INAClass     Dim routesFC As IFeatureClass          Dim naClassFieldMap As INAClassFieldMap     Dim naLoader As INAClassLoader     Dim blankstops As Object          Dim messages As New GPMessagesClass()          Dim naSolver As INASolver               Public Sub Init()              pUID.Value = "esriGeoDatabase.NetworkDatasetWorkspaceExtension"                   Dim workspaceExtensionManager As IWorkspaceExtensionManager = featureWorkspace         Dim datasetContainer2 As IDatasetContainer2 = workspaceExtensionManager.FindExtension(pUID)                   networkDataset = datasetContainer2.DatasetByName(esriDatasetType.esriDTNetworkDataset, NETWORK_DATASET)                              End Sub               'Create the analysis layer, load the locations, solve the analysis, and write to disk     Public Function SolveRoute() As Integer()         ' Open the feature workspace, input feature class, and network dataset         inputStopsFClass = featureWorkspace.OpenFeatureClass(INPUT_STOPS_FC)                   ' Create the Route NALayer         NALayer = CreateRouteAnalysisLayer("Route", networkDataset)         naContext = NALayer.Context         stopsNAClass = naContext.NAClasses.ItemByName("Stops")         routesFC = naContext.NAClasses.ItemByName("Routes")                   ' Load the Stops         naClassFieldMap = New NAClassFieldMapClass()         naClassFieldMap.MappedField("Name") = SHAPE_INPUT_NAME_FIELD              naLoader = New NAClassLoaderClass()         naLoader.Locator = naContext.Locator         naLoader.NAClass = stopsNAClass         naLoader.FieldMap = naClassFieldMap                   naLoader.Load(inputStopsFClass.Search(Nothing, True), New CancelTrackerClass(), 0, 0)         'New QueryFilterClass()         ' Message all of the network analysis agents that the analysis context has changed         'Dim naContextEdit As INAContextEdit = naContext         'naContextEdit.ContextChanged()              'Solve              messages = New GPMessagesClass()         Try             naSolver = naContext.Solver             naSolver.Solve(naContext, messages, New CancelTrackerClass())         Catch e As Exception                  TextFromFile1.errorwriter("Solver exception: " + e.Message, MakeRoute.person_no_global)                  If Not messages.Description Is Nothing Then                      TextFromFile1.errorwriter("Solver error: " + messages.Description, MakeRoute.person_no_global)             Else                 TextFromFile1.errorwriter("Solver error: ", MakeRoute.person_no_global)                       End If             Return Nothing         End Try                   'Save the layer to disk (disabled to improve execution speed         'SaveLayerToDisk(NALayer, System.Environment.CurrentDirectory + "\RouteGEN.lyr")         'naLayerStore = naLayer              'AREA FOR EXTRACTING THE FEATURES TRAVERSED BY THE ROUTE         Dim naTraversalResult As INATraversalResult = CType(naLayer.Context.Result, INATraversalResult)         Dim naTraversalResultQuery As INATraversalResultQuery = CType(naTraversalResult, INATraversalResultQuery)                        Dim traversalFClass As IFeatureClass = naTraversalResultQuery.FeatureClass(esriNetworkElementType.esriNETEdge)         Dim featureCursor1 As IFeatureCursor = traversalFClass.Search(Nothing, True)              Dim feature As IFeature = featureCursor1.NextFeature()         If feature Is Nothing Then                  MsgBox("no result found to be added to database")         End If              ' Loop through all the traversal result elements for that feature         Dim routeNums() As Integer = {0}         Dim i As Integer = 0         Dim totlength As Double = 0         Dim poly As IPolyline         While Not feature Is Nothing                       poly = feature.Shape             totlength += poly.Length             'This is the FID value of each feature                  ReDim Preserve routeNums(0 To i)             routeNums(i) = feature.Value(3)                  'MsgBox(feature.Value(0).ToString + " . " + feature.Value(1).ToString + " . " + feature.Value(2).ToString + " . " + feature.Value(3).ToString + " . " + feature.Value(4).ToString + " . ")                  feature = featureCursor1.NextFeature()             i = i + 1         End While         Marshal.FinalReleaseComObject(featureCursor1)         MakeRoute.resultDistance = totlength         Return routeNums          End Function          'Create a new network anlaysis layer and set some solver settings     Private Function CreateRouteAnalysisLayer(ByVal layerName As String, ByVal networkDataset As INetworkDataset) As INALayer         Dim naRouteSolver As INARouteSolver = New NARouteSolverClass()         Dim naSolverSettings As INASolverSettings = naRouteSolver         Dim naSolver As INASolver = naRouteSolver              'Get the NetworkDataset's Data Element         Dim datasetComponent As IDatasetComponent = networkDataset         Dim deNetworkDataset As IDENetworkDataset = datasetComponent.DataElement              'Create the NAContext and bind to it         Dim naContext As INAContext = naSolver.CreateContext(deNetworkDataset, layerName)         Dim naContextEdit As INAContextEdit = naContext         naContextEdit.Bind(networkDataset, New GPMessagesClass())              'Create the NALayer         Dim naLayer As INALayer = naSolver.CreateLayer(naContext)         Dim layer As ILayer = naLayer         layer.Name = layerName              'Set some properties on the the route solver interface         naRouteSolver.FindBestSequence = True         naRouteSolver.PreserveFirstStop = True         naRouteSolver.PreserveLastStop = False         naRouteSolver.UseTimeWindows = False              naRouteSolver.OutputLines = esriNAOutputLineType.esriNAOutputLineTrueShapeWithMeasure              'Set some properties on the general INASolverSettings interface         Dim restrictions As IStringArray = naSolverSettings.RestrictionAttributeNames              'Dim att As IStringArray = naSolverSettings.AccumulateAttributeNames              'att.Add("LEN_KM")         'naSolverSettings.AccumulateAttributeNames = att         naSolverSettings.RestrictionAttributeNames = restrictions              'CODE TO USE HIERARCHY, or not         naSolverSettings.UseHierarchy = MakeRoute.useHierarchy              ' Update the context based on the changes made to the solver settings         naSolver.UpdateContext(naContext, deNetworkDataset, New GPMessagesClass())              'Return the layer         Return naLayer     End Function          'Write the NALayer out to disk as a layer file.     Public Sub SaveLayerToDisk(ByVal layer As ILayer, ByVal path As String)         Try             Console.WriteLine("Writing layer file containing analysis to " + path)             Dim layerfile As ILayerFile = New LayerFileClass()             layerfile.New(path)             layerfile.ReplaceContents(layer)             layerfile.Save()             Console.WriteLine("Writing layer file successfully saved")         Catch err As Exception             ' Write out errors             Console.WriteLine(err.Message)             TextFromFile1.errorwriter("save layer error: ", MakeRoute.person_no_global)              End Try          End Sub      End Class Again, this worked a month ago. I think the problem is with routeclass.init() . The rest of the code should be fine even if it doesn't follow conventional methods. I have initialised the licence for arcgis in another part of my code.


title: Is there a TOC control available for Esri's Silverlight/WPF SDK?

text: I've got a requirement calling for a TOC control on a Silverlight/WPF application. Does anyone know of a TOC control (analogous to the TOCControl in ArcEngine) that I can hook into Esri's Map control?


title: legend groups and layer position/order qgis 1.8

text: I'm tring to programmatically create layer sequence also using groups. Simple python code:               layer = QgsVectorLayer("004cp000.shp", "004cp000", "ogr")     QgsMapLayerRegistry.instance().addMapLayer(layer)     legend = qgis.utils.iface.legendInterface()     index = legend.addGroup ("G1")     index = legend.addGroup ("G2",index)     legend.moveLayer (legend.layers()[0],index)        * How it is possible to remove a nested group ? qgis.utils.iface.legendInterface().removeGroup only works with top level groups.   * How a layer can be moved into a nested group? qgis.utils.iface.legendInterface().moveLayer only works with top level groups.   * How can I move a layer before or after a another layer inside the same group? After some research into QGis Source code I've discovered that it is not possible to move a layer into a sub group and also it is not possible to delete a soub group.               void QgsLegend::moveLayer( QgsMapLayer *ml, int groupIndex )     {       if ( !ml )         return;            QgsLegendLayer *layer = findLegendLayer( ml->id() );       if ( !layer )         return;            QgsLegendGroup *group = dynamic_cast<QgsLegendGroup*>( topLevelItem( groupIndex ) );       if ( !group )         return;            insertItem( layer, group );            emit itemMovedGroup( dynamic_cast<QgsLegendItem*>( layer ), groupIndex );     }          void QgsLegend::removeGroup( int groupIndex )     {       QgsLegendGroup * lg = dynamic_cast<QgsLegendGroup *>( topLevelItem( groupIndex ) );       if ( lg )       {         removeGroup( lg );       }     }


title: why are split-lines saved in the attribute table? (qgis 1.7.4)

text: A few months ago I splitted some lines in a shapefile with official roads because they are not accessible anymore. Now I loaded that file again into my project and the splitlines appears. Because I didn't made a legend yet, they show up. (see the picture below) ![enter image description here](http://i.stack.imgur.com/TEsAW.jpg) edited: I first selected the layer and turned on the toggle editing. Then I activated split features ![enter image description here](http://i.stack.imgur.com/8SWER.jpg). I drew a line across the point where the road had to be split. After doing all the splitting, I turned off the toggle editing mode. Is it normal that splitlines are saved in the attributetable? Can this be avoided? Because it is messing up my data.


title: Has MaxMind's GeoLiteCity IP geocoding database been reduced in accuracy?

text: I took a copy of MaxMind's free city-level geocoding database a year or so ago, and the file size was roughly 30MB. Checking their site again today (July 2012), the equivalent file is now only 17MB, approximately half the size. I cannot find any reference to this online. Does anyone know if they've reduced the quality or range of their free offering, or whether it's just an improvement in the binary encoding/layout of the data that resulted in the smaller file.


title: Dividing Earth into continent polygons

text: Is there a simple (low res is fine) polygonal data set that divides the **entirety** of Earth's surface into continents, with no overlap and no area excluded? Specifically: North America, South America, Africa, Europe, Asia, Oceania and Antarctica. I am working with a data set where most of the points are located in seas and oceans. An approximate sense of which continent all these points 'belong' to is needed. All continent shapefiles I've seen so far trace out the coast lines, however this excludes these marine locations. Even using country EEZ lines to map to countries and then to continents will leave much of the ocean's surface unclaimed by any continent. I understand that such a data set might not provide a valid definition of the continents, but in my case I require all points on Earth to be associated with a continent. Alternatively, if I must create such a data set myself, can anyone recommend a tool for doing so? **EDIT** lynxlynxlynx made the interesting suggestion of using tectonic plate data: ![](http://www.worldatlas.com/aatlas/infopage/tectonic.gif) However, for my purposes I would have to subdivide the Eurasian plate, and the ocean to the immediate west of North America would end up in the Pacific plate which is best matched with Oceania. **EDIT 2** I have a data set of lat/lng points of scuba diving locations. I wish to collate all of the dive sites at which one of my users has dived and produce a list of all the continents in which that user would reasonably think that they have dived. As one can dive anywhere on the surface of Earth, I would ideally like to cover the earth with such a diagram. I'm undecided about outliers such as Hawaii, which might be considered Oceania by location yet North America politically. I would tend to favour geographical location in such cases, but if a solution is easier the other way around, then that's fine too.


title: How to adjust the granularity of a heat map in QGIS?

text: Not sure I am using the correct terminology here, so I have tried to describe my problem in detail. Please edit as needed. I have a layer of points, a layer of roads and a heat map in a raster layer. I would like the heat map to reflect density of points in the point layer in such a way that I can zoom into a part of the city and see the colours for that area. To make the heat maps I am using the heat map plugin in QGIS 1.8. Unfortunately, my heat maps' resolutions don't seem to capture the level of granularity I am looking for. I have been fiddling with the settings but don't really know what settings are relevant in this case. Here are my points with the heat map overlaid: ![heat overlay](http://i.stack.imgur.com/9or6v.png) Here is what I get when I zoom in: ![zoomin](http://i.stack.imgur.com/ths2D.png) As you can see, then replacing points with roads, the result is not very illustrative: ![enter image description here](http://i.stack.imgur.com/YgeEd.png) I would like to obtain a result which gives a better idea of the point density at the neighbourhood level. These are my current settings: ![enter image description here](http://i.stack.imgur.com/cJBft.png) Update to respond to Underdark's comment below: I have tried lower values in the radius field like you suggest but that results in two problems:   1. maps are mostly "cold", it's almost as if the points don't pass a certain threshold to really register. I've thought increasing the radius makes it more likely QGIS considers them as forming clusters.   2. As you can see in the picture below, the result is highly pixelated when one zooms in. Far too much for what I am trying to do. ![enter image description here](http://i.stack.imgur.com/VyK8o.png)


title: Landsat Gap Fill

text: I am looking for a Landsat gap fill method for quantitative data analysis. The best I have found is a reference to a method called "single file gap fill triangulation" using in ENVI, but am unsure how this method was carried out and I couldn't find a tool in ENVI. I don't normally use ENVI, is there method or documentation for preforming this gap fill in ENVI v5.0? Or another method out there?


title: How to dock the default ObjectInspector in a new parent container?

text: Into my object inspector fever, I noticed that the default Object Inspector object does not dock anymore. I only have a TabPage control and two tabs, one of those hosting the default inspector. (check black marks) The other tabs works fine. Also, this is the code I'm using for the Inspect() method:                   public void Inspect(IEnumRow objects, IEditor Editor)         {             SetParent(_Inspector.HWND, DefaultPanelHandle);             ShowWindow(_Inspector.HWND, SW_SHOW);             _Inspector.Inspect(objects, Editor);         }      Why the default does not dock anymore? Any workarounds? ![alt text](http://i.stack.imgur.com/fB9Q1.png) Thanks!


title: Dynamic Display doesn't respect my annotations when my layers are in a BasemapLayer

text: I have 2 annotated layers which I load in a BaseMapLayerClass. When the map is static the annotations are visible. When I enable DynamicMap I lose all annotations. I would like a way to annotate layers in a BaseMapLayerClass while dynamic is enabled. An interesting note about the problem. If I don't have my layers in a BaseMapLayerClass, dynamic display does not affect my annotations. However I would like them to be in a BaseMapLayerClass because would make it much easier for me to save layers and restore them across sessions. I'm using ArcEngine C# VS2010 * * * I think the problem is somehow related to extents. This post: http://briannoyle.wordpress.com/2008/03/18/esri-dynamic- displaysome-limitations/ Led me to the limitations page http://edndoc.esri.com/arcobjects/9.2/NET/0b062497-a3fc-4c56-a688-0fb9b2fcb4ea.htm which states : _Labels will be cut off on custom layers that don't return the correct extent_ And actually it's weird. Sometimes the annotations will draw and then I'll move the map up a few pixels with the hand tool, and then they will no longer draw. Then move it back a few pixels and it will draw again. It's helpful there is a consistency there, but I haven't a clue how to implement a workaround.


title: Optimize the locations of 2 additional restrooms on a trail (ArcGIS 10.2)

text: In the image below, the line is a trail and the point is an existing restroom. I would like to find the rough locations of 2 additional restrooms that minimize the maximum distance along trail from anywhere on the trail to the closest of the 3 restrooms. The trail was GPS-collected in segments and connected with Edit Vertices in ArcMap. The dataset comprises 1 line shp for the trail and 1 point shp for the restroom. I have no idea how to proceed. I have posted a simplified version of this question in stackoverflow. Hints will be appreciated as well. Thanks. ![enter image description here](http://i.stack.imgur.com/In26E.jpg)


title: Enable dynamic layer with Arcgis 10

text: I'm creating a flexviewer application (arcgis API for flex 3.4) with application builder. I'm new with flexviewer. I want to add thematic widget. I was able to add the widget. But when I want to run the widget, an error occurred saying I need to have dynamic layer. I went googling all day..but I could not find a way to do it in arcgis 10.0. Most of them are guides on 10.1. So, is it possible to enable dynamic layer using arcgis server manager 10 or arccatalog 10? If it is possible, can someone direct me how? maybe some links? Thanks


title: How much time is required to get my osm data updated

text: I updated some roads in osm map. I updated using JOSM. But the updates are not visible in the map after 2 days. I am new to osm maps. Can anyone please tell me the possible reasons.


title: map composer grid annotation in DMS format

text: I have digitized a map using CRS EPSG:4326 in QGIS. In Map Composer I set the grid. I can't get the annotation to display as DMS (degrees-minutes-seconds) rather than DD (decimal degrees) format. What can I do to solve the problem?


title: Simplest way to add built-in commands to an custom add-in

text: Currently in _Config.esriaddinx_ I have               <Button id="my_UniqueID" class="MyWrapperClass" caption="My Caption"              image="Images\My_Image.png" ...> </Button>      and then I implement MyWrapperClass with:               protected override void OnClick()     {         UID uid = new UIDClass();         uid.Value = "esriArcMapUI.GraphicsToFeaturesCommand";              ICommandItem cmd = ArcMap.Application.Document.CommandBars.Find(uid);         if (cmd != null)             cmd.Execute();     }      I'd like to just have the following in _Config.esriaddinx_               <Button id="esriArcMapUI.GraphicsToFeaturesCommand" />      and have the built-in command object do all the work for me. Is this possible?


title: how to implement genetic algorithm in ArcGIS?

text: I want to implement genetic algorithm to find optimum location for placement problem with ArcGIS. Can any body explain it (step by step ) or any resources about that ?


title: How to connect QGIS web client with Geoserver?

text: i installed qgis server and qgis web client. and i can see localhost/web/site/qgiswebclient.html?map=/var/www/web/projects/helloworld.qgs&format=image/png&visibleLayers=Hello,Country ![enter image description here](http://i.stack.imgur.com/JlOi5.png) And i installed postgresql, postgis, and geoserver. i connected import shape files and i see localhost:8080/geoserver/GIStest/wms?service=WMS&version=1.1.0&request=GetMap&layers=GIStest:Country&styles=&bbox=-176.248494673366,-67.5929956881973,182.893410989597,83.6210855559975&width=783&height=330&srs=EPSG:4326&format=application/openlayers ![enter image description here](http://i.stack.imgur.com/RGt0s.png) **How to read qgis web client from geoserver wms?**


title: Syling Leaflet Maps

text: Is there anyway to style base maps using Leaflet? I know there are quite a few options for loading different pre-made styles (such as Stamen, Mapbox, etc...), but I am looking for a way to create my own style without going back to Google Map's API. Or is there a way to style one of the particular base maps? Basically, I'd like to keep using Leaflet for this project as it is quite nice and simple, but without that ability I'll have to move on. Thanks. Edit: Is there a way, through Leaflet, to use Google Maps custom styling? I like the coding aspects of Leaflet as compared to GMaps, but still would like the flexibility in terms of design that custom styling of the tiles gives me. As this is a very small project, I don't want to set up my own Tile Map Server if possible.


title: Arcgis 10.2.2 soap API Error {401 - Unauthorized}

text: I'm trying to consume secured map service via arcgis 10.2.2 soap api. I've configured my web adapter to negotiate windows domain security and everything looks fine when I check over the rest api. However I get a this exception when I try to consume the same service with soap api. I'm getting 401-unauthorized at getDefaultMapName method.                mapServer = new MapServerBindingStub(uri, username, password);      defaultMapName = mapServer.getDefaultMapName(); //<error here 401-unauthorized>      Error               Caused by: com.esri.arcgisws.runtime.exception.ArcGISWebServiceException: HTTP s     tatus code {401 - Unauthorized} received from server             at com.esri.arcgisws.runtime.transport.http.HttpClientConnectionChannelI


title: Is there a ModestMaps or other flex library for finding the distance between 2 lat/lon coordinates?

text: I'm using ModestMaps in my Adobe Flex application. Is there a flex library method to calculate the distance between two of my lat/lon coordinates, or am I stuck coding the Haversine by hand? I've looked through the API but couldn't find any distance-calculating functions.


title: Passing the value of runtime variable in WMS request

text: I'm trying to make runtime variable inside the where clause of my sql query like "select geom from table1 where sname = '%s%'. but when I'm trying to pass the value of s during WMS request (from the browser over the URL) I'm getting empty result, anyone knows how to make the correct run-time variable definition and WMS request in my case, Thank you!


title: Arcv2CAD for QGIS plugin?

text: I'm wondering if there's a QGIS plugin to do what Arcv2CAD does: export a shape file to DXF with options to include attributes as CAD object attributes or a layer of labels. My current experience with the DXF export option (Layer --> Save as… AutoCAD DXF) results in only vector lines without any data, as mentioned in thread How to export a layer to DXF from QGIS?


title: How do I insert a DATAFRAME_ELEMENT value into a dynamically generated shape file containing its extent?

text: This is for ArcGIS 10, Python 2.65. I have code which captures the extent of dataframes within all mxd's in a directory. I also want it to add the dataframe name, scale and the name of the mxd it extracted the extent polygon from. The current code that I have give the following errors. > P:\2011\Job_031_TownPlanning_SeriesProduction\Working\mxd\1 > P:\2011\Job_031_TownPlanning_SeriesProduction\Working\mxd\1\temp_extents > Created...Polygon_Extent_1 Irvinebank I Upto here the code is fine. > Executing: AddField > P:\2011\Job_031_TownPlanning_SeriesProduction\Working\mxd\1\temp_extents\Polygon_Extent_1.shp > ExtentID TEXT # # # # NON_NULLABLE NON_REQUIRED # Start Time: Thu May 12 > 08:49:57 2011 Adding ExtentID to > P:\2011\Job_031_TownPlanning_SeriesProduction\Working\mxd\1\temp_extents\Polygon_Extent_1.shp... > Succeeded at Thu May 12 08:49:57 2011 (Elapsed Time: 0.00 seconds) It creates the column but I want it to put in the x'th value from the array of dataframe names and the test print statement only produces either the full list or the x'th character in the list. > Filling in Traceback (most recent call last): File > "P:\2011\Job_031_TownPlanning_SeriesProduction\Working\mxd\1\extent_creation4.py", > line 64, in row.strField = str(x) File "C:\Program > Files\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects_base.py", line 35, in > **setattr** return setattr(self._arc_object, attr, ao) RuntimeError: ERROR > 999999: Error executing function. These are errors in writing to the row. I also want it to capture the mxd file name The updated code is below.                   import arcpy, glob, os         from arcpy import env         from arcpy import mapping         env.overwriteOutput = True              path = os.getcwd() # Script in same directory as files being processed         path2 = path + r'\temp_extents'         mxdList = glob.glob(path + "\*.mxd")         #env.workspace = path + '\extents' # Directory as files being processed              print 'Working in ' + path         #os.remove(path2 + r'\*.*')         #os.rmdir(path2)         os.mkdir(path2)         env.workspace = path2         #print env.workspace              x = 0         y = 1         z = 1              while y < 7:              print 'Processing ' + mxd          for mxd in mxdList:              mxd2 = mapping.MapDocument(mxd)              try:             dataframe = mapping.ListDataFrames(mxd2, "*")[x]             frameExtent = dataframe.extent             print frameExtent             XMAX = frameExtent.XMax             XMIN = frameExtent.XMin             YMAX = frameExtent.YMax             YMIN = frameExtent.YMin             pnt1 = arcpy.Point(XMIN, YMIN)             pnt2 = arcpy.Point(XMIN, YMAX)             pnt3 = arcpy.Point(XMAX, YMAX)             pnt4 = arcpy.Point(XMAX, YMIN)             array = arcpy.Array()             array.add(pnt1)             array.add(pnt2)             array.add(pnt3)             array.add(pnt4)             array.add(pnt1)             polygon = arcpy.Polygon(array)                  #Create Polygon Shapefile             arcpy.CreateFeatureclass_management(path2, "Polygon_Extent" + "_" + str(z), "POLYGON")             arcpy.CopyFeatures_management(polygon, "Polygon_Extent" + "_" + str(z))             print 'Created...' "Polygon_Extent" + "_" + str(z)               #Capture and Enter ExtentID, File Name and other info             for el in arcpy.mapping.ListLayoutElements(mxd2, "DATAFRAME_ELEMENT"):              ExtentIDStore = str(el.name)              #print x, y, z              print ExtentIDStore #used to see which item to add to the shape file              print ExtentIDStore[x] # should show the full entry at storage x but instead shows the character at x                   dfStore = path2 + "\Polygon_Extent" + "_" + str(z) + ".shp"              #open (dfStore, 'w')              arcpy.AddField_management(dfStore,'ExtentID','TEXT') #should only be done once per shp              arcpy.AddField_management(dfStore,'FileRef','TEXT') #should only be done once per shp              arcpy.AddField_management(dfStore,'Scale','SHORT') #should only be done once per shp              #print arcpy.GetMessages()              # code to add ExtentID value to shp file              #rows = arcpy.InsertCursor(dfStore)              #print 'Filling in ', rows              #row = rows.newRow()              #row.strField = str(x)              #rows.insertRow(row)              #print arcpy.GetMessages()                #END - Capture and Enter ExtentID, File Name and other info                  z = z + 1              except IndexError:             pass     x = x + 1     y = y + 1                 list = []            lstFCs = arcpy.ListFeatureClasses("Polygon_Extent*")       print 'Merging Polygon_Extents* to Extent.shp'            for fc in lstFCs:     list.append(fc)            arcpy.Merge_management(list, "Extent")       print 'Created Extent.shp and exiting'            #for item in list:       #arcpy.Delete_management(item)      * * * In response to the answer by Nathan. Yes extentID=name of Dataframe. I am still having problems with writing to the shape file. I thought                        rows = arcpy.InsertCursor(dfStore, mxd2)              row = rows.newRow()              print 'Filling in ', rows              row.ExtentID = ExtentIDStore              filepath = mxd2.filePath # Gets full path of MXD used in your dataframe loop              filename = os.path.basename(filepath) # Extracts the file name and extension              row.FileRef = filename               scale = el.scale              row.Scale = scale      would work but the values aren't inputed into the shape file. What am I not understanding here? best, \---- UPDATED ISSUE The code is now creating the desired output (of file name, extent, scale etc) but it seems to be overwriting the row which holds the polygon attribute -rather than appending the scale, extent etc values into the appropriate columns in the row. I've tried placing the code for this in different loops but it makes no difference. Thanks for any advice. Output > > > > Reading 2 mxd files from > Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents Processing > Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents\10_TownPlanning_B&L_Overlay_Ver2b_Kuranda.mxd > Created...Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents\t10\Polygon_Extent_1.shp > desc = Name: Shape Type : True Scale   >  : False Number of fields: 1 Filled in Koah > 10_TownPlanning_B&L_Overlay_Ver2b_Kuranda.mxd 40000.0 into > Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents\t10\Polygon_Extent_1.shp > desc = Name: Shape Type : True Scale   >  : False Number of fields: 1 Filled in Overview > 10_TownPlanning_B&L_Overlay_Ver2b_Kuranda.mxd 750000.0 into > Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents\t10\Polygon_Extent_1.shp > desc = Name: Shape Type : True Scale   >  : False Number of fields: 1 Filled in Kuranda > 10_TownPlanning_B&L_Overlay_Ver2b_Kuranda.mxd 40000.0 into > Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents\t10\Polygon_Extent_1.shp > Processing > Q:\Software\GIS\ESRI\ArcGIS10\CustomScripts\temp\Extents\9_TownPlanning_B&L_Overlay_Ver2b_MtMalloy.mxd > ... Merging Polygon_Extents* to Extent.shp Deleting identical entries > Created Extent.shp and exiting                   import arcpy, glob, os     from arcpy import env     from arcpy import mapping     env.overwriteOutput = True          path = os.getcwd() # Script in same directory as files being processed     path2 = path + r'\t6'     mxdList = glob.glob(path + "\*.mxd")     count_Mapdocs = len(mxdList)     #env.workspace = path + '\extents' # Directory as files being processed          print 'Reading '+str(count_Mapdocs)+ ' mxd files from ' + path     os.mkdir(path2)     env.workspace = path2               x = 0     y = 1     z = 1          while y < 7:              for mxd in mxdList:                  print 'Processing ' + mxd             mxd2 = mapping.MapDocument(mxd)                  try:                 dataframe = mapping.ListDataFrames(mxd2, "*")[x]                 frameExtent = dataframe.extent                 XMAX = frameExtent.XMax                 XMIN = frameExtent.XMin                 YMAX = frameExtent.YMax                 YMIN = frameExtent.YMin                 pnt1 = arcpy.Point(XMIN, YMIN)                 pnt2 = arcpy.Point(XMIN, YMAX)                 pnt3 = arcpy.Point(XMAX, YMAX)                 pnt4 = arcpy.Point(XMAX, YMIN)                 array = arcpy.Array()                 array.add(pnt1)                 array.add(pnt2)                 array.add(pnt3)                 array.add(pnt4)                 array.add(pnt1)                 polygon = arcpy.Polygon(array)                      #Create Polygon Shapefile                 arcpy.CreateFeatureclass_management(path2, "Polygon_Extent" + "_" + str(z), "POLYGON")                 arcpy.CopyFeatures_management(polygon, "Polygon_Extent" + "_" + str(z))                      dfStore = path2 + "\Polygon_Extent" + "_" + str(z) + ".shp"                 print 'Created...' + dfStore                 arcpy.AddField_management(dfStore,'ExtentID','TEXT') #should only be done once per shp                 arcpy.AddField_management(dfStore,'FileRef','TEXT') #should only be done once per shp                 arcpy.AddField_management(dfStore,'Scale','FLOAT') #should only be done once per shp               #Capture and Enter ExtentID, File Name and other info                 for el in arcpy.mapping.ListLayoutElements(mxd2, "DATAFRAME_ELEMENT"):                  desc = arcpy.Describe(dfStore)                  print "desc =" + str(desc)                  indexList = arcpy.ListIndexes(dfStore)                      # Iterate through the list of indexes                 # This is just to check on whats going on.                  for index in indexList:                     # Print index properties                     #                     print "Name: " + index.name                     print "\tType            : " + str(index.isAscending)                     print "\tScale           : " + str (index.isUnique)                     print "\tNumber of fields: " + str(len(index.fields))                  rows = arcpy.InsertCursor(dfStore, desc)                  row = rows.newRow()                  ExtentIDStore = el.name                  row.ExtentID = ExtentIDStore                  filepath = mxd2.filePath # Gets full path of MXD used in your dataframe loop                  filename = os.path.basename(filepath) # Extracts the file name and extension                  row.FileRef = filename                  scale = el.scale                  row.Scale = float(scale)                                   # This is just to check on whats going on.                  print 'Filled in ', str(ExtentIDStore), str(filename), str(scale) + " into " + dfStore                        rows.insertRow(row) # insert value into the row (should be the same row as                  z = z + 1             except IndexError:                 pass              x = x + 1         y = y + 1          list = []          lstFCs = arcpy.ListFeatureClasses("Polygon_Extent*")     print 'Merging Polygon_Extents* to Extent.shp'          for fc in lstFCs:         list.append(fc)          arcpy.Merge_management(list, "Extent")     print 'Deleting identical entries'     arcpy.DeleteIdentical_management("Extent.shp", ["SHAPE"])     print 'Created Extent.shp and exiting'          #for item in list:     # arcpy.Delete_management(item)


title: What is GDAL equivalent to ArcGIS Lookup (Spatial Analyst) tool?

text: I use the **Lookup (Spatial Analyst)** tool to create several raster layers from a combined raster layer. I know ArcPy can do the same thing. However, I would like to know whether GDAL (or another open source python package) has the same lookup function?


title: Get coordinates from address

text: I am writing a application which allow users to save a address into database (db) and in the next step it shows all addresses on a map. But in the saving step I need automatically calculate a GPS coordinates and save them into db too. But I don't know how calculate them because I don't want to show map in saving step so I can't use a Google Maps API.               Note: the Geocoding API may only be used in conjunction with a Google map;      geocoding results without displaying them on a map is prohibited.      Is there a similar service without the restriction?


title: How to convert LAS to GRID/Raster for ArcGIS input?

text: I was wondering if there is an option to convert LAS (LIDAR data) to raster or GRID format that ArcGIS could read it. Thanks.


title: How can I associate a consistent polygon area with a set of attributes using a topology?

text: Consider the following topology.               SELECT CreateTopology('topotestschema', 0, 0.01);     CREATE TABLE topotest (id SERIAL PRIMARY KEY, data TEXT);     SELECT AddTopoGeometryColumn('topotestschema', 'public', 'topotest', 'topogeom', 'polygon');     INSERT INTO topotest (data, topogeom) VALUES ('outer', toTopoGeom('POLYGON((0 0,0 10,10 10,10 0,0 0))'::geometry, 'topotestschema', 1, 0.01));      When we look at what's in the topology:               SELECT *, ST_AsText(topogeom::geometry) FROM topotest;      we get               ╔════╦═══════╦═══════════╦═══════════════════════════════════════════╗     ║ id ║ data  ║ topogeom  ║                 st_astext                 ║     ╠════╬═══════╬═══════════╬═══════════════════════════════════════════╣     ║  1 ║ outer ║ (3,1,1,3) ║ MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0))) ║     ╚════╩═══════╩═══════════╩═══════════════════════════════════════════╝      But this gets changed around when I add more shapes:               INSERT INTO topotest (data, topogeom)     VALUES     ('right', toTopoGeom('POLYGON((9 1,11 1,11 9,9 9,9 1))'::geometry, 'topotestschema', 1, 0.01))     ,('inner', toTopoGeom('POLYGON((3 3,3 7,7 7,7 3,3 3))'::geometry, 'topotestschema', 1, 0.01))     ;      Now the layer is:               ╔════╦═══════╦═══════════╦══════════════════════════════════════════════════════╗     ║ id ║ data  ║ topogeom  ║                      st_astext                       ║     ╠════╬═══════╬═══════════╬══════════════════════════════════════════════════════╣     ║  1 ║ outer ║ (3,1,1,3) ║ MULTIPOLYGON(((10 9,10 1,10 0,0 0,0 10,10 10,10 9))) ║     ║  2 ║ right ║ (3,1,2,3) ║ MULTIPOLYGON(((10 9,11 9,11 1,10 1,9 1,9 9,10 9)))   ║     ║  3 ║ inner ║ (3,1,3,3) ║ MULTIPOLYGON(((3 3,3 7,7 7,7 3,3 3)))                ║     ╚════╩═══════╩═══════════╩══════════════════════════════════════════════════════╝      How can I keep it from modifying the polygon associated with existing rows when I add new ones? Please feel free to say something if it seems like I'm on entirely the wrong track here. I'm basically experimenting with a topology and trying to figure out if it's useful for what I'm doing.


title: How to Avoid Object ID conflicts in Versioned geodatabase

text: We have just implemented a versioned Geodatabase and our schema has few feature classes and relationship classes. We do not Reconcile and Post to the parent ( Default) version, as we need some intermediate monthly job extraction and do not want to mix up jobs of past month. Our concern is , since we do not Post to default on regular base, is there a chance that Object Ids be duplicated in the monthly versions as the same as Default ? In that case wont we lose features once we post to the default at the end of the month ? Or will the Oracle/ArcSDe take care of the Object ids for all versions , so that the basetable( Default) and the delta tables( Versions) does not have features with same Object ID in a feature class ?


title: Convert 3D WKT to 2D Shapely Geometry

text: I have the WKT of many linestrings which contains three dimensions, however my script only requires a two dimensional linestring and I therefore want to discard the third dimension. My current code, using Shapely, does flatten from 3D to 2D, but it relies on the geometry's wkt method outputting as 2D (even if it is 3D). As this Shapely behaviour may change in the future, can anyone recommend a more robust method of flattening to 2D?   * Python 2.7.13   * Shapely 1.2.17   * Ubuntu 12.10                >>> from shapely.geometry import LineString       >>> from shapely.wkt import loads         >>> line_3d_wkt = "LINESTRING (1 2 0, 4 3 0, 8 9 0)"        >>> line_3d = loads(line_3d_wkt)     >>> line_3d.has_z     True     >>> list(line_3d.coords)     [(1.0, 2.0, 0.0), (4.0, 3.0, 0.0), (8.0, 9.0, 0.0)]       >>> line_3d.wkt  # outputs as 2D even if 3D geometry        'LINESTRING (1.0000000000000000 2.0000000000000000, 4.0000000000000000 3.0000000000000000, 8.0000000000000000 9.0000000000000000)'      My current technique to flatten the 3D geometry to 2D is:               >>> line_2d = loads(line_3d.wkt)     >>> line_2d.has_z     False     >>> list(line_2d.coords)     [(1.0, 2.0), (4.0, 3.0), (8.0, 9.0)]


title: HTML framework in GIS development version

text: I have a question. I need to insert an HTML frame in the QGIS developer version. That works very well. Now I will also create a title page as HTML. Now my question is which program I best create the HTML that it rejects me no formatting and lines? I've already tried it with Word and Photoshop. Both times it was me but the formatting and the lines completely discarded. Thank you for your help.


title: Adding multiple selectFeature Controls for multiple layers in OpenLayers

text:
I am trying to add two different selectFeature controls for two different layers using OpenLayers. I would like to   1. display a popup with information about the feature when **_hover_** over a feature in **vectorlayer1**   2. display a dialog with facilities to do specific operations on feature when a feature is **_clicked_** in **vectorlayer2**   3. display another dialog on right click over features in **vectorlayer2** I have added my code below and the select feature seems to work only for one layer. It was found working when i tried the same select feature for both layers but it doesnt seem to work when I tried different select feature conrol for both layers. Here is my code                       var vectorlayer1 = new OpenLayers.Layer.Vector("layer1 Overlay",{                 eventListeners:{                     'featureselected':function(evt){                         var feature = evt.feature;                         var popup = new OpenLayers.Popup.FramedCloud("popup",                         OpenLayers.LonLat.fromString(feature.geometry.toShortString()),                         null,                         "
" + feature.attributes.no+"
",                         null,                         true                     );                         feature.popup = popup;                         map.addPopup(popup);                     },                     'featureunselected':function(evt){                         var feature = evt.feature;                         map.removePopup(feature.popup);                         feature.popup.destroy();                         feature.popup = null;                     }                 }             });                  var vectorlayer2 = new OpenLayers.Layer.Vector("layer2 Overlay",{                 eventListeners:{                     'featureselected':function(evt){                         var feature = evt.feature;                         layer2_id = feature.attributes.layer2_id;                              if(layer2_id!=''){                                                  var element = document.getElementById('layer2Info');                                           var out = "";                                  out += '/*My html code with different controls for opereations on features goes here*/';                             element.innerHTML = out;                                                     $( "#layer2Info" ).dialog();                              }                          },                     'featureunselected':function(evt){                         $( "#dialog" ).dialog('destroy');                         jQuery('#layer2Info').fadeOut("slow");                         $( "#layer2Info").dialog('close');                          }                 }             });                           var hoverSelector = new OpenLayers.Control.SelectFeature([vectorlayer1],{                     hover:true,                     autoActivate:true                 });                      var clickSelector = new OpenLayers.Control.SelectFeature([vectorlayer2],{                     click:true,                     autoActivate:true                 });                                                                                map.addControl(hoverSelector);                 map.addControl(clickSelector);      Only clickSelector for vectorlayer2 is working now. I havnt written code for right click yet. But please tell me how right click is handled in openLayers. Awaiting you reply soon. Please Please Please help!


title: WFS describeFeatureType results with curl/PHP

text: sending the following request in the browser works OK for me:               http://localhost:8080/geoserver/wfs?request=DescribeFeatureType&typeName=topp:states      but I'd need its output in a PHP-variable to parse/filter the result. However, when I send the same URL using curl, I get nothing. No error-message, nothing visible in geoservers logs. just nothing:               $url = "http://localhost:8080/geoserver/wfs?request=DescribeFeatureType&typename=tanzania:Tan_regions";     ////Send  request.     $cwfs = curl_init();     curl_setopt($cwfs, CURLOPT_FORBID_REUSE, 1);     curl_setopt($cwfs, CURLOPT_FRESH_CONNECT, 1);      curl_setopt($cwfs, CURLOPT_HEADER, 0);     curl_setopt($cwfs, CURLOPT_URL, $url);     curl_setopt($cwfs, CURLOPT_RETURNTRANSFER, 1);     curl_setopt($cwfs, CURLOPT_VERBOSE, 1);          //Collect result.     $feature_description = curl_exec($cwfs) or die(curl_error($cwfs));     curl_close ($cwfs);      What am I do wrong?


title: Tableau 8 integration with Geoserver 2.3.3

text: I am using **Tableau 8** and trying to connect to Geoserver 2.3.3, but when I make a layer group on Geoserver, Tableau 8 is unable to detect WMS service of Geoserver. I have to remove the layer group in order for Tableau to detect the WMS service, but that way I can only display layers in Tableau and no layer groups. Kindly advice on how to fix this issue.


title: How to request elevation data from ERDAS Apollo?

text: I'm developing a custom application(Desktop based using uDig SDK) where user have a provision to select the line feature on map to view its elevation chart. So far I can do this using Google elevation API (as it returns elevation values in JSON format using which I can render a chart with the help of Google chart API). Now I want to do the same task by requesting the elevation data from ERDAS Apollo (i think it supports WCS) as client has hosted elevation data i.e. DEM on his ERDAS Apollo server, but I don't find any helpful documentation on "How to request elevation data from ERDAS Apollo".


title: How to add a layer group using pyQGIS?

text: Good morning, I would like to know how I can create a new group in the layer manager using python code, someone could help me?


title: Given a flight path given in latitudes and longitudes, how do I find another flight path that is parallel to the path given?

text: Given a flight route that is given in latitudes and longitudes, how do I find another route that is "parallel" to the route given (find a perpendicular line to each waypoint given)? I think the best way to do this is to transform the geographic coordinates to Cartesian and find the parallel route and then converted it back to geographic coordinates but I never get it right. I tried the following approach also but I am not 100% sure if this is correct and I never get the last point of the "parallel" flight path:


title: ogr2ogr GeoJSON GeometryCollection

text: All, I have a Postgres/PostGIS DB with two geometries in a single table. I want to export them as a GeoJSON GeometryCollection (see here) but I can't seem to get the GEOMETRY_AS_COLLECTION layer create option to work properly in ogr2ogr.               ogr2ogr -skipfailures -progress -s_srs EPSG:4326 -t_srs EPSG:4326  -f "GeoJSON" .\myFile.json PG:"host=server user=postgres dbname=name password=password" -sql "select * from table" -nln test -nlt GeometryCollection -lco GEOMETRY_AS_COLLECTION=YES      Any help on how to get this working would be greatly appreciated! A


title: Can I Turn Off or Remove Map Layers in QGIS Print Composer

text: The title pretty much says it all. I have a base map in QGIS that I am using with several Print Composers. In one of the Print Composers, I don't want one of the layers to show up. Is there any way to accomplish this?


title: What tool breaks up raster cells into individual values based upon contiguity?

text: ![The green values have a COUNT attribute of 13418](http://i.stack.imgur.com/wnSMX.jpg) The green values have a COUNT attribute of 13418. I will want to use a tool that lumps the COUNT attribute and reclassifies based on values. But first it appears I will need to break up the raster cells into individual values using a tool? I am learning raster, in the past I only used vector.


title: How to put North Arrow and Scale on map in COMPOSER in QGIS 2.2.0

text: I am new to QGIS and the questions asked and answered here about the North Arrow do not seem to help me solve my problem below. In the main mapping pane of QGIS 2.2.0 Valmiera, I can add the scale and north arrow, but as soon as I go to the Composer to print a map, those two things disappear from the map that I open in the composer. The Composer has something which allows me to add a scale thing, but it shows the actual scale on the paper that I'll be printing on and not the actual scale of the surface of the earth of the displayed map. I even tried refreshing the map but it did not work. The bottom line is this: How do I get a scale thing (sorry!) (expanding and contracting in values as I zoom in and out) in Composer, and how do I get a North Arrow in there as well.


title: Get records that are in specific area - defined by array of points

text: I have a Rails application. I have records that contain long/lat field. I want to be able query the database with array of long/lat points (together they form closed shape area), and get all the records that are in this specified area? I have MySql database. Any gem recommendations? I have tried Ruby Geocoder and alexreisner/geocoder gem but was not able to find proper functionality. Maybe I have missed something.


title: How to display a Shapefile on a browser? -offline-

text: Plot points on a **Shapefile** like this: https://github.com/gka/pyshpgeocode and then convert it to **HTML** and display it on a browser? Is this practical method to display the map on a browser?   Any cons of it? **Is there any other way of displaying a Shapefile on a browser?**


title: Strange issue with Tomcat and AGS 10 .NET

text: I have one of my applications being deploying through Tomcat of Java WebADF. It's being done by someone else, on another site. It, as it would, worls perfectly on my box, but on their application, it fails. Initially, all that displays is a empty page, with a Javascript error highlighted:               Message: 'EsriPageElement' is undefined              Line: 17         Char: 2         Code: 0         URI: http://Myserver:8080/MyApp/js/MA_startup.js      However, on investigating the Tomcat logs, I find this: `SEVERE: Unable to initialize resource 'ags1'. This resource will be removed com.esri.adf.web.ags.ADFAGSException: Cannot create MapServer instance for server object 'MyService' [hosts=[MyServer]] at com.esri.adf.web.ags.data.AGSLocalConnection.activateResource(Unknown Source)` ...and I am a bit perplexed. The user is logging on as himself. The only thing I can think of, is that this is permissions related, in that when Tomcat is fired up, it hasn't been started by anyone with the permissions to 'see' AGS? Any ideas peeps? Permissions? IT Mice under the hood? Too much alcohol abuse in my youth? Am I bad and this is Karma? I genuinely don't have a clue.


title: Stereographic projection of WGS84 ellipsoid on a plane[python]

text: I would like to ask you for help in following task.   What I need to do is to derive local x,y coordinates from latitude, longitude, applying stereographic projection of WGS84 ellipsoid on a plane, when the following parameters are known:               scale = 0.999790760649094     x0 = 41614.2107651061     y0 = 17150.1692507701     L0 = 18.4514121640352     B0 = 54.4353877827032      There's also a Q parameter but at the moment I hava no idea what is and what for it is:                 Q = -0.0203604068846608      The above parameters are given by GPS data collecting software and become from site calibration using pairs of know [local cs] and surveyed [wgs84] points. Example input coordinates:                 latitude = 54.3324735043758     longitude = 18.6188610214872      Expected result:                 x = 30395.719     y = 28273.712      Whether is this possible only with those parameters provided or I need something else?   Can it be done using Proj4 or pyproj?   If yes, what would be apropriate command?   If not, can you hint me how to do it by hand, like links to formulas etc, please?   Thanks in advance.    [edit] Adding some background to my question. In the case if my question seems absurd...   As I mentioned above the parameters become from site calibration or GPS alignmnet using pairs of know and surveyed points. I was told that after alignment is done, what I have to do to obtain x,y coordinates for any next surveyed wgs84 point, is to apply stereographic projection of WGS84 ellipsoid on a plane where:   B0, L0 are geodetic coordinates of the tangent point   X0, Y0 are grid coordinates of the tangent point   scale [k] is the scale factor   and with the least square method parameters of the transition to the local cs are determined.      Sorry if anything here is unclear but I am not quite familiar with subject.   Any help would be appreciated.      [edit]Adding additional sample points.                    latitude longitude X Y     54.438657278 18.4602869845 41989.871 17718.469     54.4192483588 18.4806069996 39857.12 19081.45     54.4124384305 18.4716474826 39087.23 18515.701     54.3891511496 18.4799499633 36506.74 19108.41     54.3410205023 18.5443421778 31239.26 23406.7     54.3498711013 18.647572105 32374.58 30096.13     54.3515138914 18.6760425965 32600.628 31942.594     54.3800402814 18.6055277359 35669.507 27287.217     54.3971889473 18.6178833183 37596.017 28046.484     54.4041298674 18.6607429298 38432.47 30811.43     54.4070050071 18.6375394617 38717.65 29297.69     54.3609834076 18.474278312 33364.349 18804.941      [edit]Final word. Lests start from the beginning.   It turned out that my question was not quite right, in its original form.   So I decided to post full context and all steps I make in order to achieve the ultimate goal which is: getting local grid coordinates.      1.Site calibration points pairs.                 latitude longitude x y     54.4353877827 18.451412164 41614.22 17150.17     54.3625019774 18.4469405743 33496.96 17024.81     54.3097484306 18.5670103377 27791.39 24957.15     54.3084348677 18.6299997771 27737.21 29059.27     54.402430173 18.6683595667 38254.91 31310.35     54.4230519006 18.6017101505 40450.41 26931.58     54.438657278 18.4602869845 41989.871 17718.469     54.4192483588 18.4806069996 39857.12 19081.45     54.4124384305 18.4716474826 39087.23 18515.701     54.3891511496 18.4799499633 36506.74 19108.41     54.3410205023 18.5443421778 31239.26 23406.7     54.3498711013 18.647572105 32374.58 30096.13     54.3515138914 18.6760425965 32600.628 31942.594     54.3800402814 18.6055277359 35669.507 27287.217     54.3971889473 18.6178833183 37596.017 28046.484     54.4041298674 18.6607429298 38432.47 30811.43     54.4070050071 18.6375394617 38717.65 29297.69     54.351702123 18.6831419293 32632.47 32403.62     54.4427340176 18.4391134326 42415.52 16335.76     54.4427368643 18.4592502382 42442.44 17641.9     54.3609834076 18.474278312 33364.349 18804.941     54.3714455799 18.7279691735 34899.27 35264.2     54.2911268152 18.6361668642 25820.32 29504.89      2.Convert above surveyed latitude and longitude to stereographic plane coordinates, taking coordinates of first calibration point as central latitude and longitude.   I don't know how important the scale factor [k_0] is so decided to force 1.0.   I also experimented with UTM instead of stereographic plane but result was not satisfactory with residuals = ~0.1m because we are talking here about obtaining survey grade coordinates not gis accuracy.               import pyproj          #input coordinates         p1 = pyproj.Proj(proj="latlong",datum="WGS84", ellps="WGS84")          #output result     p2 = pyproj.Proj(proj="stere", lat_0=54.4353877827032, lon_0=18.4514121640352, k_0=1.0, x_0=0, y_0=0, ellps="WGS84")          stereographic_coords= pyproj.transform(p2, p1, latitude, longitude)     print stereographic_coords[1], stereographic_coords[0] #keep N(=x), E(=y) order      The result:                 N E     3.50449394074e-07 2.28385549832e-06     -8113.07087572 -290.673432675     -13978.9604246 7524.00783763     -14116.6372391 11624.2036175     -3646.95639951 14088.88636     -1362.75920089 9755.69391324     363.973973779 575.837791136     -1796.13420238 1895.18404979     -2554.36871405 1313.79387579     -5146.3420395 1853.8865397     -10500.2307161 6044.01338437     -9501.32640299 12755.170692     -9312.95200167 14605.8597464     -6149.937535 10013.9390815     -4239.26761866 10812.2644791     -3459.25317782 13593.690851     -3143.43494621 12086.0443555     -9290.50857805 15067.4007146     817.802325033 -797.917280294     818.077929071 508.518662863     -8281.8734697 1486.45713731     -7082.36334017 17973.5167497     -16042.1577585 12030.6541087      3.Similarity transformation.   We can perform now similarity transformation to obtain local coordinates for any given point. One thing to remember is we firstly need to convert given latitude, longitude to stereographic coordinates. transformation.py               import math     import numpy as np               class GeometricTransform(object):     """Perform geometric transformations on a set of coordinates.          """         def __call__(self, coords):            raise NotImplementedError()              def inverse(self, coords):             raise NotImplementedError()              def __add__(self, other):             raise NotImplementedError()               class ProjectiveTransform(GeometricTransform):         _coeffs = range(8)              def __init__(self, matrix=None):             if matrix is None:             # default to an identity transform                 matrix = np.eye(3)             if matrix.shape != (3, 3):                 raise ValueError("invalid shape of transformation matrix")             self._matrix = matrix              @property         def _inv_matrix(self):             return np.linalg.inv(self._matrix)              def _apply_mat(self, coords, matrix):             coords = np.array(coords, copy=False, ndmin=2)                  x, y = np.transpose(coords)             src = np.vstack((x, y, np.ones_like(x)))             dst = np.dot(src.transpose(), matrix.transpose())              # rescale to homogeneous coordinates             dst[:, 0] /= dst[:, 2]             dst[:, 1] /= dst[:, 2]                  return dst[:, :2]              def __call__(self, coords):             return self._apply_mat(coords, self._matrix)              def inverse(self, coords):             return self._apply_mat(coords, self._inv_matrix)              def estimate(self, src, dst):             xs = src[:, 0]             ys = src[:, 1]             xd = dst[:, 0]             yd = dst[:, 1]             rows = src.shape[0]              # params: a0, a1, a2, b0, b1, b2, c0, c1             A = np.zeros((rows * 2, 9))             A[:rows, 0] = xs             A[:rows, 1] = ys             A[:rows, 2] = 1             A[:rows, 6] = - xd * xs             A[:rows, 7] = - xd * ys             A[rows:, 3] = xs             A[rows:, 4] = ys             A[rows:, 5] = 1             A[rows:, 6] = - yd * xs             A[rows:, 7] = - yd * ys             A[:rows, 8] = xd             A[rows:, 8] = yd              # Select relevant columns, depending on params             A = A[:, self._coeffs + [8]]                  _, _, V = np.linalg.svd(A)                  H = np.zeros((3, 3))         # solution is right singular vector that corresponds to smallest         # singular value             H.flat[self._coeffs + [8]] = - V[-1, :-1] / V[-1, -1]             H[2, 2] = 1                  self._matrix = H              def __add__(self, other):         """Combine this transformation with another.              """             if isinstance(other, ProjectiveTransform):             # combination of the same types result in a transformation of this             # type again, otherwise use general projective transformation                 if type(self) == type(other):                     tform = self.__class__                 else:                     tform = ProjectiveTransform                 return tform(other._matrix.dot(self._matrix))             else:                  raise TypeError("Cannot combine transformations of differing "                             "types.")                              class SimilarityTransform(ProjectiveTransform):     """2D similarity transformation of the form::                  X = a0*x - b0*y + a1 =               = m*x*cos(rotation) + m*y*sin(rotation) + a1                  Y = b0*x + a0*y + b1 =               = m*x*sin(rotation) + m*y*cos(rotation) + b1          where ``m`` is a zoom factor and the homogeneous transformation matrix is::              [[a0  b0  a1]          [b0  a0  b1]          [0   0    1]]          Parameters     ----------     matrix : (3, 3) array, optional         Homogeneous transformation matrix.     scale : float, optional         Scale factor.     rotation : float, optional         Rotation angle in counter-clockwise direction as radians.     translation : (tx, ty) as array, list or tuple, optional         x, y translation parameters.          """              def __init__(self, matrix=None, scale=None, rotation=None,                  translation=None):             params = any(param is not None                      for param in (scale, rotation, translation))                  if params and matrix is not None:                 raise ValueError("You cannot specify the transformation matrix and "                              "the implicit parameters at the same time.")             elif matrix is not None:                 if matrix.shape != (3, 3):                     raise ValueError("Invalid shape of transformation matrix.")                 self._matrix = matrix             elif params:                 if scale is None:                     scale = 1                 if rotation is None:                     rotation = 0                 if translation is None:                     translation = (0, 0)                      self._matrix = np.array([                     [math.cos(rotation), - math.sin(rotation), 0],                     [math.sin(rotation),   math.cos(rotation), 0],                     [                 0,                    0, 1]                 ])                 self._matrix *= scale                 self._matrix[0:2, 2] = translation             else:             # default to an identity transform                 self._matrix = np.eye(3)              def estimate(self, src, dst):             xs = src[:, 0]             ys = src[:, 1]             xd = dst[:, 0]             yd = dst[:, 1]             rows = src.shape[0]              # params: a0, a1, b0, b1             A = np.zeros((rows * 2, 5))             A[:rows, 0] = xs             A[:rows, 2] = - ys             A[:rows, 1] = 1             A[rows:, 2] = xs             A[rows:, 0] = ys             A[rows:, 3] = 1             A[:rows, 4] = xd             A[rows:, 4] = yd                  _, _, V = np.linalg.svd(A)              # solution is right singular vector that corresponds to smallest         # singular value             a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1]                  self._matrix = np.array([[a0, -b0, a1],                                  [b0,  a0, b1],                                  [ 0,   0,  1]])              @property         def scale(self):             if math.cos(self.rotation) == 0:             # sin(self.rotation) == 1                 scale = self._matrix[0, 1]             else:                  scale = self._matrix[0, 0] / math.cos(self.rotation)             return scale              @property         def rotation(self):                  return math.atan2(self._matrix[1, 0], self._matrix[1, 1])              @property         def translation(self):             return self._matrix[0:2, 2]                         TRANSFORMS = {         'similarity': SimilarityTransform,     }     HOMOGRAPHY_TRANSFORMS = (         SimilarityTransform,         ProjectiveTransform     )                def estimate_transform(ttype, src, dst, **kwargs):         ttype = ttype.lower()         if ttype not in TRANSFORMS:             raise ValueError('the transformation type \'%s\' is not'                          'implemented' % ttype)              tform = TRANSFORMS[ttype]()         tform.estimate(src, dst, **kwargs)              return tform               def matrix_transform(coords, matrix):         return ProjectiveTransform(matrix)(coords)      And here's my calculation:                 import numpy as np     import transformations as tf #importing above script          ##############################     ####CALIBRATION POINTS########     ##############################          #primary system coordinates obtained instep 2 from surveyed latitude,longitude.          ps1=( 3.50449394074e-07 , 2.28385549832e-06 )     ps2=( -8113.07087572 , -290.673432675 )     ps3=( -13978.9604246 , 7524.00783763 )     ps4=( -14116.6372391 , 11624.2036175 )     ps5=( -3646.95639951 , 14088.88636 )     ps6=( -1362.75920089 , 9755.69391324 )     ps7=( 363.973973779 , 575.837791136 )     ps8=( -1796.13420238 , 1895.18404979 )     ps9=( -2554.36871405 , 1313.79387579 )     ps10=( -5146.3420395 , 1853.8865397 )     ps11=( -10500.2307161 , 6044.01338437 )     ps12=( -9501.32640299 , 12755.170692 )     ps13=( -9312.95200167 , 14605.8597464 )     ps14=( -6149.937535 , 10013.9390815 )     ps15=( -4239.26761866 , 10812.2644791 )     ps16=( -3459.25317782 , 13593.690851 )     ps17=( -3143.43494621 , 12086.0443555 )     ps18=( -9290.50857805 , 15067.4007146 )     ps19=( 817.802325033 , -797.917280294 )     ps20=( 818.077929071 , 508.518662863 )     ps21=( -8281.8734697 , 1486.45713731 )     ps22=( -7082.36334017 , 17973.5167497 )     ps23=( -16042.1577585 , 12030.6541087 )          #secondary system coordinates/known local coordinates     ss1 = ( 41614.22 , 17150.17 )     ss2 = ( 33496.96 , 17024.81 )     ss3 = ( 27791.39 , 24957.15 )     ss4 = ( 27737.21 , 29059.27 )     ss5 = ( 38254.91 , 31310.35 )     ss6 = ( 40450.41 , 26931.58 )     ss7 = ( 41989.871 , 17718.469 )     ss8 = ( 39857.12 , 19081.45 )     ss9 = ( 39087.23 , 18515.701 )     ss10 = ( 36506.74 , 19108.41 )     ss11 = ( 31239.26 , 23406.7 )     ss12 = ( 32374.58 , 30096.13 )     ss13 = ( 32600.628 , 31942.594 )     ss14 = ( 35669.507 , 27287.217 )     ss15 = ( 37596.017 , 28046.484 )     ss16 = ( 38432.47 , 30811.43 )     ss17 = ( 38717.65 , 29297.69 )     ss18 = ( 32632.47 , 32403.62 )     ss19 = ( 42415.52 , 16335.76 )     ss20 = ( 42442.44 , 17641.9 )     ss21 = ( 33364.349 , 18804.941 )     ss22 = ( 34899.27 , 35264.2 )     ss23 = ( 25820.32 , 29504.89 )          ############################     ############################     ############################          from_pt = [ps1, ps2, ps3, ps4, ps5, ps6, ps7, ps8, ps9, ps10, ps11, ps12, ps13, ps14,  ps15, ps16, ps17, ps18, ps19, ps20, ps21, ps22, ps23]     to_pt =   [ss1, ss2, ss3, ss4, ss5, ss6, ss7, ss8, ss9, ss10, ss11, ss12, ss13, ss14, ss15, ss16, ss17, ss18, ss19, ss20, ss21, ss22, ss23]          fr = np.array(from_pt)     to = np.array(to_pt)     tform = tf.estimate_transform('similarity', fr, to)          #list of coordinates to transform, first converted to sterographic     l = [     [10892.8449025, -11442.6688292],     [10892.2806078, -11442.7145534],     [10880.2426997, -11440.0353547],     [10869.0233644, -11437.2787524],     [10868.4894383, -11437.3154292],     [10879.4637253, -11439.8965462],     [10879.194897, -11440.0837157],     [11661.9477535, -12607.9800965],     [11620.0600416, -12597.7915158],     [11600.121001, -13934.58708]]          for line in l:         x = line[1]         y = line[0]         res= tform(np.array([x,y]))         print res      Final result:                 my calculation result                   known control points coordinates     [[ 30395.74138597  28273.70706175]]     # 30395.716 28273.716     [[ 30395.68418212  28273.14381546]]     # 30395.659 28273.153     [[ 30398.1177287   28261.05386387]]     # 30398.092 28261.063     [[ 30400.64532909  28249.78073875]]     # 30400.620 28249.790     [[ 30400.59778907  28249.24767053]]     # 30400.572 28249.257     [[ 30398.24064822  28260.27222544]]     # 30398.215 28260.281     [[ 30398.04804427  28260.0072639 ]]     # 30398.022 28260.016     [[ 29246.33177298  29066.37582526]]     # 29246.308 29066.384     [[ 29255.66538993  29024.28939019]]     # 29255.641 29024.297     [[ 27918.74210572  29031.57195305]]     # 27918.718 29031.581      I think, residuals are acceptable.


title: OpenLayers - SelectFeature fails on a small polygon on top of a larger polygon

text: I'm using `new OpenLayers.Control.DrawFeature` to draw polygons and `OpenLayers.Control.SelectFeature` to select a specific polygon. I'm drawing one polgon and then a smaller polygon on top of the first one. PolyA is a large polygon(50 x 50) PolyB is a smaller polgon (10 x 10) drawn within (on top of) PolyA Initially I can select and unselect PolyB correctly. However once I have selected and unselected PolyA, PolyB is now hidden behind it and can no longer be selected. Why is that and what can I do to resolve it? Thanks. select code (straightforward)               select  = new OpenLayers.Control.SelectFeature(layer,                     {                         highlightOnly: true,                             clickout: true, toggle: true,                         onSelect: function(feature) {                              //do stuff               }                                     });


title: Grid vs TIF vs IMG

text: As the title suggests both of the above three formats are raster files. In ArcGIS, apparently there's no difference between them. So what's the difference between these three formats?


title: Problems getting openlayers permalink working

text: I set up an openlayers permalink function. I set the code and it shows up. But the URL does not work. Where is I get a link like this: http://babkart.no/?zoom=9&lat=7016094.01622&lon=121216.28057&layers=BFFFTTTTTF The only thing that works is the chosen layers. It does not show the correct zoom or boundingbox. What do I need in the code to get this to work? I can not see I have any variables like "lat", "lon" or "zoom". Any examples as how to get this to work? I also would like to put the permalink in the toolbar so that one can copy the URL. I am using ExtJS 2.2. But this is another issue and I have to get the permalink to work first. I add the permalink like this:               controls: [new OpenLayers.Control.KeyboardDefaults(),new OpenLayers.Control.PanZoomBar(),new OpenLayers.Control.ScaleLine(),                     new OpenLayers.Control.Navigation(),new OpenLayers.Control.Permalink()]      My app.js is here if anyone wants to take a look: http://pastebin.com/AstbpBWV The site is running here: http://babkart.no/ ~asle I got this to work this way. Don't know if it was good or bad but it seems to work: I use a function to zoom in to a location and added functionality to just zoom in to the view. I use the OpenLayers.Permalink() function.                       // If URL contains ref. to a location:             var args = OpenLayers.Util.getParameters();             if (args.x && args.y) {                 var position = new OpenLayers.LonLat(parseFloat(args.x), parseFloat(args.y));                 map.setCenter(position, 14);                 var xy = map.getViewPortPxFromLonLat(position);                 xy.x = Math.round(xy.x);                 xy.y = Math.round(xy.y);                 info.getInfoForClick({xy: xy});             } // read the permalink URL             else if (args.lat && args.lon && args.zoom) {                 var position = new OpenLayers.LonLat(parseFloat(args.lon), parseFloat(args.lat));                 map.setCenter(position, parseFloat(args.zoom));             }             else {                     map.setCenter(new OpenLayers.LonLat(528277, 7188837), 3);             }


title: Can't add map service (map cache) to mxd

text: When trying to add a map service (in this case a map cache service) to an ArcMap project, I keep getting an error that says: Could not add specified data object to the map. Could not connect to the Map Server. Even though in Catalog I can connect to the service perfectly fine and open it as an mxd project, I can't add it to another project. The goal is to use the map cache as a basemap for editing natural features from an aerial photo. Any help would be greatly appreciated. Using ArcMap 10, ArcServer 10 Thanks!


title: Tutorials for UML/GML

text: Do you know any good tutorials or any resources for learning from the beginning UML and GML ? I am asking because (later) i would like to learn how to make a application schemas in context of INSPIRE. Thank's a lot!


title: Update or destroy popup in openlayers

text: I am doing a map which refresh dinamically all markers when they get modified. I want to add popups to show some info. I create popups this way:                       var feature = new OpenLayers.Feature(markers, location); //location = lonlat             feature.closeBox = true;             feature.popupClass = OpenLayers.Class(OpenLayers.Popup.FramedCloud, {             'autoSize': true         });;             feature.data.popupContentHTML = "info";              feature.data.overflow = (false) ? "auto" : "hidden";                         var markerClick = function (evt) {                 if (this.popup == null) {                     this.popup = this.createPopup(this.closeBox);                     map.addPopup(this.popup);                     this.popup.show();                 } else {                     this.popup.toggle();                 }                 currentPopup = this.popup;                 OpenLayers.Event.stop(evt);             };             mark.events.register("mousedown", feature, markerClick);      Now, I have the problem about what to do if one of my markers update its info. How can I destroy the popup? (I will create other one) or how can I modifiy it? Can I access them using lonlat? Sometimes, I want to delete all of them, is there a way to do that? like deleting a layer or something?


title: Adding an external WMS service using REST service with Geoserver

text: I'm trying to add an External WMS layer to geoserver with a php page, using the REST service of Geoserver. I've already read the tutorials for creating the geoserver features (http://docs.geoserver.org/stable/en/user/restconfig/rest-config-api.html), but i couldn't find the way to create a new store and layer with an external wms service. in other words: i have the url of the external wms service (=store), the layer i wanto to publish (=layer) and the epgs i want to use. Is there a way to load all of these information in the geoserver schema? Thank you.


title: Dynamically change zoom level of overview map to fit with current extent in DDP

text: I am using DDP pages in a large mapset with pages ranging from 1:7500 to 1:130,000 and the overview map is zoomed out to far to see some of the smaller zones. Is it possible to set the overview map to zoom to say 10 times the current DDP page extent? so a 50k main sheet would have the overview at 500k. Map overview scale @ 1:3,000,000 - meaningless as it is to full of info. ![enter image description here](http://i.stack.imgur.com/GjkEC.jpg) Map overview scale @ 1:500,000 ![enter image description here](http://i.stack.imgur.com/ExU0j.jpg) So we can use the cartography-->ddp-->calculate adjacent fields to get the adjoining map sheets but this doesn't seem to help in changing the zoom levels of the map. Till now my workaround has been to create multiple mxd's which have all the maps in an overview zone that makes sense but this causes issues when making edits to pages etc etc and it would be ideal to automate it into one mxd. Also is it possible to get DDP to fill in the total number of maps -this should be map 14 of 38 but the 13 is carried over from a previous version as it's hard coded. \---- ONE SOLUTION I also found another workaround at 1am on a Monday morning... in the overview map data frame properties, go to Data Frame tab and choose options as per the screen capture below. ![enter image description here](http://i.stack.imgur.com/xb4On.jpg) This seems to work fine as long as the extents don't change to much. I am keeping this question "unanswered" to see what other options are available.


title: ArcGIS Flex Polygon to PostGIS Polygon

text: ArcGIS supports building polygons within the flex component that contain self- intersecting rings, as well as rings that touch at multiple points (breaking the "arbitrary" rules here: http://workshops.opengeo.org/postgis- intro/validity.html ). This causes incompatibility when passing these shapes to PostGIS. The article above mentions buffering to fix some of these issues, but this seems to only fix some cases, and often ends up with a different shape polygon than was intended (eg, only half of the poly). Is there a way to convert these to monotonic polygons within flex so that this problem is eliminated? Or a Java API to do the same? Or is there some other approach entirely to solving this problem?


title: Is there any way to fast read pixel value using gdal in C#

text: I create with C# a series of lat/lon value pair in the code(using i, j loop) and I want to extract the pixel value based on either lat/lon pair or pixel location (i, j). I know there is one gdal utility called gdallocationinfo which can read pixel value very fast. But it is written in python. I am wondering whether there is a way to use this gdal utility in my C# code, or other methods that can extract pixel value in a faster manor. Thank you in advance! **Update with final solution:**               static void Main(string[] args)         {             Console.WriteLine(DateTime.Now.ToString());             Gdal.AllRegister(); //Register all gdal data drives             List<int> resolutionId = new List<int>() { 9 };                  Dataset ds = Gdal.Open("D:/SKnoZero.tif", Access.GA_ReadOnly);   //Read raster             double[] gt = new double[6];             ds.GetGeoTransform(gt); //Read geo transform info into array             int Rows = ds.RasterYSize;  //Get the number of rows              int Cols = ds.RasterXSize;  //Get the number of columns             Band band = ds.GetRasterBand(1);    //Read band 1             double startX = gt[0];  //Upper left lon             double startY = gt[3];  //Upper left lat             double interval = gt[1];    //Cell size             double x, y;    //Current lon and lat                  var filePath = "D:/SKcsv.csv";             using (StreamWriter csv = new StreamWriter(filePath))             {                 for (int k = 0; k < Rows; k++)  //read one line at a time                 {                     y = startY - k * interval;  //current lat                     int[] buf = new int[Cols];                     //ReadRaster parameters are StartCol, StartRow, ColumnsToRead, RowsToRead, BufferToStoreInto, BufferColumns, BufferRows, 0, 0                     band.ReadRaster(0, k, Cols, 1, buf,Cols,1,0,0);                     //iterate each item in one line                     for (int r = 0; r < Cols; r++)                       {                         if (buf[r] != -32768)   //if pixel value is not NoData value                         {                             x = startX + r * interval;  //current lon                                                          var geoidpair = GetVrgId(x, y, resolutionId);   //get geoid                             var temp = geoidpair;                             var newLine = string.Format("{0},{1},{2},{3}{4}", x, y, buf[r], temp[0], Environment.NewLine);                             csv.Write(newLine);                                     }                      }              }              Console.WriteLine(DateTime.Now.ToString());              Console.WriteLine("Press any key");              Console.ReadLine();         }


title: GCS_WGS_1984 and 900913. What is the difference?

text: I use WMS Service with mymapproxy/mapproxy.yaml :               services:       demo:       kml:       tms:       wmts:       wms:         srs: ['EPSG:4326', 'EPSG:900913']         md:           title: MapProxy OSM WMS Proxy           abstract: This is the fantastic MapProxy.     layers:       - name: osm         title: Local OSM WMS         sources: [osm_cache]          caches:       osm_cache:         grids: [global_mercator_osm, global_geodetic_sqrt2]         sources: [osm_mapnik2]          sources:       osm_mapnik:         type: mapnik         mapfile: /home/oleg/src/mapnik-style/osm.xml       osm_mapnik2:         type: tile         grid: global_mercator_osm         url: http://localhost/osm_tiles/%(z)s/%(x)s/%(y)s.png          grids:       global_mercator_osm:         base: GLOBAL_MERCATOR         srs: 'EPSG:900913'         origin: nw       global_geodetic_sqrt2:         base: GLOBAL_GEODETIC         res_factor: 'sqrt2'         srs: 'EPSG:900913'         origin: nw      Then I display this WMS service in Geoserver. I **have not any problems** (if i set **Native и Declared SRS** as **900913** ) But I'd like display this service in ArcMap 10.1 This is the coordinate system of my ArcGis service. ![enter image description here](http://i.stack.imgur.com/4qGMB.jpg) So, I don't see the 900913 coordinate system. And when I add this service to TOC (left panel with layers in ArcMap) I have a lot of errrors: ![enter image description here](http://i.stack.imgur.com/30GKa.jpg) ![enter image description here](http://i.stack.imgur.com/Rq67f.jpg) Sometimes I do not see a map or this map is not good: ![enter image description here](http://i.stack.imgur.com/wqZss.jpg)


title: Can 3rd party javascript libraries be used with CartoDB?

text: I'd like to build some functionality (search, sliders, etc.) in a 3rd party javascript library (jquery mobile), and integrate it with CartoDB. Can this be done? Are there any examples available or a cookbook decribing how to do the integration ?


title: Fill no data area over ocean in orthophoto mosaic

text: We have recently received orthoimagery tiles for my county. The county is a series of islands and the imagery only covers the islands. Thus we have data gaps where there is only water. In past projects we've created a mosaic by filling in the gaps using a solid blue color that somewhat resembles the median color of the water in the tiles near the shorelines. It worked pretty well in the past, but with this project I'm finding that no single color is working best. Is there a method I can use to make the transition from imagery edge to a solid color more aesthetically pleasing?


title: How to create a symbol with graduated marker size and line width in pyQGIS?

text: I am using QGIS + pyQGIS + Python. Layer geometry is line. Each line has a magnitude and a direction (from 1st node to 2nd node). What exactly I want to perform is **"graduated width marker symbology overlaid on the graduated width line symbology"** Example: line thickness 5mm. and 5mm. width marker symbol overlaid on it. **===(5mm)===** line thickness 3mm. and 3mm. width marker symbol overlaid on it. **\-----** (3mm) **\-----** For the time being, I can perform the following; either "(i) lines with graduated thickness" or "(ii) single line symbology overlaid with single marker symbol (no graduated thickness)". Can I perform it at the same time? **Graduated width markers on graduated thickness lines.** **Below code does all same appearing markers overlaid on all same appearing lines:**               #Use the currently selected layer     layer = qgis.utils.iface.mapCanvas().currentLayer()     registry = QgsSymbolLayerV2Registry.instance()     lineMeta = registry.symbolLayerMetadata("SimpleLine")     markerMeta = registry.symbolLayerMetadata("MarkerLine")          symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())          #Line layer     lineLayer = lineMeta.createSymbolLayer({'width': '0.26', 'color': '255,0,0', 'offset': '0', 'penstyle': 'solid', 'use_custom_dash': '0', 'joinstyle': 'bevel', 'capstyle': 'square'})          #Marker layer     markerLayer = markerMeta.createSymbolLayer({'width': '0.26', 'color': '255,0,0', 'rotate': '1', 'placement': 'centralpoint', 'offset': '0'})     subSymbol = markerLayer.subSymbol()          #Replace the default layer with our own SimpleMarker     subSymbol.deleteSymbolLayer(0)     triangle = registry.symbolLayerMetadata("SimpleMarker").createSymbolLayer({'name': 'filled_arrowhead', 'color': '255,0,0', 'color_border': '0,0,0', 'offset': '0,0', 'size': '3', 'angle': '0'})     subSymbol.appendSymbolLayer(triangle)          #Replace the default layer with our two custom layers     symbol.deleteSymbolLayer(0)     symbol.appendSymbolLayer(lineLayer)     symbol.appendSymbolLayer(markerLayer)          #Replace the renderer of the current layer     renderer = QgsSingleSymbolRendererV2(symbol)     layer.setRendererV2(renderer)              QgsMapLayerRegistry.instance().addMapLayer( layer )      **Below code performs only graduated width lines (equal interval) BUT NO MARKERS OVERLAID**               def validatedDefaultSymbol(geometryType):         symbol = QgsSymbolV2.defaultSymbol(geometryType)         if symbol is None:             if geometryType == QGis.Point:                 symbol = QgsMarkerSymbolV2()             elif geometryType == QGis.Line:                 symbol =  QgsLineSymbolV2()             elif geometryType == QGis.Polygon:                 symbol = QgsFillSymbolV2()         return symbol          def makeSymbologyForRange(layer, min ,max, label ,colour, alpha, width):         symbol = validatedDefaultSymbol(layer.geometryType())         symbol.setColor(colour)         symbol.setAlpha(alpha)         symbol.setWidth(width)         range = QgsRendererRangeV2(min, max, symbol, label)         return range                              vlayer = QgsVectorLayer(str(SaveShpName), str(SaveShpName), 'ogr')     myTargetField = 'magnitude'          #For Equal Interval Classes     if dlg.ui.comboBoxSelectSymbology.currentText() == "Equal Interval":         GradSymNoOfClasses = dlg.ui.spinBoxClasses.value()          GradSymInterval = round(((int(GradSymMax) - int(GradSymMin)) / float(GradSymNoOfClasses)),0)         myRangeList = []         for i in range(GradSymNoOfClasses):             if i == 0:                 classlabel = str(GradSymMin)+" - "+str(int(GradSymMin)+GradSymInterval)                 myRangeList.append(makeSymbologyForRange( vlayer, int(GradSymMin), int(GradSymMin)+GradSymInterval, classlabel , QColor(0, 255-(255*i/GradSymNoOfClasses), 255*i/GradSymNoOfClasses), 1, 0.25 ))             elif i == (GradSymNoOfClasses - 1):                 classlabel = str(int(GradSymMin)+(GradSymInterval*i)+0.001)+" - "+str(GradSymMax)                 myRangeList.append(makeSymbologyForRange( vlayer, (GradSymInterval*i)+0.001, int(GradSymMax), classlabel , QColor(0, 255-(255*i/GradSymNoOfClasses), 255*i/GradSymNoOfClasses), 1, i/1.25 ))             else:                 classlabel = str(int(GradSymMin)+(GradSymInterval*i)+0.001)+" - "+str(int(GradSymMin)+GradSymInterval*(i+1))                 myRangeList.append(makeSymbologyForRange( vlayer, int(GradSymMin)+(GradSymInterval*i)+0.001, int(GradSymMin)+GradSymInterval*(i+1), classlabel , QColor(0, 255-(255*i/GradSymNoOfClasses), 255*i/GradSymNoOfClasses), 1, i/1.25 ))              myRenderer = QgsGraduatedSymbolRendererV2( myTargetField, myRangeList )         vlayer.setRendererV2( myRenderer )              QgsMapLayerRegistry.instance().addMapLayer( vlayer )


title: How to enable and disable IcommandItem of ArcObjects?

text: I have Commandbutton1 and Commandbutton2 in commandbar. Commandbutton1 is disabled on some condition. I want to enable Commandbutton1 on click of Commandbutton12. How do I do this?


title: Points in polygon stops processing at 0% Any ideas how to fix it?

text: I downloaded QGIS yesterday. Today, I tried to do a "Points in polygon" analysis, and it stopped working. It was working earlier today when I did it a couple times. Now, when go to do one: I click "OK", it pops up a message about the layers not having matching CRS; I click 'OK', then nothing happens. The "OK" button in the "Points in polygon" window grays out, the process bar shows 0%, the layer files are created without data (each is 1kb). I can close/cancel. If I click in QGIS outside the "Points in polygon" window, the program crashes. I've tried restarting. I've tried uninstalling and reinstalling. Any ideas?


title: OpenStreetMap XAPI request for "Member Of"

text: I'm trying to find a way to request public transportation info from the OSM XAPI. Those nodes can be best described as:               highway = *          member of = ( route = bus )      Is that even possible to do over the XAPI or do i need to query it from the planet data? this is what i tried but with no success               http://open.mapquestapi.com/xapi/api/0.6/relation[route=bus][bbox=100.63584242503583,13.659388026297398,100.65300856273045,13.66877069028344]


title: How to make the Overlays section in an OpenLayers LayerSwitcher exclusive?

text: At the simplest level I need to figure out how to create a customized OpenLayers LayerSwitcher **Overlays** section which is not using Checkboxes, but uses RadioButtons instead. For my application its not desirable for the user to be able to select multiple layers. Eventually I'm going to need to be able to create group sections in the Overlay section of the LayerSwitcher, so if anyone has any hints on that I'd appreciate it. I've seen the customized LayerSwitcher example on the OpenLayers site, but I can only get it to work in Chrome. I need to be able to support Safari on an iPad. Thanks for any advice!


title: OpenLayers.Layer.clone() not returning a unique Layer object

text: I'm taking a layer (`layer_orig`) and calling the `clone` method to create a new layer (`layer_clone`) with same features, styles, etc. I am changing the name and `map` instance of the layer with:               layer_clone.map = null;     layer_clone.setName(name);      This works exactly as expected, even giving a unique layer id to the new layer. However, if I then update the styleMap of `layer_clone`, the styleMap of `layer_orig` is updated as well - suggesting that `layer_clone` is a reference to `layer_orig` and not a completely unique Layer object. Is there a way to get a unique Layer object that is a direct clone of another layer?


title: Extract coordinates of polygon vertices in ArcMap 10?

text: I've got about a dozen polygons in a feature class loaded in ArcMap 10, all in geographic WGS 1984. How do I easily obtain the coordinates associated with each vertex of each polygon in that feature class? Ideally I would like them to be nicely tabulated in spreadsheet format.


title: How to remove small areas without reducing borders

text: I'm having a problem where I am overlaying two maps to create a new kind of map with the data conjoined. In general this is good, but I need to remove the smaller areas from it. For the interior polygons, you can see that this is fine. However, for the ones on the exterior, they are removed and area is actually lost. I don't want this to happen - I need to persist the original borders. I would like extended shape to have the data values of one of the polygons it was previously connected to. Original map: ![Original field in GRASS](http://i.stack.imgur.com/DggZu.png) Cleaned map: ![Cleaned field in GRASS](http://i.stack.imgur.com/68Xe6.png) **Additional side note...** There's an issue with `v.buffer` in Grass 6.4.3 where it will get in some sort of an infinite loop, use 100% of a CPU core, and never complete. I would like to avoid doing any operations that involve buffering the areas of the original map. However, if that's the only way please let me know. I will try to work around this problem. Thanks!


title: Multiparts and Multigeometries: checking for multigeometries with PyQGIS?

text: I have been working with **PyQGIS** and writing a small script to **check** whether the layers I have in the MapCanvas have **geometries of type "MultiXX"** because I want to run the "Multipart to Singleparts" processing algorithm on the layers that need to be processed (but not on the others because it would be too long). I need simple geometries because I am using the SPIT plugin to load them to a PostGIS database in a further step. If I used this:               for feature in layer.getFeatures():             geom = feature.geometry()             if geom.isMultipart():               print (layer.name(), ": Warning: Multipart", feature.id())     ...      I get whether every single feature is multipart or not. But I would like to determine the type of geometry on the layer itself, because I do not want to call the algorithm on the layer every time there is one feature multipart. So I tried this:               for layer in allLayers:        if layer.wkbType() > '4':           print "Warning Multigeometry in "+layer.name()         ...      But somehow my layer, in which the multipart geometries would pop up with the first test, gets a wkbType of 3 (meaning Polygon). So I am a bit confused: especially because QGIS tells me, the layer is a POLYGON layer, but the SPIT plugin tells me it is a MULTIPOLYGON layer (and indeed, there is a feature which is multipart and running the above mentionned processing algorithm resolve the problem). Could someone explain to me how the relation between the geometry type and the multiparts aspect is working (I thought it was pretty simple: polygon => singleparts and multipolygons => multiparts) or suggest to me an alternative way to reach the same result, meaning: checking accurately for layers with multigeometries in order to process them selectively.


title: Summing column values in table using arcgisscripting?

text: I need to create a python script that sums a column's values in a table. Thereafter, the script is supposed to get the percentage of each individual row in that column and populate this percentage values in a another column. I been struggling for quite some time now, I checked the net, but there is nothing of this sort. The nearest is an ArcGIS Discussion Forum thread, but its for an older version of ArcMap, or I am using it wrong.               import arcgisscripting, sys     gp = arcgisscripting.create()          intable = sys.argv[1]     field = sys.argv[2]          # Create search cursor     rows = gp.SearchCursor(intable)     row = rows.Next()     x = 0.0     # Enter while loop for each feature/row     while row:         x += row.getvalue(field)         print x         row = rows.next()     #note value can be rounded depending on field type     gp.calculatefield(intable,sys.argv[3],float(x),"PYTHON")