[gdal-dev] Define and transform a GeoTiff with a compound and projected CRS

[Hannah Besso]

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:

Horizontal Datum: NAD83(2011)(EPOCH:2010.0000)
Horizontal Projection: UTM Zone 12 North
Vertical Datum: NAVD88
Geoid Model: Geoid 12B
Horizontal Units: Meters
Vertical Units: Meters

While the ICESat-2 data is given as height above the WGS84 ellipsoid
(ITRF14 reference frame).

*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?*

Here are the steps I've taken from my most recent attempt to perform these
tasks:

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:

COMPD_CS["NAD83 / UTM zone 12N + NAVD88 height",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Transverse_Mercator (3D)"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-111],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","26912"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    VERT_CS["NAVD88 height",
        VERT_DATUM["North American Vertical Datum 1988",2005,
            AUTHORITY["EPSG","5103"],

EXTENSION["PROJ4_GRIDS","g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","5703"]],
    AUTHORITY["EPSG","5498"]]

The code I have been using to run gdal_translate is:

*!gdal_translate -of "GTiff" -a_srs '26912_5498.prj' '0115_DTM_NoSnow.asc'
'0115_DTM_NoSnow_26912_5498.tif'*

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'.
Running *!gdalinfo '0115_DTM_NoSnow_26912_5498.tif' *outputs this
information, which seems to be missing info about the UTM projection
(specifically the projection parameters):

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


2) I then used gdalwarp (again using a cell in Python) to try to convert to
WGS84 using this code:
*!gdalwarp -s_srs '26912_5498.prj' -t_srs 'EPSG:4979'
'0115_DTM_NoSnow_26912_5498.tif' '0115_DTM_NoSnow_4979.tif' *
Which runs, but when I output the bounds of the new GeoTif, the coordinates
appear to still be in UTM meters:

*src2 = rio.open('0115_DTM_NoSnow_4979.tif')*
*src2.bounds*

BoundingBox(left=579016.7999999999, bottom=5211956.0,
right=580036.2000000001, top=5212972.4)

So clearly something is not working correctly. Any help or insight is
appreciated!
Thank you,
Hannah
-- 
_____________________________________
*Hannah Besso*,  *she/her/hers*
Graduate Student, Hydrology & Hydrodynamics
Civil and Environmental Engineering
University of Washington | Seattle, WA

bessoh2 at uw.edu <lumbraca at uw.edu> | (509) 823-0946
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210819/50d12a6e/attachment.html>