area definition for a rotated pole coordinate system

[larsbuntemeyer]

Hey all, thanks a lot! I am exploring pyresample a bit for the use in resampling land cover data as input for our regional climate model. I found the BucketResampler to be exactly what i need, especially the ability to compute fractions on a corser grids! That's a great tool! However, the regional modely mainly work on rotated grid transformations and i am struggling a little bit to get the transformation right.

Code Sample, a minimal, complete, and verifiable piece of code

So this would be my area definition, it is for Europe on a 0.44 degree resolution and rotated pole and an area extent defined in the projection coordinates:

from pyresample.geometry import AreaDefinition

area_id = 'EUR-44'
description = 'Europe'
proj_id = 'rotated_pole'
projection = '+proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=39.25 +lon_0=18 +datum=WGS84 +no_defs +type=crs'
width = 106
height = 103
area_extent = (-28.21, -23.21, 17.99, 21.67)

area_def = AreaDefinition(area_id, description, proj_id, projection,
                          width, height, area_extent = area_extent)

Problem description

I think, i don't get the projection and area extent to work together. E.g. when i plot the latitude and longitude coordinates from this area definition, it doesn't look like it's getting the area extent right (e.g. ranges from -180, 180):

import matplotlib.pyplot as plt

plt.imshow(area_def.get_lonlats()[0]) # plot longitudes
plt.colorbar()

Expected Output

I usually would do the transformation myself, with something like this (and i guess pyresample does something simliar under the hood?)

from pyproj import CRS, Transformer
import numpy as np

# resolution
dlon = (area_extent[2] - area_extent[0]) / (width-1) # 0.44
dlat = (area_extent[3] - area_extent[1]) / (height-1) # 0.44
# longitude, latitude in regular coordinates (rotated)
rlon = np.arange(0, width) * dlon + area_extent[0]
rlat = np.arange(0, height) * dlat + area_extent[1]

# transform into lon lat coordinates
transformer = Transformer.from_crs(CRS(projection), CRS("EPSG:4326"),  always_xy=True)
lon, lat = transformer.transform(np.meshgrid(rlon, rlat)[0], np.meshgrid(rlon, rlat)[1])

and to compare, the longitude plot looks like expected (e.g, i covers the expected range of longitudes for my European domain)

plt.imshow(lon)
plt.colorbar()

I want to mention that all those thinks work fine when i define my domain in ESPG:4326, e.g., like this:

area_id = 'EUR-44i'
description = 'Europe'
proj_id = 'regular'
projection = 'EPSG:4326'
width = 221
height = 103
area_extent = (-44.75, 21.75, 65.25, 72.75)

area_def = AreaDefinition(area_id, description, proj_id, projection,
                          width, height, area_extent = area_extent)

I am not sure, if what i expect is actually reasonable? Thanks a lot your help and apologies for the long issue!

Versions of Python, package at hand and relevant dependencies

pyproj=3.3.1
pyresample=1.27.0

[djhoese]

I have little to no experience with ob_tran so I'll have to look at that, but based on the repeating pattern I wonder if your degrees are being treated as radians somewhere. Or I wonder if the always_xy=True isn't being used somewhere important in pyresample...but then I would have expected EPSG:4326 so show similar issues since that is also a (lat, lon) axis order.

I'll try to play with this later today (currently in a few meetings).

[larsbuntemeyer]

Thanks! No rush with this, i'll also try to figure out what's happening under the hood in that BucketResampler.

[djhoese]

It looks like the old Proj interface of pyproj does not support the ob_tran projection. Swapping usage of Proj with Transformer in pyresample/geometry.py makes things work.

[djhoese]

Hm I'm not convinced this is the issue. I do get different results even if I (I think) keep inverse/forward transformations in the correct order.

In [39]: from pyproj import CRS, Transformer, Proj

In [40]: from pyresample.geometry import AreaDefinition
    ...:
    ...: area_id = 'EUR-44'
    ...: description = 'Europe'
    ...: proj_id = 'rotated_pole'
    ...: projection = '+proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=39.25 +lon_0=18 +datum=WGS84 +no_defs +type=crs'
    ...: width = 106
    ...: height = 103
    ...: area_extent = (-28.21, -23.21, 17.99, 21.67)
    ...:
    ...: area_def = AreaDefinition(area_id, description, proj_id, projection,
    ...:                           width, height, area_extent = area_extent)

In [41]: dlon = (area_extent[2] - area_extent[0]) / (width-1) # 0.44

In [42]: dlat = (area_extent[3] - area_extent[1]) / (height-1) # 0.44

In [43]: rlon = np.arange(0, width) * dlon + area_extent[0]

