GDAL Processing for Vertical Shifts - granitehills14/Airborne-Laser-Scanning GitHub Wiki
These are Shad O'Neel's notes from a 3 hour power-session with Mike Smith.
gdalwarp -of vrt <source> <dest.vrt> -s_srs epsg:32611+5703 -t_srs epsg:32611-r bilinear
The VRT driver is a format driver for GDAL that allows a virtual GDAL dataset to be composed from other GDAL datasets with repositioning, and algorithms potentially applied as well as various kinds of metadata altered or added. VRT descriptions of datasets can be saved in an XML format normally given the extension *.vrt. https://gdal.org/drivers/raster/vrt.html?highlight=vrt
Note that using gdalwarp to build the vrt is not interchangeable with gdalbuildvrt that is useful if you want to seamlessly merge several rasters.
This is telling gdal translate (step 3) to reference the “Vertical Shift Grids”. The last line in the example below is the referenced closing tag. You paste above it.
<VerticalShiftGrids>
<Grids>g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx</Grids>
<Inverse>False</Inverse>
<ToMeterSrc>1</ToMeterSrc>
<ToMeterDest>1</ToMeterDest>
</VerticalShiftGrids>
</VRTDataset>
The line <Inverse>False</Inverse>
is telling gdal that you are going from geoid to ellipsoid. To go ellipsoid to geoid it should be <Inverse>True</Inverse>
.
gdal_translate -of <source.vrt> <output.tif>
Checking to see if it worked can be done using
GTIFF_REPORT_COMPD_CS=TRUE gdalinfo DryCreek_WGS84.tif -stats
before and after the transform, and the min and max will change, following the vertical shift you just executed.
Other notes: WGS 84 is assumed, you cant specify it for vertical datum, or gdal throws an error :
ERROR 1: SetCompoundCS() fails, vertical component is not VERT_CS.
If the WGS84 ellipsoid is your desired datum, no EPSG is specified. When editing the vrt (I used an editor called Atom) errors on the Bash prompt have line numbers that refer to the vrt. This is helpful, but I didn’t see the link right away. Setting up the environment for this to work is trickier, gdal is using proj to do the transforms. But I did not have the vertical shift grids installed. These can be installed easily using projsync:
projsync --source-id us_nga && projsync --source-id us_noaa
projsync is a program that downloads remote resource files into a local directory. This is an alternative to downloading a proj-data-X.Y.Z archive file, or using the on-demand networking capabilities of PROJ. You may have to make sure they are going to the right place. Run export
and look where the proj_lib directory is. This is where the files you just downloaded should end up. Mine, for example was:
/Users/rdcrlsro/opt/anaconda3/envs/myenv/share/proj
I had to set this variable :
XDG_DATA_HOME=/Users/rdcrlsro/opt/anaconda3/envs/myenv/share
export XDG_DATA_HOME=/Users/rdcrlsro/opt/anaconda3/envs/myenv/share
$XDG_DATA_HOME
defines the base directory relative to which user-specific data files should be stored. If $XDG_DATA_HOME
is either not set or empty, a default equal to $HOME/.local/share
should be used.Then run the projsync command again.
I still had issues with .gtx files being needed, but not grabbed during the download. These are available here.
But the naming convention has changed compared to what was entered in the vrt file. You can change the naming convention in the vrt from g2012a_conus.gtx
to g2012b#
I was given the files in the original naming, which fixed things. A good debugging step here was to :
gdalinfo $PROJ_LIB/g2012a_conus.gtx
And see that some metadata was returned.