5

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
    Commented Feb 11, 2020 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
    Commented Feb 11, 2020 at 14:49
  • @Spacedman, i guess i have not transformed them, i just used the max and min from the rlong and rlat... Commented Feb 11, 2020 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? Commented Feb 12, 2020 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
    Commented Feb 12, 2020 at 12:43

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.