title: ArcGIS Cost Path Analysis - How Is Arcpy Cost Path Analysis Tool So Much Faster Than ArcObjects

text: Although I love slinging some python to create geoprocessing scripts/services, I was under the impression that using ArcObjects to do the equivalent operation(s) will have better performance. I've made one (or two) posts over the past couple of days regarding my woes with getting geoprocessing scripts that use Spatial Analyst tools to work as geoprocessing services. My deadline is fast approaching, so I have decided to go the SOE route to achieve the desired functionality. Getting a cost path analysis in ArcObjects was relatively straight-forward using the .NET ESRI.ArcGIS.SpatialAnalyst.RasterDistanceOpClass, specifically the CostDistanceFull() and CostPath() Methods. UPDATE: Some code snippets of how I am doing things: > Python               # Get Cost Path Origin and Destination Points     inputPointsShp = 'D:/RasterStuff/test_points.shp'     arcpy.MakeFeatureLayer_management(inputPointsShp,"origin",' "TYPE" = \'ORIGIN\' ')     arcpy.MakeFeatureLayer_management(inputPointsShp,"destination",' "TYPE" = \'DESTINATION\' ')          # Check out the ArcGIS Spatial Analyst extension license     arcpy.CheckOutExtension("Spatial")          # Execute CostDistance     outCostDistance = CostDistance("origin",SOURCE_RASTER,"#","backlink")          # Execute CostPath     outCostPath = CostPath("destination", outCostDistance,"backlink")          # Convert Result to Polyline     arcpy.RasterToPolyline_conversion(outCostPath, "leastCostPath")     featSet = arcpy.FeatureSet("leastCostPath")      > C#               IDistanceOp distanceOp = new RasterDistanceOpClass();     IRasterBandCollection costDistanceRaster = (IRasterBandCollection)distanceOp.CostDistanceFull((IGeoDataset)sourceFc, (IGeoDataset)raster, true, true, false);     IRasterBand distanceRaster = costDistanceRaster.Item(0);     IRasterBand backLinkRaster = costDistanceRaster.Item(1);          IGeoDataset costPath = distanceOp.CostPath((IGeoDataset)destFc, (IGeoDataset)distanceRaster, (IGeoDataset)backLinkRaster, ESRI.ArcGIS.SpatialAnalyst.esriGeoAnalysisPathEnum.esriGeoAnalysisPathForEachCell);      A cost path analysis in ArcPy (using sa.CostDistance and sa.CostPath) takes approx 15-20 sec. Using the exact same inputs, the ArcObjects based routine takes 55-60 sec. Even using the .NET Geoprocessor is significantly slower than arcpy. So, I guess my questions here are:   1. Are the ArcPy and ArcObjects implementations pointing at the same code base (through their Python and .NET wrappers)?   2. Any tips to optimize ArcObject based Cost Path analysis?


title: How to make line width scale-dependent

text: Is there a way to tie line width to map scale in QGIS? All I found was an outside limit for display as the scale decreases. I'm using an outline line style for roads, and the roads are too big when I zoom out. What about converting them to a raster layer?


title: Is ST_GeomFromText better than providing direct geometry?

text: I have been working with postgis recently,and in my query if I use ST_GeomFromText it execute faster than running a sub-query to get geom. I thought ST_GeomFromText will be more expensive but after running many tests every time I got the result faster, my question Is there any explanation behind this? because for me getting the geom directly in sub-query is better than getting geom as text then added as GeomFromText. Thanks


title: How to calculate the average value of a raster within polygons?

text: I am using the GRASS plugin in QGIS. I am trying to calculate the average value of a raster dataset within a set of polygons. I have created a new mapset and have imported both the vector and raster layers into the GRASS mapset. I am now trying to calculate the set of statistics using the v.rast.stats tool. It starts to run and then comes up with the error _an error occurred while converting vector to raster_ and then _'raster map not found'_.               v.rast.stats raster=secant_grass@zonalstats     vector=plots_grass@zonalstats     colprefix=g/Applications/GRASS-6.4.app/Contents/MacOS/scripts/v.rast.stats:     line 190: 632 Segmentation fault v.to.rast in="$VECTORFULL"     out="${VECTOR}_$TMPNAME" layer="$GIS_OPT_LAYER" use=cat type=area     --quiet ￼An error occurred while converting vector to raster ￼Raster map not found ￼ nothing removed Finished with error      Any help would be greatly appreciated. Regards


title: Python and Data driven pages - Advanced use

text: I need to make a number of maps where only one input changes (a layer with the range of a radio tower) and I want to add the layer, apply symbology, zoom to it, change the title to include the name of the tower and then export it as a pdf with the tower name. Rather than adding the layers and applying symbology I can just put all the towers in a Group Layer and then the script can just turn all but one of them on, zoom, print and then turn it off and continue to the next tower. One long way is using: Data Driven Pages --> Page Definition option in the Definition Query tab of a Layer's Properties and specifying this but I was wondering if there was an easier work around? ![enter image description here](http://i.stack.imgur.com/s5r93.jpg) The issue with using Page Definitions is that some towers overlap and it means that I have to add in a column which will have have the tower name and signal strength repeatedly for each row. Given that we have more than 30 towers, even in a model this is laborious process that has to be then repeated when any base data changes.


title: Display problems in some layers of Geoserver 2.3.1

