ST_HexagonGrid 几何向量查找所有点

Question

我正在 PostGis 中查看此功能

https://postgis.net/docs/manual-dev/ST_HexagonGrid.html

1) 我不明白底层的几何数据是什么。如图所示，获取美国地图的来源是什么？什么是数据库架构？如果我只需要美国边界而不是每个州，我认为这可能是一个记录？

2) 结果是点列表吗？还是几何向量？

3）如果是geom向量，如何将它们转换为lat和lng的点？

4) 如何将六边形逼近到一个点的半径为 50 英里？

更新：

根据下面的 Jim Jones 示例，我使用宽度来尝试获得正确数量的六边形。不幸的是，出了点问题..

1) 长度似乎与米无关

2) 有多个大小的六边形，看起来很奇怪。

postgis_test=# WITH j AS (
postgis_test(# SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
postgis_test(#   SELECT 
postgis_test(#   generate_hexgrid(
postgis_test(#     5909968.8,
postgis_test(#     ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
postgis_test(# FROM usa_states)i) 
postgis_test-# SELECT count(j.hex) FROM j,usa_states
postgis_test-# WHERE ST_Intersects(usa_states.geom,j.hex);
 count 
-------
   119
(1 row)

postgis_test=# WITH j AS (
postgis_test(# SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
postgis_test(#   SELECT 
postgis_test(#   generate_hexgrid(
postgis_test(#     5909968.8,
postgis_test(#     ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
postgis_test(# FROM usa_states)i) 
postgis_test-# SELECT DISTINCT st_area(j.hex) FROM j,usa_states
postgis_test-# WHERE ST_Intersects(usa_states.geom,j.hex);
     st_area      
------------------
 1219.78281686003
 2089.11341619338
 2089.11341619338
 3379.93344444246
  7051.4650344734
 12076.9943663072
(6 rows)

Answer 1

根据author，以下函数应创建一个网格，其范围基于给定的 BBOX 和以米为单位的单元格大小。

SRID 3857 单位是 [非常近似] 米，使用这个 投影将创建在 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 的几何图形，您应该能够使用以下查询创建一个与美国地图重叠的网格：

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