In [44]: rlat = np.arange(0, height) * dlat + area_extent[1]

In [45]: transformer = Transformer.from_crs(CRS(projection), CRS("EPSG:4326"),  always_xy=True)

In [46]: proj = Proj(CRS(projection), always_xy=True)

In [47]: lon, lat = transformer.transform(np.meshgrid(rlon, rlat)[0], np.meshgrid(rlon, rlat)[1], direction="INVERSE")

In [48]: plon, plat = proj(np.meshgrid(rlon, rlat)[0], np.meshgrid(rlon, rlat)[1], inverse=True)

In [49]: lon[-1]
Out[49]:
array([-44.07624682, -43.7300086 , -43.38228233, -43.0330607 , -42.68233664, -42.33010331, -41.97635413, -41.62108277, -41.26428314, -40.90594944, -40.54607614, -40.184658  , -39.82169006, -39.45716765, -39.09108643, -38.72344235, -38.3542317 , -37.98345108, -37.61109745, -37.23716809, -36.86166065, -36.48457315, -36.10590394, -35.7256518 , -35.34381585, -34.96039563, -34.57539106, -34.1888025 , -33.80063069, -33.41087683, -33.01954252, -32.62662983, -32.23214125,
       -31.83607974, -31.43844873, -31.0392521 , -30.63849423, -30.23617996, -29.83231464, -29.4269041 , -29.01995469, -28.61147327, -28.2014672 , -27.78994437, -27.3769132 , -26.96238265, -26.54636219, -26.12886187, -25.70989225, -25.28946447, -24.8675902 , -24.44428168, -24.01955172, -23.59341368, -23.16588148, -22.73696962, -22.30669317, -21.87506776, -21.44210959, -21.00783544, -20.57226266, -20.13540915, -19.69729341, -19.25793446, -18.81735194, -18.37556601,
       -17.9325974 , -17.48846739, -17.04319782, -16.59681106, -16.14933003, -15.70077817, -15.25117946, -14.80055838, -14.34893995, -13.89634966, -13.44281352, -12.988358  , -12.53301008, -12.07679715, -11.61974711, -11.16188827, -10.70324937, -10.24385957,  -9.78374845,  -9.32294594,  -8.86148239,  -8.39938849,  -7.93669527,  -7.47343411,  -7.00963668,  -6.54533497,  -6.08056125,  -5.61534804,  -5.14972813,  -4.68373451,  -4.21740041,  -3.75075923,  -3.28384458,
        -2.81669018,  -2.34932992,  -1.8817978 ,  -1.41412792,  -0.94635445,  -0.47851163,  -0.01063373])

In [50]: plon[-1]
Out[50]:
array([  27.75919733,   76.4829319 ,  100.03620221,  116.12660544,  131.27075639,  148.72748004,  170.9918605 , -161.82359056, -134.67941109, -112.48461023,  -95.07367147,  -79.94006397,  -63.80957136,  -40.08294096,    9.06245631,   66.44976136,   94.6644421 ,  111.89294281,  126.91977424,  143.44773056,  164.18291791, -169.72522207, -141.89668462, -118.1341201 ,  -99.59838087,  -84.1147086 ,  -68.68580201,  -48.23418879,   -8.4031626 ,   53.76652596,   88.53498023,
        107.50949673,  122.68977075,  138.50355059,  157.82755777, -177.46953417, -149.47338964, -124.21193483, -104.36426171,  -88.30232518,  -73.21538387,  -54.9954354 ,  -22.98708503,   37.99127022,   81.30015755,  102.86449255,  118.5155618 ,  133.83729372,  151.91811117,  175.08059331, -157.29604563, -130.74208726, -109.4305581 ,  -92.57124416,  -77.52957517,  -60.79969742,  -34.57396357,   19.77038853,   72.46649533,   97.81153335,  114.32616974,  129.38905867,
        146.42236493,  168.02195323, -165.21585042, -137.71999049, -114.85339834,  -96.98524852,  -81.73049986,  -65.94979581,  -43.77739157,    1.29875539,   61.3826768 ,   92.15030905,  110.04042156,  125.09761332,  141.29310658,  161.40613159, -173.0689464 , -145.1022516 , -120.68194787, -101.60585183,  -85.90190694,  -70.65908134,  -51.2704164 ,  -15.03506915,   47.38942252,   85.59772403,  105.56034138,  120.89973522,  136.47514896,  155.24429858,  179.2981178 ,
       -152.80055915, -126.95190627, -106.49322368,  -90.116301  ,  -75.0823753 ,  -57.57751476,  -28.28561629,   30.39167549,   77.74852178,  100.76120254,  116.72818743,  131.90964963,  149.5169873 ])

