Now we are ready to start processing our dataset and make some scientific analysis.

The most common use cases of raster data are imagery, DEM, and statistical data. Right now, we'll look closer to Digital Elevation Model processing.

We can use the following:

• ST_Slope
• ST_Hillshade
• ST_Aspect
• ST_TPI
• ST_TRI
• ST_Roughness

All these functions are internally realized as MapAlgebra Callback functions. Later we will see how to do it on our way. So, calculate the slope values of our DEM data. The function syntax is like this:

ST_Slope(raster rast, integer nband=1, text pixeltype=32BF, text units=DEGREES, double precision scale=1.0, boolean interpolate_nodata=FALSE); 

We should especially pay attention to scale parameter. In case of geographical coordinate reference systems (such as our imported data), when units of distance are in degrees, and elevation in meters, use scale=111120, or in case of imperial units set it to 370400. Units parameter offers DEGREE, RADIANS, or PERCENT as a result of calculations. Because there are two variants of this function, arguments must be explicitly cast.

For calculation time reduction, I will use the clipped version of the dem table named clip. We will return to examples of extraction raster parts later in this chapter:

SELECT rid, ST_Slope(rast::raster, '1'::int,'32BF'::text,'DEGREES'::text,'111120'::double precision) AS rast FROM raster_ops.clip WHERE rid='4'; 

By analogue to slope calculations, aspect function is executed as follows, and the output raster will have aspect values in each pixel, in units specified in the fourth argument:

SELECT rid, ST_Aspect(rast::raster, '1'::int,'32BF'::text,'DEGREES'::text,'111120'::double precision) FROM raster_ops.clip WHERE rid='4'; 

More arguments must be set on ST_Hillshade functions.

ST_HillShade(raster rast, integer band=1, text pixeltype=32BF, double precision azimuth=315, double precision altitude=45, double precision max_bright=255, double precision scale=1.0, boolean interpolate_nodata=FALSE); 

We can execute a query with only raster column and band and it will work. But for LatLon datasets again, every argument must be explicitly set and cast to type.

Sometimes we need to change values of our dataset in an organized way, not in algebraic way, but as some groups-classes. There's a reclassify tool in PostGIS for that:

SELECT rid,ST_Reclass( 
(ST_SLOPE(rast::raster,'1'::int,'32BF'::text,'DEGREES'::text,'111120'::double precision))::raster, 
   '1'::int,   '0-5):1, 5-10:2, [10-15):3, 15-25):4, [25-90):5'::text,'32BF'::text) 
FROM eudem.clip; 

Pay attention to class borders and brackets. In this example, we are reclassifying the ST_slope function results to five classes of slope:

0-4.99 
5-9.99 
10-14.99 
15-24.99 
25-90 degrees 
with 1-5 values. 

For ST_TRI() and ST_TPI(), the case is more easier as the only argument for these functions is raster column.

TPI measures the relative topographic position of the central point as the difference between the elevation at this point and the mean elevation within a predetermined neighbourhood. Using TPI, landscapes can be classified in slope position classes. TPI is only one of a vast array of morphometric properties based on neighbouring areas that can be useful in topographic and DEM analysis (see Gallant and Wilson, 2000) Gallant, J.C., Wilson, J.P., 2000. Primary topographic attributes. In: Wilson, J.P., Gallant, J.C. (Eds.), Terrain Analysis: Principles and Applications. Wiley, New York, pp. 51-85.

Now, let's look into more complicated processing. As there's not much defined functions for DEM, sometimes we need to define it.

ST_MapAlgebra() is the tool that is needed.

There are two variants of this function expression and callbacfunc with one or two rasters.

The Expression version syntax is as follows:

raster ST_MapAlgebra(raster rast, integer nband, text pixeltype, text expression, double precision nodataval=NULL); 
raster ST_MapAlgebra(raster rast1, integer nband1, raster rast2, integer nband2, text expression, text pixeltype=NULL, text extenttype=INTERSECTION, text nodata1expr=NULL, text nodata2expr=NULL, double precision nodatanodataval=NULL); 

extenttype values allowed:

• INTERSECTION: The extent of the new raster is the intersection of the two rasters. This is the default.
• UNION: The extent of the new raster is the union of the two rasters.
• FIRST: The extent of the new raster is the same as the one of the first raster.
• SECOND: The extent of the new raster is the same as the one of the second raster.

These values are permitted for expression:

• [rast1] or [rast1.val]: This is the pixel value of the pixel of interest from rast1
• [rast1.x]: This is the 1-based pixel column of the pixel of interest from rast1
• [rast1.y]: This is the 1-based pixel row of the pixel of interest from rast1
• [rast2] or [rast2.val]:This is the pixel value of the pixel of interest from rast2
• [rast2.x]: This is the 1-based pixel column of the pixel of interest from rast2
• [rast2.y]: This is the 1-based pixel row of the pixel of interest from rast2

Our first expression is as follows:

CREATE TABLE ceiled AS SELECT ST_MapAlgebra(rast, 1, NULL, 'ceil([rast]*[rast.x]/[rast.y]+[rast.val])') FROM eudem.clip 

The result rendered by QGIS is as follows:

For callbackfunc version use trickier:

raster ST_MapAlgebra(raster rast1, integer nband1, raster rast2, integer nband2, regprocedure callbackfunc, text pixeltype=NULL, text extenttype=INTERSECTION, raster customextent=NULL, integer distancex=0, integer distancey=0, text[] VARIADIC userargs=NULL); 

Built-in callbackfunc:

• ST_Distinct4ma: This is the raster processing function that calculates the number of unique pixel values in a neighborhood
• ST_InvDistWeight4ma: This is the raster processing function that interpolates a pixel's value from the pixel's neighborhood
• ST_Max4ma: This is the raster processing function that calculates the maximum pixel value in a neighborhood
• ST_Mean4ma: This is the raster processing function that calculates the mean pixel value in a neighborhood
• ST_Min4ma: This is the raster processing function that calculates the minimum pixel value in a neighborhood
• ST_MinDist4ma: This is the raster processing function that returns the minimum distance (in number of pixels) between the pixel of interest and a neighboring pixel with value
• ST_Range4ma: This is the raster processing function that calculates the range of pixel values in a neighborhood
• ST_StdDev4ma: This is the raster processing function that calculates the standard deviation of pixel values in a neighborhood
• ST_Sum4ma: This is the raster processing function that calculates the sum of all pixel values in a neighborhood

As an example, use st_mean4ma:

SELECT    rid,    st_mapalgebra(rast, 1, 'st_mean4ma(double precision[][][],integer[][],text[])'::regprocedure, '32BF'::text, 'FIRST'::text, NULL) FROM raster_ops.clip 

We are not bound to these built in functions-we can prepare our own PL/pgSQL callback-just remember that it should have input arguments and should return the following:

 'sample_callbackfunc(double precision[], integer[], text[])'::regprocedure