Sf: CRS for rotated pole projection works with rgdal, but not with sf

[adrfantini]

The EURO-CORDEX climate models use a rotated pole projection described in this document (page 10, table 1)

I've been trying to find the correct CRS to use with SF, but I cannot find it. I tried for example: +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=WGS84, which is accepted by rgdal::CRS and raster::projection, but sf::st_crs fails:

st_crs("+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=WGS84")
Coordinate Reference System: NA
Warning message:
In CPL_crs_from_proj4string(x) :
  Cannot import crs from PROJ.4 string `+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180', missing crs returned

rgdal::checkCRSArgs("+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180  +ellps=WGS84")
[[1]]
[1] TRUE

[[2]]
[1] "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=WGS84"

Why is rgdal ok with this, but not sf?

[edzer]

ob_tran rings a bell, I remember @rsbivand spending a lot of time on it for rgdal: the problems seemed to be in proj.4, rgdal working around it.. I have the feeling that I am not going to do the same thing for sf.

[adrfantini]

Does this resonate with that bell?
Should be fixed now.

[rsbivand]

Two points. Unless you use the GDAL NetCDF driver (say rather than ncdf4), no changes in the GDAL NetCDF driver will be picked up in sf or rgdal. I agree that sf::st_transform() is not a reasonable place to start; lwgeom::st_transform_proj() uses proj directly, so might be the best place to handle the direction reversal shown in the rgdal diffs. And the src diffs, note the swapped DEG2RAD too.

[edzer]

Here is an example of using ob_tran with lwgeom::st_transform_proj:

library(sf)
demo(nc)
library(lwgeom)
library(maptools)
data(wrld_simpl)
x = st_transform_proj(st_as_sf(wrld_simpl), "+proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 +lon_0=-90")
g = st_graticule()
gg = st_transform_proj(g, "+proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 +lon_0=-90")
plot(x[1], graticule = gg)

[adrfantini]

Thanks for the hints, but I need to be able to set the starting projection of the data, not to remap it to something else. What is wrong with my CRS so that I get this Cannot import crs from PROJ.4 string '+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180', missing crs returned error? Does this have to do with the latest GDAL version (2.2.3-1 on ArchLinux) not including (yet) the fix I linked to above, or does it have to do with sf?

[edzer]

That is currently done by the dev version of lwgeom:

st_crs(x) = NA
orig = st_transform_proj(x, c("+proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 +lon_0=-90", "+proj=longlat"))

I'll probably try to make that somewhat more intuitive.

[adrfantini]

Oh, ok, sorry. Yep, that's quite unintuitive.

After testing, I still have problems with this, and it might be because my starting CRS is wrong.

Here is some example data.

library(sf)
library(lwgeom)
inputProj <- '+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180'
d <- st_read("data.sqlite")
st_crs(d) <- inputProj #does not work
st_crs(d) <- NA
dProj <- st_transform_proj(d, c(inputProj, "+proj=longlat"))
plot(dProj, axes=TRUE)

Outputs this:

If I instead use the CRS you provided:

dProj2 <- st_transform_proj(d, c("+proj=ob_tran +o_proj=moll +o_lat_p=45 +o_lon_p=-90 +lon_0=-90", "+proj=longlat"))
plot(dProj2, axes=TRUE)

Outputs this:

Closer, but no sugar. Is my initial CRS wrong (as reported by st_crs), or is there anything else that i'm missing? The specs document includes the following table which should contain the relevant information (row "Europe high res").

[edzer]

Do you have an image of what your result should look like?

[adrfantini]

I think in epsg::4326 it should be something along these lines:

(it might not be exactly the same: here I used another file with an LCC projection that I know for sure is correct)

[rsbivand]

By the way, released 2.2.3 does not have the NetCDF obtran fix, which is in trunk for 2.3.0 line 295. I tried with 2.2.3 and it is not present, missed the RCs by days.

[edzer]

Looks like some juggling is needed with +to_meter; inspired by the gdal hack (rgdal does sth similar):

library(sf)
# Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3
library(lwgeom)
# Linking to liblwgeom 2.5.0dev r16016, GEOS 3.5.1, proj.4 4.9.3
crs = c("+proj=longlat", "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=40 +to_meter=0.01745329")
(pole = st_point(c(0,90)))
# POINT (0 90)
st_transform_proj(pole, crs)
# POINT (-162 40.00001)
(p = st_point(c(0,80)))
# POINT (0 80)
(pp = st_transform_proj(p, crs))
# POINT (-162 30)
st_transform_proj(pp, crs[2:1])
# POINT (2.537951e-13 80)

(requires new dev install of lwgeom)

[adrfantini]

@rsbivand thanks, I guess that's why the st_crs() <- error. This will automagically work when GDAL 2.3.0 will be out, right @edzer? Any news on when that's gonna happen? (I douldn't find a GDAL release schedule).

@edzer, with +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329 I get this with the data I linked to previously:

Which clearly is still shifted and warped. Using liblwgeom 2.5.0dev r16016, GEOS 3.6.2, proj.4 4.9.3
You get something different?

[edzer]

It looks like the inverse projection needs different parameters. I looked at this change set from @rouault, but couldn't understand the difference between grid_north_pole_longitude and north_pole_grid_longitude, and how these relate to the values in your table or example above. Do you?

In any case, the example (search for "hack") points to using a different lon_0, however trying 18 (180 + -162) didn't help much when leaving everything else the same.

[rsbivand]

No magic at all. The GDAL linked to sf and rgdal would have to be built with HDF5 and NetCDF support, which is maybe OK if you build GDAL from source, but adds to the burdens of the Windows and MacOS binary maintainers. It is only going to be painful, the modellers' choices of output grids show how little they grasp about the needs of other users than they have been working for.

Do we have a complete NetCDF raster in the incrimiated projection to try to warp? I have an 80Mb tasmax 1950 daily temperatures, which I used to check that 2.2.3 doesn't recognise the metadata (I extracted the metadata from gdalInfo() to build my proj4string).

[adrfantini]

@edzer I have also tried to play around with the values to no avail so far. I do not understand much of that commit, but I'll try again.
@rsbivand sure, here is an example file in that projection from EURO-CORDEX, high res grid (~0.11).

[edzer]

The creators mindset seems to be: hey, this file contains layers with lon and lat for each pixel, why worry about describing the transformation parametrically.

[adrfantini]

@edzer
Yes. My thoughts exactly. You can plot that easily enough using... GrADS...

[rsbivand]

GDAL dev:

$ apps/gdalinfo --version
GDAL 2.3.0dev, released 2017/99/99
$ apps/gdalinfo ~/tmp/bigshape/sftlf_EUR-11_ICHEC-EC-EARTH_rcp45_r3i1p1_DMI-HIRHAM5_v1_fx.nc
Driver: netCDF/Network Common Data Format
Files: /home/rsb/tmp/bigshape/sftlf_EUR-11_ICHEC-EC-EARTH_rcp45_r3i1p1_DMI-HIRHAM5_v1_fx.nc
Size is 424, 412
Coordinate System is `'
Origin = (-28.430000000811638,21.889999083358877)
Pixel Size = (0.110000001623275,-0.109999997772440)
Metadata:
...
  rlat#axis=Y
  rlat#long_name=latitude in rotated-pole grid
  rlat#standard_name=grid_latitude
  rlat#units=degrees
  rlon#axis=X
  rlon#long_name=longitude in rotated-pole grid
  rlon#standard_name=grid_longitude
  rlon#units=degrees
