Convert regular lat/lon grid to polygons using R sf

Question

I'm trying to convert a regular 1-degree by 1-degree geographic grid to polygons while maintaining the data associated with each grid point using the R sf package.

Example code:

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))

This yields:

# Attempt #1 
#https://stackoverflow.com/questions/48383990/convert-sequence-of-longitude-and-latitude-to-polygon-via-sf-in-r


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

Which yields

# 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
# https://stackoverflow.com/questions/52669779/convert-sets-of-spatial-coordinates-to-polygons-in-r-using-sf/52671103


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

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

I'm expecting a polygon whose boundaries are defined by the regular grid but the boundaries of each polygon are not being closed. Is there a simple solution in the SF package I'm missing?

Also, when creating the polygon, I'd like for the variable "data" to be associated with each polygon. I have not made it that far, however.

What do you expect the polygon to look like? How do you want to make a polygon out of a grid? — camille, Feb 08 '22 at 14:33

mdag02 · Accepted Answer · 2022-02-08T15:15:12.127

Using 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")

You can adjust the offset depending of what you want (centers or corners)

score 0 · Answer 2 · answered Feb 08 '22 at 13:49

Probably the most convenient way to do this is to use stars, stars objects can be plotted directly with ggplot or you can make spatial polygons out of it:

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...

^{Created on 2022-02-08 by the reprex package (v2.0.1)}

Thanks. When I run the ggplot+geom_stars() command I receive the following error: "Error in if (any(names(lst) == "")) { : missing value where TRUE/FALSE needed" In addition, I'm trying to achieve polygons that are defined by the corners of the grid cells in the first image. — JFD, Feb 08 '22 at 14:04
Without a reproducible example its hard to say what the cause is, is your software up to date? — Bart, Feb 08 '22 at 14:08

Convert regular lat/lon grid to polygons using R sf

2 Answers2