Thanks Even, and appreciate pointer to the GeoTIFF 2.0 spec discussion. > "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)”
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? 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. Would be great to hear any other thoughts on best practices and options using current and future versions of GDAL/PROJ/GeoTIFF. -- David Shean Civil and Environmental Engineering University of Washington https://www.ce.washington.edu/facultyfinder/david-shean 201 More Hall, Box 352700 3760 E. Stevens Way NE Seattle, WA 98195-2700 Office: (206) 543-3105, Wilcox Hall 265 Pronouns: he, him, his > On Jan 8, 2025, at 8:56 AM, Even Rouault via gdal-dev > <gdal-dev@lists.osgeo.org> wrote: > > Scott, > >> I’d like to be able to store a custom WKT *inside GeoTiff metadata* rather >> than in an external PAM file. > 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 : https://github.com/opengeospatial/geotiff/issues/56 > <https://urldefense.com/v3/__https://github.com/opengeospatial/geotiff/issues/56__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5WklsndhfJJQDy10$> > ) > > 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) > > Note: your ending remark uses a non-ASCII ” character, which prevents the WKT > to be parsed. > > So given a test.wkt file with the following: > > 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)"]] > > $ gdal_translate in.tif out.tif -a_srs "$(cat test.wkt)" > > $ gdalinfo out.tif > > PROJCRS["WGS84 (G1150) / UTM zone 12N", > BASEGEOGCRS["WGS 84 (G1150)", > 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["Transverse Mercator", > 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]]], > CS[Cartesian,2], > AXIS["easting",east, > ORDER[1], > LENGTHUNIT["metre",1, > ID["EPSG",9001]]], > AXIS["northing",north, > ORDER[2], > LENGTHUNIT["metre",1, > ID["EPSG",9001]]]] > > $ listgeo out.tif > [... snip ... ] > Keyed_Information: > GTModelTypeGeoKey (Short,1): ModelTypeProjected > GTRasterTypeGeoKey (Short,1): RasterPixelIsArea > GTCitationGeoKey (Ascii,29): "WGS84 (G1150) / UTM zone 12N" > GeographicTypeGeoKey (Short,1): Code-7661 (WGS 84 (G1150)) > GeogCitationGeoKey (Ascii,15): "WGS 84 (G1150)" > GeogAngularUnitsGeoKey (Short,1): Angular_Degree > GeogSemiMajorAxisGeoKey (Double,1): 6378137 > GeogInvFlatteningGeoKey (Double,1): 298.257223563 > ProjectedCSTypeGeoKey (Short,1): User-Defined > ProjectionGeoKey (Short,1): Proj_UTM_zone_12N > ProjLinearUnitsGeoKey (Short,1): Linear_Meter > End_Of_Keys. > End_Of_Geotiff. > > > >> >> Here is a motivating example, using a WKT that specifies a specific WGS84 >> realization for a given UTM zone: >> ``` >> INPUT=https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh/Copernicus_DSM_10_N38_00_W109_00_DEM.tif >> >> <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$> >> 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 >> gdal-dev@lists.osgeo.org <mailto:gdal-dev@lists.osgeo.org> >> https://lists.osgeo.org/mailman/listinfo/gdal-dev >> <https://urldefense.com/v3/__https://lists.osgeo.org/mailman/listinfo/gdal-dev__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5Wklsndhflk8ZnOQ$> > -- > http://www.spatialys.com > <https://urldefense.com/v3/__http://www.spatialys.com__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5WklsndhfQhOfbpw$> > 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. > _______________________________________________ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > https://urldefense.com/v3/__https://lists.osgeo.org/mailman/listinfo/gdal-dev__;!!K-Hz7m0Vt54!l1HH4-dmmtyU4CPsIRtDt1Zf3CKzVUwtIc5gEfUTGVSyawjmrWKSW2tTag1LjZRzNGr5Wklsndhflk8ZnOQ$
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
