r - 在 R 中创建等距空间网格
问题描述
我正在尝试使用 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 来做到这一点?
注意:看起来棘手的部分是在一个非常大的区域内保持网格真正 1 公里 x 1 公里。我可以使网格保持十进制度的相等尺寸,但这不是地面上的相等距离。
由于这里有一个狡猾的答案,我已经能够使用 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
,但我对此并不过分兴奋......
解决方案
考虑这个例子;它使用布拉格市的边界作为网格的基础。
关键部分是确保您的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)
推荐阅读
- javascript - 如何在 cypress 中进行依赖 cy.request 调用
- altair - 数据通过本地 json 文件嵌入
- curl - 通过 puppet 3.7.4 中的环境变量动态设置 curl 路径返回“找不到命令''”
- numpy - 为什么向 matplotlib 图添加子图会导致 numpy 崩溃?
- python-2.7 - 将多个号码配置到 SIP 中继的简单 python 脚本
- python - 使用 pandas str.match 转义缺少的括号
- selenium - 在 Firefox 浏览器上获取 ElementClickIntercepted 异常,因为弹出窗口正在干扰(chrome 工作正常)
- javascript - 猫鼬将所有倾斜的 ObjectIds 转换为字符串格式
- c# - Coap 服务器未经授权的访问 C#
- c# - 生成素数序列,Loop Break vs LINQ TakeWhile