...
Band 1 Block=424x1 Type=Float32, ColorInterp=Undefined
  NoData Value=1.00000002004087734e+20
  Unit Type: %
  Metadata:
    coordinates=lon lat
    grid_mapping=rotated_pole
    long_name=Land Area Fraction

So still doesn't provide the metadata connection at SVN 41567. No magic here.

[edzer]

@adrfantini maybe we can discuss some of this stuff at EGU; I'll be there from Tue - Fri morning, regrettably not during your poster.

[rsbivand]

And the sample data set is less than helpful (this with GDAL 2.2.3):

> gI <- GDALinfo("sftlf_EUR-11_ICHEC-EC-EARTH_rcp45_r3i1p1_DMI-HIRHAM5_v1_fx.nc")
> attr(gI, "mdata")
 [1] "NC_GLOBAL#CDI=Climate Data Interface version 1.4.0.1"                               
 [2] "NC_GLOBAL#CDO=Climate Data Operators version 1.4.0.1 (http://www.mpimet.mpg.de/cdo)"
 [3] "NC_GLOBAL#[email protected]"                                                       
 [4] "NC_GLOBAL#Conventions=CF-1.6"                                                       
 [5] "NC_GLOBAL#CORDEX_domain=EUR-11"                                                     
 [6] "NC_GLOBAL#creation_date=2013-01-15 12:31:59"                                        
 [7] "NC_GLOBAL#driving_experiment=ICHEC-EC-EARTH,rcp45,r3i1p1"                           
 [8] "NC_GLOBAL#driving_experiment_name=rcp45"                                            
 [9] "NC_GLOBAL#driving_model_ensemble_member=r3i1p1"                                     
