Let me start by saying I’ve seen the other thread on a similar question, and even after replicating the nc maps problem they use, I’m still running into a problem with my own code.

I’m trying to create a map with 50km square grids imposed over specific countries. I can generate a map with square grids just fine, but I’m running into problems getting the cellsize argument in st_make_grid() to use sensible metrics.

I understand that I want the maps to use the same projection, and given the units I want one that’s using meters. I’ve checked this repeatedly using st_crs() and it indicates all of the objects are using meters as intended.

This first chunk of code lists the countries I want to include, and lists the CRS code I want to use. I convert each of the three maps returned by gadm() to sf objections, combine the data frames, and apply the CRS code.

library(tidyverse)
library(geodata)

countrylist <- list("DEU", "FRA", "GBR")
wgseqproj <- "EPSG:4087"

tempmap <- countrylist |>
    map(~ gadm(country = .x, level = 0, path = here::here())) |>
    map_dfr(~ st_as_sf(.x)) |>
    sf::st_set_crs(wgseqproj) |>
    sf::st_transform(wgseqproj)

The next chunk of code takes the grid and limits it to the borders of the countries in the list. I take the map from above, make the grid, and use st_intersection() to limit the grid to the country borders. The last line just creates a numeric index for each box in the grid.

  tempgrid <- st_make_grid(tempmap, cellsize = c(0.5, 0.5), square = TRUE, crs = wgseqproj) |> # grid of points
    st_intersection(tempmap) |>
    st_as_sf() |>
    mutate(names = seq(n():1))

ggplot(tempgrid) + geom_sf()

Map resulting from code chunks above

The code runs fine, and the resulting image looks normal, but the cellsize = c(0.5, 0.5) is the section where I’m running into problems. The maps are in meters, but I don’t know what this is capturing (i.e. what does 0.5 translate to?). When I try to run something like cellsize = c(50000, 50000) to get 50km square grids the resulting grids are so large they take up the entire map and I just get country border outlines.

Other code I see uses similarly "large" numbers and the results look great. Is it something about the projection or the gadm maps that’s throwing this off? When when I run more stripped down code on a single country and get rid of the map() functions to simplify things I get the same strange results. I’ve also tried multiple CRS codes. Any helps is much appreciated!

>Solution :

The problem seems to be that you are overwriting the correct crs with st_set_crs. Instead, try:

library(tidyverse)
library(geodata)
library(sf)

countrylist <- list("DEU", "FRA", "GBR")
wgseqproj <- "EPSG:4087"

tempmap <- countrylist |>
  map(~ gadm(country = .x, level = 0, path = here::here())) |>
  map_dfr(~ st_as_sf(.x)) |>
  sf::st_transform(wgseqproj)

tempgrid <- st_make_grid(tempmap, cellsize = 50000, 
                         square = TRUE, crs = wgseqproj) |> 
  st_intersection(tempmap) |>
  st_as_sf() |>
  mutate(names = seq(n():1))

ggplot(tempgrid) + geom_sf()