<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hannah,</p>
<p>there are 2 issues.</p>
<p>* First your WKT is indeed invalid. Instead of using a WKT, as
there are EPSG codes for the horizontal and vertical part, you
could just use gdal_translate -a_srs EPSG:6341+5703 (EPSG:6341 is
"<span>NAD83(2011) / </span><br>
<span>UTM zone 12N", which is slightly more appropriate given the
information you give in introduction than EPSG:26912 which uses
the older NAD83 datum)<br>
</span></p>
<p>You can get the correct WKT with:</p>
<p>$ gdalsrsinfo EPSG:6341+5703 -o WKT1<br>
<br>
COMPD_CS["NAD83(2011) / UTM zone 12N + NAVD88 height",<br>
PROJCS["NAD83(2011) / UTM zone 12N",<br>
GEOGCS["NAD83(2011)",<br>
DATUM["NAD83_National_Spatial_Reference_System_2011",<br>
SPHEROID["GRS 1980",6378137,298.257222101,<br>
AUTHORITY["EPSG","7019"]],<br>
AUTHORITY["EPSG","1116"]],<br>
PRIMEM["Greenwich",0,<br>
AUTHORITY["EPSG","8901"]],<br>
UNIT["degree",0.0174532925199433,<br>
AUTHORITY["EPSG","9122"]],<br>
AUTHORITY["EPSG","6318"]],<br>
PROJECTION["Transverse_Mercator"],<br>
PARAMETER["latitude_of_origin",0],<br>
PARAMETER["central_meridian",-111],<br>
PARAMETER["scale_factor",0.9996],<br>
PARAMETER["false_easting",500000],<br>
PARAMETER["false_northing",0],<br>
UNIT["metre",1,<br>
AUTHORITY["EPSG","9001"]],<br>
AXIS["Easting",EAST],<br>
AXIS["Northing",NORTH],<br>
AUTHORITY["EPSG","6341"]],<br>
VERT_CS["NAVD88 height",<br>
VERT_DATUM["North American Vertical Datum 1988",2005,<br>
AUTHORITY["EPSG","5103"]],<br>
UNIT["metre",1,<br>
AUTHORITY["EPSG","9001"]],<br>
AXIS["Gravity-related height",UP],<br>
AUTHORITY["EPSG","5703"]]]<br>
<br>
</p>
<p>* For the gdalwarp transformation, when wanting to apply vertical
transformation, the source and target CRS must be 3D, which is the
case in your invokation (once previous issue has been addressed),
but there's a strong limitation in the way gdalwarp operates in
that mode. The CRS definition that use orthometric heights (here
EPSG:6341+5703) must lso contain an explicit
EXTENSION["PROJ4_GRIDS"] node (if using WK1) or, a +geoidgrids= if
using a PROJ.4 style string.</p>
<p>I'm assuming you use PROJ 7 or later (to be able to use GeoTIFF
grids) and have the grids at
<a class="moz-txt-link-freetext" href="https://github.com/OSGeo/PROJ-data/tree/master/us_noaa">https://github.com/OSGeo/PROJ-data/tree/master/us_noaa</a> installed
in the resource directory of PROJ (for example by unzipping the
content of <a class="moz-txt-link-freetext" href="http://download.osgeo.org/proj/proj-data-1.7.zip">http://download.osgeo.org/proj/proj-data-1.7.zip</a> in it)</p>
<p>Then you can create a input.wkt file with the following content,
which is the same as above but by adding the explicit PROJ4_GRIDS
node:</p>
<p>COMPD_CS["NAD83(2011) / UTM zone 12N + NAVD88 height",<br>
PROJCS["NAD83(2011) / UTM zone 12N",<br>
GEOGCS["NAD83(2011)",<br>
DATUM["NAD83_National_Spatial_Reference_System_2011",<br>
SPHEROID["GRS 1980",6378137,298.257222101,<br>
AUTHORITY["EPSG","7019"]],<br>
AUTHORITY["EPSG","1116"]],<br>
PRIMEM["Greenwich",0,<br>
AUTHORITY["EPSG","8901"]],<br>
UNIT["degree",0.0174532925199433,<br>
AUTHORITY["EPSG","9122"]],<br>
AUTHORITY["EPSG","6318"]],<br>
PROJECTION["Transverse_Mercator"],<br>
PARAMETER["latitude_of_origin",0],<br>
PARAMETER["central_meridian",-111],<br>
PARAMETER["scale_factor",0.9996],<br>
PARAMETER["false_easting",500000],<br>
PARAMETER["false_northing",0],<br>
UNIT["metre",1,<br>
AUTHORITY["EPSG","9001"]],<br>
AXIS["Easting",EAST],<br>
AXIS["Northing",NORTH],<br>
AUTHORITY["EPSG","6341"]],<br>
VERT_CS["NAVD88 height",<br>
VERT_DATUM["North American Vertical Datum 1988",2005,<br>
AUTHORITY["EPSG","5103"],<br>
EXTENSION["PROJ4_GRIDS","us_noaa_g2012bu0.tif"]],<br>
UNIT["metre",1,<br>
AUTHORITY["EPSG","9001"]],<br>
AXIS["Gravity-related height",UP],<br>
AUTHORITY["EPSG","5703"]]]<br>
<br>
</p>
You can check with gdalsrsinfo input.wkt that it is correctly
understood (by default you'll see the WKT2 formulation of it, which
is more verbose, in particular in its expression of the
transformation between orthometric and ellipsoidal heights)<br>
<p>So your gdalwarp invokation would be something like:</p>
<p>gdalwarp -overwrite -s_srs input.wkt -t_srs EPSG:4979</p>
<p>If you add "--debug on", you'll see a trace like "GDAL:
GDALOpen(/path/to/us_noaa_g2012bu0.tif, this=......) succeeds as
GTiff." showing that the vertical shift grid is used.<br>
</p>
<p>Pedantic note: the relationship between NAD83(2011) and WGS 84 is
poorly known, at least not well captured in the entries of the
proj.db database. Basically the results you'll get will be
assuming that NAD83(2011) and WGS 84 are equivalent. Hard to say
which accuracy you'll really get in that transformation.</p>
<p>Note: I actually see that you mention ITRF2014. So a better match
for the target CRS would be EPSG:4912. "projinfo -s
EPSG:6341+5703 -t EPSG:4912 --spatial-test intersects" shows that
more accurate transformations are available than by using the
generic WGS 84 EPSG:4979. They are time dependent though, so for
best precision you'd need to specify a coordinate epoch, but I'm
probably digressing a bit too far...<br>
</p>
<p>Best regards,</p>
<p>Even<br>
</p>
<div class="moz-cite-prefix">Le 20/08/2021 à 01:04, Hannah Besso a
écrit :<br>
</div>
<blockquote type="cite"
cite="mid:CAC9j2qHfg4C5p_uH6uKYbbX677GrkKGJtXYOwd0xjNE_j2SQDg@mail.gmail.com">
<div dir="ltr">
<div>I have two elevation datasets I need to compare: one is an
airborne lidar-derived DTM (resolution 1m) and the other is
spaceborne altimeter-derived elevation points from NASA's
ICESat-2. These two datasets are in different datums and
projections. The lidar data has detailed metadata with this
info:</div>
<div><br>
</div>
<div><span>Horizontal Datum: NAD83(2011)(EPOCH:2010.0000)</span><br>
<span>Horizontal Projection: UTM Zone 12 North</span><br>
<span>Vertical Datum: NAVD88</span><br>
<span>Geoid Model: Geoid 12B</span><br>
<span>Horizontal Units: Meters</span><br>
<span>Vertical Units: Meters</span><br>
</div>
<div><span><br>
</span></div>
<div><span>While the ICESat-2 data is given as height above the
WGS84 ellipsoid (ITRF14 reference frame). </span></div>
<div><span><br>
</span></div>
<div><span><b>My ultimate goal is to get the airborne lidar into
the same CRS as the spaceborne ICESat-2 data. This
involves two steps: 1) converting the lidar data from
ASCII to GeoTiff with the correctly defined CRS, and 2)
transforming that GeoTiff to the same CRS as the ICESat-2
data. What is the best way to do this?</b></span></div>
<div><span><br>
</span></div>
<div><span>Here are the steps I've taken from my most recent
attempt to perform these tasks:</span></div>
<div><span><br>
</span></div>
<div><span>1) So far I have used gdal_translate in python to
convert from ASCII to GeoTiff, but am having issues with
defining the correct CRS. I tried using a custom WKT I wrote
using a text editor to combine two existing WKTs, but this
doesn't seem to be working. My custom WKT tries to combine
the info from EPSG:26912 and 4979, and reads:</span></div>
<div><span><br>
</span></div>
<div>COMPD_CS["NAD83 / UTM zone 12N + NAVD88 height",<br>
GEOGCS["NAD83",<br>
DATUM["North_American_Datum_1983",<br>
SPHEROID["GRS 1980",6378137,298.257222101,<br>
AUTHORITY["EPSG","7019"]],<br>
AUTHORITY["EPSG","6269"]],<br>
PRIMEM["Greenwich",0,<br>
AUTHORITY["EPSG","8901"]],<br>
UNIT["degree",0.01745329251994328,<br>
AUTHORITY["EPSG","9122"]],<br>
AUTHORITY["EPSG","4269"]],<br>
UNIT["metre",1,<br>
AUTHORITY["EPSG","9001"]],<br>
PROJECTION["Transverse_Mercator (3D)"],<br>
PARAMETER["latitude_of_origin",0],<br>
PARAMETER["central_meridian",-111],<br>
PARAMETER["scale_factor",0.9996],<br>
PARAMETER["false_easting",500000],<br>
PARAMETER["false_northing",0],<br>
AUTHORITY["EPSG","26912"],<br>
AXIS["Easting",EAST],<br>
AXIS["Northing",NORTH],<br>
VERT_CS["NAVD88 height",<br>
VERT_DATUM["North American Vertical Datum 1988",2005,<br>
AUTHORITY["EPSG","5103"],<br>
EXTENSION["PROJ4_GRIDS","g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx"]],<br>
UNIT["metre",1,<br>
AUTHORITY["EPSG","9001"]],<br>
AXIS["Up",UP],<br>
AUTHORITY["EPSG","5703"]],<br>
AUTHORITY["EPSG","5498"]]<span><br>
</span>
<div><br>
</div>
</div>
<div>The code I have been using to run gdal_translate is:</div>
<blockquote>
<div><i>!gdal_translate -of "GTiff" -a_srs '26912_5498.prj'
'0115_DTM_NoSnow.asc' '0115_DTM_NoSnow_26912_5498.tif'</i></div>
</blockquote>
<div>where '26912_5498.prj' is the name of the WKT text file
that I pasted in above. This is my first time trying to write
one of these WKT files, so I would not be surprised if there
is an error in this file. The command runs, outputting a file
named '0115_DTM_NoSnow_26912_5498.tif'. </div>
<div>Running <i>!gdalinfo '0115_DTM_NoSnow_26912_5498.tif' </i>outputs
this information, which seems to be missing info about the UTM
projection (specifically the projection parameters):</div>
<div>
<div class="gmail-lm-Widget gmail-p-Widget gmail-jp-Cell
gmail-jp-CodeCell gmail-jp-Notebook-cell gmail-jp-mod-active
gmail-jp-mod-selected">
<div class="gmail-lm-Widget gmail-p-Widget gmail-lm-Panel
gmail-p-Panel gmail-jp-Cell-outputWrapper">
<div class="gmail-lm-Widget gmail-p-Widget
gmail-jp-OutputArea gmail-jp-Cell-outputArea">
<div class="gmail-lm-Widget gmail-p-Widget
gmail-lm-Panel gmail-p-Panel
gmail-jp-OutputArea-child">
<div class="gmail-lm-Widget gmail-p-Widget
gmail-jp-RenderedText gmail-jp-OutputArea-output">
<pre>Driver: GTiff/GeoTIFF
Files: 0115_DTM_NoSnow_26912_5498.tif
Size is 3398, 3388
Coordinate System is:
COMPOUNDCRS["NAD83 / UTM zone 12N + NAVD88 height",
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
VERTCRS["NAVD88 height",
VDATUM["North American Vertical Datum 1988"],
CS[vertical,1],
AXIS["gravity-related height",up,
LENGTHUNIT["metre",1]],
ID["EPSG",5703]]]
Data axis to CRS axis mapping: 2,1,3
Origin = (579016.799999999930151,5212972.400000000372529)
Pixel Size = (0.300000000000000,-0.300000000000000)
Metadata:
AREA_OR_POINT=Point
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 579016.800, 5212972.400)
Lower Left ( 579016.800, 5211956.000)
Upper Right ( 580036.200, 5212972.400)
Lower Right ( 580036.200, 5211956.000)
Center ( 579526.500, 5212464.200)
Band 1 Block=3398x1 Type=Float64, ColorInterp=Gray
NoData Value=0
Unit Type: metre
</pre>
</div>
</div>
</div>
</div>
</div>
</div>
<div><br>
</div>
<div>2) I then used gdalwarp (again using a cell in Python) to
try to convert to WGS84 using this code:</div>
<div><i>!gdalwarp -s_srs '26912_5498.prj' -t_srs 'EPSG:4979'
'0115_DTM_NoSnow_26912_5498.tif' '0115_DTM_NoSnow_4979.tif' </i></div>
<div>Which runs, but when I output the bounds of the new GeoTif,
the coordinates appear to still be in UTM meters:</div>
<div><i>src2 = rio.open('0115_DTM_NoSnow_4979.tif')<br>
</i></div>
<div><i>src2.bounds</i><br>
</div>
<div>
<pre>BoundingBox(left=579016.7999999999, bottom=5211956.0, right=580036.2000000001, top=5212972.4)</pre>
</div>
<div>So clearly something is not working correctly. Any help or
insight is appreciated!</div>
<div>Thank you,</div>
<div>Hannah</div>
<div>-- <br>
<div dir="ltr" class="gmail_signature"
data-smartmail="gmail_signature">
<div dir="ltr">
<div><span>_____________________________________</span><br>
</div>
<div><b>Hannah Besso</b>, <i>she/her/hers</i> </div>
<div>Graduate Student, Hydrology & Hydrodynamics</div>
<div><span>Civil and Environmental Engineering</span></div>
<div><span>University of Washington | </span><span>Seattle,
WA</span><br>
</div>
<div><span><br>
</span></div>
<div><span><a href="mailto:lumbraca@uw.edu" rel="noopener
noreferrer" target="_blank" moz-do-not-send="true">bessoh2@uw.edu</a> |
(509) 823-0946</span></div>
</div>
</div>
</div>
</div>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<pre class="moz-quote-pre" wrap="">_______________________________________________
gdal-dev mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/gdal-dev">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="http://www.spatialys.com">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</body>
</html>