text: I deleted Geoserver 2.2-SNAPSHOT and deleted tomcat6. Set tomcat7 (from ubuntu repository) and Geoserver 2.3.1. Rewrote his workspce (store, style ..), changing the owner to tomcat7. Some layers are not displayed. The database remains exact same (PostgreSQL9.2 / PostGIS1 .5) Ubuntu server 12.04. Style               <?xml version="1.0" encoding="UTF-8"?>     <sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">         <sld:UserLayer>             <sld:LayerFeatureConstraints>                 <sld:FeatureTypeConstraint/>             </sld:LayerFeatureConstraints>             <sld:UserStyle>                 <sld:Name>Продорожная полоса</sld:Name>                 <sld:Title/>                 <sld:FeatureTypeStyle>                     <sld:Name>group 0</sld:Name>                     <sld:FeatureTypeName>Feature</sld:FeatureTypeName>                     <sld:SemanticTypeIdentifier>generic:geometry</sld:SemanticTypeIdentifier>                     <sld:SemanticTypeIdentifier>simple</sld:SemanticTypeIdentifier>                     <sld:Rule>                         <sld:Name>default rule</sld:Name>                         <sld:MinScaleDenominator>1.0</sld:MinScaleDenominator>                         <sld:MaxScaleDenominator>50000.0</sld:MaxScaleDenominator>                         <sld:PolygonSymbolizer>                             <sld:Fill>                                 <sld:GraphicFill>                                     <sld:Graphic>                                         <sld:Mark>                                             <sld:WellKnownName>hatch</sld:WellKnownName>                                             <sld:Fill/>                                             <sld:Stroke>                                                 <sld:CssParameter name="stroke">#808080</sld:CssParameter>                                                 <sld:CssParameter name="stroke-width">3.0</sld:CssParameter>                                             </sld:Stroke>                                         </sld:Mark>                                         <sld:Size>3.0</sld:Size>                                     </sld:Graphic>                                 </sld:GraphicFill>                                 <sld:CssParameter name="fill">#D95F02</sld:CssParameter>                                 <sld:CssParameter name="fill-opacity">0.5</sld:CssParameter>                             </sld:Fill>                             <sld:Stroke>                                 <sld:CssParameter name="stroke">#FFFF00</sld:CssParameter>                                 <sld:CssParameter name="stroke-dasharray">20.0</sld:CssParameter>                             </sld:Stroke>                         </sld:PolygonSymbolizer>                     </sld:Rule>                 </sld:FeatureTypeStyle>             </sld:UserStyle>         </sld:UserLayer>     </sld:StyledLayerDescriptor>      Geoserver log > org.apache.wicket.WicketRuntimeException: Error attaching this container for > rendering: [MarkupContainer [Component id = extraStyles]] at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1765) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.extensions.markup.html.form.palette.Palette.onBeforeRender(Palette.java:191) > at org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.markup.repeater.AbstractRepeater.onBeforeRender(AbstractRepeater.java:151) > at org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.extensions.markup.html.tabs.TabbedPanel.onBeforeRender(TabbedPanel.java:244) > at org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.markup.html.form.Form.onBeforeRender(Form.java:2056) at > org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > at org.apache.wicket.Component.onBeforeRender(Component.java:3946) at > org.apache.wicket.Page.onBeforeRender(Page.java:1531) at > org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.Component.prepareForRender(Component.java:2292) at > org.apache.wicket.Page.prepareForRender(Page.java:1521) at > org.apache.wicket.Component.prepareForRender(Component.java:2329) at > org.apache.wicket.Page.renderPage(Page.java:892) at > org.apache.wicket.protocol.http.WebRequestCycle.redirectTo(WebRequestCycle.java:186) > at > org.apache.wicket.request.target.component.PageRequestTarget.respond(PageRequestTarget.java:58) > at > org.apache.wicket.request.AbstractRequestCycleProcessor.respond(AbstractRequestCycleProcessor.java:105) > at > org.apache.wicket.RequestCycle.processEventsAndRespond(RequestCycle.java:1258) > at org.apache.wicket.RequestCycle.step(RequestCycle.java:1329) at > org.apache.wicket.RequestCycle.steps(RequestCycle.java:1436) at > org.apache.wicket.RequestCycle.request(RequestCycle.java:545) at > org.apache.wicket.protocol.http.WicketFilter.doGet(WicketFilter.java:484) at > org.apache.wicket.protocol.http.WicketServlet.doPost(WicketServlet.java:160) > at javax.servlet.http.HttpServlet.service(HttpServlet.java:641) at > javax.servlet.http.HttpServlet.service(HttpServlet.java:722) at > org.springframework.web.servlet.mvc.ServletWrappingController.handleRequestInternal(ServletWrappingController.java:159) > at > org.springframework.web.servlet.mvc.AbstractController.handleRequest(AbstractController.java:153) > at > org.springframework.web.servlet.mvc.SimpleControllerHandlerAdapter.handle(SimpleControllerHandlerAdapter.java:48) > at > org.springframework.web.servlet.DispatcherServlet.doDispatch(DispatcherServlet.java:923) > at > org.springframework.web.servlet.DispatcherServlet.doService(DispatcherServlet.java:852) > at > org.springframework.web.servlet.FrameworkServlet.processRequest(FrameworkServlet.java:882) > at > org.springframework.web.servlet.FrameworkServlet.doPost(FrameworkServlet.java:789) > at javax.servlet.http.HttpServlet.service(HttpServlet.java:641) at > javax.servlet.http.HttpServlet.service(HttpServlet.java:722) at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:305) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.geoserver.filters.ThreadLocalsCleanupFilter.doFilter(ThreadLocalsCleanupFilter.java:27) > at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.geoserver.filters.SpringDelegatingFilter$Chain.doFilter(SpringDelegatingFilter.java:74) > at > org.geoserver.filters.SpringDelegatingFilter.doFilter(SpringDelegatingFilter.java:45) > at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.geoserver.platform.AdvancedDispatchFilter.doFilter(AdvancedDispatchFilter.java:49) > at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.vfny.geoserver.filters.SetCharacterEncodingFilter.doFilter(SetCharacterEncodingFilter.java:109) > at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:311) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:68) > at > org.springframework.security.web.access.intercept.FilterSecurityInterceptor.invoke(FilterSecurityInterceptor.java:116) > at > org.springframework.security.web.access.intercept.FilterSecurityInterceptor.doFilter(FilterSecurityInterceptor.java:83) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:72) > at > org.geoserver.security.filter.GeoServerCompositeFilter.doFilter(GeoServerCompositeFilter.java:91) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:323) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:68) > at > org.springframework.security.web.access.ExceptionTranslationFilter.doFilter(ExceptionTranslationFilter.java:113) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:72) > at > org.geoserver.security.filter.GeoServerCompositeFilter.doFilter(GeoServerCompositeFilter.java:91) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:323) > at > org.geoserver.security.filter.GeoServerAnonymousAuthenticationFilter.doFilter(GeoServerAnonymousAuthenticationFilter.java:53) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:323) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:68) > at > org.springframework.security.web.authentication.AbstractAuthenticationProcessingFilter.doFilter(AbstractAuthenticationProcessingFilter.java:182) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:72) > at > org.geoserver.security.filter.GeoServerCompositeFilter.doFilter(GeoServerCompositeFilter.java:91) > at > org.geoserver.security.filter.GeoServerUserNamePasswordAuthenticationFilter.doFilter(GeoServerUserNamePasswordAuthenticationFilter.java:115) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:323) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:68) > at > org.springframework.security.web.authentication.rememberme.RememberMeAuthenticationFilter.doFilter(RememberMeAuthenticationFilter.java:146) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:72) > at > org.geoserver.security.filter.GeoServerCompositeFilter.doFilter(GeoServerCompositeFilter.java:91) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:323) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:68) > at > org.springframework.security.web.context.SecurityContextPersistenceFilter.doFilter(SecurityContextPersistenceFilter.java:87) > at > org.geoserver.security.filter.GeoServerSecurityContextPersistenceFilter$1.doFilter(GeoServerSecurityContextPersistenceFilter.java:46) > at > org.geoserver.security.filter.GeoServerCompositeFilter$NestedFilterChain.doFilter(GeoServerCompositeFilter.java:72) > at > org.geoserver.security.filter.GeoServerCompositeFilter.doFilter(GeoServerCompositeFilter.java:91) > at > org.springframework.security.web.FilterChainProxy$VirtualFilterChain.doFilter(FilterChainProxy.java:323) > at > org.springframework.security.web.FilterChainProxy.doFilter(FilterChainProxy.java:173) > at > org.geoserver.security.GeoServerSecurityFilterChainProxy.doFilter(GeoServerSecurityFilterChainProxy.java:103) > at > org.springframework.web.filter.DelegatingFilterProxy.invokeDelegate(DelegatingFilterProxy.java:346) > at > org.springframework.web.filter.DelegatingFilterProxy.doFilter(DelegatingFilterProxy.java:259) > at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at org.geoserver.filters.LoggingFilter.doFilter(LoggingFilter.java:75) at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at org.geoserver.filters.GZIPFilter.doFilter(GZIPFilter.java:42) at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.geoserver.filters.SessionDebugFilter.doFilter(SessionDebugFilter.java:47) > at > org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:243) > at > org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:210) > at > org.apache.catalina.core.StandardWrapperValve.invoke(StandardWrapperValve.java:224) > at > org.apache.catalina.core.StandardContextValve.invoke(StandardContextValve.java:169) > at > org.apache.catalina.authenticator.AuthenticatorBase.invoke(AuthenticatorBase.java:472) > at > org.apache.catalina.core.StandardHostValve.invoke(StandardHostValve.java:168) > at > com.googlecode.psiprobe.Tomcat70AgentValve.invoke(Tomcat70AgentValve.java:38) > at > org.apache.catalina.valves.ErrorReportValve.invoke(ErrorReportValve.java:98) > at org.apache.catalina.valves.AccessLogValve.invoke(AccessLogValve.java:927) > at > org.apache.catalina.core.StandardEngineValve.invoke(StandardEngineValve.java:118) > at > org.apache.catalina.connector.CoyoteAdapter.service(CoyoteAdapter.java:407) > at > org.apache.coyote.http11.AbstractHttp11Processor.process(AbstractHttp11Processor.java:987) > at > org.apache.coyote.AbstractProtocol$AbstractConnectionHandler.process(AbstractProtocol.java:579) > at > org.apache.tomcat.util.net.JIoEndpoint$SocketProcessor.run(JIoEndpoint.java:309) > at > java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1145) > at > java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:615) > at java.lang.Thread.run(Thread.java:724) Caused by: > java.lang.NullPointerException at > org.geoserver.wms.web.publish.StyleNameRenderer.getIdValue(StyleNameRenderer.java:21) > at > org.apache.wicket.extensions.markup.html.form.palette.component.Recorder.initIds(Recorder.java:116) > at > org.apache.wicket.extensions.markup.html.form.palette.component.Recorder.onBeforeRender(Recorder.java:87) > at org.apache.wicket.Component.internalBeforeRender(Component.java:1071) at > org.apache.wicket.Component.beforeRender(Component.java:1105) at > org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1753) > ... 140 more


title: Rendering order of layers in QGis libraries

text: I am using the QGis libraries in a external Qt App. When adding vector and raster layers (in any order) to the QgsMapLayerRegistry instance, it somehow sorts the layers so that the vector layers is first in the QgsMapLayerRegistry.mMapLayers list and then the raster layers. Which make perfectly sense otherwise at render time you wont see the vector layers if the order was the other way round. But when adding a vector point layer with a memory data provider, i would like the layer to be first in the list, or at least "sorted" with the other vector layers. But is always end up last in the QgsMapLayerRegistry.mMapLayers list, and at render time it always gets rendered "beneath" the raster layer, how can i add the vector point layer with the memory data model, to be rendered "on top" of all other layers?


title: ArcPy space in textfield

text: I got a tiny script:               def calc(value):      if value is '123':         return 'ok'      elif value is '12 3':         return 'not_ok'      else:         return value      if the value is without spaces it works perfectly, but if there is a space in it, it just returns the value. Where is my error? :)


title: How to multiple select and delete in ArcCatalog?

text: I have so many files to delete, to do that I am using ArcCatalog but I can do only one by one. Is there any way to delete many files at the same time. Also, I would like to know how to select more than files in Arc catalog as well.


title: How to count polygons in multipart features?

text: I'm trying to count the number of polygons in a multipart feature. My aim is to have a 'count' attribute to quantify fragmentation of multipart land parcels. I've looked at the field calculator and for suitable plugins, but can't see anything that would help. There are lots of ways to count points, lines and vertices, but not polygons.


title: Access to custom GeoJSON members in OpenLayers

text: I need to read a layer-wide attribute, from a GeoJSON file, set by a custom member to the object:               {"type": "FeatureCollection", "features": [...], "search_polygon": "POLYGON (...)"}      Since it seems neither the Vector Layer nor the Layer classes have any attribute or method to expose this data, I tried to use the callback from HTTP Protocol with no success. (I put an alert to the callback function, but no response.) My question is: Is there any way to access that `search_polygon` custom member using OpenLayers?


title: How to speed up an st_intersects() query?

