如何将世界分割成单元格（网格）

Question

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

我相信已经有一个库可以做到这一点，将单元格标记为S1、S2等等。

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

谢谢！

- Slackware

---

“...以便将每个纬度和经度映射到不同的单元格。”定义“每个纬度和经度”。 - Mike Sherrill 'Cat Recall'

Answer 1

PostGIS 3.1+

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

示例数据

考虑以下多边形：

创建一个以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

受到这篇文章的启发，我开始编写一个函数来实现这个功能——它仍然需要一些调整，但肯定会给你提供一个方向。下面的函数创建一个网格，其单元格大小覆盖给定几何区域：

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)

---

其他选择

迈克的渔网函数基本上是一样的，但您需要手动提供行数和列数以及左下角的坐标对：

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

您可以使用makegrid_2d函数来创建一个基于多边形的区域网格，例如，一个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函数。它创建了一个具有固定大小单元的六边形网格，覆盖给定的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);

更多阅读：

• 导入世界Shapefile
• 等待PostGIS 3.1：网格生成器，作者Paul Ramsey

- Jim Jones

Answer 2

Jim的回答很棒。然而，有些情况下您不需要实际的几何图形。就像您提到的那样，您只需要在同一个单元格中映射到相同代码的坐标。因此，您可以调用一个函数，它仅评估将坐标转换为代码的公式，而无需进行昂贵的点-多边形操作，该操作需要O(n)个多边形 - 没有索引;)。这非常方便，可以快速地对数据进行空间聚合。

就我个人而言，我喜欢Uber的H3库，用于这种类型的事情，但我相信S2也会做类似的事情，并且做得很好。 H3有一个维护良好的PostgreSQL绑定，一个简单的聚合示例可能如下所示：

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的扩展构建库和扩展本身。

---

有人在评论中要求一个点在多边形内的例子。这与问题不完全相关，但突出了如何使用H3。

多边形图层准备：

CREATE TABLE poly_h3 AS (
    SELECT id, h3_polyfill(geom_4326, 13) AS h3res13
    FROM poly
);
CREATE INDEX ON poly_h3 (h3res13);

点图层准备：

ALTER TABLE points ADD COLUMN h3res13 h3index;

UPDATE points
SET h3res13 = h3_geo_to_h3(geom_p_4326, 13);

CREATE INDEX ON points (h3res13);

计算每个多边形的点数：

SELECT poly.id, COUNT(*) AS n
FROM poly_h3 AS poly
    INNER JOIN points AS x
        ON poly.h3res13 = x.h3res13
GROUP BY poly.id;

- chris

---

H3看起来很不错...怎么我直到现在才听说它！谢谢分享+1 - Jim Jones

啊，吉姆，世界真小啊——我们可能错过了一年。抱歉，我通常不会在Stackoverflow上跟踪人，但是我在IfGI度过了非常愉快的时光！ - chris

这真是个好巧合。我几年前“离开”了ifgi，但我仍然在大学工作。我在那里度过了最美好的时光...我非常想念它。向伦敦问候 :) - Jim Jones

嘿，Chris，你能给出一个你在第一段描述的函数的例子吗？这正是我正在寻找的，但我需要一个起点 :) - four-eyes

@Stophface 我在最后添加了一个例子 - chris

谢谢 Redu，谢谢。这下明白了！ - four-eyes