Search code examples
rggplot2mapsrnaturalearth

R+rnaturalearth: Creating a Facet Plot and Removing an Inner Border


Please have a look at the reprex at the end of the plot. My goal (which I achieve step by step) is to construct a panel plot combining different color maps of European countries while removing in the map the inner border in Cyprus (no political statement on my side, I have just to do it). I somehow achieve this, but in a rather cumbersome way, so I wonder if someone can suggest some improvements.

Thanks!

library(tidyverse)
library(rnaturalearth)
#> Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE

## data I want to plot in a facet plot.

df_plot <- structure(list(expenditure_year = c(2019, 2019, 2019, 2019, 2019, 
2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 
2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 
2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 
2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 
2020, 2020, 2020, 2020, 2020, 2021, 2021, 2021, 2021, 2021, 2021, 
2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 
2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2022, 
2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 
2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 
2022, 2022, 2022, 2022), iso2 = c("IT", "FI", "FR", "NL", "CY", 
"BE", "SE", "GR", "DE", "LT", "BG", "PT", "CZ", "IE", "EE", "HU", 
"SK", "ES", "PL", "SI", "DK", "AT", "RO", "LU", "MT", "LV", "HR", 
"IT", "FI", "FR", "NL", "CY", "BE", "SE", "GR", "DE", "LT", "BG", 
"PT", "CZ", "IE", "EE", "HU", "SK", "ES", "PL", "SI", "DK", "AT", 
"RO", "LU", "MT", "LV", "HR", "IT", "FI", "FR", "NL", "CY", "BE", 
"SE", "GR", "DE", "LT", "BG", "PT", "CZ", "IE", "EE", "HU", "SK", 
"ES", "PL", "SI", "DK", "AT", "RO", "LU", "MT", "LV", "HR", "IT", 
"FI", "FR", "NL", "CY", "BE", "SE", "GR", "DE", "LT", "BG", "PT", 
"CZ", "IE", "EE", "HU", "SK", "ES", "PL", "SI", "DK", "AT", "RO", 
"LU", "MT", "LV", "HR"), percentage_share = c(0.450577656119157, 
0.92248580410076, 0.975584667105617, 0.409803641820049, 0.539818102589104, 
0.96123699068389, 0.826519143748971, 0.570729391308522, 1.51184843024544, 
1.56004652853805, 0.682394691276564, 0.529918189934815, 1.51537052534534, 
0.340282452391896, 1.23411076526779, 1.78839893466253, 0.648651854236538, 
0.474837572951868, 1.11702742041526, 0.858524606698324, 1.43464767790405, 
0.505813688224417, 0.655039437261184, 0.310317708208196, 1.6984935879487, 
1.07006270258955, 1.33269223306193, 2.02652222755559, 1.52387828834052, 
2.27427055972995, 0.875130252470089, 1.50963072632275, 1.43394831922113, 
1.02805626070114, 4.31903491607162, 3.71236026065522, 2.00271869460953, 
1.43822454206211, 1.63660567657213, 2.74250293931477, 0.407182979009171, 
1.69199343784178, 3.12677589492696, 1.76946859138244, 0.88979342454491, 
5.16611685495998, 3.54047898922094, 2.48037431884307, 1.99660745861322, 
1.23774421769849, 0.879255567282404, 4.24216627476823, 2.70305053222405, 
2.66814646113911, 2.11099250992334, 1.59620461257039, 2.69079496650438, 
1.50012485828527, 2.13084813620244, 1.29750144508215, 1.04199044469184, 
3.03512300799337, 3.17970791856142, 1.58793142828813, 1.68604027084116, 
1.42812640590373, 3.62657019057752, 0.721667603152213, 1.41210754275081, 
3.61800356104344, 2.22916452630624, 1.11319571460128, 1.81329721467997, 
3.00851645408846, 1.56776732632244, 2.46506682189924, 1.58129951074308, 
0.813434327102068, 3.99050063429371, 2.81167234901301, 2.20079792507412, 
1.41830583744773, 0.944648699957193, 1.71235447267469, 1.03518425244823, 
0.538183029124816, 1.02736324514123, 0.893949584792747, 1.37693432013489, 
2.22128961439947, 1.13396128434949, 1.01149050905179, 0.950283072433532, 
1.43451796238923, 0.487096865306793, 0.92012157362591, 2.08815527711486, 
0.884735492296972, 1.27188526690518, 0.927269440893475, 0.84679781968768, 
1.39310961513098, 1.45742166676803, 1.46526610855142, 0.665295824788144, 
1.50079360929557, 1.00760277849241, 0.839755449738563), iso3 = c("ITA", 
"FIN", "FRA", "NLD", "CYP", "BEL", "SWE", "GRC", "DEU", "LTU", 
"BGR", "PRT", "CZE", "IRL", "EST", "HUN", "SVK", "ESP", "POL", 
"SVN", "DNK", "AUT", "ROM", "LUX", "MLT", "LVA", "HRV", "ITA", 
"FIN", "FRA", "NLD", "CYP", "BEL", "SWE", "GRC", "DEU", "LTU", 
"BGR", "PRT", "CZE", "IRL", "EST", "HUN", "SVK", "ESP", "POL", 
"SVN", "DNK", "AUT", "ROM", "LUX", "MLT", "LVA", "HRV", "ITA", 
"FIN", "FRA", "NLD", "CYP", "BEL", "SWE", "GRC", "DEU", "LTU", 
"BGR", "PRT", "CZE", "IRL", "EST", "HUN", "SVK", "ESP", "POL", 
"SVN", "DNK", "AUT", "ROM", "LUX", "MLT", "LVA", "HRV", "ITA", 
"FIN", "FRA", "NLD", "CYP", "BEL", "SWE", "GRC", "DEU", "LTU", 
"BGR", "PRT", "CZE", "IRL", "EST", "HUN", "SVK", "ESP", "POL", 
"SVN", "DNK", "AUT", "ROM", "LUX", "MLT", "LVA", "HRV"), country = c("Italy", 
"Finland", "France", "Netherlands", "Cyprus", "Belgium", "Sweden", 
"Greece", "Germany", "Lithuania", "Bulgaria", "Portugal", "Czechia", 
"Ireland", "Estonia", "Hungary", "Slovakia", "Spain", "Poland", 
"Slovenia", "Denmark", "Austria", "Romania", "Luxembourg", "Malta", 
"Latvia", "Croatia", "Italy", "Finland", "France", "Netherlands", 
"Cyprus", "Belgium", "Sweden", "Greece", "Germany", "Lithuania", 
"Bulgaria", "Portugal", "Czechia", "Ireland", "Estonia", "Hungary", 
"Slovakia", "Spain", "Poland", "Slovenia", "Denmark", "Austria", 
"Romania", "Luxembourg", "Malta", "Latvia", "Croatia", "Italy", 
"Finland", "France", "Netherlands", "Cyprus", "Belgium", "Sweden", 
"Greece", "Germany", "Lithuania", "Bulgaria", "Portugal", "Czechia", 
"Ireland", "Estonia", "Hungary", "Slovakia", "Spain", "Poland", 
"Slovenia", "Denmark", "Austria", "Romania", "Luxembourg", "Malta", 
"Latvia", "Croatia", "Italy", "Finland", "France", "Netherlands", 
"Cyprus", "Belgium", "Sweden", "Greece", "Germany", "Lithuania", 
"Bulgaria", "Portugal", "Czechia", "Ireland", "Estonia", "Hungary", 
"Slovakia", "Spain", "Poland", "Slovenia", "Denmark", "Austria", 
"Romania", "Luxembourg", "Malta", "Latvia", "Croatia")), row.names = c(NA, 
-108L), class = c("tbl_df", "tbl", "data.frame"))



