Transform coordinate from rotated lat/lon to normal lat/lon WGS84

Question

I have access to an OPeNDAP service that allows me to peek into a humongous ocean wave dataset. It takes latitude and longitude as inputs (in fact even a range of the two) - but the coordinate system is a rotated one (presumably because the dataset focus is the North Sea and the Arctic):

grid_mapping_name: rotated_latitude_longitude
grid_north_pole_longitude: 140.0
grid_north_pole_latitude: 22.0
earth_radius: 6371000.0
proj4: +proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +R=6.371e+06 +no_defs

I want to query specific coordinates, but I do not know how to translate these coordinates (i.e. WGS84 lat/lon from Google Maps) into coordinates for this service. Anyone have an idea how to do this?

For example, let's say I want to transform the coordinate (57.93, 3.63) (somewhere out in the North Sea).

Any solution is welcome, using pyproj or something completely different.

EDIT: The OPeNDAP dataset can be found Here - it sadly does not offer a whole lot of information, but snippets of data can be extracted as ASCII.

EDIT2: This is the code I've tried, but it seems the proj4 string is not correctly interpreted:

from osgeo import osr
src = osr.SpatialReference()
tgt = osr.SpatialReference()
src.ImportFromEPSG(4326)
tgt.ImportFromProj4('+proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +R=6.371e+06 +no_defs')

transform = osr.CoordinateTransformation(src, tgt)
coords = transform.TransformPoint(-122, 46)
x,y = coords[0:2]

print(coords)

Returns:

(-13580977.876779377, 5780349.220256355, 0.0)
ERROR 6: No translation for an empty SRS to PROJ.4 format is known.
Traceback (most recent call last):
  File "coordtrans.py", line 24, in <module>
    coords = transform.TransformPoint(-122, 46)
  File "C:\tools\Anaconda3\envs\testEnv\lib\site-packages\osgeo\osr.py", line 1283, in TransformPoint
    return _osr.CoordinateTransformation_TransformPoint(self, *args)
NotImplementedError: Wrong number or type of arguments for overloaded function 'CoordinateTransformation_TransformPoint'.
  Possible C/C++ prototypes are:
    OSRCoordinateTransformationShadow::TransformPoint(double [3])
    OSRCoordinateTransformationShadow::TransformPoint(double [3],double,double,double)

Run in a conda environment with gdal from conda.

Answer 1

This was addressed in PROJ by Even Roualt in this issue: https://github.com/OSGeo/PROJ/issues/1525

With the fix in PROJ 6.1.1 (planned release July 1, 2019), you should be able to do this in pyproj if you install PROJ from source and install pyproj without the binary (https://pyproj4.github.io/pyproj/stable/installation.html#installing-from-source):

$ git clone https://github.com/pyproj4/pyproj.git
$ cd pyproj
pyproj$ PROJ_DIR=$HOME/scripts/pyproj/pyproj/proj_dir bash ci/travis/proj-dl-and-compile 6.1.1RC1
$ pip install -e .

Then you can run:

>>> from pyproj import show_versions
>>> show_versions()

System:
    python: 3.7.3 | packaged by conda-forge | (default, Mar 27 2019, 23:01:00)  [GCC 7.3.0]
executable: /home/snowal/miniconda3/envs/pyproj/bin/python
   machine: Linux-4.15.0-51-generic-x86_64-with-debian-buster-sid

PROJ:
      PROJ: 6.1.1
  data dir: /home/snowal/scripts/pyproj/pyproj/proj_dir/share/proj

Python deps:
    pyproj: 2.2.1
       pip: 18.0
setuptools: 41.0.1
    Cython: 0.29.10
     aenum: (2, 1, 2)

>>> from pyproj import Transformer
>>> trans = Transformer.from_crs('epsg:4326', '+proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +R=6.371e+06', always_xy=True)
>>> trans.transform(-122, 46)
(-44.37073742982609, 10.359937997567403)