Making a raster coverage that shows the average speed of points

Question

I had a layer of points of trajectory and I wanted to create a map of the average speed of the movement. I have a column of speed in the table of points. I created a big base raster coverage. Then, I used ST_PixelAsPolygons() to create a polygon of each pixel. Then, I calculated the average speed of points in each pixel. Finally, I used ST_AsRaster() to create a raster coverage of the speed map.

I used the following code (where passages is the point table and grid_pixel is the polygon of each pixel and sog(speed over ground) is the speed):

CREATE TABLE speed_map AS(

with speed_average AS (
SELECT g.id , AVG(p.sog) AS sog_avg
FROM grid_pixel g
INNER JOIN passages p
ON ST_Within(p.geom_point, g.geom)
GROUP BY g.id
),

speed_avg_geom AS (
SELECT g.id AS id, g.geom AS geom, sa.sog_avg 
FROM speed_average sa
INNER JOIN grid_pixel g
ON g.id = sa.id)


SELECT id AS id , ST_AsRaster(geom, 1, 1, '8BUI', sog_avg) AS rast
FROM speed_avg_geom
)

Now, each small polygon of grid_pixel table is transfered to a raster tile (one row of raster table) and I have a very big table that has one pixel in each row. The map is still correct and I can see the correct raster in QGIS. However, I think it is better to have tiles of 100*100 pixels.

I wanted to retile the raster and I found ST_reTile() function but I do not know how to use it to retile the raster to 100*100 rows.

Or maybe even my method is totally wrong from the beginning where I created a polygon for each pixel. Can anyone please let me know what is the best way to create a raster of average speed of points.

Answer 1

First advice would be to create a ST_MapAlgebra() custom callback function instead of converting everything to polygon and then back to raster.

ST_MapAlgebra() will iterate over all the raster pixels assigning the value your function computed from the vector coverage. You can extract any king of metric with that. The function has to return only one value compatible with the pixel type of your raster.

Have a look at the example in the documentation.

Another exmple is ST_ExtractPixelValue4ma() in the PostGIS Addons which implements many extraction methods from a vector coverage. ST_ExtractToRaster() calls ST_ExtractPixelValue4ma() with all the pixels properties necessary to compute an intersection between the extent of each pixel and a vector coverage (or anything).

If you continue the way you started, you cav recombine 1x1 rasters into 100x100 rasters with ST_Union(). Something along the line;

SELECT ST_Union(rast) rast
FROM speed_map 
GROUP BY round(ST_UpperLeftX(rast), 2), round(ST_UpperLeftY(rast), 2)