A Totally Random Map Algebra Callback in PostGIS, Part 2

Goal

In this post we’ll try to get a better understanding of the structure of plpgsql ST_MapAlgebra callback functions. Specifically, we’ll learn how to write a callback that works on one raster with many bands.

This is the second post in this series, and I encourage you to go back and read the first if you haven’t already. It gives a much more detailed account of what the parameters mean and why I’m doing things the way I am.

Purpose

It’s pretty common to operate on more than one band of a raster when performing MapAlgebra calculations. If you’re doing any serious image processing, it’s a fact of life. While the last example was bordering on overkill, we approach the limits of the Expression Version of ST_MapAlgebra very quickly.

Here we’ll provide another fairly basic example of a plpgsql callback function for ST_MapAlgebra that elucidates how one can interact with the data. This example works one-cell-at-a-time, but you can also operate on neighboring cells using this same callback format (post coming soon).

Getting Started

If you’re working with real data, carry on! If you’re worried about damaging your dataset or simply don’t want to find one, try simulating some data like I do in my post on creating rasters from scratch..

We’ll start off with our same basic structure as last time, but notice that we declare 2 variables that will hold our pixel values.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

-- The basic structure of an ST_MapAlgebra() callback
CREATE OR REPLACE FUNCTION
totally_random_2band_callback(pixel double precision[][][], pos integer[][], VARIADIC userargs text[])
  RETURNS integer
  LANGUAGE plpgsql
  IMMUTABLE  -- careful: this function is immutable, yours may not be
  AS $$
    DECLARE
      r1_val integer;
      r2_val integer;
    BEGIN
      -- Do Something
    END;
  $$;

Of course, you don’t really need to declare variables to hold your pixel values. They just help with clarity for this example.

Let’s “Do Something”

For this example, I’ll use a dummy raster I created whose first band is initialized with all values set to 5. The second band of the raster contains random integers from 1 to 3 (inclusive). Our callback will generate a random integer between 1 and 10 (inclusive), multiply it by the first band, and divide the product by the second band.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

-- The basic structure of an ST_MapAlgebra() callback
CREATE OR REPLACE FUNCTION
totally_random_2band_callback(pixel double precision[][][], pos integer[][], VARIADIC userargs text[])
  RETURNS integer
  LANGUAGE plpgsql
  IMMUTABLE  -- careful: this function is immutable, yours may not be
  AS $$
    DECLARE
      r1_val integer;
      r2_val integer;
      rando integer;
      result integer;
    BEGIN
      SELECT trunc(random() * 10 + 1)::int INTO rando;  -- Get a random value between 1 & 11, convert to int
      r1_val := value[1][1][1]::int;
      r2_val := value[2][1][1]::int;
      result :=  (r1_val * rando) / r2_val;
      RETURN result;
    END;
  $$;

Run that.

Get FUNCTION CREATED back.

Let’s play ball.

We’ll call the function and insert the result into an empty table we’ve already created.

1
2
3
4
5
6
7
8
9
10

-- Inster into the new table
INSERT INTO public.two_band_random(rast)
-- Call MapAlgebra
SELECT ST_MapAlgebra(
  ARRAY[
    ROW(rast, 1),
    ROW(rast, 2)
  ]::rastbandarg[],  -- Notice that this is pretty different from last time
  'totally_random_callback(double precision[][][], integer[][], text[])'::regprocedure)
FROM public.sim_250;

Let’s use QGIS and the DBManager plugin to check our results.

Looks pretty good!