Create a raster with fixed 10meter pixel scale in PostGIS

Question

I have a series of points that I would like to convert to a raster, and then calculate the area of the raster as a proxy for "amount of area covered". I'm trying to do this using PostGIS rather than another language so that I can share my query easily with collaborators.

One annoyance is that my points are stored as json (not geojson).

id | points
---+--------
 0 | [{'lat': 33.18424, 'lon': -100.0428, 'property': 'Group A'}, {'lat': 33.18424, 'lon': -100.0428, 'property': 'Group A'}, ...
 1 | [{'lat': 33.18424, 'lon': -100.0423, 'property': 'Group A'}, {'lat': 33.18424, 'lon': -100.0458, 'property': 'Group A'}, ...
 2 | [{'lat': 33.18436, 'lon': -100.0422, 'property': 'Group B'}, {'lat': 33.18420, 'lon': -100.0428, 'property': 'Group B'}, ...

So I have to do a little bit of work to get my points in a workable format.

SELECT ST_SetSRID(ST_Point((latlon->'lon')::text::numeric, (latlon->'lat')::text::numeric), 4326)::geometry as geom
FROM (SELECT JSON_ARRAY_ELEMENTS(points::json) as latlon
    FROM my_table) as my_table

This works fine, but now I am having some trouble getting ST_AsRaster to work for me. I want the pixels to represent a 10mx10m area however there are so many variants of ST_AsRaster I am not sure how to specify this correctly. Here is my full query:

WITH my_geom AS (SELECT ST_SetSRID(ST_Point((latlon->'lon')::text::numeric, (latlon->'lat')::text::numeric), 4326)::geometry as geom
    FROM (SELECT JSON_ARRAY_ELEMENTS(points::json) as latlon
        FROM my_table) as my_table)
SELECT ST_AsRaster(geom, scalex=10.0, scaley=10.0) as my_raster
FROM my_geom

I think I may need to transform the point geometry SRID or otherwise somehow specify units (meters) for ST_AsRaster but I am not sure this is possible.

Answer 1

I've found the best way to create a raster feature is to use a template raster. The template raster should have the same extent as the union of the features you want to rasterise, and each pixel the desired dimensions.

I'll copy some code from a project I recently worked on, where I converted very large, irregular vector features into raster objects.

First, I created a mask_raster (a raster template with pixels of 25x25 metres, and an extent equal to the extent of the union of all of my inputs). This makes for large rasters in many cases, but it is easier if you then want to do some calculations on overlapping rasters. The syntax is:

ST_MakeEmptyRaster(integer width, integer height, float8 upperleftx, float8 upperlefty, float8 scalex, float8 scaley, float8 skewx, float8 skewy, integer srid=unknown);

And my code was:

CREATE TABLE mask_raster AS (
SELECT ST_MakeEmptyRaster(5861, 4141, 1735076, 5390595, 25, 25, 0, 0, 2193) AS rast
);
UPDATE mask_raster
SET rast = ST_AddBand(rast,'1BB'::text,0);
--Set constraints
SELECT AddRasterConstraints('mask_raster'::name, 'rast'::name);

You'll need to read the documentation for these functions to determine appropriate parameters, such as the pixel type. I used 1BB because I was making simple Boolean surfaces of presence and absence, rather than holding quantitative values.

Later, I used this template as an input to the ST_AsRaster function. The template raster thereby gives its extent and pixel dimensions to your conversion of vector features. My command was:

WITH ra AS (
    WITH mb AS (
        SELECT isogeom,
        (SELECT rast FROM mask_raster) AS rast --Holds a null raster built from envelope of all MBs
        --This controls the resolution, origin, size of the output raster
        FROM mb_isochrones
        WHERE source = 'RTI'
        AND mb2013 = '%s'
        AND time = '07:15:00'::time
        AND date >= '%s'::date --start date for loop
        AND date < '%s'::date --end date for loop
        ORDER BY date
    )
    SELECT
    rast,
    ST_AsRaster(isogeom, rast, ARRAY['1BB'], ARRAY[1], ARRAY[0], true) AS rast2
    FROM mb)
SELECT
ST_MetaData(ra.rast2),
ST_MetaData(ra.rast),
ST_SameAlignment(ra.rast, ra.rast) as sm,
(ST_DumpValues(ra.rast2,ARRAY[1],false)).*
FROM ra;

I presented the whole query to give some context. (Using psycopg2, I was also then dumping the raster features to arrays to work with them in NumPy... PostGIS raster is surprisingly difficult). The important part of the above query for your purposes are:

(SELECT rast FROM mask_raster) AS rast

and

ST_AsRaster(isogeom, rast, ARRAY['1BB'], ARRAY[1], ARRAY[0], true) AS rast2

In the line just above, the ST_AsRaster function uses rast, the mask_raster as a template. That's the easiest way to handle this, in my experience. Failure to use a mask raster layer can lead to misalignment problems if you're making many raster features.