title: What is the correct way to pass a polygon to a geoprocessing service successfully

text: We have a geoprocessing service that has been in limbo for 6 months. Finally after a great decision to eliminate one element of the service, we got it working this week. I am now at the point where I can 'kill' the service call by just passing the polygons like a light bulb turned on and off. So it seems this is focused now on how I pass the polygon to the service. I've looked at it until I'm blue in the face and cannot see anything wrong with it. The symptoms are as follows:                 The server service receives the call.  It gets the bulk data just fine.       The server service locates the class to call and creates an instance.       Parameter info is called and it starts to process the parameters for calling Execute       It returns Fail.      I figured all of that out by putting the server logging in Debug mode. So now the bombardment of code. This takes the origin point, targets, and blockages (the polygons) and preps them for calling the service.               Point origin = (Point)GeometryOperations.Project(Point.CreateFromLatitudeLongitude(Convert.ToDouble(SelectedStation["lat"]), Convert.ToDouble(SelectedStation["long"])), ESRI.ArcGISExplorer.Geometry.CoordinateSystem.ProjectedCoordinateSystems.World.WGS1984WebMercatorAuxiliarySphere);     List<Point> targets = chosen.Select(location => (Point)GeometryOperations.Project(Point.CreateFromLatitudeLongitude(location.Latitude, location.Longitude), ESRI.ArcGISExplorer.Geometry.CoordinateSystem.ProjectedCoordinateSystems.World.WGS1984WebMercatorAuxiliarySphere)).ToList();     BlockagesAndDelays bandd = Functions.GetBlockagesAndDelays();     List<Polygon> blockages = new List<Polygon>();     foreach (Polygon item in bandd.Blockages)     {         blockages.Add((Polygon)GeometryOperations.Project(item, ESRI.ArcGISExplorer.Geometry.CoordinateSystem.ProjectedCoordinateSystems.World.WGS1984WebMercatorAuxiliarySphere));     }          helper.UpdateMessage("Calling service");     double[] results = Geoprocessing.GetLaunchPointDistanceToTargets(helper, origin, targets.ToArray(), blockages.ToArray());      The method called now converts the raw data into the format to send to the server and makes the call. Here is where I may be screwing up somehow.               using (CalculateDistances_GPServer service = new CalculateDistances_GPServer())     {     try     {         List<double> results = new List<double>();         progressDialog.UpdateMessage("Establishing parameters");              service.Url = Functions.GetGPServiceEndpoint("CalculateDistances");         CalculateDistances.GPToolInfo toolinfo = service.GetToolInfo("Calculate Distances");         CalculateDistances.GPParameterInfo[] parameters = toolinfo.ParameterInfo;         CalculateDistances.GPFeatureRecordSetLayer firstParameterLayer =             (CalculateDistances.GPFeatureRecordSetLayer) parameters[0].Value;         CalculateDistances.GPFeatureRecordSetLayer secondParameterLayer =             (CalculateDistances.GPFeatureRecordSetLayer) parameters[1].Value;         CalculateDistances.GPFeatureRecordSetLayer thirdParameterLayer =             (CalculateDistances.GPFeatureRecordSetLayer)parameters[2].Value;              CalculateDistances.PointN newPoint =             SoapConverter.GeometryToSoap<Point, CalculateDistances.PointN>(launchPoint);         CalculateDistances.Record[] firstParameter = new CalculateDistances.Record[1];         object[] recordValues = new object[2];         recordValues[0] = 1;                // OBJECTID         recordValues[1] = newPoint;         // SHAPE         firstParameter[0] = new CalculateDistances.Record {Values = recordValues};         firstParameterLayer.RecordSet.Records = firstParameter;              CalculateDistances.Record[] secondParameter = new CalculateDistances.Record[targets.Length];         for (int index = 0; index < targets.Length; index++)         {             CalculateDistances.PointN targetPoint =                 SoapConverter.GeometryToSoap<Point, CalculateDistances.PointN>(targets[index]);             recordValues = new object[2];             recordValues[0] = index+1;        // OBJECTID             recordValues[1] = targetPoint;    // SHAPE             secondParameter[index] = new CalculateDistances.Record {Values = recordValues};         }         secondParameterLayer.RecordSet.Records = secondParameter;              CalculateDistances.Record[] thirdParameter = new CalculateDistances.Record[blockages.Length];         for (int index = 0; index < blockages.Length; index++)         {             recordValues = new object[2];             CalculateDistances.PolygonN blockageArea =                 SoapConverter.GeometryToSoap<ESRI.ArcGISExplorer.Geometry.Polygon, CalculateDistances.PolygonN>(blockages[index]);             recordValues[0] = index + 1;        // OBJECTID             recordValues[1] = blockageArea;     // SHAPE             thirdParameter[index] = new CalculateDistances.Record { Values = recordValues };         }         thirdParameterLayer.RecordSet.Records = thirdParameter;              CalculateDistances.GPValue[] values = {firstParameterLayer, secondParameterLayer, thirdParameterLayer};         progressDialog.UpdateMessage("Calling service");         string jobId = service.SubmitJob(toolinfo.Name, values, null, null);         CalculateDistances.esriJobStatus jobStatus = service.GetJobStatus(jobId);      If anyone can see what is messed up I'd greatly appreciate it.


title: Is there a difference between projecting shapefiles then rasterizing & rasterizing shapefiles then projecting?

text: I currently have a project where I rasterize shapefiles into simple presence/absence (1 & 0) rasters and then add each one to a running total to create a total richness raster at the end. The shapefiles are in WGS84 projection, but I need the final raster to be in Behrmann Equal Area projection. What I want to know is if there is any difference in the final raster if I project all the shapefiles to Behrmann and then rasterize and merge each of them compared to rasterizing and merging, then projecting the final raster to Behrmann at the end? Kind regards,


title: How can I use py2exe with the OSGeo4W Python installation?

text: I have a default Python installation in c:\Python27. I have also installed Python through the OSGeo4W installer, which has created a "sandboxed" Python installation with all the geospatial bindings I need. Now, I need to package a script using py2exe, but I can't figure out how to install or configure py2exe to work with the OSGeo4W installation. I have configured OSGeo4W to install external packages with pip, but py2exe doesn't seem to install with pip. If I try to use the downloaded installer, it only finds the system install in the registry. Any ideas how to get it working?


title: How to query multiple features from Geoserver WFS?

text: Are there any methods to query several features by id other than one entry by a single url request? The scenario is I store around 200 point layers in PostgreGIS. Each maintains 100k entries GPS points. My boss wants to let them "move" in the browser. Obviously it's impossible to query all entries at once. So I want to query like 1k entries for one layer in every query. So basically I want to query the features by their id, like "starID<=ID<=endID". But it seems like by WFS url query, it can only query one entry every time by the featureID=ID. I'm also thinking about using the ordinary Postgre other then PostgreGIS, but just in case there might be spatial queries in the future.


title: Why is my shapefile not visible in Tilemill?

text: I have 3 layers of vectors (shapefiles) that I've loaded into Tilemill. However, nothing is visible in Tilemill. **What am I doing wrong, and what should I do to get my map to show up?** This is part of a web app to show the shop floor of a factory. I've built 3 layers in Inkscape, and used QGIS to get it into the shapefile format. I initially had some errors getting Tilemill to load the files, but that now seems to be resolved. You can see the problem here: ![Tilemill with 2 layers loaded](http://i.stack.imgur.com/Prz0s.png) I've looked at all zoom levels, and there's nothing there. When I view a layer's features, though, it seems that there's data: ![Tilemill with layer's features menu visible](http://i.stack.imgur.com/1nQQK.png) Thanks!


title: How to use -where SQL in an ogr2ogr loop in GDAL/OGR (bash)?

text: I'm using the following code:               while read line; do     for att in $line; do     ogr2ogr -where "newcode='$att'" -f "ESRI Shapefile" $OUTPUT/$att.shp $INDATA     done     done < $SCRIPT/newcodes.txt      The text file is being read correctly as shapefiles named after all names in the textfile list are being created. However, the shapefiles are empty (no attribute table or polygons). The ogr2ogr line of code works outside of the loop on a single feature within the shapefile. The newcodes.txt contains 3-5 letter codes e.g. aaa, gbe. I think the problem is getting bash to treat $att as text to search for in the attribute table of the shapefile within the `-where` switch. Any help much appreciated. **Update:** ![Attribute Table Sample](http://i.stack.imgur.com/0KpzK.png) The shapefile has 15 columns in the table, and 6506 rows/features. The 'newcode' field which I'am trying to read with the `-where` switch was joined to the table. I have tried using the loop with a pre-existing field and it still produced empty shapefiles. **Update 2:** After a bit more interrogation, I've noticed that the shapefiles are being produced with a space after their name. E.g. 'aaa .shp' rather than 'aaa.shp'. This probably means that the `$att` variable is being read with a space also, and therefore not matching anything in the attribute table. I have checked the text file, and there are no spaces in there. Still unsure how to fix.   ![enter image description here](http://i.stack.imgur.com/LkACK.png) e.g. cdg cdh cdi cej cea eig fjj fja fjf


title: How to Add a Legend in config.yaml File in GeoServer Printing

text: I am unable to add a legend in config.yaml file in GeoServer. Kindly, could any one share the code to show the legend in GeoServer print?


title: Mixed geometry types styling with Geoserver 2.2

text: I have a PostGIS table with geometry colunm. In this column i have geometries of several types (polugons, points, lines). I publish this this table on Geoserver. Now i want to create OpenLayers layer with this data.   But how to make a style for this data?   Im already read http://docs.geoserver.org/stable/en/user/styling/sld- tipstricks/mixed-geometries.html . And dont want to create `gtype` in postGIS if this possible.   I try to create new style on geoserver but i think that do something wrong               <?xml version="1.0" encoding="ISO-8859-1"?>     <StyledLayerDescriptor version="1.0.0"      xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"      xmlns="http://www.opengis.net/sld"      xmlns:ogc="http://www.opengis.net/ogc"     xmlns:xlink="http://www.w3.org/1999/xlink"     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">     <!-- a Named Layer is the basic building block of an SLD document -->              <Rule>          <ogc:PropertyIsEqualTo>           <ogc:Function name="geometryType">            <ogc:PropertyName>the_geom</ogc:PropertyName>           </ogc:Function>           <ogc:Literal>Polygon</ogc:Literal>          </ogc:PropertyIsEqualTo>          <PolygonSymbolizer>            <Stroke>               <CssParameter name="stroke">#0000ff</CssParameter>               <CssParameter name="stroke-width">1</CssParameter>             </Stroke>          </PolygonSymbolizer>         </Rule>      And get a error:               line 10: cvc-complex-type.2.4.a: Invalid content was found starting with element 'Rule'. One of '{"http://www.opengis.net/sld":Name, "http://www.opengis.net/sld":Title, "http://www.opengis.net/sld":Abstract, "http://www.opengis.net/sld":NamedLayer, "http://www.opengis.net/sld":UserLayer}' is expected.     line 11: cvc-complex-type.2.4.a: Invalid content was found starting with element 'ogc:PropertyIsEqualTo'. One of '{"http://www.opengis.net/sld":Name, "http://www.opengis.net/sld":Title, "http://www.opengis.net/sld":Abstract, "http://www.opengis.net/sld":LegendGraphic, "http://www.opengis.net/ogc":Filter, "http://www.opengis.net/sld":ElseFilter, "http://www.opengis.net/sld":MinScaleDenominator, "http://www.opengis.net/sld":MaxScaleDenominator, "http://www.opengis.net/sld":Symbolizer}' is expected.      What i doing wrong? **UPDATE** I try to show a points and make this rule for this geom type:                <Rule>             <ogc:Filter>                         <ogc:PropertyIsEqualTo>                     <ogc:Function name="geometryType">                         <ogc:PropertyName>the_geom</ogc:PropertyName>                     </ogc:Function>                     <ogc:Literal>Point</ogc:Literal>                 </ogc:PropertyIsEqualTo>             </ogc:Filter>                 <PolygonSymbolizer>             <Stroke>                 <CssParameter name="stroke">#0000ff</CssParameter>                 <CssParameter name="stroke-width">2</CssParameter>             </Stroke>             </PolygonSymbolizer>         </Rule>      But not see then on my map. Whats wrong?


How can I tell mapserver to ignore a filter that is set in a mapfile? I am displaying a wms layer with point features using mapserver. In one of the layers in my mapfile I have defined a filter similar to this:               FILTER          "mycolumn='%myvalue%'"      such that I can use runtime substitution in a mapserver request similar to this one:               http://someurl.com/mapserv?map=mymap.map&layer=mylayer&mode=map&myvalue=thisvalue      This works just fine. However, sometimes I would like to display the entire dataset and thus forcing mapserver to ignore the filter. Is it possible to deactivate the filter in a request?


title: How to display rental/property data from Zillow on ArcGIS Online?

text: I have a requirement wherein I need to plot the Zillow data on ArcGIS Online. I have looked at Zillow, and it does provide an API, but I am not sure how to integrate the API data in the ArcGIS Online template. The only head start I have is this Zillow Widget, but this is based on Flex so I am not sure how to start. I do know that I can create web services and plot the data, but the problem is, whenever I try to plot the data, the only thing that gets plotted is the State, i.e a point on the state. How do I plot the rental/property Zillow data on ArcGIS Online?


title: how to permanently add raster to satellite image in qgis

text: As mentioned in the title i need to integrate a raster to a satellite image but have no idea as to how to do this. I tried the add grid function and it works nice but is not permanent. I insert it via map composer but as this inhibits zooming I need to save it on the image. The person who will use the image as a map wishes to be able to zoom just as in any image. How is it possibly to permanently add the a raster with 10 square kilometre grid size permanently to a satellite image I have as TIFF format?


title: I'm Having trouble with the Arcpy geometry.disjoint function

text: I have points along along a line. Points that intersect the segment of the line to return attribute information. My line as three connected segments. ESRI Docs says > disjoint (second_geometry) Indicates if the base and comparison geometries > share no points in common. Two geometries intersect if disjoint returns > False.               with arcpy.da.UpdateCursor(in_point_along_line, ["SHAPE@", "Azimuth", "OID@"]) as cursor:             for row in cursor:                 azimuth = get_azimuth_intersect(row[0], row[2])                 row[1] = azimuth         def get_azimuth_intersect(in_geom, id):             for geom in in_geom_azimuth_list:                 if in_geom.disjoint(geom[0]) is False:                     print id                     print in_geom.centroid                     print "Intersect"                     print geom[0].firstPoint                     print geom[0].lastPoint                     print geom[0].length      Result:               This one is reporting False, and is correct the point does intersect the line segment.     252501.7229 5653692.1437 NaN NaN The Point     252657.6907 5653761.9575 NaN NaN The line firstpoint     251107.9884 5653068.2849 NaN NaN The line lastpoint          The one is reporting it's also False? But it does not intersect...     252501.7229 5653692.1437 NaN NaN The Point     251107.9884 5653068.2849 NaN NaN The line firstpoint     251169.5741 5652467.8243 NaN NaN The line lastpoint          The one is reporting it's also False? But it does not intersect...     252501.7229 5653692.1437 NaN NaN The Point     251169.5741 5652467.8243 NaN NaN The line firstpoint     250415.1493 5651836.5709 NaN NaN The line lastpoint      What I'm I missing?


title: Problems searching for attributes

text: I've come across this problem a few times and have not been able to find a satisfactory answer. I am new to python though so perhaps I have not been looking in the right places! Anyway, In ArcGIS I have been using the following code to select entries by attribute in the attribute table:               import arcpy          arcpy.env.workspace = "C:\Users\T0015685\Documents\WCGIS\Geog485\Lesson1\us_hydro.shp"          qry = '"COUNTRY" = \'CANADA\''     flayer = arcpy.MakeFeatureLayer_management("us_hydro", "us_hydro_Hudson6")     arcpy.SelectLayerByAttribute_management(flayer, "NEW_SELECTION", qry)     cnt = arcpy.GetCount_management(flayer)          print "The number of selected records is: " + str(cnt)      This script works fine. HOWEVER, when I alter the query entries from verbal to numeric attributes I get an error. So, for example, if I change the query to - qry = '"HYDRO_= \'1744\'' - the script fails. I think this is something to do with the SQL expression but I'm not sure. I want to be able to use code such as this to search feature ID's but cant because every time I use numbers the code doesn't work. Could someone please explain why this is the case or perhaps offer a solution.


title: Geoportal resource URL to Geoserver does not show Preview link in search result

text: I have a Geoportal catalog entry with a resource URL like: > > http://data.ilri.org:8080/geoserver/nurc/Arc_Sample/wms?request=getcapabilities However the preview link is not shown in the search result. I thought that maybe this was related to the **_resourceLinkBuilder.preview.filter_** variable in gpt.xml but even if I set the filter to empty (No filters at all)               <parameter key="resourceLinkBuilder.preview.filter" value=""/>      The preview link still does not show. Any ideas how to make it work? Please note that if I call the Geoportal Preview Window with my URL the Preview works fine. > > http://data.ilri.org/geoportal/catalog/livedata/preview.page?url=http://data.ilri.org:8080/geoserver/nurc/Arc_Sample/wms?request=getcapabilities Many thanks, Carlos


title: QGIS contours smoothing and cleaning

text: Trying to utilize the contour extraction QGIS tool but having trouble creating smoothed contours. The smooth, generalize, and dissolve tools have not worked for me. Also the plugin for contours appears to no longer be working? ![screen shot](http://i.stack.imgur.com/QcHrC.png)


title: How to work with Mobile Project Center tool in Arcgis?

text: ![Herramienta de Arcgis: Mobile Project Center ](http://i.stack.imgur.com/xR5aD.png) How is this tool structured and how to work with it? Sincerely, I can't find detailed information on this.


title: How to filter data by category, then map only maximum-value row for each?

text: I'm building a consolidated map of all transit systems in my area in QGIS. I have a GTFS dataset of bus lines, and it includes shapes for every trip on every route in the system (as it should), but that's far more detail than I need, and very slow to render. I'd like to filter the data to map only the longest possible trip on each route in each direction (inbound/outbound). What's the best way to do this? Line number, route direction, and shape length are all separate columns. I'd prefer to keep the underlying dataset intact if possible.


title: Extract an ArcSDE datasets (rasters/shape files) using python?

text: I would like to write a python script to extract/export a dataset (rasters/shape files) in ArcSDE. I saw this posts to connect to the server: Define Workspace for SDE Connection in Python But how can I extract the a dataset? Many thanks. Carlos.


title: How to read/convert arc000dat files?

text: I have a number of files that I am trying to import, but I'm just failing at the first hurdle. I've got a directory structure and file list that looks like this:                Directory of C:\blah          01/04/2011  14:41    <DIR>          info     01/04/2011  17:01               951 log     01/04/2011  14:41    <DIR>          mastertic     01/04/2011  14:41    <DIR>          site     01/04/2011  14:41    <DIR>          street     01/04/2011  17:01               575 verify.wat                    2 File(s)          1,526 bytes           Directory of C:\blah\info          01/04/2011  14:41    <DIR>          .     01/04/2011  14:41    <DIR>          ..     01/04/2011  17:01                80 arc000dat     01/04/2011  17:01               432 arc000nit     01/04/2011  17:01                80 arc001dat     01/04/2011  17:01               576 arc001nit     01/04/2011  17:01                80 arc002dat     01/04/2011  17:01               432 arc002nit     01/04/2011  17:01                80 arc003dat     01/04/2011  17:01             1,008 arc003nit     01/04/2011  17:01                80 arc004dat     01/04/2011  17:01               576 arc004nit     01/04/2011  17:01                80 arc005dat     01/04/2011  17:01               864 arc005nit     01/04/2011  17:01                80 arc006dat     01/04/2011  17:01               576 arc006nit     01/04/2011  17:01                80 arc007dat     01/04/2011  17:01             1,728 arc007nit     01/04/2011  17:01                80 arc008dat     01/04/2011  17:01               432 arc008nit     01/04/2011  17:01             3,420 arcdr9     01/04/2011  17:01                 0 arcnsp                   20 File(s)         10,764 bytes           Directory of C:\blah\mastertic          01/04/2011  14:41    <DIR>          .     01/04/2011  14:41    <DIR>          ..     01/04/2011  17:01                32 dblbnd     01/04/2011  17:01                80 dbltic     01/04/2011  17:01               169 log                    3 File(s)            281 bytes           Directory of C:\blah\site          01/04/2011  14:41    <DIR>          .     01/04/2011  14:41    <DIR>          ..     01/04/2011  17:01            17,732 arc     01/04/2011  17:01               780 arx     01/04/2011  17:01             1,504 cnt     01/04/2011  17:01               452 cnx     01/04/2011  17:01                32 dblbnd     01/04/2011  17:01                80 dbltic     01/04/2011  17:01             2,508 lab     01/04/2011  17:01             1,024 log     01/04/2011  17:01             4,424 pal     01/04/2011  17:01               260 par     01/04/2011  17:01             1,496 pat     01/04/2011  17:01               452 pax                   12 File(s)         30,744 bytes           Directory of C:\blah\street          01/04/2011  14:41    <DIR>          .     01/04/2011  14:41    <DIR>          ..     01/04/2011  17:01             1,288 aat     01/04/2011  17:01             9,268 arc     01/04/2011  17:01               324 arx     01/04/2011  17:01                32 dblbnd     01/04/2011  17:01                80 dbltic     01/04/2011  17:01           352,802 log     01/04/2011  17:01             1,330 nat     01/04/2011  17:01               260 par                    8 File(s)        365,384 bytes      There are no file extensions, but from googling, it suggests that you could convert these using Arc Workstation? Does anybody know if there is another way of loading or converting these files? I have an ArcGIS 10 desktop licence, and a variety of open source tools at my disposal!


title: Is there a way to link eVis with time manager in QGIS

text: I have a point vector layer created from way points of a GPS. I then used the eVis plugin to view geotagged photos for each of these point. I also used the time manager plugin to play through the point layer. Is there anyway that as the time manager shows each point that I can see the photo associated with that point? Thanks for any suggestions!


title: Clear or Disable Cache in ArcObjects Globe Control

text: Is it possible to clear or disable the cache for individual layers on the AxGlobeControl? If so, how? I would like to be able to do this programatically through ArcObjects.


title: Programmatically change geometry column in OpenLayers WFS layer

text: is there a way to programmatically change the geometry column on which a WFS layer is based and displayed? WFS (from GeoServer) is correctly returning me two geometry columns, I want to change "on the fly" the geometry column openLayers uses to display items, obviously having openLayers displaying items (points, in this case) at the new coordinates. I already tried to change the protocol.geometryName field, even forcing reload o refresh on the layer, but without success. Thanks in advance


title: Why can't I "Add label columns to attribute table" etc. in QGIS 1.9?

text: Why can't I "Add label columns to attribute table" and why can't I use "Letter spacing", "Word spacing" and so on? I have a shp-file. QGIS version 1.9 ![enter image description here](http://i.stack.imgur.com/UQbnn.png)


title: Select By Location in ArcGIS using C#

text: I am creating a tool which the user will select point or line features from the map and it will display their names in a list then after clicking on any of these features from the list, it will show the intersected polygons in another list. The problem that I don't know how to do a `SpatialFilter` to do this intersecting functionality, which it will be done like the one in the `Select By Location` tool. Any help will be appreciated, thanks in advance.


title: Adding multiple layers to Tilemill from a single sqlite DB

text: This is a repost of a question I had over at the MapBox support board. Since I wasn't having much luck over there, I'm hoping I would have better success here. I'm having trouble adding more than a single layer from my sqlite/spatialite DB to my project. I have already referred to the previous discussion on this topic HERE and I made sure that my layers all have unique #IDs, yet I'm still receiving an error when I try to add a layer. If I add the layer and use "Save" instead of "Save & Style" then I am able to bring the layer into my map document. However, as soon as I try to apply any styling to it will raise the error. The error:   Sqlite Plugin: 'database idx_tl_2010_51_county10_Geometry is already in use' (C:\Users\User\Desktop\QGIS\Work Progress\Project_Regions.sqlite) Full sql was: 'attach database 'C:\Users\User\Desktop\QGIS\Work Progress\Project_Regions.sqlite.index' as "idx_tl_2010_51_county10_Geometry"'


title: How to add layer with ArcPy and retain feature class metadata (like Credits) in layer properties?

text: I have multiple feature classes in a geodatabase, each feature class has metadata including a Description, Credits, and other information. If I use the 'Add Data' button in ArcMap 10.1 to add one of feature classes then that metadata is imported into the layer description. Ie, I view the layer properties and Description, Credits, etc., are filled in. However, if I use arcpy and add the feature class using `MakeFeatureLayer_management` and `AddLayer` then that metadata is not imported. Does anyone know if its possible to get that full functionality of feature classes and their description, credits, etc to load with arcpy? Am I doing something wrong? Or, must I export the metadata, read the sections I want, and add that text to the layer file? As requested, a code sample is:               mxd=arcpy.mapping.MapDocument(path_to_blankmxd)     df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]     arcpy.MakeFeatureLayer_management(featureclass_sde,activation+"_"+featureclass)     addLayer=arcpy.mapping.Layer(activation+"_"+featureclass)     addLayer.name=layername     arcpy.mapping.AddLayer(df,addLayer)      The feature class to be inserted in the MXD is 'featureclass_sde' with full metadata. However, once its inserted all metadata is lost. At present I am exporting metadata as XML, parsing XML to get metadata elements (description, source, etc), and then inserting those elements into the MXD. It works but seems so.... wrong....


title: Display raster outlines

text: I guess mine is a very simple task, however I don't know how to solve it. I m trying to display the outlines of a raster created from a polygon shp. What I would like to come up with is an image where I can visualize the outlines of each poligon within my raster image. This image has two different colours (namely two "classes") and one is switched off, so only one is visible. Here is an example of what I need (but this is a shp, my image now only displays the yellow parts, while I also need the black lines): ![enter image description here](http://i.stack.imgur.com/XetE0.jpg) I work both with ArGIS 10.1 and QGIS 1.8. Any idea?


title: How to generate a shapefile containing an UTM grid (graticule)?

text: I'm looking for a shapefile containing an 1000m or 2000m UTM grid ("graticule"). All I was about to find on the web is a 10000m UTM grid (here: http://en.giswiki.net/wiki/UTM) and I am aware that a 1000m UTM grid for the whole world might be quite big. That's why I'm now wondering if and how it would be possible to generate such a shapefile for a part of the world. I would prefer to use open source tools like gdal. Others output formats than a shapefile might be ok too. Has anyone an idea?


title: Colouring areas between (vector) contours?

text: Is it possible to colour areas between (vector) contours without making it a DEM first? Joshdoe seems to have done so in Create filled contour map with GRASS but does not say how.


title: How to create a point grid within a polygon layer?

text: Want to generate a point layer, that contains the maximum points in a polygon layer. all points should have the distance of 6 meters between each other. how can i do it?


title: How to upgrade from Windows 7 to 8 without losing ArcGIS for Desktop 10.1 Student Version License?

text: I have Arc GIS installed in my windows 7 with 1 year license key student version.I have recently installed licensed Arc GIS 1 month ago.Now i want to update my laptop to windows 8. Now after installing windows 8 does the same license works since the license key is the same? hoping for response Eliza


title: Issue with address locator in ArcGIS 10.1

text: I am attempting to geocode a table of addresses for a client. I researched into how to create your own address locator and have two separate reference sources: one that has exact address information for our county and another that has road information. I have made a locator for addresses using the US Address - Single House style, and one for roads (I experimented with one range/dual ranges style). When I make a composite locator using both of these and add it to ArcMap during geocoding, I am not presented with an "Address" field input. I only have street or intersection, city, state, zip. What setting exactly will allow me to input the physical house number? Edit: I should mention that the client data has addresses split into multiple fields such as: House_NO, ST_NAME, CITY. Should I simply combine the house_no and st_name fields?


title: QPJ files on ESRI software

text: We have some GIS data, which was converted from AutoCAD dwg to TAB using mapinfo, then converted to shapefiles as that's what we are required to create. Editing the data fields was done in the shape files using QGis. Now for some of the files, the coordinate system settings created using MapInfo were wrong, so I tried to solve this by using the 'save as' option on the edited shapefiles from within QGis, setting the correct projection. For me, all lines up as it should. However, for the one requiring the data (using Esri software) it doesn't. I _think_ this is due to QGis creating a .qpj file along the .prj file, and QGis knows to look at the qpj but esri software doesn't. First of all: **Could this interpretation be correct?** Second: if it is correct, **how do I 'fix' this** so it aligns as it should in the Esri software?


title: Python script for running ArcGIS10 Area Solar Radiation tool

text: I am trying to run the Area Solar Radiation tool using a Python script, because I want to run it on several DEMs at once. I am using a server with Windows 7, ArcGIS 10 and PyScripter. I am a beginner with Python and ArcGIS and for some reason my script doesn't work. In most of my attempts I end up with an error saying that there is a problem with executing AreaSolarRadiation. When I try running it on only one DEM it runs OK, but the outputs are empty files. My guess is that there is some problem either with the glob command, the outGlobalRad.save command, or simply with the solar radiation command. But really I have no clue. I'd be very grateful for any corrections/advice on how to deal with this! Thanks! JJ This is my script:               #Name: AreaSolarRadiation_example02.py     #Description: Derives incoming solar radiation from a raster surface.     #Outputs a global radiation raster and optional direct, diffuse and direct duration rasters     #for a specified time period. (April to July).     #Requirements: Spatial Analyst Extension     #Author: ESRI               #Import system modules     import arcpy     from arcpy import env     from arcpy.sa import *          #Check out the ArcGIS Spatial Analyst extension license     arcpy.CheckOutExtension("Spatial")          #To run on many shapefiles in one folder     import glob          #Set environment settings     env.workspace="Z:/Personal/MyDropbox/GISwork_DQ_JJS/Analysis/UK/Results/microRep/GISoutput/PolygonToRaster/"             folder="Z:/Personal/MyDropbox/GISwork_DQ_JJS/Analysis/UK/Results/microRep/GISoutput/PolygonToRaster/"               # Set local variables     inRaster = name     latitude = 51     skySize = 200     timeConfig = TimeMultipleDays(2008, 5, 160)     dayInterval = 14     hourInterval = 0.5     zFactor = 1     calcDirections = 32     zenithDivisions = 16     azimuthDivisions = 16     diffuseProp = 0.7     transmittivity = 0.4     outDirectRad = ""     outDiffuseRad = ""     outDirectDur = ""     #outDirectRad = Raster("Z:/Personal/My Dropbox/GISwork_DQ_JJS/Analysis/UK/Results/microRep/GISoutput/AreaSolRad/outDirectRad")     #outDiffuseRad = Raster("Z:/Personal/My Dropbox/GISwork_DQ_JJS/Analysis/UK/Results/microRep/GISoutput/AreaSolRad/outDiffuseRad")     #outDirectDur = Raster("Z:/Personal/My Dropbox/GISwork_DQ_JJS/Analysis/UK/Results/microRep/GISoutput/AreaSolRad/outDirectDur")               # Execute AreaSolarRadiation     outGlobalRad = AreaSolarRadiation(inRaster, latitude, skySize, timeConfig, dayInterval, hourInterval, "NOINTERVAL", zFactor, "FLAT_SURFACE", calcDirections, zenithDivisions, azimuthDivisions, "UNIFORM_SKY", diffuseProp, transmittivity, outDirectRad, outDiffuseRad, outDirectDur)          # Save the output     outGlobalRad.save("Z:/Personal/MyDropbox/GISwork_DQ_JJS/Analysis/UK/Results/microRep/GISoutput/AreaSolRad/outGlobalRad.tif")


title: How to use (make sense of) results from ArcGIS Area Solar Radiation tool?

text: I am using the tool Area Solar Radiation in ArcMap10 to analyse DEMs of parts of a city (area equivalent to a neighbourhood). My goal is to calculate the average annual solar irradiation (Wh/m2) on the roofs of buildings. When I run the tool, the output I get is a raster with the global radiation (direct + diffuse) for each raster cell of the input DEM. So far, so good! I can click on a raster cell and see the value of the radiation. But (question A) how can I 'extract' all the values in a way that I can use them (e.g. in a table per building or for the whole study area) and calculate the average solar radiation? And (question B) how can I export the pretty coloured image in a useful format? (I tried exporting it as .tiff and got a black & white image). Thanks! JJ


title: Python/GDAL: Best way to read in raster for use in A* search algortihm

text: I want to create a least cost path program with Python and GDAL since I couldn't find any open source cost path additions for use in Python. As recommended here the best approach is to use A* search algorithm using slope as the cost between generated nodes. I want to read in my Slope Raster with GDAL. What is the best way to read in the Raster in order to further process it in my (not written yet) python A* search algorithm? So far I read it in as a array but I am not sure if this is the best:               import gdal          # Open data     input_raster = "slope2.tif"     raster = gdal.Open(input_raster)          # read raster as array     bandraster = raster.GetRasterBand(1)     datamask = bandraster.ReadAsArray()


title: programmatically add field to shapefiles and insert shapefile name in every row of field

text: **situation** We have a folder containing 1 sub-folder for every participant in a study. The subfolders are named with the participants' ID, e.g., 1000, 1001, 1002, 1003. Every participant folder contains a shapefile (and associated files) of the same name (e.g., 1000.shp, 1001.shp, etc). Each shapefile has several points. **challenge** We need to add a column called "PID" to every shapefile. Then we need to add the participant's ID (i.e., the file name, 1000, 1001, 1002) to every row. For instance, in the folder "1000", there is a shapefile called "1000.shp". In the associated data table, we need to add a column called "PID". Then we need to insert the value "1000" for every row in this new PID column. **question** Is there a way to automate this process? **program/languages** I work with ArcGIS, QGIS, and R. I use R extensively, but not for spatial work. That said, I might have an easier time understanding a programmatic approach in R if one exists.


title: ArcGIS License - ArcInfo level license set during python script execution vs via Desktop Administrator

text: We have an ArcInfo floating license available. Can I check out an ArcInfo level license via arcpy. Here, from the arcgis help, it says you can import a license level for stand-alone scripts. If a user has ArcView installed/running on their machine and they execute a custom script tool that requires ArcInfo licensing, in arcpy, can we browse to a floating license on a license server somewhere and set the product there. I don't see any methods in the help. Otherwise, is my only option to have the user close down all ArcGIS instances, navigate to Desktop Administrator, browse to the server, open ArcMap with ArcInfo and then run the script tool? I would like to avoid this. Thanks.


title: Model builder question regarding iterations (newbe : ( ..... )

text: Ultimate goal is to pull data from a Feature Class, group out the different types of Lease Status (field), run a dissolve and sum calculation on each new layer/FC then merge all new layers/FC's together into one FC. This is my manual process for what I am trying to accomplish: I went about creating a table that includes one Tract_ID, the lease_status and a sum calculated on the Percent_interest field for all the owners. (probably very rudimentary LOL) _(create initial tables to work from)_   * Export DBLeaseOR8312m into a new layer/FC for each lease_status for example DBLeaseOR8312m_open, DBLeaseOR8312m_ThirdParty, etc.   * From the DBLeaseOR8312m_open FC select all other lease_statuses and delete   * Run dissolve on the new layer     * Select Tract_ID as the dissolve field     * Select percent_interest and “SUM” in the statistics area _(modify new tables to include new field)_   * This creates a new table that only has one Tract_ID, one type of Lease_status and the percent_interest calculated for all owners and included in percent interest field.   * Creat new field in this table for Lease_status and add the specific type to that field   * Do this for each Lease_status type _(merge data together into one final table)_   * Then merge all tables into one table to include all lease_status types. I'm thinking I can use Iterate Feature Selection in model builder but not quite sure how to bring it all together and set my outputs...


title: (Solved) Near/Spatial Join pick the wrong polyline segment as closest

text: I have a set of points near a street centerline. I need to find the perpendicular point on the street centerline, and use their x y instead of the GPS information of the original points. I tried Near tool and spatial join in Arcgis 10.1. They returned the same result of a segment which is far away from the right one. For example, for point (60.0001, 40.0001), I could see from the street centerline attribute table that there's a segment start from (60.0002, 40.0002). Yet the Near tool and Spatial Join tool return the same wrong result, which a segment start from(60.002, 40.002) This could be a neighborhood away. Anyone has any idea? Any other intuitive tools I could use? Thanks!!


title: Error with save path. Does not exist or is not supported

text: The final save raster file errors without message with Does not exist or is not support. But it does exist. If you can find anything wrong with this script I would be very greatfull. I had a vast array of errors with this bit of code. But much of it is already working in other parts of the script.               path ="C:/gisdata/india/imagery/"     parentdir ="off"          print "===== using Not isfile and listing folder names only ======="     for fname in os.listdir(path + parentdir):       if not os.path.isfile(os.path.join(path + parentdir, fname)):     #   pathfolders.append(os.path.join(path, fname))         dirname.append(fname)          print dirname          #Loop through image directories and list the tif raster files within each and     #process to reflectance values     for dir in dirname:         print "Setting workspace"         #Set workspace to current working directory         env.workspace = path + parentdir + "/" + dir         print "Workspace set!"                    print "Creating local variables for input bands to Dark Pixel Subtraction"     #Local variable = Raster("LE71460402010085ASN00_B1.TIF")     Rt_B1dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b1.TIF")     Rt_B2dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b2.TIF")     Rt_B3dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b3.TIF")     Rt_B4dp = Raster(path + parentdir + "/" + dir + "/" + "Rt_b4.TIF")          sample = "C:/gisdata/india/imagery/sample"     bands = ("Rt_B1.TIF", "Rt_B2.TIF", "Rt_B3.TIF", "Rt_B4.TIF")     dpRasters = (Rt_B1dp, Rt_B2dp, Rt_B3dp, Rt_B4dp)     fieldNames = ("Rt_b1", "Rt_b2", "Rt_b3", "Rt_b4")          for dpRaster in dpRasters:         for band in bands:             for fieldName in fieldNames:                 #Sampling raster values using select points from Delhi_dp.shp                 print "Sampling raster values using selected points"                 arcpy.sa.Sample( band,"C:/gisdata/india/imagery/Delhi_dp.shp","sample","NEAREST")                      queryString = '"{0}" IS NOT NULL AND "{0}" > 0'.format(fieldName)                 #Use SearchCursor to query non NULL values for creating a new variable                 with arcpy.da.SearchCursor("sample", (fieldName),queryString) as cursor:                     for row in cursor:                         dpValue = row[0]                              print "Starting dpRaster subtraction"                         if arcpy.Exists(dpRaster):                             arcpy.Delete_management(dpRaster)                         print "Subtracting " + str(dpValue) + " from " + str(dpRaster)                         raster_dp = dpRaster - dpValue                         print "Saving " + str(dpRaster) + " as image file"                         raster_dp.save(path + parentdir + "/" + dir + "/" + dpRaster + ".TIF")                         print str(dpRaster) + " has been created"


title: Cartodb: How to make dates appear in Torque time slider

text: I am trying to use my date column as a time-slider in my cartodb map. But it doesn't appear or seem to be an option. What am I doing wrong? Here's my map: Map Here's my data: ![Data structure](http://i.stack.imgur.com/1R3ic.png) Thanks.


title: How to redraw / refresh map on canvas resize?

text: I've created a web page that contains a button to set the map full screen and back to its normal form, the problem map doesn't set to full screen until I resize the output window, how to overcome this, Any guidance and help is appreciated. Here is the working example: Demo


title: How to layer statistics from Excel onto a map?

text: Would another free app or program be better and easier to do this than mapbox.com ? I simply want to layer census 2010 data by shading areas onto a map like on theatlanticcities.com/posts/map/ & developmentseed.org.


title: Layer cannot be enabled in geoserver after restart the computer

text: I use PostGIS database, mapserver is Geoserver, everytime I create a workspace to upload layers into geoserver, I can see the layers be enabled. I can see them in "Layer preview". When I restart my computer, the layers disappear, I checked in "Layers" of "Data", the layers all were not enabled, I had to upload again. I don't understand what happen! I install geoserver as a service, webserver as Apache. Please help me figure out where I did wrong! I appreciate everything yours help! (sorry my English seem to be bad)


title: What are the benefits of RGBI Geotiffs vs ECW format?

text: I just received a batch of imagery for my project area in 2 formats: Geotiff (RGBI) and ECW. The imagery looks very similar and covers the same footprint, but the file size difference is massive: 270GB for the Geotiffs vs 11GB for the ECW. Is there any added value in the RGBI Geotiffs I should be aware of? Is it even worth taking up this kind of server space to store this data? The intended uses include map making, landscape change monitoring, visualization, and supporting engineering work... but always looking for new & interesting ways to make use of the information we have.


title: Building with MinGW 4.8

text: I'm working on a project where we would like to add QGIS features to an existing Qt application. The Qt application has been built with Qt5.2.1 using MinGW 4.8 and thus I think I need to build QGIS with MinGW 4.8. If I'm wrong on this I'd be glad to hear what can be done. If I'm right, my problem is that the QGIS building-from-source tutorial (https://raw.github.com/qgis/QGIS/master/INSTALL) is indeed outdated as far as MinGW is concerned and the given msys environnement is not working with the current QGIS cmake files. I tried to edit those but some libraries are not in msys (spatialite and spatialindex) and it seems linking to Qt5 instead of Qt4 isn't working (problem with QT4_WRAP_CPP). I thought I'd just try to build spatialite myself first but the procedure seems to take ages (http://www.gaia-gis.it/gaia-sins/mingw_how_to.html) and I'm not even sure it would work in the end. If there's anything else I can do (a QGIS plugin for Qt for instance) I'd forget about building without an inch of regret so do not hesitate to share. The only thing to consider is that I'd rather keep Qt5.2.1 because it includes OpenGL and Android support which are both probably going to be needed for further developments. In case nothing can be done for QGIS to work with my project, do you know any other GIS API that has comparable features ? Thanks for your help !


title: How to analyse point data as covariate for areal data

text: I have two data sets: a) a set of polygons with areal data variables.   b) point measures of various co-variates of interest (note I do not necessarily have a point for each polygon) How do I correctly analyse this data set for correlations between an areal variable (dependent variable) and point measures (independent variable)? Unsure what techniques I need to use. Working in R. Edit on 7 July 2014: What I'm really wondering is - if you have areal data i.e. some variable that depends on an area and not points - so something related to population density for example - and you want to analyze that versus a point measurement - rainfall at specific measurement points for example. What is the correct/best way to do this ? I've thought about it since and I see two approaches: 1) Krige the point data and get an average value for each polygon (I presume by picking multiple sample points within the polygon and getting an average) - then build polygon based models using the kriged estimates as covariates. 2) Treat the areal data as point estimates of a field taking the centroid of each polygon as the location for each point estimate, then build a geostatistical model using all data as point estimates. Both these approaches strike me as "incorrect" - what I'm asking really is which is least "incorrect" ?


title: Build Openlayers 3 on Windows

text: On GitHub is a wonderful installation instructions, but most users GitHub, use Linux, it's not worked on Windows.... https://github.com/openlayers/ol3/wiki/Developer-Guide How can I build OpenLayers 3 on Windows? Can't install Pystache               Traceback (most recent call last):       File "C:\Openlayers\pystache-0.5.3\setup.py", line 111, in <module>         import setuptools as dist     ImportError: No module named setuptools      And Closure Linter               Warning (from warnings module):               File "C:\Openlayers\lib\distutils\dist.py", line 267             warnings.warn(msg)         UserWarning: Unknown distribution option: 'install_requires'


title: How to merge 2 shapefiles in ArcObjects

text: Using ArcObjects SDK for Java, I have the full path of 2 shapefiles: String in1, String in2, and I want to merge them using the geoprocessing tool Merge. I tried to several methods:               IGpValueTableObject vtObject = new GpValueTableObject();     vtObject.setColumns(2);     vtObject.setValue(0, 0, in1);     vtObject.setValue(0, 1, 1);     vtObject.setValue(1, 0, in2);     vtObject.setValue(1, 1, 2);      or               IGpValueTableObject vtObject = new GpValueTableObject();     vtObject.setColumns(2);     vtObject.addRow("'" + in1 + "' 1");     vtObject.addRow("'" + in2 + "' 2");      or               IGpValueTableObject vtObject = new GpValueTableObject();     vtObject.setColumns(2);     vtObject.loadFromString(in1 + ";" + in2);      Before executing the the gp tool:               Merge merge = new Merge();     merge.setInputs(vtObject);     merge.setOutput(out);     IGeoProcessorResult result = gp.execute(merge, null);      In each case, I get this error: AutomationException: 0x80004005 - Unspecified error I guess I am wrong in my way to set up the Merge inputs but I didn't find any samples explaining how to do it, except in C# or VB.


title: Preload tiles in OpenLayers 2 and add them to the map

text: I have an XYZ layer on an OpenLayers 2 map. Now I've written a function to preload all tiles of `map.getExtent().scale(2)`. That works, but when I pan the map, the tiles still appear after a small time gap. I think, this is because OpenLayers changes the DOM not immediately, when a new tile would become visible. Is there a way to tell OpenLayers to scale the visible area of the layer?


title: Export selected raster cells to polygon

text: I have selected certain values in a raster (roughly 1000 cells), and would like to export them to be able to create a polygon. Converting it to points first then to polygons does not work. Can anyone suggest a better method?


title: QGIS convert polyline to fixed width polygon

text: I have a shp file of polylines for trails, and I need to convert them all to be polygons that are 8 ft wide (4 ft from the center of each side of the polyline). I tried a following a few other posts that I thought were doing the same thing, but I always just ended up a big blob of color. Any help would be much appreciated.


title: Adding a custom projection to the epsg.properties

text: I have created an entry at www.spatialreference.org of a customized .prj file. http://www.spatialreference.org/ref/sr-org/7973/ I'm trying to add the generated geoserver entry to the EPSG.Properties file. Upon a restart, Geoserver does not recognize my custom projection that I have attempted to add. Reviewing the logs is not indicating any issues, its as if I didn't change the EPSG.properties file at all. I have confirmed, in the log, that I'm using the right EPSG.properties file. Any guidance would be greatly appreciated. UPDATE: The mystery has deepened, but I think I resolved my immediate problem through trial and error. My EPSG.properties file had another entry in it that I had successfully added months ago, and it was loading fine. So, copied its entry, and then replaced values and parameters from the 7973 projection I wanted. Change, Restart, Check it. Eventually, I was able to change all the parameters, and changed the description, and it shows up in my SRS List. Where the mystery deepens, I then tried to take out the junk entries between my two good entries, and it broke. Here are my EPSG Entries: / _This one works_ / 6948=PROJCS["NAD_1927_StatePlane_Wisconsin_South_FIPS_4803",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",2000040.0],PARAMETER["False_Northing",-780.0],PARAMETER["Central_Meridian",-90.0],PARAMETER["Standard_Parallel_1",42.73333333333333],PARAMETER["Standard_Parallel_2",44.06666666666667],PARAMETER["Latitude_Of_Origin",42.0],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG","6948"]] /*This one does not work, and originated from www.spatialreference.org */ 7973=PROJCS["NAD_1983_HARN_UTM_Zone_18N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",1640416.666666667],PARAMETER["central_meridian",-75.0],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG", "7973"]] /*This was my copied 6948, which I modified one parameter at a time, it worked, but then removing 7973 above broke it! */ 96948=PROJCS["NAD_1983_HARN_UTM_ZONE_18N_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-75.0],PARAMETER["Latitude_Of_Origin",0],PARAMETER["scale_factor",0.9996],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG","96948"]] /*So, then I added it a second time, as EPSG = 97973 */ 87973=PROJCS["PHI_NAD83_UTM_ZONE18N_FEET",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-75.0],PARAMETER["Latitude_Of_Origin",0],PARAMETER["scale_factor", .9996],UNIT["US_Foot",0.3048006096012192],AUTHORITY["EPSG","87973"]] /*Again, deleting the 7973 and 96948 entries broke it again! Not only that, but restoring the old file was also broke! But rather than re-add, I renumbered 97973 to 87973 and it now works. I'm leaving the junk entries in there! */ I don't know what to make of it, but its in there now.


title: OpenLayers: add WFS layer from GeoServer and OSM

text: in a OpenLayers HTML file, i'd like to request the default tasmania data of my GeoServer (runs on localhost) and add them as an WFS Layer. The second layer should be an OSM-layer. This is my code up to now:                       function init() {                      var map;                 map = new OpenLayers.Map('map');                 map.addControl(new OpenLayers.Control.LayerSwitcher({'div':OpenLayers.Util.getElement('layerswitcher')}));                 map.addLayer(new OpenLayers.Layer.OSM());                      var roads = new OpenLayers.Layer.Vector("WFS", {                     strategies: [new OpenLayers.Strategy.BBOX()],                     protocol: new OpenLayers.Protocol.WFS({                         url:  "http://localhost:8080/geoserver/ows?",                         featurePrefix: "topp", //<-- GeoServer workspace name                         featureType: "tasmania_roads",//<-- layer name                         featureNS: "http://www.openplans.org/topp", //<-- Edit Workspace Namespace URI                         srsName: "EPSG:4326",                         geometryName: "the_geom" //<-- geometry field                     })                 });                 map.addLayer(roads);                      map.zoomToExtent(roads.getDataExtent());             }      Concerning the OSM-Layer, there are no problems in the result. But the data of the WFS layer is not being displayed. When i zoom to the WFS data extent (last code row) the map centers at (0,0) in the point of origin. Could it have to do with projection problems? Thanks for suggestions. Thanks in advance


title: Use ESRI Web Adaptor to restrict Rest Services breadcrumbs

text: Am I able to restrict users from following the breadcrumb path on my rest services back to the source and seeing every rest service offered on my server? We are using the Web Adaptor on a 10.1 ESRI server.


title: QGIS 2.0 and Debian Jessie

text: I added the official QGIS repos to my installation of Debian Jessie using the instructions on the official QGIS page. As of now, there appears to be some packages missing from the repos that resolve dependencies required by the package. The following is the output for the aptitude:               $ sudo apt-get install qgis python-qgis qgis-plugin-grass           Reading package lists... Done     Building dependency tree            Reading state information... Done     Some packages could not be installed. This may mean that you have     requested an impossible situation or if you are using the unstable     distribution that some required packages have not yet been created     or been moved out of Incoming.     The following information may help to resolve the situation:               The following packages have unmet dependencies:      python-qgis : Depends: libgdal1 (>= 1.8.0) but it is not installable                    Depends: libqgis2.0.1 but it is not going to be installed      qgis : Depends: libgdal1 (>= 1.8.0) but it is not installable             Depends: libqgis2.0.1 but it is not going to be installed             Depends: qgis-providers (= 2.0.1-3+jessie1) but it is not going to be installed      qgis-plugin-grass : Depends: libgdal1 (>= 1.8.0) but it is not installable                          Depends: libqgis2.0.1 but it is not going to be installed                          Depends: grass642 but it is not installable     E: Unable to correct problems, you have held broken packages.      For those of you curious about what my sources.list looks like:               deb http://ftp.us.debian.org/debian testing main contrib non-free     deb-src http://ftp.us.debian.org/debian testing main contrib non-free          deb http://ftp.debian.org/debian/ jessie-updates main contrib non-free     deb-src http://ftp.debian.org/debian/ jessie-updates main contrib non-free          deb http://security.debian.org/ jessie/updates main contrib non-free     deb-src http://security.debian.org/ jessie/updates main contrib non-free          #QGIS     deb     http://qgis.org/debian jessie main     deb-src http://qgis.org/debian jessie main      The culprit appears to be libgdal1, does anyone have a solution for this outside of compiling all the dependencies by hand?


title: Spatio-temporal Overlay

text: I have a `SpatialPointsDataFrame` that lists satellite-detected fires (`firesutm`), including location and date. I also have a `SpatialPolygonsDataFrame` with county boundaries (`countyutm`). I would like to combine these two files to give a monthly count of fires inside of a county or within X km of its border. The package `spacetime` has an over() method that seems ideal. My first step is to expand the county boundaries by X km:               expanded <- gBuffer(countyutm,byid=T,width=rep(X,times=length(countyutm)))      Then I convert both the fire data and the county boundaries to the classes defined by `spacetime` (irregular for the fires, full grid for the counties):               firest <- STIDF(sp=as(firesutm,'SpatialPoints'),                   time=firesutm$date,                   endTime=firesutm$date,                   data=firesutm@data)     countyst <- STF(sp=as(expanded,'SpatialPolygons'),                 time=seq(from=mindt,to=maxdt,by="month"))      The fires are "instances" in the sense that they are associated with a single point in time, i.e. the time they were detected. The counties, on the other hand, are associated with intervals corresponding to the beginning and ending of each month. Finally, the overlay:               firecounty <- over(firest,countyst,returnList=T)      Since I have expanded the county boundaries, at least some of the fires should fall into multiple county-months. Unfortunately, the code above shows only one county-month for each fire, i.e. all elements of the list `firecounty` have length 1. The expected result emerges when I use `over()` on the `sp` objects:               firecountysp <- over(firest@sp,countyst@sp,returnList=T)      Here, it's clear that some fires will fall into multiple counties after expanding the county boundaries, i.e. some elements of the list `firecountysp` have length > 1. With a little more investigation, I was able to determine that `firecounty` gives only the last county-month that contains a given fire. This stuff is pretty new to me. I feel like I'm missing something obvious but I'm also feeling stuck. I'd be grateful for any insight.


title: Maps zoomed around specific locations (magnifying glass effect)

text: How do I create a static map where the scale varies based on location? Example: to drive from a specific location X in Albuquerque to a specific location Y in Chicago, I want to "zoom" around X and Y as follows:   * Near X and Y, the scale is fairly large, so you can see minor streets, buildings, etc, near X and Y.    * Far away from both X and Y, the scale is fairly small, so you might just see long lengths of highway, but no details.  In other words, it's like a regular map with a "magnifying glass effect" around X and Y. The effect should "fade off" smoothly, so there isn't a break in the map. Reasoning: when traveling between remote cities, the highway travel is often easy (in the sense it doesn't require details), but getting from X to the highway and from the highway to Y can require detailed instructions, knowing minor streets, etc. EDIT: I've attached a horrible example of what I want to do. Notes:   * This only zooms in on one city, not two.   * This is just plain old image zoom. What I would want is more detail in the Albuquerque area.   * Notice, however, that the roads/highways are continuous: the magnifying does not have a "cutoff" problem. ![alt text](http://i.stack.imgur.com/LCbvE.png)


title: MapCruncher makes maps but what scale?

text: Here is my workflow for creating tiles and linking them to a layer.               1) load shapefile     2) grab layer representing shapefile     3) init & execute mapcruncher     4) turn on dynamic display     5) create dynamiccachelayermanager and init with my layer     6) save layer as layerfile      I've managed to successfully modify this sample: http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//000100000mm8000000 It creates my png files and a file called conf.xml. In conf.xml it has the following data:               <LODInfo xsi:type="typens:LODInfo">     <LevelID>1</LevelID>      <Scale>485413.57628360949</Scale>      <Resolution>0.0011550226764674709</Resolution>      </LODInfo>      I thought if I set my mapScale to 485413.57628360949 then it would grab my generated tiles instead of regenerating them, but no luck, it still wants to regenerate them. I think its just a mapscale issue, do you know how to get the mapscale of the generated tiles so I can put it set it via the IMapControl.mapScale interface? * * * Here's what happened next:               1) I used mapcruncher to generate tiles, it also generated a conf.xml     2) I open conf.xml to get the mapscales     3) I hard coded those mapscales in as my zoom levels     4) Dynamically generate the tiles at my hard coded zoom level     5) compare the mapcruncher generated tiles to the tiles which are dynamically created      The MapCruncher files start at Level 0 - 8, and the images are numbered in an orderly fashion. C0000000.png C0000001.png... etc. But the dynamically generated pngs are in the following format C0000003d.png even if only 1 png represents the whole map. One would expect C00000000.png here.... I am truely lost here.


