[PROJ] Hodine Oblique Mercator

[Even Rouault]

Andrew,

I don't reproduce the issue with proj 8.0, 8.1.1 or latest master with 
the WKT you mention (that has a rectified_grid_angle parameter).

$ projinfo 
'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]]'
PROJ.4 string:
+proj=omerc +no_uoff +lat_0=45.3091666666667 +lonc=-86 +alpha=-22.74444 
+gamma=-22.74444 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 
+datum=NAD83 +units=m +no_defs +type=crs

I assume your real original WKT lacks rectified_grid_angle, in which 
case I do reproduce the gamma = 0.

The WKT importer of PROJ (used by GDAL 3) has logic in 
https://github.com/OSGeo/PROJ/blob/master/src/iso19111/io.cpp#L3578 to 
use the azimuth values as the angle to skew grid (+gamma) when it is 
missing, which is the case in ESRI WKT. But your above WKT, which looks 
strongly like ESRI WKT, uses "Hotine_Oblique_Mercator" which isn't a 
ESRI projection name (they use different variants), so it takes a 
different road which doesn't do the same handling of missing 
rectified_grid_angle. But it should. I've just fixed this per 
https://github.com/OSGeo/PROJ/pull/2985

 > What changed from PROJ 6.1 to PROJ 8.1 that caused this?

This isn't necessarily the good way to formulate the question. PROJ 8 is 
in the direct continuity of PROJ 6, with many bugfixes and enhancements. 
For your use case, what matters more is the GDAL version. In GDAL 2.x, 
GDAL had its own logic to parse WKT and transform it to PROJ.4 string, 
even if used with PROJ 6 (by the way GDAL 2.x with PROJ 6 isn't 
necessarily a "vetted" conf. I'd say people should stick at PROJ 5.2 if 
they used GDAL 2.x). GDAL 3 requires PROJ >= 6, which has a WKT import 
logic completely re-written from scratch, , and that GDAL uses (GDAL 3 
no longer has a built-in WKT parser), so no wonder that there might he 
differences of behavior compared to GDAL 2, although all the effort and 
testing we put in trying to minimize the gaps. You seem to have a good 
dataset to shake the implementation :-)

Even

Le 16/12/2021 à 20:01, Andrew Patterson a écrit :
> 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 <mailto:andrew at avenza.com>
> phone: 416.487.5116
>
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

-- 
http://www.spatialys.com
My software is free, but my time generally not.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20211216/d55dd03d/attachment.html>