1

I'm currently using CDO to extract data from netCDF files containing satellite meteorological data. In the netCDF file points have a resolution of 0.5 degrees, so when data for a specific coordinate x,y is requested, I use cdo -remapbil,lon=x/lat=y to generate a temporary netCDF file from which to extract the data relative to x,y. Now as new data is being added the process is getting slower and slower.

In search of new ways to accomplish the task, I see that PostGIS can save raster data (in fact I've been able to import some data, converting netCDF to GeoTIFF then using raster2pgsql to import the data).

Now I'm stuck because I see that the raster is imported as pixels and I don't know how to formulate a query for a generic x,y coordinate that would return the same data as with CDO. How shall I tackle it?

Improve this question
2
  • I don't know CDO. What is the extent and the size of the temporary netCDF raster you produce with your command? 1x1 or more?
    – Pierre Racine
    Commented Aug 22, 2019 at 18:09
  • It's 1x1, since I need the value for a single x,y point.
    – marc
    Commented Aug 23, 2019 at 10:13

1 Answer 1

Reset to default
2 
SELECT ST_Value(rast, ST_SetSRID(ST_MakePoint(X, Y), ST_SRID(rast))) val
FROM rasterTable;

If you imported your raster without tiling it (with the -t option) or:

SELECT ST_Value(rast, ST_SetSRID(ST_MakePoint(X, Y), ST_SRID(rast))) val
FROM rasterTable
WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(X, Y), ST_SRID(rast)));

If you tiled your raster when importing it.

Improve this answer
2
  • 1
    Seems to be working fine thanks. So does PostGIS automagically work out which pixel in the raster to select based on the POINT() geometry returned by ST_MakePoint()? What if that point does not match one of the fixed-spaced points in the raster (that is, it would fall between points in the original netCDF)?
    – marc
    Commented Aug 23, 2019 at 12:02
  • 1
    Yes to the first question. All the space covered by the raster is covered by at least one pixel extent. So as long as the point intersects with the raster, it will intersect with a pixel. There is no empty space. So it more like 'nearest neighbour". If you want more sophisticated interpolation I guess you can use ST_Neighborhood() instead of ST_Value() and compute a value using your own algorythm.
    – Pierre Racine
    Commented Aug 23, 2019 at 17:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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