A Totally Random Map Algebra Callback in PostGIS

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 one band.

Purpose

There are a few posts about how to do this. This one from George MacKerron’s blog gives a much better example of how MA callbacks work in PostGIS, and serves as a model of what I’ll attempt here. As Mr. MacKerron points out, the callback function in the documentation is pretty useless, but it’s also totally uninformative.

Basically, we’ll write a totally useless function which illustrates how we can manipulate raster values using ST_MapAlgebra.

Getting Started

If you’re working with some real data, great! If you’re trying to learn more about ST_MapAlgebra but don’t have a dataset, go check out my post on creating rasters from scratch.. This type of stuff is exactly what I had in mind when I wrote it.

One thing worth noting, is that there are two versions of ST_MapAlgebra: the Expression Version and the Callback Version. Obviously, we’ll be working with the latter in this post. Both are great to be familiar with, but the Callback Version let’s us do much more complicated stuff, including working with cell neighbors (post coming soon). Here however, we’ll be doing a pretty basic operation. Still, it’s one that you couldn’t perform with the Expression Version (as far as I know, anyway)!

Using MacKerron’s post as a guide, I’ve laid out a basic structure that our callback has to follow.

We need three arguments:

1. A 3 dimensional array of a type double precision pixel double precision[][][]
2. A 2 dimensional integer array pos integer[][]
3. A 1 dimensional variadic array to accept any additional user arguments variadic userargs text[]

In our case, we’re working with integer data so we’ll plan our RETURNS statement accordingly:

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

-- The basic structure of an ST_MapAlgebra() callback
CREATE OR REPLACE FUNCTION
totally_random_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
      pixval integer;
    BEGIN
      -- Do Something
    END;
  $$;

A note of caution:

I’ve seen several tutorials, including MacKerron’s and this one from Regina Obe, that use an array of a different type for the pixel value. I tried this, and it never worked for me with integer. So, I’m sticking to double precision here even though we’ll convert back to integer at the end (as you can see from RETURNS integer).

Let’s “Do Something”

We’ve initialized our raster with a value of 10 for all cells. As you can see from this picture, it’s pretty boring to look at. Let’s throw some chaos in there!

Since we are creating a random raster, but not entirely random. After all, it would be nice to learn a little bit about how to interact with our raster data since that’s the primary goal of this post. Bearing that in mind, let’s perform our operations using the values of our initial raster.

Let’s just get a random integer between 1 and 10 (inclusive) for each cell/pixel. We’ll multiply this by pixval, our initial raster value.

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

CREATE OR REPLACE FUNCTION
totally_random_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
      pixval integer;
      rando integer;
      result integer;
    BEGIN
      SELECT floor(random() * 9 + 1)::int INTO rando;  -- Get a random value between 1 & 11, convert to int
      pixval := value[1][1][1]::int                    -- pixel indices: [raster #][xdistance][ydistance]
      result :=  pixval * rando;                       -- Notice we cast the pixel value back to integer
      RETURN result;
    END;
  $$;

Okay, so we run that, and we get FUNCTION CREATED back. Time to test it out!

Here we’ll call the function on a dummy raster we saw above. We’ll take the result and stuff it in the empty table we created for the purpose.

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

-- Create a table for our new raster to live in
CREATE TABLE public.totally_random(
  rid SERIAL primary key, rast raster
);
-- Inster into the new table
INSERT INTO public.totally_random(rast)
-- Call MapAlgebra
SELECT ST_MapAlgebra(
  rast,  -- Raster data column
  1,  -- Which band to use (1-based index)
  'totally_random_callback(double precision[][][], integer[][], text[])'::regprocedure)
  -- ^^ Our callback function in a string cast as a regprocedure ^^
  -- Notice that we don't use the names of the inputs here
FROM public.test_250;  -- The table containing our initial raster

Popping on over to QGIS, we can check out our results using the DBManager plugin. You can also use whatever tools you like, but I think those ones are great and they work very well together.

Wow, look at all that uniform noise! It seems like we achieved our objective. For a look at how to do write a callback that works with multiple bands, check out the next post in this series.