0

PostgreSQL: 12.4 (running on Amazon Aurora)
PostGIS: 3.0.2

I'm attempting to clip a much larger raster (global meteorological data on a 0.25 degree resolution grid, split into 64x64 tiles, with 9 bands for different altitudes) down to a subset around the point I'm interested in.

Beyond that, for context, I want to obtain a Binlinearly interpolated point from the raster. To do that, I've been using ST_SnapToGrid to manipulate the pixel centroids to align with my desired point and interpolate the data. My intention with the clipping is that rather than manipulating the entire raster, I'm only processing the part that's relevant, which I hope will improve performance.

Data has been ingested from a TIFF file using raster2pgsql and piped to pgsql:

raster2pgsql -s 4236 -I -a -M data_20200513_12z_t06.tiff -F -t 64x64 public.data_table -Y | psql --dbname $DB_NAME --host $DB_SERVER --username $DB_USERNAME -q

I'm running this query to get my subset of the original raster:

WITH myquery AS (
    SELECT 
        ST_GeomFromText('POINT (148.1484375 -89.630)', 4236) as geom,
        'data_20200513_12z_t06.tiff' as filename
),
subset_raster AS (
    SELECT
        ST_Clip(rast, ST_Buffer(geom, ST_ScaleX(rast)), ST_BandNoDataValue(rast)) as clipped_rast, data_table.filename, rast as original_rast, geom as query_geom
    FROM
        data_table
    INNER JOIN 
        myquery ON myquery.filename=data_table.filename     
    WHERE
        ST_Intersects(myquery.geom, rast)
)
SELECT (ST_DumpValues(clipped_rast)).* FROM subset_raster

The clip seems to work, as it creates a new 3x3 pixel raster based on the expected coordinates SELECT (ST_MetaData(clipped_rast)).* FROM subset_raster:

147.875 -89.375 3   3   0.25    -0.25   0   0   4236    9

But there's no data in the raster SELECT (ST_DumpValues(clipped_rast)).* FROM subset_raster:

1   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
2   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
3   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
4   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
5   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
6   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
7   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
8   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}
9   {{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}

So what am I missing? Is this the intended behaviour?
How should I query the data if I want to constrain the raster to minimise unnecessary processing?

0

It seems that my issue was related to the point I was querying (POINT (148.1484375 -89.630)) which is very close to the South Pole.

Querying a point further away from the pole produced results more consistent with what I was expecting.

What's odd is that the point definitely fell within a pixel that had data, but it didn't appear in the clipped result (and neither did its neighbours).

For now I've put it down to "polar oddness" in lieu of a more definitive explanation

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.