Am I missing something obvious here? Unless always_xy doesn't behave the way I think it should?

[djhoese]

I wonder if the datum note at the top of this page is causing us issues with using Proj:

https://pyproj4.github.io/pyproj/stable/api/proj.html#proj

pyproj.Proj is functionally equivalent to the proj command line tool in PROJ.

The PROJ docs say:

The `proj` program is limited to converting between geographic and
projection coordinates within one datum.

[larsbuntemeyer]

Thanks for this swift support! I had now a quick glance, and i think you are right with the units. Seems like Proj requires the coordinates for ob_tran in radians. E.g., if i change to:

area_extent = (np.deg2rad(-28.21), np.deg2rad(-23.21), np.deg2rad(17.99), np.deg2rad(21.67))

this works with the resampler and also i get same results from proj and transformer. I'll add a simple example.

[djhoese]

I guess I would have just thought that Proj still handled radians versus not, but I guess if it doesn't properly support multi-datum CRS definitions maybe this translation is lost somewhere...I don't know.

I was hoping to avoid bugging @snowman2 (the pyproj maintainer), but they would know best. Alan, for a PROJ.4 definition like:

projection = '+proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=39.25 +lon_0=18 +datum=WGS84 +no_defs +type=crs'

Would you expect different handling between Proj and Transformer objects? Specifically, it looks like Proj assumes projection units are radians, while Transformer does not.

[larsbuntemeyer]

This works and give same results:

from pyproj import CRS, Transformer, Proj
import numpy as np

area_id = 'EUR-44'
description = 'Europe'
proj_id = 'rotated_pole'
projection = '+proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=39.25 +lon_0=18 +datum=WGS84 +no_defs +type=crs'
width = 106
height = 103
area_extent_rad = (np.deg2rad(-28.21), np.deg2rad(-23.21), np.deg2rad(17.99), np.deg2rad(21.67))
area_extent_deg = (-28.21, -23.21, 17.99, 21.67)

# resolution
dlon_deg = (area_extent_deg[2] - area_extent_deg[0]) / (width-1) # 0.44
dlat_deg = (area_extent_deg[3] - area_extent_deg[1]) / (height-1) # 0.44
# longitude, latitude in regular coordinates (rotated)
rlon_deg = np.arange(0, width) * dlon_deg + area_extent_deg[0]
rlat_deg = np.arange(0, height) * dlat_deg + area_extent_deg[1]

# resolution
dlon_rad = (area_extent_rad[2] - area_extent_rad[0]) / (width-1) # 0.44
dlat_rad = (area_extent_rad[3] - area_extent_rad[1]) / (height-1) # 0.44
# longitude, latitude in regular coordinates (rotated)
rlon_rad = np.arange(0, width) * dlon_rad + area_extent_rad[0]
rlat_rad = np.arange(0, height) * dlat_rad + area_extent_rad[1]

# transform into lon lat coordinates
transformer = Transformer.from_crs(CRS("EPSG:4326"), CRS(projection), always_xy=True)
lon, lat = transformer.transform(np.meshgrid(rlon_deg, rlat_deg)[0], np.meshgrid(rlon_deg, rlat_deg)[1], direction="INVERSE")

# proj
proj = Proj(CRS(projection), always_xy=True)
plon, plat = proj(np.meshgrid(rlon_rad, rlat_rad)[0], np.meshgrid(rlon_rad, rlat_rad)[1], inverse=True)

assert np.allclose(lon, plon)
assert np.allclose(lat, plat)

[snowman2]

One thing that stands out to me is that an old version of pyproj is being used. Does this issue occur with the latest version of pyproj (3.6.0)?

[djhoese]

I'm on pyproj 3.5.0 and my conda-forge environment won't let me get pyproj 3.6.* because I also have GDAL installed and that is limiting the version of PROJ that can be installed.

[djhoese]

Created a new environment and I still see this difference.

Edit: with pyproj 3.6.0

[snowman2]

This may be a helpful way to view how PROJ views the transformations in this scenario:

>>> proj = Proj(projection)
>>> proj
<Other Coordinate Operation Transformer: ob_tran>
Description: PROJ-based coordinate operation
Area of Use:
- undefined
>>> proj.definition
'proj=ob_tran o_proj=longlat o_lon_p=0 o_lat_p=39.25 lon_0=18 datum=WGS84 no_defs ellps=WGS84 towgs84=0,0,0'
>>> trans = Transformer.from_crs("EPSG:4326", projection, always_xy=True)
>>> trans
<Conversion Transformer: pipeline>
Description: unknown
Area of Use:
- undefined
>>> trans.definition
'proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step proj=ob_tran o_proj=longlat o_lon_p=0 o_lat_p=39.25 lon_0=18 ellps=WGS84 step proj=unitconvert xy_in=rad xy_out=deg'

