+ multi-res pt density lyr: H3 hexagons, pt vector tiles

[bbest]

• Tile Serving with Dynamic Geometry

  CREATE OR REPLACE
  FUNCTION tilehexagons(z integer, x integer, y integer, step integer,
                        OUT geom geometry(Polygon, 3857), OUT i integer, OUT j integer)
  RETURNS SETOF record
  AS $$
      DECLARE
          bounds geometry;
          maxbounds geometry := ST_TileEnvelope(0, 0, 0);
          edge float8;
      BEGIN
      bounds := ST_TileEnvelope(z, x, y);
      edge := (ST_XMax(bounds) - ST_XMin(bounds)) / pow(2, step);
      FOR geom, i, j IN
      SELECT ST_SetSRID(hexagon(h.i, h.j, edge), 3857), h.i, h.j
      FROM hexagoncoordinates(bounds, edge) h
      LOOP
         IF maxbounds ~ geom AND bounds && geom THEN
              RETURN NEXT;
           END IF;
       END LOOP;
       END;
  $$
  LANGUAGE 'plpgsql'
  IMMUTABLE
  STRICT
  PARALLEL SAFE;

  see animation of zoom in/out
  

• Waiting for PostGIS 3.1: Grid Generators · Paul Ramsey

  SELECT (ST_HexagonGrid(100000, ST_Transform(a.geom, 3857))).* 
  FROM admin a  
  WHERE name = 'Germany';

  

[bbest]

Woohoo, got it working! 🎉

tile.calcofi.io/public.hexagons.html

DONE:

• Upgraded pgadmin.calcofi.io to v 8.0 to show geometry columns in a map

TODO:

• add all variables (depth, time, etc) as arguments into function with joins as needed

Created spatial index on ctd_casts.geom

CREATE INDEX ctd_casts_geom_idx ON ctd_casts USING GIST (geom);

Create hexagon() function

CREATE 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;

Create hexagoncoordinates() function

CREATE 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;

Prep ctd_casts.geom_3857 column

For 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);

Create visible function 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.';

[bbest]

Next steps:

• create app to view with hexagons colored by n_ctd_casts
• use H3 hexagons instead
• index every point with H3 hexagons (all in one column possible, or need to individuate per H3 scale?)

[bbest]

Protoype of app for coloring hexagons by n_ctd_casts:

• shiny.calcofi.io/hex
  

Issues:

• don't know max value for given view

[bbest]

ctd casts as point vector tiles

This vector tile layer of CTD casts as points through our PostGIS tile server was revelatory for @evsatt .

tile.calcofi.io: ctd_casts