Raster analysis in a web application, using PostGIS Raster

A few weeks ago, I posted some techniques for using Python and agoodle to analyze raster data,clipping to a polygon and then reporting min/max/mean and pixel counts. Another technique available as of a few months ago, is PostGIS Raster.

PostGIS 2.0 brings a lot of small changes and improvements, among them support for processing raster datasets. Yes, all the wonderful GEOS goodness of PostGIS... can now be done on rasters too.

Importing The Data

PostGIS can in fact import the binary content of your raster file into the database, as a set of binary blobs. We don't do this; we keep the raster in the filesystem (the -R flag) so we can download it, run other command-line programs, etc. Our use of raster2pgsql looks like this:

raster2pgsql -s 3174 -t 1000x1000 -I -C -F -R /projects/maine/rasters/landvalue_2005.tif

You can read the help in raster2pgsql to see the details about each flag, but in summary:

• The table is created with a reference to a raster in the filesystem, and the SRS is stated.
• The table will have multiple rows, each one being a 1000x1000 slice of the raster.
• The column rast refers to the raster content.

Querying A Slice of The Data

So, we have a polygon such as county boundary, and classified raster data such as NLCD land cover codes, and we want to find the pixel values and pixel counts within that county.

First, a test of the raster and pixel counting:

SELECT ST_ValueCount(rast,1,false)
AS valuecount FROM nlcd2005;

This performs a ST_ValueCount() on the rast column, without any other sort of filtering, clipping, etc. This returns back a set of rows, each row being an array of 2 elements: the pixel value and the count of pixels having that value.

 (0,23975242)
 (1,442029)
 (2,265327)
 (3,42985)

Let's do this again, but this time clipping to a polygon:

SELECT ST_ValueCount(
ST_Clip(rast,1, 
    ST_Transform(ST_GeomFromText('POLYGON((521879 5011246,522182 5011140,521751 5010750,521751 5010966,521358 5010966,521879 5011246))',26910), 3174)
    ,1,false
)
AS valuecount FROM nlcd2005;

All this did, was replace "rast", the first parameter to ST_ValueCount(), with a call to ST_Clip() which clips a raster to a given geometry, then returns the resulting raster. Simple.

Processing The Value Counts

A lot of the raster functions return a set of rows, and these are a nuisance to work with. You could tease them apart like this:

SELECT (ST_ValueCount(rast,1,false)).value AS pixel, 
    (ST_ValueCount(rast,1,false)).count AS howmany
    FROM nlcd2005;

But this has two problems: a) It runs the value count twice, which takes twice as long; and b) those field names within the record are poorly documented at this time, so you may take some guessing or reading the source code to find them.

Eventually I gave up on this, and took the easy way out: regular expressions. Given that it's a set of numbers inside parens, it's pretty simple to just process it within PHP where we would be handling these row results anyway:

$counts = array();
$bits = null;
preg_match('/\((\d+)\,(\d+)\)/', $pxcount->valuecount, $bits );
$pixel = $bits[1];
$count = $bits[2];
$counts[$pixel] += $count;

And now you have an associative array of pixel values onto pixel counts. From here, follow your usual techniques to figure out acres, percentages, or whatever you had in mind.

Pixel Counts To Acres

The easiest way to find acres, is to multiply the pixel count by the known size in meters/feet, then by the conversion factor to get acres. Like this:

$acres_nlcd_11 = $counts['11'] * 30 * 30 * 0.000247105;

If you don't know the pixel size of this specific raster, PostGIS Raster can tell you: ST_PixelWidth() and ST_PixelHeight() It won't tell you the units, though.