[gdal-dev] GDAL 3.7.0 Proj 6.3.1 strange behavior when performing NAD83(EPSG:4152) to NAD83 HARN(EPSG:4269) realization transformation

[David Klaus]

Hello,

I apologize if this is not the correct forum for this question. I am
currently trying to determine the cause of some strange behavior I've
encountered when running  a custom build of GDAL that my company uses. To
be clear I don't think the behavior is a bug on GDAL's part, but I think
that we may be doing something wrong.

First a few details. The company I work for targets GDAL 3.7.0 for our
custom build. Pertinent to this issue: we use Proj4 6.3.1 and we include in
our installation the grid files available in the proj-datumgrid-1.7 package.

So, here is a description of the strange behavior I mentioned. When I
attempt to perform the realization transformation from EPSG:4152 to
EPSG:4269 in my code, there is no change. When I attempt to perform this
same transformation using the ogr2ogr.exe built during our Custom GDAL
3.7.0 build, I see some sort of transformation applied. That said, I can't
square the change in lat/lon that ogr2ogr.exe produces with the results of
an equivalent coordinate transformation from either the epsg.io or the
https://www.ngs.noaa.gov/NCAT/ websites.

I tried to dive a bit deeper and I took a look at the information reported
when PROJ_DEBUG is set to 3. Using this I was able to verify that Proj4 is
successfully finding the HPGN grid files needed to implement the NAD83 1986
to NAD83 HARN realization transformation in both ogr2ogr.exe and in my
code. I was also able to verify that in both places the proj.db database
was found and opened.

In fact almost all of the debug printouts are identical between my code and
ogr2ogr.exe. However once it comes time to select a coordinate
transformation, my code makes the following choice:

*"PROJ: Using coordinate operation Ballpark geographic offset from NAD83 to
NAD83(HARN)"*

ogr2ogr.exe on the other hand makes this choice:



*"PROJ: Using coordinate operation NAD83 to NAD83(HARN) (3)PROJ_TRACE: NTv2
- loading grid cn      PROJ: pj_open_lib(cnhpgn.gsb): call
fopen(c:/users/dklau/appdata/roaming/carlson
software/carlson2025/icad12_1_x64/sup/GDAL370/proj-data/\cnhpgn.gsb) -
succeeded"*

If I'm not mistaken, cnhpgn.gsb is the grid file used to implement the
NAD83 1986 to NAD83 HARN transformation for Connecticut. Now the coordinate
points I am transforming are in San Francisco, California
(lon:-122.419416000000000, lat:37.774929000000000). So, the fact that
ogr2ogr.exe is using a Connecticut grid file is strange and as I mentioned
before the result (lon:-122.419417091077321 lat:37.774930105493752) does
not seem correct.

I thought that perhaps my code was correctly detecting that there were no
good grid files to cover the San Francisco transformation. But I know that
I have grid files that should apply to Pennsylvania, "pahpgn.gsb." So, I
tried a transformation for coordinate points in Philadelphia, Pennsylvania
(lon:-75.165147 lat:39.952673). However, here too my code chose to use the
"Ballpark" geographic offset method. ogr2ogr on the other hand performed a
transformation that seems very close to the correct transformation (down to
perhaps a matter of millimeters):

   - Input:
      - lon: -75.165147
      - lat: 39.952673
   - My code results "ballpark":
      - lon: -75.165147
      - lat: 39.952673
   - ogr2ogr results:
      - lon: -75.165145360841308\
      - lat: 39.952674501581825
   - NGS results:
      - lon: -75.1651458936
      - lat: 39.9526745470


Looking at the PROJ_DEBUG output in ogr2ogr.exe seems to confirm that
ogr2ogr.exe is finding and applying the pahpgn.gsb file as I expected:


*"PROJ: Using coordinate operation NAD83 to NAD83(HARN) (46)PROJ_TRACE:
NTv2 - loading grid paPROJ: pj_open_lib(pahpgn.gsb): call
fopen(c:/users/dklau/appdata/roaming/carlson
software/carlson2025/icad12_1_x64/sup/GDAL370/proj-data/\pahpgn.gsb) -
succeeded"*

So, I'm not sure exactly what I'm doing wrong. To me it appears ogr2ogr is
doing an okay job of selecting a grid file when determining how to perform
this coordinate transformation while my code is not using any grid files at
all regardless of correctness. I suspect this is due to some GDAL or PROJ
setting that I've overlooked, but I can't determine what setting could
cause this.

Hopefully this is the correct forum for this question. Any help would be
greatly appreciated. Thank you,

P.P.S. Judging from the grid files available in the proj-datumgrid-1.7
package, it looks like realization transformations such as NAD83 HARN to
NAD83 FBN (or any other NAD83 realization shift besides NAD83 1986 to NAD83
HARN) aren't available in the package we're using. Do I understand this
correctly? If so, does anyone know which version of Proj4 is needed to
correctly implement the various NAD83 realization transformations? I know
that this is not strictly a GDAL question, but I thought someone here might
know the answer.


-- 
David Klaus
Carlson Software

Disclaimer

The information contained in this communication from the sender is confidential. It is intended solely for use by the recipient and others authorized to receive it. If you are not the recipient, you are hereby notified that any disclosure, copying, distribution or taking action in relation of the contents of this information is strictly prohibited and may be unlawful.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240715/40186a78/attachment.htm>