rggplot2r-sfmap-projectionsrnaturalearth
R, ggplot2, sf and rnaturalearth: an easy way to switch from one projection to another
First of all let me say that I am not a cartographer. I need to plot a map of Europe with some segments connecting the European capitals. Without going into the details, I provide a simplified example in the reprex at the end of the post. It is very easy to make a plot of the whole world and then change its projection (I do that at the beginning of the reprex).
It is also not difficult to select a window corresponding roughly to Europe and to plot a segment linking Rome to Paris as long as you use the default CRS (coordinate reference system).
Since I want to resort to the equal-earth projection, I did not find anything better than explicitly calculating the new coordinates of the window (and of Paris and Rome) in the new CRS. In doing so, I play with spatial objects from which I need to extract some numbers. It all feels rather cumbersome so I wonder if a real cartographer using R has a better way of achieving this and/or some functions which simplify the process.
Many thanks!
library(tidyverse)
library(rnaturalearth)
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
ww_ini <- ne_countries(scale = "medium",
type = 'map_units',
returnclass = "sf")
bb <- ne_download(type = "wgs84_bounding_box", category = "physical",
returnclass = "sf")
#> OGR data source with driver: ESRI Shapefile
#> Source: "/tmp/RtmpyDkGhC", layer: "ne_110m_wgs84_bounding_box"
#> with 1 features
#> It has 2 fields
gpl1 <- ggplot(data = ww_ini) +
geom_sf( col = "black", lwd = 0.3 )+
xlab(NULL) + ylab(NULL) +
ggtitle("Test title")+
geom_sf(data = bb, col = "grey", fill = "transparent") +
theme(plot.background = element_rect(fill = "white"),
panel.background = element_rect(fill = 'white'),
panel.grid.major = element_line(colour = "grey"),
legend.position="top",
plot.title = element_text(lineheight=.8, size=24, face="bold",
vjust=1),
legend.text = element_text(vjust=.4,lineheight=1,size = 14),
legend.title = element_text(vjust=1,lineheight=1, size=14,
face="bold" ))
ggsave("simple_world.pdf", gpl1, width=6*1.618,height=5)
gpl2 <- ggplot(data = ww_ini) +
geom_sf( col = "black", lwd = 0.3 )+
xlab(NULL) + ylab(NULL) +
ggtitle("Test title")+
geom_sf(data = bb, col = "grey", fill = "transparent") +
theme(plot.background = element_rect(fill = "white"),
panel.background = element_rect(fill = 'white'),
panel.grid.major = element_line(colour = "grey"),
legend.position="top",
plot.title = element_text(lineheight=.8, size=24, face="bold",
vjust=1),
legend.text = element_text(vjust=.4,lineheight=1,size = 14),
legend.title = element_text(vjust=1,lineheight=1, size=14,
face="bold" ))+
coord_sf( crs = "+proj=eqearth +wktext")
ggsave("simple_world_equal_earth.pdf", gpl2, width=6*1.618,height=5)
gpl3 <- ggplot(data = ww_ini) +
geom_sf( col = "black", lwd = 0.3 )+
xlab(NULL) + ylab(NULL) +
ggtitle("Test title")+
geom_sf(data = bb, col = "grey", fill = "transparent") +
theme(plot.background = element_rect(fill = "white"),
panel.background = element_rect(fill = 'white'),
panel.grid.major = element_line(colour = "grey"),
legend.position="top",
plot.title = element_text(lineheight=.8, size=24, face="bold",
vjust=1),
legend.text = element_text(vjust=.4,lineheight=1,size = 14),
legend.title = element_text(vjust=1,lineheight=1, size=14,
face="bold" ))+
coord_sf(xlim=c(-20,45), ylim=c(30, 73) )
ggsave("simple_europe.pdf", gpl3, width=6*1.618,height=5)
### now I try making a plot of Europe using the equal earth projection.
## I need to recalculate the coordinate for the zoom in the new equal-earth
## CRS (coordinate reference system)
## See https://www.r-bloggers.com/2019/04/zooming-in-on-maps-with-sf-and-ggplot2/
## Step 1: determine the target coordinate reference system
target_crs <- '+proj=eqearth +wktext'
## Step 2: I transform the world data into the new coordinate system
worldmap_trans <- st_transform(ww_ini, crs = target_crs)
### Step 3: define the display window as a latitude and longitude interval
disp_win_wgs84 <- st_sfc(st_point(c(-20, 30)), st_point(c(45, 73)),
crs = 4326)
## Step 4: transform the display window in the new coordinate system.
disp_win_trans <- st_transform(disp_win_wgs84, crs = target_crs)
#### Step 5: retrieve the window coordinates in the new coordinate system
disp_win_coord <- st_coordinates(disp_win_trans)
gpl4 <- ggplot(data = worldmap_trans) +
geom_sf( col = "black", lwd = 0.3 )+
xlab(NULL) + ylab(NULL) +
ggtitle("Test title")+
geom_sf(data = bb, col = "grey", fill = "transparent") +
theme(plot.background = element_rect(fill = "white"),
panel.background = element_rect(fill = 'white'),
panel.grid.major = element_line(colour = "grey"),
legend.position="top",
plot.title = element_text(lineheight=.8, size=24, face="bold",
vjust=1),
legend.text = element_text(vjust=.4,lineheight=1,size = 14),
legend.title = element_text(vjust=1,lineheight=1, size=14,
face="bold" ))+
coord_sf(xlim = disp_win_coord[,'X'], ylim = disp_win_coord[,'Y'],
datum = target_crs## , expand = FALSE
)
ggsave("simple_europe_equal_earth.pdf", gpl4, width=6*1.618,height=5)
####Now I want to add a segment connecting Rome to Paris to the Europe map.
### easy with the default CRS.
cities <- structure(list(city_ascii = c("Rome", "Paris"), country = c("Italy",
"France"), lat = c(41.8931, 48.8566), lng = c(12.4828, 2.3522
)), row.names = c(NA, -2L), class = c("tbl_df", "tbl", "data.frame"
))
gpl3bis <- ggplot(data = ww_ini) +
geom_sf( col = "black", lwd = 0.3 )+
xlab(NULL) + ylab(NULL) +
ggtitle("Test title")+
geom_sf(data = bb, col = "grey", fill = "transparent") +
theme(plot.background = element_rect(fill = "white"),
panel.background = element_rect(fill = 'white'),
panel.grid.major = element_line(colour = "grey"),
legend.position="top",
plot.title = element_text(lineheight=.8, size=24, face="bold",
vjust=1),
legend.text = element_text(vjust=.4,lineheight=1,size = 14),
legend.title = element_text(vjust=1,lineheight=1, size=14,
face="bold" ))+
geom_segment(aes(x = cities$lng[1], y =cities$lat[1] ,
xend = cities$lng[2], yend = cities$lat[2]), colour = "red")+
coord_sf(xlim=c(-20,45), ylim=c(30, 73) )
ggsave("simple_europe_segment.pdf", gpl3bis, width=6*1.618,height=5)
#####Now I try to do the same using the equal-earth projection
## I need to convert the coordinates of Rome and Paris into those of the
## equal earth CRS.
point_mat <- cities %>%
select(lng, lat) %>%
as.matrix
pp <- st_multipoint(point_mat)
pp_equal_earth <- st_sfc(pp,
crs = 4326)
points_equal_earth <- st_transform(pp_equal_earth, crs = target_crs)
##Now I need a way to extract the number from points_equal_earth.
## I did not find anything better than
points_simple <- as_Spatial(points_equal_earth)%>%as_tibble
gpl4bis <- ggplot(data = worldmap_trans) +
geom_sf( col = "black", lwd = 0.3 )+
xlab(NULL) + ylab(NULL) +
ggtitle("Test title")+
geom_sf(data = bb, col = "grey", fill = "transparent") +
theme(plot.background = element_rect(fill = "white"),
panel.background = element_rect(fill = 'white'),
panel.grid.major = element_line(colour = "grey"),
legend.position="top",
plot.title = element_text(lineheight=.8, size=24, face="bold",
vjust=1),
legend.text = element_text(vjust=.4,lineheight=1,size = 14),
legend.title = element_text(vjust=1,lineheight=1, size=14,
face="bold" ))+
geom_segment(aes(x = points_simple$X1[1], y =points_simple$X2[1] ,
xend = points_simple$X1[2], yend = points_simple$X2[2]), colour = "red")+
coord_sf(xlim = disp_win_coord[,'X'], ylim = disp_win_coord[,'Y'],
datum = target_crs## , expand = FALSE
)
ggsave("simple_europe_equal_earth_segment.pdf", gpl4bis, width=6*1.618,height=5)
Created on 2021-08-04 by the reprex package (v2.0.0)
Solution
My advice in your case is to transform both the points of the cities and the connecting line to spatial POINT and LINESTRING. Once that you have that, you don't need to worry about extracting numbers on different CRS, st_transform would do it for you. Also you can plot everything with ggplot2/geom_sf rather than geom_segment, as you would have everything on sf format. I think the geom_segment was the conflicting part on your code.
See a reprex, although I have to admit I am not sure I fully understood the question. Hope this helps:
library(tidyverse)
library(rnaturalearth)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
# Get countries
ww_ini <- ne_countries(
scale = "medium",
type = "map_units",
returnclass = "sf"
)
# Cities
cities <-
structure(
list(
city_ascii = c("Rome", "Paris"),
country = c(
"Italy",
"France"
),
lat = c(41.8931, 48.8566),
lng = c(12.4828, 2.3522)
),
row.names = c(NA, -2L),
class = c("tbl_df", "tbl", "data.frame")
)
# To spatial points
cities_sf <- st_as_sf(cities, coords = c("lng", "lat"), crs = 4326)
# Key point! Create an spatial line connecting the points
line <-
st_linestring(st_coordinates(cities_sf)) %>% st_sfc(crs = 4326)
# Your display window
disp_win_wgs84 <- st_sfc(st_point(c(-20, 30)),
st_point(c(45, 73)),
crs = 4326
)
# Ensure everything is in the same proj
# Clasic error when working with maps
target_crs <- "+proj=eqearth +wktext"
ww_ini_eqea <- st_transform(ww_ini, target_crs)
cities_sf_eqea <- st_transform(cities_sf, target_crs)
line_eqea <- st_transform(line, target_crs)
disp_win_eqea <- st_transform(disp_win_wgs84, target_crs)
st_coordinates(disp_win_eqea)
#> X Y
#> 1 -1793068 3764325
#> 2 2837546 7872578
# Map everything
# Start with world
ggplot(ww_ini_eqea) +
geom_sf() +
# Add cities
geom_sf(data = cities_sf_eqea, color = "blue") +
# Now the line
geom_sf(data = line_eqea, color = "red") +
# Limit the map
coord_sf(
xlim = st_coordinates(disp_win_eqea)[, "X"],
ylim = st_coordinates(disp_win_eqea)[, "Y"]
) +
theme_minimal()
sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19043)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252
#> [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
#> [5] LC_TIME=Spanish_Spain.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] sf_1.0-2 rnaturalearth_0.1.0 forcats_0.5.1
#> [4] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
#> [7] readr_2.0.0 tidyr_1.1.3 tibble_3.1.3
#> [10] ggplot2_3.3.5 tidyverse_1.3.1
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.7 lattice_0.20-44 lubridate_1.7.10
#> [4] class_7.3-19 assertthat_0.2.1 digest_0.6.27
#> [7] utf8_1.2.2 R6_2.5.0 cellranger_1.1.0
#> [10] backports_1.2.1 reprex_2.0.0 rnaturalearthdata_0.1.0
#> [13] evaluate_0.14 e1071_1.7-8 httr_1.4.2
#> [16] highr_0.9 pillar_1.6.2 rlang_0.4.11
#> [19] readxl_1.3.1 rstudioapi_0.13 rmarkdown_2.9
#> [22] styler_1.5.1 munsell_0.5.0 proxy_0.4-26
#> [25] broom_0.7.9 compiler_4.1.0 modelr_0.1.8
#> [28] xfun_0.24 pkgconfig_2.0.3 rgeos_0.5-5
#> [31] htmltools_0.5.1.1 tidyselect_1.1.1 fansi_0.5.0
#> [34] crayon_1.4.1 tzdb_0.1.2 dbplyr_2.1.1
#> [37] withr_2.4.2 grid_4.1.0 jsonlite_1.7.2
#> [40] gtable_0.3.0 lifecycle_1.0.0 DBI_1.1.1
#> [43] magrittr_2.0.1 units_0.7-2 scales_1.1.1
#> [46] KernSmooth_2.23-20 cli_3.0.1 stringi_1.7.3
#> [49] farver_2.1.0 fs_1.5.0 sp_1.4-5
#> [52] xml2_1.3.2 ellipsis_0.3.2 generics_0.1.0
#> [55] vctrs_0.3.8 tools_4.1.0 glue_1.4.2
#> [58] hms_1.1.0 yaml_2.2.1 colorspace_2.0-2
#> [61] classInt_0.4-3 rvest_1.0.1 knitr_1.33
#> [64] haven_2.4.3
Created on 2021-08-05 by the reprex package (v2.0.0)
- Arranging bars in ggplot
- How to smooth an arc
- ggplot of multiple geom_sf datasets gives mismatched legend colours
- I get an error when using {future} and {furrr} functions within a Golem Shiny App, what does it come from?
- ggpredict and quartiles: predictions with only one data point?
- How to properly create secondary y-axis with percentage scale and multiple lines in R?
- ipwtm function in R: First argument a vector
- How can I download GADM data in R?
- Use dplyr to compute "streaks" in column of data
- Count empty strings?
- How to show a ConditionalPanel when a row in a reactable table is selected?
- How to lookup and sum multiple columns in R
- Is there a way to make the bars' colors conditional to the column?
- How to assign values to a df$column from another df using conditions in R
- Is there a way to read Qlikview data object using R and perform an statistical operation and return results to Qlikview?
- R: Extracting characters from string vectors to data frame rows
- Else if statement to create a new column
- Calculate percentage in a data frame R
- How can I calculate sums of multiple columns without them being summed together into 1 column?
- How to convert a binary file with extension .GRD to .TXT in R?
- How to specify the size of a graph in ggplot2 independent of axis labels
- Using sqldf filter and group by in r dynamically
- Grouping data and reassembling rows as columns
- Coordinates of current mouse position on leaflet map with shiny
- read_csv drops 'f' from end of string
- wes_palettes: How to obtain color hex codes from a wesanderson palette when it set to continuous with more than the predefined colors
- ggplot2 change axis limits for each individual facet panel
- Mixing Rcpp modules and Rcpp::export
- How to create a forest plot that uses estimates obtained from ClubSandwich?
- how to include an abstract in a Rmd file