[gdal-dev] Define and transform a GeoTiff with a compound and projected CRS
Even Rouault
even.rouault at spatialys.com
Fri Aug 20 05:25:54 PDT 2021
Hannah,
there are 2 issues.
* 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 "NAD83(2011) /
UTM zone 12N", which is slightly more appropriate given the information
you give in introduction than EPSG:26912 which uses the older NAD83 datum)
You can get the correct WKT with:
$ gdalsrsinfo EPSG:6341+5703 -o WKT1
COMPD_CS["NAD83(2011) / UTM zone 12N + NAVD88 height",
PROJCS["NAD83(2011) / UTM zone 12N",
GEOGCS["NAD83(2011)",
DATUM["NAD83_National_Spatial_Reference_System_2011",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","1116"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","6318"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-111],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","6341"]],
VERT_CS["NAVD88 height",
VERT_DATUM["North American Vertical Datum 1988",2005,
AUTHORITY["EPSG","5103"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Gravity-related height",UP],
AUTHORITY["EPSG","5703"]]]
* 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.
I'm assuming you use PROJ 7 or later (to be able to use GeoTIFF grids)
and have the grids at
https://github.com/OSGeo/PROJ-data/tree/master/us_noaa installed in the
resource directory of PROJ (for example by unzipping the content of
http://download.osgeo.org/proj/proj-data-1.7.zip in it)
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:
COMPD_CS["NAD83(2011) / UTM zone 12N + NAVD88 height",
PROJCS["NAD83(2011) / UTM zone 12N",
GEOGCS["NAD83(2011)",
DATUM["NAD83_National_Spatial_Reference_System_2011",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","1116"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","6318"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-111],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","6341"]],
VERT_CS["NAVD88 height",
VERT_DATUM["North American Vertical Datum 1988",2005,
AUTHORITY["EPSG","5103"],
EXTENSION["PROJ4_GRIDS","us_noaa_g2012bu0.tif"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Gravity-related height",UP],
AUTHORITY["EPSG","5703"]]]
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)
So your gdalwarp invokation would be something like:
gdalwarp -overwrite -s_srs input.wkt -t_srs EPSG:4979
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.
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.
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...
Best regards,
Even
Le 20/08/2021 à 01:04, Hannah Besso a écrit :
> 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 <mailto:lumbraca at uw.edu> | (509) 823-0946
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210820/19ad53be/attachment-0001.html>
More information about the gdal-dev
mailing list