Fundamental misunderstanding with crs and map projections in R

Question

I think I have a fundamental misunderstanding with how using crs and map projections work in R. First let me show what my final goal is here. I would like to create something like the following: Using rnaturalearth. Specifically I am interested in the curved longitude lines and the diagonal latitude lines.

So far, I have used

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

gl <- ne_countries(scale = "large", country = "Greenland",returnclass = "sf")
ggplot() + 
  geom_sf(data = gl)

to create the following

The part I am struggling with is indicating these latitude and longitude lines. I know this has something to do with how I project Earth onto the 2D plane (by changing the crs) but I am struggling to tie it all together in R. I believe I have a fundamental misunderstanding of how all this works and could do with some explanation.

I should note that my final goal is to indicate (through the use of dashed lines) areas bounded by longitudes 60,62,64,66 and 68.

Allan Cameron · Accepted Answer · 2023-06-08T14:31:31.383

The CRS used in your picture looks like EPSG:3184, which is the appropriate longitudinal slice of the Universal Transverse Mercator system. As long as you tell coord_sf this is the co-ordinate reference system that you want to replicate, it will automatically convert your map to that system for plotting. You can do this as follows:

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

north <- 72.5; south <- 59; east <- -52; west <- -45

extent <- st_sf(a = 1:2, crs = 4326,
            geom = st_sfc(st_point(c(east, south)), 
                          st_point(c(west, north)))) |>
            st_transform(crs = 'EPSG:3184') |>
            st_bbox()

lat <- lapply(c(60, 62, 64, 66, 68), 
            function(x) cbind(seq(-70, -40, 0.1), x)) |>
  st_multilinestring() |> st_sfc(crs = 4326) |> st_sf() |>
  st_transform(crs = 'EPSG:3184') 

lon <- lapply(c(-60, -55, -50, -45, -40), 
              function(x) cbind(x, seq(56, 74, 0.2))) |>
  st_multilinestring() |> st_sfc(crs = 4326) |> st_sf() |>
  st_transform(crs = 'EPSG:3184') 
  
ne_countries(scale = "large", country = "Greenland",returnclass = "sf") |>
  ggplot() + 
  geom_sf(fill = 'white') +
  geom_sf(data = lat, linetype = 2, color = 'gray50') +
  geom_sf(data = lon, color = 'gray', linewidth = 0.2) +
  coord_sf(crs = 'EPSG:3184', xlim = extent[c(1, 3)], ylim = extent[c(2, 4)]) +
  scale_x_continuous(breaks = -5 * 1:10) +
  theme_bw() +
  theme(panel.background = element_rect(fill = '#f4f6f5'),
        panel.grid = element_blank())

On my original projection, using coord_sf(xlim = c(-61.2, -45), ylim = c(59, 72.5)) would crop the South-West coast. On the new projection this is not the case. Could you please explain how the coordinates for the SW coast are projected into EPSG:3184? — Fish_Person, Jun 08 '23 at 10:01
@Fish_Person the co-ordinates in EPSG:3184 are given in terms of metres of Northing and Easting from the equator at Greenwich meridian. Because you want diagonal longitude lines and curved latitude lines, setting degrees of latitude and longitude as a bounding box won't produce a rectangular final plot. However, you can easily convert between the WGS84 (lat/long) co-ordinates and UTM, so I'll update my post to show a cropped version. — Allan Cameron, Jun 08 '23 at 10:14
And if I was then to try and add black dashed lines following the longitudinal lines starting at latitudes 60,62,64,66,68? I have attempted to use geom_curve() to no avail. Thank you again. — Fish_Person, Jun 08 '23 at 11:11
@Fish_Person I have changed my answer to incorporate your request — Allan Cameron, Jun 08 '23 at 14:32

Fundamental misunderstanding with crs and map projections in R

1 Answers1