Search code examples
rmapsprojectionrnaturalearth

R + Naturalearth: ETRS89 / ETRS-LAEA projection


Please have a look at the reprex below. I would like to change the coordinate system and use, instead of the equal-earth projection, the ETRS89 / ETRS-LAEA projection.

See

https://spatialreference.org/ref/epsg/etrs89-etrs-laea/

I simply do not know if that is implemented somewhere in R and how it is called. I believe that if I am lucky, I just need to change one line in the reprex.

Any suggestion is welcome!

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

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






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



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

### how is the ETRS89 / ETRS-LAEA projection called?


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



##Now a test. I have already transformed the coordinated and
## identified the window limits in the new coordinates 



gpl <- ggplot() + geom_sf(data = worldmap_trans) +
    coord_sf(xlim = disp_win_coord[,'X'], ylim = disp_win_coord[,'Y'],
             datum = target_crs, expand = FALSE) +
    theme_bw()

gpl


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.9000 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] styler_1.10.2           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           R.oo_1.25.0            
#> [16] R.cache_0.16.0          jsonlite_1.8.7          R.utils_2.12.2         
#> [19] e1071_1.7-13            DBI_1.1.3               httr_1.4.7             
#> [22] fansi_1.0.5             scales_1.2.1            cli_3.6.1              
#> [25] rlang_1.1.2             units_0.8-4             R.methodsS3_1.8.2      
#> [28] munsell_0.5.0           reprex_2.0.2            withr_2.5.2            
#> [31] yaml_2.3.7              tools_4.3.2             tzdb_0.4.0             
#> [34] colorspace_2.1-0        vctrs_0.6.4             R6_2.5.1               
#> [37] proxy_0.4-27            classInt_0.4-10         lifecycle_1.0.4        
#> [40] fs_1.6.3                pkgconfig_2.0.3         pillar_1.9.0           
#> [43] gtable_0.3.4            Rcpp_1.0.11             glue_1.6.2             
#> [46] rnaturalearthdata_0.1.0 xfun_0.41               tidyselect_1.2.0       
#> [49] knitr_1.45              farver_2.1.1            htmltools_0.5.7        
#> [52] rmarkdown_2.25          compiler_4.3.2          sp_2.1-1

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


Solution

  • Your reprex is far more complex than it needs to be. You can specify the crs inside coord_sf. The limits can be worked out in advance using st_transform

    library(tidyverse)
    library(rnaturalearth)
    library(sf)
    
    lims <- st_sfc(st_point(c(-20, 30)), st_point(c(45, 73)), crs = 4326) |>
      st_transform(3035) |>
      st_bbox()
    
    ne_countries(scale = "medium",
                 type = 'countries',
                 returnclass = "sf") |>
      ggplot() + 
      geom_sf() +
      coord_sf(crs = 3035, xlim = lims[c(1, 3)], ylim = lims[c(2, 4)]) +
      theme_bw()
    

    enter image description here