Since introduction of resampling on ST_Value (v3.2) behaviour on inexisting band has changed

[lieven]

Hello,

Regarding the changes made since version 3.2 of PostGIS, the behavior of ST_Value with geometry points now utilizes "RASTER_getPixelValueResample". However, there have been alterations in the behavior when dealing with non-existing bands.

Here's my query:

WITH poi AS (
    SELECT ST_SetSRID(ST_Point(1.75, 46.37), 4326) AS coord
)
SELECT
    rast.time,
    ST_Value(rast, 1, coord) AS u_wind,
    ST_Value(rast, 2, coord) AS v_wind,
    ST_Value(rast, 3, coord) AS temperature,
    ST_Value(rast, 4, coord) AS humidity,
    ST_Value(rast, 5, coord) AS pressure,
    ST_Value(rast, 6, coord) AS cloud_cover,
    ST_Value(rast, 7, coord) AS precipitations,
    ST_Value(rast, 8, coord) AS ww,
    rast.model
FROM poi
LEFT JOIN weather.inspire_2023_12_raster_4326 rast ON (ST_Intersects(rast, coord))
WHERE time = TIMESTAMP WITHOUT TIME ZONE '2023-12-20 00:00:00';

Regarding the bands 6, 7, etc., potentially being empty, there's a change in behavior from prior versions where the function used to return the first band's value and set null for others. However, since version 3.2, an error message occurs and I have no result at all:

"ERROR: Could not find raster band of index 6 when getting pixel value. Returning NULL SQL state: XX000"

This change in behavior seems to be associated with the new method "RASTER_getPixelValueResample" utilized with geometry, whereas using raster coordinates with version > 3.2 employs "RASTER_getPixelValue", maintaining the previous version's functionality without exceptions.

The main issue here is that the behavior changes depending on the parameters With this constraint I'm currently stuck with version < 3.2

In all cases, thank you for your amazing job with PostGis

[lieven]

Working scenario with result, even if extra band doesn't not exist

[lieven]

Hi, Have you any alternative solution for this scenario ? I wish to upgrade my postgis version but I'm stuck with this problem. If you have any clue to help me, it will be amazing

Thanks a lot

[robe]

Replying to lieven:

Hi, Have you any alternative solution for this scenario ? I wish to upgrade my postgis version but I'm stuck with this problem. If you have any clue to help me, it will be amazing

PS: it seems to be also a problem with postgresql version. Working on 14 not on 15 version

Thanks a lot

It's postgresql version specific? T Are you testing against the same PostGIS version? If you are seeing a difference between PostgreSQL 14 and 15 using the same version of postgis, then it's most likely a GDAL version difference. Anyrate I'm seeing an issue with my latest dev instance running:

POSTGIS="3.5.0dev 3.4.0rc1-886-g03d6a1d77" [EXTENSION] PGSQL="160" GEOS="3.13.0dev-CAPI-1.18.0" PROJ="9.3.1 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/share/proj/proj.db" GDAL="GDAL 3.8.1, released 2023/11/28" LIBXML="2.9.14" LIBJSON="0.17" LIBPROTOBUF="1.4.1" WAGYU="0.5.0 (Internal)" RASTER 

Please post the

SELECT postgis_full_version();

Of the versions working as you would expect and the versions not working as you would expect.

At anyrate if this is a postgis change, it was unintentional, so I'll see if we can resolve on our end..

Regarding your query, an alternative solution that will work in both cases (this version does rely on postgresql returning null if array element item does not exist.

You could also alternatively use a CASE WHEN ST_NumBands(rast) ≥ 7 THEN ST_Value(rast, 7) ELSE NULL END but from experience I think the CASE statement might slow it down more than below

WITH poi AS (
    SELECT ST_SetSRID(ST_Point(1.75, 46.37), 4326) AS coord
)
SELECT
    rast.time,
    v.elements[1] AS u_wind,
    v.elements[2] AS v_wind,
    v.elements[3] AS temperature,
    v.elements[4] AS humidity,
    v.elements[5] AS pressure,
    v.elements[6] AS cloud_cover,
    v.elements[7] AS precipitations,
    v.elements[8] AS ww,
    r.model
FROM poi
LEFT JOIN weather.inspire_2023_12_raster_4326 r ON (ST_Intersects(r.rast, poi.coord))
    LEFT JOIN LATERAL 
      (SELECT array_agg(ST_Value(r.rast, i, poi.coord) ORDER BY i) AS elements 
          FROM generate_series(1, ST_NumBands(r.rast)) AS i ) AS v ON true
WHERE time = TIMESTAMP WITHOUT TIME ZONE '2023-12-20 00:00:00';

[Paul Ramsey <pramsey@…>]

In 595a67e/git:

Raise a notice, not an error, when unable to read a band, and return NULL, references #5649

[Paul Ramsey <pramsey@…>]

In 28d30004/git:

ST_Value should return NULL on missing band, references #5649

[Paul Ramsey <pramsey@…>]

In 373326f/git:

ST_Value should return NULL on missing band, references #5649