Volgens de auteur zou de volgende functie een raster moeten maken met de omvang gebaseerd op de gegeven BBOX en de celgrootte in meters .
SRID 3857-eenheden zijn [ongeveer] meter, en door deze projectie te gebruiken, worden hexcellen gemaakt die er "goed uitzien" op een webkaart (de meeste die een webmercatorprojectie gebruiken).
CREATE OR REPLACE FUNCTION generate_hexgrid(width float, xmin float, ymin float, xmax float, ymax float, srid int default 3857)
RETURNS TABLE(gid text,geom geometry(Polygon)) AS $$
DECLARE
b float := width / 2;
a float := tan(radians(30)) * b;
c float := 2 * a;
height float := 2 * (a + c);
index_xmin int := floor(xmin / width);
index_ymin int := floor(ymin / height);
index_xmax int := ceil(xmax / width);
index_ymax int := ceil(ymax / height);
snap_xmin float := index_xmin * width;
snap_ymin float := index_ymin * height;
snap_xmax float := index_xmax * width;
snap_ymax float := index_ymax * height;
ncol int := abs(index_xmax - index_xmin);
nrow int := abs(index_ymax - index_ymin);
polygon_string varchar :=
'POLYGON((' || 0 || ' ' || 0 || ' , ' || b || ' ' || a || ' , ' ||
b || ' ' || a + c || ' , ' || 0 || ' ' || a + c + a || ' , ' ||
-1 * b || ' ' || a + c || ' , ' || -1 * b || ' ' || a || ' , ' ||
0 || ' ' || 0 ||'))';
BEGIN
RETURN QUERY
SELECT
format('%s %s %s', width,
x_offset + (1 * x_series + index_xmin),
y_offset + (2 * y_series + index_ymin)),
ST_SetSRID(ST_Translate(two_hex.geom,
x_series * width + snap_xmin,
y_series * height + snap_ymin), srid)
FROM generate_series(0, ncol, 1) AS x_series,
generate_series(0, nrow, 1) AS y_series,
(SELECT 0 AS x_offset, 0 AS y_offset, polygon_string::geometry AS geom
UNION
SELECT 0 AS x_offset, 1 AS y_offset, ST_Translate(polygon_string::geometry, b , a + c) AS geom
) AS two_hex;
END; $$ LANGUAGE plpgsql;
Aangezien je een tabel hebt met de naam usa
en het bevat de geometrieën van dit shapefile
je zou een raster moeten kunnen maken dat de kaart van de VS overlapt met de volgende vraag:
CREATE TABLE usa_hex AS
WITH j AS (
SELECT ST_Transform((hex).geom,4326) AS hex FROM (
SELECT
generate_hexgrid(
80467,
ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
FROM usa)i)
SELECT j.hex FROM j,usa
WHERE ST_Intersects(usa.geom,j.hex);
BEWERKEN :Het is nog steeds geen antwoord, omdat het geen zeshoeken maakt met meters, maar het kan andere gebruikers helpen. De volgende functie (afgeleid van dit answer
) maakt geometrietype zeshoeken van exact dezelfde grootte in graden .
CREATE OR REPLACE FUNCTION generate_hexagons(width FLOAT, bbox BOX2D, srid INTEGER DEFAULT 4326)
RETURNS TABLE (gid INTEGER, hexagon GEOMETRY) AS $$
DECLARE
b FLOAT := width/2;
a FLOAT := b/2;
c FLOAT := 2*a;
height FLOAT := 2*a+c;
ncol FLOAT := ceil(abs(ST_Xmax(bbox)-ST_Xmin(bbox))/width);
nrow FLOAT := ceil(abs(ST_Ymax(bbox)-ST_Ymin(bbox))/height);
polygon_string VARCHAR := 'POLYGON((' ||
0 || ' ' || 0 || ' , ' || b || ' ' || a || ' , ' || b || ' ' || a+c || ' , ' || 0 || ' ' || a+c+a || ' , ' ||
-1*b || ' ' || a+c || ' , ' || -1*b || ' ' || a || ' , ' || 0 || ' ' || 0 || '))';
BEGIN
RETURN QUERY
SELECT
row_number() OVER ()::INTEGER,
ST_SetSRID(
ST_Translate(geom, x_series*(2*a+c)+ST_Xmin(bbox), y_series*(2*(c+a))+ST_Ymin(bbox)),srid)
FROM generate_series(0, ncol::INTEGER, 1) AS x_series,
generate_series(0, nrow::INTEGER,1 ) AS y_series,
(SELECT polygon_string::GEOMETRY AS geom
UNION
SELECT ST_Translate(polygon_string::GEOMETRY, b, a + c) AS geom) AS two_hex;
END;
$$ LANGUAGE plpgsql;
Overlapt met de hierboven gebruikte dataset:
WITH j (hex_rec) AS (
SELECT generate_hexagons(3.0,ST_Extent(geom))
FROM usa
)
SELECT (hex_rec).gid,(hex_rec).hexagon FROM j, usa
WHERE ST_Intersects(usa.geom,(hex_rec).hexagon);
Verder lezen:
ST_Extent
ST_Intersects