[10] "NC_GLOBAL#driving_model_id=ICHEC-EC-EARTH"                                          
[11] "NC_GLOBAL#experiment=ICHEC-EC-EARTH,rcp45,r3i1p1"                                   
[12] "NC_GLOBAL#experiment_id=rcp45"                                                      
[13] "NC_GLOBAL#frequency=fx"                                                             
[14] "NC_GLOBAL#institute_id=DMI"                                                         
[15] "NC_GLOBAL#institution=DMI"                                                          
[16] "NC_GLOBAL#model_id=DMI-HIRHAM5"                                                     
[17] "NC_GLOBAL#NCO=4.0.9"                                                                
[18] "NC_GLOBAL#product=output"                                                           
[19] "NC_GLOBAL#project_id=CORDEX"                                                        
[20] "NC_GLOBAL#rcm_version_id=v1"                                                        
[21] "NC_GLOBAL#tracking_id=9175edb8-5f0f-11e2-b420-6c626dd8513d"                         
[22] "rlat#axis=Y"                                                                        
[23] "rlat#long_name=latitude in rotated-pole grid"                                       
[24] "rlat#standard_name=grid_latitude"                                                   
[25] "rlat#units=degrees"                                                                 
[26] "rlon#axis=X"                                                                        
[27] "rlon#long_name=longitude in rotated-pole grid"                                      
[28] "rlon#standard_name=grid_longitude"                                                  
[29] "rlon#units=degrees"                                                                 
[30] "sftlf#coordinates=lon lat"                                                          
[31] "sftlf#grid_mapping=rotated_pole"                                                    
[32] "sftlf#long_name=Land Area Fraction"                                                 
[33] "sftlf#missing_value=1e+20"                                                          
[34] "sftlf#standard_name=land_area_fraction"                                             
[35] "sftlf#units=%"                                                                      
[36] "sftlf#_FillValue=1e+20"                                                             

With my tasmax file, they did add the rotated pole values:

