6

I have to analyze some CORDEX models which have a rotated pole projection. The projection is described in this document, see page 10, Table 1, row "Europe high res".

I have tried to describe the proj4 string like so: +proj=ob_tran +o_proj=longlat +o_lon_p=-162.0 +o_lat_p=39.25 +lon_0=180 +lat_0=0, but it seems that this is not correct, at least using sf in R:

st_crs(rsf) <- "+proj=ob_tran +o_proj=longlat +o_lon_p=-162.0 +o_lat_p=39.25 +lon_0=180 +lat_0=0"
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.0 +o_lat_p=39.25 +lon_0=180 +lat_0=0', missing crs returned

What is the correct CRS for this projection?

1
  • There is very valuable answer from Patrick Brown. The correct CRS is +proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=39.25 +lon_0=18 +to_meter=0.01745329 Nov 6 at 0:31

3 Answers 3

5

With the help of Edzer and others we conjured a solution that works. See the GitHub issue here. It is a workaround and is not the most beautiful approach, but it works.

Basically using the proj.4 crs "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329" works, with the caveat of using R lwgeom to project instead of sf::st_transform, since GDAL currently does not understand this projection (it will in the next version).

# d is the input dataset saved in an sf object
x <- lwgeom::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 +to_meter=0.01745329"))
sf:st_crs(x) = NA # ?! This is mandatory for some reason
sf::st_crs(x) = 4326
sf::plot(x, axes = TRUE)

enter image description here

4
  • 1
    This proj4 string isn't accepted in QGIS. Any idea why? QGIS doesn't give much detail.. just "The proj4 definition of 'CORDEX EUR-44' is not valid." ... Maybe a better question: Which version of GDAL was "the next version", when you posted this?
    – naught101
    Jan 22, 2020 at 5:13
  • @naught101 as you can see in R I had some issues as well, I have no idea what's the issue with this proj4 string. I think the version of GDAL I talk about was going to be GDAL3.0, but I might be wrong
    – AF7
    Jan 23, 2020 at 7:35
  • OK, well I'm stuck on 2.4 for the time being. Maybe I will try and upgrade.. Any idea if there's a simple way to test the validity of a proj4 string?
    – naught101
    Jan 25, 2020 at 1:29
  • @naught101 I'm sure there are several ways; google returns for example this: rdrr.io/github/lbusett/MODIStsp/man/check_proj4string.html, but try man proj for a more direct way using the command line proj program (on Unix).
    – AF7
    Jan 26, 2020 at 9:23
1

Another simpler approach is to use the "terra" package to open the nc or raster files and then assign the projection.

library(terra)

# Upload data where d is the input dataset in a nc or raster object
     x <- terra::rast("d")

# Define projection    
    crs<-"+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329"

# Change projection    
    x <- project(x, crs)

enter image description here

5
  • 1
    I tried this (using terra 1.5-12) but got the following error : Error: [project] cannot create dataset from source In addition: Warning message: The NETCDF driver does not support update access to existing datasets. (GDAL error 6). Would it be the same to do: terra::crs(x) <- "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329"?
    – User1234
    Jan 19, 2022 at 17:18
  • Yep. It should be the same. ## S4 replacement method for signature 'SpatRaster' crs(x)<-value
    – LFG
    Jan 19, 2022 at 19:11
  • Thank you! When I plot the resulting raster, I got it with the wrong extent and axis labels (the same as plotting the raw raster). Thus, looks like the terra::crs(x) is not working or I am missing something?
    – User1234
    Jan 20, 2022 at 9:31
  • It seems to be related to the version. I am using terra version 1.2-5
    – LFG
    Jan 20, 2022 at 19:03
  • 1
    Yes. It is working with terra 1.2-5. Thank you!
    – User1234
    Jan 21, 2022 at 10:21
0

My file had the following metadata gdal_metadata(file1)

[29] "rotated_pole#grid_mapping_name=rotated_latitude_longitude"
[30] "rotated_pole#grid_north_pole_latitude=31.75831245449315"
[31] "rotated_pole#grid_north_pole_longitude=87.59703130293302"           
[32] "rotated_pole#north_pole_grid_longitude=0"   

everything worked with the following proj4 string

paste0("+proj=ob_tran +o_proj=longlat +o_lon_p=", 0,
                " +o_lat_p=", 31.75831245449315, 
                " +lon_0=", 87.59703130293302-180,
                " +to_meter=0.01745329")

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.