I am new to working with rasters in PostGIS. Although I have spent couple of days trying to find solution for my problem, I couldn't seem to crack it without some help.

Here's my problem:
I have multiple overlapping rasters in a PostGIS table (about 500 records) with each raster size about 1200x1000 (30m res.) but varies since data is dumped from different sources. The table looks as follows:

-- Table "coverage_rast"

   id   |   rast
(bigint)| (raster)
   5489 | "02301ff..."
   5489 | "238000f..."
   5489 | "02301ff..."
   1234 | "238000f..."
   1234 | "02301ff..."
   4567 | "238000f..."
   4567 | "238000f..."

Each raster record represents a coverage with some DN values(ranging 0-100). When I perform CREATE TABLE sometable as SELECT ST_UNION(rast,1,'MIN') from coverage_rast; gives me resultant raster with minimum values as DN values at each grid cell and takes about 2 minutes of query time. Instead I want that the resultant raster should be such that the DN value at each grid cell should represent "id" of the raster which has minimum value at that location.

I tried using SELECT val,id FROM (SELECT ST_Value(rast, 1, geom) val,id FROM coverage_rast,(SELECT st_setsrid(ST_MakePoint(-80.85495, 35.21887),4326) geom)tbl WHERE ST_Intersects(rast, geom)) as t1 ORDER BY val LIMIT 1 which gets the id and the minimum DN value at given location. Was thinking of extending this by creating a uniform grid of points spaced 30m apart and passing geom column and then creating a raster out of it, but this would be pretty slow and sloppy.

I am hoping to find a solution which can extend st_union to assign id value as DN of raster where existing DN value is minimum of all the overlapping rasters.

  • Combine ID and DN values into a single 32-bit value, then use ST_MapAlgebra (which runs only on a single band). – CL. Mar 21 at 11:51
  • @CL. following your comment I tried to combine ID and DN values by using create table sometable as select ST_MapAlgebra(rast, 1, NULL, '[rast]*"id"') from coverage_rast but it gives column "id" not found error. I understand that column names are not permitted in MapAlgebra expression. I looked into st_reclass but it seems that reclassexpr also does not allow column names. So is there some other way to do this? – Parthesh Bulbule Mar 22 at 14:36
  • That's a complex case. You will need to develop your own ST_UnionID() aggregate function using ST_MapAlgebra() as the state function, a mix of raster value and ID as the state value and your own callback function returning the desired mix of those values (MIN or MAX). – Pierre Racine Mar 22 at 16:08
  • The combination step would have to be done by hand (outside of PostGIS, or by handling individual pixels). It might be a better idea to not use PG at all. – CL. Mar 22 at 18:57

Your Answer

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

Browse other questions tagged or ask your own question.