diff --git a/README.md b/README.md index e7ec72e..b124c1e 100644 --- a/README.md +++ b/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 ``` diff --git a/python/DEVELOP.md b/python/DEVELOP.md index 3b97be2..30d0ddc 100644 --- a/python/DEVELOP.md +++ b/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 ``` diff --git a/python/src/geo.rs b/python/src/geo.rs index ef33e71..e4611cc 100644 --- a/python/src/geo.rs +++ b/python/src/geo.rs @@ -37,7 +37,10 @@ pub(crate) struct PyGeoKeyDirectory { geog_azimuth_units: Option, #[pyo3(get)] geog_prime_meridian_long: Option, - + #[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>, #[pyo3(get)] projected_type: Option, #[pyo3(get)] @@ -117,6 +120,7 @@ impl From 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 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, diff --git a/src/geo/geo_key_directory.rs b/src/geo/geo_key_directory.rs index acf9aa2..e341b98 100644 --- a/src/geo/geo_key_directory.rs +++ b/src/geo/geo_key_directory.rs @@ -32,6 +32,7 @@ pub(crate) enum GeoKeyTag { GeogInvFlattening = 2059, GeogAzimuthUnits = 2060, GeogPrimeMeridianLong = 2061, + GeogToWGS84 = 2062, // Projected CRS Parameter Keys ProjectedType = 3072, @@ -101,6 +102,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, + pub geog_to_wgs84: Option>, pub projected_type: Option, pub proj_citation: Option, @@ -154,6 +156,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; @@ -206,6 +209,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()?), @@ -265,6 +269,7 @@ impl GeoKeyDirectory { geog_inv_flattening, geog_azimuth_units, geog_prime_meridian_long, + geog_to_wgs84, projected_type, proj_citation, diff --git a/tests/geogtowgs_tiff.rs b/tests/geogtowgs_tiff.rs new file mode 100644 index 0000000..94e3142 --- /dev/null +++ b/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 = + Arc::new(LocalFileSystem::new_with_prefix(folder).unwrap()); + let reader = Arc::new(ObjectReader::new(store, path)) as Arc; + 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(); + } +} diff --git a/tests/images/geogtowgs_subset_USGS_13_s14w171.tif b/tests/images/geogtowgs_subset_USGS_13_s14w171.tif new file mode 100644 index 0000000..1f21e29 Binary files /dev/null and b/tests/images/geogtowgs_subset_USGS_13_s14w171.tif differ diff --git a/tests/images/readme.md b/tests/images/readme.md new file mode 100644 index 0000000..94efdd0 --- /dev/null +++ b/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" +``` \ No newline at end of file