title: Batch Update Fields in Field Calculator or Model

text: I have a dataset of block-level census data that I clipped to an Illinois State Senate district. To update the data (e.g. for population) I'm applying an area multiplier, i.e taking the area of each attribute after and dividing it by the area before and applying the percentage to each field. Since I wasn't aware of how to do a geoprocessing clip that would also clip data I'm left editing the table manually My question is: how do I batch multiply columns against one column? I have 100 columns of data I need to multiply against a percentage contained in an aptly named proportion column. Can this be done through a model in ArcGis, or through a Python script in the Field Calculator? I have both QGIS and a demo license of ArcGis available.


title: Access attributes on featureCollection

text: The geoJSON-String:               {       "type": "FeatureCollection",       "features": [         {           "type": "Feature",           "properties": {"title": "Standort 1 "},           "geometry": {             "type": "Point",             "coordinates": [               13.444175720214844,               52.580522509085554             ]           }         }         ,...       ]     }      This is my code to read and create the features to markers in openlayers:               var geo = new OpenLayers.Format.GeoJSON();          var obj = JSON.parse(featureJson);     var features = geo.read(obj, "FeatureCollection");          alert(features[0].attributes.length);      This works fine and I'm able to create the markers. But what am I missing on this line?               for(var i=0;i < features.length;i++) {         alert(features[i].attributes.title);     }      The output = "undefined". Thanks! :-)