[snowman2]

The main difference I am seeing is that Proj adds towgs84=0,0,0. This is consistent with the other method for initializing a Transformer from a pipeline:

>>> pipe_trans = Transformer.from_pipeline(projection)
>>> pipe_trans
<Other Coordinate Operation Transformer: ob_tran>
Description: PROJ-based coordinate operation
Area of Use:
- undefined
>>> pipe_trans.definition
'proj=ob_tran o_proj=longlat o_lon_p=0 o_lat_p=39.25 lon_0=18 datum=WGS84 ellps=WGS84 towgs84=0,0,0'

[djhoese]

I've never used +towgs84, but I think I understand the concept. Is it correct to say that the added towgs84=0,0,0 is added because it isn't there in the user-provided definition? Also, is it an error (logically) for it to be added to this projection definition? By that I mean, by telling PROJ that towgs84 is 0,0,0 are we negating/conflicting with the nature of ob_tran with longlat (as you can tell I don't know how towgs84 works).

Is this a bug in pyproj or PROJ or is this a limitation of Proj when working with these types of projections? Or something else that I'm not understanding?

[snowman2]

It seems that proj_angular_input is not detecting the need for radians input in the inverse direction.

I added this method to the cython Transformer:

    def needs_radians(self, direction):
        cdef PJ_DIRECTION pj_direction = get_pj_direction(direction)
        return proj_angular_input(self.projobj, pj_direction)

>>> from pyproj import Proj
>>> projection = '+proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=39.25 +lon_0=18 +datum=WGS84 +no_defs +type=crs'
>>> proj = Proj(projection)
>>> proj._transformer.needs_radians("FORWARD")
1
>>> proj._transformer.needs_radians("INVERSE")
0

[snowman2]

Not sure if this has to do with it: https://github.com/OSGeo/PROJ/blob/2504995b860c3d6beacd13e1e05c519fbae30a3e/src/projections/ob_tran.cpp#L286-L290

[snowman2]

Related:

• https://github.com/OSGeo/PROJ/issues/1525
• https://github.com/OSGeo/PROJ/pull/1526

[djhoese]

@snowman2 Do you have any good suggestions for what the source/input CRS should be in a Transformer to match behavior of Proj? I've run into a failure in pyresample when switching to Transformer where a test that used EPSG:2056 are slightly different than they were when I was purely using Proj. My Transformer is using EPSG:4326 as the default lon/lat projection. What projection should be used to most closely reproduce what Proj is/was doing?

[snowman2]

From the example in the FAQ page:

from pyproj import CRS, Proj

crs = CRS("EPSG:2056")
proj = Proj(crs)
trans = Transformer.from_crs(crs.geodetic_crs, crs, always_xy=True)

[djhoese]

Thank you! I literally had the FAQ page open but got distracted by something else in the docs and didn't keep scrolling.

[djhoese]

@snowman2 I'm finally coming back to this and I've hit a bit of a snag. I implemented the crs.geodetic_crs usage so I can use Transformer while preserving the Proj behavior. However, this doesn't work for CRS's that use +pm=180 in their definition. Any ideas for an easy way to get the geodetic_crs without the prime meridian shift?

Example

In [1]: from pyproj import CRS, Transformer, Proj

In [2]: crs = CRS.from_proj4("+proj=merc +datum=WGS84 +ellps=WGS84 +pm=180 +no_defs")

In [3]: p = Proj(crs)

In [4]: trans = Transformer.from_crs(crs.geodetic_crs, crs)

In [5]: p(0, 0)
Out[5]: (-20037508.342789244, 0.0)

In [6]: trans.transform(0, 0)
Out[6]: (0.0, 0.0)

In [7]: p(-90, 0)
Out[7]: (10018754.171394622, 0.0)

In [8]: trans.transform(-90, 0)
Out[8]: (-10018754.171394622, 0.0)

[djhoese]

Or should I just hardcode 4326? Or fallback to 4326 if prime meridian isn't 0? Are current interfaces (historically) assume that a user requesting lons/lats wants the prime meridian 0 (4326-style) case. This legacy behavior basically ignores the idea that datums exist.

[snowman2]

Any ideas for an easy way to get the geodetic_crs without the prime meridian shift?

No easy ways come to mind.

[djhoese]

FYI #532 should fix this.