CODEWITHSUNDEEP
rterrargdalproj
jukiba
jukiba

Reputation: 1

Projection of EURO CORDEX data

I tried to project the EURO CORDEX data in R with the PROJ4 string. It is not working anymore (PROJ4 and RDGAL retired). Is there a more recent approach to reproject/transform the EURO CORDEX data in R based on a EPSG code or anything else? I am working with the terra package. Which libraries need to be installed in R and what is the correct command?

The EURO CORDEX data I use is downloaded from the Copernicus Data Storage CDS https://cds.climate.copernicus.eu/cdsapp#!/dataset/projections-cordex-domains-single-levels?tab=form (Domain: Europe, Experiment: Historical, Horizontal resolution: 0.11°, Temporal resolution: seasonal mean, Variable: 2m air temp, global climate model: MPI-M-MPI-ESM-LR, Regional climate model: MPI-CSC-REMO2009, ensemble member: r1i1p1, Start year: 2000, end year: 2005).

This is my code:

library(terra)
library(sf)

test <- rast(".../tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_sem_200012-200511.nc", subds="tas")

analysis_loc <- st_read("... .shp")

dput(analysis_loc)

#output: 
structure(list(id = c(1, 2, 3), geometry = structure(list(structure(c(13.4094290074978, 
        52.5208244273539), class = c("XY", "POINT", "sfg")), structure(c(13.3500999325312, 
        52.5145123927038), class = c("XY", "POINT", "sfg")), structure(c(13.2781309381902, 
        52.5050087446307), class = c("XY", "POINT", "sfg"))), n_empty = 0L, crs = structure(list(
            input = "WGS 84", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"), class = "crs"), class = c("sfc_POINT", 
        "sfc"), precision = 0, bbox = structure(c(xmin = 13.2781309381902, 
        ymin = 52.5050087446307, xmax = 13.4094290074978, ymax = 52.5208244273539
        ), class = "bbox"))), row.names = c(NA, -3L), class = c("sf", 
        "data.frame"), sf_column = "geometry", agr = structure(c(id = NA_integer_), class = "factor", levels = c("constant", 
        "aggregate", "identity")))

cor_crs <- "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329"
test2 <- project(test, cor_crs)

plot(test2[[1]])
plot(analysis_loc$geometry, add=T)

ext <- terra::extract(test2, analysis_loc)

in the extraction data frame all values are na, but on the map the points and the raster are aligned. why? rasterplot

Upvotes: 0

Views: 143

Answers (1)

jukiba
jukiba

Reputation: 1

"Extract" still does not work, but I managed to solve the problem with a workaround with "zonal":

poly <- buffer(analysis_loc, 100)
poly_zones <- zonal(test2 , poly, fun="mean", na.rm=TRUE)

Upvotes: 0

Related Questions