texts: What CRS is the USGS Global Land Cover Characteristics Data Base?

texts:
I have trouble bringing land cover data from the USGS into a "decent" coordinate system.

I obtained an .img file of the "Global Land Cover Characteristics Data Base Version 2.0" from http://edc2.usgs.gov/glcc/tabgeo_globe.php. I renamed the file to *.bsq and created an .hdr file according to the following description on the USGS website:

  9.2 Geographic Projection Parameters

  The data dimensions of the Geographic projection for the global land
  cover characteristics data set are 21,600 lines (rows) and 43,200
  samples (columns) resulting in a data set size of approximately 933
  megabytes for 8-bit (byte) images. The following is a summary of the
  map projection parameters used for the Geographic projection:

  Projection Type: Geographic

Units of measure: arc seconds

Pixel Size: 30 arc seconds

Radius of sphere: 6370997 m.

XY corner coordinates (center of pixel) in projection units (arc seconds):

    Lower left: (-647985, -323985)

    Upper left: (-647985, 323985)

    Upper right: (647985, 323985)

    Lower right: (647985, -323985)

So I have created the following .hdr file:

nrows   21600
ncols   43200
nbands  1
nbits   8
layout  bsq
skipbytes   0
ulxmap  -647985
ulymap  323985
xdim    30
ydim    30

When I now load the file into QGIS, it looks fine. But to convert the file from the original CRS to one of my choosing, I would have to know the original CRS in the first place. And this is where I am stuck: What is the CRS of the "Global Land Cover Characteristics Data Base Version 2.0" database? Can someone tell me from the aforementioned projection infos? I tried EPSG:4326 and EPSG:3857, but they do not seem to be correct.


texts: ArcGIS tiles on OpenLayers are very blurry

texts:
I'm trying to get this map on my openlayers map

When I just add the layer in the code below, I do get the tiles, correctly referenced, but very blurry.

var url = "http://services.geodataonline.no/arcgis/rest/services/Geocache_UTM33_EUREF89/GeocacheBasis/MapServer";
layer = new ol.layer.Tile({
  source: new ol.source.TileArcGISRest({
    url: url
  })
});

The ArcGISRest seems to request tiles from the wrong zoom level. Is there a way to get tiles of a higher zoom level?

I've tried to switch to ol.source.XYZ too (see below for the source definition), then it does show good quality tiles, but I can't get it to use the right projection (changing the projection option actually doesn't change the map).

source: new ol.source.XYZ({
  attributions: ["Kartverket, Geovekst og Kommuner - Geodata AS"],
  url: url + "/tile/{z}/{y}/{x}",
  projection: ol.proj.get("EPSG:25833"),
})

Here's the layers meta-information: http://services.geodataonline.no/arcgis/rest/services/Geocache_UTM33_EUREF89/GeocacheBasis/MapServer

Could someone help me with either the projection issue in the source.XYZ layer, or the quality issue in the TileArcGISRest layer?


texts: Updating large Shapefile with new polygons using QGIS?

texts:
In QGIS I have a multipolygon layer that has about 100 polygons in it (property boundaries) and I will need to continually add polygons to it as they become available.

What is the most efficient way to merge the smaller polygons into the larger file without creating a new shapefile each time?

Would it be easiest to continue using the merge tool and overwrite the original shapefile?

Each polygon was taken from a large dataset that contains all private properties for the state and each has only one ring. There is very little attribute data with each polygon and the process runs very quickly in only a few seconds.


texts: Multibeam bathymetry xyz data 53000000 features ArcGIS 9.3

texts:
I'm working with an xyz file that has 53000000 points, I have only been able to see this in Microsoft Access, because in text it says that the amount of data is too large. So what I've tried to do:

Change .xyz to .txt ---> 3D Analyst tool, conversion, from file, ASCII 3D to Feature Class: So when I try to convert to natural neighbor I can't see the z data so I cannot do the interpolation. I also try using the Tin tool, but it just say that an error happened when it try to draw it.
I've opened .txt with Microsoft Access (it's when I saw the amount of data that I have), and I saved this table in .mdb, but I have no idea on how to use this to make an interpolation.

Please I need some help to work with this data!


texts: Convert coordinates from EPSG 4326 to EPSG 28992 in PostGIS

texts: I imported a CSV file with coordinates into a PostGIS table. This table contains a column "latitude" and a column "longitude". The coordinates in these columns are in EPSG:4326. What is the best way to convert these coordinates to EPSG:28992? Is it possible to do that directly, or should I convert these coordinates to geometry first?


texts:
I would like to edit the metadata attributes of an tiff file. In particular the TIFFTAG attributes as described here.

In Python it's a piece of cake:

from osgeo import gdal
ras = gdal.Open('test.tif', gdal.GA_Update)
ras.SetMetadataItem('TIFFTAG_DATETIME', '2016:08:12 13:00:41')
ras = None

Or you just use gdal_edit.py:

gdal_edit.py -mo TIFFTAG_DATETIME='2016:08:12 13:00:41' test.tif

Up until Version 0.5.27 rgdal seemed to have functions available for this (see here).

Does anyone know whether this can still be done?


texts: Display and label a stack of Sentinel 1 images with the date

texts:
I have a stack of Sentinel 1 images and I want to display all of them one by one and want the label to be the date of corresponding image. I am trying a for loop for that but I am unable to visualize the results.
Link :https://code.earthengine.google.com/e58e28ed2fd14cc2b8cbaf4ef666b5bb


texts: Can not get the connecting line to align with label in desired place in QGIS 3

texts:
I have labelled a points layer by defining a line between points and labels, so as I could remove the label to outside the map area (Figure 1). Using additional columns in the points layer attributes table ("x" and "y"), I can move the labels manually.
To make a line from the point to the label, I have added a 'Geometry Generator' in the Style section of the layer's Properties dialog. To define the geometry generator, I have used the following expression where $x and $y are the features' coordinates and "x" and "y" are the labels' coordinates:
make_line(make_point($x,$y), make_point("x","y"))

My issue is that labels do not align with the line connecting them to features. The labels are basically where I want them, but the line heads elsewhere (see Figure 1). How do I get these lines to basically be positioned from the first letter of the town name to the point?

Image shows map with points being towns and the labels their town name. The lines are supposed to connect the town name label (e.g. Warrnambool) to the point (the town location), but the lines are fixed in odd directions


texts: Not Able To See Buffer on Click Event On Graphics Layer

texts:
Using ArcGIS JavaScript API 3.21 and This Demo Fiddle I am trying to create a bufferd area on click event but apparently this is not showing any graphic on bufferLayer graphics Layer

require([
  "esri/map",
  "esri/SpatialReference",
  "esri/graphic",
  "esri/layers/GraphicsLayer",
  "esri/Color",
  "esri/geometry/Point",
  "esri/geometry/geometryEngine",
  "esri/symbols/SimpleMarkerSymbol",
  "esri/symbols/SimpleLineSymbol",
  "esri/symbols/SimpleFillSymbol",
  "esri/renderers/SimpleRenderer",
  "dojo/domReady!"
], function(
  Map,
  SpatialReference,
  Graphic,
  GraphicsLayer,
  Color,
  Point,
  geometryEngine,
  SimpleMarkerSymbol,
  SimpleLineSymbol,
  SimpleFillSymbol,
  SimpleRenderer
) {
  var map = new Map("map", {
    basemap: "streets",
    center: [-122.4, 37.785],
    zoom: 14,
  });

  var buffSymb = new SimpleFillSymbol(
    SimpleFillSymbol.STYLE_NULL,
    new SimpleLineSymbol(
      SimpleLineSymbol.STYLE_SHORTDASHDOTDOT,
      new Color([105, 105, 105]),
      2
    ), new Color([255, 255, 0, 0.25])
  );

  var graphicsLayer = new GraphicsLayer();
  var bufferLayer = new GraphicsLayer();
  map.addLayer(graphicsLayer);

  map.on("click", function(evt){
    var buffer = geometryEngine.geodesicBuffer(
      new Point(evt.mapPoint.x, evt.mapPoint.y),
      1000,
      "meters"
    );
    bufferLayer.add(new Graphic(buffer, buffSymb));
    map.addLayer(bufferLayer);
  });
});

I am  not getting any error message but also no result back! can you please take a look at this and let me know what I am doing wrong?


texts: Connectivity Issue in Road Network - Error Graphic Pick 1

texts:
I have read around 6 various forums posts and still could not figure out where I was making the mistake so please help me.

I have downloaded the OSM shapefile and and have clipped a pretty large area which I need for some analysis I need to make.

I have attached the shapefile I was working on. Exact Coordinates I am having problem are: 55.733952 24.180805 Decimal Degrees

Referring to Attached Picture Route is unable to resolve for Point 4 to 5, while up to 4 it works fine.

The shapefile is here.


texts: Carto - Search for point by name and zoom to

texts: I have a map I've put together with Carto. The main layer is a brewery point layer. I want to set up a search box that allows the map reader to search for a brewery by name. I know I can set up a filter that does this. However, the crucial piece is that the map then zooms to the brewery when one is selected. Is this possible using Carto.js?  I'm pretty sure I could do it in Leaflet but wouldn't be able to access the data which is in Carto.


texts: PDAL Reader.las and Writer.las warning: 'found invalid value of '0' for point's return number'

texts:
I am trying to read and write LAS (*.las) file using PDAL but getting the following warning messages by PDAL's LAS reader and writer for ASPRS LAS files.

(pdal pipeline readers.las Warning) C:/denoise\lake.las:
      Found invalid value of '0' for point's return number.

(pdal pipeline readers.las Warning) C:/denoise\lake.las:
      Found invalid value of '0' for point's number of returns.

I'm trying to remove noise from las file using the code below:

 { "pipeline":[ "C:/denoise/lake1.las", { "type": "filters.outlier", "method": "statistical", "multiplier": 3, "mean_k": 8 }, { "type": "filters.range", "limits": "Classification![7:7],Z[-100:3000]" }, { "type": "writers.las", "compression": "true", "minor_version": "2", "dataformat_id": "0", "filename":"C:/denoise/clean.las" } ] }


texts: Compare workspace schemas in Modelbuilder

texts:
My model has two model parameters:

A file gdb workspace or sde workspace containing the data
An empty "template" file gdb

As the very first process in my model, I want to check if the schema of the data workspace matches the schema of the template gdb. If they don't match, the model should stop and preferably give some meaningful error message (I suppose by using some Python in a Calculate Value tool). I know this can be done simply by using the Export XML Workspace Document tool now in 10.1, but I'm still on 10 for now, so I need my models to work properly on 10.


texts: Creating local web GIS server using Node.js

texts:
This is a question related to OpenLayers/jQuery without CDN?. I was looking to copy a minimal set of files on js/css files to replace online CDNs, so that I can run a test server on localhost without internet access.

I heard from the comments that it is possible to use node.js to run a test server entirely on the localhost. This is a interesting approach. I was wondering what are the general steps, e.g. to run such a webgis server with (or that matches) typical options like OpenLayers and GeoServer?


texts: Copy Polygon from one feature class to another, but auto create relationship records

texts:
FYI, I have a very limited GIS background, but I can author code.

I have a Stands (polygon) layer and I want a user to be able to copy stand(s) from this layer and paste it to another polygon layer (say a HarvestActivity).

When the paste occurs I want to record the following in a related table of the HarvestActivity feature class that the polys were pasted into.

StandId - Id of stand from source "stand" layer