ww_ini <- ne_countries(scale = "medium",
                       type = 'countries',
                       returnclass = "sf")

bb <- ne_download(type = "wgs84_bounding_box", category = "physical",
                  returnclass = "sf") 


### select a subset of the data for a single year to plot

df_plot_sel <- df_plot |>
    filter(expenditure_year==2022)


ww_plot <- ww_ini |>
    select(adm0_a3) |>
    left_join(y=df_plot_sel, by=c("adm0_a3"="iso3"))


## start with a map plot without colors

gpl <-   ggplot(data = ww_plot) +
        geom_sf( aes(fill=percentage_share),  col = "black", lwd = 0.3 )+

        coord_sf(xlim=c(-20,45), ylim=c(30, 73) ) +
  theme_minimal()


gpl


#### The above is OK, but I need to remove the boundary inside Cyprus (no political statement whatsoever from my side).

## See https://stackoverflow.com/questions/68672575/r-rnaturalearth-and-sf-remove-a-single-border-from-map

 cyprus <- ww_ini |> 
  select(adm0_a3)  |> 
  filter(adm0_a3 %in% c("CYP", "CYN")) |>  
  mutate(adm0_a3 = "CYP")  |> ## I renamed as Cyprus my version
  group_by(adm0_a3)  |> 
  summarise()



