texts: PostGIS - Partition over distinct geometries ignoring directionality

texts:
I have a table containing geometries, some of which are spatially equal (ST_Equals returns true), but are not strictly equal due to different directionality (= operator returns false, assuming the new behavior of the = operator since PostGIS 2.4):

 gid |        geom
-----+---------------------
   1 | LINESTRING(0 0,1 1)
   2 | LINESTRING(1 1,0 0)
   3 | LINESTRING(2 2,3 3)
   4 | LINESTRING(0 0,1 1)
   5 | LINESTRING(4 4,5 5)
   6 | LINESTRING(3 3,2 2)

Note in the table above, there are 3 distinct geometry groups:

Geometries with gid 1, 2, and 4 are spatially equal (despite different directionality).
Geometries with gid 3 and 6 are spatially equal (despite different directionality).
Geometry with gid 5 is not spatially equal with any other geometry in the table.

I would like to partition over this table and assign a geometry "group ID" to each row. The group ID should be the same for all rows where the geometry is spatially equal (ignoring directionality).

The desired result should look like this:

 gid |        geom         | grpnum
-----+---------------------+--------
   1 | LINESTRING(0 0,1 1) |      1
   2 | LINESTRING(1 1,0 0) |      1
   3 | LINESTRING(2 2,3 3) |      2
   4 | LINESTRING(0 0,1 1) |      1
   5 | LINESTRING(4 4,5 5) |      3
   6 | LINESTRING(3 3,2 2) |      2

If directionality could be ignored, I could do something like this (the geoms CTE is just for sample data):

WITH
  geoms AS
  (
    SELECT 1 AS gid, ST_GeomFromText('LINESTRING(0 0, 1 1)') AS geom
    UNION ALL
    SELECT 2 AS gid, ST_GeomFromText('LINESTRING(1 1, 0 0)') AS geom
    UNION ALL
    SELECT 3 AS gid, ST_GeomFromText('LINESTRING(2 2, 3 3)') AS geom
    UNION ALL
    SELECT 4 AS gid, ST_GeomFromText('LINESTRING(0 0, 1 1)') AS geom
    UNION ALL
    SELECT 5 AS gid, ST_GeomFromText('LINESTRING(4 4, 5 5)') AS geom
    UNION ALL
    SELECT 6 AS gid, ST_GeomFromText('LINESTRING(3 3, 2 2)') AS geom
  )
SELECT
  geoms.*,
  dense_rank() OVER (ORDER BY geom) AS grpnum
FROM geoms
ORDER BY gid;

The result of which is:

 gid |      st_astext      |   grpnum
-----+---------------------+------------
   1 | LINESTRING(0 0,1 1) |          1
   2 | LINESTRING(1 1,0 0) |          2
   3 | LINESTRING(2 2,3 3) |          3
   4 | LINESTRING(0 0,1 1) |          1
   5 | LINESTRING(4 4,5 5) |          5
   6 | LINESTRING(3 3,2 2) |          4

However, that assigns a different group to 1 and 2, as well as to 3 and 6, which I would like to treat as equal. I need to use ST_Equals here somehow, which has the behavior I want (it ignores directionality), but I'm not sure how to structure the query since ST_Equals requires that you pass it in a pair of geometries.


texts: osm2pgsql. Loading data into tables with non-standard columns

texts:
I have an * .osm file with a lot of data. To load the data into postgis I use osm2pgsql.
Tell me, please, how to load in the table "my_table" data that have amenity = parking. Coordinates of polygons must be written in the "polygon" field


texts:
I have a point table and the buffer geometry is created around the point.

I have defined the if the point is within the buffer_distance then its value is =1
and if it is outside the buffer its value is =0, and when it is on buffer_polygon then =2.

Now when ever i change the location of point then it should check the following conditions and the flag value should change according to the location of point new position.

I am changing the position of the point in QGIS.

My table structure is like this:-

GID(PK) Flag(INT)   Buffer_Distance THE_Geom    Buffer_GEOM
1   1   100 point(X,Y Location) Buffer Created around the Point
2   2   200 point(X,Y Location) Buffer Created around the Point
3   2   200 point(X,Y Location) Buffer Created around the Point
4   3   300 point(X,Y Location) Buffer Created around the Point
5   3   300 point(X,Y Location) Buffer Created around the Point
6   1   100 point(X,Y Location) Buffer Created around the Point


texts: Why does my buffer take up the entire screen?

texts: I am using QGIS. I am aware there is a similar question here but the solution (changing the layer CRS) is not working for me. I cannot change my CRS as I am using information for the entire Asia-Pacific region. Is there another option to create a 1000 km buffer zone from all volcanoes in the region? It is currently projected as geographic coordinates, the default WGS:84.


texts: "ee.Algorithms.If" gives different value by NDVI

texts:
I want to use ee.Algorithms.If. It gives different value by NDVI. I write code like this：

fsl=ee.Algorithms.If(ndvi1.gt(0.2),
                     0.92,
                     ee.Algorithms.If(ndvi1.lt(0.5),
                     0.99,
                     0.8))；

I want to use this Polynomial. If NDVI is less than 0.2, the grid value is 0.9, if NDVI is greater than 0.5, the grid value is 0.92, if NDVI is between 0.2 and 0.5, the grid value is 0.91.

I get the error Cannot add an object of type <ComputedObject> to the map.


texts:
I have a very good idea of what I would like to do but I really do not know if it is possible with QGIS.
I have points that each one corresponds to a group size and I would like to use GIS to sort this data. In fact, in order to conserve independence between my group sizes I would like to keep only one group size per day.
Thus, is QGIS capable of analysing all my points and telling me how many duplicates I have within 100 meters radius for the same day? Then, if the software finds more than one point per day, can I run an application that will only keep one group size (randomly). I do not know if I am being very clear but I hope you understand my problem.

Date    Day Month   Year    Week no Season  Time    EASTING NORTHING    DATUM   Groupsize
14/03/2006  14  March       2006    11  WINTER  09:20   599707  5733317 NAD83   1
02/04/2002  2   April       2002    14  SPRING  14:34   602500  5732500 NAD83   500
05/02/2002  5   February    2002    6   WINTER  09:00   600359  5734295 NAD83   250
06/02/2002  6   February    2002    6   WINTER  08:30   599000  5733600 NAD83   600
07/02/2002  7   February    2002    6   WINTER  18:00   599200  5733400 NAD83   650
11/02/2002  11  February    2002    7   WINTER  14:30   599850  5734150 NAD83   565
12/02/2002  12  February    2002    7   WINTER  17:10   599661  5733382 NAD83   29
13/02/2002  13  February    2002    7   WINTER  09:00   599800  5734200 NAD83   650
16/02/2002  16  February    2002    7   WINTER  10:30   599010  5733850 NAD83   100
16/02/2002  16  February    2002    7   WINTER  10:30   601034  5734740 NAD83   150
16/02/2002  16  February    2002    7   WINTER  10:30   600463  5734418 NAD83   75
16/02/2002  16  February    2002    7   WINTER  10:30   600628  5734783 NAD83   150
16/02/2002  16  February    2002    7   WINTER  10:30   600147  5734358 NAD83   100
20/02/2002  20  February    2002    8   WINTER  12:00   602500  5734200 NAD83   325
26/02/2002  26  February    2002    9   WINTER  18:05   599661  5733382 NAD83   155

Thank you,


texts: Calculate sum of attribute field using Leaflet and Esri Leaflet plugin

texts:
I am adding Esri feature layer to Leaflet map from Esri REST API.
I can access the feature properties, like this

feature.properties.CasesField

How can I iterate through this CasesField and get:

All the values
Sum the values

function initMap() {
  map = L.map("map", {
    zoomControl: true,

    zoomSnap: 0.25
  }).setView([33.1625662, -96.6390040], 5); //home position
  ;

  L.esri.basemapLayer("Topographic").addTo(map);

  casese = L.esri.featureLayer({
    url: casesURL,
    onEachFeature: function (feature, layer) {
      var props = feature.properties.CasesField;

      var sum = 0;
      for (x in props) {
        sum += props;
      }

      // use the value to do something...

      console.log(sum);

    }

  }).addTo(map);

}


texts: OpenLayers React export map and vector layers

texts:
For my project I am trying to export the map I have created with all it's Vector layers. I can't seem to get the example I found on the site of OpenLayers to work in React. The export doesn't need to produce a PDF. It would be nice but if I can export a image file, that would be nice as well. It's just that I can't seem to get the code to work in React.

  exportMap() {
    this.map.once("rendercomplete", () => {
      var viewResolution = this.map.getView().getResolution();
      var mapCanvas = document.createElement("canvas");
      var size = this.map.getSize();
      mapCanvas.width = size[0];
      mapCanvas.height = size[1];
      var mapContext = mapCanvas.getContext("2d");

      this.map.getLayers().forEach((canvas) => {
        if (canvas !== undefined) {
          if (canvas.width > 0) {
            var opacity = canvas.parentNode.style.opacity;
            mapContext.globalAlpha = opacity === "" ? 1 : Number(opacity);
            var transform = canvas.style.transform;
            // Get the transform parameters from the style's transform matrix
            var matrix = transform
              .match(/^matrix\(([^\(]*)\)$/)[1]
              .split(",")
              .map(Number);
            // Apply the transform to the export map context
            CanvasRenderingContext2D.prototype.setTransform.apply(
              mapContext,
              matrix
            );
            mapContext.drawImage(canvas, 0, 0);
          }
        }
      });

      var pdf = new jsPDF("landscape", undefined, "a4");
      pdf.addImage(mapCanvas.toDataURL("image/jpeg"), "JPEG", 0, 0, 297, 210);
      pdf.save("map.pdf");
      // Reset original map size
      this.map.setSize(size);
      this.map.getView().setResolution(viewResolution);
      document.body.style.cursor = "auto";
    });
    this.map.renderSync();
  }


texts:
I am trying to find out if a certain function is available to link an excel sheet with a kmz.

What i need to happen is have a Placemark that pulls information off an excel sheet row and display's the data in the description part of the Placemark.

I know that there is a macro in excel to export KMZ file but I want the KMZ itself to perform the task of collecting the information and displaying it.

I would guess that it would be an extension of the “ad link command,” but I can’t seem to get my head around it. This would be cut down on data duplication if there was a way of doing it.

If you know of a way to do this please comment.

Thanks

Shane


texts: Search images in USGS earth explorer using Python API

texts:
Is it possible to search for images in USGS Earth Explorer using Python (API) script?

I tried the example here:

from usgs import api

# Set the EarthExplorer catalog
node = 'EE'# this indicates earth explorer website

# Set the Hyperion and Landsat 8 dataset
hyperion_dataset = 'EO1_HYP_PUB'
landsat8_dataset = 'LANDSAT_8'

# Set the scene ids
hyperion_scene_id = 'EO1H1820422014302110K2_SG1_01'
landsat8_scene_id = 'LC80290462015135LGN00'

# Submit requests to USGS servers
api.metadata(hyperion_dataset, node, [hyperion_scene_id])
api.metadata(landsat8_dataset, node, [landsat8_scene_id])

But I got:

USGSError                                 Traceback (most recent call last)
<ipython-input-4-dab0c21371c4> in <module>()
     11
     12 # Submit requests to USGS servers
---> 13 api.metadata(hyperion_dataset, node, [hyperion_scene_id])
     14 api.metadata(landsat8_dataset, node, [landsat8_scene_id])
     15

C:\Users\Ran\Anaconda2\lib\site-packages\usgs\api.pyc in metadata(dataset, node, entityids, extended, api_key)
    235     response = r.json()
    236
--> 237     _check_for_usgs_error(response)
    238
    239     if extended:

C:\Users\Ran\Anaconda2\lib\site-packages\usgs\api.pyc in _check_for_usgs_error(data)
     33     error = data['error']
     34
---> 35     raise USGSError('%s: %s' % (error_code, error))
     36
     37

USGSError: AUTH_UNAUTHORIZED: Forbidden - Valid API Key required for access to 'metadata' method.

I don't understand where my password and user name should go?

I also tried this according to this:

import landsatxplore.api

# Initialize a new API instance and get an access key
api = landsatxplore.api.API('user', 'password')

# Request
scenes = api.search(
    dataset='LANDSAT_8',
    latitude=19.53,
    longitude=-1.53,
    start_date='2015-01-01',
    end_date='2016-01-01',
    max_cloud_cover=10)

print('{} scenes found.'.format(len(scenes)))

for scene in scenes:
    print(scene['acquisitionDate'])

api.logout()

But I got:

---------------------------------------------------------------------------
EarthExplorerError                        Traceback (most recent call last)
<ipython-input-6-a48259d60bd9> in <module>()
     11     start_date='2015-01-01',
     12     end_date='2016-01-01',
---> 13     max_cloud_cover=10)
     14
     15 print('{} scenes found.'.format(len(scenes)))

