2

I am trying to figure out how to do a vertical transformation using the GDAL API. In my case I'm using a C# wrapper for the GDAL API, but if answers are given in C++ or python I should be able to translate that into my project.

As an example, if I have a point in WGS84 with ellipsoid altitude at 40N, 95W, -31.82 m, then according to this site https://www.unavco.org/software/geodetic-utilities/geoid-height-calculator/geoid-height-calculator.html the orthometric height should be around 0 m.

In the GDAL API this is how I'm setting it up.

var sourceCRS = new OSGeo.OSR.SpatialReference("");
sourceCRS.ImportFromEPSG(4979); // WGS84, but specifically 3D with ellipsoidal height defined as the Z axis.
var destinationCRS = new OSGeo.OSR.SpatialReference("");
destinationCRS.ImportFromEPSG(5498); // NAD83 + NAVD88 in meters

var transformation = new OSGeo.OSR.CoordinateTransformation(sourceCRS, destinationCRS);
var point = new double[]{40, -95, -31.82};
transformation.TransformPoint(point);

After the transform, I would expect the Z value in point[2] to be close to 0, but instead the Z value is just getting passed through and is still -31.82

I've also tried using a PROJ string to setup the destination CRS like below and I downloaded the gtx files to the share directory in the GDAL folder, but this resulted in the same transformation and the Z value was just passed through as -31.82

var destinationCRS = new OSGeo.OSR.SpatialReference("");
destinationCRS.ImportFromProj4("+proj=longlat +datum=NAD83 +vunits=m +no_defs +type=crs +geoidsgrid=g2012a_conus.gtx");
TJR
  • 1,612
  • 13
  • 21
  • Maybe you need the GTX files: https://gis.stackexchange.com/questions/191565/vertical-datum-conversion-of-a-raster-with-gdal – GBG Oct 11 '22 at 18:31

1 Answers1

3

The problem ended up being two different things.

  1. The PROJ option is +geoidgrids not +geoidsgrid. The frustrating part of this realization is that if you have a typo in your PROJ string, then it will just silently fail and act as if you have no parameter set at all (at least in the GDAL.Native C# wrapper that I'm using). Once the option was typed correctly as +geoidgrids then I started to get errors indicating that the specified grid file was not getting loaded.

  2. I did download the gtx files, but I put them into a share folder location that wasn't being used where my vertical transform test was being run. In the GDAL.Native package there are paths that can be initialized like below, and environment variables can be initialized as well. The .gtx or .tif file specified needs to be in the projSharePath as initialized in the code below.

Environment.SetEnvironmentVariable("PROJ_LIB", projSharePath);
NativeGdal.SetConfigOption("PROJ_LIB", projSharePath);
Environment.SetEnvironmentVariable("PROJ_DATA", projSharePath);
NativeGdal.SetConfigOption("PROJ_DATA", projSharePath);
NativeOsr.SetPROJSearchPath(projSharePath);
TJR
  • 1,612
  • 13
  • 21