使用 sf 围绕点（质心）创建网格

我有 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 看起来不正确。希望你能解决。