<html><head><meta http-equiv="content-type" content="text/html; charset=utf-8"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;">Thanks Even, and appreciate pointer to the GeoTIFF 2.0 spec discussion.<div><br></div><div><blockquote type="cite">"If you can sacrifice the ellipsoidal height axis to avoid using a Projected 3D CRS, which is a bit of a geodetic dubious object (don't ask me more, but minds more enlightned than me question the soundness of using an ellipsoidal heights with a projected CRS)”</blockquote></div><div><br></div><div>I won’t ask you more :) but maybe this can be discussed further with those enlightened minds, perhaps on PROJ mailing list or a GH issue? </div><div><br></div><div>We’ve been exploring approaches to embed custom 3D CRS (both Geographic and Projected, with dynamic frame epoch definition), vertical datum information (essential), and data acquisition epoch for geotiff elevation products, with the hope of eliminating metadata ambiguity for future downstream integration and analysis. But then we’ve also mistakenly failed to bundle/transfer the aux.xml sidecar files with some geotiffs, leading to downstream issues and confusion.</div><div><br></div><div>Would be great to hear any other thoughts on best practices and options using current and future versions of GDAL/PROJ/GeoTIFF.</div><div><div>
<div dir="auto" style="caret-color: rgb(0, 0, 0); color: rgb(0, 0, 0); letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div dir="auto" style="caret-color: rgb(0, 0, 0); color: rgb(0, 0, 0); letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div dir="auto" style="caret-color: rgb(0, 0, 0); color: rgb(0, 0, 0); letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none; overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br><br><br>--</div><div>David Shean<br>Civil and Environmental Engineering<br>University of Washington<br><a href="https://www.ce.washington.edu/facultyfinder/david-shean">https://www.ce.washington.edu/facultyfinder/david-shean</a><br><br>201 More Hall, Box 352700<br>3760 E. Stevens Way NE<br>Seattle, WA 98195-2700<br>Office: (206) 543-3105, Wilcox Hall 265<br>Pronouns: he, him, his</div></div></div></div>
</div>
<div><br><blockquote type="cite"><div>On Jan 8, 2025, at 8:56 AM, Even Rouault via gdal-dev <gdal-dev@lists.osgeo.org> wrote:</div><br class="Apple-interchange-newline"><div>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<div><p>Scott,<br>
</p>
<blockquote type="cite" cite="mid:2729FB7D-0BB2-45A5-B209-BEADCEA83916@uw.edu">
<pre class="moz-quote-pre" wrap=""> I’d like to be able to store a custom WKT *inside GeoTiff metadata* rather than in an external PAM file.</pre>
</blockquote><p>GeoTIFF CRS encoding is a custom one, and doesn't support WKT
(there's some hypothetical idea of allowing/mandating WKT2 for an
hypothetical GeoTIFF 2.0 spec :
<a class="moz-txt-link-freetext" href="https://urldefense.com/v3/__https://github.com/opengeospatial/geotiff/issues/56__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5WklsndhfJJQDy10$">https://github.com/opengeospatial/geotiff/issues/56</a> )</p><p>If you can sacrifice the ellipsoidal height axis to avoid using a
Projected 3D CRS, which is a bit of a geodetic dubious object
(don't ask me more, but minds more enlightned than me question the
soundness of using an ellipsoidal heights with a projected CRS),
then the driver will make a quite reasonable effort to encode it
into GeoTIFF keys (you'll loose the comments, area of use, and
DYNAMIC[] node in the process, although one could argue that the
later should be automatically reconstructed by the GeoTIFF reader)<br>
</p><p>Note: your ending remark uses a non-ASCII <span style="white-space: pre-wrap">” character, which prevents the WKT to be parsed.</span></p><p><span style="white-space: pre-wrap">So given a test.wkt file with the following:
</span></p><p><span style="white-space: pre-wrap">PROJCRS["WGS84 (G1150) / UTM zone 12N",
BASEGEOGCRS["WGS 84 (G1150)",
DYNAMIC[
FRAMEEPOCH[2001]],
DATUM["World Geodetic System 1984 (G1150)",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",7661]],
CONVERSION["UTM zone 12N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-111,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16012]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
USAGE[
SCOPE["unknown"],
AREA["Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],
BBOX[0,-114,84,-108]],
REMARK["Custom WKT2 combining UTM_12N with WGS84 (G1150)"]]
</span></p><p>$ gdal_translate in.tif out.tif -a_srs "$(cat test.wkt)"</p><p>$ gdalinfo out.tif</p><p>PROJCRS["WGS84 (G1150) / UTM zone 12N",<br>
BASEGEOGCRS["WGS 84 (G1150)",<br>
DATUM["World Geodetic System 1984 (G1150)",<br>
ELLIPSOID["WGS 84",6378137,298.257223563,<br>
LENGTHUNIT["metre",1]]],<br>
PRIMEM["Greenwich",0,<br>
ANGLEUNIT["degree",0.0174532925199433]],<br>
ID["EPSG",7661]],<br>
CONVERSION["Transverse Mercator",<br>
METHOD["Transverse Mercator",<br>
ID["EPSG",9807]],<br>
PARAMETER["Latitude of natural origin",0,<br>
ANGLEUNIT["degree",0.0174532925199433],<br>
ID["EPSG",8801]],<br>
PARAMETER["Longitude of natural origin",-111,<br>
ANGLEUNIT["degree",0.0174532925199433],<br>
ID["EPSG",8802]],<br>
PARAMETER["Scale factor at natural origin",0.9996,<br>
SCALEUNIT["unity",1],<br>
ID["EPSG",8805]],<br>
PARAMETER["False easting",500000,<br>
LENGTHUNIT["metre",1],<br>
ID["EPSG",8806]],<br>
PARAMETER["False northing",0,<br>
LENGTHUNIT["metre",1],<br>
ID["EPSG",8807]]],<br>
CS[Cartesian,2],<br>
AXIS["easting",east,<br>
ORDER[1],<br>
LENGTHUNIT["metre",1,<br>
ID["EPSG",9001]]],<br>
AXIS["northing",north,<br>
ORDER[2],<br>
LENGTHUNIT["metre",1,<br>
ID["EPSG",9001]]]]</p><p>$ listgeo out.tif<br>
[... snip ... ]<br>
Keyed_Information:<br>
GTModelTypeGeoKey (Short,1): ModelTypeProjected<br>
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea<br>
GTCitationGeoKey (Ascii,29): "WGS84 (G1150) / UTM zone 12N"<br>
GeographicTypeGeoKey (Short,1): Code-7661 (WGS 84 (G1150))<br>
GeogCitationGeoKey (Ascii,15): "WGS 84 (G1150)"<br>
GeogAngularUnitsGeoKey (Short,1): Angular_Degree<br>
GeogSemiMajorAxisGeoKey (Double,1): 6378137 <br>
GeogInvFlatteningGeoKey (Double,1): 298.257223563 <br>
ProjectedCSTypeGeoKey (Short,1): User-Defined<br>
ProjectionGeoKey (Short,1): Proj_UTM_zone_12N<br>
ProjLinearUnitsGeoKey (Short,1): Linear_Meter<br>
End_Of_Keys.<br>
End_Of_Geotiff.<br>
<br>
<br>
</p>
<blockquote type="cite" cite="mid:2729FB7D-0BB2-45A5-B209-BEADCEA83916@uw.edu">
<pre class="moz-quote-pre" wrap="">
Here is a motivating example, using a WKT that specifies a specific WGS84 realization for a given UTM zone:
```
INPUT=<a class="moz-txt-link-freetext" href="https://urldefense.com/v3/__https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh/Copernicus_DSM_10_N38_00_W109_00_DEM.tif__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5WklsndhfYx7Qu_8$">https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh/Copernicus_DSM_10_N38_00_W109_00_DEM.tif</a>
OUTPUT=/tmp/Copernicus_DSM_10_N38_00_W109_00_DEM_UTM.tif
CPL_DEBUG=ON \
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
gdalwarp -co GEOTIFF_VERSION=1.1 -overwrite -dstnodata 0 -r bilinear -tap -tr 30.0 30.0 -s_srs EPSG:7661 -t_srs utm12N_wgs1150.wkt $INPUT $OUTPUT
```
Here is the wkt, I’d like for it not to end up in a sidecar file where I might loose it, but embedded in the GTiff:
```
PROJCRS["WGS84 (G1150) / UTM zone 12N",
BASEGEOGCRS["WGS 84 (G1150)",
DYNAMIC[
FRAMEEPOCH[2001]],
DATUM["World Geodetic System 1984 (G1150)",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",7661]],
CONVERSION["UTM zone 12N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-111,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16012]],
CS[Cartesian,3],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["ellipsoidal height (h)",up,
ORDER[3],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
USAGE[
SCOPE["unknown"],
AREA["Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],
BBOX[0,-114,84,-108]],
REMARK["Custom WKT2 combining UTM_12N with WGS84 (G1150)”]]
```
Cheers,
Scott
_______________________________________________
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://urldefense.com/v3/__https://lists.osgeo.org/mailman/listinfo/gdal-dev__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5Wklsndhflk8ZnOQ$">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="https://urldefense.com/v3/__http://www.spatialys.com__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5WklsndhfQhOfbpw$">http://www.spatialys.com</a>
My software is free, but my time generally not.
Butcher of all kinds of standards, open or closed formats. At the end, this is just about bytes.</pre>
</div>
_______________________________________________<br>gdal-dev mailing list<br>gdal-dev@lists.osgeo.org<br>https://urldefense.com/v3/__https://lists.osgeo.org/mailman/listinfo/gdal-dev__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5Wklsndhflk8ZnOQ$ <br></div></blockquote></div><br></div></body></html>