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