How to make geom_sf() point northwards?

Question

I am trying to make a mapping package for R and ggplot that would allow the users to define their data using decimal degrees but would use appropriate projection depending on where the user's data are located. I decided to use the ggspatial package for this purpose since it automatically converts the coordinates. The ggspatial package, in turn, uses the sf, package, which is supported by ggplot.

I use Arctic polar stereographic projection for data originating >60 latitude. The projection is cartesian with meters from the North Pole as coordinate units. Naturally, when I try to plot this projection on the Pacific side of the Arctic, North will be pointing down (because of y-coordinates point towards the North Pole).

How do I turn the maps to point towards North using coord_sf and ggplot? If there is no reasonable way of doing this, would you be able to suggest a way around this problem?

Example data:

library(rnaturalearth)
library(sp)
library(raster)
library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(ggplot2)
library(ggspatial)

# Get coastline data a example

dt <- rnaturalearth::ne_coastline()

# Clip and reproject the data to Arctic polar stereographic

clip_boundary <- sp::SpatialPolygons(
  list(sp::Polygons(
    list(sp::Polygon(
      data.frame(lon = c(-180, 180, 180, -180), lat = c(60, 60, 90, 90)))), ID = 1)
  ), proj4string = sp::CRS(sp::proj4string(dt)))

arctic <- raster::crop(dt, clip_boundary)

arctic <- sp::spTransform(arctic, sp::CRS("+init=epsg:3995"))
arctic_sf <- sf::st_as_sf(arctic)

# Define the map limits

limits_dd <- c(-160, 160, 60, 80)
limits_dd_shp <- sp::Polygon(
  matrix(c(limits_dd[1], limits_dd[3], limits_dd[1], limits_dd[4], limits_dd[2], 
           limits_dd[4], limits_dd[2], limits_dd[3], limits_dd[1], limits_dd[3]), 
         ncol = 2, byrow = TRUE))
limits_dd_shp <- sp::SpatialPolygons(
  list(sp::Polygons(list(limits_dd_shp), ID = "clip_boundary")), 
  proj4string= sp::CRS("+init=epsg:4326"))
limits_utm_shp <- sp::spTransform(limits_dd_shp, sp::CRS("+init=epsg:3995"))
limits_utm <- unname(c(limits_utm_shp@bbox[1,], limits_utm_shp@bbox[2,]))

ggplot() + 
  ggspatial::layer_spatial(data = arctic_sf) + 
  ggspatial::layer_spatial(data = limits_utm_shp, color = "red", fill = NA)

Map showing the example data. The red polygon indicates the region which should be plotted.

 # Decimal degree north is "wrong way around" because underlying projection defines so

    ggplot() + 
      ggspatial::layer_spatial(data = arctic_sf) + 
      coord_sf(xlim = limits_utm[1:2], ylim = limits_utm[3:4])

    # A dirty trick trying to invert the y-axis does not help

    ggplot() + 
      ggspatial::layer_spatial(data = arctic_sf) + 
      coord_sf(xlim = limits_utm[1:2], ylim = limits_utm[c(4,3)])
    #> Error in st_normalize.sfc(x, c(x_range[1], y_range[1], x_range[2], y_range[2])): domain must have a positive range

    # Neither does scale_y_reverse

    ggplot() + 
      ggspatial::layer_spatial(data = arctic_sf) + 
      coord_sf(xlim = limits_utm[1:2], ylim = limits_utm[3:4]) + 
      scale_y_reverse()
    #> Error in st_normalize.sfc(x, c(x_range[1], y_range[1], x_range[2], y_range[2])): domain must have a positive range

The only way I have achieved this is to use the cartesian coordinates directly, but this approach would not allow plotting data on top using decimal degrees (geom_spatial_point) later.

ggplot(data = arctic) + 
  geom_path(aes(x = long, y = lat, group = group)) + 
  coord_cartesian(xlim = limits_utm[1:2], ylim = limits_utm[c(4,3)])

^{Created on 2020-04-17 by the reprex package (v0.3.0)}