Displaying degrees on map axes using coord_sf in R

I am making a map and wanted to use an alternative projection. When I do this using coord_sf, the x- and y-axes get labels of some other units like 0.001°W. I'm not a geographer, so I'm not sure what these are. I'd like to just see degrees latitude and longitude...

Code to generate the basic map:

mapWorld <- borders("world", colour="gray70", fill="gray70") 
ggplot() + mapWorld +  coord_sf(xlim=c(-125, 155),ylim=c(-22,40), crs = "+proj=laea") 

Resulting plot:

The simplest way of fixing your problem here is to simply remove the reprojection. The reason you are seeing funny coordinates is that you supplied a different coordinate reference system through coord_sf(...).

#install.packages('maps')
mapWorld <- borders("world",colour="gray70", fill="gray70") 
ggplot() + mapWorld + coord_sf(xlim=c(-125, 155),ylim=c(-22,40))

The more complicated answer relates to whether you want to use a different projection system. Where you have coord_sf(xlim=c(-125, 155),ylim=c(-22,40), crs = "+proj=laea") you are applying a modification to geom_sf to use a different projection of the map. Your code laea I assume refers to the Lambert Azimuthal Equal Area projection centred on Europe, right? That one has CRS code 3035.

Where you will run into a problem, though, is that in order to apply a spatial projection, you need to be working with a spatial object, not a picture of a map. When you call: borders("world",colour="gray70", fill="gray70") you are using a quick hack in ggplot2 to get simple x and y values from the data in the maps package. The problem is that borders() doesn't take projection as an argument when it retrieves data from the maps package, so it only takes simple x(long) and y(lat) values and turns it into a dataframe.

Now, you can reproject and plot a map using just the maps (and mapproj) packages. If you check out the documentation on maps::map it shows you can supply a projection = ... argument, but the projections are limited to the ones included in mapproj. For example, you could use an Azimuth Equal Area projection:

map("world", projection = "azequalarea", xlim = c(-125,155), ylim = c(-22,40))

And it would look like this, which is kind of ugly. You can modify how it is centred, etc, but it still won't look that great. Which shouldn't be that surprising as you're using a library that predates R itself, coming from the language S.

You could approach this using more modern tools that give you a lot more flexibility. geom_sf in ggplot is designed to work with spatial objects, ideally in sf format, but it is easy to convert between different kinds of spatial objects. Here's a quick example how to make the map you have in mind.

First, let's get a map of the world as a spatial data object

#install.packages('rnaturalearth')
#install.packages('rnaturalearthdata')
#install.packages('sf')
#install.packages('rgeos')

library(rnaturalearth)
library(rnaturalearthdata)
library(sf)

world <- ne_countries(scale = "medium", returnclass = "sf")
> world
Simple feature collection with 241 features and 63 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
First 10 features:

This object has a projection attached CRS:+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0. Which is essentially the default projection when working with latitude and longitude, and what most people are used to. We can plot this:

ggplot(data = world) +
  geom_sf() + 
  xlab("Longitude") + ylab("Latitude")

Now that it is a spatial data object, it's easy to change projections or to customize it based on elements within the data. For example, using the LAEA projection code (3035) and adding colours. for different continents:

ggplot(data = world) +
  geom_sf(aes(fill = continent)) + 
  coord_sf(crs = 3035)

If you still want to have boundaries like xlim=c(-125, 155),ylim=c(-22,40), you need to be mindful about what those numbers represent when dealing with other projections. Those numbers refer to lat/lon values, but when we're using other projections we have to use the units those projections are in. You would have to reproject those limits into coordinates used by that specific projection; however, based on the limits you've set, I would suggest they will not work well at all with this type of projection. I would advise you to stick with lat/long, i.e. CRS = 4326 if you want your map to look like what you had in mind.

ggplot(data = world) +
  geom_sf() + 
  xlab("Longitude") + ylab("Latitude") + coord_sf(xlim=c(-125, 155),ylim=c(-22,40))