使用 sf 在点（质心）周围创建网格答案

【问题标题】：Create grid around points (centroids) using sf使用 sf 在点（质心）周围创建网格
【发布时间】：2019-11-11 11:19:36
【问题描述】：

我有 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()

我尝试了多种方法来尝试解决此问题，但均未成功。任何关于替代解决方案的建议将不胜感激。

【问题讨论】：

标签： r grid spatial sf centroid

【解决方案1】：

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)

Example grid

此时这篇文章已经很老了，因此链接已损坏，但 0.11 的单元格大小看起来不正确。希望你解决了。

【讨论】：