Skip to content

Instantly share code, notes, and snippets.

@bitner
Created February 12, 2020 18:43
Show Gist options
  • Save bitner/1e5017997b4c3dbb6d04fdb595c0bec6 to your computer and use it in GitHub Desktop.
Save bitner/1e5017997b4c3dbb6d04fdb595c0bec6 to your computer and use it in GitHub Desktop.
Create a grid to aggregate phenomena by.
-- 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