>  gI <- GDALinfo("tasmax_EUR-11_ICHEC-EC-EARTH_historical_r1i1p1_KNMI-RACMO22E_v1_day_19500101-19501231.nc")
>  attr(gI, "mdata")[grep("^tasmax", attr(gI, "mdata"))]
[1] "tasmax#cell_methods=time: maximum"                          
[2] "tasmax#coordinates=lon lat height"                          
[3] "tasmax#grid_mapping=rotated_pole"                           
[4] "tasmax#long_name=Daily Maximum Near-Surface Air Temperature"
[5] "tasmax#missing_value=1e+20"                                 
[6] "tasmax#standard_name=air_temperature"                       
[7] "tasmax#units=K"                                             
[8] "tasmax#_FillValue=1e+20"                                    
>  attr(gI, "mdata")[grep("^rotated", attr(gI, "mdata"))]
[1] "rotated_pole#grid_mapping_name=rotated_latitude_longitude"
[2] "rotated_pole#grid_north_pole_latitude=39.25"              
[3] "rotated_pole#grid_north_pole_longitude=-162"              

giving me "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"

as the key. I could then "transform" longlat shorelines to the oddity, and they seemed to fit.

[adrfantini]

@edzer yes, I'll for sure be there wed-fri (if I ever manage to finish the work for that poster).
@rsbivand Looks great, I'll test it tomorrow (not in office and ssh server down right now), then report back.

Thanks.

[adrfantini]

@rsbivand , using your CRS I do not get different results compared to what I have shown before. See below:

#data downloaded from
#http://s000.tinyupload.com/download.php?file_id=58259885377323739503&t=5825988537732373950370534
d <- st_read("data.sqlite")
plot(st_transform_proj(d, c(st_crs(4326)$proj4string, "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180")), axes=TRUE)

As you can see the axes are wrong. If I do not reproject the data and instead reproject the shorelines on top of it, as you did... the axes are still going to be wrong, with (0,0) in central Europe.

[edzer]

x <- st_transform_proj(d, c(st_crs(4326)$proj4string, "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=175 +to_meter=0.01745329")
st_crs(x) = NA # ?!
st_crs(x) = 4326
plot(x, axes = TRUE)

comes pretty close:

[adrfantini]

Mmh, looks good. Was that trial and error, Edzer?

[edzer]

Don't tell anyone.

[edzer]

Don't use the lon_0=175, use 180. Plot a world map on top of it, e.g.

data(wrld_simpl, package = "maptools")
w = st_as_sf(wrld_simpl)
plot(st_geometry(w))
plot(x[1], axes=F, add=T)

[adrfantini]

x <- st_transform_proj(d, c(st_crs(4326)$proj4string, "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=175 +to_meter=0.01745329"))
st_crs(x) = NA # ?! (!)
st_crs(x) = 4326
bbox <- st_bbox(x)
borders <- st_as_sf(maps::map('world', plot = FALSE, fill = TRUE))

ggplot() + geom_sf(data=x, aes(fill=near.surface.air.temperature), color=NA) + geom_sf(data = borders, fill=NA, colour="black") + coord_sf(xlim = bbox[c(1,3)],ylim = bbox[c(2,4)])

Looks a bit off still :/

[adrfantini]

My bad, I used 175 instead of 180:

:+1:

[edzer]

Nice! Poster ready?

[adrfantini]

Ha! I'll start working on it in... about two weeks. Ahem. Could be worse, it could rain.

[edzer]

ready to close this one?

[adrfantini]

I'll pull the trigger. My understanding is that native sf support for this projection will happen only with the next version of GDAL, since this is what sf uses underneath.

Thanks for your help, I really appreciate the fact that you did not dismiss this as "not sf to be blamed here" and helped me out.

[rsbivand]

I doubt very much that the NetCDF driver support will work unless the NetCDF providers actually insert the proper metadata. Further, use lwgeom::st_transform_proj, not sf::st_transform to ensure direct contact with PROJ.4.

[adrfantini]

Looks like the files I am using actually are within the CF-conventions, it's just that they comply to chapter 5.6 (Horizontal Coordinate Reference Systems, Grid Mappings, and Projections) and not the (optional) 5.6.1 (Use of the CRS Well-known Text Format). I'm not sure GDAL is happy with this tho.

EDIT: link to the CF-C