在R中创建等距空间网格

Question

我正在尝试使用R在给定区域上制作等尺寸的正方形网格。我希望我的网格为1km x 1km正方形。我看到这样的示例，它们说明了相等的经/纬网格：

在空间范围内创建规则的多边形网格，旋转给定角度

但这还不算大小。看来我应该能够使用该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做到这一点？

注意：看来棘手的部分是在很大的区域内将网格保持为1km x 1km。我可以以十进制度度保持网格相等的尺寸，但是那不等于地面上的距离。

由于这里有一个巧妙的答案，我已经能够使用PostGIS做到这一点。在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，但对此我并不太兴奋...

Answer 1

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

关键部分是：
-确保您的sf对象位于公制CRS中（如果您是美国人并且感觉爱国，则应选择脚尖的对象）
-计算所需的平方数
-运行代码并享受...

library(tidyverse)
library(RCzechia)
library(sf)


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

krabicka <- st_bbox(mesto)

xrange <- krabicka$xmax - krabicka$xmin
yrange <- krabicka$ymax - krabicka$ymin

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

rozmery <- c(xrange/grid_spacing , yrange/grid_spacing) %>% # number of polygons necessary
  ceiling() # rounded up to nearest integer

polygony <- st_make_grid(mesto, square = T, n = rozmery) %>% # the grid, covering bounding box
  st_intersection(mesto) %>% # only the inside part 
  st_sf() # not really required, but makes the grid nicer to work with later

plot(polygony, col = 'white')