Creating Rasters from Scratch in PostGIS, Part 1

2017-03-23
#postgis #sql #raster #intro to postgis
database

Goal

In this post, we’ll learn how to create a single-band raster from scratch using PostGIS. This ability will be particularly useful as we begin exploring the relatively uncharted waters of the PostGIS Raster Module.

Purpose

Using ST_AddBand and ST_MakeEmptyRaster, we’ll create a single band, 8 bit unsigned integer raster. While it’s nothing fancy, it will help us better understand the mechanics of the Raster Module and it will also prove useful as a source of training data for the posts in the Raster Math Series.

Here, I’ll create a 250m resolution raster with a constant value. This method can easily create any other size and at whatever starting value the user may like. See the “All Together Now” section for some other common resolutions translated to decimal degrees.

Getting Started

First, lets create a table to hold the results of our operations.

1
2
3
4
-- Create a table for our new raster to live in
CREATE TABLE public.test_250(
  rid SERIAL primary key, rast raster
);

Creating a Raster

Now, let’s go about constructing our raster. I’m creating my raster starting somewhere in the Colorado Rockies. You don’t need to do the same place or size, it’s completely up to you.

You can adjust where and how big your raster is by changing the upper left coordinates and the number of pixels wide and tall your raster should be. The upper left coordinate is the origin, so no raster cells will appear “upper” or “lefter” of it.

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
-- Insert into our newly created table
-- Insert into that table
INSERT INTO public.test_250(rast)
SELECT ST_AddBand(
  -- Make empty raster
  ST_MakeEmptyRaster(
      1500, -- Raster width x (in pixels)
      1500, -- Raster width y (in pixels)
      -105.4330444336, -- Upper left X coordinate
      40.7170785158, -- Upper left Y coordinate
      0.00208333, -- X Cell Resolution (in degrees) (~250m)
      0.00208333, -- Y Cell Resolution (in degrees) (~250m)
      0, -- X skew
      0, -- Y skew
      4326 -- SRID (WGS 84)
    ),
    -- We're making a single band raster, but you can add
    -- as many bands as you like by adding additional rows
    -- to the array.
    ARRAY [
      ROW(
        1,  -- Band index: sets this as the first band
        '8BUI'::text,  -- Pixel Type (string rep of ST_BandPixelType types)
        5,  -- Initialized pixel value
        255  -- Nodata Value
      )
    ]::addbandarg[]
);

Et voilà! There you have it - a raster from scratch in PostGIS.

Now, If you did like me and left yourself with a massive raster, you may want to go ahead and tile it out to improve your query times. Additionally, this is a very common operation to perform on large rasers and it’s a good thing to be familiar with.

If you don’t want to tile your raster, you can go ahead and skip to the Adding Constraints section.

Tiling our New Raster

Tiling is a very common operation on large rasters. It allows you to build extremely fast spatial indexes on the data, vastly improving query times. For example, when I tried to load untiled large rasters into QGIS from PostGIS, it often crashes unless the data are tiled and indexed.

To tile your rasters, just do this:

1
2
3
4
5
6
7
8
9
-- Create a new table to hold our tiled raster
CREATE TABLE public.test_250_tiled(
  rid SERIAL primary key, rast raster
);
-- Tile it out
INSERT INTO public.test_250_tiled(rast)
  SELECT ST_Tile(rast, 1, 100, 100, TRUE, 0)
  FROM public.test_250;

Adding Constraints

There are a couple main reasons to add constraints to your raster tables after you’ve populated them. First, it improves query time by allowing the software to make certain assumptions about the geometry of your data. Second, it ensures that your data meet the standards PostGIS requires for various operations, such as alignment.

Add constraints by using the AddRasterConstraints function. See the documentation for more details on which constraints you can add and what they do.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
-- Add constraints explicitly
SELECT AddRasterConstraints(
  'test_250_tiled'::name,
  'rast'::name,
  'regular_blocking',
  'blocksize',
  'same_alignment',
  'srid'
);
-- Or, add all possible constraints
SELECT AddRasterConstraints(
  'test_250_tiled'::name,
  'rast'::name,
);

All Together Now

So, putting it all together, here’s what we’ve got:

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
45
46
47
48
49
50
51
-- Some common raster cell sizes:
-- 3 (90m) arcseconds = 0.000833333 degrees
-- 7.5 (250m) arcseconds = 0.00208333 degrees
-- 15 (500m) arcseconds = 0.00416667 degrees
-- Create a table for our new raster to live in
CREATE TABLE public.test_250(
  rid SERIAL primary key, rast raster
);
-- Insert into that table
INSERT INTO public.test_250(rast)
SELECT
  -- We're making a single band raster,
  -- though these can be chained together
  ST_AddBand(
    -- Make empty raster
    ST_MakeEmptyRaster(
      1500, -- Raster width x (in pixels)
      1500, -- Raster width y (in pixels)
      -105.4330444336, -- Upper left X coordinate
      40.7170785158, -- Upper left Y coordinate
      0.00208333, -- X Cell Resolution (in degrees) (~250m)
      0.00208333, -- Y Cell Resolution (in degrees) (~250m)
      0, -- X skew
      0, -- Y skew
      4326 -- SRID (WGS 84)
    ),
    '8BUI'::text -- Pixel Type (string rep of ST_BandPixelType types)
    ,5 -- Initialized pixel value
    ,0 -- Nodata Value
  ) AS rast;
-- Create a new table to hold our tiled raster
CREATE TABLE public.test_250_tiled(
  rid SERIAL primary key, rast raster
);
-- Tile it out
INSERT INTO public.test_250_tiled(rast)
  SELECT ST_Tile(rast, 1, 100, 100, TRUE, 0)
  FROM public.test_250;
SELECT AddRasterConstraints(
  'test_250_tiled'::name,
  'rast'::name,
  'regular_blocking',
  'blocksize',
  'same_alignment',
  'srid'
);

There you have it!

Enjoy your new raster powers! I hope this was helpful to you.