General Oblique Transformation yields unexpected results using o_lat_c, o_lon_c and o_alpha

[Jeitan]

I have tried to use the General Oblique Transformation to obtain a rotated Oblique Mercator projection, but the results of +ob_tran are unexpected. If I construct it by using a new pole (as per the documentation here), it works as expected. But if I use the "rotate the projection around a given point" option, the results are unexpected and don't match the version that uses a pole. I don't know if it is an error in my understanding or an error in PROJ. I have a question here on GIS stack exchange where someone else confirmed my results.

If the following is not a bug but is an error in my understanding, then perhaps an explanation something like #3908 is requesting would help?

Example of problem

I am using pyproj to interface with PROJ, but pyproj is not the culprit as the answer at the stack exchange question above confirmed. I'm using the following code to create a projection then test its origin, and what the geodetic coordinates of 1 approximate degree in the +X and +Y directions yield.

import numpy as np
from pyproj import CRS, Transformer

geod = CRS.from_proj4("+proj=latlon")

custom3 = CRS.from_proj4(
    "+proj=ob_tran +o_proj=omerc +lonc=0 +lat_0=0 +alpha=0 +o_lat_p=90 +o_lon_p=0"
)

cust32geo = Transformer.from_crs(custom3, geod)
print("Custom General Oblique to mimic Oblique Mercator, using pole")
print("Origin: ", np.round(cust32geo.transform(0, 0)))
print("+X: ", np.round(cust32geo.transform(111e3, 0)))
print("+Y: ", np.round(cust32geo.transform(0, 111e3)))

custom4 = CRS.from_proj4(
    "+proj=ob_tran +o_proj=omerc +lonc=0 +lat_0=0 +alpha=0 +o_lat_c=0 +o_lon_c=0 +o_alpha=0"
)

cust42geo = Transformer.from_crs(custom4, geod)
print("Custom General Oblique to mimic Oblique Mercator, using center and rotation")
print("Origin: ", np.round(cust42geo.transform(0, 0)))
print("+X: ", np.round(cust42geo.transform(111e3, 0)))
print("+Y: ", np.round(cust42geo.transform(0, 111e3)))

Problem description

The above code prints the following:

Custom General Oblique to mimic Oblique Mercator, using pole
Origin:  [0. 0.]
+X:  [1. 0.]
+Y:  [0. 1.]
Custom General Oblique to mimic Oblique Mercator, using center and rotation
Origin:  [90.  0.]
+X:  [90. -1.]
+Y:  [91.  0.]

I expected the second version to be identical to the first ... but instead, the origin is at 0 degrees latitude and 90 degrees longitude, instead of (0, 0) like I expected. Also, the X and Y coordinate axes are rotated so +X is negative latitude and +Y is positive longitude.

ETA: I also get the same exact behavior if I use +o_proj=tmerc instead.

ETA 2: If I set +o_alpha=90, I get the X and Y axes oriented as I expect, which does make sense. The origin still shifts though.

Expected Output

I expected identical outputs from the two projections above.

Environment Information

Using pyproj 3.6.1 which according to the documentation contains PROJ 9.3.0.
Environment: Windows with an installation of Miniconda

Installation method

Installed pyproj using Poetry (pip under the hood)

[jjimenezshaw]

Out of curiosity: are you trying to rotate a transverse Mercator using the parameters of the oblique Mercator, or is your initial projection already oblique?

[Jeitan]

@jjimenezshaw My true use-case is actually for plotting - I have values that are in an oblique projection with a particular lat/lon as the origin and a particular rotation from north (it's not really a true oblique mercator, but close enough for the distances involved) , and I want to display them together with some other values that are in lat/lon. I'm using Cartopy, which has plotting capability and whose transforms are based on proj. My particular example is easily solved by exposing +gamma for the PROJ oblique mercator in the Cartopy interface (I just made my own class). But I'm still confused by the results of the inputs to +proj=ob_tran.

[Jeitan]

Maybe this is a bit clearer why I suspect a bug. The following is what I get for a normal Oblique Mercator (that is mimicking a Mercator) with 10 degree rotation using gamma:
+proj=omerc +lonc=0 +lat_0=0 +alpha=90 +gamma=80

Let's say I try to reproduce this. This is what I get for the following, wherein I try to make the ob_tran section actually do nothing:

+proj=ob_tran +o_proj=omerc +lonc=0 +lat_0=0 +alpha=90 +gamma=80 +o_lon_c=0 +o_lat_c=0 +o_alpha=90

Why o_alpha needs to be 90 is unclear to me since its description in the documentation is "angle to rotate the projection with", and not the angle of the centerline. Still, let's call that a conventional misunderstanding or documentation error. However, in the resulting plot the rotation is fine, but the origin is off by 90 degrees of longitude. If I "fix" that by setting o_lon_c to 90, I do get the result I was going for:

+proj=ob_tran +o_proj=omerc +lonc=0 +lat_0=0 +alpha=90 +gamma=80 +o_lon_c=90 +o_lat_c=0 +o_alpha=90

Regardless of the use-case, this origin offset seems suspicious to me, and is causing my attempts to use ob_tran to fail (e.g. to make an oblique projection from something that doesn't have a handy alpha or gamma parameter).