3

I am trying to create a dataset that shows future climatic condition in Europe aggregates to NUT2 administration level in R. To do so, i work with cordex data (domain EU 11), that is coming withrotated coordinates.

I used the data hurs_EUR-11_ICHEC-EC-EARTH_rcp85_r1i1p1_SMHI-RCA4_v1_day_20310101-20351231.nc (data here ) which is the relative humidity Its looks as follows :

 6 variables (excluding dimension variables):
    double height[]   
        axis: Z
        long_name: height
        positive: up
        standard_name: height
        units: m
    double time_bnds[bnds,time]   
    char rotated_pole[]   
        grid_mapping_name: rotated_latitude_longitude
        grid_north_pole_latitude: 39.25
        grid_north_pole_longitude: -162
    float hurs[rlon,rlat,time]   
        grid_mapping: rotated_pole
        _FillValue: 1.00000002004088e+20
        missing_value: 1.00000002004088e+20
        standard_name: relative_humidity
        long_name: Near-Surface Relative Humidity
        units: %
        coordinates: lon lat height
        cell_methods: time: mean
    double lon[rlon,rlat]   
        standard_name: longitude
        long_name: longitude
        units: degrees_east
    double lat[rlon,rlat]   
        standard_name: latitude
        long_name: latitude
        units: degrees_north
 4 dimensions:
    time  Size:1826   *** is unlimited ***
        standard_name: time
        units: days since 1949-12-01 00:00:00
        calendar: standard
        long_name: time
        bounds: time_bnds
        axis: T
    bnds  Size:2

2 "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named bnds BUT this dimension HAS NO DIMVAR! Code will probably fail at this point" rlon Size:424 standard_name: grid_longitude long_name: longitude in rotated pole grid units: degrees axis: X rlat Size:412 standard_name: grid_latitude long_name: latitude in rotated pole grid units: degrees axis: Y

First, i created an array

hurs.array <- ncvar_get(data2, "hurs") 

Then i create a grid with the projection as suggested in this post Choosing CRS / PROJ4 string for EURO-Cordex rotated pole projection? Note that the issue raised about this post in about an simple feature and not about raster.

r_brick <- brick(hurs.array, xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon), crs=CRS("+proj=ob_tran +o_proj=longlat +o_lon_p=-162.0 +o_lat_p=39.25 +lon_0=180 +lat_0=0"))

Now my question is how to transform this data into the EU projection? with

r_brickP<- projectRaster(r_brick, crs="+init=epsg:3035")

I get

not finiteno non-missing arguments to min;

What is going on?

Is this the right approach?

Is there a better approach to unrotate the data?

If yes which one?

What is my most efficient way to then perform a zonal statics with as layer with epsg:3035 projection?

I also tried the function rotate, but i got a line suggesting that this is not this method was not appropriate.

6
  • Note to anyone contemplating this - the data is 399Mb – Spacedman Feb 11 '20 at 12:48
  • What are your lon and lat whenyou make the brick? To correctly specify the extent of your brick they need to be in transformed coordinates, not true lat-long. – Spacedman Feb 11 '20 at 14:49
  • @Spacedman, i guess i have not transformed them, i just used the max and min from the rlong and rlat... – Catherine Pfeifer Feb 11 '20 at 17:41
  • so i dont seem to be able to get at the bottom of this. How do i transform these coordinate correctly, and in one go for the whole stack? – Catherine Pfeifer Feb 12 '20 at 8:14
  • The lon and lat variables in the netCDF appear to be the correct EPSG 4326 lat-long coordinates (features in the humidity match map features from OpenStreetMap). If the rlon and rlat variables are the X and Y coords of the grid system, we need to find the transform that turns those into lon and lat. If that fails you can work with a set of {lon, lat, h} points rather than a raster (but sloooow). – Spacedman Feb 12 '20 at 12:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.