Как нарисовать линию на 66-й северной параллели в ggplot при использовании полярной (epsg: 3995) проекции и ggspatial

#r #ggplot2 #gis #ggspatial

#r #ggplot2 #гис #ggspatial

Вопрос:

Я пытаюсь нарисовать синий круг на 66-й северной параллели. на карте, которую я делаю.

Я обычно fortify() использую SpatialPolygonsDataFrame, а затем использую комбинацию geom_polygon() coord_map("ortho" , ylim = c(50, 90)) и geom_hline() для достижения этого.

Однако для нового проекта мне придется использовать ggspatial::layer_spatial() объект, у которого уже есть дата CRS, поэтому рисование прямой линии и автоматическое перепроектирование не работают.

Это демонстрация моего обычного подхода, который работает:

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

введите описание изображения здесь

и я пытаюсь воспроизвести его в этом коде, но это, очевидно, не работает, потому geom_hline() что у него нет правильной проекции:

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

введите описание изображения здесь

Ответ №1:

Посмотрите, работает ли это для вас?

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

график