text: I have a "huge" spatial database: 72GB, 200 millions documents with geometries (points) I have a GIST index on the geometries the database:               CREATE TABLE resources     (     id bigint NOT NULL,     accuracy varchar(3),     owner varchar(100),     title varchar(500),     posted timestamp,        nbtags int,     url text,     description text,     rawtags text,     tags text,     geom geometry NOT NULL,     CONSTRAINT resources_pkey PRIMARY KEY (id)     );      INDEX resources_index_rep ON resources USING GIST (geom);          CREATE RULE resources_rule AS ON INSERT TO resources     WHERE EXISTS(SELECT 1 FROM resources                 WHERE id =NEW.id)     DO INSTEAD NOTHING;      Of course I want to query this database. Nevertheless, it's very very slow ... (I did analyzes and Vacuums) for example with area of France:               EXPLAIN ANALYZE     SELECT count(*)      FROM resources     WHERE st_intersects(st_geomfromtext('POLYGON((-5.134723 41.364166,-5.134723 51.09111,9.562222 51.09111,9.562222 41.364166,-5.134723 41.364166))',4326), geom)          QUERY PLAN     Aggregate  (cost=14411182.45..14411182.46 rows=1 width=0) (actual time=17368911.185..17368911.186 rows=1 loops=1)      ->  Bitmap Heap Scan on resources  (cost=890493.27..14391392.83 rows=7915848 width=0) (actual time=2627029.235..17357552.440 rows=21381758 loops=1)         Recheck Cond: ('0103000020E6100000010000000500000077137CD3F48914C05628D2FD9CAE444077137CD3F48914C0562B137EA98B4940F52EDE8FDB1F2340562B137EA98B4940F52EDE8FDB1F23405628D2FD9CAE444077137CD3F48914C05628D2FD9CAE4440'::geometry && geom)         Filter: _st_intersects('0103000020E6100000010000000500000077137CD3F48914C05628D2FD9CAE444077137CD3F48914C0562B137EA98B4940F52EDE8FDB1F2340562B137EA98B4940F52EDE8FDB1F23405628D2FD9CAE444077137CD3F48914C05628D2FD9CAE4440'::geometry, geom)         ->  Bitmap Index Scan on resources_index_rep  (cost=0.00..888514.31 rows=23747545 width=0) (actual time=2626090.609..2626090.609 rows=21381788 loops=1)               Index Cond: ('0103000020E6100000010000000500000077137CD3F48914C05628D2FD9CAE444077137CD3F48914C0562B137EA98B4940F52EDE8FDB1F2340562B137EA98B4940F52EDE8FDB1F23405628D2FD9CAE444077137CD3F48914C05628D2FD9CAE4440'::geometry && geom)     Total runtime: 17368946.473 ms     7 row(s)      it took almost 5 hours. and with Switzerland:               EXPLAIN ANALYZE     SELECT count(*)      FROM resources     WHERE st_intersects(st_geomfromtext('POLYGON((5.96611 45.829437,5.96611 47.806938,10.488913 47.806938,10.488913 45.829437,5.96611 45.829437))',4326), geom)          QUERY PLAN     Aggregate  (cost=5623425.68..5623425.69 rows=1 width=0) (actual time=6798842.064..6798842.064 rows=1 loops=1)     ->  Bitmap Heap Scan on resources  (cost=88982.29..5621452.65 rows=789212 width=0) (actual time=783383.496..6797250.482 rows=2539057 loops=1)         Recheck Cond: ('0103000020E61000000100000005000000AF5A99F04BDD1740D28BDAFD2AEA4640AF5A99F04BDD174028F38FBE49E74740B22D03CE52FA244028F38FBE49E74740B22D03CE52FA2440D28BDAFD2AEA4640AF5A99F04BDD1740D28BDAFD2AEA4640'::geometry && geom)         Filter: _st_intersects('0103000020E61000000100000005000000AF5A99F04BDD1740D28BDAFD2AEA4640AF5A99F04BDD174028F38FBE49E74740B22D03CE52FA244028F38FBE49E74740B22D03CE52FA2440D28BDAFD2AEA4640AF5A99F04BDD1740D28BDAFD2AEA4640'::geometry, geom)         ->  Bitmap Index Scan on resources_index_rep  (cost=0.00..88784.99 rows=2367636 width=0) (actual time=782664.847..782664.847 rows=2539063 loops=1)               Index Cond: ('0103000020E61000000100000005000000AF5A99F04BDD1740D28BDAFD2AEA4640AF5A99F04BDD174028F38FBE49E74740B22D03CE52FA244028F38FBE49E74740B22D03CE52FA2440D28BDAFD2AEA4640AF5A99F04BDD1740D28BDAFD2AEA4640'::geometry && geom)     Total runtime: 6798852.955 ms     7 row(s)      it took almost 2 hours. I'm looking for how to improve the performance? About our server:               - CPU: 4     - RAM: 16 GB     - OS: Ubuntu 12.04 LTS 64bit


title: compound variables in model builder?

text: I'm trying to create a model wherein one of the parameters is a compound variable, that is made of multiple variables joined together. How can this be done? For Example: Create an insider buffer using _Buffer(analysis)_ and name the output feature class from the distance parameter. Variables:               BufferDistance = 1000     InputFC = d:\test.gdb\boundaries      Buffer tool parameters:               Input Features = %InputFC%     Output Feature Class = %InputFC%_%BufferDistance%     Distance = -%BufferDistance% Meters      A negative buffer distance is needed to create an inside buffer, but that means the %BufferDistance% variable can't be used in the output FC name. Thus the need for a compound variable.


title: Show vehicle traveled route as animation

text: > **Possible Duplicate:**   >  Animate feature in openlayers I create map web app with openlayers and yii php framework.   I want to show vehicle traveled route as animation like below site in my application.   http://mstrackweb.mobil-i.com   In the above site use demo/demo for login to site.In the history tab you can see traveled route of vehicles as animation.   In my app i can add markers to map after 2seconds with below code.In the history array i have lon/lat and dateTime of the vehicle as location history.               var idx = 0;     var interval = setInterval(function(){         if(idx >= history.length)             clearInterval(interval);         else{             markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(                 history[idx]['longitude'],history[idx]['latitude'] ).transform(             new OpenLayers.Projection("EPSG:4326"),             new OpenLayers.Projection('EPSG:900913')                 ),icon.clone()));             idx++;         }     },2000);      But if a marker is out of the map bound i don't know how to change the bounds of the map to showing added marker!.   how can i change bounds of the map to showing added marker?   is this approach is true? When a marker added to map create event and if this marker isn't in the map bounds change map bound to showing it.   Is there any library for showing vehicle traveled route as animation?


title: Qgis crash on ubuntu precise

text: ubuntu precise qgis 1.8 I open a shape file and qgis crashes hard. i.e. it completely disappears from my desktop. qgis stable version is what i am using. (1.8)


title: How do I split polygons that cross the dateline?

text: I have a dataset of country borders which I obtained from: http://thematicmapping.org/downloads/world_borders.php/ I have that dataset in both shp and SQL, but I'll put out there that I don't feel that affects my problem; however, please let me know if it does/can. Because I'll be displaying those borders in something that uses Google Maps, I've got a problem. That problem is that countries whose polygons cross 180 degrees are displaying extremely weird in Google Maps. Either it's splitting them completely wrong, taking them all the way across the map, whatever the weirdness might be. So, I thought I can simply address this by splitting the polygons of the countries that cross the line so that no country contains a polygon that presents a problem. Examples of such countries/regions would be Russia and Antarctica. I therefore invoked the following: ogr2ogr -wrapdateline -t_srs EPSG:4326 new countryBorders.sql I got: ERROR 1: TopologyException: found non-noded intersection between LINESTRING (-59.7742 -72.9003, -59.8494 -73.2331) and LINESTRING (0 -73.1221, -89.3211 -73.0542) at -59.8141 -73.0766 ERROR 1: TopologyException: found non- noded intersection between LINESTRING (-59.7742 -72.9003, -59.8494 -73.2331) and LINESTRING (180 -73.2589, -89.3211 -73.0542) at -59.8141 -73.0766 A few questions: 1: am I even doing this the right way? for example, I find it odd that I have to provide ogr2ogr the same srs as of the input. Is that translation even needed? Is that even the right SRS to use? I obtained that SRS from taking the .prj file and looking it up on this site: http://prj2epsg.org/search 2: all I want is for a few of the polygons of those few troublesome countries to get split so they don't cause any more problems. -wrapdateline seemed like the perfect solution ... maybe it's not? I'm completely open to advice, suggestions, references, whatever. I am perfectly willing to redo this dataset in another projection, run a one- time query against it to redo it so it is more compatible with displaying on Google Maps, anything at all. I'd love any advice. It doesn't have to involve my most probably incorrect approach outlined above. Thanks so much in advance!


title: Automatic Watershed Delineation

