如何将世界划分为单元格（网格）

Question

如何将世界划分为大小几乎相等的单元格，以便每个纬度、经度可以映射到不同的单元格？

我很确定我见过一个库这样做，将单元格标记为 S1、S2 等。

假设我们有 62.356279,-99.422395 ，如何将其映射到名为“FR,23”的 2km*2km 单元格？

谢谢！

Answer 1

Jim 的回答非常好。但是，在某些用例中，您不需要实际的几何形状。正如您所提到的，您所需要的只是将同一单元格中的坐标映射到同一代码。因此，您无需为 n 多边形采用 O(n) 的昂贵的多边形点操作 - 没有索引，即 ;) - 您调用一个简单评估的函数O(1) 中将坐标转换为代码的公式。在空间上快速聚合数据非常方便。

就我个人而言，我喜欢 Uber 的 H3 库，但我确信 S2 做了类似的事情并且做得很好。 H3 有一个维护良好的 PostgreSQL binding，一个简单的聚合示例如下所示：

SELECT h3_geo_to_h3(geom_4326, 9) AS h3res09, SUM(pop_19) AS pop_19
FROM uk_postcode_population
GROUP BY 1;

阅读：为每个分辨率 9 六边形总结英国邮政编码级别的人口。

您仍然可以在需要时创建实际的六边形几何形状（甚至是整个六边形网格）。但根据我的经验，一旦您采用网格方法，最终您将只需要多边形来进行可视化。

我应该注意，你不能使用这个库自己划分世界 - Uber 已经为你划分了它。因此，如果 2 公里方格是硬性要求，那么这不适合您。

安装 H3 扩展并不像CREATE EXTENSION postgis 那样简单，但您也不必成为命令行勇士。您至少必须安装PGXN，很可能是PostgreSQL's extension build library 和extension 本身。

Answer 2

PostGIS 3.1+

PostGIS 3.1 引入了非常易于使用的网格生成器，即ST_SquareGrid 和ST_HexagonGrid。将这些函数与表中的数据一起使用的一种简单方法是使用LATERAL 来执行它，例如大小为 0.1° 的单元格::

样本数据

考虑以下polygon：

在给定几何体上创建一个单元格大小为 0.1° 的方格

SELECT grid.* FROM isle_of_man imn, 
LATERAL ST_SquareGrid(0.1,imn.geom) grid;

如果您只想要与几何体相交的单元格，只需在 WHERE 子句中调用函数 ST_Intersects：

SELECT grid.* FROM isle_of_man imn, 
LATERAL ST_SquareGrid(0.1,imn.geom) grid
WHERE ST_Intersects(imn.geom,grid.geom);

同样的原则也适用于ST_HexagonGrid：

SELECT grid.* FROM isle_of_man imn, 
LATERAL ST_HexagonGrid(0.1,imn.geom) grid;

SELECT grid.* FROM isle_of_man imn, 
LATERAL ST_HexagonGrid(0.1,imn.geom) grid
WHERE ST_Intersects(imn.geom,grid.geom);

旧的 PostGIS 版本

受post 的启发，我开始编写一个函数来做到这一点——它仍然需要一些调整，但它肯定会给你一个方向。以下函数创建一个网格，其中包含给定大小的单元格，覆盖给定几何图形的区域：

CREATE OR REPLACE FUNCTION public.generate_grid(_size numeric,_geom geometry)
RETURNS TABLE(gid bigint, cell geometry) LANGUAGE 'plpgsql'
AS $BODY$
DECLARE
  _bbox box2d := ST_Extent(_geom);
  _ncol int := ceil(abs(ST_Xmax(_bbox)-ST_Xmin(_bbox))/_size);
  _nrow int := ceil(abs(ST_Ymax(_bbox)-ST_Ymin(_bbox))/_size);
  _srid int DEFAULT 4326;
BEGIN  
  IF ST_SRID(_geom) <> 0 THEN
    IF EXISTS (SELECT 1 FROM spatial_ref_sys crs
               WHERE crs.srid = ST_SRID(_geom) AND NOT crs.proj4text LIKE '+proj=longlat%') THEN
      RAISE EXCEPTION 'Only lon/lat spatial reference systems are supported in this function.';
    ELSE
      _srid := ST_SRID($2); 
    END IF; 
  END IF;
  
  RETURN QUERY 
   SELECT ROW_NUMBER() OVER (), geom FROM (
    SELECT 
     ST_SetSRID(
      (ST_PixelAsPolygons(
        ST_AddBand(
          ST_MakeEmptyRaster(_ncol, _nrow, ST_XMin(_bbox), ST_YMax(_bbox),_size), 
          '1BB'::text, 1, 0), 
         1, false)).geom,_srid)) j(geom);
END;
$BODY$;

注意：此功能依赖于扩展PostGIS Raster。

SELECT cell FROM isle_of_man, 
LATERAL generate_grid(0.1,geom);

...如果您只对与多边形重叠的单元格感兴趣，请在查询中添加ST_Intersects：

SELECT cell FROM isle_of_man, 
LATERAL generate_grid(0.1,geom)
WHERE ST_Intersects(geom,cell)

---

其他选择

Mike 的fishnet function 基本相同，但需要手动提供行数和列数，以及左下角的坐标对：

SELECT ST_SetSRID(cells.geom,4326)
FROM ST_CreateFishnet(4, 6, 0.1, 0.1,-4.8411, 54.0364) AS cells;

您可以使用此makegrid_2d function 在使用多边形的区域上创建网格，例如一个包含 5000 米大小的单元格的网格：

CREATE TABLE grid_isle_of_man AS
SELECT 'S'||ROW_NUMBER() OVER () AS grid_id, (g).geom 
FROM (
  SELECT ST_Dump(makegrid_2d(geom,5000))
  FROM isle_of_man) j(g)
JOIN isle_of_man ON ST_Intersects((g).geom,geom);

同样的逻辑适用于这个hexagrid function。它在给定的 BBOX 上创建具有固定大小单元的六边形网格。您可以手动提供 BBOX（函数的第二个参数）或从给定的多边形中提取它。例如，要创建一个与多边形范围相匹配的六边形网格，并将其存储在带有所需标签的新表中 - 单元格大小为 0.1°：

CREATE TABLE hexgrid_isle_of_man AS 
WITH j (hex_rec) AS (
  SELECT generate_hexagons(0.1,ST_Extent(geom)) 
  FROM isle_of_man
)
SELECT 'S'||ROW_NUMBER() OVER () AS grid_id,(hex_rec).hexagon FROM j
JOIN isle_of_man t ON ST_Intersects(t.geom,(hex_rec).hexagon);

进一步阅读：

• Import World Shapefile
• Waiting for PostGIS 3.1: Grid Generators，作者：保罗·拉姆齐