0

I'm looking to divide the world into grids using the Behrmann projection (CRS 'ESRI:54017'). I came across some R code that seems to do the job:

world_single_cell <- st_polygon(list(
    rbind(c(-180, -90),
          c(180, -90),
          c(180, 90),
          c(-180, 90),
          c(-180, -90)
    )
))
world_single_cell_sfc <- st_sfc(world_single_cell)

# Define the CRS for the world as WGS84 (EPSG:4326)
st_crs(world_single_cell_sfc) <- 4326
worldGrids <- world_single_cell_sfc %>% 
    st_transform("ESRI:54017") %>%
    st_make_grid(n = 3426)

However, I'm unsure if this code is correct. Additionally, I'd like to adjust the number of grids based on the length of each grid cell. Specifically, I'm aiming for resolutions of 96.5 km × 96.5 km, 193 km × 193 km, and 385.9 km × 385.9 km, using equal-area Behrmann projection grids with dimensions of 1° × 1°, 2° × 2°, and 4° × 4° longitude at 30°N latitude.

How can I achieve this in R?

1 Answer 1

1

I think you should be able to make use of the cellsize argument when using st_make_grid().

?st_make_grid:

cellsize
numeric of length 1 or 2 with target cellsize: for square or rectangular cells the width and height, for hexagonal cells the distance between opposite edges (edge length is cellsize/sqrt(3)). A length units object can be passed, or an area unit object with area size of the square or hexagonal cell.

So let's try this:

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(units)

edge <- units::as_units(385.9, "km")

x <- world_single_cell_sfc |>  
  st_transform("ESRI:54017") |> 
  st_make_grid(cellsize = c(edge, edge)) |> 
  st_as_sf()

x <- dplyr::mutate(x, "area" = st_area(x))
x
#> Simple feature collection with 3549 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -17367530 ymin: -7342230 xmax: 17749370 ymax: 7707870
#> Projected CRS: World_Behrmann
#> First 10 features:
#>                                 x               area
#> 1  POLYGON ((-17367530 -734223... 148918810000 [m^2]
#> 2  POLYGON ((-16981630 -734223... 148918810000 [m^2]
#> 3  POLYGON ((-16595730 -734223... 148918810000 [m^2]
#> 4  POLYGON ((-16209830 -734223... 148918810000 [m^2]
#> 5  POLYGON ((-15823930 -734223... 148918810000 [m^2]
#> 6  POLYGON ((-15438030 -734223... 148918810000 [m^2]
#> 7  POLYGON ((-15052130 -734223... 148918810000 [m^2]
#> 8  POLYGON ((-14666230 -734223... 148918810000 [m^2]
#> 9  POLYGON ((-14280330 -734223... 148918810000 [m^2]
#> 10 POLYGON ((-13894430 -734223... 148918810000 [m^2]

Created on 2024-04-24 with reprex v2.1.0

2
  • Thank you dimfalk! The code is very helpful!
    – Emma
    Commented Apr 24 at 17:00
  • @Emma: If this answered your question, feel free to accept my approach to mark this as "completed".
    – dimfalk
    Commented May 5 at 20:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.