1

I have a four band raster, and I would like to generate a new raster that has the argmax of the bands for each pixel. For example, if band 3 holds the maximum value at that pixel, the new raster should have a value 3 at that pixel.

I feel that this is a simple operation, but I am having a hard time with it. There is no argmax function in Postgres that I can find. And the callbackfunc syntax for ST_MapAlgebra() defeats me.

Is there a simple way to do this with an existing function?

Or are there examples of callback functions (along similar lines to this question, although that question wasn't really answered)?

Improve this question
  • Yes, I have really struggled with ST_MapAlgebra also. I am sure this has to be possible using a rasterbandarg and something like ST_Max4ma, but I haven't quite figured out how to control things so that the neighbourhood is 1x1 and the number of bands. It is possible to do this using ST_Value, and the select the max of an unnested array for each pixel, eg, – John Powell Sep 15 '17 at 8:55
  • Using dummy_rast from the docs, WITH pixels (x, y, ps) AS (SELECT x, y, array[ST_Value(rast, 1, x, y), ST_Value(rast, 2, x, y), ST_Value(rast, 3, x, y)] FROM dummy_rast CROSS JOIN generate_series(1, 5) As x CROSS JOIN generate_series(1, 5) As y WHERE rid = 2) SELECT min(p), x, y FROM (SELECT unnest(ps) AS p, x, y FROM pixels) foo GROUP BY x, y order by x, y; , but this is laughbably inefficient (and you still have to reconstruct the new raster from the return values). I will persist, because ST_MapAlgebra is so useful, but the callbacks are seriously hard to grok. – John Powell Sep 15 '17 at 8:56
  • 1
    I willl have another go at this later. I really need to understand these damn ST_MapAlgebra callbacks better, unless someone puts me out of my misery and answers the question. – John Powell Sep 15 '17 at 9:46
1

Here is an example that works with two bands. Up to you to generalize to four bands.

WITH raster AS (
  SELECT ST_AddBand(ST_AddBand(ST_AddBand(ST_SetValue(ST_AddBand(ST_MakeEmptyRaster(2, 2, 0, 0 ,1), '8BUI'), 1, 1, 2), '8BUI'::text, 2), '8BUI'::text, 2), '8BUI'::text, 4) rast
), bands AS (
  SELECT 1  band, ST_Band(rast, 1) rast FROM raster
  UNION ALL
  SELECT 2  band, ST_Band(rast, 2) rast FROM raster
)
, max AS (
  SELECT ST_Union(rast, 'MAX') rast FROM bands
), index1 AS (
  SELECT ST_MapAlgebra(b1.rast, max.rast, 'CASE WHEN [rast1] = [rast2] THEN 1 ELSE NULL END') rast
  FROM (SELECT rast FROM bands WHERE band = 1) b1, 
       max 
), index2 AS (
  SELECT ST_MapAlgebra(b2.rast, max.rast, 'CASE WHEN [rast1] = [rast2] THEN 2 ELSE NULL END') rast
  FROM (SELECT rast FROM bands WHERE band = 2) b2, 
       max
), indexbands AS (
  SELECT rast FROM index1
  UNION ALL
  SELECT rast FROm index2 
), indexmax AS (
  SELECT ST_Union(rast, 'MIN') rast FROM indexbands
) SELECT (ST_PixelAsPolygons(rast)).* FROM indexmax;

The idea is to create a temporary raster being the union(max) and to compare each band to this raster to build indexes rasters and then to union again these indexes rasters to find the first matching index.

Might be slow on big rasters given the number of operations...

I agree this is overwhelmingly complex for a task that seems easy. We would need a ST_Union() taking a user custom aggregate function. That's how ST_Union() was designed at the very beginning but it's not the way it was implemented. The Pl/pgSQL version still exists if you want it. Might need a little refresh.

Improve this answer
  • Thanks, @Pierre. This is a global 30 arc-second raster so it is insanely slow, but I don't see a better way to do it. – amball Sep 23 '17 at 17:26

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.