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
bess...@uw.edu <lumbr...@uw.edu> | (509) 823-0946
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
