Based on the input given in r-spatial/sf#651, several examples exist on how to bring EURO-CORDEX rotated pole projection data to e.g. WGS84, using vector (sf) or raster (terra) objects as input, as described in this answer.
I'm interested in the exact opposite case here. Given some vector administrative boundary dataset in WGS84, I'd like to be able to reproject this to EURO-CORDEX rotated pole projection to overlay the datasets for visualization and processing purposes without the need for cost-intensive resampling of raster data altering original grid values.
Not sure if the challenge here is because of the use of proj4strings for CRS definition, which I'm not very used to, or because the initial solution is somewhat "hacky", but I can't wrap my head around how to approach this.
Since you need to register to obtain the original NetCDF4 data, I created a small subset for demonstration purposes which is available here.
Here is the preparation code chunk:
library("terra") #> terra 1.7-83 fname <- "tasmax_EUR-11_ICHEC-EC-EARTH_rcp45_r12i1p1_SMHI-RCA4_v1_day_20510101-20551231.nc" r <- terra::rast(fname)[[1]] r <- terra::flip(r) terra::writeRaster(r, "tasmax_EUR-11_ICHEC-EC-EARTH_rcp45_r12i1p1_SMHI-RCA4_v1_day_20510101-20551231_subset.tiff") So this is basically the point of entry (and obviously the crs definition is wrong):
fname <- "tasmax_EUR-11_ICHEC-EC-EARTH_rcp45_r12i1p1_SMHI-RCA4_v1_day_20510101-20551231_subset.tiff" x <- terra::rast(fname) x #> class : SpatRaster #> dimensions : 412, 424, 1 (nrow, ncol, nlyr) #> resolution : 0.11, 0.11 (x, y) #> extent : -28.43, 18.21, -23.43, 21.89 (xmin, xmax, ymin, ymax) #> coord. ref. : lon/lat WGS 84 (EPSG:4326) #> source : tasmax_EUR-11_ICHEC-EC-EARTH_rcp45_r12i1p1_SMHI-RCA4_v1_day_20510101-20551231_subset.tiff #> name : tasmax_1 #> min value : 241.0844 #> max value : 303.5063 #> unit : K #> time (days) : 2051-01-01 terra::plot(x) 
In the end, I'd like to know which crs argument I need to use in order to bring e.g. NE data to this projection via sf::st_transform():
library("sf") #> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE ne_countries <- rnaturalearth::ne_countries(scale = 50) |> sf::st_geometry() ne_countries #> Geometry set for 242 features #> Geometry type: MULTIPOLYGON #> Dimension: XY #> Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961 #> Geodetic CRS: WGS 84 #> First 5 geometries: #> MULTIPOLYGON (((31.28789 -22.40205, 31.19727 -2... #> MULTIPOLYGON (((30.39609 -15.64307, 30.25068 -1... #> MULTIPOLYGON (((53.08564 16.64839, 52.58145 16.... #> MULTIPOLYGON (((104.064 10.39082, 104.083 10.34... #> MULTIPOLYGON (((-60.82119 9.138379, -60.94141 9... 