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