AreaOfStand - this is the area of the polygon of the HarvestActivity layer that covers a stand from the stands layer.

AreaPercent - the percentage ratio to the original stand

Another important note is that if someone edits the polygons in the activity layer, i.e. cuts or extend it. then the related table needs to be updated (changing the ratios, adding records, deleting records) automatically.

E.g. If I extend a polygon in the HarvestActivity Layer and that polygon overlaps with 3 stands in the stands layer then we would need to have 3 records created.

What is the best way to handle this in the GIS world? We will be using Portal and ArcGis Online, thus support of offline sync needs to be also considered.


texts: Listing all overlapping polygons using PostGIS?

texts:
Assume a table of overlapping polygons where the number of polygons overlapping a single area varies from 1 to n. I am trying to find a PostGIS/SQL query that will return a list of all of the polygons that contribute to a particular overlap.

I am currently using a two pass approach that is 1/2 in SQL and 1/2 in Python. I believe that this approach could be significantly improved, but am unsure exactly how.

SQL: This query returns a table of polygons that represent all of the overlaps between all of the input polygons.

SELECT ST_AsEWKB(geom) AS geom FROM ST_Dump((
    SELECT ST_Polygonize(the_geom) AS the_geom FROM (
        SELECT ST_Union(the_geom) AS the_geom FROM (
            SELECT ST_ExteriorRing(footprint).geom)  AS the_geom
            FROM study) AS lines
    ) AS noded_lines
  )
)

This picture might help to describe what is occurring above. The polygons 'contributing' to the intersection are noted in red. The above query computes the polygons, but does not result in an array of 'contributing members'.

Python: The second pass is to take the original input layer and intersect it with the results of the above SQL query.

Question: Does a one pass solution in pure SQL exist for this spatial problem? The posted SQL gets the polygon boundaries, but I am unsure how to identify which input polygons created that boundary in a way that does not presume the number (or max number) of contributing polygons.

This question is similar to: Find and list all polygons that overlap with another polygon except answers along the lines of:

SELECT *
FROM table as a, table as b
WHERE ST_Intersects(a.Geometry, b.Geometry) AND a.id > b.id

result in a selection of pairwise overlaps.


texts:
In my flask application I have a folium map displayed inside an iframe. I defined a folium html popup with a form.  How can I route the form?
Here is my popup html code:
"""<!DOCTYPE html>

                {}	""".format(row['properties.NTAName']) + """

                <form id="restForm" method='POST' action='/exp_neigh' style="border: 1px solid #333333;padding:20px; margin-bottom:20px;">
                    <input type="hidden" id="lat_boc" name="latitude" value="{}" """ .format(lat) + """>
                    <input type="hidden" id="lon_box" name="longitude" value="{}" """ .format(lon) + """>
                    <button type="submit" class="button">Select this neighbourhood</button>
                </form>

            """

Here is my Flask route:
@app.route('/exp_neigh', methods=['GET', 'POST'])
def exp_neigh():

    return render_template('neighbour.html')


texts: How can I list all of the sites on my CARTO map on the side of the map for easy access?

texts: I have created a map in CARTO that I would like to have all of the different sites in a list on the side of the map (similar to what the current widget does, but with more than 5 sites, and ideally zooming to each when clicked). Is there an easy way to do this?


texts: How can I plot WorldPopProject data (GeoTIFF) in R with an appropriate scale?

texts:
I am a complete beginner to using any GIS tools and struggling through some very basic concepts.
I am trying to load data from the WorldPopProject (https://www.worldpop.org/geodata/summary?id=6365) and though I'm successful in converting the file into a raster that I can plot, the scale is from 0 to 1. The site says "The units are number of people per pixel" and I have to imagine more than 1 person can live in a pixel when the resolution is 3 arc/100m at equator.
library(raster) #Required for raster files. Also loads 'sp' package
library(rgdal) #Required to manipulate GEOTIFF and advanced tools
library(tiff) #Required for reading in TIFF
library(ggplot2) #Plotting
#Seemingly all the GIS functions were built under 3.6.2 (?)

ghanaPopLocation = "filepath_inserted_here"

# Loading and conversion
ghanaPop = readTIFF(ghanaPopLocation) #Get initial GEOTIFF loaded
ghanaPopRaster = raster(ghanaPop) #Convert GEOTIFF to raster

# Plotting
plot(ghanaPopRaster)

Am I misinterpreting this plot, the data itself, or have I done something wrong?


texts: Arcscene export issue: only upper layer showing up

texts:
I have created a 3d visualization with ArcScene 10.1.

I have 1 hillshade of a portion of an area, which take elevation ifno from DEM of the same area. Also, the DEM is coloured and overlayed onto the hillshade. Plus, vectors data are added (i.e., rivers and springline).

Now, I have tried to export the scene as 2D; I have used various format like PDF, TIFF and other formats. The issue is that the exported scene is just the upper layer, i.e., the DEM. The Hillshade as well as the vertoc layers do not show up.

Any insight in to issue?


texts:
I have a raster layer of some roads and I want to interpolate it such that on the new layer, each cell has a value corresponding to how far it is from a road cell. The end result would look like a set of ridges (high values) for the roads, sloping away to low values as you move away from the roads.

What would be a suitable command for this processing? I believe 'interpolating' is the correct term but I am not certain as I am new to GIS.


texts:
I have been trying to mosaic two Worldview images together using the QGIS merge tool. However, the result produces a very noticeable offset.

Having visually checked the original images, they do have a small offset from each other.

Have people encountered this in the past. If so, how did you solve it?


texts: Using Make Grids and Graticules Layer in ArcMap?

texts:
I'm trying to use the Make Grids and Graticules Layer tool in ArcMap  (have a ArcInfo license).   I'm looking for a basic example like in North America.

eg.

Start with empty map
Set data frame to  NAD83 UTM Zone 1N
Create File geodatabase
Create GridExample dataset in FGDB in NAD83 UTM Zone 1N
Open; Make Grids and Graticules Layer tool

Select grid template C:\Program Files (x86)\ArcGIS\Desktop10.0\GridTemplates\Quad_24K_NAD83.xml
Select Extent, and Same as Display from drop down list
this sets Left: -5620263.767724, Top: -5620263.767724, Right: -5620263.767724, Bottom: -5620263.767724
Select the GridExample dataset in FGDB for Input Feature Dataset
Output layer is set automatically to Quad_24K_NAD83

Run the tool

Output is as follows, but fails to create the grid.  Creates the feature classes ok.   Can anyone hint on improving this example to create the grid?

  Executing: MakeGridsAndGraticulesLayer "C:\Program Files
  (x86)\ArcGIS\Desktop10.0\GridTemplates\Quad_24K_NAD83.xml"
  "-5620263.767724 -5275836.261749 5621263.767724 5276836.261749"
  C:\data\FGDB\Grids.gdb\Grid1 Quad_24K_NAD83 Quad_24K_NAD83 24000
  -3.18645767674729E-04 "3 Centimeters" "2 Meters" PROJCS['GCS NAD 1983 CORS96 UTM Zone 60N
  (Calculated)',GEOGCS['GCS_NAD_1983_CORS96',DATUM['D_NAD_1983_CORS96',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',177.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]];IsHighPrecision
  NO_CONFIGURELAYOUT
  "PROJCS['NAD_1983_StatePlane_Hawaii_1_FIPS_5101',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-155.5],PARAMETER['Scale_Factor',0.9999666666666667],PARAMETER['Latitude_Of_Origin',18.83333333333333],UNIT['Meter',1.0]];-5123000
  -12084900 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision"
  "GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400
  -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision" # # Start
  Time: Tue May 08 16:46:56 2012

  The area-of-interest is much larger than the grid's intended area.
  Processing time may be long. Failed to create features in GRD_Grid1.
  Failed to execute (MakeGridsAndGraticulesLayer).


texts: Routing on rail network with interim points (pgrouting)

texts:
I want to use pgrouting on a rail network, stored in PostGIS. For each service number I have all the station it stops at. Until now, I set up start and end nodes for the rail edges and created the topology in Postgis.

Unfortunately, the coordinates of the stations don't correspond with the positions of the vertices of the edges, although they overlap with the edges. What should be achieved in the end, is routing along the train lines, e.g.

  "show pathway of the line X, which stops at stations [5565,5619,7719]"

Because the underlying data changes more or less often and I need their original IDs, I want to avoid to change it too much by splitting the lines. Is there a way in pgrouting to use multiple interim points for calculating the shortest path?


texts: Adding raster Data from Geodatabase to QGIS

texts:
I have raster data in a geodatabase and want to open them with QGIS. When I navigate with QGIS to the Geodatabase I get to this window.

Unfortunately, I can't open These atx and GDBTABLX files and I don't recognize which of layers are the files. In ArcGIS the window just Shows the layer names.

Is there a way to recognize the layer name and open the files in QGIS?


texts:
I am using ArcGIS Web AppBuilder (Version 2.13 August 2019).

I have a widget, and on button click, I need to place a symbol at the centre of the map.

I tried using :

selectFeatureService: function (value) {
    value = Number(value);
    if (this.overviewLayer) {
        this.map.removeLayer(this.overviewLayer);
        this.overviewLayer = null;
    }
    this.graphicExists = false;
    //this.allCoverage.set("checked", false);
    this.map.graphics.clear();
    this.imageDialog.hide();
    ***var ctrPt = this.map.center;
    console.log("center Point found(X) :"+ctrPt.x);
    console.log("MinX: " + this.map.extent.XMin + "\n");
    console.log("MinY: " + this.map.extent.YMin + "\n");
    console.log("MaxX: " + this.map.extent.XMax + "\n");
    console.log("MaxY: " + this.map.extent.YMax + "\n");***
    if (value !== "select") {
        var url = this.config.oic[value].serviceUrl;
        if (url.indexOf("ImageServer") === -1) {
            var query = new Query();
            query.where = "1=1";
            query.returnGeometry = false;
            query.outSpatialReference = this.map.extent.spatialReference;
            var queryTask = new QueryTask(url);
            queryTask.executeForExtent(query, lang.hitch(this, function (response) {
                if (response.extent) {
                    this.map.setExtent(response.extent);
                }
            }));
        } else {
            var request = new esriRequest({
                url: url,
                content: {
                    f: "json"
                },
                "handleAs": "json",
                "callbackParamName": "callback"
            });
            request.then(lang.hitch(this, function (response) {
                if (response && response.extent) {
                    if (response.extent.spatialReference.wkid === this.map.spatialReference.wkid)
                        this.map.setExtent(new Extent(response.extent));
                    else {
                        var params = new ProjectParameters();
                        params.geometries = [new Extent(response.extent)];
                        params.outSR = new SpatialReference(this.map.spatialReference.wkid);
                        this.geometryService.project(params, lang.hitch(this, function (geometry) {
                            this.map.setExtent(geometry[0]);
                        }));
                    }
                }
            }));
        }
        if (this.config.oic[value].overviewUrl) {
            this.overviewLayer = this.layerModuleSelector(this.config.oic[value].overviewUrl);
            this.map.addLayer(this.overviewLayer);
        }
    }
},

Where am I going wrong?


texts: QGIS regular points around each polygon in a layer

