Fundamental misunderstanding with crs and map projections in R

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.

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