ww_ini_cyp <- ww_ini  |>  
  select(adm0_a3)  |> 
  filter(!adm0_a3 %in% c("CYP", "CYN"))  |>  
  bind_rows(cyprus)



gpl_cyp <-   ggplot(data = ww_ini_cyp) +
        geom_sf( col = "black", lwd = 0.3 )+

        coord_sf(xlim=c(-20,45), ylim=c(30, 73) ) +
  theme_minimal()

### and this is OK. I now have a map without the inner border in Cyprus
gpl_cyp


### Now let us replot a colormap with Cyprus without border




ww_plot_cyp <- ww_ini_cyp |>
    ## select(adm0_a3) |>
    left_join(y=df_plot_sel, by=c("adm0_a3"="iso3"))

gpl_cyp <-   ggplot(data = ww_plot_cyp) +
        geom_sf( aes(fill=percentage_share),  col = "black", lwd = 0.3 )+

        coord_sf(xlim=c(-20,45), ylim=c(30, 73) ) +
  theme_minimal()


gpl_cyp   ## and this works


## Now let us try a facet plot



ww_plot_cyp_facet <- ww_ini_cyp |>
    left_join(y=df_plot, by=c("adm0_a3"="iso3")) 

gpl_cyp_facet <-   ggplot(data = ww_plot_cyp_facet) +
        geom_sf( aes(fill=percentage_share),  col = "black", lwd = 0.3 )+
facet_wrap( ~expenditure_year, nrow = 2, scales = "fixed")+
        coord_sf(xlim=c(-20,45), ylim=c(30, 73) ) +
  theme_minimal()


gpl_cyp_facet   ## and this is a disaster. How do I fix it?

## The problem is that the expenditure_year is missing for adm0_a3 codes
## I just need to create it for the all the adm0_a3 codes






years <- df_plot|>
    pull(expenditure_year) |>
    unique()

iso_list <- ww_ini$adm0_a3 |>
    unique()

df <- crossing(years, iso_list)


ww_plot_cyp_facet_fix <- ww_ini_cyp |>
    left_join(df, by=c("adm0_a3"="iso_list")) |>
    left_join(df_plot, by=c("adm0_a3"="iso3",
                            "years"="expenditure_year"))



gpl_cyp_facet_fix <-   ggplot(data = ww_plot_cyp_facet_fix) +
        geom_sf( aes(fill=percentage_share),  col = "black", lwd = 0.3 )+
