著者によると、次の関数は、指定されたBBOXとメートルのセルサイズに基づいた範囲のグリッドを作成する必要があります。 。
SRID 3857ユニットは[非常に]メートルであり、この投影法を使用すると、Webマップ上で「正しく見える」六角形のセルが作成されます(そのほとんどはWebメルカトル図法を使用します)。
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;
usa
というテーブルがあるとします。 このshapefile
のジオメトリが含まれています 次のクエリを使用して、USAマップとオーバーラップするグリッドを作成できるはずです。
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);
編集 :メーターを使用して六角形を作成しないため、まだ答えではありませんが、他のユーザーには役立つ可能性があります。次の関数(このanswer
から派生) )ジオメトリタイプを作成します 度でまったく同じサイズの六角形 。
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;
上記で使用したデータセットとの重複:
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);
さらに読む:
-
ST_Extent
-
ST_Intersects