texts:
I have an arcpy script that works on one machine, but fails very early on another.

import arcpy, os

arcpy.env.overwriteOutput = True

xlFile = r"D:\GIS\mt\data\ExcelTable.xlsx"
outGDB = r"D:\GIS\mt\data\Data.gdb"
arcpy.env.workspace = outGDB
mapTable= "xxExcelOutput"
inMemMapTable= r"in_memory\xxExcelOutput"

sheetName = r"Map Sheet"

arcpy.ExcelToTable_conversion(xlFile, inMemMapTable, sheetName)

It fails on the arcpy.ExcelToTable_conversion() with the following error message:

Traceback (most recent call last):
  File "D:\GIS\mt\readexcel.py", line 13, in <module>
    arcpy.ExcelToTable_conversion(xlFile, inMemMapTable, sheetName)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\conversion.py", line 60, in ExcelToTable
    retval = convertArcObjectToPythonObject(gp.ExcelToTable_conversion(*gp_fixargs((Input_Excel_File, Output_Table, Sheet), True)))
ExecuteError:
Traceback (most recent call last):
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 252, in <module>
    arcpy.GetParameterAsText(2))
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 203, in excel_to_table
    out_fields = gen_out_fields(workbook, sheet, out_path, is_gdb)
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 175, in gen_out_fields
    out_fields.append(clsField(f, out_path, is_gdb, out_fields))
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 33, in __init__
    if not name[0].isalpha():
IndexError: string index out of range

Failed to execute (ExcelToTable).

The machine it works on is running 10.3.1, and it fails on 10.5.  I thought maybe there was a change/difference in parameter requirements in 10.5, so in ArcMap I ran just the Excel to Table tool from ArcToolbox

and I get pretty much the same error:

Executing: ExcelToTable D:\GIS\mt\ExcelTable.xlsx in_memory\xxExcelOutput "Map Sheet"
Start Time: Tue Mar 21 09:30:15 2017
Running script ExcelToTable...
Failed script ExcelToTable...

Traceback (most recent call last):
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 252, in <module>
    arcpy.GetParameterAsText(2))
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 203, in excel_to_table
    out_fields = gen_out_fields(workbook, sheet, out_path, is_gdb)
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 175, in gen_out_fields
    out_fields.append(clsField(f, out_path, is_gdb, out_fields))
  File "c:\program files (x86)\arcgis\desktop10.5\ArcToolbox\Scripts\ExcelToTable.py", line 33, in __init__
    if not name[0].isalpha():
IndexError: string index out of range

Failed to execute (ExcelToTable).
Failed at Tue Mar 21 09:30:16 2017 (Elapsed Time: 1.20 seconds)

What could be causing my arcpy to not be able to read my Excel file, and how do I fix it?

The Excel file was originally an export from Google Sheets, so I re-saved it in Excel with no effect.  I removed spaces from the file name and the sheet name, also with no effect.

I created another Excel file with a similar filename and sheet name but with dummy data, and this has worked successfully, so I suspect now there's something in my original file that the Excel to Table tool doesn't like, but I am unsure what to look for, particularly because the exact same file works fine on another machine.


texts: Access-Control-Allow-Origin in geoserver and ol2

texts:
I am trying to get a pop up from my map with the getfeauture code. I tried to followed the following example. When i tried to run the code i get the following error in my console:
XMLHttpRequest cannot load http://localhost:8080/geoserver/wms?LAYERS=test%3Aput_25_vlak_1_vulling&QUE…04&FORMAT=image%2Fpng&INFO_FORMAT=text%2Fhtml&SRS=EPSG%3A28992&X=731&Y=197. Response to preflight request doesn't pass access control check: No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'null' is therefore not allowed access.

I found some simular problems online like Access-Control-Allow-Origin in geoserver and ol3 but when i add a proxy i get the following error:
XMLHttpRequest cannot load file:///C:/Program%20Files%20(x86)/GeoServer%202.8.1/data_dir/www/proxy.cgi…g%26INFO_FORMAT%3Dtext%252Fhtml%26SRS%3DEPSG%253A28992%26X%3D751%26Y%3D189. Cross origin requests are only supported for protocol schemes: http, data, chrome, chrome-extension, https, chrome-extension-resource.

Is there somebody that can help me?
My code is as following:
<html>

  OpenLayers Example

    OpenLayers.ProxyHost = "proxy.cgi?url=";

 var bounds = new OpenLayers.Bounds(
                        189796.953125, 427589.40625,
                        189836.1875, 427625.5625
                 );
             var options = {
                             controls: [],
                             maxExtent: bounds,
                             maxResolution: 0.02741796875,
                             projection: "EPSG:28992",
                             units: 'm'
             };

             map = new OpenLayers.Map('map', options);

             var ccounties = new OpenLayers.Layer.WMS(
                             "test", "http://localhost:8080/geoserver/wms",
                             {
                                     srs: 'EPSG:28992',
                                     layers: 'test:put_25_vlak_1_vulling',
                                     format:'image/png'
                             },
                             {singleTile: true, ratio: 1}
             );
             map.addLayer(ccounties);
             // build up all controls
             map.addControl(new OpenLayers.Control.PanZoomBar({
                     position: new OpenLayers.Pixel(2, 15)
             }));
             map.addControl(new OpenLayers.Control.Navigation());
             map.addControl(new OpenLayers.Control.Scale());
             map.addControl(new OpenLayers.Control.MousePosition());
             map.zoomToExtent(bounds);
    info = new OpenLayers.Control.WMSGetFeatureInfo({
        url: 'http://localhost:8080/geoserver/wms',
        title: 'Identify features by clicking',
        queryVisible: true,
        eventListeners: {
            getfeatureinfo: function(event) {
                map.addPopup(new OpenLayers.Popup.FramedCloud(
                    "chicken",
                    map.getLonLatFromPixel(event.xy),
                    null,
                    event.text,
                    null,
                    true
                ));
            }
        }
    });
    map.addControl(info);
    info.activate();
  </script>


texts: QtDesigner is not starting for qgis bat file

texts:
The qt-designer.bat file did not launch the QtDesigner. I have tried the file available with OSGeo4W64 installation but the qt-designer did not starts.I also try to start the qtdesigner from the QtDesginer with QGIS 3.10 custom widgets but also did not start.
I have tried this solution in this link: No custom widgets in Qt Designer for QGIS 3

Note Qgis is v3.4


texts: Creating curves using the Parcel Editor tool in ArcGIS 10.1

texts:
We are migrating parcel editing environments from AutoCAD to ArcGIS and I'm having trouble creating curves within a Parcel Fabric that are not tangent to lines.  For example, what would be some methods for creating the curve (a=9.618, r=190.000) for this image.

I have yet to find good documentation for performing this task.  Anyone know of good sources?


texts:
NASA makes available a Digital Elevation Model for the Moon in GeoTIFF format. I am trying to use QGIS to create a map of a specific region, bounded by certain values of latitude and longitude. From what I can gather (I am not a GIS person; this is my first time using QGIS), the GeoTIFF in question is encoded as x, y, elevation, but I need to work with latitude and longitude.
I load the file in question into a new project by adding a new raster layer and selecting the file. However, when I add a grid decoration, it is displayed as x, y. In addition, if I want to select a point or region, I must input x, y.
I tried extracting elevation contours (using Raster->Extraction->Contour), add that layer, and try to add fields via the Attribute Table (specifically the $x and $y attributes, as suggested in this answer). However, the computed values are NULL.
From what I can gather, my CRS looks correct: +proj=ecq +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=1737400 +b=1737400 +units=m +no_defs. It looks correct to me (again, not my wheelhouse) because it shows that the map is in a simple cylindrical projection, with lat/lon zero at xy zero, over a sphere with radius 1,737,400 m (the correct radius of the Moon).
I also tried converting the raster into a vector (using Raster->Conversion->Polygonize) and work with the resulting shapefile to no avail.
How would you go about loading the linked GeoTIFF in QGIS and displaying a map of, say, upper-left corner in lat/lon (30 deg, 0 deg) and lower-right corner in lat/lon (-30 deg, 30 deg). The map would have a grid separated in, say, 5 degree increments.


texts:
I'm using OpenLayers and OpenStreetMap tiles on my map. I'm drawing a circle somewhere in Britain. Then I call getRadius() method on this circle. Map is set to projection: epsg3857
I don't understand why it gives me that strange results. Could you please help me to understand why getRadius() function returns value that is not represented by drawn circle? As a result I suppose that point and radius will represent different circle from what I can see on the map (it is used to later on find some objects inside that circle).


texts:
I have a 3D line (e.g. flightpath) and want to create a 2D polygon buffer projected onto the ground (like a shadow), which depends on the altitude - i.e. this would look like a long thin polygon with variable width.
The naïve way would seem to be to create a series of points along the line at suitable spacing, buffer each point based on the exact interpolated altitude of the specific point and then dissolve these circles into one polygon.
However, the above approach is quite expensive if you don't want "gaps" between the circles, and I feel like there must be a more efficient way of doing this directly. Note that this is not quite the same as buffering a line via attribute, as each line segment will have a gradient (i.e. it might start lower or higher than it ends).
Can anyone provide insight to the approach here - is this a known subclass of problem? Note that I'm happy to use any GIS platform or derive something from a suitable algorithm, but not sure where to start here.


texts:
After the opening of Google Satellite through the OpenLayers plugin in QGIS, it shows the error

  "For Development process only"

How to resolve this issue?


texts: How to query two tables using ST_Within?

texts:
I am stuck trying to find out how to see if points from my parking table exists in the polygons in my community districts table. So basically the the community districts table contains a borocd(borough ID) and the_geom(polygon geometry). The parking table has a zipcode column and a thegeom(point geometry). I am trying to figure out how to see if thegeom(from parking) is within the_geom(community district number 101). So basically my community district table looks like:

boroid   the_geom
101    | 0106000020D7080000040000000103000000010000003803000000728A4...
102    | 0106332020D7080000040000000103000000010000003093000000998A4...
101    | 0107112020D7080000040000000103000000010000003803000000087A4...

And my parking table looks like:

 zipcode |                      thegeom
---------+----------------------------------------------------
  11203 | 0101000020E6100000000000C0C351444000000000817C52C0
  14634 | 0101000020E6100000000000C0D0564440000000A0597152C0
  13401 | 0101000020E6100000000000804D584440000000C0B37E52C0

The query I am trying to run is:

select a.thegeom from parking a, cd b  where b.boroid = '101' AND  ST_WITHIN(a.thegeom,b.the_geom);

After I run this query nothing is displayed.


texts: Keep the options of the cells of each column of the attribute table of a layer in QGIS when exporting the layer in a format

texts: Preserve the options of the cells of each column of the attribute table of a layer in QGIS when exporting the layer in a format or when importing the layer in PostGIS with those options in the cells of each column without leaving those options of the cap


texts: Unable to edit shp files in saved QGIS project

texts:
Yesterday I was working on a new project in QGIS. I loaded some shp files, edited them with the "edit" option and saved the project in qgs format.

I've opened the same project today and noticed that the "edit" option is greyed out.

Am I missing some option to re-enable editing?


texts: Reprojecting raster data from a rotated pole grid

texts:
I have some raster data (ASCII grid files) that use a custom rotated-pole projection, that I'd like to reproject into something more useful, such as lon/lat WGS84. A description of the grid projection is here, but I can't see how I can use this to create a projection string for use in GIS.

