2

I am working with a EUR-11 CORDEX dataset of daily temperature values (nc file) from ESGF database and EUROSTAT NUTS data (from here: https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts) as my shp file for the EU-countrie's borders.

My goal is to extract raster values from the nc file based on the polygons in the shp file.

    library(sf)
    library(ncdf4)
    library(raster)
    library(rgdal)
    library(lwgeom)

    # vector data    
    NUTS_data <- st_read("./data/NUTS_RG_60M_2016_4326/NUTS_RG_60M_2016_4326.shp")
        NUTS2 <- subset(NUTS_data,LEVL_CODE == 2)
        idx <- c("FRY1","FRY2","FRY3","FRY4","FRY5","ES70", "PT30","PT20") #getting rid of some small islands outside my study area
        for(i in 1:length(idx)){
          NUTS2 <- subset(NUTS2,NUTS_ID != idx[i])
        }
        EUstates = NUTS2 %>%
          group_by(CNTR_CODE) %>%
          summarize(
            urban  = sum(URBN_TYPE)) 
        EUstates <- EUstates[,c(1,3)]
    
    crs(EUstates)
        
    Coordinate Reference System:
    Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
    Warning:
    In wkt(x) : CRS object has no comment

    # raster data
    file <- "tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_SMHI-RCA4_v1a_3hr_199801010000-199812312100.nc" 
    ncfile <- nc_open(file)
    lon <- ncvar_get(nc, varid ="lon")
    lat <- ncvar_get(nc, varid="lat")
    TdmEUR98 <- ncvar_get(ncfile, varid="tas") 
    TymEUR98 <- rowMeans(TdmEUR98,na.rm = T,dims=2) # Calculate yearmean
    data <- as.data.frame(TymEUR98)
    r <- raster(t(data),xmn=min(lon),xmx=max(lon),ymn=min(lat),ymx=max(lat),
                              crs=CRS("+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329")) 
    r <- flip(r, direction="y")
    r.clip <- mask(r,EUstates)

The CRS for the raster is from the Q&A at Choosing CRS / PROJ4 string for EURO-Cordex rotated pole projection? Before I load my data in R, I preprocess the CORDEX data with CDO, calculating daily means from those 3h values and remapping the rotated pole projection to a regular lonlat grid.

When I plot my raster and vector, it looks good so far.

plot(r.clip)
plot(EUstates,col="transparent",add=T)

I am aware, that the CRS of my raster and vector differ up to this point, but they seem to match visually, so I just tried to extract...which did not work. ext contains only Values of NULL and NA

ext <- extract(r, EUstates, method='simple')
# Warning:
In .local(x, y, ...) :
  Transforming SpatialPolygons to the crs of the Raster

Trying to reproject both CRS as follwing did not get me anywhere either, because the raster and shp now don't align and the longitude and latitude seem off.

CRS("+init=epsg:4326") # as CRS 4326 seems to already be the projection of shp file
v <- st_transform_proj(EUstates,"+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
r.rep <- projectRaster(r.clip,crs = crs(v))
plot(r.rep)
plot(v,col="transparent",add=T)

enter image description here

I assume I have to change something about the projection or reproject my raster and/or my shp file, but I tried reprojecting both an I just can't get it right. So at this point I am just lost and don't really know what to do.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.