I have used postgis for vector operations for years, but have recently begun a project mixing vector and raster data and I am bit confused as to what is the best way to do what I want to do. I have the trajectory of a boat as lon,lat points that I want to rasterize, counting the number of points in each grid cell. The trick is that I want to groups the points according to time, counting a set of points only once if they all occur in the same grid cell within the same time window (e.g., within an hour of each other). Efficiency is important because in reality I have many many such trajectories each with a large number of points.
I can think of several ways of attacking this problem, a couple of which I detail below, but in all cases there is a step that I am not sure how to do. More generally, I am wondering if there is a standard approach to these types of problems.
- Group the points of a trajectory by time using ST_Union grouping over a truncated version of the timestamp for each point
- Use ST_AsRaster on each of these groups of points to produce one raster layer for each group with 1's wherever the points intersect the reference raster layer.
- Sum together these individual raster layers to get the final result.
My problems with this approach are (1) that I am not sure how to sum all the rasters together and (2) it seems inefficient to me to create all these raster layers. Regarding the first problem, I believe I can use ST_SummaryStatsAgg to do this (though I would need to update my version of postgis to have this function), but an example would help. Regarding the second problem, I will have many many time periods and trajectories and calculating all of these rasters and then computing a full suite of summary statistics over these rasters seems daunting.
This approach essentially passes through vector calculations only returning to raster maps at the end.
- Convert the raster layer to a vector layer using ST_PixelAsPolygons (adding a spatial index to these polygons)
- Intersect these polygons with the points to associate a grid cell with each point.
- Group the points by time period and grid cell so that we have a table with one row per time period and intersected grid cell combination.
- Aggregate across time periods so that we count the number of "positive" (boat present) time periods per grid cell.
- Convert these counts per grid cell back into a raster layer.
The problem here is that I am not sure how to do the last step. I can image using ST_AsRaster on the centroid of each grid cell to generate a raster with one non-zero pixel, then using ST_MapAlgebraExpr to multiply this 1-pixel raster by the counts, an then using ST_SummaryStatsAgg to put all these 1-pixel rasters back together, but there must be a better way.