facet_wrap( ~years, nrow = 2, scales = "fixed")+
        coord_sf(xlim=c(-20,45), ylim=c(30, 73) ) +
  theme_minimal()


gpl_cyp_facet_fix   ## and this works, but if feels very cumbersome. Any

## suggestion is appreciated.




sessionInfo()
#> R version 4.3.2 (2023-10-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 12 (bookworm)
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
#>  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
#>  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Europe/Brussels
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] sf_1.0-14           rnaturalearth_0.3.4 lubridate_1.9.3    
#>  [4] forcats_1.0.0       stringr_1.5.0       dplyr_1.1.3        
#>  [7] purrr_1.0.2         readr_2.1.4         tidyr_1.3.0        
#> [10] tibble_3.2.1        ggplot2_3.4.4       tidyverse_2.0.0    
#> 
#> loaded via a namespace (and not attached):
#>  [1] s2_1.1.4                utf8_1.2.4              generics_0.1.3         
#>  [4] class_7.3-22            KernSmooth_2.23-22      stringi_1.7.12         
#>  [7] lattice_0.22-5          hms_1.1.3               digest_0.6.33          
#> [10] magrittr_2.0.3          evaluate_0.23           grid_4.3.2             
#> [13] timechange_0.2.0        fastmap_1.1.1           jsonlite_1.8.7         
#> [16] e1071_1.7-13            DBI_1.1.3               httr_1.4.7             
#> [19] fansi_1.0.5             scales_1.2.1            cli_3.6.1              
#> [22] rlang_1.1.2             units_0.8-4             munsell_0.5.0          
#> [25] reprex_2.0.2            withr_2.5.2             yaml_2.3.7             
#> [28] tools_4.3.2             tzdb_0.4.0              colorspace_2.1-0       
#> [31] vctrs_0.6.4             R6_2.5.1                proxy_0.4-27           
#> [34] lifecycle_1.0.4         classInt_0.4-10         fs_1.6.3               
#> [37] pkgconfig_2.0.3         pillar_1.9.0            gtable_0.3.4           
#> [40] glue_1.6.2              Rcpp_1.0.11             rnaturalearthdata_0.1.0
#> [43] xfun_0.41               tidyselect_1.2.0        knitr_1.45             
#> [46] farver_2.1.1            htmltools_0.5.7         labeling_0.4.3         
#> [49] rmarkdown_2.25          wk_0.9.0                compiler_4.3.2         
#> [52] sp_2.1-1

Created on 2023-11-24 with reprex v2.0.2


Solution

  • If you can switch from Natural Erath data to Eurostat GISCO shapes (accessible through {giscoR} package), you'd get just one shape for Cyprus.

    When joining sf object and a data.frame, consider using right_join() to get sf object as a result while keeping all records of df_plot.

    library(sf)
    #> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
    library(dplyr) 
    library(ggplot2)
    
    # using Eurostat GISCO shapes and ETRS89 / ETRS-LAEA projection
    countries_sf <- giscoR::gisco_get_countries(year = "2020", epsg = "3035")
    bb_3035 <- st_bbox(c(xmin = -20, ymin = 30, xmax = 45, ymax = 73), crs = "WGS84") |>
      st_as_sfc() |> 
      st_transform(st_crs(countries_sf))
    countries_sf <- st_crop(countries_sf, bb_3035)
    
    # we want to have sf object as a first argument while joining 
    # to df_plot so all df_plot records would remain, hence the right_join
    countries_sf |>  
      right_join(df_plot, by = join_by(ISO3_CODE == iso3)) |>
      ggplot() +
      geom_sf(data = countries_sf) +
      geom_sf(aes(fill = percentage_share), col = "black", linewidth = 0.3) +
      facet_wrap(~expenditure_year, nrow = 2) +
      coord_sf(expand = FALSE) +
      theme_minimal() +
      theme(legend.position = "bottom")
    

    Created on 2023-11-24 with reprex v2.0.2