<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<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://github.com/opengeospatial/geotiff/issues/56">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://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh/Copernicus_DSM_10_N38_00_W109_00_DEM.tif">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://lists.osgeo.org/mailman/listinfo/gdal-dev">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="http://www.spatialys.com">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>
</body>
</html>