text: I am trying to delineate watersheds automatically using Arcgis/Achydro. I got the drainage area contribution of stream segment (part of stream) as shown in figure stream. The corresponding drainage area as shown in figure area. But I am looking for drainage areas/catchment areas as shown in figure Water Shed. In watershed figure , at any point on the stream, it calculates drainage area by considering entire upstream length of the stream. Suppose, if you choose point "8", its catchment/drainage area covers the area of points "9" and "10". Similarly for point "5" covers all points upstream of stream. I want to get similar automatic watershed delineations. Please guide me to the software/method to get the result as shown in figure watershed for my entire stream network of an area. ![Stream](http://i.stack.imgur.com/z8CBY.jpg) ![Area](http://i.stack.imgur.com/2VE8t.jpg) ![WaterShed](http://i.stack.imgur.com/hENfm.jpg)


title: Using M 'Measure' value in SQL Server 2008

text: MS Sql Server 2008 geometry type supports adding Z and M values to geometry objects, including linestrings. However, I cannot find any way to make use of these values, since methods such as `.STIntersection()` return a 2D point. My end goal is to find where a given point intersects with a linestring and then find the equivalent Z and M values. Googling and scouring Technet has yielded no useful information other than to make me wish I was using PostgreSql instead, which seems to include the functionality I require (unfortunately I am constrained to Sql Server). Does anyone know how to actually make use of the Z and M values within linestrings, or of any toolsets that provide this functionality?


title: Remove lock on feature class

text: When i try to delete a feature class , after or during the use of an Arcgis engine app which uses the same feature class i got an error concerning the lock , so i can't get ride of this lock until i restart the computer. i would like to know a way to remove locks using arcpy or arcobjects 10.1 , it would be perfect if i can get ride of the lock inside the arcengine app which also locks the geodatabase , i have seen different posts in SE about locks but no one of them gives the solution to free locks using Arcobjects. Any help please


title: ArcMap 10.1 Python split inconsistent string

text: I converted a kml file to a feature class. All of the fields from the kml were combined into two fields in the feature class. One of them contains information delimited by "/" and the other contains information delimited by "br". I need to extract all of the information into separate fields. I met with some luck using the methods described here: ArcGIS 10 Python split string The problem is that the nothing is consistent about the information in the fields except for the delimiter. Some rows only contain one piece of information (no delimiter) and others contain up to five. When I use .split with a [0] to return the first piece of information it works fine, when I try to return the second piece [1] I get an error. If I get rid of the data that does not contain a delimiter then [1] works, but I don't know of a way to get rid of all of it. I am very new to ArcMap, and I am aware that I'm probably not using the correct terminology. I have read a lot of forums, looked at python dictionaries ( for example: http://docs.python.org/library/string.html) and tried to make this work over the course of many hours. The bottom line is that I still don't understand Python well enough to come up with something that works on my own! Please help me!! I tried to post images showing what I have tried and what my data is like, but I can't because I am a new user. If you think you have any ideas I would love to share more info. Thanks!


title: Query a polyline layer using arcGIS's Javascript API

text: I'm attempting to query a polyline layer using the ArcGIS Javascript API I all ready have a lot of experience doing this with Point layers. So I set my Polyline layer up the same way and stopped and started my service again to get every thing going. When I run the query (specifically when I run PolylineQueryTask.execute()) I get an error that just says the operation can't be completed. There is some code below that shows how I'm setting my query layer up. Any ideas on what I can do to fix this issue or things I should try?               PolylineQueryTask = new esri.tasks.QueryTask("http://" + DomainAddress + "/ArcGIS/rest/services/SNND/MapServer/1");          PolylineQuery = new esri.tasks.Query();     PolylineQuery.returnGeometry = true;     PolylineQuery.outFields = ["ObjID"];          PolylineQuery.where = 'where 1 = 1';     PolylineQueryTask.execute(PolylineQuery, ShowPolylineResults, FatlError);      So as you can see from my code I should be getting all my polylines back. Any thoughts are very appreciated.


title: Can anyone point me to good development environment in linux for client side search with flex or javascript?

text: I want to create a client in mapfish, and/or opengeo, using open layers, maybe cartaro and would like to play in a development environment using flex or javascript. I need to take my data, make it searchable from table columns and features, and be able to search strings for geolocation content, any good tuts, links, or linux development environments that you could recommend?


title: How can I fix QGIS 1.8 startup issues due to missing qwt5.dll library?

text: I installed the QGIS 1.8 a couple of weeks ago and it run ok. Today the qwt5.dll library is missing apparently. I deinstalled and also deinstalled the microsoft visual c++ 2005 distributables to do a clean new install but nothing seems to work. Windows 7 64 bit .exe


title: How can I programmatically display the "Create Features" dockable window in ArcMap 10?

text: When I start editing in ArcMap 10, by default the "Create Features" window is automatically opened and docked. Stopping the edit session hides this window. When I use the IEditor:StartEditing method this does not display the "Create Features" window however if it is already displayed the IEditor:StopEditing method hides it. How can I Start Editing programmatically and display the "Create Features" window? ![alt text](http://i.stack.imgur.com/wUgWw.jpg)


title: I have several layers (polylines and polygons) that are of the same area but they don't overlap in ArcMap as they should, what do I do?

text: For my GIS project for my class, I have downloaded several vector layers from different websites (polyline of the continental divide, a boundary polygon of New Mexico, a rivers and water bodies layer of the United States, etc.). Obviously, these features in real life are associated with each other but when I import each of these shapefiles into ArcMap they are not overlapping as they should (i.e the continental divide does not cut through the U.S. polygon or the New Mexico polygon). I noticed that the scales for each shapefile are completely different from each other (i.e. 1:143 vs. 1: 4,203,394), is this the problem and if so does anyone know who to fix this, and if not, does anyone know what the problem I am having actually is and how to fix it.


title: Intergraph TerramapServer and PostGIS

text: Has anyone ever used the two of them together, and if so, any special caveat to be aware of ? Furthermore, does it even make sense to use PostGIS / do TerramapServers have their own storage format instead ?


title: Register as versioned with move edits to base - when not to?

text: What limitations do I get, when I chose to register as versioned with the option to move edits to base?


title: Why doesn't the gml 3.1.1 schema validate with xmllint?

text: Running xmllint on the gml 3.1.1 schema fails, even though the release notes claim that this release was created to get the schema to validate using 'common tools'. There's questions about this on the web, but no satisfying answers.               xmllint --schema http://schemas.opengis.net/gml/3.1.1/base/gml.xsd \             --noout foo.xml      returns something like: > element complexType: Schemas parser error : local list type: A type, derived > by list or union, must have the simple ur-type definition as base type This is not the same problem as the one described here, which was a Xerces C++ issue. ~~If there is some other xmllint issue here, what xml validator is known to work with GML?~~ I have this problem on a few different flavors of Ubuntu, CentOS, and RHEL 5.x **Update:** The 3.1.1 and 2.1.2 GML schemas validate using the following tools:   1. The W3C XML Schema Validator   2. Sun's msv:  running `java -jar msv.jar -warning gml\base\gml.xsd` returned 1 warning about "anyType" in gmlBase.   3. Using the validator bundled in eclipse worked with no errors or warnings.


title: 'GIS' grade mobile mapping app that isnt ArcPad?

text: I am a GIS analyst/geologist with a mining company based in Canada. However most of my fieldwork is in remote places with minimal access to the internet. I am after a reliable alternative to ArcPad and ESRI running on either ipad or android tablets. I have already tried several apps like Mapitfast and cartomobile but Ive been unimpressed with either their UI and/or inability to share/access data off the device. Some tools I need are:   * Creating rock outcrop maps   * Editing/importing geology polygons   * GPS tracking that can be exported for later use   * Photo-tagging would be a plus   * Bing or Google as a background image Thanks in advance. -David-


title: Extract pixel value for corresponding gps location in R

text: I have animal gps data with several locations for each day and regular gaps of several days when no location of the animal was recorded. Additionally, I have satellite data in intervals of 16 days. I would now like to **extract the pixel value** that corresponds to the **specific point** and to the **specific time**. This means if the location of the animals was recorded e.g. 2 days before a satellite image was taken, I want to extract the pixel value of this image (that comes after) rather than from the images that was taken 14 days before the animal location was recorded. I always want to extract from the image that lies closer according to the time. I created some test data that hopefully illustrates the problem:               library(sp)     library(raster)          ### Create test data     # create first raster     edc2012001_m <- raster(ncol=36, nrow=18)     edc2012001_m[] <- sample(1:ncell(edc2012001_m))          # create second raster     edc2012017_m <- raster(ncol=36, nrow=18)     edc2012017_m[] <- sample(1:ncell(edc2012017_m))          rasters<-stack(edc2012001_m,edc2012017_m)          # Create xy coordinates     time<-c("2012-01-01", "2012-01-01", "2012-01-01", "2012-01-02", "2012-01-02", "2012-01-02", "2012-01-12", "2012-01-12", "2012-01-13", "2012-01-13")     x <- rep(-50,10)     y <- sample(c(-80:80), 10)     data<-data.frame(x,y,time)          # Convert data to spatial points data frame     coordinates(data) <- c("x","y")          ### Extract all data from raster stack     extract(rasters, data)      The number string in the name of the first raster indicates that the image was taken on the first day in 2012, the second image was taken on the 17th day of this year. The 7th location in the test data for example shall now be extracted from the second raster file as it lies closer according to the time. All in all I have 87 raster files and 600 observations. I am really not sure how to program this. I guess I could use `substr()` to retrieve the date info from the raster names. But besides that... I am thankful for every hint I can get, also about functions that might be helpful in this context.


title: Using python in ArcGIS, create a line of varying distance from an origin point

text: What I'm trying to do is create a line feature from a single point, using a set distance and angle. Using ArcGIS and Python (I hope). Here's the scenario: Say I have a point at these coordinates: X = 400460.99, Y = 135836.76 From this point, I want to create a 800 Meter long line at a 15 degree angle from that point. I do not know what the resulting endpoint will be, which is frustrating. (My data are projected in Maryland State Plane South - Meters) Thanks for any help and suggestions. James


title: How to display my country's map?

text: I am newbie in Geomajas and I tried to create new Geomajas application using the documentation in http://files.geomajas.org/maven/trunk/geomajas/docbook- gettingstarted/html/master.html. I don't know how to display my country's map (Mongolia) in my created Geomajas application. Can you tell me how to do that or do you have any step by step guide/tutorial on this issue? Please help me


title: Converting tax map polygons from Shapefile to table of map number and corner coordinates

text: Fairfax County, Virginia provides a Shapefile here which shows the boundaries of the tax maps as polygons. Each polygon simply consists of four corners, and what I'd like to do is create a table showing the map number and the coordinates of the corners, such as the following: Map number, bottom-left, bottom-right, top-right, top-left   77-1, -77 32, -76 32, -76 33, -77 33 (though they aren't quite square like this) I don't have access to ArcGIS software, but do have QGIS and OpenJUMP, and can install any other software that is freely available. I believe this should be very easy but I'm just missing something. Thanks!


title: Joining lines when direction is not known

text: I am struggling with this problem. I have a polyline made up of several segments. Each segment is made up of points, each point identified in 3D space by coordinates and elevation. Together, if plotted, the segments make up a more or less continuous line (there could be breaks between segments), but the segments themselves are not consecutive, nor do the points of all segments follow the same direction of travel.   The question is: how can I create, using Python preferably, a single line from these non-consecutive segments so that I can measure the distance between points and the total length of the line.   I don't even know which of the segments is the first one or the last one in the line, but somehow I have to put them in the right sequence and make sure they all point in the same direction so I can measure them.   thank you for any assistance. I'll be happy to supply additional info, data, etc.. I stress I am not asking for the actual python code (not that I would refuse it...), just the strategy. Bob


title: How to generalize road network to be exact angle to each other?

text: I need to generalize a road network to be an exact angle to each other, we are talking about parallel streets like New York. The thing is that when you generalize, you do it per feature. I am trying dissolve and unsplit, but then it goes past the parallel streets and the generalization destroys curved streets. I do not want to do it manually, it is a line feature. think of one road that is straight, but the data is not. all line comes as splited in intersections. when generalizing, the output is not generalizing the hole straight road, so the road will still be full of vertexes. but when I dissolve the generalization gets better, but then if one street is changing direction at the end of the block, the topology disappears . look at the link :)https://docs.google.com/file/d/0B8sJ0k30rFr2SkNWVllWMDMxY1k/edit?usp=sharing https://docs.google.com/file/d/0B8sJ0k30rFr2TTNuLWpaNlJFY1k/edit?usp=sharing so on the brown lines you can see that within the red area the block is basically rectangular, but the end curves, so when generalising, the curve gets vertexes removed, and now it is really not straight –


title: Calculating slope in percentage

text: I'm working with qgis, v. 2.0.1. I have several soil types. They are (amongst others) described with slope indications. I already calculated the slope from the dtm in degrees, but as de description of the soil types is in percentage, I need to calculate the slope in percentage. I tried two things: first I tried the analysis tool from the raster section, but it didn't work. The result was just black with values from 0 to 0. The same happened when I tried to calculate the percentage in the raster calculator but the result was the same. Maybe I need to say that the terrain is pretty plain. What else can I do? Thanks. UPDATE: I tried the raster calculator and did what radouxju suggested. It gave me a black picture with values from 1756.64 to 1756.64. But as soon as I click in the picture the values are "normal", say 5.88 or 13.224 or something like that. The same happens when I try the gdaldem tool but with values from 0 to 0. I also tried to change the contrast enhancement (like strech to min/max) but it didn't change anything. Update 2: I just found out that maybe i have to set the scale to 111120.(http://learninggis.com/working-with-dem-models-qgis-quantum-gis) but that doesn't work neither (result: black, values 0-0).


title: Error with existing OSM QGIS style

text: I'm trying to use existing OSM style file like those of AsgerPetersen (https://github.com/3liz/osm-in- qgis/tree/master/styles/openstreetmap/geofabrik) or anitagraser's (https://github.com/anitagraser/QGIS- resources/tree/master/qgis2/osm_spatialite) but always got an error = * > " unexpected character à la ligne 31 colonne 72 " I'm on mac osX QGIS 2.2


title: SQL Server 2012 spatial check constraints disabled after migration from 2008

text: Currently we are using SQL Server 2008 geometry datatyped columns in our databases. These are registered with ESRIs enterprise geodatabase repository, so every spatial table has a check constraint which checks if the SRID of every geometry instance is equally               ([Shape].[STSrid]=(31256))      After restoring the 2008 database backup file to SQL Server 2012 for migration purposes every single of these check constraints is disabled (because of changes MS did on the spatial CLR types). As far as I can see, the interface to the CLR types is the same in this case, so I should simply be able to reenable the constraints and everything should work fine again. Are there any major issues I could run into, if I do this? At this point I can't see any obstacles.


title: Multivalue iteration in Vehicle Routing Problem (VRP)?

text: I'm currently building a model that would simulate three different set of orders (shapefile format) by using the **Vehicle Routing Problem (VRP)** tool of network analyst. The model I built runs correctly but shows me only the first delivery solution in the excel file. My aim is to have all the different route results into an Excel file (or even multiple). Does anybody know how to do that or fix the model? I also add that one issue is related to the fact that I cannot specify the input data in the delivery "Add location" during different iterations (only the first one is considered). ![enter image description here](http://i.stack.imgur.com/vWzXx.jpg)


title: Zero filling attributes using python

text: Being the AML dinosaur that I am, I am just learning python and I need to convert a 2 digit numeric attribute in a table to a 2 digit zero-filled attribute in a script. Example 9 needs to be 09. Any help would be appreciated.


title: How to convert GTFS to XML or JSON?

text: Do you know any tool that provides that? I need to manipulated GTFS file but I would like to have data available in any well-known format like XML or JSON. EDIT GTFS are Google Transit Feed Specification


title: I exported a DXF from QGIS, but it won't open in AutoCAD

text: I'm new in QGIS and here is my problem: I exported a DXf from QGIS by right-clicking on the layer. Then tried to open in AutoCAD 2000. That did not work. AutoCAD gave a failure-message: `Unknown Value "AC1018" found in drawingversion.` Then I tried to view the file in DWG-Truview2013 - which didn't work either giving the following error: `Invalid database handle 0 on line 1548. Invalid or incomplete DXF input -- drawing discarded.` Does anyone know a solution? Or is there anywhere in QGIS change the DXF- Version to export? Any Help will be useful. Thank you! Regards, Frank


title: Convert Shp to Dxf with QGIS

text: 8.0. I have asked this before and seen similar but not the same (all of which I Have not got an answer to so not a duplicate). But is there a way to export a shp as a dxf so that it will not crash AutoCAD when I open it.


title: Fixing "spatial index grid size is invalid" without Recalculate button

text: I've got a python script that is connecting a squence of points to make lines, and putting them into a feature class (based on answers to my previous question). It works wonderfully on a small area that I was using to debug and test (~25 meters, 20 points), but then when I put it into a real-world scenario (~1000 meters, ~500 points), it breaks.               The spatial index grid size is invalid.      Looking into the ArcGIS help suggestions to fix/prevent that, this page suggests recalculating the spatial index grid (directions here), but my polyline feature class doesn't have that option when I am in the "Indexes" tab -- there's no "Recalculate" button at all, and the "Delete" button is grayed out. So how can I fix the error? (I tried simply deleting my original feature class and making a new one, on the theory that my original "little" test case was setting the spatial index grid too small; however, even with a brand-new feature class, I have the same error.) The problem is more specifically the number of dots I'm trying to connect during one run of the script -- less than ~50 is fine, more than ~100 is not...* ArcGIS for Desktop with Advanced level license, data stored in feature class in a file geodatabase.


title: How to reclassify raster data on a scale of 0-1?

text: I have various raster layers generated after using euclidean distance tool.   How can I reclassify them on a scale of 0-1 instead of 1-10, etc?


title: GeoServer and Dynamic ColorMap Extension

text: I am running GeoServer 2.4.3 on Centos 6 using Tomcat6. Java version 1.7.0_51 OpenJDK I am trying to install and use the DynamicColorMap community extension. I extracted the jar in to WEB-INF/lib and restarted the server as usual. On the restart I am receiving this error:               2014-04-10 21:14:42,655 INFO [geoserver.wps] - Blacklisting process ras:DynamicColorMap as the input colorRamp of type interface org.geotools.styling.ColorMap cannot be handled      Thoughts? I have cross posted this to the GeoServer user list.


title: How do you change starting point on a polyline? (split tool)

text: There must be an easy fix for this. I am using the split tool and I want it to split a polyline into equal parts. For some reason, it is starting in the middle of my polyline. How do I force it to start where I want it to? I am using ArcGIS v. 10.1.


title: Why do I get an empty shapefile when converting this KMZ in QGIS?

text: I have downloaded a KMZ file containing global elevation information. I have tried to convert the data to a shapefile using qgis (save as shape option) and using the zonum solutions website but both times the conversion has not worked. When using qgis to convert, the new shape file was empty and all attributes were 'NULL'. Any idea on how to get around this? The KML/KMZ files I am using is from the Marine Geoscience Data system and is available from their website on the following link http://www.marine- geo.org/link/download.php?data_uid=14255&client=DataLink.


title: OpenLayers.Format.GeoJSON throws an exception "SyntaxError: JSON.parse: unexpected character" for valid JSON

text: In my GeoExt client, I have the following line that retrieves JSON data from a web service into a FeatureStore:               var store = new GeoExt.data.FeatureStore({             layer: route_layer,             fields: [                 {name: "length"}             ],             proxy: new GeoExt.data.ProtocolProxy({                 protocol: new OpenLayers.Protocol.HTTP({                     url: "./php/pgrouting.php",                     format: new OpenLayers.Format.GeoJSON({                         internalProjection: epsg_900913,                         externalProjection: epsg_4326                     })                 })             })         });      When this code is executed, Firebug reports the error `SyntaxError: JSON.parse: unexpected character`, even though the JSON returned by the web service is valid (I've confirmed this against http://jsonlint.com). The response header has been configured accordingly with the line `header('Content-type: application/json',true);`. Anyone who might have encountered this issue before, please assist. I've waded around, but no luck yet. **UPDATE:** Sample JSON output is at http://pastebin.com/LPLjBDgp. It's valid JSON, but the error `SyntaxError: JSON.parse: unexpected character` persists. I've even updated my html page by adding the following line `<meta http- equiv="Content-Type" content="text/html; charset=UTF-8"/>` in an attempt to restrict character encoding.


title: Normalize with gdal_translate

text: I do not understand the command line when I try to work with gdal_translate in QGIS. I found this at the forum: Normalize raster values to 0-255 And tried this: * **"If not, from within QGIS choose Raster -> Conversion - > Translate. Choose your files and options. At the bottom of the window the full gdal_translate command will be displayed. Click edit and add the -scale option."*** When I type the following command into the line I get an error message "wrong parameters"??? Maby the commands are wrong. ![enter image description here](http://i.stack.imgur.com/kHhvt.jpg) My current grid file has a range from 1 to 9. I need to have a range from 1 to 5 (5 classes).


title: Importing attributes from Excel file into new Polygons

text: I have a constantly updated Excel file containing address information which I am currently copying into the attributes of new polygons which I am creating on a daily basis. I am currently doing this my manually by copying the address across from the Excel file each time a add a new polygon. I am looking for a way of automating this process. After producing my polygon I would like a drop down box/ or search function of a unique ID which once selected will fillout the correct attributes by pulling them from the Excel file. I am new to arcgis 10 and would have previously used the VBA editor to have a go at this (no expert). However, I see that unlike the previous versions there is not longer and VBA editor and has been replaced by python, which I am unfamiliar with. I have also spotted the model builder tool, not sure if this could be used to solve my problem? So really I'm just asking for a point in the right direction. Is there any scripts available to solve a similar problem which I can modify? Many thanks Matthew


title: From gpx to line segment 1Km using ArcGIS 10

text: I have gpx files containg points with attribute data such as time and date. I'd like to convert them into lines, and that is not a problem, but what I want to do is create 1Km segments with time and date attribute. I'm using Arcgis 10.2 Thanks ![enter image description here](http://i.stack.imgur.com/UcAsg.png)


title: Leaflet - how can I call custom js after drawing polygon?

text: I want to draw a polygon on a leaflet map and afterwards send it to the DB and retrieve all objects that intersect that polygon. In Openlayers I did it by using custom function in "featureAdded:" when declaring OpenLayers.Control.DrawFeature. How can I do the same in leaflet ? Also , Is there an easy way to add custom controls from top right to top left as in openlayers control panel ? thanks.


title: Geoserver and SRS

text: I have a very basic problem, and please bear with me since I'm new to GIS. I have geoserver installed in a Ubuntu Box and installed a few maps in shapefile. The maps were originally in ArcGIS format .mdb, but I converted them to .shp file using ogr2ogr in Windows. The maps, as known from the .prj file is using LCC projection. and geoserver correctly identifies them in `native srs` as `lcc`. But it mandates me to use `EPSG:3034` in `Declared SRS`, which again forces me to use the same SRS via OpenLayers. Isn't it the duty of WMS Server to convert projections. Shouldn't I be allowed for asking geoserver to give me map tiles in say `EPSG:4326` projection?? Please give me an insight..


title: Can ArcPy's AsShape be used with GeoJSON FeatureCollection of Polygons with Labels?

text: The AsShape documentation provides examples that work very well. However, I have been given what I believe to be valid GeoJSON that is a FeatureCollection of two polygons each with a Label. It validates fine with the JSON Formatter & Validator but when I try AsShape on it, as per the code below, it throws an AssertionError on FeatureCollection. I'm assuming that is because AsShape does not support FeatureCollections and so am wondering if there is another approach that I could use?               import arcpy          polyX2_string = {"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[151.669111722222,-28.8913181547568],[151.669712422222,-28.890733765868],[151.670252222222,-28.8896441436458],[151.670847433333,-28.8891966436458],[151.670976977778,-28.8865389769791],[151.670112977778,-28.8846608880903]]]},"properties":{"Label":"National Park"}},{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[151.711278388889,-28.8201784658692],[151.710854355556,-28.8192543325359],[151.7100674,-28.8187800769803],[151.709037344444,-28.8184233992026],[151.707838755556,-28.8182234992026],[151.705821266667,-28.8180901992026],[151.705597933333,-28.8194629992025],[151.705198088889,-28.8212781880914]]]},"properties":{"Label":"National Park"}}]}          polygon = arcpy.AsShape(polyX2_string)


title: How do you produce a output table containing distances between mineral grains at various azimuthal directions?

text: I am using **ArcGIS 10** , the newest version and I am a senior at Western Carolina University and working on my senior thesis and the data set in which I have are images of mineral thin sections all separated into different groups of grains, magnetite, quartz, etc.. either being a vector layer or a raster layer. I am wanting to calculate the distance between each grain or the grain spacing, boundary to boundary, **NOT** center to center and also have those distances given every 10 degrees on a 0 to 360 degree directional scale. Is there any code out there that does this very exact method? Or is there another program, R or Fragstats, that would do a better job, ArcGIS 10 preferably. An example of the grains is shown here: ![enter image description here](http://i.stack.imgur.com/iIsQl.jpg) Thanks, -Michael Stillwagon


title: Access Control List for web GIS

text: I need to implement Access Control List (ACL) for my web application and I need to know common practice of using permissions in GIS systems. What levels of acces exists: layer, feature, attribute, extent, maybe anything else? What types of activities can be used by users: editing, browsing? Any response would be appreciated.


title: ArcGIS10 Split Line at Point Failing due to Dataset size

text: Im having an Issue with my Split Line at Point tool. Seems to run if lucky for 700 entries before arcmap/arccatalog falls over. Debugging of the errorlog has produced nothing. The total size of the line layer is 400,000 and the point layer 450000 covering the size of most states. Is there anyway I could set a python script to pull a row from current mxd, run the split line to point, to the point layer held in ram then append it to a shape/fgdb file, then loop back until complete. I have tried the various logical issues of confirming projections the same, repairing geometry, even splitting the data into 10% chunks/FCs. Also tried FME to the same result that it appears to run out of memory. Cheers. EDIT> At the moment Im tilling using Iteration by rowid using a dummy field that is populated with suburb names (spatial join population). This keeps the rows selected under 2000 each pass and at this stage seems to be working. 26/03/2012


title: How to find wmsURL? And how to use QueryTask method?

text: I have couple of questions.   1. var wmsLayer = new esri.layers.WMSLayer(wmsURL); with this line of code i can display wms layer over any base Layer. Now From where I get wmsURL ? I have published my shapefile on GeoServer, but i dont know how i get the URL (wmsURL) to use as parameter in above line of code ?    2. Now, Let's suppose I am successful to display wmsLayer. Now i want to query data against this wmsLayer. Like code below:  queryTask = new esri.tasks.QueryTask("wmsURL"); Can I use wmsURL as a parameter in QueryTask method ? From ArcGIS API reference, I come to know that the URL parameter of QueryTask method is the URL to the ArcGIS Server REST resource. But as u know i am using GeoServer, so how it could be fixed. Need some explanation ?


title: OD cost matrix Arcpy

text: I'm trying to create an OD cost matrix python script. I'm doing an OD cost matrix analysis. I'll find out the average time with patients to hospitals over 10 years. When I do this analysis, the result is that each patient visit every hospital, but in reality a patient only visit one hospitals. In my orgins table I have an attribute value HOSPITALS with the names of hospitals that can be linked to attribute value HOSPITALS in my destinations table. This means that the analysis should understand that each patient only visit a hospital. Does anyone have a tip on how I do this in arcpy?


title: Why does arcpy.ListFields list GUID field as SmallInteger?

text: I want to get a list of all GUID fields in a feature class, but it seems that arcpy.ListFields is listing my GUID field as a SmallInteger. What am I doing wrong?               global_id_fieldnames = [f.name for f in arcpy.ListFields(layer) if f.type.lower() == 'guid']      after this list was empty, I called arcpy.AddMessage() for each of the fields:               for f in arcpy.ListFields(layer):         arcpy.AddMessage('{},{}'.format(f.name, f.type))      ![enter image description here](http://i.stack.imgur.com/vyago.png)


title: How to correctly set the Independent Cities of Virginia in Open Street Map

text: I would like to use the data from Open Street Map to identify which county I am in, but I noticed something unusual with the Independent Cities of Virginia. Essentially, there are a few cities that are independent of counties, even if they are completely contained within a county. When I do a query (Via the Map Quest Geocoding API), it returns that I am in the surrounding county, and not the city. For instance, here is Augusta County, which contains two independent cities, Staunton and Waynesboro. I've seen this for a half dozen other independent cities, BTW. ![enter image description here](http://i.stack.imgur.com/G4Vni.png) It's a bit tricky to see, but Staunton has a county boundary approximating the city boundaries. When I call Staunton, it correctly shows Staunton City as the county. When I look at Waynesboro, it stays it is in Augusta, not Waynesboro City. What I'm trying to do is to figure out how to make the appropriate correction. I can see a couple of possibilities, including:   1. Somehow trace an interior boundary that matches the city of Waynesboro, and identify that county as Waynesboro City.   2. Somehow make the existing Waynesboro think it's a county. I'm not quite sure how either of those could be done, or which is the best, or even if there's a better solution all together. How can I make this correction?


title: How to Geocode use GeoPy API combined with GDAL Python bindings?

text: My question is related to the 'How to Geocode 300,000 addresses on the fly?' topic that has been asked by Mehul Nahar. I'm interested in running the processes explained by capooti in response to mehul question. I've done everything that was explained, but I was unsuccessful. Can you please let me know about the tools or components that I need to have on my computer to run the script? I have already ArcGIS 10.2, QGIS and GRASS installed. Is there any documentation related to the script?


title: Set Grass Region in Python grass.script

text: How can I set the Grass Region in Python without already having a map in the Mapset? Basically I am looking for the equivalent to this window ![enter image description here](http://i.stack.imgur.com/4MRKP.png) Currently I am only able to extent an existing region, which leads in my case to very big regions.


title: Bentley DGN format to Oracle Spatial

text: Hello did anyone ingested the DGN (all versions) to Oracle Spatial? -Rao.


title: how to add wms to QGis via cpp-plugin

text: I try to add a wms-Layer to QGIS Desktop from my plugin(c++)               QString uri= "crs=EPSG:28992&layers=natura2000:natura2000&styles=&format=image/png&url=http://geodata.nationaalgeoregister.nl/natura2000/ows?          QString myProviderName      = "wms";     QgsRasterLayer* lyr= mQGisIface->addRasterLayer(uri ,"test", myProviderName);      But in QGis i get the Message, the extent of the layer could not be calculated because of unknown dataprovider. What do i wrong? In QGIs Desktop manually adding this wms-layer is possible without problems. I also tested several other wms with same result. Thanks for all help.


title: How to get X and Y coordinate values of the map?

text: I want the `upper right` and `lower left` _coordinate_ of the exposed portion of the map, and for the whole map. What is the way to do this?


title: How to create cohesive Spatial Pixels from Spatial Points Dataset

text: I have a Spatial Point DF `spo` (covering an irregular shaped area of interest). The data are not on a regular grid due to crs transformation. My goal is a raster with predefined resolution and extent of the area of interest ( more spatial point data are to be mapped on this master raster) . Problems start when I               rasterize(spo, raster(ncol, nrow, extent, crs), spo$param)      I need to adjust `nrow`and `ncol` in a way so that I wont get moire patterns of NAs within my area of interest. I cant use a predefined (higher) resolution. As a solution to this, I thought I would need some kind of Spatial Pixel DF `spi`, that covers my whole area of interest (just like _meuse.grid_ in `library(raster); data(meuse.grid)`), and serves as a master grid. Then, I can use it to interpolate my data, e.g.               idw(param~1,spo,spi)      and by this, get full cover of my area of interest at my chosen resolution. But how can a SpatialPixelsDataFrame be produced from the point data? So in my view, the question boils down to: How to produce _meuse.grid_ from _meuse_ dataset? Maybe I'm taking the wrong approach here, so please let me know if more easily can achieved what I'm after, using a different way.


title: ArcGIS-Location "Location 1" in "Stops" is unlocated in Online Routing

text: I was trying to run a sample project "Routing" provided by ArcGis Runtime SDK. I tried the code but when the method `solve()` of `RouteTask` class is executed, it throws the following exception without returning any routes: > com.esri.core.io.EsriServiceException: Unable to complete operation. > Location "Location 1" in "Stops" is unlocated. Location "Location 2" in > "Stops" is unlocated. Need at least 2 valid stops. "Stops" does not contain > valid input for any route. I am not able to figure out what the actual problem is. Here is the link for the sample code and the code looks like this:                       map.setOnLongPressListener(new OnLongPressListener() {                  private static final long serialVersionUID = 1L;                  public boolean onLongPress(final float x, final float y) {                      // Clear the graphics and empty the directions list                 routeLayer.removeAll();                 hiddenSegmentsLayer.removeAll();                 curDirections = new ArrayList<String>();                 mResults = null;                      // retrieve the user clicked location                 final Point loc = map.toMapPoint(x, y);                      // Show that the route is calculating                 dialog = ProgressDialog.show(RoutingSample.this, "",                         "Calculating route...", true);                 // Spawn the request off in a new thread to keep UI responsive                 Thread t = new Thread() {                          @Override                     public void run() {                         try {                             // Start building up routing parameters                             RouteParameters rp = mRouteTask                                     .retrieveDefaultRouteTaskParameters();                             NAFeaturesAsFeature rfaf = new NAFeaturesAsFeature();                             // Convert point to EGS (decimal degrees)                             Point p = (Point) GeometryEngine.project(loc, wm,                                     egs);                             // Create the stop points (start at our location, go                             // to pressed location)                             StopGraphic point1 = new StopGraphic(mLocation);                             StopGraphic point2 = new StopGraphic(p);                             rfaf.setFeatures(new Graphic[] { point1, point2 });                             rfaf.setCompressedRequest(true);                             rp.setStops(rfaf);                             // Set the routing service output SR to our map                             // service's SR                             rp.setOutSpatialReference(wm);                                  // Solve the route and use the results to update UI                             // when received                             mResults = mRouteTask.solve(rp);                             mHandler.post(mUpdateResults);                         } catch (Exception e) {                             Log.w(ESRI_EXCEPTION_TAG, e);                             mException = e;                             mHandler.post(mUpdateResults);                         }                     }                 };                 // Start the operation                 t.start();                 return true;             }      I get exception on the following line: `mResults = mRouteTask.solve(rp);`


title: Error using con in Map Algebra in Python script

text: I am experiencing a painful problem in the middle of a Python script I am writing. When I get to a Map Algebra function, I have tried two different options and both return a different error. At first I wrote:               facteur_Pente = Con((surf_Drainee > 12)  &  (surf_Drainee <= 80)  &  (pente_Rclass <= 24),pente_Rclass + 1, Con((surf_Drainee > 80)  &  (pente_Rclass < 24),pente_Rclass + 2,Con((surf_Drainee > 80)  &  (pente_Rclass == 24),pente_Rclass + 1,pente_Rclass)))      The error I get from this is: " **TypeError: cannot concatenate 'str' and 'int' objects** " This seems normal, as I have strings (surf_Drainee, pente_Rclass) so I convert the values to str using the following:               facteur_Pente = Con(("surf_Drainee" > str(12))  &  ("surf_Drainee" <= str(80))  &  ("pente_Rclass" <= str(24)),"pente_Rclass" + str(1), Con(("surf_Drainee" > str(80))  &  ("pente_Rclass" < str(24)),"pente_Rclass" + str(2),Con(("surf_Drainee" > str(80))  &  ("pente_Rclass" == str(24)),"pente_Rclass" + str(1),"pente_Rclass")))      The problem is that it doesn't work, and produces the following error: " **RuntimeError: ERROR 000732: Input Raster: Dataset pente_Rclass1 does not exist or is not supported.** " I don't even have a pente_Rclass1 in my script! Does anybody know how I could fix this problem in order to continue my script? Regards,


title: Activate layers that are in a checkbox list

text: I'm trying to do a code to activate layers that are in a checkbox list, the code that I'm trying to use is this:                   MxDoc = My.ArcMap.Document         Maps = MxDoc.Maps                   If chxActivated.SelectedItem = True Then             GxLayer.Layer.Visible = True [COLOR="#FF0000"]HERE IS THE PROBLEM[/COLOR]             MxDoc.UpdateContents()             MxDoc.ActiveView.Refresh()          End If      I add the layers with this code                   MxDoc = My.ArcMap.Document         Map = MxDoc.FocusMap         'new layer         GxLayer = New GxLayer         GxFile = GxLayer         'Layer file path         GxFile.Path = "C:\Layers\Aditional Information.lyr"         'Add the new layer to ArcMap         Map.AddLayer(GxLayer.Layer)         'Add the layer invisible         GxLayer.Layer.Visible = False         'Refresh the view         MxDoc.ActiveView.Refresh()     End Sub      Thank you!


title: Merging multiple polygon shapefiles

text: I want to merge two polygon shapefiles in one without losing the attributes using ArcMap 10.1. To be exact I am creating a geological map and have some geologic formations in one shapefile and some other formations in another shapefile. Thank you for any answer!


title: ArcMap Model Builder Iterator to create selections

text: This is a sequel to my question about aggregating polygons in order to create neighborhoods. I want to automate the geoprocessing of the above task. When I use Model Builder's "Iterate Feature Selection" iterator I get the below error: ![enter image description here](http://i.stack.imgur.com/CRefK.png) I originally tried to link the "Iterate Feature Selection" iterator directly to the Aggregate Polygon tool. But as per this thread I inserted the Make Feature layer tool in between. (Now removed, see comments below) ![enter image description here](http://i.stack.imgur.com/rLzPo.png) My problem is I'm unsure how to link the iterator's selection to the input of the Aggregate Polygon tool. Really, **why does the GROUP BY function of "iterate feature selection" iterator not create a complete SQL condition statement??** The following threads have proved un-useful as well: Model builder question regarding iterations (newbe : ( ..... ) Model builder: select feature based on attribute, iterator? thanks in advance


title: How can I retrieve the features in a layer without adding it to a map?

text: I am using the ArcGIS API for JavaScript. My layer of `Point` features is exposed through a REST service. I want to retrieve all of these features when the application starts up, then show different subsets of features at different times, depending on user input (such as text search) and client-side calculations (such as clustering calculations). If I just create a `FeatureLayer` normally, it won't actually retrieve any features before it's added to a map -- its `graphics` property is just `[]` when I look at it in Firebug. I tried calling `FeatureLayer.refresh`, but it just complains that `map` is undefined. How do I retrieve the features without having to attach them to a map?


title: Problem with merging and dissolving a shapefile in R

text: I am migrating from ArcMAP to opensource GIS. I've succesfully installed R and many packages in a Ubuntu 12.04 machine. I am used to linux and programming in C, but not to R and its syntax. I am trying to "dissolve" (fuse, merge) the spatial polygons of a shapefile by a certain code that is loaded from a csv file (e.g., to join counties' polygons by its FIPS State code to obtain a single polygon(s) for each State). I began composing my script with pieces from different places, but mostly this questions and its answers (that I will upvote as soon as I get enough reputation): http://stackoverflow.com/questions/19791210/r-ggplot2-merge-with-shapefile- and-csv-data-to-fill-polygons http://stackoverflow.com/questions/19613004/what-is-the-proj4-string- associated-with-my-shapefile The following code is the initial, incorrect version of the script, that has been improved thanks to the answer by SlowLearner. I leave it here (for now) until the actrual question and the answers get its final form.               # In this script I'm trying to perform the next steps:     # - Load a shapefile that contains some base polygons ("zones").     # - Load a table that contains a possible grouping of the zones into "regions".     # - Merge ("join" in ArcMap) the zones' shapefile with the grouping table     #   (common code "PROVMUN").     # - Create the regions' polygons by merging the zone polygons with same region     #   identifier ("CODER").     # - Merge the newly created regions' polygons with another table that contains     #   some region-level statistics.     # - Plot into a PDF a map with regions filled with different colors, borders     #   of zones in grey, and borders of regions in black          library(gpclib)     library(rgdal)     library(maptools)     library(ggplot2)     gpclibPermit()     #args <- commandArgs(TRUE)     ## Input parameters (fixed in this example)     work.dir <- "."     zones.shp.name <- "SP0111"     zones.shp.id <- "PROVMUN"     zones.csv.name <- "statsz.csv"     zones.csv.id <- "PROVMUN"     zones.csv.region <- "CODER"     regions.csv.name <- "statsr.csv"     regions.csv.id <- "CODER"     regions.shp.name <- "SP01_R_diss"     plot.pdf.name <- "SP01_R_plot.pdf"          # Read base polygon (zones) shapefile     zones.sp <- readOGR(dsn = work.dir, layer = zones.shp.name)          # Read regionalisation and zone-level statistics     zones.data <- read.csv(zones.csv.name)     # and merge (join) with the base polygons     zones.spwd <- merge(x=zones.sp, y=zones.data, by.x=zones.shp.id, by.y=zones.csv.id)          # Create region polygons by dissolving zone polygons     zones.IDs <- as.numeric(zones.spwd[[zones.csv.region]])     regions.sp <- unionSpatialPolygons(SpP = zones.spwd, IDs = zones.IDs, avoidGEOS = TRUE)     # (add threshold parameter to remove slivers on tricky shapefiles)          # Now regions.sp is an SpatialPolygons object, without data slot, so it seems     # the codes of regions are lost. But shooting in the dark I did this:     names(regions.sp)     # and I get a list of the codes I need (BTW: why this object has no dimension?)          # I need to (re?)insert the codes of the regions so that it allows to perform another     # merge (see below in the code)          # For that, I found this answer:     # http://stackoverflow.com/questions/15201875/merge-polygons-and-plot-using-spplot     # But here I'm shooting in the dark, as I don't really understand the linked example     # and could not make it work (no dimensions where there should be, etc.).          # So I try this (from the examples about plotting with ggplot2):     regions.spdf <- fortify(regions.sp)     # And the codes are now in the column named "id", so I use that          # Read region-level statistics (to merge with the dissolved regions' polygons)     regions.data <- read.csv(regions.csv.name)     # Merge that data with the regions' SpatialPolygons object     regions.spdf.wd <- merge(x=regions.spdf, y=regions.data, by.x="id", by.y=regions.csv.id)          # Write to shapefile the region polygons (with the stats already atached?)     writeOGR(obj = regions.spdf.wd, dsn = work.dir, layer = regions.shp.name, driver="ESRI Shapefile")      The writeOGR() command fails with error message: > Error: inherits(obj, "Spatial") is not TRUE               # Finally, the map plotting:     # Make PDF document     pdf(plot.pdf.name)     # Make empty map plot     mapplot <- ggplot()     # Add regions' polygons filled by region id (unique color for each region)     mapplot + geom_polygon(data = regions.spdf.wd, aes(fill = "id"))     # Add zone borders in grey     mapplot + geom_path(data = zs.df, aes(color = "grey40"), size=0.4)     # Add region borders in black     mapplot + geom_path(data = rs.df, aes(color = "black"), size=0.8)     # Add labels and title     mapplot + labs(x = "Longitude", y = "Latitude", fill = "FR") + opts(title = "Test map")     dev.off()      If I disregard the error in writeOGR and go through the plotting, I get several error like this: > Error in as.environment(where) : where' is lost" I have uploaded the shapefile and reduced version of the tables to Google Drive: https://drive.google.com/file/d/0B-TKzP- _yJNMdDhlbi1HbjFqNDg/edit?usp=sharing Any help will be greatly appreciated. Finally, this is my first question in stackexchange (any of them), don't hesitate to comment/criticize anything, editing title and whatnot. * * * EDIT: After trying SlowLearner's answer I have improved the script and hopefully the replicability of the problem: it now uses rgeos and a free shapefile from http://www.gadm.org/country (country Spain, complresed file ESP_adm.zip, shapefile ESP_adm4.*). However, it still fails to perform the desired dissolve operation. This is the content of file grouping.csv necessary for the dissolve operation:               MUN_NAME,MA_NAME     Alicante,ALICANTE     El Campello,ALICANTE     San Vicente del Raspeig,ALICANTE     Sant Joan d'Alacant,ALICANTE     Busot,ALICANTE     Mutxamel,ALICANTE     Agost,ALICANTE     Crevillent,ELCHE     Elche,ELCHE     Santa Pola,ELCHE     Benidorm,BENIDORM     La Nucia,BENIDORM     Villajoyosa,BENIDORM     L'Alfàs del Pi,BENIDORM     Polop,BENIDORM     Finestrat,BENIDORM     Callosa d'En Sarrià,BENIDORM      Following is the current state of the script (with most of the intended plotting removed, for now):               library(rgdal)     library(rgeos)     library(ggplot2)          # ESP_adm4 contains the municipalities of Spain (polygons)     # You can download the same shapefile from http://www.gadm.org/country     mun.spdf <- readOGR(dsn = ".", layer = "ESP_adm4")          # From this answer I learned to subset my SpatialPolygonsDataFrame easily:     # http://stackoverflow.com/a/13443741     mun.spdf <- mun.spdf[mun.spdf$ID_1 == 944,c("ID_1","NAME_1","ID_2","NAME_2","ID_4","NAME_4")]     # Now z.spdf contain only municipalities in Comunidad Valenciana (COM==944)     #head(mun.spdf@data)     #mun.spdf@data[mun.spdf$ID_1 != 944]          # grouping.csv contains a grouping of some of the municipalities into three regions     grouping.df <- read.csv("grouping.csv")     # For the dissolve operation, we need to add it to the polygons     # In ArcGIS we would use Merge (Data Management), a.k.a. "Join" in the interactive menu     # http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000055000000     # R scripting is beautiful, once you understand it:     mun.spdf.wd <- mun.spdf     mun.spdf.wd@data <- merge(x = mun.spdf@data, y = grouping.df, by.x = "NAME_4", by.y = "MUN_NAME", all.x = TRUE, all.y = FALSE)          # Now we should be ready to perform the dissolve     # In ArcGIS 10 we would use Dissolve (Data Management)     # http://resources.arcgis.com/en/help/main/10.1/index.html#//00170000005n000000     # and set the dissolve_field to "AM_NAME", and it would generate a new shapefile with     # where each polygon would be identified by the corresponding "CODER" value          # Here in R:     regions.ids <- mun.spdf.wd[["MA_NAME"]]     # z.IDs contains the values of the column MA_NAME, where rows with AM_NAME != NA belong to a region     # and polygons with the same MA_NAME should be fused together (there can be many distinct regions,     # in the example data there are three: "ALICANTE", "BENIDORM", "ELCHE")     # If I understand rgeos documentation, the next command should dissolve all of the polygons with     # same value in MA_NAME into one polygon, discarding the ones with AM_NAME == NA     regions.sp <- gUnionCascaded(mun.spdf.wd, id = regions.ids)          # http://stackoverflow.com/a/13666703 : this answer from SlowLearner was very useful to allow me     # plot with x11(type='Xlib'), that seems to be the default.     set_Polypath(FALSE)     plot(regions.sp)      The results make no sense: there should appear three contiguous clusters, instead of this spread. However, we can check that there are three regions in the new shapefile with the correct names: names(regions.sp) Maybe the merge of the loaded CSV file with the shapefile is not performed correctly? Lets try now another dissolve to obtain the borders of the provinces (NAME_2), using the original shapefile               prov.ids <- mun.spdf[["ID_2"]]     provinces.sp <- gUnionCascaded(mun.spdf, id = prov.ids)     # This should plot the borders of the three Spanish provinces Alicante, Castellón y Valencia:     par(new = TRUE)     plot(provinces.sp)      And again this fails, as my understanding of what is wrong here. So this happens with a shapefile downloaded from gadm as well as with the shapefile initially uploaded with this question. It makes me think we can rule out a problem in the shapefile. Therefore, what is wrong in this code? why the union of polygons does not produce the expected provinces or metropolitan areas? Any ideas? Thank you very much. Edit: I have finally figured out why this was not working: you can not merge only the data slot because it reorders the rows of that slot while leaving untouched the polygons slot, therefore, the union of polygons is performed over corrupted indices. The solution is to merge the whole object, R is really awesome and can handle it, so the code from @SlowLearner should look like this:               #zones.spwd@data <- merge(zones.sp@data, zones.data, by = "PROVMUN") ## <- Wrong, it corrupts the corresponce between slots @polygons and @data     zones.spwd <- merge(zones.sp, zones.data, by = "PROVMUN") ## This works! the changes in parameter all are


title: How to make traffic time contour map in ArcMap?

text: I want to make travel time contour map which shows the area people can go from an origin (e.g. downtown, airport) within certain time interval. I have a roadway linear shapefile with the information of travel time and some example maps found online, but lack of the idea how to make the map. Thanks for your help! ![enter image description here](http://i.stack.imgur.com/7998r.jpg) ![enter image description here](http://i.stack.imgur.com/wluBL.jpg)


title: WmsGetFeatureInfo not working in larger scales

text: I use Openlayers/Geoserver on ubuntu server. Geoserver runs at tomcat6. I have grouplayer of several layers all shown at certain scales. When I click at feature at small scale, getfeatureinfo works. When I zoom in and then click at same feature, response from server is empty. No errors at geoserver log or apache log. Also when I send that response manualy via web browser, response is empty, too. At one layer, there are even two zoom levels working, only the last one is empty. BBox and pixels in request look OK. ![](http://i.imgur.com/YMYpLMF.png) >>> zoom >>> ![](http://i.imgur.com/tWsMZXn.png) Most probably, the problem isn't in code, but here it is:               info = new OpenLayers.Control.WMSGetFeatureInfo({             url: 'http://158.196.109.46:8080/geoserver/diplomka/wms',              title: 'Identify features by clicking',             queryVisible: true,             maxFeatures: 1,             eventListeners: {                 getfeatureinfo: function(event) {                     map.addPopup(new OpenLayers.Popup.FramedCloud(                         "chicken",                         map.getLonLatFromPixel(event.xy),                         null,                         event.text,                         null,                         true                     ));                 }             }         });         map.addControl(info);         info.activate();      Proxy.cgi is also set properly I think. I also tried raising tomcat heap memory, it didn't help. So I really don't know, what should be causing this. Anybody experienced this, or knows what to do? I would really appreciated that.


title: How to extract Z data from a point layer

text: I have created a point layer with ~150 points. I have a DEM beneath it, which is where I have taken the X and Y values from. Does anybody know how I can retrieve the Z values, without using the info button! Using ArcGIS 10. I have licences for the spatial analyst etc tools.


title: How to get QGIS to play with pandas (Windows)?

text: has anybody been successful getting pandas to work with QGIS's Python installation, on Windows? I've used the OSGeo4W installer to install QGIS2.2, which comes with Python 2.7 and Numpy 1.5.1. Any reasonably new version of pandas requires Numpy 1.6.1 or newer. So that seems like my main problem. I've tried: \- Building an older version of pandas (e.g. 0.7.3 and 0.4.1) using easy_install, but that gives me a very long list of errors, at least the last of which are related to distutils\msvc9compiler.py \- Using binaries of old pandas versions, but they are all built for Numpy 1.6.1 \- Upgrading Numpy, which also gives me an error message. In case anybody is interested, the message ends with "File “C:\OSGeo4W\apps\Python27\lib\distutils\msvc9compiler.py”, line 299, in query_vcvarsall” newline “raise ValueERror(str(list(result.keys())))” newline “ValueError: [u’path’]" I'm wondering if anybody else has ever been successful in getting pandas to work with a QGIS install on Windows?? It would be so easy if QGIS would ship with, say, a version of Numpy that's less than three years old. Does anybody know if such an upgrade is in the works for upcoming QGIS versions?


title: Does an open source alternative to ESRI's Image Server exist?

text: I noticed the default way to store imagery is by using a database in the open source world. I like ESRI's Image Server software because it removes the need for a database and makes backing up a lot easier since its just files. Is there an Image Server type software in the open source world?


title: What User Testing techniques do you use?

text: In producing dynamic maps for students/users I have found the time invested in user testing well worth it. I use a varient of think aloud protocol called hallway testing. I grab any 2 or 3 likely users for 20 minutes each and set them a sample task (its important to get them to **do** something rather than just comment). They explain what they are doing when they get stuck or find something confusing and I prompt them if they aren't saying anything. Its a great way to discover consistent issues users have with your maps, most major problems will pop up with 3 users. I use software to record the screen (silverback) but just taking notes was pretty effective. With a big project I also advocate 'test early, test often'. I've found clients to be unhappy about funding user testing but its so quick and low tech that I tend to do it anyway to make sure there aren't unknown issues with a project Do you use user testing when producing dynamic maps?   If so, what techniques have you found to be effective and why?


title: Batch processing within geocommons/geocoder

text: I've just installed geocommons/geocoder from github. The only example they provide is how to geocode a single address using interactive ruby. I am looking for a way to batch-geocode thousands of addresses. Is there a way to run this geocoder in batch mode, using an input file of addresses? Thank you! natalia


title: Arcmap: a different delimiter when symbolizing "Unique values, many"?

text: In Arcmap 10.1, is it possible to use something other than a comma to separate the value fields when symbolizing by "Unique values, many"?


title: Monster Spatial database server suggestions

text: I'm going to setup a open Mapserver w/ Openlayers, Geojango, Postgis using Postgresql. And I need some suggestions from people that have been shopping for what I need longer than I have. What would be the best configuration for such a db server and from where? Im a big mac fan and was thinking of a Xserve w/ 42g ram and two quad cores (2.66). What does postgresql like best as far as server software/hardware? Can it fully utilize multi cores? Im on somewhat of a tight budget so price is a big factor. But what would be the dream machine and its cheaper inequivalent? sockets; i'd probably have 20 to 30 ppl at the lease almost 100 at most, maybe more. Serving up more than 100 layers with some common canned geoprocessing type stuff going on.


title: clicking button in QueryAttributesTask by code

text: I need help, I have a website developed using Argis Server 9.3 and Visual Studio 2008. The map is inside the website via an iframe and it receives the input to be used in a custom query via javascript, this is working fine, except that, the users still have to clic the button of the QueryAttributesTask component to execute the query task. Does anyone knows how can I call the button clic code from another button (like an asp button)? Any ideas on how to accomplish this?


title: modeling of spatial objects

text: I have a line and polygon in this image. If i find the intersection of line exterior and polygon boundry, what will be the graphical representation of both objects intersection.. its my concept.. ![http://imm.io/2xTh](http://i.stack.imgur.com/LBUJw.png)


title: How to get spatial reference of a shapefile using python?

text: I always use arcpy.GetParameterAsText to get the spatial reference from the user. There should be a way to decide what is the spatial reference of a shapefile but I can't find it. Something like Describe....


title: Abstract Data Store

text: I'm trying to create a new extension for data store in geoserver that can take data from in-memory arrays.But as I'm kind of new to geotools I'm not sure how to do it.Any help would be really appreciated.Thanking you in advance