我有 EURO-CORDEX 气候数据,该数据位于 11 度旋转的极网格上。我已通过将投影转换为 WGS84 来预先准备好这些数据。数据以点的形式出现,代表方形网格的质心。我需要创建围绕这些点的方形网格。我已经导出了实现此目的的通用方法,但网格单元的最终区域显示的误差高达 50%。
我的代码如下。之前我曾被告知要以 tidyverse 表示法提供代码,因此我的目标是尽可能将其删除。数据和代码在github上:https://github.com/avisserquinn/exampleData
首先,我从 csv 加载质心的经度和纬度,并使用 WGS84 投影转换为空间特征数据框。这些点应代表 11 x 11 度或 12 x 12 公里的网格。
> library(tidyverse)
> library(sf)
>
> data <- read_csv("stackExample.csv", col_types = cols())
> data <- st_as_sf(data, coords = c("lon", "lat")) # Spatial feature data frame
> data <- st_set_crs(data, 4326) # Set projection
> data
Simple feature collection with 2221 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: -9.996 ymin: 50.051 xmax: 1.965 ymax: 61.938
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
# A tibble: 2,221 x 2
grid geometry
<dbl> <POINT [°]>
1 1 (-9.996 51.768)
2 2 (-9.979 53.544)
3 3 (-9.96 52.013)
4 4 (-9.931 51.666)
5 5 (-9.924 52.258)
6 6 (-9.912 53.442)
7 7 (-9.906 54.034)
8 8 (-9.895 51.91)
9 9 (-9.875 53.687)
10 10 (-9.869 54.278)
# ... with 2,211 more rows
>
> ggplot(data) + geom_sf() + theme_bw()
我通过应用 st_make_grid 两次(来自 sf 空间特征包)来创建网格。第一次,我找到了点之间的中心。第二次,我找到网格角,这样这些点现在就是质心。
> cellsize = .11
> dataGrid <- st_make_grid(data, cellsize = cellsize, what = "centers")
> dataGrid <- st_make_grid(dataGrid, cellsize = cellsize)
> dataGrid <- dataGrid %>% as_tibble %>% st_as_sf
> dataGrid
Simple feature collection with 11772 features and 0 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: -9.941 ymin: 50.106 xmax: 1.939 ymax: 62.096
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 10 features:
geometry
1 POLYGON ((-9.941 50.106, -9...
2 POLYGON ((-9.831 50.106, -9...
3 POLYGON ((-9.721 50.106, -9...
4 POLYGON ((-9.611 50.106, -9...
5 POLYGON ((-9.501 50.106, -9...
6 POLYGON ((-9.391 50.106, -9...
7 POLYGON ((-9.281 50.106, -9...
8 POLYGON ((-9.171 50.106, -9...
9 POLYGON ((-9.061 50.106, -8...
10 POLYGON ((-8.951 50.106, -8...
接下来,我将网格数据与质心聚合,以仅查找匹配的网格单元。
> dataGrid <- aggregate(data, dataGrid, FUN = mean)
> dataGrid <- as_tibble(dataGrid)
> dataGrid <- dataGrid[!is.na(dataGrid$grid),]
> dataGrid$area_sqm = st_area(dataGrid$geometry)
> dataGrid$area_sqkm = as.numeric(unlist(dataGrid$area_sqm * 10^-6))
> dataGrid$area_deficit = (12*12) - dataGrid$area_sqkm
> dataGrid
# A tibble: 2,175 x 5
grid geometry area_sqm area_sqkm area_deficit
<dbl> <POLYGON [°]> [m^2] <dbl> <dbl>
1 656 ((-5.651 50.106, -5.541 50.106, -5.541 50.216, -5.651 50.216, -5.651 50.106)) 96173304 96.2 47.8
2 678 ((-5.431 50.106, -5.321 50.106, -5.321 50.216, -5.431 50.216, -5.431 50.106)) 96173304 96.2 47.8
3 702 ((-5.211 50.106, -5.101 50.106, -5.101 50.216, -5.211 50.216, -5.211 50.106)) 96173304 96.2 47.8
4 730 ((-5.101 50.106, -4.991 50.106, -4.991 50.216, -5.101 50.216, -5.101 50.106)) 96173304 96.2 47.8
5 693 ((-5.321 50.216, -5.211 50.216, -5.211 50.326, -5.321 50.326, -5.321 50.216)) 95954257 96.0 48.0
6 720 ((-5.101 50.216, -4.991 50.216, -4.991 50.326, -5.101 50.326, -5.101 50.216)) 95954257 96.0 48.0
7 762 ((-4.881 50.216, -4.771 50.216, -4.771 50.326, -4.881 50.326, -4.881 50.216)) 95954257 96.0 48.0
8 1044 ((-3.891 50.216, -3.781 50.216, -3.781 50.326, -3.891 50.326, -3.891 50.216)) 95954257 96.0 48.0
9 712 ((-5.211 50.326, -5.101 50.326, -5.101 50.436, -5.211 50.436, -5.211 50.326)) 95734844 95.7 48.3
10 746. ((-4.991 50.326, -4.881 50.326, -4.881 50.436, -4.991 50.436, -4.991 50.326)) 95734844 95.7 48.3
# ... with 2,165 more rows
当我绘制最终输出时,您可以看到问题。网格单元尺寸有很大偏差(赤字);由于投影是弯曲的,我预计与 12 x 12 公里会有一些偏差,但这种程度的赤字是极端的。另外,网格中存在间隙;我推测这是因为网格单元的大小不正确意味着并非所有点都被捕获?
> ggplot(dataGrid) +
+ geom_sf(aes(geometry = geometry, fill = area_deficit)) +
+ theme_minimal() +
+ scale_fill_viridis_c(direction = -1) +
+ scale_colour_viridis_c(direction = -1) +
+ coord_sf()
我尝试了多种方法来尝试解决此问题,但没有成功。任何有关替代解决方案的建议将不胜感激。
最佳答案
st_make_grid()
仅使用第一个参数的边界框,因此您只需在每个方向上将边界框扩展单元格大小的 1/2:
# Generate some centroids
centroids <- st_make_grid(what="centers") %>% st_sf()
# Make a new grid from them
cellSize <- 10
grid <- (st_bbox(centroids) + cellSize/2*c(-1,-1,1,1)) %>%
st_make_grid(cellsize=c(cellSize, cellSize)) %>% st_sf()
ggplot() + geom_sf(data=grid) + geom_sf(data=centroids)
这篇文章目前已经很旧了,因此链接已损坏,但单元格大小 0.11 看起来不正确。希望你能解决这个问题。
关于r - 使用 sf 围绕点(质心)创建网格,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/58800327/