rggplot2gisggspatial

How to draw a line on the 66th parallel north in ggplot when using a polar (epsg:3995) projection and ggspatial

I am trying to draw a blue circle on the 66th parallel north. on a map I am making.

I usually fortify() a SpatialPolygonsDataFrame and then use the combination geom_polygon(), coord_map("ortho" , ylim = c(50, 90)) and geom_hline() to achieve this.

However, for a new project I will have to use ggspatial::layer_spatial() on an object that already has a CRS date, so this straight line drawing and automatic re-projection does not work.

This is a demonstration of my usual approach, which works:

library(tidyverse)

map_data <- rnaturalearth::ne_countries()
map_data_df <- fortify(map_data)

ggplot() + 
  geom_polygon(data = map_data_df, 
               aes(x = long,
                   y = lat,
                   group = group),
               col = "black",
               fill = "lightgrey") + 
  coord_map("ortho" , ylim = c(50, 90))  +
  geom_hline(yintercept = 66, col = "blue")

enter image description here

and I am trying to replicate it in this code, but that obviously does not work because the geom_hline() does not have the right projection:

map_data_transformed <- sp::spTransform(map_data, sp::CRS("+init=epsg:3995"))
# needs some fixing to close all the polygons
map_data_transformed <- rgeos::gBuffer(map_data_transformed, byid=TRUE, width=0)


ggplot() + 
  ggspatial::layer_spatial(data = map_data_transformed, 
                           fill = "lightgrey", col = "black", size = .2)+
  coord_sf( xlim = c(-4500000 , 4500000),
            ylim = c(-4500000 , 4500000)) +
  geom_hline(yintercept = 66, col = "blue")

enter image description here

Solution

  • See if this works for you?

    # define parallel north line as a spatial polygon with latitude / longitude projection
new_line <- sp::Polygon(coords = matrix(c(seq(-180, 180),
                                          rep(66, times = 361)), # change this for other
                                                                 # intercept values
                                        ncol = 2))
new_line <- sp::SpatialPolygons(Srl = list(sp::Polygons(srl = list(new_line), 
                                                        ID = "new.line")),
                                proj4string = sp::CRS("+proj=longlat +datum=WGS84"))

# change projection to match that of data
new_line <- sp::spTransform(new_line, sp::CRS("+init=epsg:3995"))

# plot
ggplot() + 
  ggspatial::layer_spatial(data = map_data_transformed, 
                           fill = "lightgrey", col = "black", size = .2) +
  ggspatial::layer_spatial(data = new_line,
                           fill = NA, col = "blue") +
  coord_sf( xlim = c(-4500000 , 4500000),
            ylim = c(-4500000 , 4500000))

    plot