[PROJ] Hodine Oblique Mercator

[Andrew Patterson]

Hello again!

I'm afraid I'm back with another GDAL 2.x -> GDAL 3.x mystery. I've
unraveled some of it, but I'm hoping someone can illuminate the rest and
maybe even suggest a course of action.

We have existing maps in our map store that are using the
Hotine_Oblique_Mercator projection. They worked fine in GDAL 2.2 (PROJ 6.1)
but are not working correctly in GDAL 3.3 (PROJ 8.1). They load and
transform, but whereas the map is of Michigan, with the WKT below it thinks
the map is in New Mexico. The WKT in question is:

PROJCS["NAD_1983_Michigan_GeoRef_Meters",GEOGCS["GCS_NAD83(1986)",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER["false_easting",2546731.496],PARAMETER["false_northing",-4354009.816],PARAMETER["latitude_of_center",45.30916666666666],PARAMETER["longitude_of_center",-86],PARAMETER["azimuth",-22.74444],PARAMETER["rectified_grid_angle",-22.74444],PARAMETER["scale_factor",0.9996],UNIT["Meter",1]]

Having learned something from the last time this list assisted me, I
checked the proj4 string in 2.4 vs 3.3 and it revealed a difference:

GDAL 2.4/ PROJ 6.1:
+proj=omerc +no_uoff +lat_0=45.30916666666666 +lonc=-86 +alpha=-22.74444
+k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +datum=NAD83 +units=m +no_defs

GDAL 3.3/PROJ 8.1:
+proj=omerc +no_uoff +lat_0=45.3091666666667 +lonc=-86 +alpha=-22.74444
+k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +datum=NAD83 +units=m +no_defs
+gamma=0

The difference is, clearly, the "+gamma=0". If I remove that and import the
string, instead of the above WKT, it works. If I then export it to WKT, it
adds a parameter:

PARAMETER["rectified_grid_angle",-22.74444]

This appears to be set to the same value as the "azimuth" parameter. So
currently I have a hack where if I see a "Hotine_Oblique_Mercator"
projection that has an "azimuth" parameter but no "rectified_grid_angle", I
add "rectified_grid_angle" and set its value to the same as the "azimuth"
parameter. This feels a bit sketchy, but maybe that is the best solution.

So I guess what I'm asking is: IS that the best solution? What changed from
PROJ 6.1 to PROJ 8.1 that caused this? I need to support these older maps
somehow, so I'd love some guidance as to what a 'safe' hack is if possible.

Any help would be greatly appreciated!

..............................
Andrew Patterson
Lead Software Architect
Avenza Systems Inc.

email: andrew at avenza.com
phone: 416.487.5116
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211216/557390cf/attachment.html>