2

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.

  • 1
    Possible duplicate of gis.stackexchange.com/questions/226200/… – nmtoken Feb 13 '19 at 17:18
  • That does not really solve it, as I do not have the ESPG for the projection but rather a proj4 string. Or am I confused? – dingobar Feb 13 '19 at 19:58
  • Are you sure that the proj4 string is valid? Gdalsrsinfo gives errors about that. Do you have a link for test data or to documentation of that dataset? – user30184 Feb 13 '19 at 20:20
  • You don't need the EPSG code, just the parameters. – nmtoken Feb 13 '19 at 20:55
  • see gdal.org/gdaltransform.html ... The coordinate systems that can be passed are anything supported by the OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes (i.e. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prj file containing well known text. – nmtoken Feb 13 '19 at 21:08
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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