http://ukclimateprojections-ui.defra.gov.uk/ui/docs/grids/prob_land_25km_rotated/index.php

This SE question seems to discuss using rotated pole data, but I can't work out how to apply it to my problem: Lon/lat transformation

I found a thread here asking exactly same question (about the same dataset), but with no replies: http://lists.osgeo.org/pipermail/gdal-dev/2011-March/028174.html

The header of the ASCII grid data looks like this:

ncols 39
nrows 52
xllcorner -0.50
yllcorner -0.50
cellsize 1.00


texts:
I am trying to get the polygon of a basin boundary using R. I got a 7.5 arcsec DEM that I download from this page:  https://topotools.cr.usgs.gov/gmted_viewer/viewer.htm

I look through a lot of websites and papers where I found that it is possible but not a with a clear example of how to implemented. In this paper: http://confluence-jwsm.ca/index.php/jwsm/article/view/2/1 ; they use RSAGA package to get the boundary but I have found no function to that in the package.

Does anybody know how to do this ??


texts:
I am trying to convert feature collection GeoJSON into image. I want to use geotools to convert. I tested the soultion from this post here. But I am really noob that I ran into some errors like this:

  Exception in thread "main" java.lang.IllegalArgumentException: Width (1000) and height (0) cannot be <= 0.

When I print the features.getBounds() it produces this result:

  features ReferencedEnvelope[0.0 : -1.0, 0.0 : -1.0]
  Why is features.getBounds is not producing the correct result instead of 0? What should be the correct result?

This is the complete code:

    URL states = new URL("http://geojson.xyz/naturalearth-3.3.0/ne_110m_admin_1_states_provinces.geojson");
    FeatureJSON featureJSON = new FeatureJSON();
    FeatureCollection features = featureJSON.readFeatureCollection(states.openStream());
    MapContent mapContent = new MapContent();
    mapContent.setTitle("Quickstart");
    File outputFile = new File("states.png");
    ImageOutputStream outputImageFile = null;
    FileOutputStream fileOutputStream = null;
    try {
      fileOutputStream = new FileOutputStream(outputFile);
      outputImageFile = ImageIO.createImageOutputStream(fileOutputStream);
      int w = 1000;
      ReferencedEnvelope bounds = features.getBounds();
      int h = (int) (w * (bounds.getHeight() / bounds.getWidth()));
      BufferedImage bufferedImage = new BufferedImage(w, h, BufferedImage.TYPE_INT_ARGB);
      Graphics2D g2d = bufferedImage.createGraphics();

      mapContent.getViewport().setMatchingAspectRatio(true);

      mapContent.getViewport().setScreenArea(new Rectangle(Math.round(w), Math.round(h)));
      mapContent.getViewport().setBounds(bounds);

      g2d.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);

      Rectangle outputArea = new Rectangle(w, h);

      GTRenderer renderer = new StreamingRenderer();
      LabelCacheImpl labelCache = new LabelCacheImpl();
      Map<Object, Object> hints = renderer.getRendererHints();
      if (hints == null) {
        hints = new HashMap<>();
      }
      hints.put(StreamingRenderer.LABEL_CACHE_KEY, labelCache);
      renderer.setRendererHints(hints);
      renderer.setMapContent(mapContent);
      renderer.paint(g2d, outputArea, bounds);
      ImageIO.write(bufferedImage, "png", outputImageFile);
    } catch (IOException ex) {
      ex.printStackTrace();
    } finally {
      try {
        if (outputImageFile != null) {
          outputImageFile.flush();
          outputImageFile.close();
          fileOutputStream.flush();
          fileOutputStream.close();
        }
      } catch (IOException e) {
          log.error("cannot produce image");
      }
    }
  }

Where am I getting this wrong? I just need to convert to get image. I would be happy if there is any other solution too. I need to use Java library so geotools is the best I found.


texts:
I am trying to use the network analysis toolset to build a route. I made the the tool in model builder initially and have tried to edit it to make it a script.

When I try to build the network I get the following error:

The following is the current script:

import arcpy
from arcpy import env
from arcpy.sa import *
import os

arcpy.CheckOutExtension("Network")
arcpy.env.overwriteOutput = True

#getting geodatabase
geoBase = arcpy.GetParameterAsText(0)
arcpy.env.workspace = geoBase
arcpy.env.scratchWorkspace = geoBase
path = geoBase

#getting network database
network = arcpy.GetParameterAsText(1)
#getting stops
Stops = arcpy.GetParameterAsText(2)

Route = "Route"
Route__2_ = Route
Route__3_ = Route__2_
Routes = "Route\\Routes"
Current_Best_Route = geoBase + "\\current_bestRoute"

# Process: Build Network
arcpy.na.BuildNetwork(network)

# Process: Make Route Layer
arcpy.na.MakeRouteLayer(network, "Route", "Length", "USE_INPUT_ORDER",
                        "PRESERVE_BOTH", "NO_TIMEWINDOWS", "", "NO_UTURNS",
                        "", "NO_HIERARCHY", "", "TRUE_LINES_WITH_MEASURES", "")

# Process: Add Locations
arcpy.na.AddLocations(Route, "Stops", Stops, "Name Name #", "5000 Meters", "",
                      "Export_Output SHAPE;Trial_ND_Junctions NONE",
                      "MATCH_TO_CLOSEST", "APPEND", "NO_SNAP", "5 Meters",
                      "INCLUDE", "Export_Output #;Trial_ND_Junctions #")

# Process: Solve
arcpy.na.Solve(Route__2_, "SKIP", "TERMINATE", "", "")

# Process: Copy Features
arcpy.CopyFeatures_management(Routes, Current_Best_Route, "", "0", "0", "0")

arcpy.CheckInExtension("network")

I am not sure how to fix the error. When I run the tool without the arcpy.na.BuildNetwork I get the error saying:

  ERROR 030024: Solve returned a failure. The network dataset has no
  edge elements. It may be unbuilt. Failed to execute (Solve).

So it seems I must build it, but that won't work for some reason.


texts: Selecting granule from ASTER Global DEM V2

texts:
I am following this tutorial in order to download ASTER Global DEM V2 data:
http://www.echo.nasa.gov/reference/reverbastergdem_tutorial.htm
on http://reverb.echo.nasa.gov/reverb/

Step 7 is to "Select the shopping cart button for the granules you desire to order."
Here there is an unordered list of 22'702 granules, such as ASTGTM2_N68W154.zip, ASTGTM2_N62E087.zip and so on...

I want to get the ASTGTM2_N41W002.zip (around Barcelona). However, I cannot search for this in this unordered list of 22'702 rows, the list being downloaded and shown bit by bit as I scroll down the list.

I tried the filter feature ("More search options"). But as soon as I add a filter option, whatever it is, such as the given example of "Equatorial Crossing Longitude: -45.0", then it reports "No granules found".

I tried also the map view, which I can zoom and center to Barcelona. But once there, there is not option to choose that granule.

Step 7 of the tutorial says "Select the shopping cart button for the granules you desire to order.".

How do I select the ASTGTM2_N41W002.zip granule?


texts: Openlayers Vector Layer Features are empty even after DrawFeature executed

texts:
I am new to Openlayers and Javascript.

I am creating a web map that will allow users to draw shapes on top of a base map. I'm trying to figure out how to access the coordinates of the shapes that are drawn on the map, but the layer that I thought would contain the drawings always appears empty.

My current code:

function getFeatures() {
    var annotation_layer = map.getLayersByName("Annotations")

    window.alert(annotation_layer.features)
    for(var a = 0; a < annotation_layer.features.length; a++ ){
        window.alert(annotation_layer.features[a].geometry)
    };
}

function init() {

    map = new OpenLayers.Map('map', {
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG:4326"),
        controls: [
            new OpenLayers.Control.PanZoom(),
            new OpenLayers.Control.Navigation()
        ]
    });

    var gphy = new OpenLayers.Layer.Google(
        "Google Streets"
    );

    annotation = new OpenLayers.Layer.Vector("Annotations", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        projection: new OpenLayers.Projection("EPSG:4326"),
        protocol : new OpenLayers.Protocol.HTTP({
                url : 'polygons.geojson',
                format : new OpenLayers.Format.GeoJSON()
            })
    })

    map.addLayers([gphy, annotation]);

    var panel = new OpenLayers.Control.Panel({
        displayClass: 'customEditingToolbar',
        allowDepress: true
    });

    var draw = new OpenLayers.Control.DrawFeature(
        annotation, OpenLayers.Handler.Polygon,
        {
            title: "Draw Feature",
            displayClass: "olControlDrawFeaturePolygon",
            multi: true
        }
    );

    panel.addControls([draw]);
    map.addControl(panel);

}

The function getFeatures() is linked to a button in the html. When pressed, I would like it to create popups for the features that have been drawn.

The problem is that the features array inside the annotations vector layer always appears empty to the getFeatures() function even after several shapes have been drawn on the map. How do I get annotations vector layer to correctly contain the drawn features?


texts: Load CSV in SQLite and then create layer in Tilemill

texts:
I have a big CSV file with addresses and their lon/lat. I am not able to load that in Tilemill directly. I thought to load that CSV in SQLite Db and then create a layer with SQLite DB.
I am not able to progress on that. I could not find a tutorial where I can load the lon/lat to SQLite and create a layer. I am getting an error saying no column as 'geometry_field' type. I think, it need some SRS value to get that worked. Can someone please tell what is the issue and where to look?


texts: Converting DEM to Floating Point Raster (FLT) using QGIS?

texts:
I am looking to convert a .DEM exported from AutoCad to a floating point raster file type (.FLT) utilizing QGIS to implement in HEC-RAS Mapper.

This is accomplished easily in Arc using the Raster>Float option.  Referencing ESRI Supported Formats it does not appear the .FLT format is proprietary. It does mention  "Write—Only via Raster To Float tool, or developer code" which could suggest that it could be incorporated into QGIS.  Although, the QGIS Raster>Translate does not offer the .FLT format.

Is there open-source workaround to convering DEM to .FLT?


texts: Feature set buffer script is not working in ArcGIS Pro?

texts:
I've created a script that takes in a user selected point through a FeatureSet parameter then creates a buffer of that point. It worked fine in ArcMap but now our work is migrating everything to ArcGIS Pro. I made the necessary changes with arcpy.mp instead of arcpy.mapping, etc. but am now getting an rutime error message when I try and add my buffer layer. The script is also created using a toolbox that I cam accessing in ArcGIS Pro.

This is my code:

arcpy.env.workspace = "E:/GIS/Data/Sinks.gdb"
arcpy.env.overwriteOutput = True
arcpy.env.addOutputsToMap = False
inPts = arcpy.GetParameterAsText(0) #Empty featureSet that the user clicks
anywhere on the map to create, the schema is set as an x/y point previously
created
distance = arcpy.GetParameterAsText(1)
arcpy.CopyFeatures_management(inPts, "sourcePt_1")
arcpy.MakeFeatureLayer_management("sourcePt_1", "point1_lyr")

p = arcpy.mp.ArcGISProject("CURRENT")
m = p.listMaps("Layers")[0]

#automatically clearing out previous buffer
if arcpy.Exists("Sink_Search_Radius"):
    arcpy.Delete_management("Sink_Search_Radius")
else:
    pass
