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@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.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
