Skip to content

Commit

Permalink
update examples
Browse files Browse the repository at this point in the history
  • Loading branch information
vincentsarago committed Nov 17, 2022
1 parent e81e371 commit 6e0b627
Show file tree
Hide file tree
Showing 6 changed files with 275 additions and 2 deletions.
140 changes: 140 additions & 0 deletions data/functions/hexagon.sql
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
-- Custom "hexagon" layer
-- Adapted from https://github.com/CrunchyData/pg_tileserv
--
--
-- Given an input ZXY tile coordinate, output a set of hexagons
-- (and hexagon coordinates) in web mercator that cover that tile
CREATE OR REPLACE FUNCTION tilehexagons(
bounds geometry,
step integer,
epsg integer,
OUT geom geometry(Polygon, 3857), OUT i integer, OUT j integer)
RETURNS SETOF record AS $$
DECLARE
maxbounds geometry := ST_TileEnvelope(0, 0, 0);
edge float8;
BEGIN
edge := (ST_XMax(bounds) - ST_XMin(bounds)) / pow(2, step);
FOR geom, i, j IN
SELECT ST_SetSRID(hexagon(h.i, h.j, edge), epsg), 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;

-- Given coordinates in the hexagon tiling that has this
-- edge size, return the built-out hexagon
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;

-- Given a square bounds, find all the hexagonal cells
-- of a hex tiling (determined by edge size)
-- that might cover that square (slightly over-determined)
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;

-- Given an input tile, generate the covering hexagons Step parameter determines
-- how many hexagons to generate per tile.
CREATE OR REPLACE FUNCTION hexagon(
-- mandatory parameters
xmin float,
ymin float,
xmax float,
ymax float,
epsg integer,
-- additional parameters
query_params json
)
RETURNS bytea AS $$
DECLARE
result bytea;
step integer;
bounds geometry;
BEGIN
-- Find the bbox bounds
bounds := ST_MakeEnvelope(xmin, ymin, xmax, ymax, epsg);

step := coalesce((query_params ->> 'step')::int, 4);

WITH
rows AS (
-- All the hexes that interact with this tile
SELECT h.i, h.j, h.geom
FROM TileHexagons(bounds, step, epsg) h
),
mvt AS (
-- Usual tile processing, ST_AsMVTGeom simplifies, quantizes,
-- and clips to tile boundary
SELECT ST_AsMVTGeom(rows.geom, bounds) AS geom, rows.i, rows.j
FROM rows
)
-- Generate MVT encoding of final input record
SELECT ST_AsMVT(mvt, 'default')
INTO result
FROM mvt;

RETURN result;
END;
$$
LANGUAGE 'plpgsql'
STABLE
STRICT
PARALLEL SAFE;

50 changes: 50 additions & 0 deletions data/functions/landsat_poly_centroid.sql
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
CREATE OR REPLACE FUNCTION landsat_poly_centroid(
-- mandatory parameters
xmin float,
ymin float,
xmax float,
ymax float,
epsg integer,
-- additional parameters
query_params json
)
RETURNS bytea
AS $$
DECLARE
bounds geometry;
tablename text;
result bytea;
BEGIN
WITH
-- Create bbox enveloppe in given EPSG
bounds AS (
SELECT ST_MakeEnvelope(xmin, ymin, xmax, ymax, epsg) AS geom
),
selected_geom AS (
SELECT t.*
FROM public.landsat_wrs t, bounds
WHERE ST_Intersects(t.geom, ST_Transform(bounds.geom, 4326))
),
mvtgeom AS (
SELECT
ST_AsMVTGeom(ST_Transform(ST_Centroid(t.geom), epsg), bounds.geom) AS geom, t.path, t.row
FROM selected_geom t, bounds
UNION
SELECT ST_AsMVTGeom(ST_Transform(t.geom, epsg), bounds.geom) AS geom, t.path, t.row
FROM selected_geom t, bounds
)
SELECT ST_AsMVT(mvtgeom.*, 'default')