arcpy.AddMessage("Removing previous buffers.")
for lyr in m.listLayers():
    if lyr.name == "Sink_Search_Radius":
        m.removeLayer(lyr)
    arcpy.AddMessage("Sink selection buffer removed.")
else:
    pass

if arcpy.Exists("E:/GIS/Data/Sources.gdb/Source_Search_Radius"):
   arcpy.Delete_management("E:/Data//Sources.gdb/Source_Search_Radius")
else:
    pass

for lyr in m.listLayers():
    if lyr.name == "Source_Search_Radius":
        m.removeLayer(lyr)
        arcpy.AddMessage("Source selection buffer removed.")
    else:
        pass

for lyr in m.listLayers():
    if lyr.name == "nPoint":
        m.removeLayer(lyr)
        arcpy.AddMessage("Point removed.")
    else:
        pass

 #Buffering point
 arcpy.AddMessage("Drawing a search radius around your point.")
 arcpy.Buffer_analysis("point1_lyr", "Sink_Search_Radius", distance)
 Bufflyr = ("E:/GIS/Data/Sinks.gdb/Sink_Search_Radius")
 m.addDataFromPath(Bufflyr) #where I get my errror
 arcpy.ApplySymbologyFromLayer_management("Sink_Search_Radius",
 "E:/GIS/Data/Sinks/BuffFc.lyr")
 arcpy.AddMessage("Beginning to perform your selection.")


texts: Extracting primary and secondary roads from OpenStreetMap data

texts:
After importing my country roads from OSM data, I need to create a network ( roads connecting big cities ) for use with pgRouting.

