How to add a band with conditional value to a raster with PostGIS

Question

I'm trying to convert a vector polygon into a raster, with pixel values that indicate the area of the intersection. I have a query that gets the area of intersection between an empty raster and the polygon:

SELECT ST_Area(ST_Transform(geometry, 26986)) / 1000000 as area, x, y
FROM (
    SELECT ST_Intersection(ST_Transform(geom, 4326), wkb_geometry) as geometry,
    x, y
    FROM (
        SELECT (ST_PixelAsPolygons(rast)).*
        FROM empty_raster
    ) raster_pixels,
    vector_polygon
    WHERE region_name = 'Calgary-Buffalo'
    AND ST_Intersects(ST_Transform(geom,4326), wkb_geometry)
) intersects

Which gives me the output:

  | area  | x  |  y  |
1 | 4.998 | 79 | 203 |
2 | 1.535 | 80 | 203 | 

The empty_raster table is an empty raster with width of 140 and height of 248.

What I need to do is to add a band that sets values to 0 except for the above x-y values.

I've thought to try something like this:

SELECT ST_AddBand(rast, '32BF'::text, area)
FROM empty_raster, see_above_query

But I don't know what my WHERE clause would be.

---

I have found that I should create the band with 0 values first:

UPDATE empty_raster
SET rast = ST_AddBand(rast, 1, '32BF'::text, 0)

And update the values with ST_SetValues or ST_SetValue.

Now, when I try to update the values, it only updates the first value:

UPDATE empty_raster
SET rast = ST_SetValue(rast, 1, x, y, area)
FROM see_first_query

Am I using this function incorrectly? Do I need to do a for loop to generate the proper SQL?

Answer 1

I found that ST_SetValues can work with a collection of geometries and values as an aggregated array. So my SQL solution is this:

UPDATE empty_raster
SET rast = ST_SetValues(rast, 1, geomval)
FROM (
    SELECT array_agg((geom, (ST_Area(ST_Transform(intersected_geometry, 3400)) / 1000000))::geomval) as geomval
    FROM (
        SELECT geom, ST_Intersection(ST_Transform(geom, 4326), wkb_geometry) as intersected_geometry
        FROM (
            SELECT (ST_PixelAsPolygons(rast)).geom
            FROM empty_raster
        ) raster_pixels,
        vector_polygon
        WHERE region_name = 'Calgary-Buffalo'
        AND ST_Intersects(ST_Transform(geom,4326), wkb_geometry)
    ) intersects
) values;

Starting from the inner SELECT statement, it gets the polygon geometries from the raster, intersects them with the vector polygon (where they intersect), and puts the raster polygon geometry in an aggregated array with the area of the intersected geometry in km2. Those values are passed to ST_SetValues in band 1.