texts: ArcGIS javascript api - layer reordering - how to add new layers?

texts:
I am learning ArcGIS javascript API and wanted to make a reordering panel for my webapp. I looked through the ArcGIS sample.

but I'm stuck at addLake() function. If I'm understand the sample correctly, to add the new lake layer

1. get the existing DynamicLayerInfos from the current map
2. create a new DynamicLayerInfo for the new layer
3. push to the DynamicLayerInfos
4. set the new DynamicLayerInfos to the map

The example uses LayerDataSource that reference a table as the new data source for the layer.

And how do you add a new layer from a public layer on a MapServer/FeatureServer. Can someone give me an example?


texts:
I am using SAR image and SNAP tool to extract water bodies. Followed the below link ( http://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-flood-mapping/step-by-step) and could extract the water bodies successfully.

However, the problem what I am facing is along with the water bodies, standing water with crop fields are also extracted.

How can that be avoided?


texts: Finding overlap of three polygons using ArcGIS for Desktop?

texts:
I have three polygon areas showing different land types and I want to create a new layer based on where the three layers overlap each other.

I thought this was a subtraction query but havent been able to achieve this in ArcGIS for Desktop.


texts: TopoToRaster tool execution gives exceptions with ArcObjects?

texts:
I am trying to execute TopoToRaster tool (ArcToolbox>3D Aanalyst Tool>Raster Interpolation>Topo To Raster) in ArcGIS Desktop 10.0 with .NET.

Topo to Raster Input

But get exceptions.  My code is here :

  private ESRI.ArcGIS.Analyst3DTools.TopoToRaster InitTool(out object beginningEnvironmentValue)
    {
        _gp = new Geoprocessor { OverwriteOutput = true, AddOutputsToMap = false };

        beginningEnvironmentValue = _gp.GetEnvironmentValue("workspace");

        _gp.SetEnvironmentValue("workspace", _workspace.PathName);
        _gp.SetEnvironmentValue("scratchWorkspace", _workspace.PathName);

        var topoToRaster = new ESRI.ArcGIS.Analyst3DTools.TopoToRaster
        {
            in_topo_features = _inputFeature,//new InputFeatureTopoToRaster() { FeatureLayer = _inputFeature.ToString(), Field = "HGL", Type = "CONTOUR" },
            out_surface_raster = _outputRaster
        };
        return topoToRaster;
    }

and

private readonly object _inputFeature;
    private readonly string _outputRaster;

    _inputFeature:  is the input feature class path.<<

what is the problem of my input?

I have also tried with class of InputFeature

 public class InputFeatureTopoToRaster
{
    public string FeatureLayer { get; set; }
    public string Field { get; set; }
    public string Type { get; set; }
}


texts:
I'm trying to merge or append layers with identical schemas. When I run the append arcpy function in a stand alone script (only the arcpy tool) it works but when I put it my complete script, the script completes but the output from the tool is empty.

I've tried using shapefiles, converting them to feature layers, and using the append and merge tool (both have the same output).

I'm not sure how to get around this and can't seem to find what could be causing this.


texts: Is there a way to change the display of a raster file into a Robinson projection in Qgis?

texts:
I had asked and answered about the way to make a world map with an angle.

Now I am trying to make the same thing with a raster file.

I had changed the raster CRS, but I did not succeed to make him appear with the shp files

.

When the raster CRS is set on WGS 84 it's display separately from the other layers

.

Actually I am not sure if this is even possible to do.


texts: JavaScript: Adding/removing layers from ArcGIS Server

texts:
I've developed an application which loads the services from the ArcGIS Server via REST API.

Each service has multiple layers, and I can check/uncheck any layer. The problem is the layers are loaded via one URL, and each time I have to remove all layers of the service and load again.

Is there a better way?


texts: Splitting polygon using another polygon in QGIS?

texts:
I have two polygon layers (yellow and orange in the attached image).

I am wondering if there is a way to get a count of the number of isolated areas that the orange polygon has divided the yellow polygon into using QGIS?

I have tried using various vector processing tools in QGIS but have yet to find one that will produce the desired output. I am using QGIS version 3.4.


texts: Change both border and fill colors in one place in QGIS 3.0

texts:
This seems somewhat basic, but is there a way to set the colors for both fill and stroke of a polygon layer at the same time?  I find myself frequently having to change the colors in multiple places, especially when I am using line fill types. Changing the highest level color changes the fills in subordinate layers, but doesn't affect the stroke.

This dialogue allows you to change fill colors in all subordinate styles

but doesn't affect the strokes

Am I missing a setting somewhere?


texts: How to subset point data by outside of polygon data in R

texts:
How can I subset point data from outside of polygon data in R?
Currently, I can subset point data inside polygon data,
now I want to get point data outside polygon data.

library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
nc_point <- st_sample(x = nc, size = 100)
nc1 <- nc[1:50,]

nc_point[nc1]  %>% st_geometry() %>% plot(add=T,pch=16,col="blue")

I want to get the data easily as follows.

nc_point[!nc1]


texts: Pixel size when changing from geographic to projected coordinate system?

texts:
I have a large raster dataset spanning ~5 degrees of latitude N-S. I know the geographic projection has a spatial resolution of 0.00025 degrees, meaning that the dimensions of each pixel in meters  varies depending on latitude.
I have counted the number of pixels in the raster representing different land cover classes. How would I go about converting these pixel counts to area in m2 if the pixel dimensions are variable?
Would simply converting to a projected coordinate system (ie in ArcGIS) standardize the pixel dimensions so that I can multiply pixel counts by a constant area?


texts: Style OpenLayers.Feature.Vector in OpenLayers.Strategy.Cluster()

texts:
I'm using the clustering feature in my project. I collect list of points like so:

cluster_features.push(new OpenLayers.Feature.Vector(
     new OpenLayers.Geometry.Point(lon,lat),
     {heading:deg}
));

Set style:

var clusterstyle = new OpenLayers.Style({
     pointRadius: 12,
     rotation: "${heading}",
     externalGraphic: "ship.png"
});

And set the layer:

clusters = new OpenLayers.Layer.Vector("Clusters", {strategies: [strategy],
     styleMap: new OpenLayers.StyleMap({
          "default": clusterstyle
     })
});

And the problem is rotation in styling. When all the points are unclustered in high zoom, the rotation works just fine, BUT with clustered points it doesn't work at all (the icons/images are all rotation=0 / same direction / no rotation).

P.S. I'm relatively new to OpenLayers.


texts:
I'm using a QGIS project as template for other projects. This template has got around 5000 rows. Because I do several updates to the file, it's getting larger and more complex by the time. I know, that some plugins leave there codes in the file, even if they are not installed to QGIS anymore. I would like to get rid of that residual codes because they slow down the project load time and could cause errors.

Is there a way to automatically clean up the QGIS-project-file? I'm looking for something similar to the "SVG-Cleaner" whit which you can clean up the XML-Code of a SVG.

Until now I was randomly recreating the project by copy and pasting the "good" stiles and codes from one QGIS window with the old project to another window with the new project.


texts: How to Interpolate very large datasets in QGIS?

texts:
I am working with a very large set of point data, (XYAttribute) which I wish to display as a raster.

I've been using Interpolate, this is successful but very slow.

I have other programmes that simply display the data according to a colourbar, is there a faster/less data intensive means of converting my dataset to a raster?


texts: Use Arcpy to Zoom to multiple layers with Def Queries

texts:
I have a project where we create multiple mxd's per geographic area. Each mxd contains a unique definition query for 3 layers. I then select the 3 layers in the Table of Contents view and right-click for "Zoom to Layers". None of the 3 layers are consistently the largest, so I need to approximate this behavior using Arcpy.

Is there a way to do so?

Here are the steps I have taken so far:

Group the layers, then zoom to layer on the group. I have been getting strange behavior from this and the zoom goes out much further than the actual extent when I zoom to the layers individually selected.
Try to use ListLayers to select multiple layers, but I cannot determine if this is possible.
Use the following code to zoom to the extents:

def GoToLayer(LayerNum):
   ...     lyr = arcpy.mapping.ListLayers(mxd, "", df)[LayerNum]
    ...     ext = lyr.getExtent()
    ...     df.extent = ext
     GoToLayer(0)

I get inconsistent zooms from this, where it will zoom out VERY far (58 mil), and then zoom in incrementally (35 mill, 6 mill, etc) as I apply the function again until it finally zooms to the actual layer extent.

My next solution will be to try selecting the individual features of each layer and using "Zoom to selected features".


texts:
I have been digitising a historical map which contains large areas called townlands, filled by smaller areas called fields. I digitised both as two separate vector layers, both polygons.

When digitised the field layers I snapped them to the townland boundaries where the townlands and fields shared a boundary. Unfortunately, I managed to shift the field layer by accident so now the fields and the townlands don’t overlap. Is there a topological function to fix this? To snap one layer of polygons back to another?

I tried v.clean in the GRASS module but this only appears to work for cleaning within one layer.


texts: Force ArcMap to release unused Excel file

texts:
This is not a scripting issue or even a lock file issue.
I have an Excel file with one sheet.
I drag that one sheet into the ToC of ArcMap.
I look at the table.
I remove the sheet from the project.
Now the Excel sheet is locked for editing.
I return to ArcMap and close the project by starting a new one.
I disconnect the folder.
The Excel file is STILL locked.
I have to close ArcMap completely in order to edit the Excel file.
This results in having to go through the full 2 minute start up (yes, 2 minutes on a 2 month old machine).
Is there any other way of forcing ArcMap (10.8) to release its grip on a file it is no longer using?


texts: What is the bounding box of filtered features?

texts:
I'm using geoserver to store a number of gps tracks that are recorded during test drives with a number of vehicles. I have a layer defined that contains the following columns in the database:

date
point <-- postgis geometry
vehicleID
sensorID

That means that all tracks are recorded in the same layer. In geoserver, there is a feature to calculate the bounding box of the entire table. I'm looking for a way to calculate the bounding box for all tracks that fit a number of filters, say vehicleID and sensorID.

I would like to use this information to center the map over the middle of the selected points and to zoom in as close as possible (with all points still being visible).

The filtering is so far done using OpenLayers.Filter on a OpenLayers.Layer.WMS.

I'm looking for any way to achieve this behavior using geoserver.


texts: Possible to resample tifs in numpy array?

texts: I need to resample a large number of tif files. Is it possible to put the list of rasters in a numpy array and resample the array? I'm thinking this would be faster than looping through the list, but since resampling even one raster takes from 5 - 10 seconds normally, maybe it would cause a chokepoint.


texts: Turning on raster painting toolbar of ArcMap?

texts:
My aim is to delete raster cells using the erase tool in the raster painting toolbar. However, the raster painting tool bar won't turn on in ArcMap (it is greyed out), as the screenshot below.

How do I turn the toolbar on and use the erase tool to erase raster cells?


texts: ArcGIS query "dash" character multiple times in a string & replace/remove dash

texts:
Using ArcMap 10.2.2 and have a string field in an attribute table that has values such as:

xxxx-xxxx-x
xxxx-xxxx-x
xxxx-xxxx-x
xxxx-xxxx
xxxx-xxxx
xxxx-xxxx

I want to query for all records that have two dashes in a string (three records in the sample above).
I want to remove the second dash from the string.

I would prefer to do this in two separate steps where I can first query/identify identify the records w/two dashes and then run an update to remove the second dash.


texts:
I'm developing a QGIS Plugin with python.
I need to do a request to a REST server. But when I try to use request module "No module named requests" error displayed.

requests module is installed in my system, but I don't know why QGIS doesn't detect it.

What can I do?

My SO is Windows 7 Ultimate, and QGIS version installed is 2.14.3


texts: 3d arcgis raster .tiff image to CFD software

texts:
I am trying to export 3d raster image prepared in ArcGIS to cfd package such as Fluent. I was wondering if anyone has experience to convert raster image to autocad 3d format, stl, vrml or any other format of its likes which could be imported to cfd package.

I could import format likes IGES, .stl, .DWG, .DXF or even text files.


texts:
I want users to interact with my CartoDB map by interacting with filters. The CartoDB Editor allows to filter data but there doesn't appear to be an option to "expose" the filters to the end-user. Is this something that requires some custom CartoDB.JS scripting or am I missing something ?

The only available user interface options on published maps are toggling through layers and a "search box" in which I am not sure what I am searching.

Lastly, any place where I can find CartoDB driven interfaces i.e. a CartoDB map with custom JS scripting on top of the map content ?


texts: Blank crs attributes in fiona: Can't get proj4 parameters from OpenFileGeodatabase

texts:
I’m having trouble obtaining crs values for geodatabase feature classes in fiona (version 1.7.5). Specifically, the crs attribute returns empty for many of the geodatabases I try to read. This is not an isolated incident -- it occurs frequently with OpenFileGeodatabase types.

For example...

$ fio info --indent 2 /Users/felix/Data/VilasTransportation.gdb
/usr/local/Cellar/python3/3.6.0/Frameworks/Python.framework/Versions/3.6/lib/python3.6/importlib/_bootstrap.py:205: ImportWarning: can't resolve package from __spec__ or __package__, falling back on __name__ and __path__
  return f(*args, **kwds)
{
  "driver": "OpenFileGDB",
  "schema": {
    "properties": {
      "Name": "str:100",
      "NameAlt": "str:100",
      "Type": "str:20",
      "Ownership": "str:30",
      "Municipality": "str:20",
      "Source": "str:30",
      "Comment": "str:250",
      "EditDate": "datetime",
      "Editor": "str:5",
      "SHAPE_Length": "float"
    },
    "geometry": "MultiLineString"
  },
  "crs": "",
  "crs_wkt": "PROJCS[\"NAD_1983_HARN_WISCRS_Vilas_County_Feet\",GEOGCS[\"GCS_North_American_1983_HARN\",DATUM[\"NAD83_High_Accuracy_Reference_Network\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4269\"]],PROJECTION[\"Lambert_Conformal_Conic_1SP\"],PARAMETER[\"False_Easting\",441000.0],PARAMETER[\"False_Northing\",165147.666],PARAMETER[\"Central_Meridian\",-89.48888888888889],PARAMETER[\"Scale_Factor\",1.0000730142],PARAMETER[\"Latitude_Of_Origin\",46.07784409055556],UNIT[\"Foot_US\",0.3048006096012192],AUTHORITY[\"Esri\",\"103463\"]]",
  "bounds": [
    298925.8670191318,
    83800.42979672551,
    583149.6328880489,
    246478.92202439904
  ],
  "name": "TransportationLines",
  "count": 3946
}

As demonstrated below, GDAL has no problem generating a proj4 string for the same dataset.

$ gdalsrsinfo /Users/felix/Data/VilasTransportation.gdb

PROJ.4 : '+proj=lcc +lat_1=46.07784409055556 +lat_0=46.07784409055556 +lon_0=-89.48888888888889 +k_0=1.0000730142 +x_0=134417.0688341377 +y_0=50337.10927101854 +ellps=GRS80 +units=us-ft +no_defs '

OGC WKT :
PROJCS["NAD_1983_HARN_WISCRS_Vilas_County_Feet",
    GEOGCS["GCS_North_American_1983_HARN",
        DATUM["NAD83_High_Accuracy_Reference_Network",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["False_Easting",441000.0],
    PARAMETER["False_Northing",165147.666],
    PARAMETER["Central_Meridian",-89.48888888888889],
    PARAMETER["Scale_Factor",1.0000730142],
    PARAMETER["Latitude_Of_Origin",46.07784409055556],
    UNIT["Foot_US",0.3048006096012192],
    AUTHORITY["ESRI","103463"]]

I need proj4 parameters to use in a pyproj.Proj() instance, so simply grabbing the fiona.crs_wkt value doesn’t solve my issue. Similarly, I want to avoid conversions between proj4 and wkt formats with ogr/osr, or parsing the wkt for parameters.

Why does fiona miss the crs info, but GDAL doesn't? Is it possible to obtain proj4 parameters in fiona with the same consistency as GDAL command line?

The gdb example can be downloaded here (transportation layer) --> http://vcgis.co.vilas.wi.us/vcom/Download.html


texts: How do I count the matching 'many' entries in a one-to-many relationship?

texts:
I have a simple SHP line file - where each line has a unique reference. I have a related text csv file using the same reference, in a one-to-many situation. I need a map with the lines (in the 'one' side of the relationship) colour coded according to the number of matching entries (from the 'many' side of the relationship).

Practically this represents sections of route and repeated dated route inspection reports - and the map should show the total number of inspections for each bit of route. As it happens I also have the same route reports available with the addition of exactly matching geometries in a SHP or TAB file if I need them, but I don't need the geometry information as this is duplicated.

How can this map be produced (given that the format that the data is supplied in is fixed as Mapinfo TAB files which I'd assume to convert to SHP anyway - and that this needs to be an easily repeated process as data is updated).


texts: Using nested for loops to iterate through features in two layers

texts:
I have two layers with 20+ polygon features, layer1 and layer 2. I would like to extract attributes from all of the features in layer 2 that are within 1km of layer 1, for each feature in layer 1. I am doing this by iterating through the features in layer 1, and buffering each feature by 1km. Within this loop, I am iterating through the features in layer 2, and checking whether these features intersect with the 1km buffer:

siteboundary=None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "@@@@@":
        layer 1 = lyr
        break

targetlayer=None
for lyr2 in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr2.name() == "@@@@2":
        targetlayer = lyr2
        break

siteboundaryshapes = siteboundary.getFeatures()
targetlayershapes = targetlayer.getFeatures()

for count,i in enumerate(siteboundaryshapes):
    siteboundarygeom = i.geometry()
    siteboundarybuffer = siteboundarygeom.buffer(1000,25)
    siteboundaryatt = i.attributes()
    print (siteboundaryatt[1])

    for j in targetlayershapes:
        targetlayergeom = j.geometry()
        targetlayeratt = j.attributes()

    if targetlayergeom.intersects(siteboundarybuffer):
        print (targetlayeratt[1])
    else:
        continue

I am trying to return the attributes of the features in layer 2 that intersect with the buffer on layer 1 in the following format:

layer 1 feature 1 attribute

layer 1 feature 2 attribute
layer 2 feature 1 attrbute
layer 2 feature 2 attribute

layer 1 feature 3 attribute
layer 2 feature 1 attribute

layer 1 feature 4 attribute
layer 2 feature 1 attribute

etc etc

layer 1 feature 20 attribute

The above currently does not work; it loops through all of the features in layer 1, but only returns attributes for layer 2 in the first instance that layer 2 intersects with the buffer in layer 1. I have tried moving the indent of the if statement, but this just loops through the features in layer 2 for the first feature in layer 1 and then stops looping through layer 2 thereafter.


texts:
I want to shorten a feature (a polyline representing a road) in QGIS 2.0.

I tried the 'cut features'-tool that look like a pair of scissors. First I chose the right layer, and then clicked 'Toggle'. Then I selected the road and chose the cutting tool.
But, when I selected it, my polyline "disappears"? I never got a chance to start cutting.

I don't remember having any problems with this in QGIS 1.8.
Is there a step that I am missing here?

I really want to get this to work, without having to download QGIS 1.8 again.


texts: Missing roads of San Francisco and extra parts on the map

texts:
I'm using OSM dataset of San Francisco that I have imported using osm2po to PostgreSQL and I think that there are some data missing, because when I use pgr-astar to find the shortest path, there is sometimes no result, so there are some roads between some points that are not available on my dataset.
In other words there are some segments of roads unavailble on my dataset.

Additionally, when I compare the map of San Francisco generated by this dataset with the one generated with the osm file without any convertion or manipulation, it's visible that there are some roads missing.
Does osm2po had eliminate some roads? how to avoid this problem please

The map of San Francisco(osm2po dataset)

The map of San Francisco(osm dataset)


texts: How do you launch QGIS (Mac) without any plugins enabled?

texts: On a Mac, how can you launch Qgis without any plugins enabled?


texts: GeoServer create layer from PostGIS (data type point)

texts:
(Several days ago a user on gis stackexchange suggest that I should use OpenGeo (I was using geojson before). I succeed at what I was doing at that time, but now I have a different problem.)

I was using MySQL, but for easier work with GeoServer I decidec to migrate to PostgreSQL. I spend a few days to migrate a PHP project.

I have a table events with fields: id(integer), name (varchar), ... , and coords(point). I am trying to display some point/icon/whatever on layer, where the event is taking place. I went through different tutorials on OpenGeo, GeoServer, Google ...
For now I manage to create a SQL layer with SQL statement: SELECT id, coords FROM realestate WHERE coords IS NOT NULL.

But when I go to LayerPreview and click OpenLayers->Go, the browser download file (with wms filename) with content

  ?xml version="1.0" encoding="UTF-8" standalone="no"?>!DOCTYPE
  ServiceExceptionReport SYSTEM
  "http://localhost.just_for_validation:8080/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
  ServiceExceptionReport version="1.1.1" >   ServiceException>
  /ServiceException>/ServiceExceptionReport>

(I had to remove all < characters.)

I am not sure if I took the right approach. To sum up: I have PostgreSQL table with several fields, one of them is coords(point). I want to display icon/dot/image/something on map with OpenLayers through GeoServer WMS. I don't know how to create layer which draw dot on lat,lon based on coords field in PostgreSQL.


texts: Mapproxy as transparent Openstreetmap cache

texts:
I followed the guide from https://wiki.openstreetmap.org/wiki/MapProxy to create a locale tile proxy for OpenStreetMap, because I need to access the same tiles over and over again. However, there seem to be some problems with the grid definitions.

As a example if I try to access http://localhost:8080/tms/1.0.0/osm/EPSG900913/6/12/11.png the URL is translated to http://c.tile.openstreetmap.org/7/12/116.png

What has to be changed in the config that mapproxy acts as a purely transparent proxy?


texts: I am working with population data organised in parish units for the whole of Denmark. Many parishes are very small (ie <500 inhabitants compared to a median of around 1500) so I plan to merge these with the next largest adjacent parish. I am working with shapefiles and can manually merge the polygons using QGIS, but with over 1000 small parishes, it would be nice to be able to automate this process! Any suggestions as to how - in QGIS or other software (I also have access to ArcInfo)?


texts: Strip Map Index Feature cuts off map parts?

texts:
I am using ArcMap 10.2 on Win 7 SP1.
When I use the SMIF-Tool for creating a mapbook at a fixed scale it seems to have an issue with overlapping the bounding areas correctly in order to get a proper connection to the following map.
In the example I have a DIN A4-Layout (297 x 210 mm, Portrait) and I have these settings:

Given the track I get these results:

Within the North-South direction the overlapping works quite fine, but when it comes to West-East, connection parts of the following map are missing because ArcMap starts to rotate the data frame within the process of Data Driven Pages. And once you start taking the value of "length along the line" for "length perpendicular to the line" it must get screwy because you still want to populate a regular DIN A4-Layout having 180 mm usable width with now 260 mm "length along the line". So, therefore there has to be a missing part of approx. 4 cm on each side of the A4-page.

My usual workaround at this point is: instead of creating a SMIF use the Grid Index Feature Tool (GIF) based on any given line feature. Within this tool the rectangles have no rotation and they keep in a constant portrait position that is useful for exporting maps within the A4-Layout. The only thing is that I have to bring them in order manually by drag&drop, so I lose quite a bit of automation functionality (sometimes 1000+ rectangles).

No missing parts here.

My question: is there any way either to get that SMIF-Tool working properly or to get the GIF-Tool working without manually draging or anything else?


texts: Georeferencing .hdf file in ArcGIS for Desktop?

texts:
Using a Noggin 250 GPR and Smart Cart. I generated data for a grid. I did not have a gps unit on the smart cart. I do have gps x'y for the corners of each grid. After post processing the data in EKKOmapper. I exported out 3D data. I have two different types of data. One is a .hdf file.

Is it possible to georeference a .hdf file in ArcGIS?


texts: File not found error using Contour in QGIS 3.4

texts:
CSV file load failure in "Contour" procedure, even though the same file will load and display from Layers>Add Layer>.  This has not previously been true.

The problem could be related to the use of forward slash vs. back slash.

The source window for the Layer load uses back slashes and does work.

The source window for the "Contour" load window uses forward slashes and receives the "can't find file" error.  The slashes are not entered by me.  It is entirely a program generated syntax.  Selection for both methods uses the standard windows file selection window.

I previously noticed a similar error when using GRASS through QGIS, but when I used stand alone GRASS, I had no problem.

Is there a setting that can be changed?  Or does it have a different cause?


texts: Export map to PDF using Esri JavaScript API

texts:
I've searched a lot and tried many ways to export a web map to PDF.
I need to export a web map and all active layers.
Can anyone help me?


texts: Couldn't load plugin 'processing' on starting QGIS3 on a mac

texts:
I'm getting an error message when I start QGIS3. I've tried all sorts but don't know what else to do?

 Couldn't load plugin 'processing'

ImportError: cannot import name 'ogr'
Traceback (most recent call last):
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 312, in loadPlugin
    __import__(packageName)
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 674, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/__init__.py", line 28, in
    from processing.tools.dataobjects import *          # NOQA
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 674, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/tools/dataobjects.py", line 51, in
    from processing.algs.gdal.GdalUtils import GdalUtils
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 674, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/gdal/GdalUtils.py", line 38, in
    from osgeo import ogr
ImportError: cannot import name 'ogr'

Python version: 3.6.5 (default, Apr 25 2018, 14:22:56) [GCC 4.2.1 Compatible Apple LLVM 8.0.0 (clang-800.0.42.1)]
QGIS version: 3.2.1-Bonn Bonn, exported

Python Path:
/Applications/QGIS3.app/Contents/MacOS/../Resources/python
/Users/Lauh/Library/Application Support/QGIS/QGIS3/profiles/default/python
/Users/Lauh/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins
/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins
/Library/Frameworks/Python.framework/Versions/3.6/lib/python36.zip
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/lib-dynload
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages
/Users/Lauh/Library/Application Support/QGIS/QGIS3/profiles/default/python


texts:
I am looking for a textbook for ArcGIS Online that is similar to the common textbook, "Getting to know ArcGIS for Desktop".

A text book that simplifies and organizes all the information for Online and provides lessons to practice working with it.


texts:
I've calculated covariance between 7 rasters using R and ArcGIS. It gives me different results. Each times covariance values from R are approximately 1.874 times higher than ArcGIS values.

R - package raster, function layerStats(stack, 'cov', na.rm=TRUE)

ArcGIS - Band Collection Statistics tool (http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000q3000000)

Why covariance values are different in R and ArcGIS?


texts: How do I use the SearchCursor where_clause to get a max value from a column?

texts:
This is my first time asking a question on StackExchange.  I hope I phrase my question and describe my attempts to understand my problem well enough!

While working with SearchCursor, I'm trying to specify the SQL where_clause from a table called "Parcel" within the featureClass, and I want the MAX and MIN value from the column LAND_VALUE.  Below is what I try just for MAX.

import arcpy

feature_class = "C:\foo.ext"
where_clause = "SELECT MAX(LAND_VALUE) AS MaxLandVal from Parcel"

rows = arcpy.SearchCursor(feature_class, where_clause, "", "", "LAND_VALUE D")

I'm told an invalid SQL statement was used.  However, I'm not sure if the SQL statement is invalid, or if perhaps I'm stating the SQL query within SearchCursor incorrectly.  Despite using the reference material for SearchCursor, SQL MAX() and SQL MIN(), and despite referencing a few different examples, I can't get this to work!  Please advise!

Finally, how would I display MaxLandVal if the query were to be correct?  Perhaps something like...?

for blah in rows:
    print blah.MaxLandVal

Thank you for your help!


texts:
I want to do an Ordered Weighted Averaging (OWA) analysis using SAGA in QGIS (not by using pyqgis) as a suitability assessment for recreation/restoration of specific habitats.

There's a number of both positive and negative factors to include in the overlay for which I have .shp polygon files, clipped to my study area. I've converted everything to raster as 10m squares ready for the analysis.

When I tried to run the OWA I got an error saying:

  input layers do not have the same grid extent.

When I checked I saw my rasters didn't line up with the source vectors properly (see below).

The misalignment doesn't seem to be consistent, it's out in different places for each raster/vector pair (despite the original vectors all being aligned) & isn't just a shift in one direction so guessing I've missed a setting somewhere?  I checked the CRS is consistent so it isn't that. I tried using the align rasters tool  (it messes up all the values).

I then wondered if because the vectors cover different extents within the study area (see below) I need a complete grid including areas with no data for each raster for the OWA to work anyway.  So I tried the Grid (interpolation) tool on two vectors using inverse distance to power and making sure the extents were exactly the same to fill all the empty grid squares.  When I run the OWA again just on these two rasters I get a 'list index is out of range' error.


texts: Reducing point numbers in a certain area

texts:
We have a list of poi data in our database(Access mdb db), and we have a web service which serve data by a specified BoundingBox.

However we found that the returned points are too close to another at some area, so I wonder if it is possible to reducing the numbers of the points to make the result more even distributed.

Is this possible?

update:

Data schema:

id name address detail x y type

Where the type means if this point is a hotel or something else. We may take this field as the weight, for example, a point of hotel may have a high priority than that of a untyped point. So when reduce the number of the points, we may take this into account.

For our service, we get data like this:

http://server/ws?minx=10&maxx=20&miny=10&miny=20&size=50 //value of the parameter just for example

And the number of the items to be returned is specified by the paramter size from the client.

Please me know if further information is required.


texts: Where can I find a list of states, districts and tehsil names in India?

texts:
I need only names for making a chained dropdown list for users to select location correctly and submit the data.
Please provide any resources.
I emailed censusindia.gov.in but no reply from them.


texts: Passing multiple arguments to callback function in Leaflet layer method bindTooltip()

texts:
I have a few GeoJSON layers on a Leaflet map that I am using the bindTooltip method on.  Specifically, I am trying to use a function as the argument to bindTooltip so that I can produce some custom html for each of my tooltips.

The documentation on this method states : "Binds a tooltip to the layer with the passed  content and sets up the necessary event listeners. If a Function is passed it will receive the layer as the first argument and should return a String or HTMLElement."

For example:

function simpleTooltip(layer){
    return 'simple tooltip';
}
geojsonLayer.bindTooltip(simpletooltip);

returns 'simple tooltip' as a tooltip for each feature of my GeoJSON layer as expected.

But, if I want to pass another variable to simpleTooltip with the following syntax:

function simpleTooltip(layer, name){
    return 'simple tooltip'+ name;
}
geojsonLayer.bindTooltip(simpletooltip(layer, "steve" );

I get no tooltips and the rest of my <script> doesn't run.

I'm assuming that I'm having these errors because bindTooltip is expecting to receive the name of a function and not the return value of the callback function (in this case simpleTooltip()).

Is it possible to make the callback function be a function of both the layer and an additional variable that helps customize the html(in my case)?

Maybe there is some way to overload or edit the bindTooltip() function to additionally accept a list of arguments to the callback function.

Using Leaflet 1.5.0


texts: Changing color of GeoTIFF using GDAL and Mapbox?

texts:
I have a geotiff that I am importing to a mapping software that I want in the color pink.  It is currently a white image, and I would like a GDAL option to change the color.

I am using Mapbox, which requires an 8-bit .tiff format.  I had a 16-bit tiff that I "rescaled" using gdal_translate -of GTiff -ot Byte -scale 0 65535 0 255 src_dataset dst_datase, which worked.  So now, I have an 8-bit geotiff that I've imported as a layer into Mapbox, but it visualizes as white, and I wish to change to pink.


texts:
I have this projection in PROJ4:

+proj=tmerc +lat_0=41.836944 +lon_0=-87.684722 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

What's the WKT translation?

I saw this site: http://cfconventions.org/wkt-proj-4.html but that doesn't help if you don't know the underlying syntax. For me, it's like being given the definition of 5 Chinese words and being told to construct a sentence out of them.


texts: Clipping topo-data takes too much time and file become too big

texts:
I was following http://www.qgistutorials.com/en/docs/working_with_terrain.html This tutorial.

I added hgt file from SRTM3 as a raster layer, and tried to clip a small area in 1:150000 scale. But this clipping takes eons of time, and the result file seems more than 100gb. I tried the same process with tif file from GMTED data, but the same symptom occurred.

How can I escape from this situation?

Is there anyway to shorten the clipping time, and create result file in small size?

edit 1) I'm not sure how the files looks like, but their size were very small. Looks like they were compressed.

edit 2) Following is the command in command line
gdal_translate -projwin 13984505.7052 4710091.33247 14010189.9721 4674886.87201 -of GTiff -co "COMPRESS=LZW" C:/Users/myname/Documents/QGIS_test/NxxExxx.hgt/NxxExxx.hgt C:/Users/myname/Documents/QGIS_test/NxxExxx.hgt/cliptest.tif

-co "COMPRESS=LZW"  line helped to reduce the size of the file, but still it takes really really long time to clip


texts: GdalInfo and Mapserver : Difference between data

texts:
I'm trying to create a generic XML response from DescribeCoverage response. For my tests, i used Mapserver 7.0.7
With gdalinfo, i got this :

Driver: GTiff/GeoTIFF
Files: geostationary_ir-120_20170619083000_c7142d3e384e192461debbe8c32751c1___b8aedb3f.
tiff
Size is 2288, 2288
Coordinate System is:
PROJCS["Geostationary_Satellite",
    GEOGCS["GCS_unnamed ellipse",
        DATUM["D_unknown",
            SPHEROID["Unknown",6378160,298.2539162964669]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",86.5],
    PARAMETER["satellite_height",35785993.3373],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (-5732339.950191651470959,5732339.950191651470959)
Pixel Size = (5010.786669747946689,-5010.786669747946689)
Metadata:
    AREA_OR_POINT=Area
    TIFFTAG_ARTIST=""
        TIFFTAG_DATETIME=2017:06:19 08:30:00
    TIFFTAG_DOCUMENTNAME="TIFF-MF CMS" 126 155 0
    TIFFTAG_HOSTCOMPUTER=REANIM02
        TIFFTAG_IMAGEDESCRIPTION= 126 155 0
    TIFFTAG_RESOLUTIONUNIT=1 (unitless)
    TIFFTAG_SOFTWARE=Pnm2TiffMF v1.2
    TIFFTAG_XRESOLUTION=1
    TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
    COMPRESSION=JPEG
    INTERLEAVE=BAND
Corner Coordinates:
ERROR 1: tolerance condition error
Upper Left  (-5732339.950, 5732339.950)
ERROR 1: tolerance condition error
Lower Left  (-5732339.950,-5732339.950)
ERROR 1: tolerance condition error
Upper Right ( 5732339.950, 5732339.950)
ERROR 1: tolerance condition error
Lower Right ( 5732339.950,-5732339.950)
Center      (   0.0000000,   0.0000000) ( 86d30' 0.00"E,  0d 0' 0.01"N)
Band 1 Block=2288x1 Type=Byte, ColorInterp=Gray

With Mapserver DescribeCoverage response, i got this :

...
    <gml:origin>
        <gml:Point gml:id="geostationary_ir-120" srsName="http://www.opengis.net/def/crs/EPSG/0/4326">
            <gml:pos>5729834.556857 -5729834.556857</gml:pos>
        </gml:Point>
    </gml:origin>
...

Why there is a difference between these data (Origin Point) ?

Moreover, in XML response, the size of the TIFF is 2287 2287 while in gdalinfo is 2288 2288.


texts: GEE - Applying Savitzky-Golay filter after LS8 NDVI calculation for smoothing curve generating

texts: I'm new in GEE platform, thereby trying to generate a smoothed LS8 NDVI time series by using the Savitzky-Golay filter. I'm following https://groups.google.com/d/msg/google-earth-engine-developers/A7wKQ8WMLIs/kNGt8gwQGQAJ which posted by Guido Lemoine to the GEE Developers Board. I'm almost done the image collection but I am not able to generate the chart. If anyone helps me, for sure I will be glad, let check my code https://code.earthengine.google.com/fc4cfcf05f409e2cc6205ecbf52ec2b3


texts: How do you export a temporal KML from QGIS?

texts:
If I have data with a timestamp attribute, how do I export that from QGIS into a time-aware KML that invokes the GoogleEarth time slider?

I do not see any sort of options in the Save vector layer as dialogue and internet searches have turned up zero.


texts: Given a terrain, how to draw the stream flow path?

texts:
Assuming I have a terrain, as usual the terrain has ridges, creeks and all the characteristics that you can find on a real life map. Water flows from the top of the mountain into lower area, the path that water flows is termed stream flow path.

The terrain is given in terms of triangular irregular network ( TIN), which each point p(x,y) has a z value. How to use this information to construct the stream flow path? What is the physics behind this?

From what I know, steepest descent method can be used to solve this problem. I am thinking about writing my own stream flow algorithm, so I am interested in the theoretical background rather than using the existing tools.


texts: Generate Height map for unity-3D from Lidar point data

texts:
I have Downloaded a data set from udacity open source self driving car program
from here. It is a ROS (Robotic operating system) BAG file more info.The BAG file contains data collected from Velodyne HDL-32E LIDAR and stored as velodyne_packets that I converted to pointcloud Velodyne_pointcloud and stored it in another BAG file.The new bag file I have converted to .csv

But the tutorial I am following to create a terrain in unity from lidar dataset uses .las file in ARCGis . I want to know how to convert the said csv file to .las or is there any other way I could create a heightmap using this lidar data.


texts: What expression needed for universal naming of parts polygon feature class?

texts:
For example, i have n count of polygons and little parts of it.

I want in field [Name] get string: "PNum1(1)".

where "1" - number of main polygon [Num], "(1)" - number of part of main polygon [Num_app].

In field calculator i use this expression: "PNum" & [Num] & "(" & [Num_app] & ")".

On the screenshot is Pnum 13(),Pnum 10()
What i need write in expression if [Nomer_app] is Null and "(", ")" must do not shown (when main polygon don`t have any parts the string looks: Pnum2())?


texts: Collapsing multi-row output of one-to-many spatial join to single-row in ArcPy?

texts:
I'm working with TIGER data, looking at the relationships among incorporated places and counties. Not all places are completely contained by counties; some do span multiple counties.

In the final output I'm preparing, I need to have a "location reference memo" for each place; for a city like Cleveland that is completely contained, this might say something like

"Incorporated Place in Cuyahoga County"

for a city like Birmingham, AL, this would instead need to say something like

"Incorporated Place in Jefferson County, Shelby County"

or

"Incorporated Place in Jefferson, Shelby Counties"

-- or something like that.

I know that I can do a spatial join to get a multirow output like

Birmingham | Jefferson
Birmingham | Shelby
[city in AL] | [County in AL]

Is there any facility, perhaps after the fact, for generating a single row equivalent of this? E.g.,

Birmingham | Jefferson; Shelby


texts: How to convert kmz raster with rotation to GeoTiff?

texts:
I have a map that I have georeferenced in Google Earth and I'd like to easily use it in QGIS. The problem is that this map does not have north at the top, so the georeferencing involves rotating the image in Google Earth.
There are multiple answers here, elsewhere, and web sites for converting a kmz (saved from Google Earth) into a GeoTiff, but I haven't found one that works if the raster file has been rotated in Google Earth. gdal_translate seems to ignore the rotation parameter that Google Earth saves in the kmz or kml file. Also, if you just import the kmz as a layer into QGIS it also comes in unrotated.
Clearly a solution is to do the georeferencing in QGIS, but I was wondering if anyone knows a faster solution.


texts:
I would like to use special pie chart for my map, where the angle is not important, but area or radius of the arc carries information. I found, that this type of pie chart is calling Coxcomb/polar/rose (in Czech it is known as comparative) chart (see http://blogs.sas.com/content/jmp/files/2012/03/spending1.png or http://visual.ly/learn/coxcomb-chart), basicaly it is bar chart masked in pie :). Do you have any suggestions how to create this chart automaticaly or semiautomaticaly in gis (I prefer QGIS)?

Thanks for help,
Daniel


texts: How to use Matplotlib in Python 2.6 installed by default ArcGIS 10.0 installation?

texts:
I am wondering if anyone has successfully run Matplotlib with the default installation of Python 2.6 that comes Arcgis 10.0. It does not seem to work with the python version installed by Arcgis 10.0. I was able to download and successfully run other packages such as Scipy. I have looked at many blogs where they talk about using other installations of python such as Python(x,y) etc. I prefer to run it without changing the version of python. Please post the details of the steps if anyone was successful in running Matplotlib. I have installed Matplotlib version 1.1.0. For example, the following command to check version works.

      import matplotlib

      print matplotlib.version shows
      1.1.0

I think the problem is with other things that comes with Matplotlib such as pylab, pyplot. Every time I try to import pylab or pyplot, the python IDLE environment is RESTARTED. Please give me the details of steps (from download, installation to running it) if anyone of your are successful in using with the default python installation of ArcGIS 10.0.


texts: How do I use 'model builder' in QGIS 2.2?

texts:
I am trying to learn how to use the 'model builder' in QGIS 2.2. I have used model builder in ArcGIS, but I don't understand how to do it in QGIS 2.2?

What I want to do is to make buffers, clip vector layers and, if it's possible, add columns and data to attribute tables?

I want to do the same process for quite a few different layers, so it would be nice to get help from a model.

Is there any instructive tutorials to be found online that can show me how to make a model?
I tried to google it, but didn't find what I needed.


texts: Recover lost polygons from shapefile in QGIS

texts:
I've committed a newbie error it seems in GQIS with a polygon shapefile.

I wanted to remove about 15 of the 210 polygons, and put those 15 in another file.

I made a duplicate of the polygon shapefile.
In the duplicate file I removed everything from the attributes table, started cutting and pasting lines from the original to duplicate file attribute tables.
I saved the edits in duplicate with the 15 polygons, gave it a new name.
I tried to save the edits in the original now with 195 polygons, wouldn't let me.
I made another copy of that original, but it only had the 15 removed polygons.
Started closing error messages, thinking it would allow me to save.
I tried closing edits of the original, I selected discard changes. All remaining 195 polygons dissapeared, all 15 polygons removed are now present.
Closed project without saving changes, but alas the shapefile had been saved. The cpg, prj and qpg files seem unchanged with the original date, so I copied those for safekeeping. Loading the project again it comes back with the 195 polygons missing (about 8 hours work).

Is this a lost cause or is there a way to recover them?


texts:
Is there any document for using NVS vector stream tool? I found difficulties when using it for my data (shape file). It says  "Streams Layer must contain a TNODE_ field denoting the arc's (stream's) destination node". The question is how to calculate TNODE?

Regards,

Anas


texts: Why do features in WGS84 disappear from the (ArcMap) display when zooming in really close?

texts:
I routinely import centerlines and other data from various counties that are in different projections, so I've begun to convert everthing to WGS84 so all is in the same projection.  I want to match up topology where centerlines from different counties meet, so in order to do that, I have to zoom in (really) close to see if centerline endpoints match up.  However, since I've converted stuff to WGS84, when I zoom in close to features, they disappear from the display, and I can't zoom in close enough for what I need.  Anyone else experience this, and is there a fix?  Using ArcGIS 10.1 Advanced...

EDIT

By "really close" I mean sub-centimeter at times, maybe even sub-millimeter... essentially I want to do is be able to create a network dataset so I need the lines to connect, and I need a way to identify places where they don't connect, even from small distances...

If WGS84 is not high enough precision, what would you recommend for a "universal" projection (mainly for United States) that is "high precision"?  Projections have always been kind of confusing for me (as well as others I'm sure...)  Thanks for your response/comments...


texts: Writing a Batch File (.bat) Script to run a QGIS Python Script without opening QGIS

texts:
I have saved a python script (named as 'Task3.py'). Currently, I will open QGIS and then run this Task3.py . This python script is saved in C:\Users\Sreeraj.qgis2\processing\scripts\Task3.py .

But, I want to write a Windows Batch File (Task3.bat), such that when I run this Task3.bat file, the Task3.py will get automatically executed (even without opening QGIS). This batch file script must be saved in E:\Sreeraj\Task Global Data CONVERSION\Split\Task3.bat .

How do I write Task3.bat (so that I can run Task3.py without even opening QGIS) ?


texts:
I can't stop getting this error message when I try to run the following code in Qgis 2.14.10 :
from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtCore import QSettings
from PyQt4.QtCore import QSettings
from qgis.core import QgsVectorLayer, QgsDataSourceURI

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "name of my database", "my user name", "my password")
uri.setDataSource("database schema", "table name", "geom")
couche = "PM"
layer = QgsVectorLayer(uri.uri(), couche, "postgres")
if not layer.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayers(layer)

TypeError: QgsMapLayerRegistry.addMapLayers(list-of-QgsMapLayer, bool addToLegend=True): argument 1 has unexpected type 'QgsVectorLayer'

Could someone tell me what I'm doing wrong please ? :)


texts: How to only display active layers in legend?

texts:
In QGIS, in the Print Composer's Legend function, how do I make it so that the Legend automatically displays only the active layers (i.e. the layers that are checked and therefore displayed on the map)?

If I check the "Auto-Update" box, it just generates a list of ALL the layers (including the inactive ones), and I have to manually go back and remove the inactive layers from the Legend which gets to be quite a pain after a while.


texts: How to highlight polygons within another polygon?

texts:
I want to symbolize public buildings in OSM data.

For Example: There is a polygon and in column "amenity" it has the value "school". But i don't want to highlight this polygon, but the buildings inside this polygon.

How can I do this?

Thanks!


texts: How to import a GeoJSON file from GeoServer using OpenLayers

texts:
I have been trying to import the layer that I created in GeoServer. I need to use OpenLayers. This is what I have so far.

<!-- language: js -->

<script type="text/javascript">
var serviceUrl = 'https://gis.arlingtonva.us:6443/arcgis/rest/services/public/Bike_Routes/MapServer/';
var layer = '0';
var otherURL = 'Specialty/ESRI_StateCityHighway_USA/MapServer';

var raster = new ol.layer.Tile({
      source: new ol.source.OSM()
    });

var map = new ol.Map({
  target: 'map',
  layers: [raster],
  view: new ol.View({
    center: ol.proj.fromLonLat([-77.0910, 38.8816]),
    zoom: 13.58
  })
});

var bikes = new OpenLayers.Layer.WMS(
 "Bike Routes of Arlington", "http://localhost:8080/geoserver/ArlingtonBikePath",
 {
         srs: 'EPSG:2283',
         layers: 'Bike_Route',
         styles: '',
         format:'image/png'
 }
);

map.addControl(new OpenLayers.Control.LayerSwitcher());
var myswitch = new OpenLayers.Control.LayerSwitcher();
map.addControl(myswitch);

</script>


texts:
Does anyone know what the current state of the SQL Server extension for GeoServer is?

When it came out in 2008, there were a lot of complaints of poor performance, as the comments in this blog show: http://blog.geoserver.org/2008/11/10/146/.

Has any further work been done on it since then?


texts: Can I access Worldclim 2.1 CMIP6 data in 30-second resolution?

texts:
I am trying to access the Worldclim 2.1 CMIP6 future climate data at 30-second resolution. The Worldclim website says 30-second data will be available by then end of March 2020. This is long overdue and I am not sure if maybe there is another way of accessing it or similar downscaled CMIP6 data?
I have tried emailed their support email a few times over the last several months with no response.
If not possible, I'm looking for similar downscaled data for different CMIP6 SSPs (for the Arctic in particular) if anyone has directions.


texts: Iterate Feature Selection Error -Export grouped features

texts:
I am trying to group values inside my shapefile using Iterate Feature Selection and then Copy selected feature to new shapefile.

I don't know why but my copy features don't allow me to using %Value% expression. What am i doing wrong?

I am using Arcgis Pro Basic 1.3.


texts: Custom QGIS 2.6 map tool - mouse click and release issue

texts:
I have create a custom map tool for QGIS 2.6, implementing a canvasReleaseEvent and canvasPressEvent. I can't though get the signals to work on one single mouseclick (I have to click twice).

My map tool was implemented like this:

class IdentifyGeometryTool(QgsMapToolIdentify, QgsMapTool):
    def __init__(self, canvas):
        self.canvas = canvas
        QgsMapToolIdentify.__init__(self, canvas)

    geomIdentified = pyqtSignal(['QgsVectorLayer', 'QgsFeature'])
    def canvasReleaseEvent(self, mouseEvent):
        results = self.identify(mouseEvent.x(), mouseEvent.y(), self.ActiveLayer, self.VectorLayer)
        if len(results) > 0:
            self.geomIdentified.emit(results[0].mLayer, results[0].mFeature)

    canvasClicked = pyqtSignal(['QgsPoint', 'Qt::MouseButton'])
    def canvasPressEvent(self, mouseEvent):
        pnt = self.toMapCoordinates(mouseEvent.pos())
        self.canvasClicked.emit(pnt, mouseEvent.button())

It was then created and connected to a couple of functions in a plugin:

class myPluginClass:

    def __init__(self, iface):
        self.clickTool = IdentifyGeometryTool(self.canvas)

    def initGui(self):
        self.clickTool.canvasClicked.connect(self.firstFunction)
        self.clickTool.geomIdentified.connect(self.secondFunction)

    def firstFunction(self, point):
        do stuff

    def secondFunction(self, layer, feature):
        do more stuff

Clicking the canvas (on a feature on the active layer) the first time connects to firstFunction and clicking one more time connects to a secondFunction. What can I do to call each function in one mouse click and release?


texts:
In QGIS (2.10.1 on a Mac), I am trying to make some multiline labels on a map I have with a bunch of trees in them. In the attribute table I have columns with the Latin scientific name, the English name (when one exists), the year the tree was planted etc.

I am labelling the trees in my map with an expression like this:

genus_species + '\n' +
case
    when  "english_species" is not null then  "english_species" + '\n'
end
+ to_string( "planted_year")

This works well for the rows where english_species is not null, however whenever I have a row with no English name, the entire expression just seems to break and nothing is printed, not even the scientific name or the year (which always exists).

I have tried adding a

when  "english_species" is null then ''

But that doesn't work either. The expression just breaks for rows with no value in the field english_species.


texts: Searching for detailed shapefiles in Vietnam: Tam Ky and Quy Nhon cities

texts:
I'm currently interning with UN-Habitat Vietnam in Hanoi, and am helping them produce some basic maps for community development strategies in the coastal cities of Tam Ky and Quy Nhon. This is my first time working in GIS internationally, and I'm realizing quite quickly the differences between here and North America in terms of reliable data acquisition. I'm having a really hard time tracking down shapefiles for the study areas, and am hoping someone with experience doing GIS work in Vietnam could help out with some pointers.

One thing I should mention, is that I have been provided fairly detailed DWG files for both cities. I personally don't have any experience converting CAD into GIS, but I have begun to poke around in AutoCAD to see what I can come up with. However, the file layers are not organized in any way I can comprehend, and many of the polylines are disconnected, and therefore don't translate well to GIS.

Does anyone know where to look for shapefiles in Vietnam?


texts: Calculating distance of sites to riverways using ArcGIS Desktop?

texts:
How do you calculate the distance of about 10 sites scattered across a region to bodies of water (rivers and lakes) using ArcGIS 10.1?

I have the lat/lon coordinates of the 10 sites as a point feature class and I have downloaded the Surface Water Area (polygon shapefile) for the study  area from the UK Ordnance Survey webpage.

I want to determine the easiest route-ways that are available in the study area by taking into account slope and proximity to rivers. I understand that I can achieve this by running a cost path analysis if I have a dem raster, slope raster and distance rasters to rivers and valley floors. But I am not sure how to create distance rasters to rivers and valley floors using Spatial analyst tools.


texts:
I am trying to set the progress bar. I have  a list contain 1000 numbers but progressbar showing up to 100.

a=range(1000)
for i in a:
    arcpy.SetProgressor("step", "Copying shapefiles to geodatabase...",0, len(a), 1)
    arcpy.SetProgressorPosition()


texts: Problem creating a raster mosaic to run a viewshed analysis

texts:
I am using the freely available dtm data from the OS UK website (Land-form panorama and its in the ASCII format). when you download the data from the website, the dtm is found in different subfolders (labelled according to the national grid squares, e.g. NN, NO, HY etc.). within each folder there are number of ASC files (which looks like tiles when you view on ArcGIS - I am not sure what the correct technical term is). My sites are spread out across an area of about 23 national grid squares - so there is a large number of ASC files that needs to be added. When I try to open the dtm data on ArcGIS 10.1 I encountered a problem, I can only open the ASC files from about 4 national grid squares before the program freezes and shuts down. The primary focus of my study is to conduct viewshed analysis for my sites. when I add the various ASC files and try to conduct a viewshed analysis it looks like this:

The problem is my sites are in clusters and some sites are found spread out into the neighbouring tiles. it looks like when I try to do the viewshed analysis it doesn't take into account the sites that are found in the nearby tiles - which is a huge problem for me.

I have been told to create a raster mosaic but when I tried doing that it still comes up looking like separate tiles:

Am I doing something wrong? I am also having trouble performing the viewshed analysis here (although I can perform a hillshade), I was able to do it before I created the mosaic dataset. The other problem is I still have to do the hillshade and viewshed analysis one asc file at a time. is there a way around this?

one other thing - when I added the raster data to create a mosaic dataset - it comes up as empty green squares - the terrain only becomes visible if I calculate the statistics for each individual asc file - it remains blank when I try to do a batch statistics. I am not sure what I am doing wrong, if you could help me it will be greatly appreciated!


texts: Is there a default setting to save ArcGIS 10.3.1 .mxd to 10.1

texts: Does anyone know if there is a setting in ArcGIS Desktop 10.3.1 where we can permanently set the save function to version 10.1? I know that there is Save a Copy function, but I don’t want to create duplicate mxds for different versions. I have read through this site but it seems to just talk about the Save a Copy function rather than having a permanent change to the save function.


texts: Finding which points from dataset within specified boundary using CartoDB?

texts: I'm a non-coder using Cartodb.  I've uploaded a dataset that has longitudinal and latitudinal coordinates, and overlayed that with a map with neighborhood boundaries.  Is there a way to easily compare (or sort) these two layers?  I want to be able to quickly find out which locations reside inside a certain area, without having to look at the visualization and figure it out manually.  Thoughts?


texts: How to update the Z value using an UpdateCursor?

texts:
I have created the following script to mimic the ArcGIS 3D Analyst tool, FeatureTo3DByAttribute, which creates a 3D copy of a 2D feature class and sets the Z to a specified attribute.

import arcpy, os
from arcpy import env

# Parameters
in_features = arcpy.GetParameterAsText(0)   #Feature Layer
out_features = arcpy.GetParameterAsText(1)  #Feature Class
height_field = arcpy.GetParameterAsText(2)  #Field

# disable overwrite protection
arcpy.env.overwriteOutput = True

# enable environment outputZFlag
arcpy.env.outputZFlag = "Enabled"

# copy input to output
arcpy.CopyFeatures_management(in_features, out_features)

# Use cursor to update geometry with z value based on height_field
with arcpy.da.UpdateCursor(out_features, ("SHAPE@Z", height_field)) as cursor:
    for row in cursor:
        row[0] = row[1]
        cursor.updateRow(row)

Unfortunately, the z value is 0 for all of the output features.  I can't find any documentation as to why this is.  An alternative would be read every point and change the z value.  Does anyone know why it is failing to update the z?

EDIT:
Additional discussion

https://geonet.esri.com/thread/93322

https://geonet.esri.com/thread/82123


texts: Coverting TIFF file by Raster Calculator in QGIS

texts:
My ultimate goal is same as here
Here is my steps:

The image is created using Python and saved as PNG format:

Import the image into QGIS using Georeferencing tool and saved as TIFF:

Reclassified the image value using Raster -> Raster Calculator in order to use processes such as r.thin and r.to.vect later according to here:

But the result is weird:

Not sure what's wrong here.
Or are there any other solutions?


texts: Changing extent of a raster based on a different raster

texts:
I two rasters with differen't dimensions, but they overlap.  I want to change the extent of the larger raster so that it is 'clipped' to the smaller raster so I can stack them.

I found out how to determine where the intersection from this post, How to add different sized rasters in GDAL so the result is only in the intersected region

but I think I am getting something wrong in the final step or two.

So far I have this:

#read in the files
larger_raster = r'C:\large.tif'
larger_raster = gdal.Open(larger_raster)

smaller_raster = r'C:\small.tif'
smaller_raster = gdal.Open(smaller_raster)

#get the geotransformation information
gt = smaller_raster.GetGeoTransform()
larger_raster = lc.GetGeoTransform()

#get the xmin, xmax, ymin, ymax of each
r1 = [gt[0], gt[3], gt[0] + (gt[1] * smaller_raster.RasterXSize), gt[3] + (gt[5] * smaller_raster.RasterYSize)]

r2 = [gt2[0], gt2[3], gt2[0] + (gt2[1] * larger_raster.RasterXSize), gt2[3] + (gt2[5] * larger_raster.RasterYSize)]

#find where they intersect
intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

#only select the locations where they intersect from the larger raster
new_extent = larger_raster.ReadAsArray(intersection)

this last line is where I go wrong though and this error is returned:

TypeError: in method 'BandRasterIONumPy', argument 3 of type 'double'

I have also tried this (picking up after finding the intersection):

xoffset = int(abs(intersection[0] - gt2[0]) / gt2[1])
yoffset = int(abs(intersection[1] - gt2[3]) / gt2[5])

 new_extent = lc.ReadAsArray(xoffset, yoffset, int(gt2[1]), int(gt2[5]))

but this returns:

ValueError: negative dimensions are not allowed

Also,

gt2 returns:

(-11119505.203637687,
 463.3127136230469,
 0.0,
 8895604.157333,
 0.0,
 -463.3127136230469)

intersection returns:

[-8895604.157333, 7783653.637667, -7783653.637667, 6671703.118]

xoffset returns:

4800

and yoffset returns:

2400

The code runs when I change the new extent calculation to this:

new_extent = large_raster.ReadAsArray(xoffset, yoffset, 1, 1)

but the output shape is (1L, 1L) whereas the shape of smaller_raster is (2400L, 2400L).


texts: Unable to clip by mask a raster in QGIS 3.2

texts:
OS: Windows 8 with QGIS 3.2

when attempting the following clip-by-mask exercise, the result works on 2.18 but not 3.2.

The data I am using is the QGIS training data with the Alaska Shapefile and the SR_50M_alskaa_nad.tif raster.

The Alaska shapefile needs cleanup before using as a mask and here are the steps followed:

Load the Alaska shapefile
Enable Check Geometry and Topology plugins
Run Vector/Check Geometries
Run Check Topology with no overlaps, no dups, and no gaps
Run Grass/Vector/v.clean setting sliver to 10,000
Load SR_50M_alaska_nad.tif
Processing/Toolbox/GDAL/Raster Extraction/clip raster by a mask

The error message is
GDAL command output:

Warning 1: Ring Self-intersection at or near point -2632521.5885887709 5137728.2501277877

ERROR 1: Cutline polygon is invalid.

Is there another cleanup tool I can run on this file in QGIS 3.2 to correct this issue?


texts: Is it possible to configure the popup window to display only visible elements

texts:
I currently have a map that contains a MapImageLayerwith several sublayers, which I can toggle visibiliy with a LayerList. The problem is that, when clicking on the map, the popup opens even for invisible elements.

Is there a way to tell the popup engine to show results for visible elements only?


texts:
Is there any built-in function or quick hack to efficiently return the closest element (index) to a point for a given geopandas GeoSeries?

Something like

>> geopandas.GeoSeries(
>>     {
>>         0:Point(1,1),
>>         1:Point(3,2),
>>         2:Point(2,3),
>>     }
>> ).closest_element(Point(0.99,0.99))
>  0

I could calculate the distance to every single element and return the sorted list index's first element, but this solution is really slow and I feel like there should be some better way to do it.


texts:
Does anyone know of a software or service that can define whether a location is inside or outside of city and county jurisdictions?

Our company is getting hit with huge tax audits because we did not collect sales tax in certain cities that impose a 1-2% tax rate as well as counties. I'd like to be able to imput an address, have it converted to GPS coordinates for accurate location and provide location information on whether a location is inside or outside of city and county jurisdictions so the correct rate can be calculated.


texts: Generating grid that matches Google Earth Engine polygons

texts:
I'm trying to create a grid in R over the whole continent of Africa (I received some help in this GIS StackExchange post), but I'm having trouble generating a grid that matches the CRS used in the Google Earth Engine, which I'm going to use later to download satellite imagery.
The discrepancies are higher when I depart from the Equator line. For instance, here is some code to generate a grid of 10 km by 10 km patches over South Africa, using the stars library:
library(sf)
library(stars)
library(giscoR)
library(tidyverse)

# country shapefile
south_africa <- gisco_get_countries(country="ZAF", resolution = 10)  %>% st_transform(3857)

# method 1: stars
grid = st_as_stars(st_bbox(south_africa), dx = 10000, dy = 10000)
grid = st_as_sf(grid)
grid = grid[south_africa, ]

Here I used CRS 3857, so I have my units in meters to define the grid. I can compute the area of the pixels, and get the geometry of a pixel so we can see that the area is correct (1, because I normalized it) under this coordinate reference system:
grid$area = st_area(grid$geometry)/(10000)^2
print(paste('area:',grid[[3]][[10000]]))
grid[[2]][[10000]]

[1] "area: 1"
POLYGON ((2331832 -3517624, 2341832 -3517624, 2341832 -3527624, 2331832 -3527624, 2331832 -3517624))

I also checked the grid using st_make_grid, and I obtain exactly the same results when checking the same patch in the grid:
grid_2=st_make_grid(south_africa, cellsize = 10000)
grid_2 = st_as_sf(grid_2)
grid_2 = grid[south_africa, ]
grid_2$area = st_area(grid_2$geometry)/(10000)^2
print(paste('area:',grid_2[[3]][[10000]]))
grid_2[[2]][[10000]]

[1] "area: 1"
POLYGON ((2331832 -3517624, 2341832 -3517624, 2341832 -3527624, 2331832 -3527624, 2331832 -3517624))

The problem is that if I take this polygon and import it to Google Earth Engine, I see that the dimensions of the polygon are not right. In this case, I get that the square root of the area is 8,637 meters, and the average of the length of the sides of the polygon is 8,632 meters, but both measures should be 10,000 meters if the geometries where properly defined.
I also take a Landsat-8 Satellite image around the centroid of the polygon, covering 10km by 10km, and I obtain an image that is bigger than my polygon. Here you can see the satellite image, and the polygon in black:

Here is the code to visualize this in the code editor of Google Earth Engine.
Is there any way in which I could fix this? I probably need to define a CRS that matches the one in Google Earth Engine. I tried working with CRS  4326, but I couldn't define a grid in meters with that CRS.


texts:
I download this data from the Argo website.
When I tried to open it in ArcGIS 10, it produced only one set of LAT/LONG but I know there should be more than that.

I tried using R.

 library (ncdf4)
 nc_open("E:\\Python\\20150123_prof.nc")

So I opened it. But, I do not know how to proceed because I do not know how to plot this in R. I want to get the sea surface temperature of all points during the date I downloaded and interpolate them. Has anyone get an experience on this in R?

EDIT:

 This was the result of the code.


texts: Editing with auxiliary features (lines, points etc.) in ArcGIS Desktop?

texts:
In order to digitize exact features in ArcGIS Desktop 9.3.1, our draftsman would like to use auxiliary features (lines and points mostly) that don't have to be saved later on. I told them the easiest way was to digitize the auxiliary features into a seperate shapefile - is there a simpler or better way?

(in case I didn't explain what I want very good I'll give an example. Suppose someone has to digitize a lighting post (point feature) 4 meters from a curb. In Autocad our draftsman would have digitized a line four meters long, using snapping he'd put a lighting post feature in the end of the line (in a different layer), and in the end he'd erase the line.

Is there any way to mimck the process without actually digitizing into a shp/gdb?)


texts: Using point geometry object in Spatial Join with ArcPy?

texts:
I'm trying to do something that is conceptually very straightforward.  I'm pulling a center point from a raster (geotiff, etc) and then doing a CONTAINS operation with that point against a UTM zone grid.  The intent is to obtain some field values for the polygon feature the raster center point falls within.  Then, I can use those values to initiate a Project Raster method using the correct UTM zone.  I have already done this with GDAL/QGIS but need to try and make it work in the Esri world.

Here is my code:

import arcpy

def getRasterCenter(raster):
    xmin = arcpy.GetRasterProperties_management(raster, 'LEFT')
    xmin = arcpy.GetRasterProperties_management(raster, 'RIGHT')
    xmin = arcpy.GetRasterProperties_management(raster, 'BOTTOM')
    xmin = arcpy.GetRasterProperties_management(raster, 'TOP')

    center_x = (float(xmax.getOutput(0)) + float(xmin.getOutput(0))) / 2
    center_y = (float(ymax.getOutput(0)) + float(ymin.getOutput(0))) / 2
    center_coords = center_x, center_y

    return center_coords

#stash intermediate data in a scratch db
tempData = arcpy.CreateScratchName(workspace=arcpy.env.scratchGDB)

#write json to file since arcpy doesn't read json string
utm_grid_geojson = '{ validated Esri json content here }'
file_json = r'C:\Temp\utm_zones.json'
with open(file_json, 'w') as outfile:
    outfile.write(utm_grid_geojson)

# set up the raster
terrain_file = r'C:\Temp\raster.tif'
terrain_file_center = getRasterCenter(terrain_file)
x = float(terrain_file_center[0])
y = float(terrain_file_center[1])
point = arcpy.Point(x, y)
terrain_center = arcpy.PointGeometry(point)

# sanity check the point
print(terrain_center.WKT)

# now get the utm grid again and do the join
utm_grid = arcpy.JSONToFeatures_conversion(file_json, tempData)

zone = arcpy.SpatialJoin_analysis(utm_grid, terrain_center, tempData, "#", "KEEP_COMMON", "#", "CONTAINS")

# check if valid result
result_desc = arcpy.Describe(zone)
print('result type: {}'.format(result_desc.featureType))

This code results in the following error:

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid. ERROR 000670 Output Feature Class is same as input Target Features.

I am trying to avoid creating intermediate data sets if possible as this script is only one small part of a much larger model.  Looking for a bit of hand-holding to understand how to do spatial joins via arcpy, I guess.


texts:
How do I manage floating point calculations in QGIS?

I have carried out 3 consecutive calculations on landsat bands. I end up with 32 bit floating point rasters on which I can't carry out a 4th calculation. The calculation result gives a substantial file but the values of the pixels are shown as "Inf".

I am not able to find a syntax guide for the float() function although I use it successfully in 1 of the preceeding calculations. And I do not have an in-depth understanding of how to manage integer/floating point or how to convert between them in QGIS.

I found these but no big clues

How to reclassify a float raster?
Floating point: understanding their inaccuracy


texts: Is possible to merge landsat 8 with High Resolution Satellite Imagery?

texts:
i have landsat8 and High Resolution Imagery (1m) in the same study area, is possible to merge them to get one imagery so that i can take advantage of the accuracy satellite imagery.

some_text http://rghost.net/8WGsT7Kgy/image.png

some_text http://rghost.net/6N9clcVYQ/image.png


texts: Write a layer in a specific table of a spatialite file

texts:
I'm manipulating several layers using python and qgis and I would like to save all those layers as separate tables in the same spatialite file.

At the moment I'm using :

QgsVectorFileWriter.writeAsVectorFormat(layer, output_file_path_str, 'utf-8',
  crs, 'SQLite', False, None, ['SPATIALITE=YES', ])

However, this only enables me to store one layer into one file. I can't figure out how to store several layers into the same file.

In a qgis unit test file i found this :

options = QgsVectorFileWriter.SaveVectorOptions()
options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer

But it seems to be outdated, i get errors such as AttributeError: type object 'QgsVectorFileWriter' has no attribute 'SaveVectorOptions'.


texts: Remote sensing Classification and texture analysis compared with QGIS zonal statistics

texts:
I am Using Geomatica to analyse suburban spectral signature in Las Vegas. Ran a texture analysis on all band images. Combined texture bands with original image bands to run an unsupervised ISODATA classification. I aggregated the classification to 3 classes (urban, non-urban, and suburban). I now want to compare the results of both quantitatively thus run zonal statistics. I presently at school only have access to QGIS. I assume that both my PIX files (texture analysis combined image bands, as well as classification image bands) are raster based.

Presently converting my PIX files to TIFF to make it usable in QGIS. Seeing as QGIS is my only present option to use, is it possible for me to compare two rasters here?

I am a NOOB with GIS and have no scripting knowledge. Opening the zonal statistics plugin I see only a vector vs raster option.

Does this mean I would have to convert one of my files?

I'm also not sure what would be correct in this case, do I treat my image file that contains the classification aggregation as the vector, and the texture image as the raster, or the other way around?

Lastly, both images contain heaps of data, Texture analysis is combined with 6 EMR bands, while the classification image has 14 bands (original 6, texture bands, plus 2 additional channels for the classification output and class aggregation results).

Prior to attempting zonal statistics, do I extract the classification data itself and compare solely that against the texture analysis bands?


texts:
I'd like to know if there are any GPS emulators available?
What I'd like to do is the following: Given a series of straight lines and a vehicle moving at a certain velocity on the lines, I'd like the GPS emulator to write data about the vehicle's position to a computer port or somewhere in /dev/.
I'd like to  query this port using Python and read information. I googled but I'm not sure what to make of the information available. GPSFake is one of the results. Can it do what I'm asking?


texts:
I work with QGIS 3.6.0 Noosa and closed the WMS-window in the log messages panel by accident (see image).

How can I get that window back?


texts: How to get Borough names using ArcGIS API for Python through geocoding?

texts:
I'm new to ArcGIS API for Python. Is it possible to do do geocoding of an area, specifically New York City and get all the boroughs' names?

i.e.: Geocode request for all boroughs in New York city.

in[x]: geocoding.geocode('New York, NY')
out[x]:{'address': 'New York',
'location': {'x': -74.00713999999994, 'y': 40.71455000000003},
 'score': 100,
 'attributes': {'Loc_name': 'World',
  'Status': 'T',
  'Score': 100,
  'Match_addr': 'New York',
  'LongLabel': 'New York, NY, USA',
  'ShortLabel': 'New York',
  'Addr_type': 'Locality',
  'Type': 'City',
  'PlaceName': 'New York',
  'Place_addr': 'New York',
  'Rank': 2.5,
  'AddBldg': '',
  'Nbrhd': '',
  'District': '',
  'City': 'New York',
  'MetroArea': 'New York City Metro',
  'Subregion': 'New York County',
  'Region': 'New York',
  'RegionAbbr': 'NY',
  'Territory': '',
  'Country': 'USA',
  'LangCode': 'ENG',
  'Distance': 0,'},
 'extent': {'xmin': -74.25713999999994,
  'ymin': 40.46455000000003,
  'xmax': -73.75713999999994,
  'ymax': 40.96455000000003}}

and get  in the 'District' element: 'Bronx, Brooklyn, Manhattan, Queens, Staten Island'


texts: Compute Moran's I on more than one attribute

texts:
I am trying to compute global Moran's I to check for spatial autocorrelation among areas of a city. I am using python, in particular geopandas and pysal.

My geodataframe looks like:

District_name|District_code|Attribute_1|Attribute_2|...|Attribute_n|geometry

my areas are polygons, and each attribute is a mean value related to the inhabitants of the different areas (e.g. mean age, average years of schooling, average income, etc.).

I have no issue in computing the Moran's I for a single attribute. I first compute a spatial weights matrix for a rook neighborhood, then I calculate the index:

import geopandas as gps
import pysal as ps
urban_data=gpd.read_file(gdf_file_path)
w=ps.weights.Rook.from_dataframe(urban_data)
mi=ps.Moran(urban_data.Attribute_1, w)
print(mi.I)
print(mi.EI)
print(mi.p_sim)

and I get

0.1225081380916777
-0.002932551319648094
0.001

which looks pretty good to me. But, beside the interpretation of this result, I really do not understand how I can check for spatial correlation among ALL the different attributes. I'd like to discover if there issome kind of pattern considering all the different features, not a single one each time.
Is that possible? Or am I completely losing the point of spatial autocorrelation?


texts: In QGIS what files do I need to recreate this map?

texts:
http://gaez.fao.org/Main.html#

I am trying to use the maps on this website in my QGIS project. When I download the zip it comes with a few different types of files but I do not know which ones to use to make it look like the preview.

Here is a image of the zip:


texts: Converting XY coordinate points to Lat Lon

texts:
I have a set of points representing the nodes making up a polygon from a city-border. These points are projected in RD New EPSG:28992 (The Netherlands). I want to convert these coordinates to Lat Lon.
I have tried to save the layer as a new layer with the WGS:84 EPSG:4326 projection and opened this new layer in a new project with the same CRS.
However, this either shows a completely warped view of the points in a wrong place or it simply shows nothing.

Original shape in RD New

Warped points in Lat Lon/WGS:84
Am I missing something obvious here? I am not new to QGIS and have done this trick with other CRS's lots of other times.


texts: Datatype and filter for a SDE connection parameter in a Python toolbox?

texts:
I have a tool in a Python toolbox. As a parameter, I want to let the user select an existing (or possibly create a new) SDE connection. What value for datatype should I use for this?

I have tried DERemoteDatabaseFolder, definied by the documentation as:

  The database connection folder in ArcCatalog.

However, it only lets me select files of mysterious "Basic Types" and does not show me any SDE files. I also tried DEFile, setting param.filter.list = ['sde']. But that options required the user to find the SDE files on her own, which is not so easy if they are located somewhere in the AppData\Roaming folder.

Is there any better alternative?


texts: ArcCatalog: Can you add attributes to a relationship class at a later date?

texts:
I have a relationship class that could use a new attribute (or two or three).  I have looked for a way to do this in ArcCatalog and have not discovered it yet.  I have:

Opened the Relationship Class Properties by right clicking on the Relationship class.
Previewed the Relationship Class and right-clicked in the Column Headings.
Looked in ArcToolbox>Data Management Tools>Relationship Classes.

None of these attempts have resulted in a new attribute.  If anyone knows how to do this, I would appreciate an answer!


texts: Extract pixel values at given point dataset and if masked find the nearest pixel value

texts:
I need to extract the values at each pixel foe a feature collection from an image and if they are masked the value at nearest pixel,
my feature collection is named table and the image collection I am using is Landsat 7 Collection 1 Tier 1 Annual NDVI Composite

var mapfunc = function(feat) {
 var geom = feat.geometry()
 var addProp = function(img, f) {
   var newf = ee.Feature(f)
   var date = img.date().format()
   var value = img.reduceRegion(ee.Reducer.first(), geom, 30).get("NDVI")
   return ee.Feature(ee.Algorithms.If(value,
                                  newf.set(date, ee.String(value)),
                                  newf.set(date, ee.String('No data'))))
  }
  var newfeat = ee.Feature(ndvi.iterate(addProp, feat))
  return newfeat
 };

 var newftNDVI = table.map(mapfunc);

This code extracts values at said pixels but does not extract the nearest pixel values for the masked ones.


texts: Is my WKT projection definition wrong or is there a bug in PROJ?

texts:
I'm using PROJ with a WKT projection string to try and match a coordinate system being employed in a 3rd party solution.
I seem to be able to match the 3rd party coordinates when projecting EPSG:3857 coordinates to my WKT definition, but when I try to convert back again the Y value is correct but the X values is completely wrong.
I'm pretty new trying to use WKT projection strings so I'm wondering if there's something wrong with my definition?
This is my custom WKT definition:
PROJCS["PacificIslands",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4167"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",-21.2],
    PARAMETER["central_meridian",-159.75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","100001"]]

If I take the following EPSG:3857 coordinates:
X: 22295668.360
Y: -2414760.379

Projecting them using the WKT given above:
X: 3674.005
Y: 909.963

Then if I project them back to EPSG:3857 again I get:
X: -17779348.325
Y: -2414760.379

Any idea what's going on?
I've tried this in both the C++ PROJ 8.0.0 library and using the C# ProjNet library and get the same result.
Here's my C++ PROJ example code in case it's of any use:
#include <stdio.h>
#include <proj.h>

int main(void) {
    PJ_CONTEXT* C;
    PJ* P;

    /* or you may set C=PJ_DEFAULT_CTX if you are sure you will     */
    /* use PJ objects from only one thread                          */
    C = proj_context_create();

    P = proj_create_crs_to_crs(C,
        "EPSG:3857",
        "PROJCS[\"PacificIslands\",GEOGCS[\"NZGD2000\",DATUM[\"New_Zealand_Geodetic_Datum_2000\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6167\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4167\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",-21.2],PARAMETER[\"central_meridian\",-159.75],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"100001\"]]",
        NULL);

    PJ_COORD coordFrom, coordTo, coordReversed;
    coordFrom = proj_coord(22295668.3602, -2414760.3791, 0, 0);

    coordTo = proj_trans(P, PJ_FWD, coordFrom);
    printf("ORIGINAL x: %.3f, y: %.3f\n", coordFrom.enu.e, coordFrom.enu.n);
    printf("CONVERTED x: %.3f, y: %.3f\n", coordTo.enu.e, coordTo.enu.n);

    // Reverse again
    coordReversed = proj_trans(P, PJ_INV, coordTo);
    printf("REVERSED BACK x: %.3f, y: %.3f\n", coordReversed.enu.e, coordReversed.enu.n);

    /* Clean up */
    proj_destroy(P);
    proj_context_destroy(C); /* may be omitted in the single threaded case */
    return 0;
}


texts: Aligning data in GCS_North_American_1983 and GCS_North_American_1927 using ArcGIS for Desktop?

texts: I have tried everything that I have found on this site so far and nothing has worked!  I have used ArcGIS for 3 years but had never had to mess with coordinate systems as the 1st company I worked for had it down. However the new company does not and I am getting data from the county GIS which is correct, and trying to match the data form the 2nd company with the counties data. The data from the county is GCS_North_American_1983 the company data is GCS_North_American_1927. Not sure of what information is needed but I can get the projection, limits, what ever is needed. Please help as the is really confusing!


texts: Slow rendering of 25 MB raster in ArcScene

texts:
I am loading one raster layer that is 25 MB in size. It is a 1926x2042 1m DEM. I need to render it at a cell size of one for optimal display of the layer. When I open a saved .mxd with the only the raster, it takes ~ 10 minutes. Every time I go into the layer properties and hit OK, it takes about 5 minutes to load, even if I don't change any settings. If I change the color ramp it takes even longer.

I'm running Arc from a 2 TB hard drive on Parallels. The virtual machine has 2 GB of RAM designated, and 2GB of virtual memory. I allocated 4 GB of virtual memory according to ESRI's help page (http://desktop.arcgis.com/en/arcmap/latest/extensions/3d-analyst/optimizing-arcscene.htm), but still have the same problem. I've come across people with issues with 100-200MB files, but not with a 25 MB file.

When the layer does load, I can move around fine and this is great. However I have more similar sized rasters to load into ArcScene, and am in a crunch for time.

If this is a normal amount of time for a file of this size, does anyone know if it is possible to speed up rendering times?


texts: WMS layer details not clearly visible in Leaflet

texts:
I cannot figure out where the problm is. I use GeoServer and Leaflet. My problem is that when I visualise some layers (in particular that with the names of the clients) I cannot see them all at each zoom level as I in the layer preview in Geoserver. Why??
The SLD must be correct otherwise I wouldn't have a correct preview right? Then when I load the site in Leaflet and I zoom in and out I cannot see all the names of the clients. I also have a similar problem with pie-chart: in layer preview they are perfect, while in leaflet are cut or not visible. Where is the problem? Moreover, is there a way to solve it??


texts:
I would like to put some images from Earth Engine (EE) in my report. So far I screenshot it and paste PNG on MS Word.

I see EE can export images to Google Drive. But its geotiff, I can't just paste it on Word, even to see them is tricky.

Curious how you guys do it or what's the recommended way to go? And interested the most in a solution that works on macOS as well.

Perhaps you are doing the following:
GeoTiff convert to PNG with transparency?


texts: Determine mudflats area from remote sensing images

texts:
I need to determine the mudflat area for a group of small islands in the Red Sea. I have read some publications and it seems to be possible using remote sensing, for example Landsat. But I am not exactly sure how to go about it, I know I should try to get images at low tide of my area of interest but I am not sure how to go about selecting the mudflats areas.

The most intuitive way seems to be to select pixels with a specific RGB combination, since the mudflats seem to have a very specific color. I am just interested in mapping the extent of the mudflats.

Someone has suggested to me "multispectral classification using a standard method like the maximum likelihood classifier" and "ratio analysis or decision-tree classifications" but I am not a remote sensing person and these seem to be way too technical.

edit: I have available argic 10.2 full license and might be able to get Imagine 2014


texts:
I was hoping someone could help me get Spyder working with QGIS. I am finding the python editor limiting within QGIS as I would like to debug my code, use the variable explorer, and be able to run highlighted code. I have tried two routes:

Route 1: Install Spyder within the OSGeo4W QGIS python instance

I installed OSGeo4W Network Installer (64 bit) v3.2 and opened OSGeo4W Shell and typed:

py3-env
python -m pip install spyder
spyder3

But I am getting errors in my IPython console:

SyntaxError: invalid syntax
File "C:\OSGEO4~1\apps\Python36\lib\site.py", line 177
file=sys.stderr)

Also, when I put in

from qgis.core import *
import qgis.utils

into my python script, it complains with unable to detect undefined names

I noticed all the QGIS files are in C:\OSGeo4W64\apps\qgis\

Route 2: Create new conda environment and install QGIS via conda

I used the install instructions at this url: https://anaconda.org/conda-forge/qgis

I create a new conda environment and then typed:

conda install -c conda-forge qgis
conda install spyder
start spyder

My IPython console is fine (i.e. no errors as route 1), but I get errors when trying to import qgis.core etc into my script as in route 1.

I noticed that conda installed all the QGIS files to C:\Users\<user_name>\AppData\Local\Continuum\miniconda3\envs\<my_env>\Library\

Would anyone have any idea how to get this working?


texts: Will .sbn and .sbx files (spatial index) speed up operations in GDAL/OGR?

texts:
I'm trying to intersect large datasets in OGR and python using SQL queries like ST_Intersects and ST_Contains. These operations run very fast in ArcMap but I'm designing a system to automate them, hence using OGR and python. My code is taking forever to run (days). I've read recently that spatial indexing will help speed up these operations. I can generate these files by indexing the shapefiles in ArcMap, which generates .sbx and .sbn files. If those files accompany my .shp .dbf .shx and .prj files in the corresponding directories, will the speed of the operations benefit?

Alternatively, will the operation speed increase if the shapefiles are loaded in to a postGIS database instead of read in via the file system?


texts:
With PostgreSQL I can query an entire row for a string:

SELECT * FROM table_name WHERE table_name::text LIKE '%string%';

Is there an equivalent (CQL?) for a Geoserver WFS serving PostGIS data?

I've had a look here and here but found nothing.


texts:
I would like to use custom styled tiles for quite a small area. I know there are general solutions like mapnik or mapbox but since I am really only interested in showing tiles for one forrest, is there a way to pre-render those and but the resulting png's on a static assets-server?

How would you go about this? Setup a local OSM-Server, crawl mapnik for that area and dump the data yourself? Or is there a tool/cli-option for this particular use-case?


texts: I am looking for the best way in either Drone2map or ArcGIS orthoimagery tool to mask parts of individual photos from being projected onto the final mosaic ortho. Like moving traffic so we can use projection from image that doesn’t have the vehicle. Like how pix4d allows drawing mask and changing which image is used.


texts: Iterate through each raster for executing an equation to calculate index using Python

texts:
I am working with multiple rasters, and in need to execute an equation on each raster, and save them on a date-wise for each raster in an out specified folder. Input Rasters are as follows

"S1A_IW_GRDH_1SDV_20150226T003942_20150226T004007_004787_005F2A_8D84.VV" # out name for this raster 20150226

"S1A_IW_GRDH_1SDV_20150230T003942_20150226T004007_004787_005F2A_8D84.VV" # out name for this raster 20150230.

Below code, I am trying to execute, but its giving error
dates   = []                        # empty lists to start
Rasters = []

for thisFile in os.listdir(r'D:\C1'):
    fN, fE = os.path.splitext(thisFile) # break into name and extension
    if fE.upper() == '.TIF':
        Rasters.append(thisFile)

for ThisRaster in Rasters:
    RasterNameSplit = ThisRaster.split('_')
    if len(RasterNameSplit) > 5: # omit non-conforming names
        ThisDateString  = '_'.join(RasterNameSplit[:9]) # Will take only 9th from raster e.g 20150226
        if not ThisDateString in dates:
            dates.append(ThisDateString)

# validate the dates
ValidDates=[]
for ThisDateString in dates:
    L8_red = 'S1A_IW_GRDH_1SDV_{}_20150226T004007_004787_005F2A_8D84.VV.tif'.format(ThisDateString)
    if os.path.exists(os.path.join(r'D:\C1',L8_red)):
        ValidDates.append(ThisDateString) # both red and nir bands exist

for ThisDateString in ValidDates:
    L8_red = ''S1A_IW_GRDH_1SDV_{}_20150226T004007_004787_005F2A_8D84.VV.tif''.format(ThisDateString)
    outFile = "out_{}.tif".format(ThisDateString) # make the name unique

# open
redds = gdal.Open(os.path.join(r'D:\C1',L8_red))

redband = redds.GetRasterBand(1)

red = redband.ReadAsArray()

red = red.astype(np.float64)

#The actual calculation
def SR_pqa(red):
    SR = Abs((red - ( -22.97)) / (- 0.01 - (-22.97)))
    return SR

SR = SR_pqa(red)

#Write the out file
driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create("out.tif", redds.RasterXSize, redds.RasterYSize, 1, gdal.GDT_Float32)
CopyDatasetInfo(redds,dsOut)
bandOut=dsOut.GetRasterBand(1)
BandWriteArray(bandOut, SR)

Error message:
AttributeError                            Traceback (most recent call last)
<ipython-input-4-4787ad7a0460> in <module>
     30     # rest of your script
     31
---> 32 redband = redds.GetRasterBand(1)
     33
     34 red = redband.ReadAsArray()

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'


texts: Does setting a precise position of the base matter when using RTK and you only care about relative position?

texts: I'm using RTKLIB with my own rover and base station. Most of the time I have a "float" solution. Would setting a more accurate position of the base station in the RTKLIB settings help to have more a "fix" solution? I don't care about an absolute solution, only relative.