postgresql - PostGIS可以用来创建一个国家的网格地图吗？

Question

我正在尝试将一个区域的地图划分为 1 公里 x 1 公里平方的网格，并最终计算和检索有多少点（给定它们的纬度和经度值）落入该地图的每个方形网格中。这个操作可以在 PostGIS 中完成，如果可以，我该怎么做？

更新：Mike Toews在这里有一个详细的答案：

https://gis.stackexchange.com/questions/16374/how-to-create-a-regular-polygon-grid-in-postgis

Answer 1

正如我的评论中提到的， 制作一个常规网格。要为整个国家制作 1 公里的网格，这可能具有挑战性，因为地球不是平坦的，无法划分为完美的 1 公里网格。

要制作 1 公里的网格，您需要一个投影坐标系，长度单位为米。WGS84 (EPSG:4326)不能这样做，因为它的单位是纬度/经度。要找到合适的投影系统，您需要找到“等面积”投影，例如兰伯特方位角等面积投影（LAEA）。例如，整个欧洲都可以使用ETRS-LAEA (EPSG:3035)，尽管某些部分可能存在一些失真。或者如果在新西兰，新西兰横轴墨卡托 2000。每个区域通常都有一个很好的投影可以使用。

要运行 PostGIS 查询，您需要使用ST_Transform(geom, 3035)（例如，对于 ETRS-LAEA）将几何图形投影到网格上。

Answer 2

我编写了一个 PostGIS 函数来在另一层之上生成十六进制网格。

DO $$
DECLARE
  _curs   CURSOR FOR SELECT geom3857 FROM nrw;
  _table  TEXT    := 'nrw_hx_10k';
  _srid   INTEGER := 3857;
  _height NUMERIC := 10000;
  _width  NUMERIC := _height * 0.866;
  _geom   GEOMETRY;
  _hx     TEXT    := 'POLYGON((' || 0 || ' ' || 0 || ',' || (_width * 0.5) || ' ' || (_height * 0.25) || ',' ||
                 (_width * 0.5) || ' '
                 || (_height * 0.75) || ',' || 0 || ' ' || _height || ',' || (-1 * (_width * 0.5)) || ' ' ||
                 (_height * 0.75) || ',' ||
                 (-1 * (_width * 0.5)) || ' ' || (_height * 0.25) || ',' || 0 || ' ' || 0 || '))';
  _hx_g   GEOMETRY := ST_SetSRID(_hx::GEOMETRY, _srid);

BEGIN
  CREATE TEMP TABLE hx_tmp (geom GEOMETRY(POLYGON));

  OPEN _curs;
  LOOP
    FETCH
    _curs INTO _geom;
    EXIT WHEN NOT FOUND;

    INSERT INTO hx_tmp
      SELECT
        ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON) geom
      FROM
        generate_series(
          (st_xmin(_geom) / _width)::INTEGER * _width - _width,
          (st_xmax(_geom) / _width)::INTEGER * _width + _width,
          _width) x_series,
        generate_series(
          (st_ymin(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) - _height,
          (st_ymax(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) + _height,
          _height * 1.5) y_series
      WHERE
        ST_Intersects(ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON), _geom);

    INSERT INTO hx_tmp
      SELECT ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON) geom
      FROM
        generate_series(
          (st_xmin(_geom) / _width)::INTEGER * _width - (_width * 1.5),
          (st_xmax(_geom) / _width)::INTEGER * _width + _width,
          _width) x_series,
        generate_series(
          (st_ymin(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) - (_height * 1.75),
          (st_ymax(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) + _height,
          _height * 1.5) y_series
      WHERE
        ST_Intersects(ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON), _geom);

  END LOOP;
  CLOSE _curs;

  CREATE INDEX sidx_hx_tmp_geom ON hx_tmp USING GIST (geom);
  EXECUTE 'DROP TABLE IF EXISTS '|| _table;
  EXECUTE 'CREATE TABLE '|| _table ||' (geom GEOMETRY(POLYGON, '|| _srid ||'))';
  EXECUTE 'INSERT INTO '|| _table ||' SELECT * FROM hx_tmp GROUP BY geom';
  EXECUTE 'CREATE INDEX sidx_'|| _table ||'_geom ON '|| _table ||' USING GIST (geom)';
  DROP TABLE IF EXISTS hx_tmp;
END $$;

Input parameters must be set for: _curs: The geometry field name and the table name of input geometries. _table: The name of the output table. _srid: The geographic projection code of the input (and output) geometries. _height: The height of a hexgon grid cell in projection units.

I generate a scaled hexagon geometry in the declare block and then loop through the input geometries. In the loop, I generate series for the x and y extent plus some for each input geometry. The hexagon is translated and inserted into a temporary table if the two geometries intersect. A second pair of series is generated alternative offset rows. Finally I group the hexagon grid cells by their geometries to remove duplicates.

There is a more detailed description and some background on medium.