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 '20 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 '20 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 '20 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 '20 at 18:57

With the help of CL's comments I combined the ID and DN values using psycopg2 such that e.g id=5489 and DN=144.5601 then the DN of new raster becomes 144564.5489. So I take float value with 2 decimal precision and in units place I keep the number of digits in the id(since id could be 4 or more digits long), and in the decimal places I keep the id value. So when I run ST_Union on the new raster I would get the ID also embedded in the DN. Then I run mapalgebra to use only the ID values and construct a new raster. Following is how I did it:

import psycopg2

tbl = "coverage_rast"
conn = psycopg2.connect("dbname="+dbname+" port="+port+" host="+host+" password="+password+" user="+user)
cur = conn.cursor()

def create_ID_DN_combo(tbl,conn,cur):
    cur.execute("ALTER TABLE "+tbl+" ADD COLUMN no_digits integer; update "+tbl+" set no_digits = (floor(log(id))+1);")
    cur.execute("Select id,no_digits from "+tbl+";")   # no_digits is no. of digits in id
    cols = cur.fetchall()
    for i in range(len(cols)):
        cur.execute("update "+tbl+" set rast = st_mapalgebra(rast,1,'64BF','(floor([rast]*100)*10)+"+str(cols[i][1])+"+"+str(cols[i][0]/(10**len(str(cols[i][0]))))+"') where id="+str(cols[i][0])+";")
    cur.execute("create table combo_tbl as select st_union(rast,1,'MIN')rast from "+tbl+";")
    cur.execute("update combo_tbl set rast = st_mapalgebra(rast,1,'32BUI','round(([rast]-floor([rast]))*(10^(((floor([rast])/10)-floor(floor([rast])/10))*10)))');")


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.