title: Finding coordinates from a nested grid as a numpy array

text: So this is a pretty basic question, but I'm trying to figure it out. I've written a script that allows a user to read in a raster and provide a cell resolution, and outputs the amount of cells necessary to resample the grid (to provide a pre-processing tool to run the WRF weather model). I am still trying to figure out a test to make sure the cell size actually fits the input grid. Anyways, I'm having some difficulty figuring out how to nest grids within a 'master' domain. IDEALLY this would be via user-IO, where multiple grids are selected within themselves. Haven't added this function yet, but the nesting is easy enough with a for loop (though I need to make sure that the nested grids actually fit within the master) But my main question is obtaining the origin points for each of the nested grids in terms of cells of the master domain (shown in the link above). For instance, finding the x,y coordinates of the 0,0 origin of a nested domain WITHIN the larger master-domain. Does anybody have any ideas? Here's my code:               #!/bin/python     import numpy as np     from sys import argv          file = argv[1]     grid = np.loadtxt(file)     array_dims = grid.shape     x = array_dims[0]     y = array_dims[1]          cell_size = raw_input("What is the cell size you want? (x,y) \n >")     cell_size = (cell_size.split(","))     cell_size = [int(i) for i in cell_size] #goes through list and turns them to array     print cell_size     cell_x = cell_size[0]     cell_y = cell_size[0]          print "The shape of the grid is %r by %r \n" %(x,y)          e_we = int(round(x/cell_x))     e_sn = int(round(y/cell_y))          print "With that cell size, the WRF-domain is %r by %r cells" %(e_we,e_sn)          center_x = x/2     center_y = y/2          print "The center point is at %r, %r" %(center_x,center_y)     print "This point's value is", grid[center_x,center_y]


