Open bbest opened 10 months ago
Woohoo, got it working! 🎉
tile.calcofi.io/public.hexagons.html
DONE:
TODO:
ctd_casts.geom
CREATE INDEX ctd_casts_geom_idx ON ctd_casts USING GIST (geom);
hexagon()
functionCREATE OR REPLACE
FUNCTION
hexagon(i integer, j integer, edge float8)
RETURNS geometry
AS $$
DECLARE
h float8 := edge*cos(pi()/6.0);
cx float8 := 1.5*i*edge;
cy float8 := h*(2*j+abs(i%2));
BEGIN
RETURN ST_MakePolygon(ST_MakeLine(ARRAY[
ST_MakePoint(cx - 1.0*edge, cy + 0),
ST_MakePoint(cx - 0.5*edge, cy + -1*h),
ST_MakePoint(cx + 0.5*edge, cy + -1*h),
ST_MakePoint(cx + 1.0*edge, cy + 0),
ST_MakePoint(cx + 0.5*edge, cy + h),
ST_MakePoint(cx - 0.5*edge, cy + h),
ST_MakePoint(cx - 1.0*edge, cy + 0)
]));
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE
STRICT
PARALLEL SAFE;
hexagoncoordinates()
functionCREATE OR REPLACE
FUNCTION hexagoncoordinates(bounds geometry, edge float8,
OUT i integer, OUT j integer)
RETURNS SETOF record
AS $$
DECLARE
h float8 := edge*cos(pi()/6);
mini integer := floor(st_xmin(bounds) / (1.5*edge));
minj integer := floor(st_ymin(bounds) / (2*h));
maxi integer := ceil(st_xmax(bounds) / (1.5*edge));
maxj integer := ceil(st_ymax(bounds) / (2*h));
BEGIN
FOR i, j IN
SELECT a, b
FROM generate_series(mini, maxi) a,
generate_series(minj, maxj) b
LOOP
RETURN NEXT;
END LOOP;
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE
STRICT
PARALLEL SAFE;
ctd_casts.geom_3857
columnFor faster querying in same projection
-- Create index
CREATE INDEX ctd_casts_geom_3857_idx ON ctd_casts USING GIST (geom_3857);
-- Update geom_3857
UPDATE ctd_casts
SET geom_3857 = ST_Transform(geom, 3857);
-- Add a spatial column to the table
SELECT AddGeometryColumn ('public','ctd_casts','geom_3857', 3857, 'POINT', 2);
hexagons()
-- Create visible function at https://tile.calcofi.io
CREATE OR REPLACE
FUNCTION public.hexagons(z integer, x integer, y integer, step integer default 4)
RETURNS bytea
AS $$
WITH
bounds AS (
-- Convert tile coordinates to web mercator tile bounds
SELECT ST_TileEnvelope(z, x, y) AS geom
),
hex AS (
-- All the hexes that interact with this tile
SELECT h.i, h.j, h.geom
FROM TileHexagons(z, x, y, step) h
),
counts AS (
-- count number of ctd_cast points per hexagon
SELECT h.i, h.j, count(c.geom) AS n_ctd_casts
FROM hex h
LEFT JOIN ctd_casts c ON st_intersects(h.geom, c.geom_3857)
GROUP BY h.i, h.j
),
hex_counts AS (
SELECT h.i, h.j, c.n_ctd_casts, h.geom
FROM hex h
LEFT JOIN counts c USING (i, j)
),
mvt AS (
-- Usual tile processing, ST_AsMVTGeom simplifies, quantizes,
-- and clips to tile boundary
SELECT ST_AsMVTGeom(hc.geom, b.geom) AS geom,
hc.i, hc.j, hc.n_ctd_casts
FROM hex_counts hc, bounds b
)
-- Generate MVT encoding of final input record
SELECT ST_AsMVT(mvt, 'public.hexagons') FROM mvt
$$
LANGUAGE 'sql'
STABLE
STRICT
PARALLEL SAFE;
COMMENT ON FUNCTION public.hexagons IS 'Hex coverage of the number of ctd_casts dynamically generated. Step parameter determines how approximately many hexes (2^step) to generate per tile.';
Next steps:
Protoype of app for coloring hexagons by n_ctd_casts
:
Issues:
This vector tile layer of CTD casts as points through our PostGIS tile server was revelatory for @evsatt .
Tile Serving with Dynamic Geometry
see animation of zoom in/out
Waiting for PostGIS 3.1: Grid Generators · Paul Ramsey