<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 style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Horizontal Datum: NAD83(2011)(EPOCH:2010.0000)</span><br style="box-sizing:inherit;color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Horizontal Projection: UTM Zone 12 North</span><br style="box-sizing:inherit;color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Vertical Datum: NAVD88</span><br style="box-sizing:inherit;color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Geoid Model: Geoid 12B</span><br style="box-sizing:inherit;color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Horizontal Units: Meters</span><br style="box-sizing:inherit;color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Vertical Units: Meters</span><br></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><br></span></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">While the ICESat-2 data is given as height above the WGS84 ellipsoid (ITRF14 reference frame). </span></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><br></span></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><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 style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><br></span></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">Here are the steps I've taken from my most recent attempt to perform these tasks:</span></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><br></span></div><div><span style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)">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 style="color:rgb(29,28,29);font-family:Slack-Lato,appleLogo,sans-serif;font-size:15px;font-variant-ligatures:common-ligatures;background-color:rgb(248,248,248)"><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"]]<font color="#1d1c1d" face="Slack-Lato, appleLogo, sans-serif"><span style="font-size:15px;font-variant-ligatures:common-ligatures"><br clear="all"></span></font><div><br></div></div><div>The code I have been using to run gdal_translate is:</div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><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 class="gmail-lm-Widget gmail-p-Widget gmail-jp-CellFooter gmail-jp-Cell-footer"></div></div><div class="gmail-lm-Widget gmail-p-Widget gmail-jp-Cell gmail-jp-CodeCell gmail-jp-Notebook-cell"><div class="gmail-lm-Widget gmail-p-Widget gmail-jp-CellHeader gmail-jp-Cell-header"></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 style="color:rgb(0,0,0)">_____________________________________</span><br></div><div><font color="#000000"><b>Hannah Besso</b>, <i>she/her/hers</i> </font></div><div><font color="#000000">Graduate Student, Hydrology & Hydrodynamics</font></div><div><span style="font-size:12.8px"><font color="#666666">Civil and Environmental Engineering</font></span></div><div><font color="#666666"><span style="font-size:12.8px">University of Washington | </span><span style="font-size:12.8px">Seattle, WA</span></font><br></div><div><span style="font-size:12.8px"><font color="#000000"><br></font></span></div><div><span style="font-size:12.8px"><font color="#000000"><a href="mailto:lumbraca@uw.edu" rel="noopener noreferrer" target="_blank">bessoh2@uw.edu</a> | (509) 823-0946</font></span></div></div></div></div></div>