[postgis-users] PostGIS Raster: Two Questions

[Regina Obe]

*  First question: how can one recode gbr100x100d so that all pixels with value 0 are NODATA?

Use ST_SetBandNoDataValue https://postgis.net/docs/manual-3.1/RT_ST_SetBandNoDataValue.html

UPDATE gr100x100d SET ST_SetBandNoDataValue(rast,1,0);

 

General side note, using ST_Reclass for what you are doing below should be much faster than ST_MapAlgebra

 

*  Second question:  I haven’t had a chance to look that closely.  What is the purpose of the ST_SnapToGrid?

 

 

From: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] On Behalf Of Simon SPDBA Greener
Sent: Tuesday, March 16, 2021 8:15 PM
To: postgis-users at lists.osgeo.org
Subject: [postgis-users] PostGIS Raster: Two Questions

 

I have a 1 band raster, srid 4326, (0.001 x 0.001 roughly 100m x 100m) , called gbr100x100d, whose pixels are coded 0 and 1 where 1 is water deeper than 100m and 0 water shallower than 100m. 

It was created as follows:

create table gbr100x100d
as
SELECT row_number() over (order by 1) as rid,
       ST_SetSrid(ST_MapAlgebra(a.rast, 1, NULL, '([rast.val]<-100.0)::integer'),4326) as rast
  from gbr100x100c as a;

First question: how can one recode gbr100x100d so that all pixels with value 0 are NODATA?

Next. 

Polygons of size 0.02x 0.01 need to be extracted from the 1 valued 0.001 pixels.

I think this should be done as follows.

First we resample to get a raster with cells of size 0.02 x 0.01:

create table gbr2kmx1kmd1
as
SELECT row_number() over (order by 1) as rid,
       ST_SetSrid(ST_Rescale(a.rast,0.02,-0.01,'Bilinear',0.125),4326) as rast
  from gbr100x100d as a;

Then we extract the polygons of these 0.02 x 0.01 pixels.

create table gbr2kmx1kmd1p
as 
SELECT row_number() over (order by 1) as id, (gv).val, ST_SnapToGrid((gv).geom,0.001)::geography as geog
  FROM (SELECT ST_PixelAsPolygons(a.rast,1,true) as gv 
          FROM gbr2kmx1kmd1 as a
       ) gv
 WHERE (gv).val = 1;

The result is close as the polygons are 0.02 x 0.01 but they are organised as a set of striped polygons where the columns don't touch. 

See attached.

What was expected was a complete coverage of abutting polygons. 

Can anyone spot what I am doing wrong? 

regards

Simon


2021-03-17_10-27-06.png 



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210316/e98df005/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 18942 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210316/e98df005/attachment.png>