J'essaie de dessiner un cercle bleu sur le 66e parallèle nord. sur une carte que je fais.

J'ai généralement fortify() un SpatialPolygonsDataFrame, puis j'utilise la combinaison geom_polygon(), coord_map("ortho" , ylim = c(50, 90)) et geom_hline() pour y parvenir.

Cependant, pour un nouveau projet, je devrai utiliser ggspatial::layer_spatial() sur un objet qui a déjà une date CRS, donc ce dessin en ligne droite et cette re-projection automatique ne fonctionnent pas.

Ceci est une démonstration de mon approche habituelle, qui fonctionne :

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

Et j'essaie de le reproduire dans ce code, mais cela ne fonctionne évidemment pas car le geom_hline() n'a pas la bonne 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

0
L Smeets 18 oct. 2020 à 20:43

1 réponse

Meilleure réponse

Voir si cela fonctionne pour vous ?

# 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

2
Z.Lin 19 oct. 2020 à 07:17