使用 R sf 将常规 lat/lon 网格转换为多边形

Question

我正在尝试将常规的 1 度乘 1 度地理网格转换为多边形，同时使用 R sf 包维护与每个网格点关联的数据。

示例代码：

library(tidyverse)
library(sf)
library(sfheaders)
library(sfheaders)

# latitude
lldata <- tibble(lat=as.numeric(rep(seq(40,49.5),10)))

# Longitude + add data
lldata <- lldata %>%  
          group_by(lat) %>%  
          mutate(lon = -90 + row_number()-1,
                 data = runif(1))

#Plot the grid as points
lldata %>% 
  ggplot() + 
  geom_point(aes(x=lon,y=lat,color=data))

这会产生：

# Attempt #1 
#


 
lldata %>% arrange(lat,-lon)  %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON") %>%
  ggplot() + 
  geom_sf()

哪个产量

# Also tried 
poly <- st_sf(st_sfc(st_make_valid(st_polygon(list(as.matrix(lldata))))), crs = 4326)
# error that the polygon is not closed

# Attempt #2
# 


sf <- sfheaders::sf_polygon(
  obj = lldata
  , x = "lon"
  , y = "lat"
)

# Same result as above
sf %>% 
  ggplot() + 
  geom_sf()

我期待一个边界由规则网格定义的多边形，但每个多边形的边界都没有闭合。我缺少的SF包中是否有简单的解决方案？

此外，在创建多边形时，我希望变量“数据”与每个多边形相关联。然而，我还没有做到这一点。

Answer 1

可能最方便的方法是使用星星，星星对象可以直接用 ggplot 绘制，或者你可以用它制作空间多边形：

library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(sfheaders)
require(stars)
#> Loading required package: stars
#> Loading required package: abind
# latitude
lldata <- tibble(lat=as.numeric(rep(seq(40,49.5),10)))

# Longitude + add data
lldata <- lldata %>%  
  group_by(lat) %>%  
  mutate(lon = -90 + row_number()-1,
         data = runif(1))
ggplot()+geom_stars(data=st_as_stars(lldata, dims  = c('lon','lat')))

st_as_stars(lldata, dims  = c('lon','lat')) %>% st_as_sf(as_points=F)
#> Simple feature collection with 100 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -90.5 ymin: 39.5 xmax: -80.5 ymax: 49.5
#> CRS:           NA
#> First 10 features:
#>         data                       geometry
#> 1  0.5560062 POLYGON ((-90.5 49.5, -89.5...
#> 2  0.5560062 POLYGON ((-89.5 49.5, -88.5...
#> 3  0.5560062 POLYGON ((-88.5 49.5, -87.5...
#> 4  0.5560062 POLYGON ((-87.5 49.5, -86.5...
#> 5  0.5560062 POLYGON ((-86.5 49.5, -85.5...
#> 6  0.5560062 POLYGON ((-85.5 49.5, -84.5...
#> 7  0.5560062 POLYGON ((-84.5 49.5, -83.5...
#> 8  0.5560062 POLYGON ((-83.5 49.5, -82.5...
#> 9  0.5560062 POLYGON ((-82.5 49.5, -81.5...
#> 10 0.5560062 POLYGON ((-81.5 49.5, -80.5...

^{由 reprex package (v2.0.1)}

于 2022-02-08 创建

Answer 2

使用sf::st_make_grid：

library(tidyverse)
library(sf)

lldata <- tibble(lat = as.numeric(rep(seq(40, 49.5), 10))) %>%  
  group_by(lat) %>%  
  mutate(lon = -90 + row_number() - 1,
         data = runif(1)) %>% 
  ungroup() %>% 
  st_as_sf(coords = c("lon", "lat")) 

lldata_poly <- lldata %>% 
  st_make_grid(cellsize = 1, offset = c(-90.5, 39.5)) %>% 
  st_as_sf() %>% 
  st_join(lldata)

# plot
ggplot() +
  geom_sf(data = lldata_poly, aes(fill = data)) +
  geom_sf(data = lldata)

# export
write_sf(lldata_poly, "~/grid.shp")

您可以根据需要调整偏移量（中心或角）

使用 R sf 将常规 lat/lon 网格转换为多边形

Convert regular lat/lon grid to polygons using R sf

r

sf