dian m fay

Let's say you have two rasters, and you want to combine them with some extra value math -- perhaps you want to grade grassland around 50°N 100°E by latitude and elevation, so blue to green to red the more northerly or higher up the point. ST_MapAlgebra to the rescue!

select
  c.rid,
  st_mapalgebra(
    e.rast,
    c.rast,
    '
      case when [rast2.val] in (13, 14, 16, 17, 18)
        then -[rast1.y] + [rast1.val]
        else 0
      end
    '
  ) as rast
from diva.coverage as c
join diva.elevation as e on e.rid = c.rid;

Oops. (effect exaggerated for visibility)

If they are tiled it's a bit more complicated...

— Pierre Racine

Problem: instead of two big aligned rasters, you have a bunch of tiled aligned rasters, so instead of a smooth gradient south to north you have lots of little smooth south-north gradients. The rest of it works, but only within each tile, or here band since we're ignoring x.

How to fix this? The expression argument only has a handful of available values, and all of them are internal to a specific raster value, that is, the tile. There's a function-callback version of ST_MapAlgebra that could add conditional logic based on external factors like tile latitude; however, I don't want to write and maintain a whole function for this really rather straightforward calculation.

But! ST_MapAlgebra executes per raster, and that means the expression argument is passed into the invocation each time. This means we can use format to pass in external variables -- here, converting pixel y-value to latitude within the SRID. The y-value of any given pixel relative to the SRID as a whole is given by subtracting its y-value within the tile from the latitude of the tile's top edge (ST_UpperLeftY) divided by the height a pixel represents in the SRID (ST_PixelHeight). Mileage may vary in the southern hemisphere, but this too is tractable.

select
  c.rid,
  st_mapalgebra(
    e.rast,
    c.rast,
    format(
      '
        case when [rast2.val] in (13, 14, 16, 17, 18)
          then %1$s - ([rast1.y] * %2$s) + [rast1.val]
          else 0
        end
      ',
      st_upperlefty(c.rast),
      st_pixelheight(c.rast)
    )
  ) as rast
from diva.coverage as c
join diva.elevation as e on e.rid = c.rid;