Creating Rasters from Scratch in PostGIS, Part 3

Goal

In this post, we’ll learn how to create a multi-band raster from an existing raster in PostGIS.

Purpose

This technique is useful for more than the simulation I discussed in the first and second posts in this series. It allows us to combine the data of many rasters into a single raster, which can make the data a little easier to work with.

In this post, I’ll demonstrate how I created the initial raster I used in my post on writing multi-band ST_MapAlgebra callback functions.

Getting Started

What I wanted when I set out was a raster made up of two 8-bit Unsigned Integer bands: the first with all cell values set to 5 and the second with values all ranging between 1 and 5 roughly based on a normal distribution.

In practice, we rarely modify existing raster tables when doing analysis or data management (but you can). We typically create new ones. Bearing this in mind, it didn’t make sense to begin with a two band raster. Especially considering that we can’t initialize a raster band with an existing value. We have to modify the cells with an ST_MapAlgebra callback, and ST_MapAlgebra only returns single band rasters.

Using Map Algebra

So, we begin with a single band raster:

1
2
3
4
5
6
7
8
9
10
11

-- Create a new table for our new raster
CREATE TABLE public.starter(
  rid SERIAL primary key, rast raster
);
-- Insert into our newly created table
INSERT INTO public.starter(rast)
SELECT ST_AddBand(
  ST_MakeEmptyRaster(1500,1500,-105.4330444336,40.7170785158,0.00208333,0.00208333,0,0,4326),
  1, '8BUI'::text, 5, 0);
-- Want to learn more about this query and its parameters?
-- Check out the first post in this series.

With our single band raster created, let’s see what we can do about adding another. As I said above, we’ll need an ST_MapAlgebra callback function to do the job. Here’s what that function looks like:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

CREATE OR REPLACE FUNCTION
generate_random_int_raster(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL::text[])
  RETURNS double precision
  LANGUAGE plpgsql
  IMMUTABLE  -- careful: this function is immutable, yours may not be
  AS $$
    DECLARE
      rando integer;
      non_zero boolean;
      mean integer;
      sd integer;
    BEGIN
      non_zero := TRUE;    -- Used in the conditional below
      mean := 3;           -- Mean of my normal distribution
      sd := 1;             -- Standard deviation of my normal distribution
      -- I've set up this conditional statement so that we can
      -- control whether or not we want to return 0 as a value.
      -- In my case, 0 is my nodata value, so I don't want that
      -- to happen.
      IF non_zero THEN     -- Make sure that our random value is greater than 0
      rando := 0;          -- Initialize our random value at 0
      LOOP
        IF rando > 0 THEN  -- If rando is non-zero...
          EXIT;            -- exit the loop.
        ELSE               -- Otherwise,
          rando := trunc(normal_rand(1, mean, sd))::integer;  -- get a new value.
        END IF;
      END LOOP;
                           -- If 0 is allowed, just get a value.
      ELSE
      rando := trunc(normal_rand(1, mean, sd))::integer;
      END IF;
      RETURN rando;        -- Notice that we only return the random number,
                           -- and don't do anything with the existing band.
    END;
    $$;

When we created our starter raster, we added a band to the empty raster in the torast position of the ST_AddBand function (see docs, Varaiant 5).

This time, we’ll specify starter as the torast (instead of an empty raster) and we’ll specify the results of our ST_MapAlgebra callback function as the fromrast. We also stipulate that this new band should be band two, although this isn’t strictly necessary as omitting the band value will automatically append the new band to the end of the raster’s band list.

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

CREATE TABLE public.sim_250(
  rid SERIAL primary key, rast raster
);
INSERT INTO public.sim_250(rast)
SELECT
  ST_AddBand(b1.rast,  -- torast
    ST_MapAlgebra(b1.rast,  -- Starter raster
      1,  -- First band of our starter raster is the target band (not that we use it)
      'generate_random_int_raster(double precision[][][], integer[][], text[])'::regprocedure
    ),  -- fromrast
    1,  -- Target FROM band (band one of fromrast)
    2   -- Target TO band (band two of torast)
  )
AS rast
FROM starter AS b1;

Just like that, you’ve got a new two band raster! This technique is quite useful if you want to perform an operation on a few bands of a raster, then store the result in the same raster in a new band.

Using Existing Rasters

Of course, you can do the same sort of thing using only existing rasters. Perhaps you need to combine several rasters about climate into a single one to use in a model. You can do so like this:

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

CREATE TABLE public.new_raster(
  rid SERIAL primary key, rast raster
);
SELCT INTO public.new_raster(rast)
SELECT ST_AddBand(
          ST_AddBand(b1.rast, b2.rast, 2, 2), -- Set band 2 of our second raster
                                              -- as band 2 of the first raster
            b3.rast, 1, 3) -- Set band 1 of the third raster
                           -- as band 3 of the first raster
AS rast
FROM test_b1 AS b1, test_b2 AS b2, test_b3 AS b3;

I hope this post was useful to you. To learn more about ST_MapAlgebra and its callbacks, make sure to check out my Raster Math Series.