using lasvdatum to apply an arbitrary correction raster

[Andrew Ritchie]

TLDR; to vertically shift a laz pointcloud using a given input raster, do something like this:

(1)> gdalwarp -of GTX -r bilinear -overwrite -s_srs EPSG:6340 -t_srs EPSG:6318 <input_raster> <adjustment_raster.gtx>

(2)> lasvdatum -i <input_laz>* -vgrid <adjustment_raster.gtx> -olaz -odix _adj_lasvdatum

*note that the input laz file can be in a projected CS like EPSG:6340 and the output will be too, but lasvdatum works on the lat/lon coords.

The gdalwarp resampling was important (in my case) because otherwise the transformation from a projected CS to unprojected resulted in aliasing artifacts that introduced an unwanted undulation into the correction surface. I spent most of a day figuring out what was going on.

I'm correcting data using an empirically derived error model. This results in an error surface of around +/- 10cm that I want to apply to my .laz files. I can do it pretty easily on 1m gridded tiff data using gdal_calc, but I was having all sorts of problems trying to use lasvdatum to apply it to the original pointcloud, partly because the .gtx format isn't documented within gdal.

I'm posting this here so hopefully it helps others trying to do the same thing, and so I can find it in ten years when I'm stuck again, and because I've seen a few posts struggling to do this.

For me the hardest thing was figuring out that in order to get a smooth interpolation surface I need to specify a resampling method for the .gtx format raster (undocumented format in gdal).

(in my case I started with an EPSG:6340 text file so I specified the source CS (-s_srs), the result needs to be WGS84 or NAD83 lat lon, and according to NOAA the difference between the two isn't significant enough to matter given the typical size of adjustment surfaces - true in my case)

BTW, NOAA uses biquadratic interpolation for their vdatum geoid translations. The function is described here (in fortran) in a nice technical memorandum.