Created
February 12, 2020 18:43
-
-
Save bitner/1e5017997b4c3dbb6d04fdb595c0bec6 to your computer and use it in GitHub Desktop.
Create a grid to aggregate phenomena by.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
-- Procedural Function that takes in an extent and a cellsize (in extent srid units) | |
-- And creates a grid bounding that extent with given cellsize with an origin that is in | |
-- a multiple of the cellsize + zero | |
CREATE OR REPLACE FUNCTION ST_GRID( | |
IN _extent geometry, | |
IN _cellsize float8, | |
OUT id bigint, | |
OUT r integer, | |
OUT c integer, | |
OUT geom geometry | |
) | |
RETURNS SETOF record AS | |
$$ | |
WITH e AS ( | |
SELECT st_snaptogrid(st_expand(_extent,_cellsize), _cellsize) as extent | |
) | |
, dims AS ( | |
SELECT | |
st_xmin(extent) as xmin, | |
st_ymin(extent) as ymin, | |
st_xmax(extent) - st_xmin(extent) as width, | |
st_ymax(extent) - st_ymin(extent) as height | |
FROM e | |
), rowcol AS ( | |
SELECT | |
xmin, | |
ymin, | |
floor(width / _cellsize)::int as ncol, | |
floor(height / _cellsize)::int as nrow | |
FROM dims | |
) | |
SELECT | |
row_number() over () as id, | |
i + 1 AS row, | |
j + 1 AS col, | |
st_setsrid( | |
ST_Translate( | |
cell, | |
j * _cellsize + xmin, | |
i * _cellsize + ymin | |
), | |
st_srid(_extent) | |
) AS geom | |
FROM | |
rowcol, | |
generate_series(0, rowcol.nrow - 1) AS i, | |
generate_series(0, rowcol.ncol - 1) AS j, | |
st_makeenvelope(0,0,_cellsize, _cellsize) AS cell | |
$$ LANGUAGE sql IMMUTABLE STRICT; | |
-- Create an index on the bbox in phenomena using a transform to the srid that the | |
-- grid will be in | |
-- the st_contains is there since we can't index data that falls outside of the bounds | |
-- of the 3857 projection | |
CREATE INDEX ON phenomena USING GIST (st_transform(bbox, 3857)) WHERE | |
st_contains(st_makeenvelope(-180,-85, 180, 85, 4326), bbox); | |
-- Query to create an aggregated grid with grid created on the fly | |
SELECT g.id, count(*), geom | |
FROM | |
phenomena, | |
ST_GRID( | |
st_transform( | |
st_makeenvelope(0,0,15,15,4326)::geometry, -- envelope in epsg:4326, in practice would be the view bounds | |
3857 | |
), | |
100000 -- grid size in srid units | |
) g | |
WHERE | |
st_contains(st_makeenvelope(-180,-85, 180, 85, 4326), bbox) -- needed right now to filter bogus data | |
AND | |
st_intersects(st_transform(bbox,3857), g.geom) -- the magic sauce | |
GROUP BY g.id, geom; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment