在 R 中创建等距空间网格答案

【问题标题】：Creating an equal distance spatial grid in R在 R 中创建等距空间网格
【发布时间】：2019-05-16 06:58:24
【问题描述】：

我正在尝试使用 R 在给定区域上制作一个等尺寸的正方形网格。我希望我的网格为 1 公里 x 1 公里正方形。我看到这样的例子说明了相等的纬度/经度网格：

Creating a regular polygon grid over a spatial extent, rotated by a given angle

但这不是大小。似乎我应该能够使用st_make_grid 函数并创建它，但我不知道如何使网格成为 1km x 1km。

https://r-spatial.github.io/sf/reference/st_make_grid.html

例如，我想从 (37,-89.2) 开始，到 (36.2,-86.8) 结束，并创建一个 1km x 1km 的等间距网格。我将如何使用 R 来做到这一点？

注意：看起来棘手的部分是在一个非常大的区域内保持网格实际上为 1 公里 x 1 公里。我可以使网格保持十进制度的相等尺寸，但这不是地面上的相等距离。

感谢 crafty answer here. 在 PostGIS 我创建了一个函数：

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  width_step integer,
  height_step integer
)
RETURNS public.geometry AS
$body$
DECLARE
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  NextX DOUBLE PRECISION;
  NextY DOUBLE PRECISION;
  CPoint public.geometry;
  sectors public.geometry[];
  i INTEGER;
  SRID INTEGER;
BEGIN
  Xmin := ST_XMin(bound_polygon);
  Xmax := ST_XMax(bound_polygon);
  Ymax := ST_YMax(bound_polygon);
  SRID := ST_SRID(bound_polygon);

  Y := ST_YMin(bound_polygon); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
      NextX := ST_X(ST_Project(CPoint, $2, radians(90))::geometry);
      NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);

      i := i + 1;
      sectors[i] := ST_MakeEnvelope(X, Y, NextX, NextY, SRID);

      X := NextX;
    END LOOP xloop;
    CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
    NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);
    Y := NextY;
  END LOOP yloop;

  RETURN ST_Collect(sectors);
END;
$body$
LANGUAGE 'plpgsql';

然后我可以调用它并传递一个多边形：

SELECT (
    ST_Dump(
      makegrid_2d(
        ST_GeomFromText(
          'Polygon((-75 42, -75 40, -73 40, -73 42, -75 42))',
          4326
        )  ,
         1000, -- width step in meters
         1000  -- height step in meters
       ) 
    )
  ) .geom AS cell into test_grid_cell;

但正如您所见，即使使用 PostGIS，这也不是固定的例程。通过一些努力，我想我可以将它移植到 sf，但我并没有对此感到过度兴奋......

【问题讨论】：

我猜想matrix(c(37, 36.2, -89.2, -86.8), 2) %>% st_multipoint() %>% st_sfc(crs = 4326) %>% st_transform(3857) %>% st_make_grid(cellsize = 1000*1000) %>% plot 之类的东西，但cellsize 的论据并不十分清楚。如果覆盖的区域真的很大，请注意有些正方形不会很正方形。
是的，我的区域很大，所以它最终没有恒定的尺寸（以米为单位）。
是的。针对您关心的内容进行优化；等面积投影可能是个好主意。

标签： r spatial

【解决方案1】：

考虑这个例子；它使用布拉格市的边界作为网格的基础。

关键部分是确保您的 sf 对象在公制 CRS 中（如果您是美国人并且感到爱国，则为尺码）；只要它在预计的 CRS 中就没有关系（即st_is_longlat(x) 返回FALSE）。

library(dplyr)
library(RCzechia) # a package of Czech administrative areas
library(sf)


mesto <- kraje() %>% # All Czech NUTS3 ...
  filter(NAZ_CZNUTS3 == 'Hlavní město Praha') %>% # ... city of Prague
  st_transform(5514) # a metric CRS 

grid_spacing <- 1000  # size of squares, in units of the CRS (i.e. meters for 5514)

polygony <- st_make_grid(mesto, square = T, cellsize = c(grid_spacing, grid_spacing)) %>% # the grid, covering bounding box
  st_sf() # not really required, but makes the grid nicer to work with later

plot(polygony, col = 'white')
plot(st_geometry(mesto), add = T)

【讨论】：