C:\Users\Ran\Anaconda2\lib\site-packages\landsatxplore\api.pyc in search(self, dataset, latitude, longitude, bbox, max_cloud_cover, start_date, end_date, months, max_results)
    111             params.update(months=months)
    112
--> 113         response = self.request('search', **params)
    114         return response['results']
    115

C:\Users\Ran\Anaconda2\lib\site-packages\landsatxplore\api.pyc in request(self, request_code, **kwargs)
     39         response = requests.get(url, params=params).json()
     40         if response['error']:
---> 41             raise EarthExplorerError('EE: {}'.format(response['error']))
     42         else:
     43             return response['data']

EarthExplorerError: EE: Dataset is currently offline.

Do they have an official API for Python?
I would like to search (and download) images from the Earth Explorer database (not only Landsats) using Python, is it possible?


texts: GDAL for Python: how can I get XYZ surface from a Geotiff?

texts:
I want to get an XYZ-like surface from a Geotiff with GDAL in Python.

I do the following:

1) Get some geotiff from here: http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

2) Write the script:

from osgeo import ogr, osr, gdal

file_path = os.path.join(os.path.dirname("/home/felix/dump/srtm_62_05.tif")
gdal_data = gdal.Open(file_path, gdal.GA_ReadOnly)
raster_band = gdal_data.GetRasterBand(1)
raster = raster_band.ReadAsArray()

Okay, it works and even gives me some data in [latitude, longitude --> some_strange_elevation_units] format.

I want to convert this raster to some X[i,j], Y[i,j], Z[i, j] surface. And I need the units to be scaled in meters. How can it be done?

P.S. I use GDAL 2.1.3


texts: Changing double precision raster file using ArcGIS Desktop?

texts:
I am using ArcMap 10.2.

I’m trying to create an attribute table for an NDVI tif file. The pixel type of my file is double precision, so I tried first multiplying the whole thing by 10^n, however I have 16 decimal places. For some reason when I use the INT tool afterwards it just gives me one category (which is 0). I can try multiplying by 10^9, and then I have no problem using the INT tool, but I feel like I’m losing resolution by using only 9 places.

Is there any other method to change from double precision into signed integer?


texts: How to compile QGIS from source with custom pyspatialite 4.0.0?

texts: I have turned WITH_PYSPATIALITE ON, but I have a custom build of SpatiaLite (4.0.0), not the packaged version, so I set WITH_INTERNAL_SPATIALITE to OFF. How do I tell QGIS to configure PYSPATIALITE for 4.0.0? Or should I turn WITH_PYSPATIALITE OFF and install it myself?


texts: Transforming Datum in Lat/Long to Projection in Feet

texts:
I am having issues projecting my csv file I have imported in lat/long to our state government-standard projection (EPSG 102121 NAD 1983 Michigan GeoRef (feet)), which is based on Oblique Mercator.

When I pull in the csv, I load the lat/long points, I select 'EPSG 102121- NAD 1983 Michigan GeoRef Feet US' for the "Geometry CRS". When converted to points, the layer misaligns with others in the same projection. My points are still in lat/long, but the other layers have been successfully reprojected to feet.

How do I transform the lat/long points to be reprojected in feet in QGIS? I understand it is still a datum at this point, and that another step is needed, but I guess I assumed since the NAD 1983 MI GeoRef is based on Oblique Mercator, the projection was built in or something.

Is there a step I am missing that requires the geometries of the lat/long to be recalculated in feet?  I've only had this issue working with csv files, and it's frustrating.


texts: How to run a custom Geoprocessing tool using no parameters?

texts:
I created a custom geoprocessing tool called "CreateCoverLayer". It has no input or output parameters.  I have the following snippet of code, but get an error message due to the fact that the Execute method is looking for an array of parameters as the second input argument. Can I use this tool with no parameters, or is that not possible?  Thanks!

IGeoProcessor2 gp = new GeoProcessorClass();
gp.AddToolbox(@"C:\TestToolbox.tbx");
gp.Execute("CreateCoverLayer", null, null);


texts:
One of the main reasons I use groups in Arc and QGIS is to be able to group similar data and/or associated layers which I will want to display or turn on/off together.

In Arc if I have a group with say 5 layers, and I have 2 of those turned on, I can simply turn the entire group off and both those layers will turn off. Importantly though if I then turn the group back on, those two layers turn back on. The 'on/Off' status of each layer within the group is not changed by changing the whole groups On/Off status.

In QGIS when I click to turn a group off QGIS first turns all the layers in that group on, and then on clicking again will turn all the layers off, as per the image below. This is frustrating as I just want to turn the group off to collectively turn those two layers off, and most commonly then back on again after I have published a map or achieved what ever I needed to without those layer displayed.

My question is, is there a way in QGIS 2.18 to turn off a group  while retaining the On/Off status of the layers in the group such as Arc does?


texts: How can I export images with their metadata from Google Earth Engine to Google Cloud Storage?

texts:
I want to transfer some images from Earth Engine to a Google Cloud Storage bucket. I know that I can use the following code to accomplish this task:

Export.image.toCloudStorage({
            image: image,
            bucket: bucket_name,
            fileNamePrefix: file_name_prefix,
            fileFormat: 'GeoTIFF'
        })

This stores my image in the form of several chunks inside my bucket. However, I cannot see the image metadata in any of these chunks. Can someone advise me as to how to save the metadata along with the images?


texts: Select / Drag or by rectangle CMV Application

texts: Anyone able to suggest to this noob if there is a way to change the selection option in my CMV application from individual polygons (one at a time) to an option of selecting multiple polygons at once and having attributes display in the attribute table widget, in this case "parcels", by a shape of rectangle. I've notice documentation online also references it as "Select / Drag". I can't seem to find any documentation on how i would write or change the existing cmv script. I've looked at other JavaScript api maps that have this method of selection and cannot figure out how to implement this into mine.


texts: Changing ArcSDE connection credentials for map document?

texts:
Is it possible to change the ArcSDE credentials used by an MXD, so that users with the same connection can open the MXD and use their own username and password to access that data?

With our current setup, if we share an MXD, the person attempting to open it will be prompted for username and password, however since it's not the same as the person who created it, the connection is unable to be made and the account of the opening person becomes locked (if attempted 3 times).

Using ArcMap 10.3, ArcSDE is on Oracle database


texts: How can I convert a .tif DEM to XYZ format using QGIS?


texts: GDAL contour line generation at DEM edges

texts:
In a Python application I use gdal.ContourGenerate() to generate contour lines (GDAL version 1.11.4 and 2.2.1) from DTED tiles. DTED tiles overlap at the top and the right with neighboring tiles.

Comparing contour lines build from several individual DTED tiles to these, build from a VRT show discrepancies at the tile edges. That is because contour lines are extended perpendicular to the tile borders. That is the reason, why contour lines from neighboring tiles do not meet at the same position.

According the documentation this shouldn't be the case:

  This means that contour lines' nodes will not actually be on pixel
  edges, but rather along vertical and horizontal lines connecting the
  pixel centers.

Is there a way to restrict contour line generation to the area between raster cell centers?


texts: How can I optimize this point-linking query?

texts:
I have the following query, and am wondering how it could be optimised to run faster:

WITH
gps_nodes_1 AS ( SELECT "NODE_ID", "geog" FROM gps_nodes WHERE ST_Distance("geog",ST_GeogFromText('POINT( -1.54790 53.79424 )')) < 40000 ),
gps_nodes_2 AS ( SELECT "NODE_ID", "geog" FROM gps_nodes WHERE ST_Distance("geog",ST_GeogFromText('POINT( -1.54790 53.79424 )')) < 40000 )
SELECT
    row_number() over() as LINK_ID,
    gps_nodes_1."NODE_ID" as source,
    gps_nodes_2."NODE_ID" as target,
    ST_Distance(gps_nodes_1."geog",gps_nodes_2."geog")::numeric(7,1) as length
INTO
    gps_node_links_neighbouring
FROM
    gps_nodes_1
INNER JOIN
    gps_nodes_2
ON ( ST_Distance(gps_nodes_1."geog",gps_nodes_2."geog")::numeric(7,1) BETWEEN 0 AND 20 )

A description of what it does:

Table "gps_nodes" consists of ~2 million GPS positions. It has two columns, "NODE_ID" (an auto incrementing id) and "geog" which is a PostGIS geography point.

For each GPS node, I want to create a link to all other GPS nodes within 20 metres, and store that link (with start id, end id, and length) in a new table.

To start with, rather than try doing this for all GPS nodes, I want to limit it to those within 40 kilometres of a set point - which is around 90,000 records in the gps_nodes table.

So I started by creating two CTEs of the GPS points within 40km of required position.
Then creating an inner join between the two CTEs, where the join condition is that the points are between 0 and 20 metres apart (need to be greater than zero distance, otherwise you always get a point linked to itself).

I freely admit that I'm an SQL novice and there may be far better ways of doing this!

Can anyone suggest a more optimised version of this query ?

Many thanks.


texts: How to find a list of every REAL residential address in an area?

texts:
I'm trying to get demographic data in the US at a level more specific than census block groups, so am going to try out using the zillow api and other similar apis to get home price and rental data for individual homes.

So far, I can't find a way to get the zillow api to return a list of every datapoint it has in its database (within a given area) in the form of a very large array, so my strategy currently is to make specific address requests to the API.

So what I need now is a list of every home address in a given area. My first guess is low hanging fruit, just get them from a phonebook. That's a good start but I don't know an open api or data dump for a phonebook currently. Will post sources as I learn them. What would be GREAT would be something like Open Street Map for homes. Not sure if anything like that exists though


texts:
I am using ArcGIS 10.5 and I am unable to set the destination format when I try to export an image using the Export Raster option on ArcMap.

The steps:

But the destination format is locked on "File geodatabase".

On later versions I was able to choose the format in this field.

What can I do in order to activate this field and choose a different format (TIFF, JPEG, etc.)?


texts: Legend and Attribute Table change throughout map composer

texts:
I'm using QGIS 2.18.3 (Las Palmas)
For a project I have about 15 maps with the same extent, where each map displays a different attribute. I have followed Automating Atlas to create maps for all Attributes and it works perfectly.

Two things are not working as they should:
- I would like to add a legend to each map (I wouldnt mind doing so manually), but every time I go back and forth though the maps, the last legend created is set for all the maps. Since all maps display something different, I really need to make sure that each map has its own legend.

I'm adding the attribute table with the relevant parameters I want to show, but - as with the Legend - all 15 maps show the same attribute Table.

Do you have any suggestions how I can make sure that each map in the Atlas have its own Legend and Attribute Table?


texts: Data Store error when registering a database in ArcGIS

texts:
I am using ArcMap 10.3.1 and have downloaded the PostgreSQL 9.3.5 files from MyESRI.com. PostgreSQL is installed locally for testing purposes. I have everything working smoothly and am attempting to publish the data for use on a tablet in the Esri Collector app. In order to do this I have to register the database with ArcGIS Server.

I edited the pg_hba.conf to allow remote connections, since ArcServer is installed on another machine. When I attempt to register the database I get he error in the following screenshot. I have attempted an exhaustive Google sreach for a resolution with no success.


texts: Backup GIS data on postgis database from pgadmin3

texts:
I am using a postgis database with qgis and managing it with pgadmin3

I have already investigated the backup/restore options within pgadmin 3 and done one test run on one database

I am wondering if this .backup file is the only/best way to do backups?

Ideally I would want to back up several postgis databases, daily and probably send the file to a remote server for offsite storage or at least email it.

Does anyone have any suggestions on backing up postgis databases beyond my current experience?


texts: How to use current layer name in FME PythonCaller?

texts:
In a PythonCaller (FME 2014 SP4 in Data Interop in ArcGIS 10.3), I'd like to do some conditional processing based on the layer that I'm currently processing. I'm feeding in a text file of layers dynamically at runtime, so the layer(s) being processed can vary. I'd like to do something like the psuedocode below:

import fmeobjects
def processFeature(feature):
    # Get which layer is currently being processed
    layer_name = somehow get current layername via fmeobjects??
    # do some logical processing of fields based on what layer we are processing
    pass

Not sure if there is something in FME_MacroValues that exposes this (looked and couldn't find much) or if I need to expose it somehow from my source and then pipe it into my PythonCaller perhaps.


texts: Are the downloaded map legend units different than those in the online data explorer?

texts:
I have just recently downloaded the 0.05m depth CEC map from SoilGrids and noticed that the map units do not match those in the WMS version, or those in the file explorer online (even if divided by 10 per the mmol to cmol conversion). The downloaded data has a legend that goes from 0-392 and the WMS version goes from 0-852 (viewed in QGIS).
Are the values for band 1 actually in cmol(c)/kg or is there something I have to do to the data after loading into QGIS.
I am a bit new to this.


texts:
I am trying to understand why the output includes None in tracepoints. What do these mean?

  [None, None, None, {u'waypoint_index': 0, u'name': u'West Washington Street',

And then the algorithm would return a coordinate sometimes without a street name. Are these erroneous output?

  {u'waypoint_index': 26, u'name': u'', u'hint': u'eWK6g____382AAAAjwAAAAAAAAB-AAAADQAAACQAAAAAAAAAHwAAAJo5AAB55FH5-qL-AY3kUfkRo_4BAAAvE3jrEvY=', u'matchings_index': 0, u'location': [-112.073607, 33.465082], u'alternatives_count': 0}


texts: Translating multiband netCDF to GTiff using gdalwarp and Python

texts:
How do I translate (extract) multiband netCDF to individual raster geotiff?

I am new to scripting, so I am hoping you could help me with how I should write my gdal translate script starting from setting my working directory.

Instructions that come with the data:

The maps can be transformed to the GTiff format using GDAL as follows, replacing Year, Band_Name and Output_Filename with the parameters of interest:

gdalwarp -of Gtiff -co COMPRESS=LZW -co TILED=YES -ot Byte -te -180.0000000 -90.0000000 180.0000000 90.0000000 -tr 0.002777777777778 0.002777777777778 -t_srs EPSG:4326 NETCDF:xxxYEAR.nc:BandName OutputFilename.tif

So far:

import os, sys # the standard imports

BaseFolder = r'C:\Users\Dlan\Desktop\netCDF_TIFF' # change this to match your data

for FullPath, dirs, files in os.walk(BaseFolder):
for ThisFile in files:                   # iterate the files
    fN,fE = os.path.splitext(ThisFile)   # separate file name and extension
    if fE.lower() == 'LCCS-2018.nc':
        # this is a netCDF file, though other extensions might exist
        # you will need to provide these details or ascertain them from
        # the input NetCDF file name.

        Year = '2018'
        Band = 'lccs_class'
        OutputName = 'LCCS-2018'

        GDAL_Command = 'gdalwarp -of Gtiff -co COMPRESS=LZW -co TILED=YES -ot Byte -te '\
                                     '-180.0000000 -90.0000000 180.0000000 90.0000000 -tr 0.002777777777778 0.002777777777778 '\
                                     '-t_srs EPSG:4326 NETCDF:{}YEAR.nc:{} {}.tif'.format( Year, Band, os.path.join(FullPath,OutputName))
        os.system(GDAL_Command) # run the command

There are some other instructions:

The maps can be transformed to the GTiff format using GDAL as follows, replacing Year, Band_Name and Output_Filename with the parameters of interest:

with Year the year : 2016, 2017 or 2018; and BandName one of the following: lccs_class; processed_flag; current_pixel_state; observation_count; change_count.

Format:
    gdalwarp -of Gtiff -co COMPRESS=LZW -co TILED=YES -ot Byte -te -180.0000000 -90.0000000 180.0000000 90.0000000 -tr 0.002777777777778 0.002777777777778 -t_srs EPSG:4326 NETCDF:XXX-YEAR.nc:BandName OutputFilename.tif

Example code:
gdalwarp -of Gtiff -co COMPRESS=LZW -co TILED=YES -ot Byte -te -180.0000000 -90.0000000 180.0000000 90.0000000 -tr 0.002777777777778 0.002777777777778 -t_srs EPSG:4326 NETCDF:LCCS-2018.nc:lccs_class LCCS-2018.tif


texts: Error: Intersect 2 polygons in Python. Record's geometry not match: 'GeometryCollection' != 'Polygon'

texts:
I'd tried replicating intersecting two shapefiles from Python or command line by anna-magdalena which requires intersecting two shapes. I used @gene answer for that:

import os
from shapely.geometry import mapping, shape
import fiona

os.chdir('/media/dir')
# Create schema of values that new shape will contain
schema =  {'geometry': 'Polygon','properties': {'area': 'float:13.3', 'id1': 'str', 'id3': 'str'}}

with fiona.open('13fiona.shp', 'w',driver='ESRI Shapefile', schema=schema) as output:
    for c1 in fiona.open('1.shp'):
        for c2 in fiona.open('3.shp'):
           if shape(c1['geometry']c1).intersects(shape(c2['geometry'])):
              area = shape(c1['geometry']).intersection(shape(c2['geometry'])).area
              prop = {'area': area, 'id1' : c1['id'],'id3': c2['id']}
              geom = mapping(shape(c1['geometry']).intersection(shape(c2['geometry'])))
              output.write({'geometry':geom,'properties': prop})

And I got this error:

ValueError: Record's geometry type does not match collection schema's geometry type: 'GeometryCollection' != 'Polygon'

After check the structure of objects c1, c2 and geom I got this:

Structure of c1 polygon:

{'geometry': {'type': 'Polygon', 'coordinates': [[(-74.71213794799996, 1.2098971930000744), (-74.71 ...

and c2 polygon:

{'geometry': {'type': 'Polygon', 'coordinates': [[(-74.69988807399994, 1.2222...

I found that intersection polygon assigned to geom objetc have this structure:

{'type': 'GeometryCollection', 'geometries': [{'type': 'Point', 'coordinates': (-74.71039521399996, 1.2098918560000698)}, {'type': 'Point', 'coordinates': (-74.70929715499994, 1.2102729000000636)}, ...

The c1 and c2 shapes cand be found here
I'm using shapely 1.6b2 version

Any advice?


texts: Desktop ArcGIS Training Sample Manager cannot merge training samples

texts:
I have done the following procedure which has resulted in a troubling behavior:

I have a Shapefile Feature Class with multiple different areas under specific distinct description field. I need to make a supervised classification and construct a training sample from e.g, 10 samples to 1 class.

By querying I have made as many Feature Classes as there are different descriptions. Then in order not to make 10 polygons for each one, I have extracted them as Shapefiles (X_description.shp) which are fine indeed on the map.

Then I open the Training Sample Manager which allows me to Load Training Samples and I open a Shapefile, delete all but 10 sample ROIs and try to merge them in one class.

And here is the issue : ArcGIS does not respond to "Merge Training Samples" at all.

However when I draw polygons it is just fine.

Does anyone have an idea why this happens and how to overcome it in order NOT to make hundreds of polygons manually?

An explanatory screenshot of the procedure and the issue:


texts:
I downloaded GDAL on my Mac (running Mojave 10.14.1) through Homebrew and I'm having trouble getting OGR to work with Postgresql. When I try the following to create a geojson from gpx files I have, I get the following error:
    ERROR 1: Unable to find driver 'PostgreSQL'.

I have tried a bunch of commands to install gdal with postgresql including:

$brew install gdal --with-postgresql

and the command is being recognized and thus ignored. Everytime
 gdal 2.3.2, released 2018/09/21 is downloaded.

I've also tried d/ling GDAL Complete from kyngchaos and while it installs, I have no idea where it is on my Mac HD and it's not on the local driver so that Terminal can read it.

Can anyone help?

It seems to be a common problem from what i've researched but there is not good solution to it.


texts: Creating feature collection after importing different shapefiles from assets in Google Earth Engine

texts:
I am importing 2 different shapefiles from assets to plot the time series chart of MODIS

but I am unable to make the feature collection out of these shapefiles
how to do that?

// 1st Shapefile
var Haryana_state = ee.FeatureCollection('users/abhilashaanu92/HaryanaBoundary');

Map.addLayer(Haryana_state);

// 2nd Shapefile
var Punjab_state = ee.FeatureCollection('users/abhilashaanu92/punjab_state_boundary');

Map.addLayer(Punjab_state);

// Combine both shapefile features into a single feature collection.
var both_states = ee.FeatureCollection([Haryana_state, Punjab_state]);

Map.addLayer(both_states);

// Load MODIS vegetation indices data and subset annual images.

var vegIndices = ee.ImageCollection('MODIS/006/MOD13A1')
                     .filter(ee.Filter.date('2019-01-01', '2020-01-01'))
                     .select(['NDVI', 'EVI']);

// Define the chart and print it to the console.
var chart =
    ui.Chart.image
        .seriesByRegion({
          imageCollection: vegIndices,
          band: 'NDVI',
          regions: both_states,
          reducer: ee.Reducer.mean(),
          scale: 500,
          seriesProperty: 'label',
          xProperty: 'system:time_start'
        })
        .setOptions({
          title: 'Average NDVI Value by Date',
          hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
          vAxis: {
            title: 'NDVI (x1e4)',
            titleTextStyle: {italic: false, bold: true}
          },
          lineWidth: 5,
          colors: ['f0af07', '0f8755', '76b349'],
        });
print(chart);

if I put the name of a single shapefile variable in the regions then it is working but when I combine the 2 shapefiles into a feature collection then it is showing the error.
Error generating chart: Collection.geometry: Unable to use a collection in an algorithm that requires a feature or image. This may happen when trying to use a collection of collections where a collection of features is expected; use flatten, or map a function to convert inner collections to features. Use clipToCollection (instead of clip) to clip an image to a collection.

How to rectify this?
After adding .flatten() as in here
var both_states = ee.FeatureCollection([Haryana_state, Punjab_state]).flatten();

both the polygons merged into a single layer and the chart drawn was of the Punjab region only, as the values I cross-checked after exporting the CSV of Chart.

but I want a chart something like this, where the time-series is of both Haryana and Punjab on the same chart.


texts: ESRI JS API - Get Data in Infowwindow when scrolling through multiple selected features

texts:
I would like to be able to grab attribute data from a popupwindow when its associated feature is clicked. When there is only one feature in the window this is no problem by accessing the features by:

map.infowindow.features[0].attributes["desired attribute"].

If the info window contains more than one point, however, this cannot tell be which one is visible in the popup. My goal is to be able to detect when the window is scrolled through, identify what feature is shown, loop through the features attributes for the desired data, and then pass that data to another function.

I desire to do that so a graph can be synced with whatever is in the info window, but I do not want the graph in the info window itself.

Ive noticed that the infowindow object has a selectNext and Select Previous that I dont know if its useful or not.

EDIT**

        map.infoWindow.on("selection-change", function () {
            // get the selected feature/graphic in the popup window.
            var feature = map.infoWindow.getSelectedFeature();
            var featureMeter = feature.attributes["Meter"];
            if (feature) {
                console.log(featureMeter);
                // do stuff with the feature.
            }
        });

        meterLayer.on("click", function (evt) {
            map.infoWindow.setFeatures([evt.graphic]);

        });

The initial click on the map returns three values the last of which is the correct value. The "next" icon will only return one and the correct value.


texts: How to create a turn restriction in ArcObjects

texts: I am having difficulties creating turn restrictions in ArcObjects. I have found a way to Add a a network turn to restrictions through the INetworkForwardStarEx.AddTurnRestriction() interface. But how to I create the INetworkTurn that this method takes.


texts:
I have QGIS with proj 4.8.0 installed. I'd like to add Winkel Tripel projection which seems to be defined as

+proj=wintri

But that definition string does not work. Is this a bug in my proj version, or what's the correct definition string?


texts: Arcgis Javascript Api Identity Manager with proxy

texts:
I'm trying to setup an authentication architecture with a web app using arcgis Javascript api with a configured .NET proxy and a web adaptor.

When i load my app the identity manager shows up for locked resources.

Do i really need the proxy?
Why identity manager shows up on every page refresh when the validity is 60 minutes?

proxy.config

<ProxyConfig allowedReferers="*" logFile="proxylog.txt" mustMatch="false">

    <serverUrls>

        <serverUrl url = "http://arcgiserver.com/"
            matchAll="true"
            tokenServicesUri="http://arcgiserver.com/my_web_adaptor/tokens"
        />

    </serverUrls>
</ProxyConfig>


texts: Generating TIN, manually modified level curves?

texts:
I'm having trouble while generating a TIN on a manually modified topography file.

I'm new to ArcGIS, but I have intensively tried to fix this bug, with no result.

As you can see, the bug deals about continuity in the generated model, the level curves seems to be logically disposed but the generated model proposes a big hole in the ground.

I noticed that the bug (unwanted changing of color = uncontinuity) was triggered by uncontinuous lines, probably created during the manual modification of this curves.

I tried to merge, union, dissolve ans unsplit the lines with no effect on the generated model.

Thus, I'm currently lost and I would like to know if anyone already noticed this bug, and maybe have a piece of solution to give me ?

I'm using on ArcGIS 10.1 for Desktop.


texts: Extracting values from .ovr pyramid files provides NA in R? getValues?

texts:
I'm trying to extract values from many rasters to specific points and have those values be in columns that are labeled as the raster name. I can do it all here except get the values. All projections match. Is this due to the fact that is an ovr type raster?

m <- "E:/./"
allFiles <- list.files(m, full.names = T,pattern = "_ex0315.ovr",
                       no.. = TRUE, recursive=TRUE)
filenamez <- str_split_fixed(allFiles, "/", 6)[,5]
allFiles <- stack(allFiles)

names(allFiles) <- filenamez
projection(allFiles) <- CRS("+proj=utm +zone=14 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

Choruspts<-read.csv(".")

corproj<-CRS("+proj=utm +zone=14 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

points <- SpatialPointsDataFrame(Choruspts[,2:3], Choruspts,
                                 proj4string = corproj)

points.mrc <- spTransform(points, CRS(projection(allFiles)))

output<-as.data.frame(matrix(NA, nrow(Choruspts), length(filenamez)))

dim(output)
for(i in 1:length(filenamez)) {
    output[, i] <- extract(allFiles[[i]], points, method="bilinear")
}
names(output) <- filenamez

fulloutput <- cbind(Choruspts, output)

The output has the 3 columns of the point files and column headers of the 399 rasters but the values are all NAs.

dim(output)
[1] 200 399

For example:

head(fulloutput)
  Location Easting Northing A1_2013-04-06_650 A1_2013-04-06_710 A1_2013-04-06_730 A1_2013-04-06_750 A1_2013-04-20_640
1 A1 E 100  427572  4676098                NA                NA                NA                NA                NA
2 A1 E 200  427672  4676098                NA                NA                NA                NA                NA
3 A1 E 400  427872  4676097                NA                NA                NA                NA                NA
4  A1 E 50  427522  4676099                NA                NA                NA                NA                NA
5 A1 E 800  428073  4676090                NA                NA                NA                NA                NA
6 A1 N 100  427465  4676203                NA                NA                NA                NA                NA

I have determined that the stack is lacking values, which should be the last entry on the bottom:

class       : RasterStack
dimensions  : 1436, 801, 1150236, 399  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0, 801, 0, 1436  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
names       : A1_2013.04.06_650, A1_2013.04.06_710, A1_2013.04.06_730, A1_2013.04.06_750, A1_2013.04.20_640, A1_2013.04.20_700, A1_2013.04.20_740, A1_2013.04.20_800, A1_2013.04.20_820, A1_2013.04.20_840, A1_2013.04.20_900, A1_2013.04.20_920, A1_2013.04.27_617, A1_2013.04.27_644, A1_2013.04.27_708, ...


texts: Error message exporting a raster in QGIS

texts:
I have created a raster using interpolation-TIN.
I tried to export it but I get the error message:

  Invalid creation option:
  PREDICTOR=2 only valid for 8/16/32 bits per sample (using 64)
  Click on help button to get valid creation options for this format

How can I solve this?


texts: When I try to load this WMS: https://wms.geo.admin.ch/?VERSION=1.0.0 in QGIS, I get the following error message:"mswmsloadgetmapparams() image handling error. unsupported output format". In ArcMap the WMS works without issues. Any idea how to fix this?


texts: Extract derived data fields from Shapefiles

texts:
I've uploaded shapefiles with professionally supplied survey data into QGIS and it all plots fine, then when I export the attributes table into Excel I see that the XYZ values are all zeros.

I know they can't be zero as all the data points are plotting on my map. then when I select each point and view attributes table, sure enough, the XYZ fields are zero, however I do see a tab labelled derived, which when selected does show all of the XYZ attributes that I require.

Does anyone know how to extract these derived data fields from shapefiles so I can merge the values with the previously exported attributes Excel file and have all the survey information I need within a single file.

Extra information added below on 13/03/19.

the data i have is shape files extracted from civil cad that contains survey information of drainage pipes and pits. it has x,y,z and other comments. everything appears in the attributes tables except the xyz information which can be viewed in the derived tab, but not the attribute table. i would be happy to show you screenshots if i could post them.

the attributes table is produced in excel by the following steps
1. right click on the shapefile in layers list
2. select export-save feature as
3. then select excel file from list and follow the prompts

the result is a file with copy of the attributes from the shapefile. but the upstream and downstream E and N coords are zero.

the x,y coords i found when i use identify and right click - then select a pipe or pit from the displayed map. the identify results show feature and value columns where the attributes are displayed and also a derived tab.
under the derived tab when expanded, there is the x and y coords for pit or start and end of pipe for whichever feature is selected.

im trying to get the coords under the derived tab extracted and shown in the attributes table as upsteam or downstream E and N coords.

so far ive only been able to do this by hand one at a time, but there must be a better way.


texts: Show popup on marker hover / mouseover, hide on mouseout, and keep open on click

texts:
I'm trying to make a marker popup in Leaflet (1.3.1) that:

Appears on mouseover ("hover")
Disappears on mouseout
Persists (stays open) on click

Here's what I have so far:

var marker1 = L.marker([51.505, -0.09]).addTo(mymap)
  .bindPopup("Hello world!
I am a popup.");

marker1.on('mouseover', function (e) {
  this.openPopup();
});
marker1.on('mouseout', function (e) {
  this.closePopup();
});
marker1.on('click', function (e) {
  this.openPopup();
  //disable mouseout behavior here?
});

Plunker: https://embed.plnkr.co/393lgE49Ndqsj7O6dg19/

Mouseover and mouseout work as expected. When I click on the marker, the popup briefly disappears before reappearing. When I mouseout, the marker disappears.

I'd like it if the popup didn't briefly disappear when I click on it, and instead just stayed open until the user clicks the close 'x' button or clicks somewhere else on the map.


texts: Converting coordinates into decimal using SQL Server Spatial?

texts:
I am looking to convert a pair of coordinates to decimal, and am unsure how to go about doing it.

The following script

select Centroid.STY as Longitude,
Centroid.STX as Latidude
from
(    select Shape.STCentroid() as Centroid
     from gis.USA_10_CENSUS_TRACT
)U

generates these results:

-1125114.67826638   -2502243.88901991
-1054751.90111223   -2514282.59232887
-1048731.88354633   -2500496.88304733
-1046266.95217362   -2499166.70651986

When I am in fact looking for results in the following format:

34.54270600 -112.46910900
33.37959400 -111.76312100
31.64200000 -110.33700000
33.75689700 -111.99271000

(note that the above and bottom examples are not intended to match)


texts: How to merge same long lat points with different data into a unique one?

texts:
I have a layer where there is more than 1 feature for a same object. Imagine i have  a layer called signals. There can be up to 10 features in the same place (long/lat) for a unique signal but the features can have different fields or not. I want to traspose (ms access term) those fields from vertical to horizontal (1 unique row). Is this possible in QGIS?

Imagine the following geographic features:

Feature 1
Long    Lat   Field1    Field2  Field3
x        y    a         b       c

Feature 2
Long    Lat   Field1    Field2  Field3
x       y     a         b       d

What i need is to obtain:

Feature
Long    Lat   Field1    Field2   Field3_1   Field3_2
x       y     a         b        c          d

Is this posiible in QGIS?


texts: QGIS Points to Path Plugin Date/time error

texts:
I would like to use the Points to Paths plugin to plot the tracks of collared animals, exactly as specified on this page http://spatialgalaxy.net/2011/12/30/qgis-plugin-of-the-week-points-to-paths/

However, my date/time is in this format '2013/06/06 00:25' and the plugin reports this error:

ValueError: Order field is not an integer type or date format is invalid

I can get the plugin to work using time in HH:MM:SS only (so it's not my order field at fault) but as the animals are tracked over several days this is not practical. I cannot find in any of my three spreadsheet programs a date/time format option that uses '-' as the date separator as shown in the box? Does the plugin support my date/time format, and if so, what syntax should I use in the plugin's date format box?
Or can anyone tell me how to generate a CSV file with date/time in the 2013-06-25T00:25 format that the plugin shows?
Thanks!
PS Am generating the date/time field by adding together separate date and time columns in the raw data


texts: Using OGR_GEOM_AREA to add area attribute to shapefile?

texts:
I've been trying to batch convert raster files to polygons then delete areas <100m2 before clipping the raster to the filtered polygons.

I have referred to both the below posts:

Is there any way to add an area column to new shapefile created from raster with gdal's polygonize in the same step?
https://stackoverflow.com/questions/42355370/sqlite-query-on-path-as-directory-ogr2ogr

I have tried to use the desired command in its simplest form:

ogr2ogr ouput.shp input.shp -sql "SELECT *, OGR_GEOM_AREA AS area FROM input"

Despite this I still get the error "ERROR 1: SELECT from table input failed, no such table/featureclass."

I know that the input.shp exists and is in the same directory as my batch file, so I can only assume that the OGR_GEOM_AREA can not be found.

I notice that the gdal_polygonize step I use prior writes the output in GML format: "Creating output input.shp of format GML."

Could this be the source of the issue?

If so I have not been able to find how to change the polygonize output format to something more ogr friendly.

I did some tests and can confirm that the issue is definitely the GML format being output by gdal_polygonize


texts: ST_Covers fails on GEOGRAPHY column with GiST index, PostGIS 3.0

texts:
I'm trying to upgrade to Postgres 12 with PostGIS 3.0. Unfortunately, I can't get my spatial queries to work. There seems to be an issue using ST_Covers on a geography column with a GiST index.

Proof of Concept:

CREATE TABLE test (id SERIAL, the_geog GEOGRAPHY);
CREATE INDEX test_geog_idx ON test USING GIST(the_geog);
INSERT INTO test (the_geog) values (
    ST_Buffer(ST_GeogFromText('POINT(-0.1 51.5)'), 10)
);
SELECT id FROM test WHERE ST_Covers(the_geog, ST_GeogFromText('POINT(-0.1 51.5)'));

ERROR:  no spatial operator found for opfamily 19745 strategy 7

I imagine this is because GiST indexes work in a different way in PostGIS 3.0. How can I solve this issue?

Here is my Postgres version:

test=# select version(),postgis_full_version();
                                                   version                                                   |                                                                        postgis_full_version
-------------------------------------------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------
 PostgreSQL 12.1 on x86_64-apple-darwin16.7.0, compiled by Apple LLVM version 8.1.0 (clang-802.0.42), 64-bit | POSTGIS="3.0.0 r17983" [EXTENSION] PGSQL="120" GEOS="3.8.0-CAPI-1.13.1 " PROJ="6.2.1" LIBXML="2.9.10" LIBJSON="0.13.1" LIBPROTOBUF="1.3.2" WAGYU="0.4.3 (Internal)"
(1 row)

Here is the equivalent code with geometry:

CREATE TABLE geom (id SERIAL, the_geom GEOMETRY(Polygon, 4326));
CREATE INDEX test_geom_idx ON geom USING GIST(the_geom);
INSERT INTO geom (the_geom) values (
    ST_Buffer(ST_GeomFromText('POINT(-0.1 51.5)'), 10)
);
SELECT id FROM geom WHERE ST_Covers(the_geom, ST_GeomFromText('POINT(-0.1 51.5)',4326));

 id
----
  1

EDIT: Here is a snippet with a table named test1:

sx=# CREATE TABLE test1 (id SERIAL, the_geog GEOGRAPHY);
CREATE TABLE
sx=# CREATE INDEX test1_geog_idx ON test1 USING GIST(the_geog);
CREATE INDEX
sx=# SELECT id FROM test1 WHERE ST_Covers(the_geog, ST_GeogFromText('POINT(-0.1 51.5)'));
ERROR:  no spatial operator found for opfamily 19505 strategy 7


texts: Bilinear resampling with GDAL leaves holes

texts:
I'm attempting to down-sample a raster (resampling to a coarser/ larger pixel size) with continuous, floating-point values. I generally use bilinear resampling for the up-sampling case (resampling to a finer resolution), so I thought I would use it here, too. I'm finding that whether I use gdal.ReprojectImage() or gdal.Warp (which is a wrapper for gdalwarp), I get large NoData chunks in the output image wherever there was at least one NoData value in the input image. I would have assumed that as long as there was some data in an output pixel (which subsumes multiple smaller pixels), there would be a valid output value.

The effect is quite striking visually. The original raster image (9-km resolution on an equal-area grid) looks like this:

The resampled raster image (1-degree equirectangular projection) looks like this:

Sample Python code used to resample the original equal-area raster is below; there is an equivalent with gdal.ReprojectImage() but it is much more verbose.

import gdal
from gdal import gdalconst

ds = gdal.Warp(
    'temp.file', original_raster,
    format = 'MEM', xRes = 1, yRes = 1,
    dstSRS = "EPSG:4326",
    outputBounds = (-179.5, -90.5, 180.5, 89.5),
    resampleAlg = gdalconst.GRA_Bilinear,
    outputType = gdalconst.GDT_Float32)

In contrast, a nearest-neighbor resampling (changing resampleAlg argument) produces a continuous 1-degree grid, with no "holes" in the output raster. Nearest neighbor is much simpler, of course, but is there an option to calculate weights for bilinear resampling in a way that NoData pixels don't contaminate output pixels in this way?


texts: QGIS 3.0 - MMQGIS - Geocode. Inconsistent results

texts:
I have a .csv file I am trying to geocode using the "Geocode CSV with Google/Open Streetmap" tool in MMQGIS and I am getting inconsistent results using the exact same file with the exact same parameters.  In any case, the results aren't what they should be.  My .csv file only has 11 rows, so it is not a huge file.

I ran the tool 3 times using these parameters:

and I received these results...
First time:

Second time:

Third time:

Finally I tried to change the "web service" to Open Streetmap and received an error with no results. The results seem completely arbitrary and none of them are the desired result.  Can anyone explain this?


texts: interpolate raster values to the extent of a bigger vector layer in QGIS 2.18.16

texts:
I have raster data (1km resolution) of nationwide multiannual precipitation. Unfortunately, the raster layer does not completely cover the extent of an underlying vector layer (see picture below).

I would like to interpolate the original raster values to a slightly bigger raster layer, which completely covers the vector layer.

My approach (conversion of raster values to points and interpolation of point values to a new raster layer with a bigger extent) has not been very straightforward. I guess, there is a better way.


texts: Checking whether a vector layer is empty or not

texts: I am working on intersections of a number of buffers with one another. When two particular buffers have no intersection, applying intersection function on them produces an empty layer which I want to delete. But how do I check if the intersection result produced is an empty layer or not? Is there any inbuilt function which I can use in my code?


texts:
I am trying to use XYZ Esri tile source with Openlayers for a cached service in Palestine Grid Coordinate System (EPSG:28191) as in this example. However openlayers fails request incorrect tiles which are always blank.
My cached ArcGIS service is https://portal.geomolg.ps/server/rest/services/AerialPhoto/AerialPhoto2019/MapServer
What might be the issue here?
Here is my code:
import 'ol/ol.css';
import Map from 'ol/Map';
import TileLayer from 'ol/layer/Tile';
import View from 'ol/View';
import XYZ from 'ol/source/XYZ';
import {fromLonLat} from 'ol/proj';

var map = new Map({
  target: 'map',
  layers: [
    new TileLayer({
      source: new XYZ({
        attributions:
          'Tiles © ArcGIS',
        url:
          'https://portal.geomolg.ps/server/rest/services/' +
          'AerialPhoto/AerialPhoto2019/MapServer/tile/{z}/{y}/{x}',
      }),
    }) ],
  view: new View({
    center: [-4839400.0, 6616800.0],
    zoom: 7,
  }),
});

With this configuration, Openlayers requests tiles that are not available like:
https://portal.geomolg.ps/server/rest/services/AerialPhoto/AerialPhoto2019/MapServer/tile/6/21/23


texts:
I am having problems adding postgis to my database.  I'm following the instructions from here: http://postgis.net/docs/postgis_installation.html#idp32328352  (That said, I have combined some previous instructions from multiple blogs/advices.)

After creating my sample gis db, this is what I'm seeing:

/c gis
gis=# CREATE EXTENSION postgis;
ERROR:  PostGIS is already installed in schema 'public', uninstall it first
gis=# CREATE EXTENSION postgis_topology;
ERROR:  required extension "postgis" is not installed
gis=# CREATE EXTENSION postgis_tiger_geocoder;
ERROR:  required extension "postgis" is not installed

I think that I may have originally done a CREATE EXTENSION postgis; on the postgres database.  I'm not sure.  I have tried DROP EXTENSION postgis; on postgres, but that returns: ERROR:  extension "postgis" does not exist.

I have gone so far as to purge postgres and related packages, then reinstalled it.  (Though it seemed to still find the gis db after reinstalling -so I'm not sure it was thorough.)
Am I misunderstanding the message that postgis is installed in schema public?  The sites I have read have advised to install it separately on each gis-enabled DB.

Add'l information: SELECT version(); Returns: PostgreSQL 9.4.4 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu 4.8.2-19ubuntu1) 4.8.2, 64-bit

And:  SELECT postgis_full_version();  Returns: POSTGIS="2.1.8 r13780" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.9.1" LIBJSON="UNKNOWN"
This is a development setup.


texts: Problem with labels on export/print a map (qgis master)

texts:
As you can see on the images below, the labels on the exported map (right) are much bigger then in the original map (left). Same by pdf, image, and vector - export.

Whats the problem? Can I fall back in the qgis master version to the version of about 1 month ago?


texts: import gdal, return "ImportError: DLL load failed:"

texts:
I came up against this problem during the installation of gdal, I tried several solutions already read on the site but none helped me. Here's what I did:

1) uninstall gdal
2) uninstall Anaconda
3) install Anaconda
4) conda create new env("GDAL")
5) (GDAL) conda install conda-forge gdal

What am I supposed to do?


texts: GeoServer on cascading WMTS sets wrong CRS

texts:
I have one GeoServer 2.12.5 with a layer which CRS is set to EPSG:3857, and another GeoServer 2.14.2 where I want to cascade WMTS from the first one. But when I add this layer from cascaded source its native CRS is detected as EPSG:4326, not 3857. I even deleted 4326 gridset from tile options of original layer, but is didn't help. How to fix cascaded CRS?

Original layer group CRS-settings:

Original layer group gridset settings:

In cascaded layer its CRS is detected as EPSG:4326


texts: Rasterio read function producing uniform array or near-zero values for a non-uniform raster

texts:
I have been looking all over the internet for a possible cause of this, but without luck, and was hoping it's some standard Rasterio quirk that I just can't track down. I am trying to load the Tanzania poverty map from World Pop, which is a 1km resolution raster grid of population rates (available here). Values for each grid square can be between 0 and 1, and the raster looks like this:

Using the following Python code I've been loading the raster with Rasterio:

import os
import rasterio

file_path = os.path.abspath('path/to/tza10povcons125.tif')
raster = rasterio.open(file_path)
band = raster.read(1)

Based on my understanding of Rasterio, I would expect the band variable to contain a 2 dimensional array of values between 0 and 1 (or, rather, 0.664994 and 0.950991, which are the min and max values according to QGIS). However, when I print out the band variable I see:

array([[ -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38, ...,
         -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38],
       [ -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38, ...,
         -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38],
       [ -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38, ...,
         -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38],
       ...,
       [ -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38, ...,
         -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38],
       [ -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38, ...,
         -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38],
       [ -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38, ...,
         -3.39999995e+38,  -3.39999995e+38,  -3.39999995e+38]], dtype=float32)

The recurring value, -3.39999995e+38, is also suspiciously close to the minimum value for numpy's float32 datatype:

>>> import numpy as np
>>> np.finfo('float32').min
-3.4028235e+38

I am just wondering if anyone has seen this kind of behaviour for Rasterio before, or if there is something obvious that I am missing?


texts: Is there a way to open BSQ file format in IDL?


texts: Creating groups containing slopes of a specific range

texts:
I have a basic knowledge of QGIS and can do everything I need to do to achieve my requirements except the following.

I'm trying to graphically show areas of land that contain a set of specific slope ranges.

I'm working in relatively small areas of up to 50ha.

I have the DEM and can create contours with no problem.

What im trying to do is shown in the attached image, I know it can be done but I just don't know how to do it or what the process is called.


texts: How do I cast shadows from 3D Objects in ArcScene?

texts: I would like to cast shadows from my 3D objects in ArcScene onto the surface and other objects. (Much like I can in SketchUp).  Is it possible? If so how?


texts: Geoprocessing with ArcGIS Server 10.0 - %scratchworkspace%

texts:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Key_concepts_for_geoprocessing_services/002v00000002000000/

I'm trying to have a new GP Service arcpy.AddMessage() the scratchworkspace that is created on that run.  I keep getting a invalid syntax error.

But real aim is to write to this scratchworkspace but the %scratchworkspace%/ doesn't seem to be working for me.

line = arcpy.GetParameterAsText(0)
arcpy.AddMessage(line)
try:
    arcpy.AddMessage(str(%scratchworkspace%))
except:
    arcpy.AddMessage(traceback.format_exc())

This is a gp service published through the arcmap process, I've opened and run the tool layer once already and the Tool is creating the scratchworkspace in arcgisjobs directory

arcgisjobs\geoprocessing\gpservicetest_gpserver\jf93e76ce2f0f4967b7b9028c384b19f8\scratch

So how do I get a hold of this in python to write to?

Thanks.


texts: Changing symbology of all layers in table of contents using ArcPy and ArcMap?

texts:
I am attempting to put together a DEM of the U.S. I have thousands of contour elevation files that I have loaded into an ArcMap document. I need to change the symbology for all of the layers in order to represent variations in elevation. I have an ArcPy script that will change the symbology of a list of layers, based on the symbology of a layer assigned as the template.

How do I change the symbology of all of the layers in the table of contents at once?


texts: How can I save a column from the attribute table of a vector layer to a variable?

texts:
I am trying to write a python plugin for QGIS that can convert three vector layers into an object File for an urban climate model. Unfortunately I am stuck again with a possibly easy question:

How can I save a column from the attribute table of a vector layer to a variable (for later use)?

For example, I have a layer with buildings (build) with the hight of the buildings in field 35. I want to save all the building hights into a variable (theINPolyHField).

I tried the following:

build = canvas.layer(2)
build_prov = build.dataProvider()

theINPolyHField = build_prov.select([35], build.extent(), False, False) # but returns "None"
build.selectedFeatureCount() # returns 0

What am I doing wrong?


texts: Connecting ordered MultiPoint features with lines/arrows

texts:
I'm working through an XY problem, with the end goal of effectively visualising data representing different journeys.

Given I have ordered data of the various positions (each with a time stamp) along the journey, I am considering using a MultiPoint layer, with each feature representing a journey, and each point a position along that journey.

I have written minimal (but verbose to convey intent) code below to generate a MultiPoint layer and visualise some dummy journeys as points.

My question: Is it possible to visually connect each pair of adjoining points with a line (or better still an arrow)?  If so, how would I extend the following code to do this using PyQGIS?

I need to automate this, at the time the layer is created, and I don't want to rely on external plugins.

    from collections import namedtuple

    from qgis.core import QgsFeature, QgsField, QgsGeometry, QgsPointXY, QgsProject, QgsVectorLayer
    from qgis.PyQt.QtCore import QVariant

    # Create a new multipoint layer
    layer = QgsVectorLayer('MultiPoint?crs=EPSG:4326', "Journeys", 'memory')
    layer.startEditing()

    # Define attribute table headers
    field_headers = [
        QgsField('journey_id', QVariant.String),
        QgsField('positions [(lat, long, time_step), ...]', QVariant.String),
    ]
    layer_provider = layer.dataProvider()
    layer_provider.addAttributes(field_headers)

    # Define sample data for 2 journeys
    Position = namedtuple('Position', 'lat long time_step')
    journeys = [
        ('journey_1', [Position(10, 0, 1), Position(1, 1, 2), Position(3, 0, 3)]),
        ('journey_2', [Position(-5, -5, 1), Position(2, 2, 2), Position(-4.5, 9, 3)]),
    ]

    # Create feature for each journey
    features = []
    for journey in journeys:
        journey_id, positions = journey

        # Create points list from positions
        points = [QgsPointXY(p.lat, p.long) for p in positions]

        # Create feature from points list
        feature = QgsFeature()
        feature.setGeometry(QgsGeometry.fromMultiPointXY(points))

        # Add entry to attributes table
        feature.setAttributes([
            journey_id,
            ", ".join(f"({p.lat}, {p.long}, {p.time_step})" for p in positions),
        ])

        features.append(feature)

    # Add features to layer, update, and commit
    layer_provider.addFeatures(features)
    layer.updateExtents()
    layer.commitChanges()

    # Add layer to legend
    QgsProject.instance().addMapLayer(layer)


texts: Can ArcGIS mapping be moved to Google Maps?  Looking to migrate from ArcGIS to Google Maps and wondering about the challenges/risks.


texts:
Consider the following query:

WITH geoms AS(
    SELECT
    ST_Buffer(ST_MakePoint(10,10),5,200) as circle,
    ST_MakeEnvelope(5,5,15,15) as square
)
SELECT ST_Snap(square,circle,15) geom FROM geoms;

I would expect two possible outcomes:

The square being made completely into a circle (points are being introduced to the segments)
The 4 vertices of the square are snapped to their closest neighbours on the circle (the square would fit into the circle)

but the actual outcome is this:

The other way around makes more sense (ST_Snap(circle, square,15)):

Here, only 4 vertices are snapped to their corresponding nearest square neigbours.

I don't understand what happens in the first case. Somebody else does?
(my real use-case is more complex, but seems to boil down to this problem)


texts: Leaflet: flyToBounds location of bounding box returned by PostGIS

texts:
I want to simply center my map on the search result of a polygon using leaflet. I thought to return the bounding box, as this is closer to what I actually need.

I want to use flyTo()

I see a possible way to do this by using Leaflet?

I have a thought to do this

SELECT ST_AsGeoJSON (geom) AS geom, Box2D(geom) as bounds, ST_YMax(Box2D(geom)) as ymax, ST_YMin(Box2D(geom)) as ymin, ST_XMax(Box2D(geom)) as xmax, ST_XMin(Box2D(geom)) as xmin FROM nation_shapes WHERE (condition)

and then pass it to

mymap.flyToBounds([[ymax, xmax],[ymin, xmin]])

Though this seems inefficient. I would like to be able to get this information from the object I am already passing back if possible, or convert something like Box2D(geom) as bounds which would be more efficient than calling it 4 times, I think.


texts: Transforming between epsg:3857 and Google Maps Tile coordinates

texts:
Hio,
I used to use a simple python project (https://github.com/hrldcpr/mercator.py) to convert between wgs84 and google maps tile coordinates. This project doesn't have a license, and looks like it may be inaccurate (assuming a spherical earth?), so I decided to use pyproj.

I understand that google maps uses a projection called epsg:3857 to draw their tiles. When I convert from wgs84 to epsg:3857, I get numbers that I just don't understand:

>> import pyproj
>> epsg3857 = pyproj.Proj(init='epsg:3857')
>> wgs84 = pyproj.Proj(init='EPSG:4326')
>> latlng = (45.618, -73.604)
>> pyproj.transform(wgs84,epsg3857,latlng[0],latlng[1])
(5078172.5310075525, -12357511.560866801)

I know that the lat,lng coordinate I took is in the tile http://b.tile.openstreetmap.org/11/605/731.png, But I just can't figure out how to go from the returned numbers to the tile coordinates. Firstly, the numbers shouldn't be negative (tile coordinates are positive). Trying different scales based on the 11 zoom levels, I just can't figure this out. Searching the internet, I can't figure out what the epsg:3857 coordinates actually denote.


texts: Import JavaScript Google Earth Engine function scripts in a Python API script

texts:
I have external scripts in the JavaScript API for the purpose of keeping the functions needed for a mass export of data separate. For example, when I would like to reference these functions in other scripts in the main script, I would use something like:
var codes = require('users/mincej20/path:To/Script');
var cloud = require('users/mincej20/path:To/CloudScript);

The problem is that the JavaScript API is not well managed when it comes to exporting many different data tables. I have about 100,000 lakes that I need to create a .csv for, and it seems that the Python API method ee.batch.export in a loop would suit this much better.
Am I able to reference my already-made function scripts in a Python Google Earth Engine API script similar to referencing an asset? If so, how would I do this? Or will I need to re-create these scripts locally in Python scripts?


texts: Projection for Germany that preserves distance?

texts:
I have a point A (EPSG 4326) in Berlin, Germany:

{
   "type": "Point",
   "coordinates": [
     13.394994735717773,
     52.514849775310154
   ]
}

I want to create point B that is 500m to the east of point A and point C that is 500m to the North.

For this, I have to project point A to a projection that preserves distance.

Can you recommend a projection that works for this operation in entire Germany?

I am working with GeoPandas and Shapely for this project.


texts: mapbasic - how to link my database in my mapbasic coding

texts:
I'm creating a tool in MapBasic that identifies the PCI clash in 4G
I need help:

How do I add database data into MapBasic?
How can I calculate mod values and colour the output in MapInfo based on these calculations?


texts: List fields in parameters from Fc list in script validation using arcpy in arcgis

texts:
I am trying to list fields from a feature class list in script tool validation. First user selects the data frame from the list then FC populates in next box and fields of the FC in next box. I have code up to populating the FC list and working fine. but after that i'm trying list fields from the FC and getting error.

here is the code

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    try:   # this works
        mxd = arcpy.mapping.MapDocument('CURRENT')
        dflist = [df.name for df in arcpy.mapping.ListDataFrames(mxd)]
        self.params[1].filter.type = "ValueList"
        self.params[1].filter.list = dflist
        del mxd
    except:
        arcpy.AddError("\nMust run this tool in ArcMap, not ArcCatalog")

    self.params[2].enabled = True

    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    mxd = arcpy.mapping.MapDocument('CURRENT')

    if self.params[1].value:
        self.params[2].enabled = True

        df = arcpy.mapping.ListDataFrames(mxd, self.params[1].value)[0]
        # creates an empty dictionary
        if df:
            layerList = [lyr.name for lyr in arcpy.mapping.ListLayers(df) if not lyr.isGroupLayer]
            uniqueList = list(set(layerList))
            uniqueList.sort()
            self.params[2].filter.list = uniqueList

# Up to this it working and problem is from below code

    if self.params[2].value:
        self.params[3].enabled = True

        df = str(self.params[1].value)
        lyrList = arcpy.mapping.ListLayers(mxd, "", df)
        for lyr in lyrList:
            if lyr.supports("dataSource") and lyr.name == str(self.params[2].value):
                for field in arcpy.ListFields(lyr.dataSource):
                    uniqueList = list(set(field))
                    uniqueList.sort()
                    self.params[3].filter.list = uniqueList

    return

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return


texts: My QGIS project consists of two adjacent Canadian NTS map sheets. Each sheet consists of multiple shape files. Each set of shape files has a border.  How do I eliminate the border?


texts: I am having a problem with adding a layer group into QGIS, where upon making a connection to the server in QGIS, the layer group does not appear. I see only the layers that make up that layer group. Can someone tell me what might be the problem or what I might be missing?


texts: Summing mulitple GridCoverage Objects using GeoTools?

texts:
I'm working in Java and GeoTools to process some ASCII rasters. At present, I'm trying to get the rasters summed onto one large grid. I've tried a couple of methods but I seem to be having some performance degradation.

(These are manually copied so apologies for typos)

I initially tried:

Operatons op = new Operations(null);
GridCoverage2D sum = initialise();
for (GridCoverage2D grid : grid){
    sum = op.add(sum, grid);
}

This rattled through seemly quickly but when I try to access sum is takes ages, so I assume these operations are conducted lazily. I then tried to add a writer in between to access the grid so it's calculated on the fly.

Operatons op = new Operations(null);
GridCoverage2D sum = initialise();
for (GridCoverage2D grid : grid){
    sum = op.add(sum, grid);
    GridCoverageWriter writer = new ArcGridWriter(File.createTempFile(...))
    writer.write(sum);
}

Which improved things a little. I then thought doing a recursive sum might help, but I'm getting continuously slower results still and hammering a single processor. Are there any ways of increasing the speed of these operations or am I approaching this in the wrong way?

The grids are approximately 4000x4000. At the moment summing 6 grids takes 5 minutes, 8 takes about 15, 10 takes approaching half and hour. The grids are generated using a fortran model that I don't have access to and there could be up to approaching 100 grids being added together. I have resampled the grids so they all fall in the same envelope. The only problem with the first approach is the runtime, but I was wondering if there would be a more efficient way of doing it?


texts:
I Want to Add Database Server in ArcCatalog. My Database is Sql Server 2014 enterprise;When I try, It returns that "Connections to this edition of the database Server are not Supported". I Searched and found that i should install proper version of Arcsde. But I couldnt find ArcSde for Sql Server 2014 enterprise. How can solve this problem?
(ArcGis 10.2.2 ; Sql Server 2014 enterprise; Sql Client is installed;)


texts: KeyboardInterrupt doesn't work when arcgisscripting module is loaded

texts:
The KeyboardInterrupt exception -- thrown when Ctrl-C, or Ctrl-Break on some machines, are pressed in a Python console window -- is commonly used to gracefully break out of long-process loops, allowing your script to perform cleanup, logging etc, before exiting.

Under my ArcGIS 9.3.1 SP2/Windows XP SP3 environment, there is an issue with catching the KeyboardInterrupt exception if the arcgisscripting module is loaded. Instead of hitting the KeyboardInterrupt catch statement, the script aborts abruptly, spitting out the following message:

forrtl: error (200): program aborting due to control-C event

If the arcgisscripting module is not loaded (removed from the import statement), the KeyboardInterrupt exception is thrown and processed properly.

From what I've read this is due to the use of Fortran in the arcgisscripting module which somehow overrides Python's interrupt processing.

Here's an example script you can use to test and reproduce the issue. Save it as a Python script, run it from the Windows command prompt and press Ctrl-C (Ctrl-Break on some machines) to test whether the KeyboardInterrupt exception is processed. Simply toggle the commenting for the import arcgisscripting line.

#import arcgisscripting
import time

try:
    while(True):
        time.sleep(1)
except KeyboardInterrupt:
    print "KeyboardInterrupt detected!"

Anyone know of a workaround to this problem? Thanks in advance for any insights you may have.

edit: here is my test script using Jason's suggestion which unfortunately doesn't work:

import arcgisscripting
import signal
import time

def ctrlc(sig, frame):
    raise KeyboardInterrupt("CTRL-C!")

print "Waiting 65s ..."
time.sleep(65)
print "Installing SIGINT handler"
signal.signal(signal.SIGINT, ctrlc)

try:
    while(True):
        time.sleep(1)
except KeyboardInterrupt:
    print "KeyboardInterrupt detected!"


texts: Selecting only 10000 records from layer using PyQGIS

texts:
I would like to join attributes by location of 2 layers. My INPUT layer has 446k records and it is too many. I would like to divide my layer and select for example 10000 records.
For example in Python if I have a list I can select by data[:10000]. How can I select x records from layer?
My code:
from qgis.core import *
from qgis.utils import *
from qgis.analysis import QgsNativeAlgorithms
from PyQt5.QtCore import QVariant
from qgis.core import QgsApplication, QgsProcessingFeedback, QgsRasterLayer
import sys

sys.path.append('/usr/lib/qgis')
sys.path.append('/usr/share/qgis/python/plugins')
os.environ["QT_QPA_PLATFORM"] = "offscreen"

QgsApplication.setPrefixPath(r'/usr/bin/qgis', True)
qgs = QgsApplication([], False)
qgs.initQgis()

import processing
from processing.core.Processing import Processing

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
feedback = QgsProcessingFeedback()

processing.run("native:joinattributesbylocation", {'INPUT':r'/home/gis/budynki_opolskie/budynki_opolskie/bdot10k_bubd_a_opolskie.gpkg|layername=bubd_a',
                                                   'JOIN':r'/home/gis/adresypolska/PRG_PunktyAdresowe_POLSKA.shp',
                                                   'PREDICATE':[1,5],
                                                   'JOIN_FIELDS':['PNA','SIMC_id','SIMC_nazwa','ULIC_nazwa','Numer'],
                                                   'METHOD':0,
                                                   'DISCARD_NONMATCHING':True,
                                                   'PREFIX':'',
                                                   'OUTPUT':r'/home/gis/zlaczoneopolskie.shp'})


texts: Why is my function hiding all of my features?

texts:
I'm making a map, using ArcGIS javascript API, that has a filter tool that filters some polygons by habitat type. When the user selects a habitat type, only features of that habitat type should be left on the map, but for some reason, my function is hiding all features when I select different types. I'm making this map based off of a tutorial on the ArcGIS for javascript API website. When I test the value of "value" and "newValue" they have the correct value, but it's not translating into selective hiding.

Why is my filter function not changing features properly?

Here are the functions I'm using:

  require([
    "esri/map",
    "esri/layers/FeatureLayer",
    "esri/dijit/PopupTemplate",
    "esri/dijit/Legend",
    "dojo/domReady!"
  ], function(
    Map,
    FeatureLayer,
    PopupTemplate,
    Legend
  ) {

Here's is the code that creates the filter upon a change event and the function (updateDefinitionExpression) that uses a conditional statement to set the definition expression:

  var filter = document.getElementById("filter");
  // filters the layer using a definitionExpression
  // based on option selected by the user
  filter.addEventListener("change", function(event){
    var newValue = event.target.value;
    updateDefinitionExpression(newValue);
    alert(newValue);
  });
});

function updateDefinitionExpression(value){

  var definitionExpression = value ? "habitat type = '" + value + "'" : null;
  layer.setDefinitionExpression(definitionExpression);
  alert(definitionExpression);
  alert(value);
  map.infoWindow.hide();

}

And the HTML:

    Filter by habitat type:
    <select id="filter">
      <option value="">All</option>
      <option value="Mangrove">Mangrove</option>
      <option value="Pelagic">Pelagic</option>
      <option value="Reefs">Reefs</option>
      <option value="Salt Marsh">Salt Marsh</option>
      <option value="Sediment">Sediment</option>
    </select>


texts:
I'm trying to develop a plugin which grabs a layer from the project and shows its attribute table in a tableView. My code is as follows:

def getVectorLayer(self):
    layer = None
    layername = self.dlg.CBLayerInputs.currentText()
    for lyr in QgsProject.instance().mapLayers().values():
        if lyr.name == layername:
            layer = lyr
            self.dlg.tableView.QgsAttributeTableView(layer)
            break
        return layer

However, I do not get any error - the code compiles correctly without error, but no attribute table information is shown in the table view.


texts: merge two DEMs of different resolutions?

texts: I have two DEMs, the first one is in 30 meters of résolution and the second: 5 meters, and will like if there is a possibility of merging them to extend the area where I do my research?


texts: QGIS: How do I specify the layer for an get_feature expression in the model builder?

texts:
I'm trying to calculate a field (Field calculator 2) in the model builder with an get_feature expression which should look something like this attributes(get_feature('LandUse_TD_HYBAS_dissolved_calc', 'NEXT_DOWN', "HYBAS_ID"))['Calc'] . The problem is that I don't know how to specify the layer. Normally I just enter the layer name 'LandUse_TD_HYBAS_dissolved_calc' which worked as I did it outside of the model builder, but seemingly isn't working in the model builder because the layer I want to use is calculated by another field calculation one step ahead (Field calculator 1).

I have also tried to enter '@Field_calculator_1_OUTPUT' to specify the layer which did not work either(Picture 2).

What am I supposed to enter for the layer specification? Or is there another way to do it in the model builder?


texts: How to tile Mars themis raster for Cesium

texts:
I am trying to tile the Mars themis raster that has the following crs:
PROJCS["SimpleCylindrical Mars",
    GEOGCS["GCS_Mars",
        DATUM["D_Mars",
            SPHEROID["Mars",3396190,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",180],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

The --profile option in gdal2tiles can be geodetic, mercator or raster. Given the projection of the input raster above, I tile the raster as follows:
gdal2tiles.py --zoom=0-7 --profile raster input.tif output/

However, Cesium.js is unable to consume these tiles because it only supports a tiling scheme of either geodetic or mercator.
Is there a workaround? I was wondering if someone could kindly point me in the right direction.


texts:
I have a simple SpatiaLite DB containing one table with point features. I would like to generate a Minimum Convex Hull around that features. I could do that easily with one of the Convex Hull tools from GDAL/SAGA/GRASS/R-Script... from processing toolbox. But I would like to have a dynamic SQL solution. So every time a new point feature is added, the Convex Hull should automatically be updated. So I thought about a SQL-View:

I tried around following this instructions: http://postgis.net/docs/ST_ConvexHull.html

CREATE VIEW birdsmcp AS
SELECT
  ST_ConvexHull(ST_Collect(geom))
FROM birds GROUP BY bird;

I get a table, but no geometry. What do I have to do to get shapes as a result? I'm not very familiar with SQL queries.

EDIT: Thats my result in QGIS DB-Manager:


texts: Exporting time series of multiple polygons to .csv

texts:
I defined (by drawing) a significant number of polygons and need to export the time series of the average values inside each one. Following this example, that allows to display charts for each region, how can I export the data to a .csv? Ideally I will have in the first row the names of the regions.

// Define a FeatureCollection: regions of the American West.
var regions = ee.FeatureCollection([
  ee.Feature(    // San Francisco.
    ee.Geometry.Rectangle(-122.45, 37.74, -122.4, 37.8), {label: 'City'}),
  ee.Feature(  // Tahoe National Forest.
    ee.Geometry.Rectangle(-121, 39.4, -120.8, 39.8), {label: 'Forest'}),
  ee.Feature(  // Black Rock Desert.
    ee.Geometry.Rectangle(-119.15, 40.8, -119, 41), {label: 'Desert'})
]);

// Load Landsat 8 brightness temperature data for 1 year.
var temps2013 = ee.ImageCollection('LANDSAT/LC8_L1T_32DAY_TOA')
    .filterDate('2012-12-25', '2013-12-25')
    .select('B11');

// Create a time series chart.
var tempTimeSeries = ui.Chart.image.seriesByRegion(
    temps2013, regions, ee.Reducer.mean(), 'B11', 200, 'system:time_start', 'label')
        .setChartType('ScatterChart')
        .setOptions({
          title: 'Temperature over time in regions of the American West',
          vAxis: {title: 'Temperature (Kelvin)'},
          lineWidth: 1,
          pointSize: 4,
          series: {
            0: {color: 'FF0000'}, // urban
            1: {color: '00FF00'}, // forest
            2: {color: '0000FF'}  // desert
}});

// Display.
print(tempTimeSeries);


texts: Why is Shapely's snapping (GEO snaps) not working as expected?

texts:
I'm trying to snap two lines to each other using Shapely/Geopandas but the result of the snapping is very odd.
I tried :

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import snap

lines1 = gpd.GeoDataFrame.from_file('lines1.shp')
lines1 = lines1.to_crs({'init': 'epsg:2227'})
lines2 = gpd.GeoDataFrame.from_file('lines2.shp')
lines2 = lines2.to_crs({'init': 'epsg:2227'})
res = lines1
lines2_union = lines2.geometry.unary_union
res.geometry = res.geometry.apply(lambda x: snap(x, lines2_union, 14))
res.to_file('result.shp', driver="ESRI Shapefile")

And got this result :

lines1 = red lines

lines2 = black lines

After the snapping(with 14 as a tolerance) :
the blue lines are the result of the snapping

In this case the lines are correctly snapped

Another example where it didn't work as expected: (before the snapping)

And here's the result after the snapping. Only a part is snapped to the black line (the south side). Although the original lines are pretty close and within the 14 feets

If I increase the tolerance I get a wrong output, something like this (after defining 20 as the tolerance of the snapping, the green line is the result):

Any ideas on why the snapping is not working properly?
Any suggestions on how to solve this problem?


texts: Deleting rows in grid attribute table not allowed

texts:
I want to delete records in a raster attribute table. The window shows me an error message while deleting record.

  Grid attribute table is blocked from deletion.

I have attached a picture of the error.


texts: How to change color of points dynamically from a geoJson in openlayers

texts:
I am new in Openlayers and what I am trying to do is to dynamical give a color to some points read by a GeoJSON file. I am taking inspiration by this post

How can I add an infowindow to an OpenLayers GeoJSON layer?

but I can't get what is wrong in my code:

    <script type="text/javascript">

        var map;
        function initMap() {
            map = new OpenLayers.Map("info");
            var mapnik = new OpenLayers.Layer.OSM();
            map.addLayer(mapnik);

            map.addControl(new OpenLayers.Control.ScaleLine());
            map.addControl(new OpenLayers.Control.OverviewMap());
            map.addControl(new OpenLayers.Control.MousePosition());
            map.addControl(new OpenLayers.Control.LayerSwitcher());

             var myStyles = new OpenLayers.StyleMap({
                    "default": new OpenLayers.Style({
                        fillColor: "${getColor}",
                        strokeWidth: 1,
                        strokeColor: "#000",
                        fillOpacity: 0.8,
                        graphicZIndex: 5

                    },
                    {
                    context: {
                       getColor : function (feature) {
                          return feature.attributes.colore > 2 ? '#ffc000' :
                                 feature.attributes.colore > 1 ? '#00317c' :
                                                                 '#FFEDA0' ;
                         }
                       }
                  })
              });

            var geojson_layer = new OpenLayers.Layer.Vector("GeoJSON", {
                    styleMap: myStyles,
                    projection: new OpenLayers.Projection("EPSG:4326"),
                    strategies: [
                        new OpenLayers.Strategy.Fixed()],
                    protocol: new OpenLayers.Protocol.HTTP({
                        url: "colori.json",
                    format: new OpenLayers.Format.GeoJSON()
                })
            });

            map.addLayer(geojson_layer);

            map.setCenter(
                new OpenLayers.LonLat(8.692, 49.412).transform(
                    new OpenLayers.Projection("EPSG:4326"),
                    map.getProjectionObject()
                ), 10
            );

        }

    </script>

and here is the GeoJSON file colori.json I am trying to use:

{"type":"FeatureCollection",
    "features":[
        {"type":"Feature",
            "properties":{
                "colore": 7
            },
            "geometry":{
                "type":"Point",
                "coordinates":[8.692, 49.412]
            }
        }
    ]
}


texts: How do I import UK Environment Agency Lidar data (.asc) into QGIS

texts: I have downloaded Lidar data from the UK Environment Agency as ESRI ASCII files but being a complete novice as far as GIS is concerned I do not know how to view them. Quantum GIS was recommended to me so I have installed this. Any help gratefully received.


texts: Using ArcGIS Maximum Likelihood Classification gives Error 010108?

texts:
I am trying to find trees of a single species (Cedar) in my orthophoto.  My assumption is that the spectral bands of cedar trees are different to other vegetation (they definitely look different).  I have created about 8 small training areas and united them in one class which I have saved as a signature file.  I get two errors - one that the number of bands in my signature file is less than the number of bands in the raster, and the other is an error executing the grid expression.

Error 010108: The number of bands in the signature file is less than the number of bands in the raster
Error 010067: Error in executing grid expression

One factor that may be causing this is that my orthophoto is actually composed of a large amount of smaller tiles.  I have created my training areas and signature file within one of the tiles but I want to apply the classification to about 30 of my tiles.  I was getting an error message about extent so I changed the extent values to union of inputs in the environments options but I can't seem to get rid of the other error messages and I am clueless.


texts:
Hopefully a simple question.
I'm concatenating some fields in a featureCollection by running this .map function:
var fc = ee.FeatureCollection('users/xxx/featurecollection')
    .map(function (feature) {
        return feature.set('combined',
          feature.getString('ColA')
          .cat(' ')
          .cat(feature.getString('ColB'))
          .cat(' ')
          .cat(feature.getString('ColC'))
          );
    });

It works fine.  But, when I go to export to an Asset, I get the error "Error: Unable to export features with empty geometry"
Export.table.toAsset({
  collection: fc,
  description:'exportToTableAssetExample',
  assetId: 'MyExportedFC',
});

Question:  How do I skip features with Null/empty geometry in the above .map function?


texts: Combining grid of points to closest distance of several subqueries results using PostGIS

texts:
Say I have a table "table1" with a bunch of points (grid), and another table, "table2", with a bunch of lines (say they are roads), with a column "col1" (say it's a type of roads). So:

table1:

id;geography
1 ;POINT(X1,Y1)
2 ;POINT(X2,Y2)
...

table2:

geography  ; col1
LINE(Xi,Yi); 1
LINE(Xi,Yi); 1
LINE(Xi,Yi); 2
LINE(Xi,Yi); 1

I want to create a table which gives me, for each point in my grid, the distance to the closest "road" of type 1 and the distance to the closest "road" of type 2.

table3 I want (col1 is "distance to closest type 1", col2 is "distance to closest type 2"):

id ; col1  ; col2
1  ; 23.45 ; 0.234

To create this table, I thought of using lateral joins. However, doing so returns, in my case, zero rows, which is not true. I suspect that it joins both subqueries and that the intersection of the two is null. I'm just stumped as to how to get to what I want (which is basically the "sum" of the two subqueries instead of the intersection). Here is the query I used:

CREATE TABLE table3 AS (

SELECT t1.id,
ST_DISTANCE(t1.geography, t21.geography) AS col1,
ST_DISTANCE(t1.geography, t22.geography) AS col2
FROM table1 t1
CROSS JOIN LATERAL
    (SELECT t2.geography AS geography
     FROM table2 AS t2
     WHERE ST_DWITHIN(t2.geography, t1.geography, 10000) AND t2.col1 = 1
     ORDER BY ST_Distance(t2.geography, t1.geography)
     LIMIT 1
     ) AS t21
CROSS JOIN LATERAL
    (SELECT t2.geography AS geography
     FROM table2 AS t2
     WHERE ST_DWITHIN(t2.geography, t1.geography, 10000) AND t2.col1 = 2
     ORDER BY ST_Distance(t2.geography, t1.geography)
     LIMIT 1
    ) AS t22
)

I can do what I want for one lateral join at a time (i.e, have a table3 that has only col1), but I want to do all the columns at the same time. In this example I have only two columns, but I plan on having a dozen or so.


texts: Viewing relationship classes ties for multiple features

texts:
I've been working with a complex ArcGIS 9.3.1 geodatabase with multiple feature classes (FC's) interconnected by many relationship classes.

My question: supposing I have a FC with features a,b,c and another FC with features 1,2,3 - I want to see which featuers of the first FC are to connected to which  features of the second FC (of course in reality, the geodatabase is much more complex with more FC's, more features and more relations).

Now, I know that the IDENTIFY tool does it, but it does it to one feature at a time, and I can't export the ties to an outer file (it is only visual). Is there any tool I can use to do that for whole layers at a time?

(I've seen Geodatabase Diagrammer very briefly, but as far as I understood, it is useful for the planning of the structure of the geodatabase, but not it's content. correct me if I'm wrong).