1

I am trying to use PostGIS to create something similar to the attribute table produced when two raster files are combined using the ArcGIS Spatial Analyst function "Combine". This shows the count of unique combinations of values associated with each pixel across each raster.

I have a functioning query for a very small (8 pixel by 8 pixel) multiband raster which replicates this function where merged_rasters is my PostGIS raster table with 2 bands.

SELECT b1val, b2val, COUNT(*)
FROM(
    SELECT ST_Value(rast, 1, x, y) As b1val, ST_Value(rast, 2, x, y) As b2val
    FROM merged_rasters CROSS JOIN
    generate_series(1, 10) As x CROSS JOIN generate_series(1, 10) As y
    WHERE x <= ST_Width(rast) AND y <= ST_Height(rast)
    ) as merged_values
GROUP BY b1val, b2val;

However, when I try using the ST_Value on larger rasters, even tiled rasters with a tile size of 256x256, this becomes very slow.

How would I improve performance?

0

Update: Tested and verified!

Dumping the pixel values of those bands/rasters to be combined into a table using ST_DumpValues will be a lot more performant (one order of magnitude on average during my tests):

SELECT  ROW_NUMBER() OVER() AS "Value",
        b1val, b2val,
        COUNT(*) AS cnt
FROM    UNNEST(  -- parallel unnest of both arrays into a table, keeping corresponding pixels in one row
          (
            SELECT dmp.*
            FROM   merged_rasters,
                   LATERAL ST_DumpValues(rast, 1) AS dmp  -- create a table of values for each pixel of band 1
          ),
          (
            SELECT dmp.*
            FROM   merged_rasters,
                   LATERAL ST_DumpValues(rast, 2) AS dmp  -- create a table of values for each pixel of band 2
          )
        ) AS b(b1val, b2val)
WHERE   b1val IS NOT NULL
  AND   b2val IS NOT NULL  -- exclude NULL values; the dump excludes no data values by default, but UNNEST inserts NULL for missing column values to keep equal row count between arrays
GROUP BY
        b1val, b2val  -- grouping by (tuples of a) row, effectively getting distinct combinations
;

This will create the tabular output expected from the ArcGIS Combine tool.

4
  • When I try this on raster of, eg. 5000x5000 pixels it fails or just taked forever
    – Pemburu
    Nov 26 '20 at 10:36
  • @Pemburu yeah, you're attempting to dump 25 million compressed raster values into a virtual table twice, before grouping over both and aggregating without index or statistics; there is a reason why many GIS keep raster data in the file system with native OS access...
    – geozelot
    Nov 26 '20 at 10:55
  • Yes but the solution from above is actually not mine, it's that from Vince, I just tried it and it did not work for me, for reasons you explained correctly. Of course you are right, the native access is faster, but sometimes it also makes sense to integrate things in your DB workflows without depending on external software.
    – Pemburu
    Nov 27 '20 at 9:09
  • I was trying to find a solution based on gdal tools (see my link below [link] (gis.stackexchange.com/questions/377965/…) ) but without success. PostGIS just is such a wonderful tool for designing your workflows without requiring any external software packages. I'll always accept small or moderate performance drawbacks, if I get a good workflow in return. The solution I outlined below for me worked well, perhaps because one raster has a lot of NULL values. I still need to try this.
    – Pemburu
    Nov 27 '20 at 9:09
0

This corresponds to a question I just posted recently, without actually seeing your entry. Have a look here https://gis.stackexchange.com/questions/377965/combine-two-raster-with-postgis-or-gdal-calc.

When I try your solution on rasters of, eg. 5000x5000 pixels it fails or just takes forever.

I found a solution that works for me (one of the rasters has a lot of NULL pixels, though.

My solution looks like this:

 with hw_cells as (
    SELECT val::integer as hw,geom
        FROM 
        (
          SELECT dp.* 
          FROM  hw
             , LATERAL ST_PixelAsCentroids(hw.rast, 1) AS dp 
         ) foo
    )  
    ,
    stats_combine as 
    (
        select (split_part(combi,'_',1))::integer as hw
              ,(split_part(combi,'_',2))::integer as beam_id
              , n 
          from 
          (
    
                select   combi, count(*) as n 
                from 
                  (
                    select  hw_cells.hw::TEXT || '_' || 
                      st_value(beam.rast,hw_cells.geom)::TEXT as combi
                     from hw_cells, r_beam_echazerms_5 beam 
                     where st_intersects(beam.rast,hw_cells.geom)
                 
                   ) b
                   group by combi
            ) g
       )
    
select beam_id, hw, n    
from stats_combine

So, I basically create points from one raster and intersect them with the second raster. Works for me, although Combine tool in ArcGIS is still a lot faster. I admit it's a bit ugly in that I have to split the text. The query takes about 30 sec for a resulting table with ~ 16000 rows. Unlike the combine tool in ArcGIS it does not create a new raster.

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.