-- Put the query result into the result variale.
INTO result FROM mvtgeom;

-- Return the answer
RETURN result;
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE -- Same inputs always give same outputs
STRICT -- Null input gets null output
PARALLEL SAFE;

COMMENT ON FUNCTION landsat_poly_centroid IS 'Return Combined Polygon/Centroid geometries from landsat table.';
65 changes: 65 additions & 0 deletions data/functions/squares.sql
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
CREATE OR REPLACE FUNCTION squares(
-- mandatory parameters
xmin float,
ymin float,
xmax float,
ymax float,
epsg integer,
-- additional parameters
query_params json
)
RETURNS bytea AS $$
DECLARE
result bytea;
sq_width float;
bbox_xmin float;
bbox_ymin float;
bounds geometry;
depth integer;
BEGIN
-- Find the bbox bounds
bounds := ST_MakeEnvelope(xmin, ymin, xmax, ymax, epsg);

-- Find the bottom corner of the bounds
bbox_xmin := ST_XMin(bounds);
bbox_ymin := ST_YMin(bounds);

-- We want bbox divided up into depth*depth squares per bbox,
-- so what is the width of a square?
depth := coalesce((query_params ->> 'depth')::int, 2);

sq_width := (ST_XMax(bounds) - ST_XMin(bounds)) / depth;

WITH mvtgeom AS (
SELECT
-- Fill in the bbox with all the squares
ST_AsMVTGeom(
ST_SetSRID(
ST_MakeEnvelope(
bbox_xmin + sq_width * (a - 1),
bbox_ymin + sq_width * (b - 1),
bbox_xmin + sq_width * a,
bbox_ymin + sq_width * b
),
epsg
),
bounds
)

-- Drive the square generator with a two-dimensional
-- generate_series setup
FROM generate_series(1, depth) a, generate_series(1, depth) b
)
SELECT ST_AsMVT(mvtgeom.*, 'default')

-- Put the query result into the result variale.
INTO result FROM mvtgeom;

-- Return the answer
RETURN result;
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE -- Same inputs always give same outputs
STRICT -- Null input gets null output
PARALLEL SAFE;
2 changes: 1 addition & 1 deletion demo/leaflet/index_3035.html
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
}
);

var mapboxUrl = "http://0.0.0.0:8081/tiles/EuropeanETRS89_LAEAQuad/public.countries/{z}/{x}/{y}";
var mapboxUrl = "http://127.0.0.1:8000/tiles/EuropeanETRS89_LAEAQuad/public.countries/{z}/{x}/{y}";

var mapboxVectorTileOptions = {
rendererFactory: L.canvas.tile,
Expand Down
2 changes: 1 addition & 1 deletion demo/leaflet/index_4326.html
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
continuousWorld: true,
}).addTo(map)

var mapboxUrl = "http://0.0.0.0:8081/tiles/WGS1984Quad/public.countries/{z}/{x}/{y}.pbf";
var mapboxUrl = "http://127.0.0.1:8000/tiles/WGS1984Quad/public.countries/{z}/{x}/{y}";

var mapboxVectorTileOptions = {
rendererFactory: L.canvas.tile,
Expand Down
18 changes: 18 additions & 0 deletions docker-compose.yml
Original file line number Diff line number Diff line change
Expand Up @@ -37,3 +37,21 @@ services:
command: postgres -N 500
volumes:
- ./.pgdata:/var/lib/postgresql/data

# pg_tileserv:
# image: pramsey/pg_tileserv:latest
# environment:
# - DATABASE_URL=postgresql://username:password@database:5432/postgis
# ports:
# - "7800:7800"
# depends_on:
# - database

# martin:
# image: maplibre/martin
# environment:
# - DATABASE_URL=postgresql://username:password@database:5432/postgis
# ports:
# - "3000:3000"
# depends_on:
# - database

0 comments on commit 6e0b627

Please sign in to comment.