ggplot map displaying inaccurate coordinate scale

Question

I am trying to to plot a projected world map using ggplot. So far I have plotted said map, but the coordinate axis scale of the map is off by a factor of 10^5. I am unsure what interaction is causing this as the map data is in the correct scale according to the data in the geometry column. The following is my current code:

library(ggplot2)
library(rnaturalearth) #map data
library(tidyverse)
library(sf)
mapData = ne_countries(returnclass = 'sf')
#4087 is the EPSG code for the coordinate reference system used
ggplot(mapData)+geom_sf(aes(geometry = mapData$geometry))+coord_sf(crs = st_crs(4087))

Which plots the following map:

I initially thought that the map was generating at the wrong scale so I have tried shifting the position of the data into geom_sf() which had no effect on the graph. Removing geometry = mapData$geometry or the set crs value so the method will default generates a map with no coordinate labels. I would prefer to have the labels, and could not ascertain whether this fixed the issue, so this was not an option.

However, after using a vector of ordered points to add a scatter plot to the map I determined that the map was plotting at the right scale but the labels were displaying inaccurate values.

geoRand = data.frame(longitude = c(0:100, -50),latitude = c(0:100, -50))
ggplot(mapData)+geom_sf(aes(geometry = mapData$geometry))+coord_sf(crs = st_crs(4087))+
geom_point(data = GeoDist,mapping = aes(x = longitude, y = latitude))

I am unsure how to proceed from here as nothing I have attempted to change has had an apparent effect on this. I suspect it is some issue with the type/format of the data and what the method was expecting. However, this is pure speculation as to my understanding the sf package expects sf objects which the data is, and the geometry of the graph is determined by the geometry column in the data which is formatted as expected.

Answer 1

I think here the mistake is in geom_sf(aes(geometry = mapData$geometry)). You don't need to do that, since geom_sf() detects automagically the geometry column of your object.

My guess here is that you are passing mapData$geometry in the original CRS (probably EPSG:4326) but you are projecting to EPSG:4087, and that the internal mechanism of ggplot() is not capturing properly the CRS of mapData$geometry since you passed it as an additional object. For this reason the map is messed up (i.e. mapping EPSG:4326 on a EPSG:4087 space).

What you should do with sf object is just let geom_sf() does his thing with the geometry, and not use the aes(geometry) at all.

Additionally, with world maps it is a good idea to use coord_sf(expand = FALSE) to display properly the labels. See:

library(ggplot2)
library(rnaturalearth) # map data
library(tidyverse)
library(sf)
mapData <- ne_countries(returnclass = "sf")
# 4087 is the EPSG code for the coordinate reference system used

# Don't use the aes and coord_sf together!
ggplot(mapData) +
  geom_sf(aes(geometry = mapData$geometry)) +
  coord_sf(crs = st_crs(4087))

This is a replica of your map (notice the gap on Antarctica, left-bottom corner).

If we just remove the aes and use coord_sf(expand = FALSE):

# See, preserving the crs of the object
ggplot(mapData) +
  geom_sf(aes(geometry = mapData$geometry)) +
  coord_sf(expand = FALSE)

This is correct, but it is en EPSG 4326.

Putting all together

# See, also use expand = FALSE
ggplot(mapData) +
  geom_sf() +
  coord_sf(crs = st_crs(4087), expand = FALSE)

This is the right map and labels in EPSG:4087. Note that there are no gaps in the corner any more

Now we check but converting geoRand to an sf object with the right CRS as well:

# Check
geoRand <- data.frame(longitude = c(0:100, -50), latitude = c(0:100, -50)) %>%
  # As sf object  
  sf::st_as_sf(coords = c("latitude", "longitude"), crs = 4326)


ggplot(mapData) +
  geom_sf() +
  geom_sf(data = geoRand) +
  coord_sf(crs = st_crs(4087), expand = FALSE)

^{Created on 2024-02-21 with reprex v2.1.0}