Add support for GeogTOWGS84GeoKey

README.md

@@ -24,4 +24,5 @@ Download the following file for use in the tests.
 
 ```shell
 aws s3 cp s3://naip-visualization/ny/2022/60cm/rgb/40073/m_4007307_sw_18_060_20220803.tif ./ --request-payer
+aws s3 cp s3://prd-tnm/StagedProducts/Elevation/13/TIFF/current/s14w171/USGS_13_s14w171.tif ./ --no-sign-request --region us-west-2
 ```

python/DEVELOP.md

@@ -2,4 +2,5 @@
 uv sync --no-install-package async-tiff
 uv run --no-project maturin develop
 uv run --no-project mkdocs serve
+uv run --no-project pytest --verbose
 ```

python/src/geo.rs

@@ -37,7 +37,10 @@ pub(crate) struct PyGeoKeyDirectory {
     geog_azimuth_units: Option<u16>,
     #[pyo3(get)]
     geog_prime_meridian_long: Option<f64>,
-
+    #[pyo3(get)]
+    /// Note: This is not part of the official GeoTIFF specification but is a proposed extension.
+    /// See https://trac.osgeo.org/geotiff/wiki/TOWGS84GeoKey.
+    geog_to_wgs84: Option<Vec<f64>>,
     #[pyo3(get)]
     projected_type: Option<u16>,
     #[pyo3(get)]
@@ -117,6 +120,7 @@ impl From<PyGeoKeyDirectory> for GeoKeyDirectory {
             geog_inv_flattening: value.geog_inv_flattening,
             geog_azimuth_units: value.geog_azimuth_units,
             geog_prime_meridian_long: value.geog_prime_meridian_long,
+            geog_to_wgs84: value.geog_to_wgs84,
             projected_type: value.projected_type,
             proj_citation: value.proj_citation,
             projection: value.projection,
@@ -169,6 +173,7 @@ impl From<GeoKeyDirectory> for PyGeoKeyDirectory {
             geog_inv_flattening: value.geog_inv_flattening,
             geog_azimuth_units: value.geog_azimuth_units,
             geog_prime_meridian_long: value.geog_prime_meridian_long,
+            geog_to_wgs84: value.geog_to_wgs84,
             projected_type: value.projected_type,
             proj_citation: value.proj_citation,
             projection: value.projection,

src/geo/geo_key_directory.rs

@@ -31,6 +31,7 @@ pub(crate) enum GeoKeyTag {
     GeogInvFlattening = 2059,
     GeogAzimuthUnits = 2060,
     GeogPrimeMeridianLong = 2061,
+    GeogToWGS84 = 2062,
 
     // Projected CRS Parameter Keys
     ProjectedType = 3072,
@@ -100,6 +101,7 @@ pub struct GeoKeyDirectory {
     /// defined by its longitude relative to the international reference meridian (for the earth
     /// this is Greenwich).
     pub geog_prime_meridian_long: Option<f64>,
+    pub geog_to_wgs84: Option<Vec<f64>>,
 
     pub projected_type: Option<u16>,
     pub proj_citation: Option<String>,
@@ -153,6 +155,7 @@ impl GeoKeyDirectory {
         let mut geog_inv_flattening = None;
         let mut geog_azimuth_units = None;
         let mut geog_prime_meridian_long = None;
+        let mut geog_to_wgs84 = None;
 
         let mut projected_type = None;
         let mut proj_citation = None;
@@ -205,6 +208,7 @@ impl GeoKeyDirectory {
                 GeoKeyTag::GeogPrimeMeridianLong => {
                     geog_prime_meridian_long = Some(value.into_f64()?)
                 }
+                GeoKeyTag::GeogToWGS84 => geog_to_wgs84 = Some(value.into_f64_vec()?),
                 GeoKeyTag::ProjectedType => projected_type = Some(value.into_u16()?),
                 GeoKeyTag::ProjCitation => proj_citation = Some(value.into_string()?),
                 GeoKeyTag::Projection => projection = Some(value.into_u16()?),
@@ -264,6 +268,7 @@ impl GeoKeyDirectory {
             geog_inv_flattening,
             geog_azimuth_units,
             geog_prime_meridian_long,
+            geog_to_wgs84,
 
             projected_type,
             proj_citation,

tests/geogtowgs_tiff.rs

@@ -0,0 +1,31 @@
+#[cfg(test)]
+mod test {
+    use std::env;
+    use std::sync::Arc;
+
+    use async_tiff::metadata::{PrefetchBuffer, TiffMetadataReader};
+    use async_tiff::reader::{AsyncFileReader, ObjectReader};
+
+    use object_store::local::LocalFileSystem;
+
+    #[tokio::test]
+    async fn tmp_towgs84() {
+        let folder = env::current_dir().unwrap();
+        let path =
+            object_store::path::Path::parse("tests/images/geogtowgs_subset_USGS_13_s14w171.tif")
+                .unwrap();
+        let store: Arc<LocalFileSystem> =
+            Arc::new(LocalFileSystem::new_with_prefix(folder).unwrap());
+        let reader = Arc::new(ObjectReader::new(store, path)) as Arc<dyn AsyncFileReader>;
+        let prefetch_reader = PrefetchBuffer::new(reader.clone(), 32 * 1024)
+            .await
+            .unwrap();
+        let mut metadata_reader = TiffMetadataReader::try_open(&prefetch_reader)
+            .await
+            .unwrap();
+        let _ = metadata_reader
+            .read_all_ifds(&prefetch_reader)
+            .await
+            .unwrap();
+    }
+}

tests/images/readme.md

@@ -0,0 +1,35 @@
+`geogtowgs_subset_USGS_13_s14w171.tif` was created from "s3://prd-tnm/StagedProducts/Elevation/13/TIFF/current/s14w171/USGS_13_s14w171.tif" using these commands:
+
+```bash
+gdal_translate USGS_13_s14w171.tif tiny.tif -srcwin 0 0 1 1 -co COMPRESS=DEFLATE
+listgeo USGS_13_s14w171.tif > metadata.txt  # Then modify to remove information related to spatial extent
+cp tiny.tif geogtowgs_subset_USGS_13_s14w171.tif
+geotifcp -g metadata.txt tiny.tif geogtowgs_subset_USGS_13_s14w171.tif
+listgeo geogtowgs_subset_USGS_13_s14w171.tif
+```
+
+and this workspace definition:
+
+```toml
+[project]
+name = "gdal-workspace"
+version = "0.1.0"
+description = "workspace for using gdal via pixi"
+readme = "README.md"
+requires-python = ">=3.12"
+dependencies = []
+
+[tool.pixi.workspace]
+channels = ["conda-forge"]
+platforms = ["osx-arm64"]
+
+[tool.pixi.pypi-dependencies]
+gdal-workspace = { path = ".", editable = true }
+
+[tool.pixi.tasks]
+
+[tool.pixi.dependencies]
+gdal = ">=3.11.5,<4"
+libgdal = ">=3.11.5,<4"
+geotiff = ">=1.7.4,<2"
+```