title: Mysterious Data loss after compressing SDE database

text: I suffer sometimes from data loss in SDE database exactly after compressing database. i mean by data loss that some data remains in its prior state before edits, but this problem do not affect the totality of the data, it may touch just a portion of a feature class , and it may not touch all the feature classes. i have made the compressing script to take edits data in priority. i remarked also problems when i try to close to application that i'm using for editing, when i make some edits it takes too much time to close after i close it with close (X) button,so in the most cases i kill the process with task manager. i wonder if in this time the Arcobjects or Arcsde make some back up for edits. may be this is the cause for the Data loss. i still wonder if this problem is the cause of this data loss. **My question :** why my app take so long time to close after i've made some edits ? i'm available for any additional informations. I'm using Arcgis 10.1 , VS 2010 , C#


title: How to hide grid's labels when zooming out?

text: I have used polygon grid layer on my map, and I used to label the layer with ID automatically, and the result was convenient for me: ![enter image description here](http://i.stack.imgur.com/RHknf.jpg) But when I zoomed out, the labels didn't hide at all! and they keep appearing in frustrating way! ![enter image description here](http://i.stack.imgur.com/uTi98.jpg) I would like to get the result when I zooming out like this: ![enter image description here](http://i.stack.imgur.com/uYffa.jpg) What to do in this case to get this result? How to hide grid's labels when I zooming out? My vector grid is 1000 meter polygon.


title: Spatial Join in CartoDB

text: I need to update a table of points with the FIP code of the nearest county polygon from a separate table. I've managed to update the column for those points that Intersect a county polygon with this:               update points p set fips = c.fips      from counties c where ST_Intersects(c.the_geom, p.the_geom)      but I've got about 60 coastal points that, because of generalization in the county polygons, don't overlap a polygon feature. I need to update the fips of these with the NEAREST county polygon's fips code, and I've been futzing around with the <-> Indexed NN, but I just can't get it to work for me. So, where p.fips IS NULL I need to set p.fips = c.fips of the nearest county feature. Any help is appreciated. I'm sure it's mostly my complete lack of understanding of SQL subqueries that is preventing me from doing this successfully.


title: Can't generate leaders with EasyCustomLabeling

text: I'm using EasyCustomLabeling with QGIS 2.0.1 and can't get it to generate the leader lines it is supposed to in Action 2. I create a new layer with EasyCustomLabeling, move labels to new locations, and then save the layer. I'm not sure how to get it to generate the lines to connect the points and labels after this however. Any help would be much apprecaited


title: Use field calculator in QGIS so get number of points in polygone and store it in attribute table

text: I have to layers, one point one polygon (with 6 areas) stored as shapefiles. Now I want to count the number of points in each area and write it into the attribute table of the polygons. A column sum_poi already exists. Can I use the field calculator to populate the column in the polygon layer? Thanks for your help.


title: How to export a MXD to a PDF but only switching the folder

text: Let say I have two variables.                   InputMXD = parameters[0].valueAsText         PDF_Folder = parameters[1].valueAsText      I have my inputMXDs to have multivalues so I'll have around 10 mxds each time the script is run. I want to put the exported pdf's into different directory, but I want to use whatever the mxd name is. The input could be p:/mxds/oneexample.mxd and p:/mxds/onemoreexample.mxd and the output would be z:/pdfs/oneexample.pdf and p:/mxds/onemoreexample.pdf.                   MXDList = arcpy.ListFiles(InputMXD)         for MXDPath in MXDList:             MXD = arcpy.mapping.MapDocument(MXDPath)             arcpy.mapping.ExportToPDF(MXD, PDF_Folder + )      Here's what I have to get me to the last line. I don't know how to add the original file name of the mxd to the outputted PDF. * * * Updated: It will run through it, but it doesn't output anything. The output folder is empty. I think there's something wrong with the loop. Here's what I now have in my loop.                   MXDList = arcpy.ListFiles(InputMXD)         for MXDPath in MXDList:             MXD = arcpy.mapping.MapDocument(MXDPath)             outputPDF = os.path.basename(MXDPath).split('.')[0] + ".pdf"             arcpy.mapping.ExportToPDF(MXD, PDF_Folder + "\\" + outputPDF)      I tried the code below and still no dice. Here's what the results kicked back. ![enter image description here](http://i.stack.imgur.com/MYivV.png) I'm not sure why it's saying the output is empty? * * * Update 2: I updated the code to this:                   arcpy.env.overwriteOutput = True         arcpy.env.workspace = PDF_Folder         MXDList = arcpy.ListFiles('*.mxd')         for MXDPath in MXDList:             MXD = arcpy.mapping.MapDocument(MXDPath)             outputPDF = os.path.basename(MXDPath).split('.')[0] + ".pdf"             arcpy.mapping.ExportToPDF(MXD, PDF_Folder + "\\" + outputPDF)      I changed the workspace to PDF_Folder and it didn't kick back an error. It also didn't output any PDF files. Also, here's how my parameters are setup at the top.                   # First parameter         InputMXD = arcpy.Parameter(             displayName="Input MXDs",             name="InputMXD",             datatype="DEMapDocument",             parameterType="Required",             direction="Input",             multiValue=True)              # Second Parameter         PDF_Folder = arcpy.Parameter(             displayName="Output Folder for PDF",             name="PDF_Folder",             datatype="DEFolder",             parameterType="Required",             direction="Output")      Does my output need to be something besides a DEFolder datatype? * * * Update 3 Here's what my input looks like:   ![enter image description here](http://i.stack.imgur.com/X1PBn.png) Input MXDs is a MapDocument that is a multiValue input. Output Folder for PDF is a DEFolder output. I've tried strings, workspaces, folders. Nothing has worked. I know its something small making this really frustrating. I'm going basic to see if I can get just one mxd outputted to a pdf. Crawl before I can walk. * * * Update 4 I'm trying to run the command directly in ArcMap from the python window. Here is the kickback it's giving me.   ![Python Window in Arc](http://i.stack.imgur.com/Dvm4U.png) * * * Update 5 I have tried just running it in a small script to see if it worked with a file name. It does and it includes spaces, commas, periods, dashes and underscores. Here's the simple testing code that worked. I have no idea what isn't working.                import arcpy           doc = arcpy.mapping.MapDocument("Q://Field Data//6th P.M//COLORADO//6N65W//PROJECTS//6N-65W-12//HOLTON 6-65-12//BEP.mxd")      arcpy.mapping.ExportToPDF(doc, "C://temp//test2.pdf")


title: Plugin Tutorial issue: trying to work around changes in QGIS API

text: I'm currently following the tutorial to create my first python plugin for QGIS: http://www.qgisworkshop.org/html/workshop/plugins_tutorial.html and have similar issues to this thread: QGIS 2.0 Plugin: Select and Identify a Single Feature I'm attempting to list attributes of features in a dialog window when selecting a single feature. I don't receive an error message, but the feature isn't selected, nor is attribute information displayed. I'm only wetting my feet in terms of programming and using this forum, so please do excuse any newbie errors. Any suggestions would be greatly appreciated. I'm posting my entire program as I'm unsure where the problem is coming from:               # Import the PyQt and QGIS libraries     from PyQt4.QtCore import *     from PyQt4.QtGui import *     from qgis.core import *     from qgis.gui import *          # Initialize Qt resources from file resources.py     import resources_rc     # Import the code for the dialog     from zconsodialog import zConsoDialog     import os.path          class zConso:          def __init__(self, iface):         # Save reference to the QGIS interface         self.iface = iface         # reference to map canvas         self.canvas = self.iface.mapCanvas()         # out click tool to emit a QgsPoint on every click              self.clickTool = QgsMapToolEmitPoint(self.canvas)         # initialize plugin directory         self.plugin_dir = os.path.dirname(__file__)         # initialize locale         locale = QSettings().value("locale/userLocale")[0:2]         localePath = os.path.join(self.plugin_dir, 'i18n', 'zconso_{}.qm'.format(locale))              if os.path.exists(localePath):             self.translator = QTranslator()             self.translator.load(localePath)                  if qVersion() > '4.3.3':                 QCoreApplication.installTranslator(self.translator)              # Create the GUI dialog (after translation) and keep reference         self.dlg = zConsoDialog()          def initGui(self):         # Create action that will start plugin configuration         self.action = QAction(             QIcon(":/plugins/zconso/icon.png"),             u"Calculer consommation d'une zone", self.iface.mainWindow())              # connect the action to the run method         QObject.connect(self.action, SIGNAL("triggered()"), self.run)         #self.action.triggered.connect(self.run)              # Add toolbar button and menu item         self.iface.addToolBarIcon(self.action)         self.iface.addPluginToMenu(u"&zConso", self.action)              # connect our custom function to a clickTool signal that the canvas was clicked              result = QObject.connect(self.clickTool, SIGNAL("canvasClicked(const QgsPoint &, Qt::MouseButton)"), self.handleMouseDown)         #QMessageBox.information( self.iface.mainWindow(),"Info", "connect = %s"%str(result) )              # connect our select function to the canvasClicked signal              result = QObject.connect(self.clickTool, SIGNAL("canvasClicked(const QgsPoint &, Qt::MouseButton)"), self.selectFeature)         #QMessageBox.information( self.iface.mainWindow(),"Info", "connect = %s"%str(result) )               def unload(self):         # Remove the plugin menu item and icon         self.iface.removePluginMenu(u"&zConso", self.action)         self.iface.removeToolBarIcon(self.action)          def handleMouseDown(self, point, button):              self.dlg.clearTextBrowser()         self.dlg.setTextBrowser( str(point.x()) + " , " +str(point.y()) )         #QMessageBox.information( self.iface.mainWindow(),"Info", "X,Y = %s,%s" % (str(point.x()),str(point.y())) )          def selectFeature(self, point, button):         QMessageBox.information( self.iface.mainWindow(),"Info", "in selectFeature function" )         cLayer = self.canvas.currentLayer()              radius = 5         rect = QgsRectangle(point.x() - radius,                             point.y() - radius,                             point.x() - radius,                             point.y() - radius)         cLayer.select(rect, False)         selected_features = cLayer.selectedFeatures()              for i in selected_features:             attrs = i.attributeMap()             for (k,attr) in attrs.iteritems():                      print "%d: %s" % (k, attr.toString())              #code from tutorial         # setup the provider select to filter results based on a rectangle         #pntGeom = QgsGeometry.fromPoint(point)           # scale-dependent buffer of 2 pixels-worth of map units         #pntBuff = pntGeom.buffer( (self.canvas.mapUnitsPerPixel() * 2),0)          #rect = pntBuff.boundingBox()         # get currentLayer and dataProvider                        #selectList = []         #if cLayer:         #        provider = cLayer.dataProvider()         #        feat = QgsFeature()         #        # create the select statement         #        provider.select([],rect) # the arguments mean no attributes returned, and do a bbox filter with our buffered rectangle to limit the amount of features         #        while provider.nextFeature(feat):         #            # if the feat geom returned from the selection intersects our point then put it in a list         #            if feat.geometry().intersects(pntGeom):         #                selectList.append(feat.id())              #        # make the actual selection         #        cLayer.setSelectedFeatures(selectList)         #else:         #        QMessageBox.information( self.iface.mainWindow(),"Info", "No layer currently selected in TOC" )                    # run method that performs all the real work     def run(self):         # make our clickTool the tool that we'll use for now              self.canvas.setMapTool(self.clickTool)              # show the dialog         self.dlg.show()         # Run the dialog event loop         result = self.dlg.exec_()         # See if OK was pressed         if result == 1:             # do something useful (delete the line containing pass and             # substitute with your code)             pass      Mahalo!


title: Can I limit the drawing size of icons when my map is scaled?

text: I have created some custom icons that I'm using to represent points on my map. For the case in point, I'm using a skyscraper icon to indicate construction projects underway. It looks cool at certain scales, but when I zoom out farther to get the big picture, my icons just overlap into a blobby mess and all meaning is lost. So, is there a way to limit the scaling at a certain point so that when I zoom out farther than, say, a 15 zoom, my icons will no longer scale any larger?


title: WKT for local mine grid

text: I have been trying to create a WKT entry to use in gdalwarp etc. I need an affine transform from coordinates in MGA94 Zone 54 to our local mine grid. All efforts have so far failed, including using the bits and pieces I have found here. I have a MapInfow.prj file entry that works beautifully but I need to convert our imagery from MGA to mine grid to supply to mining consultants. This entry is below with the affine parameters highlighted: "SRG/AISTR", 3008, 116,7, 141, 0, 0.9996, 500000, 10000000, 7, **0.999770876952, -0.031210618094, -193025.578125, 0.031224025413, 0.999768821523, -7589075.9140625** , -10000, -10000, 500000, 500000 The standard scale, rotate parameters as supplied by our surveyors are here: Rotation: 358.2111111 degrees Scale: 0.999741893 Local East: 16000 Local North: 20000 MGAz54 E: 446231.629 MGAz54 N: 7596899.020 My latest attempt using the affine parameters from MapInfo are here: FITTEDCS["SRG", PARAM_MT["Affine", PARAMETER["num_row, 3], PARAMETER["num_row, 3], PARAMETER["elt_0_0", -193025.578125], PARAMETER["elt_0_1", -7589075.9140625], PARAMETER["elt_1_0", 0.999770876952], PARAMETER["elt_1_1", 0.031224025413], PARAMETER["elt_2_0", -0.031210618094], PARAMETER["elt_2_1", 0.999768821523], PARAMETER["Central Meridian", 140.4800753], PARAMETER["Latitude_Of_Origin", -21.73061814]],   PROJCS["GDA_1994_MGA_Zone_54", GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994", SPHEROID["GRS_1980",6378137.0,298.257222101]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Transverse_Mercator"], PARAMETER["False_Easting",500000.0], PARAMETER["False_Northing",10000000.0], PARAMETER["Central_Meridian",141.0], PARAMETER["Scale_Factor",0.9996], PARAMETER["Latitude_Of_Origin",0.0], UNIT["Meter",1.0],AUTHORITY["EPSG",28354]] Hope one of you out there can help. Cheers, Tony As requested by Andre here are some control points: > >     Local East    Local North MGAz54 East MGAz54 North >     15896 27965   446376.2852 7604861.329 >     12458 28452   442956.0464 7605455.262 >     35698 27412   466146.2673 7603690.743 >     15698 29654   446231.1445 7606555.249 >     16547 23456   446886.0792 7600335.373 >


title: Add PostGIS spatial functions to a custom schema other than "public" in PostgreSQL

text: I am new to PostgreSQL and PostGIS. Recently, I created a PostGIS 2.0.3 database on a PostgreSQL 9.1 database server using pgAdmin. The "PostGIS" extension was found installed in "Extensions". Al spatial functions were added to the "public" schema. That's fine. Now I want to store all my data into a new scheme called "gc". However, how can I make all spatial function installed in that "gc" schema? There is no single function in the schema. However, when I tried to import/new a feature class from ESRI ArcCatalog 10.1 to this schema, it worked! The feature class could be imported and displayed in QGIS. Could anyone give me any tip or idea about it? Thanks in advance!


title: Could you suggest software to draw a cross section, companion to QGIS?

text: If you use GIS systems, what software do you turn to when there is a requirement to draw a cross section? Ideally a workflow that includes QGIS. Typically there are some points which will coexist in the plan and section view. But often I find there is much more detail to be input in the section which is redundant to have in x y coordinate system. Related to this question How to draw a cross section in QGIS but what is your experience? Perhaps you do street plans and need to represent the road and verge profile from time to time. Not Arc addons please.


title: Replace google map MARKER with HTML element

text: I want to add a custom HTML element (that can become transparent dynamically) instead of a marker. Basically my html object would be a cirle (or rounded pin) with a picture of users in the middle. I have imagined different solutions:   * place a label (containing HTML) covering the marker which visibility would be set to none   * work with overlays   * hard way: replace the marker with a PNG image that I should compose dynamically with some sort of library superimposing two images (one geometric and the other user picture). Anyway I can't figure out what to do and how to do it. Thanks, Stefano


title: Is it possible to integrate (Cluster layer + Time slider) in ArcGIS JSAPI?

text: I am using ESRI suit 10.1 (ArcGIS Desktop + Server) and JSAPI 3.5/3.6 I am trying to integrate time aware data and cluster layers in single map but both are unable to display. According to web help: -Cluster layers are created using feature layers and generated graphics -Searching sample/parallel example but not able to find on google/ESRI samples So any help/direction well be great!! Thanks in advance


title: Converting multilayer KML files to GeoJSON using ogr2ogr

text: I have a kml file that I'm trying to convert to GeoJSON. However, it turns out there's multiple layers, and `ogr2ogr` isn't happy.               $ ogr2ogr -f GeoJSON shapes.json shapes.kml     ERROR 6: GeoJSON driver doesn't support creating more than one layer     ERROR 1: Terminating translation prematurely after failed     translation of layer Set 2 (use -skipfailures to skip errors)      Here's a sample of one such multi layer KML file:               <?xml version="1.0" encoding="UTF-8"?>     <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">     <Document>         <name>shapes.kml</name>         <Folder>             <name>Shapes</name>             <Folder>                 <name>Set 1</name>                 <Placemark>                     <name>Shape 1-1</name>                     <Polygon>                         <tessellate>1</tessellate>                         <outerBoundaryIs>                             <LinearRing>                                 <coordinates>                                     -81.97209580036224,40.70902818974323,0 -81.97211148995324,40.70705812260233,0 -81.97204657426639,40.70698429917707,0 -81.97170065320428,40.70689762261313,0 -81.97128823249281,40.70678591732805,0 -81.97012759802162,40.70683180630405,0 -81.97011100460675,40.70907869772486,0 -81.97209580036224,40.70902818974323,0                                  </coordinates>                             </LinearRing>                         </outerBoundaryIs>                     </Polygon>                 </Placemark>                 <Placemark>                     <name>Shape 1-2</name>                     <Polygon>                         <tessellate>1</tessellate>                         <outerBoundaryIs>                             <LinearRing>                                 <coordinates>                                     -81.9678256998772,40.70699525577179,0 -81.96710537786173,40.70714109117483,0 -81.9671339172961,40.70964816832036,0 -81.9681094085136,40.70953487937989,0 -81.96849455277558,40.70941143819032,0 -81.9687122721291,40.70922702833104,0 -81.96869523424586,40.70693622825318,0 -81.9678256998772,40.70699525577179,0                                  </coordinates>                             </LinearRing>                         </outerBoundaryIs>                     </Polygon>                 </Placemark>                 <Placemark>                     <name>Shape 1-3</name>                     <Polygon>                         <tessellate>1</tessellate>                         <outerBoundaryIs>                             <LinearRing>                                 <coordinates>                                     -81.97575255412974,40.70901137563123,0 -81.97498009801741,40.70910891904678,0 -81.97366868385656,40.70909911614361,0 -81.97156793031792,40.70919907550063,0 -81.96927940567542,40.70932415940178,0 -81.96979533529213,40.70965324506344,0 -81.97175061906975,40.70968631132654,0 -81.9733120082854,40.70970666237672,0 -81.97576756497214,40.70956450955584,0 -81.97575255412974,40.70901137563123,0                                  </coordinates>                             </LinearRing>                         </outerBoundaryIs>                     </Polygon>                 </Placemark>                 <Placemark>                     <name>Shape 1-4</name>                     <Polygon>                         <tessellate>1</tessellate>                         <outerBoundaryIs>                             <LinearRing>                                 <coordinates>                                     -81.96434152341726,40.70702598013048,0 -81.96356422472108,40.70713349505342,0 -81.96297679852724,40.70696883424459,0 -81.96309807638089,40.70812857198734,0 -81.96343567507532,40.70928783377635,0 -81.96436703230279,40.70914919861734,0 -81.96460689239605,40.70921100214736,0 -81.9650923971045,40.70894923436187,0 -81.96555478608285,40.70876312686659,0 -81.96552441635622,40.70713152110199,0 -81.96434152341726,40.70702598013048,0                                  </coordinates>                             </LinearRing>                         </outerBoundaryIs>                     </Polygon>                 </Placemark>             </Folder>             <Folder>                 <name>Set 2</name>                 <Placemark>                     <name>Shape 2-1</name>                     <Polygon>                         <tessellate>1</tessellate>                         <outerBoundaryIs>                             <LinearRing>                                 <coordinates>                                     -81.98834021695214,40.69452103418363,0 -81.98812651732483,40.68978476540717,0 -81.98732201929137,40.68982207193159,0 -81.98748028381483,40.69454809335561,0 -81.98834021695214,40.69452103418363,0                                  </coordinates>                             </LinearRing>                         </outerBoundaryIs>                     </Polygon>                 </Placemark>                 <Placemark>                     <name>Shape 2-2</name>                     <Polygon>                         <tessellate>1</tessellate>                         <outerBoundaryIs>                             <LinearRing>                                 <coordinates>                                     -81.98570709620523,40.69457415615653,0 -81.98548489260064,40.68973423650695,0 -81.98495179849937,40.6899203348901,0 -81.98523118119185,40.6946028701411,0 -81.98570709620523,40.69457415615653,0                                  </coordinates>                             </LinearRing>                         </outerBoundaryIs>                     </Polygon>                 </Placemark>             </Folder>         </Folder>     </Document>     </kml>      Are there any suggestions for successfully processing the shapes.kml using existing libraries, or must it be split apart manually? What is your suggested strategy for doing so?


title: Crashing trying to warp map files

text: I have been trying to warp ECW and Ozfx3 map formats to GeoTiff format. I have navigated to `Raster>Projections>Warp(reproject)` I have then selected the file to warp the output that I want and the process starts but each time the process fails within 5 seconds with the dialogue "The process crashed some time after starting successfully". So far I have had no success could someone assist me please. Regards Bart.


title: How to include highways type track or service in osm2po?

text: What must I do for osm2po to take into consideration highways type track or service? I have processed one file and in the results with QGIS I can't see those types of highways. Thanks!


title: Why are detailed model processes and tools being shown in the progress window of Model Builder when it is double clicked / opened?

text: ArcGIS Desktop Advanced 10.1 The ESRI help notes state that when a model is run from"edit" in Model Builder all results are shown in the progress window - e.g. the progress of every tool and process. The help notes also state that when the model tool is run from "open" or double clicked, only the overall progress of the model is shown in the progress window, not the progress of the individual tools or processes within the model. http://resources.arcgis.com/en/help/main/10.1/index.html#/Running_a_model_tool/002w0000004q000000/GUID-17578415-7F89-4BD4-BB9C-33B8F2220FEF/ This is the behavior I normally see for my Models. However I have some models that when they are double clicked in ArcCatalog / or "open" when they run they show the full progress of every tool and process in the dialogue box. Any ideas or experience of this behavior? thanks for any tips / insight


title: Copy Z-coordinate from S57 SOUNDG to Shp

text: I'm trying to create a shapefile from an S57 SOUNDG Layer, where i want to put in in a separate column the Z-coordinate of the SOUNDG Layer in the resulting Shapefile....any ideas on how to achieve that?.


title: Does ST_GEOMETRY support Linear Referencing

text: I am wondering if the ArcGIS ST_GEOMETRY provides Linear Referencing functions, such as creating a point based on a given measurement, like Oracle Spatial. I couldn't find any function like that in the ArcGIS documentation. But there are pretty functions in ArcPy. Thanks, Fred


title: Where can I find osm2po javadoc

text: I would like to build an offline routing engine to compute route distance between two coordinates. I found osm2po but I didn't found any sort of documentation. As far as i understood osm2po can be used as java library. Where can I find it? Because I didn't found any way to use osm files to find a path. There's others valid offline engines that I can use? Or there's some sort of input paramenters to use with osm2po to get the distance? Thanks


title: OpenLayers CQL filter gets info about NOT VISIBLE features

text: I have a WMS layer and a control WMSGetFeatureInfo.   If I apply a CQL Filter to the layer and I redraw it, the WMSGetFeatureInfo control continue querying also NOT VISIBLE features.   Why in your opinion? How can I do to fix the problem? **EDIT**               info = new OpenLayers.Control.WMSGetFeatureInfo({      url: 'myserver:8080/geoserver/gis/wms';,       title: 'Identify features by clicking',       layers: [layer_1,layer_2,layer_3],       infoFormat: 'application/vnd.ogc.gml',       queryVisible: true     });          info.events.register("getfeatureinfo", this, callFunct);