texts:
I have a number of buffered points (somewhere polygons intersect & are dissolved.
I'd like to create regular spaced points grid around EACH polygon -- not ALL polygons.

When I do Vector -> Research Tools -> Regular Points --- I have to choose a layer as an Input Boundary, but it doesn't seem to loop through each shape to create a GRID just for that shape, and instead creates one large grid around ALL shapes. The end result:

How can I get a point grid for EACH individual shape vs ALL?


texts: ogr2ogr *.tab to postgres: Failng on date format

texts:
I'm importing some MapInfo tab files into a PostgreSQL database running the following command in FWTools Shell:

ogr2ogr -f "PostgreSQL" PG:"dbname='<the_database>' host='<the_host>' port='5432' user='<the_username>' password='<the_password>'" -lco SCHEMA=<the_schema> -lco PRECISION=NO -nlt GEOMETRY -t_srs EPSG:25832 C:\<the_local_path>\file.tab

Some of the files run through just fine, but some of the files have date values like that 0000/00/00  0:00:00 which ogr2ogr, justifiably, cannot interpret.

ERROR 1: INSERT command for new feature failed.
ERROR:  date/time field value out of range: "0000/00/00  0:00:00"
LINE 1: ...ning', 'Tilslutning af stik', 'I brug/drift', '', '0000/00/0...

Is there a flag which will make ogr2ogr ignore those, and just insert a null?


texts: Extracting data in region with Earth Engine

texts:
I need to extract 256x256 tiles from Landsat images using Google Earth Engine. Usually,
this is done by the following function calls:
image_tiles = my_image.neighborhoodToArray(kernel_size) #kernel size is 256x256 in this case
sample = image_tiles.sample(
                     region=my_roi.geometry(),
                     scale=30,
                     numPixels=1,
                     tileScale=8
)
# sample export omitted for brevity

However, contrary to the usual use case, I need to 256x256 tiles to be sampled at a specific
location. The code above samples the 256x256 tile from a random point within the ROI, but I need to know exactly where the tile is sampled. To this end, I've created a grid across my entire study area. Each grid cell (in an equal area projection) is size 256 * 30 x 256 * 30 (30m is Landsat resolution, and 256 is the tile size I want to export). When I call
task = ee.batch.Export.image.toCloudStorage(
        image=my_image,
        bucket=my_bucket,
        description='myDescription',
        fileNamePrefix='myPrefix',
        fileFormat='TFRecord',
        region=my_roi.geometry(),
        crs='EPSG:5070',
        scale=30,
        formatOptions={'patchDimensions':256,
                       'compressed':True,
                       'maskedThreshold':0.99}
        )

The image won't export, as a 256x256 tile doesn't fit within the my_roi's geometry. I'm not sure
why this is happening; I make the ROI so it should be exactly 256x256 pixels at Landsat scale.
Is there a simple way to export 256x256 tiles at the exact location I'm interested in?


texts: PostgreSQL table corrupt not display in QGIS

texts:
I have a problem with a PostgreSQL table. It is the result of a query that joins 77 tables with 4 columns (geom, inicio, fuente, temporada). I used ST_SRID and ST_Transform because I could not join tables with different SRID.

When I want to export it as shapefile QGIS can not display the file, it does not recognize it. I tried to insert a column ID_0 but the problem continues.

Consulting on this table resulting some sub tables and I could export as a shape, but not others. I do not know if the table is missing some metadata, or did something wrong with the reprojection.


texts:
I encounter the problem more and more often and it begins seriously to become embarrassing :
I open a project of yesterday and it is empty or I just have the title groups, groups emptied of their layers.
And the layers of the project are shp and virtual layers resulting from queries.
So
The .qgs has hundreds of kilos and when I edit it, I find my layers in the  tags, the paths are good, I have no messages like what tags are badly closed and file is not truncated.
Is project finally lost ?


texts: How to populate a div with hover text in bokeh?

texts:
I am developing a geographical dashboard in bokeh. It is fairly easy to use hover tooltip in bokeh, however I want to populate a div when I hover on a field.

Python version : 3.5.2
Bokeh version : 0.13.0

1. This is the main view of the dashboard.

2. Required output

I want to add two interactions,

When I hover on a district, it should only populate the district name field.
When I click on a district, it should populate both the district name and info field(area, population, density).

3. Data:

4. Code:
    from bokeh.io import show, output_notebook, output_file
    from bokeh.models import (
        GeoJSONDataSource,
        HoverTool,
        LinearColorMapper,
        ColorBar,
        BasicTicker,
        PrintfTickFormatter,
        LogColorMapper,
        Range1d,
        Plot,
        Text
    )
    from bokeh.plotting import figure
    import geopandas as gpd
        with open('/home/drogon/Desktop/Rescue-1122-project/punjab_districts(area_pop_den).geojson', 'r') as f:
            geo_source = GeoJSONDataSource(geojson=f.read())

        df = gpd.read_file('/home/drogon/Desktop/Rescue-1122-project/punjab_districts(area_pop_den).geojson')
        print(df.density)
        density = df['density']

        colors = ['#000003', '#3B0F6F', '#8C2980', '#DD4968', '#FD9F6C']
        colors.reverse()
        color_mapper = LogColorMapper(palette=colors, low=density.min(), high=density.max())

        TOOLS = "pan,wheel_zoom,box_zoom,reset,hover,save"

        p = figure(title="Punjab Districts", tools=TOOLS, x_axis_location=None, y_axis_location=None, width=800, height=800)
        p.grid.grid_line_color = None
        p.title.text_font_size = '30pt'

        p.patches('xs', 'ys', fill_alpha=0.9, fill_color={'field': 'density', 'transform': color_mapper},
                  line_color='white', line_width=1, source=geo_source)

        hover = p.select_one(HoverTool)
        hover.point_policy = "follow_mouse"
        hover.tooltips = [("District", "@district"),
                          ("Density", "@density"),
                          ("Area", "@area")]

        color_bar = ColorBar(color_mapper=color_mapper, major_label_text_font_size="10pt",
                             ticker=BasicTicker(desired_num_ticks=8),
                             formatter=PrintfTickFormatter(format="%d"),
                             label_standoff=10, border_line_color=None, location=(0, 0))
p.add_layout(color_bar, 'right')
show(p)


texts: I am running an IDW on a data set that represents point values of water current speed. I would like to run the IDW only within the water area of another shapefile that shows the shape of the water body. So far I have only been able to keep the IDW to the XY extents (a square) of the data points, not within the waterbody shape. Any ideas about how to do this would be very helpful


texts: Is there a way to interpolate lost data/lost pixels caused by cloud and cloud shadow masks in Google Earth Engine?

texts:
I'm sure cloud and cloud shadows have been, time and time again, a frustrating problem. Usually after masking the two, I composite multiple images (in a set time interval e.g. 2015-2016, months 6-8) and that solves the missing pixels caused by the masking.

Sometimes though it doesn't quite work.

Is there a way to interpolate lost data/lost pixels caused by cloud and cloud shadow masking in Google Earth Engine?


texts: Relative paths in ArcScene document are broken if opened programmatically

texts:
When I open ArcScene document with relative paths programmatically from ArcMap, all those paths are broken. The location in source tab for broken layers is set to some different path like "C:\Windows\System32" or "C:\Windows".

Below is the code sample I use to open document. Both ArcMap and ArcScene projects are located in "C:\test" and so are files added to ArcScene doc.

Sub test_relpaths()
  Dim pDoc As IDocument
  Set pDoc = New sxDocument
  Dim iApp As IApplication
  Set iApp = pDoc.Parent
  iApp.OpenDocument ("C:\test\test_relpaths.sxd")
End Sub

When I use similar code to open ArcMap doc from ArcScene, everything is fine, relative links are kept:

Sub test_relpaths()
  Dim pDoc As IDocument
  Set pDoc = New mxDocument
  Dim iApp As IApplication
  Set iApp = pDoc.Parent
  iApp.OpenDocument ("C:\test\test_relpaths.mxd")
End Sub

Can anyone please help me to fix this issue or tell if it's a bug.

Update: I've installed 9.3.1 with SP2 and some patches, but it didn't help.
Update2: I've also tested these samples in ArcGIS Desktop 10 SP3 and results are the same. It seems like an ArcScene bug.


texts:
Say I have a latitude of 38.802610 and a longitude of -116.419389 and I want to be able to add a distance to that, 20 miles for example.

So far I have came up with a rough idea on how I would do this:

Firstly, work out how many miles are in 1° of latitude, say for example it was 30
Then, divide one by it:

1 / 30 = 0.033333

Add it to my original latitude to get my maximum ° of latitude:

38.802610 + 0.033333 = 38.8355943

Subtract it to my original latitude to get my minimum ° of latitude:

38.802610 - 0.033333 = 38.769277

But this is flawed because there seems to be no direct conversion for longitude as from what I've read the calculation varies. I have looked around and found some resources but I can't find one that applies directly to my question. I came across Haversine Formula but I can't find a way of applying it to my situation as, from what I've read, it's used for calculating the distance between two points.

Ultimately, I need to be able to find out the:

maximum latitude (my current latitude + given distance e.g 20 miles)
minimum latitude (my current latitude - given distance e.g 20 miles)
maximum longitude (my current longitude + given distance e.g 20 miles)
minimum longitude (my current longitude - given distance e.g 20 miles)


texts: Creating fishnet grid polygon layer from set of ortho rasters in ArcGIS for Desktop?

texts:
I have a set of ortho slides (SID format) for my city. The slides do not overlap, and cover the entire city and immediate environs. What I want to do is create a grid-like polygon layer with separate polygons for each raster slide. All are equal in size (6k feet by 4k feet), so I would assume that this would involve creating a fishnet, right? The problem I am encountering is how to create the polygon in the first place. Would I create a larger polygon that would encase all the slides, and then create the fishnet from there?

Working with ArcMap Desktop v10.3.


texts: What kind of datasets are necessary for geomarketing applications?

texts: If you would have to create a data base for geomarketing applications yourself, what kind of datasets would you purchase/collect? The first goal would be to find new store locations. Another would be to assess the efficiency of advertising efforts in different areas.


texts: Choropleth of calculated Voronoi polygons?

texts:
I have calculated a voronoi layer for a set of points, I then joined said layer with counts so (almost) each point that generated that polygon has a number associated to it.

My question is: how do I color each polygon according to its count in QGIS?

My data is simply the voronoi layer (I don't exactly know where it gets the data for drawing the polygons), an ID for each point that generated that polygon, and a count associated to it. Like this:

   ID, COUNT
 ABC1,    12
 ZXI6,   234
 ...

But that's all I have. In any case, each ID has a lat/lon pair, again the point that generated the polygon.

I'm sorry I cannot upload a raw screencapture of the table, the data is sensitive.


texts: Error while combining multiple partially overlapping rasters into a single raster in R

texts:
I have multiple (17) raster files which I want to process. They are from the Sentinel 2 satellite System (ESA) and already transformed in .tif format.

I wanted to build a mosaic image of all these 17 images (Portugal) for each Band I need (Blue, Green, Red, NIR, SWIR). After a short search I found following program:
Combine multiple partially overlapping rasters into a single raster in R

When I try:

national_layer <- mosaicList(raster_files )

I get error:

  Error in (function (classes, fdef, mtable)  :    unable to find an
  inherited method for function ‘raster’ for signature ‘"numeric"’

File specifications:

class       : RasterLayer
dimensions  : 10980, 10980, 120560400  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 699960, 809760, 4590240, 4700040  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\###########\Desktop\winter\T29TQG\T29TQG_B04.tif
names       : T29TQG_B04
values      : 0, 23739  (min, max)

So if I crop the raster files manually it works, but would be nice if I use a loop.

complete code:

mosaicList <- function(rasList){

  #Internal function to make a list of raster objects from list of files.
  ListRasters <- function(list_names) {
    raster_list <- list() # initialise the list of rasters
    for (i in 1:(length(list_names))){
      grd_name <- list_names[i] # list_names contains all the names of the images in .grd format
      raster_file <- raster::raster(grd_name)
      raster_file <- projectRaster(raster_file, snap, method = "ngb")
    }
    raster_list <- append(raster_list, raster_file) # update raster_list at each iteration
  }

  #convert every raster path to a raster object and create list of the results
  raster.list <-sapply(rasList, FUN = ListRasters)

  # edit settings of the raster list for use in do.call and mosaic
  names(raster.list) <- NULL
  #####This function deals with overlapping areas
  raster.list$fun <- sum
  #raster.list$tolerance <- 0.1

  #run do call to implement mosaic over the list of raster objects.
  mos <- do.call(raster::mosaic, raster.list)

  #set crs of output
  crs(mos) <- crs(x = raster(rasList[1]))
  return(mos)
}

setwd("C:/Users/sebastian.r/Desktop/winter")
raster_files <- lapply(list.files(pattern = "02.tif$", recursive = TRUE), raster, values=TRUE, full.names=TRUE,extent.values=TRUE, mean=TRUE)
snap <- raster(resolution = c(10,10), xmn = 399960, xmx = 809760, ymn = 4090200, ymx = 4700040, crs = "+proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
national_layer <- mosaicList(raster_files )


texts:
I am using a NUTITEQ API key and now I want to use CARTO mobile in my App.

It says the app bundle is already used, but this is via my NUTITEQ account which I think is probably now inactive.

Can I transport the NUTITEQ account to CARTO?


texts: When I query, for example, the Germany Boundary Commission, and the Austrian Boundary Commission, for boundary information, I am getting conflicting data sets. Is there a resolution method for this?


texts:
I am looking for a way to get a map from QGIS into QT Designer. I have the QT Designer version that comes with QGIS, meaning i have all the standard QGIS Widget. I don´t see a map widget, but many things belonging to a map e.g. a mapScaleWidget.

How can I get a few functions and a map into QT?

I can programm a little bit of python, but I have not yet made any plugins myself, because I personally think I can´t programm good enough for that yet.

I also do not know how to install plugins.


texts: Changing layer order of vector tile layer

texts:
I am creating a simple web map using vector tiles of type pbf.
The vector tiles are created with tippecanoe like:
tippecanoe -pC -P -z 19 -Z 15 -e builds --drop-densest-as-needed --extend-zooms-if-still-dropping piers.GeoJSON ... water.GeoJSON ... green.GeoJSON w --force

When I display the layers in leaflet using vectorgrid, I can see, that the water polygons are covering the piers.
The leaflet example code goes like this:
L.vectorGrid.protobuf('http://localhost:8080/{z}/{x}/{y}.pbf', {
    maxZoom: 19,
    maxNativeZoom: 16,
    vectorTileLayerStyles: {
    greenGeoJSON: {
        color: "#d4dad6",
        weight: 1,
        fillColor: "#ADD19E",
        opacity: 1,
        fillOpacity: 1,
        fill: true,
    },

    waterGeoJSON: function(properties, zoom) {
        var lineWidth = 0;
        if (properties.area > 10000000) {
            if (zoom == 15) {
                lineWidth = 3;
            }
            if (zoom == 16) {
                lineWidth = 3.25;
            }
            if (zoom == 17) {
                lineWidth = 3.5;
            }
            if (zoom >= 18) {
                lineWidth = 4;
            }
        }
        return ({
            weight: lineWidth,
            fillColor: "#99ccff",
            color: "#bfc7c8",
            opacity: 0.1,
            fill: true,
            fillOpacity: 1
            });
        },
        pierpolygonsGeoJSON: {
            weight: 1,
            fillColor: "#e3e3dc",
            color: "black",
            opacity: 1,
            fill: true,
            fillOpacity: 1
        }
    }
}).addTo(mymap);

Here is the current layer order:

and this is the same map without the water layer:

How can I reorder the layers, that are coming from a pbf? Can I change the pbf creation process so the water-layer will be at the right "zlevel"?


texts: Map display not filling window in QGIS on Mac?

texts:
In this screenshot of QGIS, you can see that the map is not filling the window. I have tried resizing the window in various ways.

Am I missing something obvious, or is this a bug?

Versions:

MacBook Pro MacOS High Sierra 10.13.5
QGIS version 3.2.1-Bonn
QGIS code branch Release 3.2
Compiled against Qt 5.9.6
Running against Qt 5.9.6
Compiled against GDAL/OGR 2.2.4
Running against GDAL/OGR 2.2.4
Compiled against GEOS 3.6.2-CAPI-1.10.2
Running against GEOS 3.6.2-CAPI-1.10.2 4d2925d6
PostgreSQL Client Version 9.6.6
SpatiaLite Version 4.4.0-RC0
QWT Version 6.1.3
QScintilla2 Version 2.10.7
PROJ.4 Version 493
Python 3.6.6


texts:
I want to select the features that are not equal. So, I made a selection by location where the geometric predicate is equal and then I thought to invert the selection. Is there any way to invert the selection using PyQGIS?

My code:

import processing
layer1 = processing.getObject('polygon1')
layer2 = processing.getObject('polygon2')

processing.runalg('qgis:selectbylocation', layer2, layer1, u'equals', 0.005,0)

#INVERT SELECTION


texts: How to change the Transparency of an ILayer?

texts:
I've created a function which sets the IColor to a specific ILayer, the value of transparency is between 0 and 255. When I set the value to 0, it gives me a transparent layer which is OK, but when I set it to 1,2,3,4,5,....255 it gives me the color as if it has 255 transparency. I don't understand, what's wrong with my code:

public static void setLayerColor(ILayer layer, IColor fillColor, IColor outlineColor, short width = 1)
    {if (layer == null || fillColor == null || outlineColor == null) return;
     IGeoFeatureLayer geoFeatureLayer = (IGeoFeatureLayer)layer;
     ISimpleRenderer simpleRenderer = (ISimpleRenderer)geoFeatureLayer.Renderer;

     ISimpleFillSymbol symbol = new SimpleFillSymbolClass();
     ILineSymbol lineSymbol = new SimpleLineSymbolClass();

     lineSymbol.Width = width;
     lineSymbol.Color = outlineColor;

     symbol.Color = fillColor;
     symbol.Outline = lineSymbol;
     simpleRenderer.Symbol = (ISymbol)symbol;}

My function has 4 parameters:

ILayer, the layer of which I want to change its color.
IColor, fillColor.
IColor, outlineColor.
short, the width of the outline.


texts: ArcPy saveACopy method saving copy of MXD to wrong path?

texts:
I have a script that looks through datasets in a geodatabase and copies a mxd document to a folder. The programs works, however it is saving the map documents to my desktop and not the specified folder on my desktop (PATH2).
This is the code I have so far:

import arcpy
import os
PATH1 = r"C:\Users\Daimon Nurse\Desktop\DFMPROJECT\DFM.gdb"
PATH2 = r"C:\Users\Daimon Nurse\Desktop\DFMPROJECT\"
arcpy.env.workspace = PATH1
arcpy.env.overwriteOutput=True

zones = arcpy.ListDatasets("Zone*")

for zone in zones:
    mxd = arcpy.mapping.MapDocument(r"C:\Users\Daimon Nurse\Desktop\DFMPROJECT\Complete Final Capstone Map\Complete Z2 Campus Map_Landscape1.mxd")
    print zone
    mxd.saveACopy(PATH2 + zone)


texts: Way to set a label leader line length (not using annotation) in ArcGIS Desktop?

texts:
How to set a label leader line length without using annotation in ArcMap 10.5?

I am using label leaders mainly for labelling features - pointing at the polygon/line/point that the label refers to.

I use Labelling Toolbar - Label Manager, than go to Symbol (/Edit Symbol/3rd tab - Advanced Text/Text Background - Properties) .... Line Callout/Leader Symbol (/Edit Symbol/Cartographic Line Symbol) ..., using simple line or arrow line mainly

example:

In this example I have set an offset that determine how far the label from the polygon border would be, but I cannot set internal zones (going from the dialog on picture bellow to properties and on the 1st tab internal zones) to influence where the leader will have it's end... which would I think answer that question - so I am looking for another way to do it.

illustration


texts: Copying raster data from QgsRasterDataProvider to OpenCV matrix?

texts:
I'm using C++ (QtCreator) QGIS 2.0 und openCV2 on an 64 bit Workstation running Debian Wheezy. The Application I try to assemble is a standalone QGIS one and will be a mix of gis and image processing software.

Problem description:

I've a QgsMapLayerCanvas with an associated QgsRasterlayer based on an GeoTiff with UTM32 CRS and a resolution of 2 cm. My objective ist to point on an location, determin the pixel data and cut out the data buffer in a 2 x 2 meter window . In the next step this buffer should go to an opencv image processing routine. I' havent found a ready to use solution in the qgis library.

My guess is, that the QgsRasterDataProvider can do this job but. But I cannot find any useful snippets in the web for this task.

Is there somebody who has experience with this matter?

Here the snippet of my first try:

void MainWindow::doMapImgCanvasClicked(const QgsPoint &point,
                                  Qt::MouseButton button) {

    if (!(qgsMapImgLayer && qgsMapImgProvider) ) return;

    // Size of the window for the image process
    int winWidth  = 10; int winHeight = 10;

    // Calculate the pixel position
    QPoint pixPos(-1,-1);
    doMapImgCalcRasterPos(point, qgsMapImgLayer, pixPos);

    // Give me a response for coords and pixel pos
    QString str;
    str.sprintf("COORD %f %f PIX %d %d",
                       point.x(), point.y(),
           pixPos.x(), pixPos.y());
    outLog(str);

    // Create a region of interest
    QgsRectangle roi(point.x()-1, point.y()-1, point.x()+1, point.y()+1);
                // Make a copy of Data 2x2 m for the first band
    QgsRasterBlock *block = qgsMapImgProvider->block(
      1, roi, winWidth, winHeight);

     // Block data is OK?
    if (block->isValid()) {
        outLog("VALID");
    }

                // Blockdata empty?
    if (block->isEmpty()) {
      outLog("EMPTY");
    } else {
     return;
    }

    // Is this the buffer?
    char* bits = block->bits();
    str.sprintf("DB.TYPE %d", block->dataType()); outLog(str);

    // OK This is the size of the Pointer
    str.sprintf("DB.NUM %d", sizeof(bits) ); outLog(str);

    // Test Output
    for (int i=0; i<winWidth*winHeight; i++) {
            char c = bits[i];
        str.sprintf("PIX[%d]=%d", i, bits[i] ); outLog(str);
    }
}


texts: Create different CSS for Two Buttons in One Control Panel

texts:
I am having a problem applying CSS to two custom OpenLayers buttons. I create a PANEL toolbar and then add 8 buttons to it. The ZoomBox, DragPan, ZoomOut, ZoomIn, and nav.next and nav.previous buttons all style the way they should. The issue is with the two BUTTON items.

As near as I can tell the OpenLayers BUTTON object only allows one CSS style called .olControlButtonItemInactive.

THE PROBLEM: I have two BUTTON objects in the PANEL and I want each one to have a different style. But, if I change the name of the CSS class from olControlButtonItemInactive the CSS no longer is picked up by the button.

I thought I could just assign a unique ID attribute to each BUTTON (for instance button #1 has an ID of 'btnMaxExtent' and button #2 has an ID of 'btnHiLite'), then create a CSS style for each one called, respectively, olControlbtnMaxExtentItemInactive and olControlbtnHiLiteItemInactive and then assign these two button styles to the olControltoolPanel class but it does not work.

Here is the script I use to create the buttons on my panel:

    var panelControls = [
            new OpenLayers.Control.ZoomBox({title: "Draw a box to ZOOM IN", text: 'Click Me'}),
            new OpenLayers.Control.DragPan({title: "Click to PAN around map"}),
            new OpenLayers.Control.ZoomOut({title: "Click to ZOOM OUT one level"}),
            new OpenLayers.Control.ZoomIn({title: "Click to ZOOM IN one level"}),
            new OpenLayers.Control.Button({title: "Click to go to MAX EXTENT",
                id: 'btnMaxExtent',
                trigger: zoomToExtent}),
            nav.next, nav.previous,
            new OpenLayers.Control.Button({title: "Click once for hover info, click it again to turn off the functionality",
                id: 'btnHiLite',
                trigger: activateHiLite})
            ];

Just for completeness here is the script I use to add the controls to my panel:

    var toolbar = new OpenLayers.Control.Panel({'div':  document.getElementById('tlbr')});

    toolbar.addControls(panelControls);
    map.addControl(toolbar);

And here is the CSS I have now that works for the BUTTON objects but it makes both BUTTONs look the same.

I want each button to have its own CSS. If I take this same CSS script and create one for each of the BUTTONS with the names .olControlbtnMaxExtentItemInactive and .olControlbtnHiLiteItemInactive they do not work.

  .olControltoolPanel .olControlButtonItemInactive {
border-top: 1px solid #96d1f8;
background: #6ed665;
background: -webkit-gradient(linear, left top, left bottom, from(#3e9c50), to(#6ed665));
background: -webkit-linear-gradient(top, #3e9c50, #6ed665);
background: -moz-linear-gradient(top, #3e9c50, #6ed665);
 background: -ms-linear-gradient(top, #3e9c50, #6ed665);
background: -o-linear-gradient(top, #3e9c50, #6ed665);
padding: 5px 10px;
-webkit-border-radius: 8px;
-moz-border-radius: 8px;
border-radius: 8px;
-webkit-box-shadow: rgba(0,0,0,1) 0 1px 0;
-moz-box-shadow: rgba(0,0,0,1) 0 1px 0;
box-shadow: rgba(0,0,0,1) 0 1px 0;
text-shadow: rgba(0,0,0,.4) 0 1px 0;
color: #9c494e;
font-size: 12px;
font-family: Georgia, serif;
text-decoration: none;
vertical-align: middle;
background-image: url('Img/zoom-world-mini.png');
background-repeat: no-repeat;
    }

Any suggestions would be appreciated. Thanks for taking the time to look at this!


texts:
Using QGIS 3.10 - new data points are disappearing after clicking on save. Layer is digitising points (park assets) which have been identified from field capture. Up until yesterday April 17, 2020, all data points were appearing, and then after changing the style (colour) for one asset type, all have now disappeared.

The point appears in the attribute table, but not symbolized on the map. I am using categories for my symbols.

Example - Each of the grey squares should have a point symbol in the centre - as per the Attribute table. Note - disregard the brown dots this is our field capture layer

I have tried following above responses to saving a new shapefile - this has not changed the issue.

Why are all points not rendering?

re @she-weeds screenshot of styling details - attached, these have all been "classified" with symbols for each asset type which relate to the park asset type


texts: GIS data for Grocery store sales in the US

texts: I am looking for sources that would provide sales data for grocery stores in the US. I have looked at ESRI's Business Data that they receive from InfoGroup, but was wondering if anyone know of any other sources.


texts: How to quantify projection distortion using ArcGIS

texts: I am trying to figure out a workflow to quantify the amount of distortion for different projections throughout a given area of interest.  I am thinking in the form that you would do a relative accuracy assessment.  The difference between the known (direction, distance, area, etc) vs the projected.  I envision the area of interest covered in a grid of points, lines and polygons (or whatever combination of vector i need for the calcs) and calculating the values both ways.  Is there something like that out there?  What would be the steps needed to take for each distortion?  I would probably only care about distance, area, and direction.


texts:
I am trying to join three data files to a raster grid:

I have environmental data in the centre of each grid cell (delimited text layer), Sightings (delimited text layer), effort points (delimited text layer), and track lines (points2one). The grid is a raster layer.
How can I join the data which is within each grid cell to form a single data set to export?


texts: FME: why do my readers sometimes not get added to the canvas?

texts: Sometimes when I generate a workspace in fme 2014, the reader and writer get created but not added to the workspace. Can anyone explain why?


texts:
I am trying to flatten an xml file with fme 2014. I´m having problems understanding how to manage the following xml schema. One of the XML child elements is called "aspects" and contains the following ...

<aspects>
<a id="1" l="de" v="Alter Schafstall und Eichen beim Südsee-Camp"/>

How do I go about flattening the attributes? FME seems to flatten most of the table, but the above gets defined as an "XML Fragment". In the table the field "XML_Fragment" then contains the above xml structure with tags etc...ie is not being flattened.


texts: ArcGIS 10.3 Address Locator - Partial Searches Not Working

texts:
I'm using an address locator service based off an address point feature class in ArcGIS 10.3.  My ultimate goal is to write an app which will parse a list of addresses in a database and determine the X/Y coordinates for them.

So far the service is working correctly, however it is not returning results when search terms are only partially entered.

For example:

"2150 Main St" returns data
"2150 Main" returns data
"2150 Ma" does not return data

I have tried searching using wildcards, but that has not made a difference.  I know that ArcGIS 10.3 has the option of doing a Suggest operation, but that will not work since all of these searches will be happening behind the scenes.

Does anyone know how I can enable the locator to work with partial search terms?


texts: What's the time zone of the timestamp of LANDSAT Images?

texts:
I am using google earth engine just to print the timestamp here of LANDSAT Image. I am using the geometry of a part of Chennai which is in India.

// Filtering Data with Date and Geometry
var Data2 = LANDSAT.filterBounds(geometry).filterDate('2013-04-09', '2014-04-10');
// Printing Timestamp
print('LANDSAT Metadata', Data2.first().date());
Map.centerObject(geometry, 8);
Map.addLayer(geometry,{color: 'FF0000'},"Chennai");

My question what exactly the printed time does indicate (with respect to UTC (GMT) or Local Time Zone (IST) ).

just in case if you need the import

IMPORTS:

var LANDSAT = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA"),
geometry = /* color: #d63000 */ee.Geometry.Polygon(
    [[[80.26491619216495, 13.08801893465061],
      [80.26422954665713, 13.077652181027775],
      [80.28002239333682, 13.08166516962143]]]);

I am using google earth engine only to print the timestamp here of LANDSAT Image and not anyway related to my question.


texts: Getting H2GIS functions to work in Spherical Space?

texts:
I've been playing around with H2GIS (http://www.h2gis.org/) for a project at work and am struggling to get the H2 functions to operate in spherical space (lat/ lon), rather than cartesian space.

The problem is basically this. I want to store points in a geometry column of an H2 table using the WKT format, where the coordinates represent a lat/ lon pair (i.e: 'POINT (0 0)' would be the point where the prime-meridian crosses the equator).

So I create my table and insert a point into it like:

CREATE TABLE GeoTable(location GEOMETRY);
INSERT INTO GeoTable (location) VALUES ('POINT (45 45)');

I did some calculations by hand to figure out how far 'POINT (45 45)' would be from 'POINT (0 0)' in spherical and cartesian space using the Haversine and Euclidean distance methods respectively to devise a test that should fail in Cartesian space but succeed in Spherical space. The numbers I used work out to a Cartesian distance of 7076.77 km (assuming a distance of 111.2km between degrees of latitude, which I understand is only accurate at the equator). The distance computed using the Haversine formula (calculated using http://andrew.hedges.name/experiments/haversine/) is roughly 6673.79 km. When measuring the distance on Google Earth, the haversine distance is within 20km, however the cartesian distance is more than 400km off, so given a difference of 20x in the error over large distances we want to make sure predicates supported by H2GIS can operate in Spherical space.

The test I've devised to confirm we can use H2GIS predicates in spherical space uses the ST_DWithin method to select all points within 6700km of 'POINT (0 0)', which should return 'POINT (45 45)' in Spherical space, but not in Cartesian space.

CREATE TABLE IF NOT EXISTS GeoTable(location GEOMETRY);
INSERT INTO GeoTable (location) VALUES ('POINT (45 45)');
SELECT * FROM GeoTable WHERE ST_DWITHIN(ST_GeomFromText('POINT (0 0)', 4326), ST_GeomFromText(location, 4326), 63);

If that final select returns 'POINT (45 45)', then I know that the distance between the points is being calculated in spherical space, as the required distance in cartesian space is > 69 degrees (calculated by dividing the cartesian distance in km between the points by the number of km per degree, again using 111.2km/degree). If we are in spherical space, then the distance should be ~ 60 degrees. I can't get this test to pass, regardless of how many ST_Transforms into WGS-84 (using SRID=4326) I apply to location or the input POINT, which leads (finally) to my question: How can I modify/ transform geometry values in H2 so H2GIS predicates function treat coordinate pairs as lat/ lon pairs and compute distances in spherical space?

I've seen that this is supported by tools like PostGIS (with a geography type that is differentiated from geometry types). Is this something that H2GIS simply doesn't support?


texts: Using database connections in ArcPy with ArcGIS Pro (arcgispro-py3)?

texts:
Since Pro stores database connections in a project, how do we connect to a database from within a standalone python script?

For example below is how I would have previously connected in ArcGIS 10.5

import arcpy

GDB_directory = r"C:\Python Testing"
targetEntGDB = r"Database Connections/espr1sql14.sde/"
extractedGDB = "600150be54a84779a461d4fc37247f91.gdb"
extractedGDB = GDB_directory + "\\" +extractedGDB

arcpy.Copy_management(extractedGDB, targetEntGDB)


texts: Why does my hillshade .ovr display in lower resolution in QGIS?

texts: I have viewed a lidar terrain model hillshade .ovr file on QGIS and Arc and Arc has significantly higher resolution than QGIS.  Any thoughts as to why this would happen and ways to enhance resolution in QGIS?


texts: OpenLayers getting OSM from local server

texts:
I have successfully installed OSM tile server on a local server. I verified operation using direct calls

http://server_ip/hot/z/x/y.png

I setup a page with OpenLayers 5.2.0. The layer init :

osm = new ol.layer.Tile({
    source: new ol.source.OSM({
        opaque: false,
        url: 'http://server_ip/hot/{z}/{x}/{y}.png'
        //url: 'http://a.tile.openstreetmap.org/{z}/{x}/{y}.png'
    })
});

When I use the a.tile.openstreetmap.org it works ok, but shows blank tiles on my local server call.

I verified the call to my server using wireshark, http request and reply. The return of HTTP are the required tiles (verified by saving them to disk). But OpenLayers is showing blank.

Does OpenLayers expect other parameters (in layer init code) when not specifically using osm server? (or timeout problem since my server is slower?)


texts: Best Bluetooth GPS for marine applications

texts:
What is the best GPS to get as an add on for an Android (ASUS tf300)? This is for a marine application I am writing.

I need sub 5 meter accuracy.
Would like decent accuracy at moderate speeds.
I don't want to spend more than 120 USD.
Must be Bluetooth since device doesn't have USB

Battery life is important since I'll be on a boat for 12+ hours at a time - the GPS will be in standby mode some of this time - battery should be able to run at least 6 hours in active mode without charging

I'm assuming that most of the devices on the market today are waas enabled (this is what I want right?). I've looked at some globalsat and dual models but not sure what is the best.

Please be gentle as I'm a noob when it comes to GPS technology.


texts: gdal_merge_raster.py ; Names of huge number of input images

texts:
gdal_merge_raster.py
  I am trying to make a stack including huge number of images. I am just wondering whether there is a way to avoid writing all the input images names !!!!!!

Thanks


texts:
In USA, I'm looking for outline of the Division 9 Pacific which consists of the following states:

Alaska
California
Hawaii
Oregon
Washington

I have got a map with the KML to generate for these states: https://drive.google.com/open?id=1lYAaPc-U3sv6hUlbhBsQJ0vzXfw&usp=sharing

But I just want the outline of the whole area and NOT the state borders. How can I get this?


texts: Calculate statistics on overlapping polygons with turf.js

texts:
In my Mapbox application, I have added draw tools to calculate the area of a polygon according to the example provided in the documentation.

var draw = new MapboxDraw({
    displayControlsDefault: true
});

map.addControl(draw);

var calcButton = document.getElementById('calculate');
calcButton.onclick = function() {
    var data = draw.getAll();
    if (data.features.length > 0) {
        var area = turf.area(data);
        // restrict to area to 2 decimal points
        var rounded_area = Math.round(area*100)/100;
        var answer = document.getElementById('calculated-area');
        answer.innerHTML = '
' + rounded_area/1000000 + '

square km
';
    } else {
        alert("Use the draw tools to draw a polygon!");
    }
};

After adding a polygon layer (census tracts) from a geojson, my objective is to now calculate the sum of population in all of the census tracts that are overlapped by the drawn polygon - for example:

It is possible using turf.js to find out whether features overlap with overlap(poly1, poly2), however according to the documentation this seems to only be possible only with  turf polygons as opposed to geojson polygons added as layers using the standard Mapbox library.

how would one read in a geojson file to turf and convert it to a multipolygon (with turf.multiPolygon()?) and then get specific properties (such as Pop2010) from those features that are overlapping?

The code below adds censusTracts, but how to get this layer to interact with turf.js?

map.addSource('censusTracts', {
    'type': 'geojson',
    'data': 'wgs84.geojson'
});

        map.addLayer({
            'id': 'density',
            'type': 'fill',
            'source': 'censusTracts',
            "paint": {
              "fill-color" : {

                "property": "Pop2010",
                "stops": [
                    [1, "#FFFFFF"],
                    [82200, "#404040"]
                ]
            },
            'fill-opacity': 0.9,
            'fill-outline-color': 'rgba(0, 0, 0, 1)'
          }

        });


texts: Select records that intersect more than 3 polygons

texts:
I am analyzing parcels in land use zones.

Some parcels overlap 3 or more zones = I am interested in selecting  only  those records (i.e. the purple-colored parcels)

At this point I have a query that joins zones in zoning to parcels in parcels when an intersection occurs

SELECT t.*, m.*
FROM parcels AS t
LEFT JOIN zoning AS m
ON ST_Intersects(t.geom, m.geom)

However, ST_Intersection() only returns true or false .  Is there a way to "count" the intersections and specify 3 or more?


texts: Using R for analysing relations between administrative borders (neighborhood of polygons) and voting results (or void votes)

texts:
What would be the relatively simple way to verify if there is a relationship between administrative borders (polygons)  and voting results (like: void votes). This would be used to verify if voting proces was tampered with.

I have:
 - shp file with administrative borders,
 - data with populations (eligible voters), vote results for particular parties, void votes within these polygons.

So far, by reading various sources I realize it is important to offset for the population in various  polygons and to account for neighboring polygons.

I am wondering what package and what model should be used for such analysis. Would GAM be appropriate? How to account for neighboring polygons as a variable?


texts: Added columns in SpatiaLite won't appear in QGIS (attribute table)

texts:
I have a SpatiaLite Database with OSM-Data from planet OSM (imported with the OSM import tool from QGIS 2.0)

Everything is fine so far. I can create, export, delete tables and draw them in QGIS... but when I create a new column in my table with:

ALTER TABLE planet_22_points
ADD xlabel integer DEFAULT 0

It will appear in the "DB Manager" and in "QspatiaLite" but NOT in QGIS! So I tried to add a column with QGIS using the attribute table. This actually seems to work, but after saving the table, zack! the columns I was creating before using the SQL command just appears! And that massage comes up:

Could not commit changes to layer planet_22_points

Errors: SUCCESS: 1 attribute(s) added.
ERROR: the count of fields is incorrect after addition/removal of fields!

QGIS 2.0.1, SQLite version: 3.7.9; SpatiaLite Library: 4.1.0


texts: Positioning Overview Map On ArcGIS API for Javascript

texts:
I am adding an overview Map to the main map by this way

  var overviewMap = new OverviewMap({
          map: map,
          visible: true
        });
        overviewMap.startup();

which the overview map is poping up at Top-Right

Can you please let me know how I can I re-position it to left-bottom


texts:
How to automatically load all my attribute values into a layer after manually adding an attribute field? Or how to achieve the effect of PointConnector plugin in a more flexible that retains attributes?

I have a point table

  City | Lat | Long

and a route table with an attribute

  From | To | Journeys

After loading both files in QGIS, I use the PointConnector plugin to generate lines connecting the cities, however the resulting layer lacks the Journeys Source Field (which I would use to determine a graduated line color in the Symbology pane).

I can add the Source Field manually, and enter each and every Journeys value in the Attribute Table of the layer. This becomes a challenge when I have 1000s of data points.

Can you advise a viable solution?

I suppose the PointConnector plugin doesn't support importing attributes. How can I automatically load all my Journeys attribute values after manually adding the attribute field?

How can I achieve the same effect offered by PointConnector (joining two points with a line) in another more flexible way that can make use of attributes in the import tables?


texts: Does anyone know where I can get data like this for north America (for use in SSRS). It does not have to be exploded. If I add Bing maps my 36 page PDF is 250MB when exported!


texts: Can arcpy modify Annotation TextString ?

texts:
Annotations that are linked to the source entity possess an attribute table within which you find a "TextString" column, containing the text displayed by the annotation.

If the content of the column is changed manually within ArcGIS Pro GUI, the displayed text changes accordingly.
But, when done with arcpy, the "TextString" cell content is updated but the text displayed in the map stays the same.
This problem is described here too :
https://community.esri.com/thread/206173-can-i-change-annotation-expression-using-arcpy

Has anyone worked his way around this ?

(Specifics of my case : I need them to be annotations and not labels, and I need to update the annotation content, and can not just get the text right from the label formula.)


texts: ArcGIS Server is truncating the WHERE clause when sending query to SQL Server

texts:
When I make a query request to a map layer with some WHERE clause to filter the response, the ArcGIS Server truncates the WHERE condition while executing the query on SQL Server.

I am dealing with very large data set of approx 600 millions. The ArcGIS version is 10.4.1

The map layer is bound with a SQL view.

I am querying map layer with following condition:

where: ReadRate30Days <= 100 AND ReadRate30Days >= 0
orderByFields: ReadRate30Days DESC

Ideally in SQL Sever the query should be genarted like below by ArcGIS Server

SELECT *
FROM OAWarehouse.Map_PM_Electric_Service_Point_View) a
WHERE ( (  ( ReadRate3Days <= 100 AND ReadRate3Days >= 0 ) ) AND ( ReportingDate >= '2019-10-09 00:00:00' and ReportingDate <= '2019-10-09 00:00:00' ))
  AND (Shape.STIntersects(geography::STGeomFromWKB(@P1,@P2)) = 1)
ORDER BY  ReadRate3Days DESC

but it is generated like this in SQL Server

SELECT * from (SELECT * FROM OAWarehouse.Map_PM_Electric_Service_Point_View) a

I am not able to understand this behaviour ?


texts: Change symbology of layer from "Single symbol" to "unique value" with python

texts:
How can I change the symbology of a layer with single symbol to unique values with Python.  I  don't know how to update it.

I used this code:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
lyr = arcpy.mapping.ListLayers(mxd,"region")[0]
if lyr.symbologyType == "UNIQUE_SYMBOLE":
  lyr.symbology.valueField = "REGION"
  lyr.symbology.addAllValues()
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd

The code runs but doesn't change anything (no error message)
the problem is that I want to change the symbology of the layer displaid with single symbol to category "uniques values"?


texts:
I'm trying to insert points using C# into a PostGIS / PostgreSQL database. I'm using the Npgsql driver. I suspect this is less a driver issue and more a syntax issue due to me being unfamiliar with ST_GEOMETRY.

Here's the code:

NpgsqlConnection conn = new NpgsqlConnection("Server=127.0.0.1;Port=5432;User Id=XXX;Password=XXX;Database=XXX;");

NpgsqlCommand command = new NpgsqlCommand("INSERT INTO points(name, point) VALUES(2, ST_GeomFromText('POINT(:longitude :latitude)', 4326));", conn);

var longitude = command.CreateParameter();
longitude.Direction = ParameterDirection.Input;
longitude.DbType = DbType.Double;
longitude.ParameterName = "longitude";
longitude.Value = fix.Longitude;
command.Parameters.Add(longitude);

var latitude = command.CreateParameter();
latitude.Direction = ParameterDirection.Input;
latitude.DbType = DbType.Double;
latitude.ParameterName = "latitude";
latitude.Value = fix.Latitude;
command.Parameters.Add(latitude);

conn.Open();

var dr = command.ExecuteNonQuery();

conn.Close();

When I run the code, I get a Npgsql.NpgsqlException with an error message ERROR: XX000: parse error - invalid geometry" The Hint property has a message "POINT(:l" <-- parse error at position 8 within geometry.

How do I construct a proper SQL statement to insert ST_GEOMETRY using parameters?


texts: Point to line tool to link all pairs of points in ArcGIS Pro

texts: I have a series of points which I want to connect together as a web, with lines between all possible pairs of points. I'm sure I've done this before using the points to line geoprocessing tool, but I'm struggling to make it work.


texts: QGIS crashes loading raster layer (.tif) through GDAL provider

texts:
I am working on a project that uses QGIS within another application using QT Creator.  When I try loading my raster data (using C++ API), I get a SIGSEGV error that traces back to my GDAL provider plugin library.
However, when testing the validity of my data just using the QGIS standalone application, it loads just fine.
I am using QGIS 3.8 and Qt 5.11.3 on a Fedora Linux OS (v28) with GDAL 2.4. I noticed several warnings related to my libproj.so version possibly conflicting with gdal and sqlite when building QGIS from source (Fedora 28 is EOL and I needed newer versions for the application I'm working on; OS choice is out of my hands, unfortunately) so it might be due to that.  But here is my main.cpp:
#include "mainwindow.h"
#include "qgsapplication.h"
#include "qgsproviderregistry.h"
#include <QApplication>
#include <QString>

int main(int argc, char *argv[])
{
    qputenv("QT_STYLE_OVERRIDE",0);
    QApplication a(argc, argv);
    QgsApplication b(argc, argv, true);
    QString qgis_pfx = "/home/me/qgis";
    QString gdal_path = "/usr/local/lib/gdalplugins";
    QgsProviderRegistry::instance(gdal_path);
    b.setPrefixPath(qgis_pfx,true);
    b.setPluginPath(qgis_pfx+"/lib");
    b.initQgis();
    MainWindow w;
    w.show();
    return a.exec();
    return b.exec();
    b.exitQgis();
}

and the snippet of code that segfaults in mainwindow.cpp between '#':
QString raster_path = "/home/me/data/raster_01.tif";
QString raster_name = "Surface Temperature (C)";
QString raster_provider = "gdal";
#QgsRasterLayer *r_l = new QgsRasterLayer(raster_path,raster_name,raster_provider);#

and the backtrace:
...received signal SIGSEGV, Segmentation fault.
0x00007fffc9bf0ef3 in int QAtomicOps<int>::load<int>(std::atomic<int> const&) ()
   from /home/me/qgis/lib/qgis/plugins/libgdalprovider.so

I've narrowed it down to either
a) bad/incomplete QGIS install or
b) improper use of the C++ API.
I have the needed libraries set in my LD_LIBRARY_PATH and in the .pro file along with needed include paths and have set a few other variables such as GDAL_DATA prior to launching the app.  The next debug step is to try an example of this in Python and see if that results in similar problems.  Will also show for the record that the raster does indeed load in QGIS GUI, so it's probably not the data itself:


texts: Getting smoother contours using ArcScene

texts:
I am trying to get nice smooth contours using ArcGIS. I have already smoothed out the DEM using focal statistics, simplified and smoothed the contours. The current result looks like this, which is way too pixel-y whenever the line is not straight:
And zoomed in on my peeve:
What I am going for is something less raster-y and more vector-looking like this:
I don't know if something like that is even possible within ArcGIS or whether I need a different application. Nevertheless I hope someone can guide me onto the right path towards perfect contours.


texts: How can I change vector grid extents in QGIS?

texts: I created many vector grids, but later I discovered that I need to modify some of their extents, hence, How can I change vector grid's extents (xmin, xmax, ymin and ymax) in QGIS?


texts: Using svg symbols for styling elements in ms4w (mapserver)?

texts:
I have MapServer installed in a windows machine. I have imported into my PostGIS db osm data and want to style certain elements (nodes) with svg symbols.

The problem is that the version of MapServer that I have (I am using ms4w) does not have precompiled the libsvg-cairo library which is needed in order to use SVG symbols for styling elements.

Does anyone know how to accomplish this?


texts: problem with resizing Polygons with gBuffer in R

texts:
I'm working on something and having a little problem.

library(tmaptools)
library(rgeos)
library(raster)
library(leaflet)
library(dplyr)

USA <- getData(country="USA",level=1)
Texas <- USA[USA@data$NAME_1=="Texas",]
texas_projected <- spTransform(Texas, CRS("+init=epsg:3857"))

places <- texas_projected %>%
spsample(50, type = "random") %>%
  # build random poisson buffers around them
  gBuffer(width = 4000*rpois(length(.), 4), byid = TRUE) %>%
  gEnvelope(byid = TRUE) %>% # take their envelope
  gUnaryUnion()

places <- spTransform(places, CRSobj = crs(Texas))

map<- leaflet() %>% addTiles() %>%
  addPolygons(data = Texas, weight = 1) %>%
  addPolygons(data = places, weight = 1,
              fillColor = "yellow",
              color = "darkblue")
map

places <- spTransform(places, CRS("+init=EPSG:32633"))
places <- gBuffer(places,width=5)
places <- spTransform(places, CRSobj = crs(Texas))

map2<- leaflet() %>% addTiles() %>%
  addPolygons(data = Texas, weight = 1) %>%
  addPolygons(data = places, weight = 1,
              fillColor = "yellow",
              color = "darkblue")
map2

map and map2 are identical, isn't that supposed to be a gain on the width of those green polygons when we plot map 2 (reproduce de code to better understandment)

why is gbuffer not resizing my polygons?


texts: I have a .dwg/.dxf file that needs to be georeferenced so I can draw some boundaries on a site in Indianapolis, Indiana, USA. I'm using GQIS 3.18.2, only the .dxf file will be imported. Now when imported, the data shows up as the map figure as a whole (reading the legend etc. as data). This is floating in space. Under layer properties, I can view the extent as (-362.3947145426504903,-998.1682020468679184 : 2891.9386173585958204,1472.8011332957680679). I've tried playing around with the "brute force" method but cannot seem to land on the right one.


texts: Dragging multiple polygons on Canvas rerender in React Leaflet

texts:
I am using React Leaflet to render Leaflet map and Leaflet.Path.Drag to drag multiple polygons together. Here is a demo.
I have a lot of polygons and would like to drag multiple polygons with preferCanvas option. If I try to drag multiple polygons with preferCanvas, it doesn't work. Dragging makes duplicate polygons and polygons are not visible during dragging.
codesandbox.io
Cross-post from SO.


texts:
I have the following PyQGIS script (see below) which works fine except from one issue: The input vector layer I would like to use is a temporary file. How can I call a temporary layer (which does not have a defined location) in PyQGIS?

##importing power plants vector layer
cpp = [call temporary layer here]

##importing population raster layer
pop = "C:/AIRPOLIM-ES/popKE.tif"

from qgis.analysis import QgsZonalStatistics

polygonLayer = QgsVectorLayer(cpp, 'zonepolygons', "ogr")
QgsProject.instance().addMapLayer(polygonLayer)

rasterLayer = QgsRasterLayer(pop, "Raster Layer")
QgsProject.instance().addMapLayer(rasterLayer)

zoneStat = QgsZonalStatistics (polygonLayer, rasterLayer , 'pre-', 1,QgsZonalStatistics.Sum)
zoneStat.calculateStatistics(None)


texts:
I have binary (0-1) rasters.

I'm trying to do an operation equivalent to:

z = arcpy.sa.Con((Raster(a) | Raster(b) | Raster(c) | ... | Raster(y))

where the raster files a-y are dynamically generated from a list. I tried using python's built-in 'any' function like so:

z = arcpy.sa.Con(any(Raster(r) for r in rasters_a_to_y))

but to no avail, I get an error Truth value of a raster is ambiguous.  Is there any equivalent to 'any' function in the ArcPy toolbox or another way to do this that I haven't considered?


texts: "User memory limit exceeded" using Random Forest classification

texts:
I am trying to classify using randomForest classification in GEE . However, I always get the error User memory limit exceeded, so I assume the number  of the feature are too large.
Nevertheless, I saw many developers can classify even bigger dataset, so I am wondering if there is a better option there?
The code can be accessed via
https://code.earthengine.google.com/23538cf3ed6b667d6e5975cb2864a436
var fire= ee.FeatureCollection('users/spatola/recovery_rf_1712').map(function(feature){
    var num = ee.Number.parse(feature.get('CLC_311'));
    return feature.set('CLC_311', num);
  })
var indices= ['EVI_mean','Dnbr','NBR2_mean','NBR2_min','NBR2_max','EVI_min','EVI_max'];
var fire1=  fire.filter(ee.Filter.notNull(indices));
var fire2= fire1.select('Class_R','Dnbr','EVI_max','NBR2_max','imageId')
var bands= ['NBR2_max', 'EVI_max'];
var classProperty= 'Class_R';
var random= fire2.randomColumn('random');
var split_train= 0.8;
var split_test= 0.2;
var trainPartition= random.filter(ee.Filter.lt('random',split_train));
var testPartition= random.filter(ee.Filter.gte('random',split_test));

var RFtrained_c= ee.Classifier.smileRandomForest({
 numberOfTrees: 100,
 minLeafPopulation:2,
 seed: 123,
 })
.setOutputMode('CLASSIFICATION')
.train({
    features: trainPartition,
    inputProperties: bands,
    classProperty: classProperty,
    })
print(RFtrained_c);
var dict =  RFtrained_c.explain();
var test = testPartition.classify(RFtrained_c)
print(test.first())


texts: Spatial joining/merging two spatial points data frames in R via buffer?

texts:
I'm trying to spatial join/merge two spatial points data frames in R. One data frame is lacking IDs, which is the information that needs to be joined to the other data frame. To be specific, any point from data frame 1 that falls within a 1500 m buffer of a point in data frame 2 should be joined with the respective ID from data frame 2. Here's some code to illustrate my problem:

library(mapview)
library(sp)
library(rgdal)
library(rgeos)

lat1 <- c(4417391,4517826, 4435680,4509372,4449390)
long1 <- c(5557780,5358439,5328731,5323168,5519670)
attrib1 <- seq(6,10,1)
attrib2 <- seq(11,15,1)
df1 <- data.frame(lat1,long1,attrib1,attrib2)

lat2 <- c(4394822,4488830,4417257,4517995,4435679)
long2 <- c(5293795,5418630,5557927,5358272,5328084)
ID <- seq(1,5,1)
df2 <- data.frame(lat2,long2,ID)

coordinates(df1) <- ~lat1+long1
proj4string(df1) <- CRS("+init=EPSG:31468")

coordinates(df2) <- ~lat2+long2
proj4string(df2) <- CRS("+init=EPSG:31468")

buff <- gBuffer(df2,width=1500)
df1_buff <- df1[buff,]

mapview(buff) + mapview(df1,color="red") + mapview(df1_buff,color="green")

I have now selected the points from df1 that lie within the buffers of df2. However, the buffer obviously does not contain IDs and I'm at a loss about how to merge this infromation. I have looked into over and point.in.poly, but with no success. This is the result I'm looking for:

attrib1 attrib2 lat1     long1    ID
6       11      4417391  5557780  1
7       12      4517826  5358439  2
8       13      4435680  5328731  3


texts: Can't add elements from an attribute value to a set

texts:
I have duplicate elements in one field (11223334) and would like them to bu unique in another (1234). I made a new field and am trying to use update cursor:

>>> with arcpy.da.UpdateCursor("SD_Pattern_CLSCOM",["Day_Pattern","Day_Pattern2"]) as cursor:
...     for row in cursor:
...         a = set()
...         for i in row[0]:
...             a.add(i)
...         row[1] = a
...         cursor.updateRow(row)

...
And I am getting this error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 7, in <module>
TypeError: value #1 - unsupported type:


texts: Looping with a range of integers for 'numeric' inputs in QGIS Python Console

texts:
I have been trying to run r.sun.incidout in a loop of integers (for the 'Day' option, type=integer) in QGIS Python Console. Since the 'day' parameter takes an integer number for its value, I get the error as:
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS 3.10\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 6, in <module>
  File "C:/PROGRA~1/QGIS 3.10/apps/qgis-ltr/./python/plugins\processing\tools\general.py", line 106, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/PROGRA~1/QGIS 3.10/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 137, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Unable to execute algorithm
Incorrect parameter value for day

I am using the following code. I know the error lies in dys inside the paramms in my code. How can I parse the variable as a numeric in the loop?
import processing
from processing.core.Processing import Processing

Processing.initialize()

import os
import numpy as np

dayss = np.arange(1, 367, 1)

outdir = r'C:\***\swiss_rad_tiffs'

for dys in dayss:

    day_str = str(dys).zfill(3)

    print('Doing Day %s' %day_str)

    outfile = os.path.join(outdir, 'rad_%s.tif'%day_str)

    paramms = {'elevation': r'C:\***\dem_wg32_500m.tif', 'aspect':r'C:\***\asp_wg32_500m.tif', 'aspect_value': 270, 'slope': r'C:\***\slp_wg32_500m.tif', 'slope_value':0.0, 'linke': None, 'albedo': None, 'albedo_value': 0.2, 'lat': None, 'long': None, 'coeff_bh': None, 'coeff_dh': None, 'horizon_basemap': None, 'horizon_step': 1, 'day': dys, 'step': 1.0, 'declination': None, 'distance_step': 0.5, 'civil_time': None, 'time': 12.00, 'beam_rad': outfile, 'GRASS_REGION_PARAMETER': r'C:\***\dem_wg32_500m.tif', 'GRASS_REGION_CELLSIZE_PARAMETER': 0  }

processing.run('grass7:r.sun.incidout', paramms)


texts: How can I split my spatial layer every 250 records?

texts:
I have a shapefile table with more or less 12000 sectors. I want to use this as a zonal feature layer, however the tool I want to use it with can only handle 255 zones. Therefore I want to divide the shapefile in 1200/155 individual files.

I thought to do it through model builder with the Select (analysis) tool. I created a list of all values (0 - 255 - 510 - etc). I want to run Select with each time using different values of my list in the Expression. I have no clue how to iterate through the values in the expression.

There is probably a better method for doing this.


texts: Need to remove small areas, but only if they do not have adjoining boundary

texts:
I'm generalizing some lat/lon data in grass64 and as a precursor to the generalization, I am removing small areas from the data using 'v.clean'.  The following command is working well:

$ v.clean input=foo output=foo_rmarea_100km tool=rmarea thresh=100000000 type=area --overwrite

However, the one problem I have is that it is removing a small area that has a shared border with the rest of the larger data set.

Is it possible to specify that only areas which DO NOT have a shared border with another area AND which are smaller than the defined threshhold be removed?

Thanks


texts: Is population density considered continuous data and why?

texts: Population counts are typically referred to as discrete or quantitative data. Why is population density a continuous data type when it is typically measured for aggregate areas such as census tracts or districts/neighbourhoods (ie, it can't be measured at any point on a surface like gradient or temperature).


texts: PyQGIS: QgsVectorFileWriter.create() doesn't flush even after deleted

texts:
If QgsVectorFileWriter constructor is directly used, the output is flushed when the reference is deleted.  I'm trying to switch to QgsVectorFileWriter.create() because QGIS warns that the former is being deprecated.  But as far as I tested, QgsVectorFileWriter.create() does not flush the output until QgsApplication.exitQgis() is called.

I'm not familiar with Python, QGIS, and C++.  I see that QgsVectorFileWriterTask does something asynchronous, but I'm so weak to follow the code.

QGIS version: 3.10.3
OS: macOS 10.15.3


texts: GEE Select only one band in my imagecollection

texts:
I want that all the images in my imagecollection will contain only one band ('VH_Filtered').
The problem is that  whenever I ask GEE to display the images in the imagecollection, I can see that the images display the original band and not the VH_Filtered and I don't understand why.

This is my code:

//Filter the collection for the VH product from the descending track

var Sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
        .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .select('VH')
        .filterDate('2019-01-01','2019-11-12')
        .filterBounds(geometry);

//Clip the image collection

var clippedVH= Sentinel1.map(function(im){
  return im.clip(geometry);
});

//apply Filter

var filterSpeckles = function(img) {
  var vh = img.select('VH') ;//select the VH polarization band
  var VH_smoothed = vh.focal_median(35,'circle','meters').rename('VH_Filtered'); //Apply a focal median filter
  return img.addBands(VH_smoothed); // Add filtered VH band to original image
};

var Clipped_Speckles = clippedVH.map(filterSpeckles);
print(Clipped_Speckles);

//check that speckle filter worked
var list_before=clippedVH.toList(clippedVH.size());
var list_after=Clipped_Speckles.toList(Clipped_Speckles.size());

var before=ee.Image(list_before.get(0));
var after=ee.Image(list_after.get(0)).clip(geometry).select('VH_Filtered');

//Map.addLayer(before,{min: -30, max: 1},'Before');
//Map.addLayer(after,{bands: 'VH_Filtered',min: -30, max: 1},'after');

// var clippedVHsize=clippedVH.size();
// print('SAR Size:',clippedVHsize);
// print('SAR images data:',clippedVH);

var listOfImagesSAR =(Clipped_Speckles.toList(Clipped_Speckles.size()));

var NumberOfImagesSAR=listOfImagesSAR.size();

var listOfNumbersSAR = ee.List.sequence(0, NumberOfImagesSAR.subtract(2));
listOfNumbersSAR = listOfNumbersSAR.map(function(n) {
  return ee.Number(n).add(1);
});
print(typeof(listOfNumbersSAR));
listOfNumbersSAR=listOfNumbersSAR.getInfo();

for (var i in listOfNumbersSAR) {
  var image = ee.Image(listOfImagesSAR.get(listOfNumbersSAR[i]));
  var toexport=image.visualize({min: -30, max: 1}).addBands(image);

  //Add the layer- here I can see it didnt work
  Map.addLayer(image,{min: -30, max: 1},'SAR'+ i);
  // Export.image.toDrive({
  // image: toexport.toFloat(),
  // description: i,
  // scale:10,
  // crs:'EPSG:4326',
  // maxPixels:1310361348,
  // region:geometry.geometry().bounds()

// });

}

The area is in Thailand so any polygon there wouold work.

my end goal: to have imagecolelction with only one band.


texts:
I'm trying to customize QGIS 3.4 interface to allow data display only. None of the default QGIS geospatial data editing functionalities should be exposed to the end user. Also, the users should not be allowed to mess with the interface, like changing layer properties, layer groups or layer ordering.

Using the standard "Interface Customization" option I could do most of the work, but I can still see a few unwanted options around:

The Attribute Table window still shows a few buttons allowing users to edit vector layer attributes
The "Layers" panel still shows a few options related to interface change (creating groups, removing layers etc.)
If I right-click a layer, I can still access the context menu

See image below indicating what elements I'm mentioning:

I have gone through the widget customization option in the "Interface Customization", but none of the widgets seem to be related to these particular items.

I found a few posts mentioning PyQGIS for interface customization, but they did not provide much information that I could use in my scenario.

Other posts suggest that we develop a completely customized interface from QGIS libraries, which I would avoid at this stage as it doesn't fit in this project's budget.

Is there a easier way to hide/remove the buttons/options I have highlighted in the attached image?


texts: Loop over attribute table of selected layer

texts:
Here's a stab at my first practical ArcPy script.

I want to grab a layer by name (rather than by layer order or by file path) in the current MXD, then loop over it's attributes (with a view to modifying them). However I don't know the syntax of connecting the desired layer to the attribute table

Here's what I've got so far:

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
needle = "1km_Grid"
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd, "", df):
    if lyr.name == needle:
    # curObj = this layer # pseudo code!

columns = ["MID", "Shape"]

for row in curObj:
    print ("column 1 is: %s" % (row[0]))


texts: Making ArcGIS Python Add-In Extension immediately functional after install?

texts:
I have a Python add-in that includes an extension class.  The only purpose of this extension is for some buttons in the add-in to be automatically enabled/disabled depending on if there are relevant layers in the current MXD's table on contents (by using the extension's itemAdded() and itemDeleted() functions).

Therefore I need this extension to be loaded automatically AND enabled when this add-in is installed.

In the config.xml, I have autoLoad="true".  In the Python class, the __init__ function for the extension class includes:

self.enabled = True

I thought this should be enough for it to work as I wanted.

However, when I install the add-in and re-launch ArcMap, the extensions dialog shows the extension but NOT checked.

Now, I have found that if I run something like:

myExtension.enabled = True

Within the __init__ method of some other toolbar class (eg, one of the buttons, NOT the extension class), then the extension IS CHECKED in the extensions dialog when ArcMap is re-launched!  However, it doesn't actually work - it's functions are not called.  But, if I then quit ArcMap and restart ArcMap, the extension is both checked in the dialogue AND fully functional (presumably because it remained enabled from the previous run of ArcMap).

So I can get the desired behaviour ONLY if I do ALL of the following ODD things:

enable the extension in the class of some other add-in button (NOT just in the extension class itself)
launch ArcMap
quit ArcMap
launch ArcMap

After that it works fine on subsequent launches of ArcMap.

Is it possible to get a Python add-in extension to automatically load and be enabled by default when the add-in is installed by any other means, and to have it work the first time ArcMap is launched after the add-in is installed?

Ie, is there a correct way to get an extension to be fully functional immediately after an add-in is installed (and ArcMap re-launched once) without any other manual interaction from the user?  Or is this behaviour simply not supported?


texts: QGIS - graphically depict movement patterns

texts:
I want to make a graphical depiction in QGIS of two different datasets showing pedestrian movement patterns in a city center. One data set is observed, and the other is projected.
My data (CSV) consists of the following columns: 1. LinkID, 2. Pedestrians_Observed, 3. Pedestrians_projected
The LinkID column must be "connected" to a line layer I have in QGIS which covers the entire city center.
Ideally the graphical depiction should look something like this, in a way that makes it possible to visualize the differences between projected and the observed :

I am working in QGIS 3.4.13-Madeira.


texts:
I cannot find any solid info on this and its a pretty rudimentry GIS process.

I have 2 shapefiles.
One shapefiles is a single large polygon.
The other is a larger multipolygon.

I need to clip the multipolygon to the extents of the smaller polygon. None of the multipolygons are completely outside of the clipping polygon.

I can do this quite easily in Arc/QGIS/MsPaint/Whatever but this has to be done in python and using OGR.

Code Snippet:

import os, sys, csv
from osgeo import ogr

## Input
inDataSource = driver.Open('OAs_Temp.shp', 0)
inLayer = inDataSource.GetLayer()

## Clip
inClipSource = driver.Open('_convexhull.shp', 0)
inClipLayer = inClipSource.GetLayer()

## Clipped Shapefile... Maybe???
outDataSource = driver.CreateDataSource('FINAL.shp')
outLayer = outDataSource.CreateLayer('FINAL.shp', geom_type=ogr.wkbMultiPolygon)

featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
ogr.Layer.Clip(inLayer, inClipLayer, outLayer)
outFeature.SetGeometry(outLayer)
outLayer.CreateFeature(outFeature)

Issue is with outFeature.SetGeometry(outLayer)

Traceback:
TypeError: in method 'Feature_SetGeometry', argument 2 of type 'OGRGeometryShadow *'

The problem is that I just can't figure out how OGR is doing things and what it wants. I am finding that cookbook thing pretty lacklustre and incomplete and the documentation is not very clear on how any of this works. It is some of the most convoluted open source libraries I have ever worked with and I regret taking this project on.

This is my last option. As always.


texts: Overlaying point data with NDVI layer with raster output (not graph)

texts:
I have point data, which are animal burrows, and a NDVI raster layer that shows plant productivity. I want to overlay these layers and determine if vegetation productivity is higher close to burrows (or in a certain radius around burrows). The output must be a raster and not a correlation graph, something like this image (Bonachela et al., 2015):
showing how productivity changes close to burrows and the further you move away from burrows. Can this be done with QGIS, R or another program?


texts:
I'm wondering if it is possible to include some las files in a LAScatalog only for the purpose of providing a buffer to other las files.

For example, I may only want to normalize 2 las files and use buffer zones of neighbouring files to avoid edge artefacts, but never fully process the neighbouring files.

I don't think lasclip() will do what I want since lasclip() essentially creates a single las file from a region of interest.