For this, I need the green, blue and red ones (http://www.openstreetmap.org/#map=11/-29.8582/-51.5238) but I'm unsure how to select them in database, since OSM is not clear about how to find them.

I just try http://taginfo.openstreetmap.org/keys/highway#values using PRIMARY and SECONDARY tag values, but the result sometimes gives me the orange and yellow ones too.


texts: Error: Exception Raised on the "for row in rows" command line?

texts:
I am comparing fields in layers. I have the AADT layer and AADTAnno layer and I want it to search the fields FLAG and TFLAG. For Example: If there is a 29C9 in AADT layer in FLAG fields then there needs to be the same thing in the ANNo layer in the TFLAG field. My code stops running on my for row in FRows command line. I think its because I'm not letting it know that FLAG is the field. Please help!!

New Code!!

import arcpy, traceback

mxd = arcpy.mapping.MapDocument(r"K:\TASS\2 - GEO-DATA PROCESSING SUPPORT\MICHELLE'S WORK_ENTER NOT!!\Work Folder\Python Programming\Wayne's Tools\Abilene_Base_Map.mxd")
lstLayers = arcpy.mapping.ListLayers(mxd)

flayer = arcpy.mapping.ListLayers(mxd, "AADT")[0]
alayer = arcpy.mapping.ListLayers(mxd, "AADTAnnoLabel")[0]

FRows = arcpy.SearchCursor(flayer)
ARows = arcpy.SearchCursor(alayer)

ffields = arcpy.ListFields(flayer, "", "FLAG")
afields = arcpy.ListFields(alayer, "", "TFLAG")

FList = []
AList = []

for row in FRows:
    Fvalue = row.getValue("FLAG")
    FList.append(str(Fvalue))

for row in ARows:
    Avalue = row.getValue("TFLAG")
    AList.append(str(Avalue))

matched = set(FList) & set(AList)

for x in matched:
    exp = "ID = " + x
    arcpy.SelectLayerByAttribute_management(flayer, "ADD_TO_SELECTION", exp)
    arcpy.SelectLayerByAttribute_management(flayer, "SWTCH_SELECTION")


texts: Arrow head at the last path segment of Bezier curve

texts:
I am making a flow map. For most flows, I used the curved arrow symbology associated with lines.
For my Mediterranean flow, I have to draw a path that ends with an arrow.
I use the Bezier plugin to get curves.
Question: How to keep this path and keep only one arrow?
Issue: When I uncheck in the arrow symbology repeat the arrow for each segment, I only get a "big" curved arrow and I allow my path.


texts: Changing title on legend in order for them to be in row and not in column of QGIS print composer?

texts:
The generic order of the layers in the print composer are in a column.

I need them to be in a row.

How do I it?


texts: Unable to zoom beyond scale of 1:2,256 in QGIS?

texts:
In QGIS, whatever might be the source, I couldn't zoom beyond the scale of 1:2256.

Whether it's a raster or vector, result is the same.

Does anyone know how to fix this ?

On @Chris' advice I've uninstalled "Tile Map Scale Plugin" which fixed the problem (at least for me).


texts: How to convert negative NMEA latitude and longitude values ('S' and 'W')?

texts:
I'm getting from GPS data in the form of NMEA.  There I'm getting positive latitude and longitude values('N' and 'E').

How do I convert the negative latitude and longitude values('S' and 'W') into decimal degree format.


texts: Is it possible to put a raster layer on a shifted Chinese basemap using leaflet API?

texts:
Recently, we want to develop an application with leafleat and gaode basemap because some features of these products fullfill our basic needs. We want to overlay some raster layers caculated by spatial models by the means of adding Tilelayer or imageOverlay.

However, the trouble is that many basemaps including gaode Map have positions shifted (encoded with GCJ-02 ) in the China area. And other basemaps which do not have shifted encoding such as OSM lose many street details in the same region.

The question is how to overlay a raster (WGS84) onto an shifted basemap in leaflet?


texts: What is happening when I register a "Feature Class" with Geodatabase?

texts:
I'm registering a table w/ native geometry with the geodatabase and it is taking a long, long time.  Actually I'm not sure if ArcCatalog is frozen or crashing.  The table does have several million features (points).

The database is a SQL Server 2012 database and there is nothing happening on the database when I run SQL Profiler so whatever is going on, I think it's on the client. What is ArcCatalog (10.2) doing that is making it slow.   Is there a way to make it go faster?


texts: How to change expression used to determine symbology of feature-linked annotation?

texts: I created and adjusted a large feature-linked annotation feature class in ArcMap and realize that I symbolized the annotation classes based on an incorrect field. Is there a way to change the symbology expression so that I don't have to regenerate and move hundreds of annotations?


texts: How to count classified points in polygons?

texts: is there a way to merge the attribute table of a csv file into a the attribute table of a vector? Join is not the answer for me because in this way the plugin GroupStats that you suggested me to use for counting the elements in the table doesn't work properly; in fact I need to count accommodations for every municipality in my province and using Join feature between municipality vector and accommodation file doesn't count the number of the elements that I need.


texts:
I am trying to determine the relationship between multiple factors (slope(from dem-raster), aspect(from dem-raster), proximity to settlements (buffers-polygon)) with fire events(point data, with date attribute) spatially and temporally.

I have classified the multiple factors as a ranking between 1-10 (10 being the most likely to increase fire risk) but would like to quantitatively and spatially assess the relationship between the ranking and fire events.

I don't know which technique will provide me with what I want. Does anyone have any ideas of how I can investigate and represent this relationship in ArcGIS?


texts: Tool to delete too many points in one area, ArcGIS, ArcMap

texts:
I have two layers: points and areas (buildings).

In many places I have two or three points in one building. I would like to leave only one point in one building. How can I do that?


texts:
First of all, I am new to python and rasterio. I am trying to open a simple raster with 5 bands using the rasterio module in python.
This is what I am trying to achieve:

import rasterio

with rasterio.open('/Users/martinsiklar/Downloads/LC81330242015176LGN00/baikal_subset.tif') as src:
    print(src.width, src.height)
    print(src.count)
    print(src.indexes)
    print(src.dtypes)
    print(src.crs)
    print(src.transform)

Unfortunately I am getting the following Error and I can't really figure out what's wrong. Google didn't help either:

  Traceback (most recent call last):   File
  "/Users/martinsiklar/Martin/Python/arcpy/arcpy101", line 3, in

      with rasterio.open('/Users/martinsiklar/Downloads/LC81330242015176LGN00/baikal_subset.tif')
  as src:   File
  "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/rasterio/init.py",
  line 105, in open
      from rasterio._io import RasterReader   File "init.pxd", line 155, in init rasterio._io (rasterio/_io.c:49668)

  ValueError: numpy.dtype has the wrong size, try recompiling

    Process finished with exit code 1

When I take a look into my modules, rasterio seems to be installed properly (as far as I can tell.) The path of the .tif is also correct.  I have numpy 1.8 installed.


texts: Extracting one band from three band raster in ArcMap?

texts:
I'm fairly new to ArcGIS and I get into trouble extracting a single band layer from a raster image .img in ArcMap that I want to use after in zonal statistics. So far, all help I found is to do the operation in the other way around and merge bands and some explanation by using the python code. But I would like to use ArcMap interface.

I tried to create a single layer from ArcCatalog and import it to ArcMap but I get a layer and not a raster. I also tried to import my image.img in ArcMap but then I don't find anything to extract only a specific band.


texts: Convert land-based camera image pixel x,y to ground lat/long or x,y,z

texts:
I have a number of cameras at various known locations (x,y,z) and directions that are above ground about 20 feet pointed at different ground locations about 100 feet away from each camera, and whose furthest fields of view (FOV) are about 300 feet. None of the cameras are pointing directly N-S or E-W. Because of the relatively short distance, I don't need to take into consideration the curvature of the Earth.

I am detecting moving objects within the camera FOV and obtaining their image pixel x,y contours and centroids. How can I convert the pixel x,y locations into (initially 'flat') lat long coordinates or relative ground x,y locations? I am obtaining elevation data, so will be adding the ground height (z).

Am trying to do this in Java as much as possible, so any pointers (e.g., algorithms, code libraries, etc) are welcome


texts: Creating Route Event layer with Right Angle?

texts: I have collected GPS points along the road. and got Right angle Tape measurements of Boundary walls and other objects from that respective points. I have those records in Excel Sheet and GPS observation is in Shp file. Now I want to plot those points. I Have Tried with Linear referencing Tool in Arcmap. The events were created as point layer but that is not in Right Angle... How to Solve this?


texts: Unable to reorder after layers from multiple ArcGISDynamicMapServiceLayers

texts:
I used the ArcGIS Javascript layer reordering sample as template.

Initially I have 4 layers added to the map by ArcGISDynamicMapServiceLayer where reordering was working. Then, instead of using add lake from a table, I add 5 other layers using ArcGISDynamicMapServiceLayer again and pushed the layers into the first ArcGISDynamicMapServiceLayer with DynamicLayerInfos.push(newDynamicLayerInfo).

The drag and drop list would update to include all the old + new layers. But the new layer wouldn't reorder.

I drag the new layers all the way to the top and move it around. What ended up happen was that no matter name in the list was the first four item, it would change the order of the older layers. Example:

old1 => old layer1
old2 => old layer2
new4 => old layer3
new3 => old layer4

so if I move it around this would happen

old1 => old layer1
new3 => old layer4
new4 => old layer3
old2 => old layer2

It doesn't seem to be registering the new layer in the older ArcGISDynamicMapServiceLayer which I push the new layers into. Can anyone tell me what I did wrong?

Allen

function addProperty(){
    var propertyLayer, propertyDynLayerInfos, dynamicLayerInfo;

    // disable the "Add Property" button
    dijit.byId("property").set("disabled", true);
    dijit.byId("property").destroy();

    // show a loading icon
    //dojo.style(dojo.byId("loading"), "display", "inline-block");
    dojo.style(dojo.byId("loading"), "visibility", "visible");

    urlDyn = "http://sampleserver6.arcgisonline.com/arcgis/rest/services/USA/MapServer";
    propertyLayer = new esri.layers.ArcGISDynamicMapServiceLayer(urlDyn, { "id": "Property2" }); //only work for mapserver not featureserver
    app.map.addLayer(propertyLayer);

    var addLayerhandler = dojo.connect(
        propertyLayer,
        "onUpdateEnd",
        function() {
            dojo.disconnect(addLayerhandler);
            dojo.style(dojo.byId("loading"), "visibility", "hidden");

            // create a new dynamic layer info object for the property layer
            if ( ! app.hasOwnProperty("dynLyrInfos") ) {
                app.dynLyrInfos = app.map.getLayer("Kings County").createDynamicLayerInfosFromLayerInfos();
            }

            // add new layers
            //propertyDynLayerInfos = propertyLayer.hasOwnProperty("dynamicLayerInfos") ? propertyLayer.dynamicLayerInfos : propertyLayer.layerInfos;
            propertyDynLayerInfos = propertyLayer.createDynamicLayerInfosFromLayerInfos();
            alert(propertyLayer + " " + propertyLayer.hasOwnProperty("dynamicLayerInfos"));
            dojo.forEach(
                propertyDynLayerInfos,
                function(dynamicLayerInfo, i)
                {
                    app.dynLyrInfos.push(dynamicLayerInfo);
                }
            );

            app.map.getLayer("County").setDynamicLayerInfos(app.dynLyrInfos);
        }
    );


texts: QGIS 3.10.7 Buffering a line by 1m with flat end cap, produces "cracked" buffer

texts:
I need to buffer a line 1 meter, however when running the buffer tool in QGIS with the end cap style as "flat" I get the result shown below.
What causes this, and how can I fix it?


texts: QGIS 3.10 Populate attribute table with ellipsoidal distances along lines

texts:
I have a GIS with a logistics layer that contains line features. The lines have been digitised from aerial photography and each line represents a path from A to B, following roads, rail and waterways.

My goal is to write the real world distance along each path to a column within the attribute table.

So far I have tried using:

[Toolbox > Vector geometry > Add geometry attributes]

The results are as follows:

My result of 16.5 is confusing, because the actual distance in this case is about 1650km. I followed a tutorial to do this and their results were in metres, so I guess my questions are:

Is my result in degrees or some other units?
Is there a better way to get the ellipsoidal distance along each line into the attribute table?

(Edited to note I am using QGIS 3.10)


texts: Working with tiled datasets in QGIS 2.18?

texts:
I have datasets (OS Opendata to be precise) which are broken down into tiles, for some data as small as one square km.

I'm creating maps from these datasets and I wanted to know if there was a simpler way to manage their appearance / styles when instead of one layer for (for example) terrain, or roads, you have a hundred or more.

While grouping them allows me to turn them off and on as a single layer and navigate the list of layers more easily, it doesn't help when it comes to managing the symbology, you can't edit the symbology of a whole group and so still have to set the style for each layer individually.

Being able to copy/paste styles helps a bit, though it's still more time-consuming than I'd like, especially if you want to tweak the symbology, each time you want to adjust the colour or thickness of a line or somesuch you have to then go layer by layer pasting the style.

Another 'solution' is to merge the tiled layers into a single shapefile, however I would like to stick to non-destructive editing if at all possible. The OS Data is updated around every 6 months and I would rather just be able to download the new data, place it in the folder overwriting the old data and hey presto, the map now uses the new data, rather than have to create my new shapefiles each time the data is updated.

I may be missing something, but is there a better way of doing this, for example of setting the symbology for an entire group, or of having a layer, the definition for which is a union query so that I can include several individual layers in a single layer?


texts:
Now,I know how to insert a shapefile to PostGIS by using shp2pgsql.exe.
If I have an updated shapefile and I only want to update the updated rows, how can I know the table being updated?


texts:
2 years ago, this question was posted here at Gillies's blog,  on date of 2009-05-15.
Future for geospatial on App Engine?
At that time there was no news relating to any apps of GIS web services using Datastore (Google's App Engine).
Is there any new information on GIS and App Engine over last 2 years?
In particular, I am looking for a WFS that uses Datastore on the backend and a browser and a javascript API on the client that will allow people around the world to collaboratively make maps, which can be represented in structured way such as with geoJSON or a shapefiles.
In general, what cool map making and spatial querying web services use Google's App Engine Datastore instead of typical relational database ?


texts: Opening saved QGIS layout in layout designer with Python

texts:
Is it possible to open saved layout in layout designer(print composer) with Python?
QGIS Python API:

openLayoutDesigner(self, layout: QgsMasterLayoutInterface) → QgsLayoutDesignerInterface

Opens a new layout designer dialog for the specified layout, or brings
an already open designer window to the foreground if one is already
created for the layout.
New in version 3.0.

I have tried iface.openLayoutDesigner("layoutname") in python console, but I get error:

TypeError: QgisInterface.openLayoutDesigner(): argument 1 has
unexpected type 'str'

I don't understand what argument 1 should be.


texts: Dissolving polygons with common attribute within certain proximity to each other using QGIS?

texts:
I have a land parcel with a lot of polygons with the same attribute (owners). I would like to dissolve all the polygons with the same owners but only if the polygons are within 100m of each other. How would I go about this on a large scale (ie >10000 polygons)?

Using QGIS 1.8

The purpose is to group together large properties (Farms) into their actual operating unit rather than based on property title.


texts: how to prevent a mixture of "old" and "new" Bing Maps tiles

texts:
The Bing Maps blog states that the "new" Bing Maps style will be implemented today.

I'm seeing a strange mixture of old and new tiles when using the ArcGIS Server JavaScript API (use this sample to test, choosing the Roads option). http://www.bing.com/maps does not exhibit this problem.

Questions

1) am I just unlucky in trying this exactly as ESRI/Microsoft are changing the tiles over? (ie, does it work differently for you?)

2) if not, how can I force my users to download new Bing Maps tiles (but only once, allowing them to cache them hereafter)?

As luck would have it, I have a product launch scheduled for tomorrow so I'm hoping desperately that this problem will resolve itself overnight, Australian time...


texts: How to create for each point feature an Euclidian Distance raster in ArcGIS

texts:
I need some comment on how to create an Euclidian distance raster for a point dataset of about 100 features in ArcGIS. I need to create a raster for each single point.
My question is, if there is a tool that allows me to interpolate the already existing Eucledian distance ratser to each single point OR if I need to create a Model Builder to cue it.
In the picture you see an example. The green feature are points.


texts: How can I make geotiffs of individual depth layers in the Hycom ocean circulation models?

texts:
I am trying to pull layers from the glocal ocean circulation model Hycom.  Before I go writing code to try to do this, is there a way to pull a single z depth as a geotiff?

e.g. I have hycom_glb_911_2015093000_t000_uv3z.nc.  gdalinfo (gdal pre 2.1.0) reports:

gdalinfo hycom_glb_911_2015093000_t000_uv3z.nc

Driver: netCDF/Network Common Data Format
Files: hycom_glb_911_2015093000_t000_uv3z.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#classification_authority=not applicable
  NC_GLOBAL#classification_level=UNCLASSIFIED
  NC_GLOBAL#Conventions=CF-1.0 NAVO_netcdf_v1.0
  NC_GLOBAL#distribution_statement=Approved for public release. Distribution unlimited.
  NC_GLOBAL#downgrade_date=not applicable
  NC_GLOBAL#field_type=instantaneous
  NC_GLOBAL#history=archv2ncdf3z
  NC_GLOBAL#institution=Naval Oceanographic Office
  NC_GLOBAL#source=HYCOM archive file
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"hycom_glb_911_2015093000_t000_uv3z.nc":water_u
  SUBDATASET_1_DESC=[1x40x2001x4500] eastward_sea_water_velocity (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"hycom_glb_911_2015093000_t000_uv3z.nc":water_v
  SUBDATASET_2_DESC=[1x40x2001x4500] northward_sea_water_velocity (16-bit integer)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)

ncdump:

ncdump -h hycom_glb_911_2015093000_t000_uv3z.nc
netcdf hycom_glb_911_2015093000_t000_uv3z {
dimensions:
    lat = 2001 ;
    lon = 4500 ;
    depth = 40 ;
    time = UNLIMITED ; // (1 currently)
variables:
    double time(time) ;
        time:long_name = "Valid Time" ;
        time:units = "hours since 2000-01-01 00:00:00" ;
        time:time_origin = "2000-01-01 00:00:00" ;
        time:calendar = "gregorian" ;
        time:axis = "T" ;
        time:NAVO_code = 13 ;
    double tau(time) ;
        tau:long_name = "Tau" ;
        tau:units = "hours since analysis" ;
        tau:time_origin = "2015-09-30 00:00:00" ;
        tau:NAVO_code = 56 ;
    double depth(depth) ;
        depth:long_name = "Depth" ;
        depth:standard_name = "depth" ;
        depth:units = "m" ;
        depth:positive = "down" ;
        depth:axis = "Z" ;
        depth:NAVO_code = 5 ;
    double lat(lat) ;
        lat:long_name = "Latitude" ;
        lat:standard_name = "latitude" ;
        lat:units = "degrees_north" ;
        lat:point_spacing = "even" ;
        lat:axis = "Y" ;
        lat:NAVO_code = 1 ;
    double lon(lon) ;
        lon:long_name = "Longitude" ;
        lon:standard_name = "longitude" ;
        lon:units = "degrees_east" ;
        lon:modulo = "360 degrees" ;
        lon:axis = "X" ;
        lon:NAVO_code = 2 ;
    short water_u(time, depth, lat, lon) ;
        water_u:long_name = "Eastward Water Velocity" ;
        water_u:standard_name = "eastward_sea_water_velocity" ;
        water_u:units = "m/s" ;
        water_u:_FillValue = -30000s ;
        water_u:missing_value = -30000s ;
        water_u:scale_factor = 0.001f ;
        water_u:add_offset = 0.f ;
        water_u:NAVO_code = 17 ;
    short water_v(time, depth, lat, lon) ;
        water_v:long_name = "Northward Water Velocity" ;
        water_v:standard_name = "northward_sea_water_velocity" ;
        water_v:units = "m/s" ;
        water_v:_FillValue = -30000s ;
        water_v:missing_value = -30000s ;
        water_v:scale_factor = 0.001f ;
        water_v:add_offset = 0.f ;
        water_v:NAVO_code = 18 ;

// global attributes:
        :classification_level = "UNCLASSIFIED" ;
        :distribution_statement = "Approved for public release. Distribution unlimited." ;
        :downgrade_date = "not applicable" ;
        :classification_authority = "not applicable" ;
        :institution = "Naval Oceanographic Office" ;
        :source = "HYCOM archive file" ;
        :history = "archv2ncdf3z" ;
        :field_type = "instantaneous" ;
        :Conventions = "CF-1.0 NAVO_netcdf_v1.0" ;


texts: NetCDF with separate lat/lon bands to GeoTiff with Python/GDAL

texts:
I know this has been asked before, but none of the suggested solutions works for me, and Converting Sentinel-3 data (netcdf) to GeoTIFF  is for Java. I tried the accepted answer from Map project a raster having separate latitude and longitude raster bands, but got the error message

ERROR 1: Unable to compute a transformation between pixel/line
    and georeferenced coordinates for Oa01_radiance.tif.
    There is no affine transformation and no GCPs.

although I made sure that my latitude.vrt and longitude.vrt have the same dimensions as Oa01_radiance.vrt and also point to the correct files.

So I have some Sentinel-3 data, which come as a folder of a number of NetCDF's, one per radiance band, together with another NetCDF holding tie points for every x-th pixel (every 64th, in that case). I managed to interpolate the lats and lons and create single tifs with the same dimensions as the radiance bands, so basically the lat/lon information for each pixel (I checked them, they are correct). From the other questions around the web, I tried this:

# create vrt for data, lat and lon
gdal_translate -of VRT lat.tif latitude.vrt
gdal_translate -of VRT lon.tif longitude.vrt
gdal_translate -of VRT NETCDF:"Oa01_radiance.nc" Oa01_radiance.vrt

Then I edited data.vrt to only contain the following:

<VRTDataset rasterXSize="4865" rasterYSize="4091">
  <metadata domain="GEOLOCATION">
    <mdi key="X_DATASET">MY_FILE_PATH\longitude.vrt</mdi>
    <mdi key="X_BAND">1</mdi>
    <mdi key="Y_DATASET">MY_FILE_PATH\latitude.vrt</mdi>
    <mdi key="Y_BAND">1</mdi>
    <mdi key="PIXEL_OFFSET">0</mdi>
    <mdi key="LINE_OFFSET">0</mdi>
    <mdi key="PIXEL_STEP">1</mdi>
    <mdi key="LINE_STEP">1</mdi>
  </metadata>
  <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="0">NETCDF:Oa01_radiance.nc</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4865" RasterYSize="4091" DataType="Float32" BlockXSize="4865" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="4865" ySize="4091" />
      <DstRect xOff="0" yOff="0" xSize="4865" ySize="4091" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

and then tried

gdalwarp -geoloc -t_srs epsg:4326 Oa01_radiance.vrt Oa01_radiance.tif

which resulted in the above mentioned error. So

is there something I'm missing and
is there some Python code that can convert the data from the given structure to GeoTiff?

One more comment: using min/max of the lat/lon bands (even the interpolated ones) and calculating the GeoTransformObject from them results in unacceptable offsets and distortions of several hundred km.


texts:
I try to change the position of the title created by qgis2web to the top left above the zoom control. Right now it is always displayed on the top right side.

And is it somehow possible to make these changes just in the qgis2web template file?


texts:
I am trying to QGIS 3.16.2 and I am getting an error message saying:

SSL errors occurred accessing URL

The certificate has expired.
How can I rectify this?
Here's an image for more detail:


texts: How to create complex legends in Mapserver?

texts:
I want to create a color legend on mapserver.
It is working but now I have to group several classes together like a heading with a group of color in it and so on. I have 10 such groups of classes.

I tried using the concept of classgroup and group in a class, but now it only returns the last classgroup with no label/heading.

I want an ordered  hierarchy showing group and contained classes.

Please, it would be really helpful if somebody could help me as to how to display all the groups of classes.


texts: Focal Statistics produces corrupted/ partial Output ArcMap 10.4.1

texts:
I recently encountered an issue with ArcMap that I am not able to solve by a workaround:

I have a grid file covering the EU plus a couple of other countries inside a *.gdp (cell size 1000 x 1000 m). The file contains a population grid, meaning different values for all cells where population is situated but in areas where nobody lives, I have 0 values resulting in around 7.1 million individual cells.

I need to use the tool 'Focal Statistics' in order to calculate the population potential within a given area around each point.

Everytime I try, I receive a corrupted output (please see image attached). It does not only produce a partial output but reutrns wrong values (negative values).
I already tried to change the environmental settings (m to map units, etc.) it did not work.


texts:
I want to deploy QGIS 3.4.7.1 to several Windows 7 x64 machines. It would be nice, if I could get the installation logs somehow.
I have tried to use /V 4 as parameter, and /O to define log file, but there is no result. Can somehow please send me a document about available CLI installation parameters, or suggest me what should I use?


texts: How to sort layers within groups with Python?

texts:
After running a longer script (Qgis Console), all I want to do is sort the layers within my groups;
since some of them are buffers and the 'biggest' buffer is now on top.

I've tried stuff based on:
Post on how to sort Layers

And Martin Dobias's post concerning the tree layer api:
Layer Tree API 1 ; Layer Tree API 2 ; Layer Tree API 3

It sounds simple enough but for some reason I can not figure it out...

Among many things I've tried this:

root = QgsProject.instance().layerTreeRoot()
for child in root.children():
  if isinstance(child, QgsLayerTreeGroup):
      lyr = child.children()
      lyr.sortItems(0,Qt.AscendingOrder)

But .children() appearantly returns a list(?)

The iface.legendInterface() group-options didn't help either...


texts:
Has anyone else had any issues using the random points tool in QGIS?

When I try using it I'm selecting the option to use a value input field. The tool starts but then freezes once it gets to the 10% processing level as indicated by the status bar. Is this a known bug?


texts: Does OpenLayers 4.6.5 (or latest 5.3.0) support reprojection of VectorTiles in MVT format?

texts:
I have been setting up my own MVT-service, which takes some GeoJSON input and convert it to MVT tiles. Then I've used Open Layers 4.6.5 to render the tiles.

The testdata I have been playing with is for the Norwegian Continental Shelf.

It all works great in EPSG:3857 projection, my tiles render nicely where they are supposed to be (I'm at zoom-level 2 for overview purposes). See the image below:

However, I then try to render the same layer in WGS84 projection, and I get this display:

This leaves two questions:

The first thing is the location, ending up in Africa could mean that I have an axis problem. How do I change that?
Secondly, despite of the location being wrong, I clearly see that the features still render in EPGS:3857 projection. So they are not reprojected according to the projection of the view object.

These are the relevant code snippets:

    //Define view projection to WGS 84
var viewProj = ol.proj.get('EPSG:4326');

//Use a baselayer with vector data
var countries = new ol.layer.Vector({
      source: new ol.source.Vector({
        url: 'https://openlayers.org/en/v4.6.5/examples/data/geojson/countries-110m.geojson',
        format: new ol.format.GeoJSON({
          defaultDataProjection: ol.proj.get('EPSG:4326'),
          featureProjection: viewProj
                })
              })
        });

//Define a "custom" tilegrid so that we match the projection (and can get proper resolutions for the view later on)
var tileGrid = ol.tilegrid.createXYZ({
        extent: viewProj.getExtent(),
        tileSize: 256
    });

//Set up the Vector Tile layer with a custom tileLoadFunction
var layer = new ol.layer.VectorTile({
    source: new ol.source.VectorTile({
    attributions: '...',
    format: new ol.format.MVT(),
    tileGrid: tileGrid,
    tilePixelRatio: 16,
    tileLoadFunction: function(tile, url) {
      tile.setLoader(function() {
            Promise.all([
              fetch(url,{
              headers: {
                    "my.custom.header": "my.custom.value"
                }
              }),
            ])
            .then(function(responses) {
              return Promise.all([
                responses[0].arrayBuffer()
              ]);
            })
            .then(function(results) {
              var format = tile.getFormat();
              var features = format.readFeatures(results[0]);
              tile.setFeatures(features);
              tile.setProjection(format.readProjection(results[0]));
            });
          });
        },
    url: "http://my/service/that/returns/mvt/tiles/"+"{z}/{x}/{y}.mvt?"
    })
});

//Create the map object
var map = new ol.Map({
target: 'map',
view: new ol.View({
      center: ol.proj.fromLonLat([8.0, 65.822], viewProj),
      projection: viewProj,
      zoom: 2,
      //Set the extent and resolutions to match the defined projection/tileGrid
      extent: viewProj.getExtent(),
      resolutions: tileGrid.getResolutions()
    })
});

//add the layers
map.addLayer(countries);
map.addLayer(layer);


texts: ArcPy function works in ArcMap Python window but not from standalone script?

texts:
So I think I'm running into a data/schema lock issue, but I have no idea how to tell, or how to fix it for that matter. Here's the setup:

I have this chain of gp functions:

#extract elevations to stations
arcpy.gp.ExtractValuesToPoints_sa(station_locations, elevation_raster, scratchGDB + "\station_elevations", "NONE", "VALUE_ONLY")
#calculate the mean air temperature for each station over the n-hour time period
arcpy.Statistics_analysis(data_table, scratchGDB + "\station_statistics", "air_temperature MEAN", "site_key")
#join elevations to the statistics table created above
arcpy.JoinField_management(scratchGDB + "\station_statistics", "site_key", scratchGDB + "\station_elevations", "Site_Key", "RASTERVALU")
#select rows in "station_statistics" that have a positive elevation (non-null)
arcpy.TableSelect_analysis(scratchGDB + "\station_statistics", "stats1", "RASTERVALU IS NOT NULL")

scratchGDB is set from:

scratchGDB = arcpy.env.scratchGDB

When I run these four gp functions one at a time in the ArcMap python window, they run successfully. But when I run the debug in PyScripter on the stand-alone script, I get this error from the "arcpy.TableSelect_analysis" function:

Looking into "ERROR 000210" is what makes me think it's a data lock issue, but I'd like to know what's making it work in the python window and not work in the PyScripter IDE?


texts: Why doesn't Ordnance Survey street view raster layer appear on qgis map?

texts: I've downloaded the relevant 100 km square of OS Street View raster data and have added the relevant 5 km square to a QGIS 2.8.2 map of scale 1:1500.  The layer can be seen in the layers menu but not on the map.  I've tried adding the georeferencing data supplied by OS, but a message came up saying these are supported file types.  I've also checked the CRS and it appears correct.


texts: Does arcpy have a tool for converting file geodatabase Table to JSON?

texts:
arcpy has:

arcpy.FeaturesToJSON_conversion()

which can handle feature layers, but what about file geodatabase tables?

Not too hard to implement custom, but I want something with the same API as FeaturesToJSON and will be maintained by Esri as updates are made.  Here is a custom solution, but I don't like it:

def feature_table_to_json(tablePath, fields='*'):
    rows = []
    with arcpy.da.SearchCursor(tablePath, field_names=fields) as cursor:
        print cursor.fields
        colnames = cursor.fields
        for row in cursor:
            rows.append(dict(zip(colnames, row)))
    date_converter = lambda r: r.isoformat() if isinstance(r, datetime.datetime) else None
    result = {}
    result['features'] = rows
    return json.dumps(result, default=date_converter)


texts: GeoAlchemy2: Cannot insert row because Function st_geomfromewkt(unknown) "does not exist"

texts:
i've been stuck on this problem for quite a while now and I can't figure out what could be the issue here.

I have a Postgres DB with postgis extension installed

CREATE EXTENSION IF NOT EXISTS postgis WITH SCHEMA 'testschema'

I set the search path to 'testschema' and create the following table:

class StationGroup(Base):
    __tablename__ = 'station_Group'

    id = Column(Integer, primary_key=True)
    name = Column(String)
    center_lat = Column(Float)
    center_lon = Column(Float)
    radius = Column(Float)
    geom = Column(Geometry(geometry_type='POINTZ', dimension=3, srid=4326))

    def __init__(self, name: str, center: Coordinate):
        self.name = name
        self.center_lat, self.center_lon = center.lat(), center.lng()
        self.radius = 0.0
        self.geom = 'POINTZ(%f %f %f)' % (
            center.lng(), center.lat(), center.alt()
        )

Now I create the object and want to add it to the table with

session.add(item)

Which results in the following query:

2019-01-15 14:04:26,495 INFO sqlalchemy.engine.base.Engine INSERT INTO "testschema"."station_Group" (name, center_lat, center_lon, radius, geom) VALUES (%(name)s, %(center_lat)s, %(center_lon)s, %(radius)s, ST_GeomFromEWKT(%(geom)s)) RETURNING "testschema"."station_Group".id
2019-01-15 14:04:26,496 INFO sqlalchemy.engine.base.Engine {'name': 'name', 'center_lat': 38.3389025, 'center_lon': 15.23485, 'radius': 0.0, 'geom': 'POINT(15.23485 38.3389025 0.000000)'}

Now I get the following result:

  sqlalchemy.exc.ProgrammingError: (psycopg2.ProgrammingError) FEHLER:  Function st_geomfromewkt(unknown) does not exist

BUT:
When I execute the query manually with

INSERT INTO "testschema"."station_Group"  (name, center_lat, center_lon, radius, geom) VALUES ('name', 38.3389025, 15.23485, 0.0, st_geogfromtext ('POINT(15.23485 38.3389025 0.000000)')) RETURNING "testschema"."station_Group".id

it works perfectly!

I can't figure out what I'm doing wrong, any ideas?


texts: Debugging error 000539 from ArcGIS Field Calculator?

texts:
What am I doing wrong in this field calculator logic?

Pre-logic

def re_score(elev_t,f2, contour):
 f2=contour
 if elev_t.endswith('.5'):
  f2 = "Midi"
 return f2

Expression for field "Contourtyp"

 re_score(!Elev_t!,f2,!Contourtyp!)

I just want to rename the contour type where the contours are at half meter intervals. The 'Elev_T' is a text column with all the elevations (converted to text as I couldn't think of a way to do it with numerical elevation data which has .25, .5, 75 and 1m intervals). I want all .5 called 'midi' and .25/.75 remain as 'minor' and 1m remains as 'basic'

  Executing (Calculate Field (3)): CalculateField 328000_8123000
  Contourtyp " re_score(!Elev_t!,f2,!Contourtyp!)" PYTHON_9.3 "def
  re_score(elev_t,f2, contour):\n f2=contour\n if
  elev_t.endswith('.5'):\n  f2 = "Midi"\n return f2" Start Time: Thu Apr
  12 17:10:40 2012 ERROR 000539: :
  unexpected indent (, line 1) Failed to execute (Calculate
  Field (3)). Failed at Thu Apr 12 17:10:40 2012 (Elapsed Time: 0.00
  seconds)

I have tried different indentation and no luck...


texts: Calculating KML Lengths from Google Earth in QGIS - results in different values?

texts:
I've been trying to use QGIS to calculate the length of polylines found in KML files, after having difficulty getting ArcGIS to work properly.
I used this process:
https://www.qgistutorials.com/en/docs/3/calculating_line_lengths.html
Essentially doing:

Starting default project in QGIS (using EPSG:4326)
Using QuickMapServices to set a base map layer (say, google hybrid)
Dragging and dropping the KML into QGIS. The placement of the KML aligns with the basemap, as well as the KML's placement in other programs like google earth.
Either using the tutorial listed above (creating a geometry attributes layer via 'ellipsoidal' CRS) or using the field calculator to create a length attributes table.

Then using the lengths as I please.
However, there seems to be a difference in value between google earth pro and QGIS in some cases. In a test KML file I made in google earth pro, the length of a 15km line is almost exactly the same value (difference of about 10 metres) but with other files, the difference is almost double, despite the fact that they look exactly the same in KML
What could be causing this difference? Is it something to do with the KML files themselves (could the imported files be interpreting the lines as rectangles)? The process of importing and/or calculating? Something else?
EDIT: Ok, something is up - it seems that only ONE of the kml files i have been using previously are acting strangely. Not sure what's causing this specific one to have massive differences in distance, but the others I have used before all seem to mostly match the distances that QGIS report.


texts:
Recently I've been using QGIS to create fracture tracemaps as line shapefiles. I need to export these lines as an SVG that uses the tag polylines rather than paths. I've used the SimpleSVG plugin but this seems to export as paths.
I've also tried importing this into Adobe Illustrator and exporting again but this still only recognizes paths.

Does anyone have any advice on how to handle this?

I'm very new to any kind of coding so I'd preferably like to not have to manually edit the SVG file myself, and don't have enough experience to create code to edit the files.


texts:
i made changes in the render order of my map canvas using a :

qgis.utils.iface.mapCanvas().setLayerSet(a)

"a" being a list of QgsMapCanvasLayer

My code is working as expected (i see result on map canvas) but the "layer render order" widget (dunno the right english name...) is not refreshing ...

How can i programmatically do that ?

Thanks !


texts: How can I extract only Polygons from Google Earth Engine's Feature Collection?

texts:
I have here a shapefile which contains some limits of productive areas from the company I work for.
When I open it with ArcMap it seems there are only polygons, and indeed the source tab from layer properties shows it as I said (Geometry Type: Polygon).
In the other hand, when I upload the shapefile (with all possibles complements - shp, shx, sbn, etc), I see there are a bunch of features (about 13 among thousands of areas) which is written as Linestring.
Why does it occur and is there any way of selecting only the Features written as Polygon?
PS: Maybe my script is failing in filtering bounds based on this Feature Collection because of this issue, because when I try to use filterBounds() method with other Feature Collection (with no Linestring type) it works well, but in this particular case it doesn't.


texts: Join text file with lat/long in a shapefile without a field in common with QGIS 1.8

texts:
I'm using QGIS 1.8. I have 3 files: a shapefile from a state, a CSV file with a common field with this shape file  and a third file (CSV, with long/lat columns).

In this third file. I have a column 'id_customer'. The shapefile is a map of a state. And the file with a common column with the shapefile, has the income for various neighbourhoods. I want to link the customer file with the shape, and after link the customer to the income field.

Can anyone help me?


texts: Using lucene query parser syntax in GeoNetwork query

texts:
On my GeoNetwork instance I can query search results by using the qi endpoint:

  https://catalog-url/geonetwork/srv/eng/q?any=temperature+wms&fast=index&sortOrder=&resultType=details&sortBy=relevance&from=1&to=10&facet.q=

I wonder, how I can use the lucene query parser syntax for customizing search results?

It is possible to use the ? or * operators for wildcard searches, but the other operators does not seem to work as expected.

If I just search for e.g. 'temperature wms', it returns all results that contain both, temperature and wms keywords. The lucene website on the contrary says, this would work as the OR operator.

Also, the AND operator leads to no search results in my case, where it should give the result from the 'temperature wms' example instead (i.e. only datasets containing both keywords), according to the lucene page.

Is there a way to use the lucene query syntax with this endpoint at all? Or any other way to customize search for datasets by using advanced query syntax?

UPDATE:

I tried to dig a bit into the search classes of GeoNetwork. According to LuceneQueryBuilder, boolean AND and OR operators should actually be supported. Not yet sure, why it does not take effect in the q-interface search.


texts: QGIS: Intersect rasterdata with vectorline

texts:
I am new to QGIS.
My aim is to be able to get information of the height at any point on my vector line. My raster is a map with height information in the format geotiff. The vector layer is in the shape format.

I tried to intersect them but it wouldn't work, I guess because they are not the same format. I also tried to vectorize the raster layer and use the plugin refFunction and geomwithin but it wouldn't work either.


texts: gdalcontour output not aligned with the DEM

texts:
I use gdalcontour to generate the contour lines of a DEM in shapefile format. However the shapefile is not aligned with the DEM. Someone have a solution for this? GDAL version is 1.10

im using thins command:

gdal_contour -a height 012_merged.tif 012_contour.shp -i 25.0

gdalinfo:

Driver: GTiff/GeoTIFF
Files: 012_merged.tif
       012_merged.tif.aux.xml
Size is 9601, 4801
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,0,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-7124459.006549800700000,6446293.880868392100000)
Pixel Size = (23.191010033770905,-35.712191611463290)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-7124459.007, 6446293.881) ( 64d 0' 0.37"W, 50d 0' 0.37"N)
Lower Left  (-7124459.007, 6274839.649) ( 64d 0' 0.37"W, 48d59'59.54"N)
Upper Right (-6901802.119, 6446293.881) ( 61d59'59.80"W, 50d 0' 0.37"N)
Lower Right (-6901802.119, 6274839.649) ( 61d59'59.80"W, 48d59'59.54"N)
Center      (-7013130.563, 6360566.765) ( 63d 0' 0.09"W, 49d30' 9.16"N)
Band 1 Block=9601x110 Type=Int32, ColorInterp=Gray
  Description = Mosaicked Band
  Min=1.000 Max=305.000
  Minimum=1.000, Maximum=305.000, Mean=66.113, StdDev=78.056
  NoData Value=-32767
  Metadata:
    STATISTICS_MAXIMUM=305
    STATISTICS_MEAN=66.11288340405
    STATISTICS_MINIMUM=1
    STATISTICS_STDDEV=78.056016483986


texts: Collapsing sublayers in legend in Print layout

texts:
Does anyone know how can I collapse the legend items (subgroups) in the printing layout in QGIS?

I found some link here: https://docs.qgis.org/3.4/en/docs/user_manual/print_composer/composer_items/composer_legend.html
but nothing is about it, what I want. Is it possible at all?


texts: QGIS view bounding box on Add vector layer GeoGSON protocol?

texts:
Our server can serve GeoJSON with an optional bounding box HTTP parameter. e.g. http://our_gis_server.org.au/layername?bbox=xmin,ymin,xmax,ymax

In QGIS I can add this vector layer as a Protocol and point directly at the server. Is it possible to pass the current view bounding box dynamically as part of this query to reduce load on the server? I.E. can QGIS pass in the map bounding box as parameters to the protocol call?


texts: Feature table: Mapnik, GeoServer, CKAN, GeoNode, GeoNetwork?

texts:
I'm trying to get a clear view of the roles of different components here, this is not a "which is better" comparison, these components have different roles, but their functions overlap.

(Comparing different open source GIS servers? would be an answer if it included GeoNode, CKAN and Geonetwork, it links to this excellent table: http://www.geotests.net/cours/sigma/webmapping/2013/tableau6_webmapping2013.pdf)

So, can anyone help me fill in the holes / correct this table:

(source: greygreen.org)
https://gist.github.com/tbnorth/3be5b91fe10243898017

Y = Yes, N = No, < = handled by included subcomponent

"Catalog data": provide searchable web / web-map front end view.

GeoNode and GeoNetwork use GeoServer for rendering and nominally WMS etc.

CKAN uses pycsw for metadata federation.


texts: Can I identify if a point is inside a volume in PostGIS?

texts:
We are currently investigating PostGIS as an option to look up three-dimensional blocks which contain a point (and later also lines). At this point in time I cannot see any functionality that operates on the volume of a three dimensional shape, all functions PostGIS offers seem to operate on the surface.

Here are a number of functions we tried:

SELECT
    ST_Contains(p , mpoly),      ST_Contains(p , ph),      ST_Contains(p , box),
    ST_3DIntersects(p , mpoly),  ST_3DIntersects(p , ph),  ST_3DIntersects(p , box),
    ST_3DDWithin(p , mpoly, 0),  ST_3DDWithin(p , ph, 0),  ST_3DDWithin(p , box, 0)
FROM
  (SELECT
    'POINT Z (0.5 0.5 0.5)'::geometry AS p,
    'MULTIPOLYGON Z (
        ((0 0 0, 0 1 0, 1 1 0, 0 1 0, 0 0 0)),
        ((0 0 1, 0 1 1, 1 1 1, 0 1 1, 0 0 1)),
        ((0 0 0, 0 0 1, 0 1 1, 0 0 1, 0 0 0)),
        ((1 0 0, 1 0 1, 1 1 1, 1 0 1, 1 0 0)),
        ((0 0 0, 0 0 1, 1 0 1, 1 0 0, 0 0 0)),
        ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0))
    )'::geometry AS mpoly,
    'POLYHEDRALSURFACE (
        ((0 0 0, 0 1 0, 1 1 0, 0 1 0, 0 0 0)),
        ((0 0 1, 0 1 1, 1 1 1, 0 1 1, 0 0 1)),
        ((0 0 0, 0 0 1, 0 1 1, 0 0 1, 0 0 0)),
        ((1 0 0, 1 0 1, 1 1 1, 1 0 1, 1 0 0)),
        ((0 0 0, 0 0 1, 1 0 1, 1 0 0, 0 0 0)),
        ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0))
    )'::geometry AS ph,
    ST_3DMakeBox('POINT Z (0 0 0)'::geometry, 'POINT Z (1 1 1)'::geometry) AS box
  ) AS foo
;

All of these return false. Given that we are using polygons and polyhedral surfaces this seems correct, but is there a way to have PostGIS apply volumetric calculations?


texts: eval and cats is not recognized as an internal or external command, operable program or batch file in GRASS GIS

texts:
I'm trying to do the r.lfp and r.accumulate examples for multiple longest flow paths, except I keep getting the "eval" or "cats" "is not recognized as an internal or external command, operable program or batch file.". Not sure what I'm doing wrong, not used to working with GRASS I usually work on QGIS I just didn't find a way do get lfp from QGIS. Using GRASS on windows.
Here are the two codes I'm trying to put in GRASS:
For r.lfp:
# get nsres
eval `r.info -g map=elevation`

# delineate streams using a threshold
r.stream.extract elevation=elevation threshold=50000 stream_vector=streams

# populate stream lengths
v.db.addcolumn map=streams columns="length real"
v.to.db map=streams option=length columns=length

# create points along the streams starting from downstream
v.to.points -r input=streams output=stream_points dmax=$nsres

# find outlets (downstream-most less nsres points)
cats=`db.select -c sql="select stream_points_2.cat from stream_points_2 \
    inner join stream_points_1 on stream_points_1.cat = stream_points_2.lcat \
    where length-along > 0.5*$nsres and length-along < 1.5*$nsres"`
cats=`echo $cats | tr " " ,`
v.extract input=stream_points layer=2 cats=$cats output=stream_outlets

# create the longest flow paths for all outlets
r.lfp input=drain_directions output=lfp id_column=id outlet=stream_outlets layer=2 outlet_id_column=lcat

And this one for r.accumulate:
# set computational region
g.region raster=elevation -p

# calculate drainage directions using r.watershed
r.watershed elevation=elevation drainage=drain_directions -s -a

# get nsres
eval `r.info -g map=elevation`

# delineate streams using a threshold
r.accumulate direction=drain_directions threshold=50000 stream=streams

# populate stream lengths
v.db.addtable map=streams
v.to.db map=streams option=length columns=length

# create points along the streams starting from downstream
v.to.points -r input=streams output=stream_points dmax=$nsres

# find outlets (downstream-most less nsres points)
cats=`db.select -c sql="select stream_points_2.cat from stream_points_2 \
    inner join stream_points_1 on stream_points_1.cat = stream_points_2.lcat \
    where length-along > 0.5*$nsres and length-along < 1.5*$nsres"`
cats=`echo $cats | tr " " ,`
v.extract input=stream_points layer=2 cats=$cats output=stream_outlets

# calculate longest flow paths and delineate subwatersheds for all outlets
r.accumulate direction=drain_directions lfp=lfp id_column=id \
    outlet=stream_outlets outlet_layer=2 outlet_id_column=lcat \
    subwatershed=subwatersheds

# convert subwatersheds to vector
r.to.vect input=subwatersheds type=area output=subwatersheds

Not sure what I'm missing. Trying to execute this code from the GRASS Console tab maybe I need to do it from somewhere else?


texts: apache vs lighttpd vs nginX as OpenLayers / GIS Server

texts:
How is using one of the above technologies to host an Openlayers-based GIS application?
My investigation with apache and lighttpd showed that apache is not the best for heavy load environments.

In some cases there have been suggestions to try out the ngin-X route. Has anyone tried the three against each other and made a benchmark or a test?

Thanks!


texts: Esri ArcObjects SDK - ArcMap Desktop Add-In - Upgrade from .NET 3.5 to .NET 4.5 (or 4.5.1)

texts:
We have an ArcMap Desktop add-in, created a few years ago using ArcObjects in .NET framework 3.5. I'm looking into updating it to framework 4.5 (or 4.5.1) to be compatible with a newer project that it needs to use now.

I've bumped it up to 4.5 in a test copy and I'm not seeing any errors or erratic behavior (in the few minutes I tested it) but I'm wondering if anyone has been through this process and has had any issues down the road or has any suggestions on what I should be on the lookout for.

Also, when I bump the test copy up to 4.5.1 I do get the following error. I'll research this one myself as well, but in case this is related the experience anyone else has had so I'm adding it in.

Error:

  Error 2   The "ValidateAddInXMLTask" task failed unexpectedly.
  System.FormatException: Input string was not in a correct format.
  at System.Number.ParseDouble(String value, NumberStyles options,
  NumberFormatInfo numfmt)    at
  ESRI.ArcGIS.AddIns.SDK.Utilities.ToAddInCLRLanguage(String
  strTargetFramework)    at
  ESRI.ArcGIS.AddIns.MSBuild.ValidateAddInXMLTask.Execute()    at
  Microsoft.Build.BackEnd.TaskExecutionHost.Microsoft.Build.BackEnd.ITaskExecutionHost.Execute()
  at
  Microsoft.Build.BackEnd.TaskBuilder.d__26.MoveNext() NDOTGeolocator


texts:
Is there any API or web service to create statistical world maps where the values are different shades of a colour?

Something like http://www.visualizingeconomics.com/wp-content/uploads/sachs-gdpdensity.png


texts: I am wondering is it possible to set up two google earths? I have a project that contains all the jobs our company done and I want this to act as an aid to estimators. It is a kml file. The problem is that this project will be forgetting about if it is used in the estimators everyday google earth as it will get lost in the mess of all there other kml files they open. Can a separate google earth be set up so that only my kml will be located there and nothing else?


texts:
I use OpenStreetMap along with OSMand on my mobile device.

I have personal points of interest on my map that I do not want to upload to the community. Currently I use Favorites in OSMand to save them.

Now I want to make a print of that map with my Favorites. How can I bring my Favorites to be displayed in OpenStreetMap on my PC?

Related view-waypoints-in-osm


texts: How can I merge geotiffs with a vrt file?

texts:
I have heard that I can use a vrt file to merge mulitple geotiffs into one?

what is the syntax?  gdal_tranlsate -of VRT file1.tif file2.tif file2.tif

Something like this?


texts: Find distance between point layer and lines layer in QGIS?

texts:
All the answers were for one point and one line, But I'm talking about one point and multiple lines, and for older versions of QGIS.

I have two shapefiles including one point layer and multiple lines layer about 5000 lines shapefile.

Does anybody know a solution to this?

I'm not familiar with PyQGIS and Python and I have "QGIS 3.0.1 Girona".


texts:
I have a need to cut tiles from a raster with bounds that result in no partially filled tiles at a particular zoom level (no transparency or black/white areas).  First thought was to overlay an index grid so I could see how to cut the raster.  Is it possible to find shapefile or other description of the tiling grid for google maps or other tiling schemes.

Also interested in other idea for achieving this goal.  One is to simply programmatically discard tiles with empty areas after the fact.    Interested in other suggestions as well.


texts: Unable to intersect a shapefile with a set of points

texts:
I have a set of 1000 urban points that I'm trying to intersect with a classification result to find the accuracy but I keep getting an error
"GEOS geoprocessing error: One or more input features have invalid geometry."


texts: Crop and reproject shapefile based on GeoTIFF sample

texts:
I have GeoTIFF with these attributes:

Size is 495, 309
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6370000,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",4],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
PROJ.4 string is:
'+proj=lcc +lat_1=47.5 +lat_0=47.5 +lon_0=4 +k_0=1 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs '
Origin = (-2963997.870569999800000,-1848004.200800000000000)
Pixel Size = (11975.754948303031000,11961.191901488673000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-2963997.871,-1848004.201) ( 24d36'23.04"W, 26d22' 5.88"N)
Lower Left  (-2963997.871, 1848004.097) ( 45d39'38.75"W, 55d16'35.82"N)
Upper Right ( 2964000.829,-1848004.201) ( 32d36'23.13"E, 26d22' 5.85"N)
Lower Right ( 2964000.829, 1848004.097) ( 53d39'38.88"E, 55d16'35.76"N)
Center      (   1.4791350,  -0.0520200) (  4d 0' 0.07"E, 47d30' 0.00"N)
Band 1 Block=495x4 Type=Float32, ColorInterp=Gray
  NoData Value=9.969209968386869e+036

and I want to overlay Natural Earth shapefile.

So I have to crop and reproject global Natural Earth shapefile to this GeoTIFF sample, however I can't understand how to do this with such projection and was wondering if I can instruct somehow ogr2ogr to use my GeoTIFF as reference to crop and reproject the shapefile?

Also I'm interested in other approaches if ogr2ogr does not have this feature.


texts: What's the easiest way to generate many attribute-based depictions of the same layer?

texts:
I'm trying to show (in QGIS) how the distribution of 7 or 8 different attributes of a given  shapefile vary across polygons within the layer.

For a given attribute, this is simple (well, the basics are): use a categorized or graduated style to show which polygons are high- and which are low-intensity in the attribute.

Because I've no formal background in GIS-ing, I'd like to know what the most common/straightforward way to do this several times.

Right now my approach is to bring in 7-8 copies of the layer into the ?layer space?, and repeat the process on each of the layer copies--so that in the end I have 7 (identical except for aesthetic style) copies of the original layer, each colored according to one attribute.

Then, I bring them all in to the print composer one-by-one and use subtitles to distinguish them.

This seems somewhat inefficient to me, so I wanted to reach out to SE to see if there are any recommendations for a speedier approach to scaling up the process of illustrating attribute distributions in a shapefile.

Compounding the issue is that the base layer doesn't have the attributes--I'm bringing them in via an attribute join, so I have to re-join the table to each layer copy.

I thought of re-saving the layer once the table's joined, but I'm hoping there's a better way--the table being joined is still sort of dynamic (early in development) so it's likely to evolve (if slightly) as analysis progresses. This is a slight improvement on my current approach (only re-save the layer once if the table updates, as opposed to re-join the table 7 times), but I'm on the look out for something more drastic.


texts:
In this example from the ArcGIS for Flex samples page:
Clip and Ship

I am curious how from this function how the label property of each check box that is passed as input for the Layers_to_Clip parameter of the GPService are interpreted by the GPService as actual pieces of data.

private function getLayersToClip():Array
        {
            var layersToClip:Array = [];

            if (incidentPointsCheckBox.selected)
            {
                layersToClip.push(incidentPointsCheckBox.label);
            }
            if (incidentLinesCheckBox.selected)
            {
                layersToClip.push(incidentLinesCheckBox.label);
            }
            if (incidentAreaCheckBox.selected)
            {
                layersToClip.push(incidentAreaCheckBox.label);
            }

            return layersToClip;
        }

I have looked at the Clip and Ship data from the ArcTutor files and don't understand how the GPString values - which in this case are "Input Incident Points", "Input Incident Polylines", and "Input Incident Polygons" - correspond to any of the map layers.

Part 2) Ideally I would like to change my Flex application to take items selected from checkboxes and pass them into my model in the same fashion, except I won't be extracting them but using them in analysis. If preprocessing is done, every input would be a raster dataset regardless of the original geometry and the output would be a raster as well. How can I have those checkboxes refer to selected input rasters in a database location or within the map document/gp service?

Thanks in advance, I've been banging my head against the wall with this one,

Sean


texts:
How can data be organized in tabs within Qfield?
The "Advanced Bee Farming Demo" that comes with QField displays data organized in tabs (see picture below: General, Picture, Issues, Review, Consumption). I attempted to recreate this, but I have been unsuccessful. I can only manage to get my data to be displayed in a liner list.
What configurations need to be set in order to display tabbed data?


texts: Conflicting results for spatial autocorrelation tests

texts:
I am working with an imbalanced (5%/95%) presence/absence dataset. I've
created 25x25 m. raster cells that are categorized '1' or '0'. My Geary's C
result (~0.993) and Moran's I result (~0.0045) are each nearly ideal- the
expectation for each, respectively, would be 1 and 0 in the absence of
spatial autocorrelation. Previous research into these results tells me that this is pretty darn close to random enough.

I also created a spherical semivariogram which displays pretty constant
semivariance across all distances, save for a few outliers near the nugget
and sparsely dotted across a few other sections. These outliers near the
nugget do suggest there may be some autocorrelation present, but I'm really
not sure- we're talking perhaps 2-3 dozen dots in relation to about 10,000
others.

My question is, why would my Moran's I and Geary's C statistics indicate
very strongly against extant autocorrelation if there truly is some
remaining in the data? Should I weigh each result (semivariogram, MI, GC)
equally, or trust one result more than the others?


texts: Exploding multipart polygons within shapefile in QGIS?

texts:
I've got a shapefile with polygons, some of which are multipart polygons.  I need to create individual polygons out of the multipart polygons.

Is there a tool or plugin that can do this in QGIS 1.8?


texts: Rasterize (R package raster) fails to rasterize island polygons?

texts:
I am having problem in rasterizing a polygon shapefile: "island" polygons (I mean polygons which are totally included into a larger one) are ignored and do not appear in the raster. My code is:

seg <- readOGR(dsn="my path", layer="seg_10_9_5s")

map <- raster(map2007.tif) # import of raster of same extension than shape
# the shape has only a field with a unique id for each polygon
length(unique(seg@data$id))  # give me the right number of polygons (1568)
print(seg@data)  # everything seem OK
nrow(seg@data)   # 1568
idpol <- unique(seg@data$id)  # give me a vector of 1568 id values

# rasterize the polygons
segr <- rasterize(seg,map,"id")
writeRaster(segr, filename="seg.tif", mat="GTiff", overwrite=TRUE,datatype='INT2S')

# I  look for the values of grid cells
freq <- freq(segr)  # freq has only 1492 rows!
f <- freq[,1]
f
# I compare the list of value in vector and raster
setdiff(idpol,f)  # id of polygons which are absent from the raster

In QGIS I display the shape and seg.tif, the polygons identified as lacking in raster were not rasterized, the cells get the value of the surrounding polygon.

Maybe the problem is with the original shapefile? but it seems OK.

I put the shape and the map in
https://dl.dropboxusercontent.com/u/22467480/misterio.zip

Can anyone see what I'm missing here?


texts: How to export from ArcMap to tab-delimited text file, while preserving formatting and values?

texts:
I'm working with a complicated fGDB featureclass which contains hundreds of fields, in a variety of text and numerical formats. This is proving difficult to export into tab-delimited text file format.

The problem

After opening the table in ArcMap:

Export to DBF results in the numbers being handled correctly, but the field names being truncated after 8 characters
Export to Text results in the numbers being placed in quotations with commas separating the thousands (eg "5,181,613") and commas (not tabs) separating the values. This makes it difficult to replace the commas with tabs

My current workaround

export the table to DBF in ArcMap
export it again to CSV
edit the CSV file and remove the values, retaining just the headers
replace the commas in the headers with tabs
open the DBF file in Excel
export from Excel to a new tab-delimited text file (which has truncated headers)
copy/paste the headers from the edited headers file

Please tell me there's an easier way!


texts: Query similar Linestrings (running workouts) using PostGIS

texts:
I track my running workouts using GPS and I have the coordinates for each of them

I create a table to save each workout

CREATE TABLE runnings (id SERIAL PRIMARY KEY,
                       userId integer,
                       activityId integer,
                       geom geometry(LINESTRING,4326)
                      );

I inserted two similar workouts (workout 1 with red, workout 2 with green)

I want to find a way to query workouts that are similar to each other. As you can see, the two workouts (red, green) haven't exact points but they are on the same road and has the same direction. I want to find the matching workouts and do further analysis on these matched running workouts.

UPDATE

I did some tests with the ST_HausdorffDistance function. The task is to scan the table and find similar workouts using the userId. Is there any better way than the following query? It gets 40 seconds to finish. Now the table has only 70 rows in total. I have created a spatial index for the geometry column.

SELECT a.activityId, b.activityId, ST_HausdorffDistance(a.geom, b.geom) as dist
FROM public.runnings a, public.runnings b
WHERE a.id < b.id
AND a.userId = '123'
AND b.userId = '123'
AND ST_DWithin(a.geom, b.geom, 1.0)
AND ST_DWithin(ST_StartPoint(a.geom)::geography, ST_StartPoint(b.geom)::geography, 200.0)
AND ST_DWithin(ST_EndPoint(a.geom)::geography, ST_EndPoint(b.geom)::geography, 200.0)
AND ST_HausdorffDistance(a.geom, b.geom) < 0.001


texts: gdalinfo documentation includes -json flag but using it returns "FAILURE: Unknown option name '-json'"

texts:
Reading the documentation for gdalinfo describes the option of formatting the output as JSON.

  -json

    Display the output in json format.

However, when I use the -json tag it appears not to be implemented.

  > gdalinfo -json translated_image_overlay.png
    Usage: gdalinfo [--help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]
                              [-norat] [-noct] [-nofl] [-checksum] [-proj4]
                              [-listmdd] [-mdd domain|all]*
                              [-sd subdataset] datasetname

  FAILURE: Unknown option name '-json'

I have GDAL 1.11.1 installed.

Can someone explain why this is happening?


texts:
I am working on land-use classification and I would like to merge some classes from my training dataset. I've a FeatureCollection with over a hundred polygons (features), with all of them having a property 'class' to them with a number from 1 to 17 representing the different classes. I would like to replace from all Features the property which have a class 8 to be replaced by class 9.

I've been trying different things, and I think it should look something like this:

var conditional = function(image){
  return ee.Algorithms.If(ee.Number(image.get('class').eq(8)),
  image.set({class: 9}),
  image)

var test = table.map(conditional)

In which:

table contains my FeatureCollection with polygons
class is the property on which the conditional function must be based
8 and 9 are linked to the 'class' properties


texts: How to protect the loaded json file name and location using html5 and javascript?

texts:
The data of markers (not the map) is stored in a json/csv/js file. We want to show all data as markers on the map. But how to protect the raw data file from public access?

Following the example on leaflet website, the name and location of the data file is currently visible as plain text in the index.html. Any people can download them after viewing the source code of the html.

Is there any way to secure our raw data from public access?

for example, by encoding the file name and location and there is no way to decode it?

Or use some authorization method, tokenization, oAuth, RestFul API?

Or using PHP instead of html?

Any tips and suggestions are greatly welcome!


texts: False positive results while mapping water bodies

texts:
I am a veterinary epidemiology researcher interested to benefit livestock in a desert. As water is an important resource for herders therefore I have become interested to map water bodies.I learned GEE and used sentinel 2 data . After pansharpening , I calculated water indices and found ranges of e.g., mndwi for some water bodies  to decide a threshold range. To some extent I was successful but I get false positive results toward periphery of desert probably at inference between desert and agriculture or green land is not well remarked  (shapefile of desert is a bit rough).
Water bodies (tobas in local language) in desert are are often small about 20 m*20m , rain-fed , muddy to light green in colour , and shallow . Waterbodies are located in compact flat soil areas (dahar in local language )total barren . There is some green vegetation around around water bodies which is source of palatable forage for livestock (like ring of green around water)
What could be possible reasons for false positive results ? Any how to overcome this problem.
I burned out but attempted object based image analysis and even supervised classification too
https://code.earthengine.google.com/ce6eae289b72af1fb96e365afbad720d


texts: What does the values of Watershed Accumulation, Basin, In QGIS indicates?

texts: While performing watershed module in QGIS each output gives some value range while changing color palette. what does it actually means in accumulation , basin, stream flow. What does the value actually means?


texts: Slow attribute table loading after join table via custom plugin

texts:
I am new to QGIS and I am currently learning pyqgis.
My problem is I tried to create a custom plugin to join csv and shapefile with 20,000+ features. I was successful in joining csv and shapefile via custom plugin that I created using this code.

from PyQt4.QtCore import *
import qgis
from qgis.core import QgsMapLayerRegistry, QgsVectorJoinInfo

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    # Change to your shapefile name
    if layer.name() == "MyShapefile":
        qgis.utils.iface.setActiveLayer(layer)
        shp = qgis.utils.iface.activeLayer()

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    # Change to your csv name
    if layer.name() == "MyDatabase":
        qgis.utils.iface.setActiveLayer(layer)
        csv = qgis.utils.iface.activeLayer()

# Set up join parameters
shpField='PIN'
csvField='UPIN'
joinObject = QgsVectorJoinInfo()
joinObject.joinLayerId = csv.id()
joinObject.joinFieldName = csvField
joinObject.targetFieldName = shpField
shp.addJoin(joinObject)

But when i try to open the attribute table it loads very slow and counting features by 8/10/15 at max. But when I join table using properties>join of shapefile and check the "cache data in virtual memory" my attributes loads perfectly fine even if I have 20,000 features.

My question is how can I make the loading of attribute table faster like in joining via properties>join using a custom plugin?

I am using QGIS 2.16


texts: different popup window styles for different tiered locations

texts: New to cartodb, how can I have 1 map with multiple locations, and they can all have different amounts of popup data depending on the plan they purchased to get featured on our map?