How to make an equal area grid in r using the terra package

Question

I want to make an equal area grid (400 square miles per grid cell) over Wisconsin. I am doing this using the code from this link: Creating an equal distance spatial grid in R.

But, this code isn't very flexible, and I also need the grid to be more than just polygons. I need it to be a shapefile. I like the Terra package, but am unable to figure out how to do this in the terra package. The WI shapefile can be downloaded from https://data-wi-dnr.opendata.arcgis.com/datasets/wi-dnr::wisconsin-state-boundary-24k/explore.

My code looks like this:

library(sf)
library(terra)
library(tidyverse)

wi_shape <- vect('C:\\Users\\ruben\\Downloads\\Wisconsin_State_Boundary_24K\\Wisconsin_State_Boundary_24K.shp')

plot(wi_shape)
wi_grid <- st_make_grid(wi_shape, square = T, cellsize = c(20 * 1609.344, 20 * 1609.344))
plot(wi_grid, add = T)

How do I define a grid that is centered on a lat/lon point, where the output is a shapefile that contains attributes for each grid cell? I'm not sure why this is so confusing to me. Thank you.

Answer 1

If your goal is to make a raster based on the extent of another spatial dataset (polygons in this case) you can do

library(terra)
wi <- vect('Wisconsin_State_Boundary_24K.shp')
r <- rast(wi, res=(20 * 1609.344))

You can turn these into polygons and write them to a file with

v <- as.polygons(r)
writeVector(v, "test.shp")

To define a lon/lat center for the grid, you could do the following.

Coordinates of an example lon/lat point projected to the crs of your polygons (Wisconsin Transverse Mercator).

center <- cbind(-90, 45) |> vect(crs="+proj=longlat")
cprj <- crds(project(center, wi))
res <- 20 * 1609.344

Create a single cells around that point and expand the raster:

e <- rep(cprj, each=2) + c(-res, res) / 2 
x <- rast(ext(e), crs=crs(wi), ncol=1, nrow=1)
x <- extend(x, wi, snap="out")

The result

plot(as.polygons(x), border="blue")
lines(wi, col="red")
points(cprj, pch="x", cex=2)

I should also mention that you are not using an equal-area coordinate reference system. You can see the variation in cell sizes with

a <- cellSize(x)

But it is very small (less than 1%) relative to the average cell size

diff(minmax(a))
#       area
#max 1690441

global(a, mean)
#           mean
#area 1036257046

Answer 2

Let's try to tidy this a little bit.

[...] and I also need the grid to be more than just polygons. I need it to be a shapefile.

It's exactly the other way around from my point of view. Once you obtained a proper representation of a polygon, you can export it in whatever format you like (which is supported), e.g. an ESRI Shapefile.

I like the Terra package, but am unable to figure out how to do this in the terra package.

Maybe you did not notice, but actually you are not really using {terra} to create your grid, but {sf} (with SpatVector input from terra, which is accepted here).

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.6.33

wi_shape <- vect('Wisconsin_State_Boundary_24K.shp')
class(wi_shape)
#> [1] "SpatVector"
#> attr(,"package")
#> [1] "terra"

wi_grid <- st_make_grid(wi_shape, square = T, cellsize = c(20 * 1609.344, 20 * 1609.344))
class(wi_grid)
#> [1] "sfc_POLYGON" "sfc"

It's a minor adjustment, but basically, you can cut this dependency here for now. Also - although I'm not sure is this is the type of flexibility you are looking for - I found it very pleasing to work with {units} recently if you are about to do some conversion stuff like square miles in meters. In the end, once your code is running properly, you can substitute your hardcoded values by variables step by step and wrap a function out of this. This should not be a big deal in the end.

In order to shift your grid to be centered on a specific lat/lon point, you can leverage the offset attribute of st_make_grid(). However, since this only shifts the grid based on the original extent, you might lose coverage with this approach:

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

wi_shape <- read_sf("Wisconsin_State_Boundary_24K.shp")

# area of 400 square miles
A <- units::as_units(400, "mi^2")

# boundary length in square meters to fit the metric projection
b <- sqrt(A)
units(b) <- "m"

# let's assume you wanted your grid to be centered on 45.5° N / 89.5° W
p <- c(-89.5, 45.5) |> 
  st_point() |> 
  st_sfc(crs = "epsg:4326") |> 
  st_transform("epsg:3071") |> 
  st_coordinates()
p
#>          X        Y
#> 1 559063.9 558617.2

# create an initial grid for centroid determination
wi_grid <- st_make_grid(wi_shape, cellsize = c(b, b), square = TRUE)

# determine the centroid of your grid created
wi_grid_centroid <- wi_grid |> 
  st_union() |> 
  st_centroid() |> 
  st_coordinates()
wi_grid_centroid
#>          X        Y
#> 1 536240.6 482603.9

# this should be your vector of displacement, expressed as the difference 
delta <- wi_grid_centroid - p 
delta
#>           X        Y
#> 1 -22823.31 -76013.3

# `st_make_grid(offset = ...)` requires lower left corner coordinates (x, y) of the grid,
# so you need some extent information which you can acquire via `st_bbox()`
bbox <- st_bbox(wi_grid)

# compute the adjusted lower left corner
llc_new <- c(st_bbox(wi_grid)["xmin"] + delta[1], st_bbox(wi_grid)["ymin"] + delta[2])

# create your grid with an offset
wi_grid_offset <- st_make_grid(wi_shape, cellsize = c(b, b), square = TRUE, offset = llc_new) |> 
  st_as_sf()

# append attributes
n <- dim(wi_grid_offset)[1]

wi_grid_offset[["id"]] <- paste0("A", 1:n)
wi_grid_offset[["area"]] <- st_area(wi_grid_offset) |> as.numeric()

# inspect
plot(st_geometry(wi_shape))
plot(st_geometry(wi_grid_offset), border = "red", add = TRUE)

If you wanted to export your polygon features ("grid") in shapefile format, simply make use of st_write(wi_grid_sf, "wi_grid_sf.shp").

PS: For this example you need none of the tidyverse stuff, so there is no need to load it.