6

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.

Improve this question
6
  • Note to anyone contemplating this - the data is 399Mb 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. 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). Commented Feb 12, 2020 at 12:43

1 Answer 1

Reset to default
0

It's been some time since you asked this but since I was playing around with CORDEX data anyway (and GDAL seems to be able to handle the projection in the meantime), I thought I'd give it a try.

Definitely not the most elegant approach, but so is the original approach handling vector data.

Since your sample data is not available anymore, I had to use my own, but this should be pretty similar from a structural point of view:

library("terra") #> terra 1.7.83 library("sf") #> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE library("rnaturalearth") fname <- "tasmax_EUR-11_ICHEC-EC-EARTH_rcp45_r12i1p1_SMHI-RCA4_v1_day_20060101-20101231.nc" # let's continue with the first layer only for demo purposes r <- terra::rast(fname)[[1]] # per default, Europe seems to be upside down r <- terra::flip(r) # crs definition provided in r-spatial/sf#651 crs <- "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329" # reproject to WGS84, override output crs definition x <- terra::project(r, crs) terra::crs(x) <- "epsg:4326" x #> class : SpatRaster #> dimensions : 339, 734, 1 (nrow, ncol, nlyr) #> resolution : 0.1496611, 0.1496611 (x, y) #> extent : -44.74486, 65.10636, 21.90072, 72.63582 (xmin, xmax, ymin, ymax) #> coord. ref. : lon/lat WGS 84 (EPSG:4326) #> source(s) : memory #> name : tasmax_1 #> min value : 224.3909 #> max value : 299.2951 #> unit : K #> time (days) : 2006-01-01

Not a big fan of reprojection casades, but this was my "brilliant" idea. Once you arrive in a well-known CRS, just continue to whatever system you need in the end:

x2 <- terra::project(x, "epsg:3035") x2 #> class : SpatRaster #> dimensions : 619, 1056, 1 (nrow, ncol, nlyr) #> resolution : 10189.62, 10189.62 (x, y) #> extent : -1042792, 9717446, -95413.04, 6211962 (xmin, xmax, ymin, ymax) #> coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) #> source(s) : memory #> name : tasmax_1 #> min value : 224.7869 #> max value : 299.2935 #> unit : K #> time (days) : 2006-01-01 # inspect visually terra::plot(x2) ne <- rnaturalearth::ne_countries(scale = 50) |> sf::st_geometry() |> sf::st_transform("epsg:3035") plot(ne, add = TRUE)

EURO-CORDEX data overlayed with NE countries in EPSG:3035

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.