Fixing common map projection issues in R

February 06, 2021

It's easy to run into projection issues when you using spatial data from R packages without paying attention. In this quick post I will outline a few easy solutions for fixing common projection issues.

dggridr

Let's start with the dggridR package, which generates discrete global grids as SpatialPolygons or data frames. Here I'm generating an Icosahedron Snyder Equal Area (ISEA) grid of resolution 5. I'm using a Robinson projection, but EPSG:4326 will have the same issues.

library(dggridR)

dggs <- dgconstruct(res = 5)
grid <- dgearthgrid(dggs, frame = FALSE) %>%
  st_as_sf()

ggplot() +
  geom_sf(data = grid) +
  coord_sf(crs = "ESRI:54030") +
  theme_void()

dggridr robinson

There's clearly an issue with polygons crossing the date line. Using an orthographic projection we can take a better look at the polygons causing the artefacts on our map.

ggplot() +
  geom_sf(data = grid) +
  geom_sf(data = grid[1153,], color = "#00ccff", fill = "#ccf5ff") +
  geom_sf(data = st_sfc(st_linestring(rbind(c(180, 0), c(180, 90))), crs = 4326), color = "#cc3300") +
  coord_sf(crs = "+proj=ortho +lat_0=90 +lon_0=0") +
  theme_void()

dggridr ortho

To fix this issue we can use st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=90")), which will split polygons that cross the date line. Note the DATELINEOFFSET option which has value 10 by default, and will determine the range around the date line within which polygons will be split. Let's reproduce the previous map, but now with the st_wrap_dateline() applied to the highlighted grid cell:

ggplot() +
  geom_sf(data = grid) +
  geom_sf(data = st_sfc(st_linestring(rbind(c(180, 0), c(180, 90))), crs = 4326), color = "#cc3300") +
  geom_sf(data = grid[1153,] %>% st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=90")), color = "#00ccff", fill = "#ccf5ff") +
  coord_sf(crs = "+proj=ortho +lat_0=90 +lon_0=0") +
  theme_void()

dggridr ortho fixed

After applying this to the entire grid, our Robinson projection map will look a whole lot better:

grid_fixed <- grid %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=270"))

ggplot() +
  geom_sf(data = grid_fixed) +
  coord_sf(crs = "ESRI:54030") +
  theme_void()

dggridr robinson fixed

We could opt to also use st_segmentize() to make the map a little more accurate (i.e. with slightly curved cell boundaries instead of straight ones), but that would make the rendering considerably slower.

ggplot2::map_data() and rnaturalearth

The ggplot2 package has a function map_data() which wraps the maps package and provides shapes of countries and states ready for use in ggplot2 maps. However, if you play close attention you will see that there's an issue with Antarctica when using the WGS84 coordinate reference system (see bottom left).

world <- map_data(map = "world") %>%
  filter(region == "Antarctica")

ggplot() +
  geom_polygon(data = world, aes(x = long, y = lat, group = group), color ="#595959", fill ="#e5e5e5") +
  coord_sf(crs = 4326) +
  theme_void()

mapdata

An easy solution is to just switch to the rnaturalearth package:

world <- ne_countries(scale = "large", returnclass = "sf") %>%
  filter(admin == "Antarctica")

ggplot() +
  geom_sf(data = world) +
  coord_sf() +
  theme_void()

rnaturalearth

However, in contrast to the ggplot2::map_data() Antarctica shape, the one from rnaturalearth does not work well with an orthographic projection:

ggplot() +
  geom_sf(data = world) +
  coord_sf(crs = "+proj=ortho +lat_0=-90 +lon_0=0") +
  theme_void()

rnaturalearth ortho

The solution here is to apply a minimal buffer using st_buffer() after transforming the shape to the orthographic projection:

world_fixed <- world %>%
  st_transform("+proj=ortho +lat_0=-90 +lon_0=0") %>%
  st_buffer(1)

ggplot() +
  geom_sf(data = world_fixed) +
  coord_sf(crs = "+proj=ortho +lat_0=-90 +lon_0=0") +
  theme_void()

rnaturalearth ortho fixed