[gdal-dev] Is it possible to write custom CRS WKTs within a GeoTiff?
Even Rouault
even.rouault at spatialys.com
Wed Jan 8 08:56:21 PST 2025
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 )
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
> 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 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.
Butcher of all kinds of standards, open or closed formats. At the end, this is just about bytes.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20250108/119cc83a/attachment.htm>
More information about the gdal-dev
mailing list