From 1b8fac8b1b6030d48f126d130b66b2538a484fa8 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Wed, 13 Mar 2024 17:42:31 +1300 Subject: [PATCH] :sparkles: Get affine transformation from GeoTIFF (#8) * :heavy_plus_sign: Add geo Geospatial primitives and algorithms! * :sparkles: Get affine transformation from GeoTIFF Implement transform method in CogReader struct to get affine transformation matrix from GeoTIFF via ModelPixelScaleTag and ModelTiepointTag. Added a unit test to check that the x/y pixel size and top left origin pixel coordinates is correct. * :memo: Document transform method's internal properties Explicitly mentioning the internal fields of the AffineTransform struct, specifically the pixel resolution, rotation and origin, and included a reference to the GeoTIFF standard. Also renamed origin_x/origin_y to x_origin/y_origin. * :construction: Raise unimplemented when getting non-zero rotation values Supposed to be parsed from the ModelTransformationTag, but haven't found a sample GeoTIFF online with non-zero rotation to implement this properly with a unit test, so marking as unimplemented for now. --- Cargo.lock | 190 +++++++++++++++++++++++++++++++++++++++++++++- Cargo.toml | 1 + src/io/geotiff.rs | 80 ++++++++++++++++++- 3 files changed, 267 insertions(+), 4 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 87b1319..4b095ec 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -17,6 +17,24 @@ version = "1.0.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f26201604c87b1e01bd3d98f8d5d9a8fcbb815e8cedb41ffccbeb4bf593a35fe" +[[package]] +name = "ahash" +version = "0.8.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e89da841a80418a9b391ebaea17f5c112ffaaa96f621d2c285b5174da76b9011" +dependencies = [ + "cfg-if", + "once_cell", + "version_check", + "zerocopy", +] + +[[package]] +name = "allocator-api2" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0942ffc6dcaadf03badf6e6a2d0228460359d5e34b57ccdc720b7382dfbd5ec5" + [[package]] name = "android-tzdata" version = "0.1.1" @@ -32,6 +50,15 @@ dependencies = [ "libc", ] +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + [[package]] name = "async-trait" version = "0.1.77" @@ -88,6 +115,12 @@ version = "3.15.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8ea184aa71bb362a1157c896979544cc23974e08fd265f29ea96b59f0b4a555b" +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + [[package]] name = "bytes" version = "1.5.0" @@ -124,6 +157,7 @@ name = "cog3pio" version = "0.1.0" dependencies = [ "bytes", + "geo", "ndarray", "numpy", "object_store", @@ -165,6 +199,16 @@ version = "0.3.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "fea41bba32d969b513997752735605054bc0dfa92b4c56bf1189f2e174be7a10" +[[package]] +name = "earcutr" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "79127ed59a85d7687c409e9978547cffb7dc79675355ed22da6b66fd5f6ead01" +dependencies = [ + "itertools 0.11.0", + "num-traits", +] + [[package]] name = "either" version = "1.10.0" @@ -212,6 +256,12 @@ dependencies = [ "miniz_oxide", ] +[[package]] +name = "float_next_after" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8bf7cc16383c4b8d58b9905a8509f02926ce3058053c056376248d958c9df1e8" + [[package]] name = "fnv" version = "1.0.7" @@ -316,6 +366,44 @@ dependencies = [ "slab", ] +[[package]] +name = "geo" +version = "0.28.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f811f663912a69249fa620dcd2a005db7254529da2d8a0b23942e81f47084501" +dependencies = [ + "earcutr", + "float_next_after", + "geo-types", + "geographiclib-rs", + "log", + "num-traits", + "robust", + "rstar", + "spade", +] + +[[package]] +name = "geo-types" +version = "0.7.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ff16065e5720f376fbced200a5ae0f47ace85fd70b7e54269790281353b6d61" +dependencies = [ + "approx", + "num-traits", + "rstar", + "serde", +] + +[[package]] +name = "geographiclib-rs" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6e5ed84f8089c70234b0a8e0aedb6dc733671612ddc0d37c6066052f9781960" +dependencies = [ + "libm", +] + [[package]] name = "getrandom" version = "0.2.12" @@ -352,11 +440,34 @@ dependencies = [ "tracing", ] +[[package]] +name = "hash32" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "47d60b12902ba28e2730cd37e95b8c9223af2808df9e902d4df49588d1470606" +dependencies = [ + "byteorder", +] + [[package]] name = "hashbrown" version = "0.14.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "290f1a1d9242c78d09ce40a5e87e7554ee637af1351968159f4952f028f75604" +dependencies = [ + "ahash", + "allocator-api2", +] + +[[package]] +name = "heapless" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0bfb9eb618601c89945a70e254898da93b13be0388091d42117462b265bb3fad" +dependencies = [ + "hash32", + "stable_deref_trait", +] [[package]] name = "heck" @@ -503,6 +614,15 @@ version = "2.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8f518f335dce6725a761382244631d86cf0ccb2863413590b31338feb467f9c3" +[[package]] +name = "itertools" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b1c173a5686ce8bfa551b3563d0c2170bf24ca44da99c7ca4bfdab5418c3fe57" +dependencies = [ + "either", +] + [[package]] name = "itertools" version = "0.12.1" @@ -539,6 +659,12 @@ version = "0.2.153" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9c198f91728a82281a64e1f4f9eeb25d82cb32a5de251c6bd1b5154d63a8e7bd" +[[package]] +name = "libm" +version = "0.2.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" + [[package]] name = "linux-raw-sys" version = "0.4.13" @@ -650,6 +776,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "da0df0e5185db44f69b44f26786fe401b6c293d1907744beaa7fa62b2e5a517a" dependencies = [ "autocfg", + "libm", ] [[package]] @@ -699,7 +826,7 @@ dependencies = [ "futures", "humantime", "hyper", - "itertools", + "itertools 0.12.1", "parking_lot", "percent-encoding", "quick-xml", @@ -974,6 +1101,23 @@ dependencies = [ "windows-sys 0.52.0", ] +[[package]] +name = "robust" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cbf4a6aa5f6d6888f39e980649f3ad6b666acdce1d78e95b8a2cb076e687ae30" + +[[package]] +name = "rstar" +version = "0.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "133315eb94c7b1e8d0cb097e5a710d850263372fd028fff18969de708afc7008" +dependencies = [ + "heapless", + "num-traits", + "smallvec", +] + [[package]] name = "rustc-demangle" version = "0.1.23" @@ -1195,12 +1339,30 @@ dependencies = [ "windows-sys 0.52.0", ] +[[package]] +name = "spade" +version = "2.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61addf9117b11d1f5b4bf6fe94242ba25f59d2d4b2080544b771bd647024fd00" +dependencies = [ + "hashbrown", + "num-traits", + "robust", + "smallvec", +] + [[package]] name = "spin" version = "0.9.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6980e8d7511241f8acf4aebddbb1ff938df5eebe98691418c4468d0b72a96a67" +[[package]] +name = "stable_deref_trait" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3" + [[package]] name = "syn" version = "1.0.109" @@ -1433,6 +1595,12 @@ dependencies = [ "percent-encoding", ] +[[package]] +name = "version_check" +version = "0.9.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "49874b5167b65d7193b8aba1567f5c7d93d001cafc34600cee003eda787e483f" + [[package]] name = "walkdir" version = "2.5.0" @@ -1734,3 +1902,23 @@ dependencies = [ "cfg-if", "windows-sys 0.48.0", ] + +[[package]] +name = "zerocopy" +version = "0.7.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "74d4d3961e53fa4c9a25a8637fc2bfaf2595b3d3ae34875568a5cf64787716be" +dependencies = [ + "zerocopy-derive", +] + +[[package]] +name = "zerocopy-derive" +version = "0.7.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ce1b18ccd8e73a9321186f97e46f9f04b778851177567b1975109d26a08d2a6" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.50", +] diff --git a/Cargo.toml b/Cargo.toml index c07223b..9108a39 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,6 +11,7 @@ crate-type = ["cdylib"] [dependencies] bytes = "1.5.0" +geo = "0.28.0" ndarray = "0.15.6" numpy = "0.20.0" object_store = { version = "0.9.0", features = ["http"] } diff --git a/src/io/geotiff.rs b/src/io/geotiff.rs index f266c5e..d843321 100644 --- a/src/io/geotiff.rs +++ b/src/io/geotiff.rs @@ -1,8 +1,10 @@ use std::io::{Read, Seek}; +use geo::AffineTransform; use ndarray::Array2; use tiff::decoder::{Decoder, DecodingResult, Limits}; -use tiff::{TiffFormatError, TiffResult}; +use tiff::tags::Tag; +use tiff::{TiffError, TiffFormatError, TiffResult}; /// Cloud-optimized GeoTIFF reader struct CogReader { @@ -18,6 +20,55 @@ impl CogReader { Ok(Self { decoder }) } + + /// Affine transformation for 2D matrix extracted from TIFF tag metadata, used to transform + /// image pixel (row, col) coordinates to and from geographic/projected (x, y) coordinates. + /// + /// ```text + /// | x' | | a b c | | x | + /// | y' | = | d e f | | y | + /// | 1 | | 0 0 1 | | 1 | + /// ``` + /// + /// where (`x'` and `y'`) are world coordinates, and (`x`, `y) are the pixel's image + /// coordinates. Letters a to f represent: + /// + /// - `a` - width of a pixel (x-resolution) + /// - `b` - row rotation (typically zero) + /// - `c` - x-coordinate of the *center* of the upper-left pixel (x-origin) + /// - `d` - column rotation (typically zero) + /// - `e` - height of a pixel (y-resolution, typically negative) + /// - `f` - y-coordinate of the *center* of the upper-left pixel (y-origin) + /// + /// References: + /// - + fn transform(&mut self) -> TiffResult> { + // Get x and y axis rotation (not yet implemented) + let (x_rotation, y_rotation): (f64, f64) = + match self.decoder.get_tag_f64_vec(Tag::ModelTransformationTag) { + Ok(_model_transformation) => unimplemented!("Non-zero rotation is not handled yet"), + Err(_) => (0.0, 0.0), + }; + + // Get pixel size in x and y direction + let pixel_scale: Vec = self.decoder.get_tag_f64_vec(Tag::ModelPixelScaleTag)?; + let [x_scale, y_scale, _z_scale] = pixel_scale[0..3] else { + return Err(TiffError::FormatError(TiffFormatError::InvalidTag)); + }; + + // Get x and y coordinates of upper left pixel + let tie_points: Vec = self.decoder.get_tag_f64_vec(Tag::ModelTiepointTag)?; + let [_i, _j, _k, x_origin, y_origin, _z_origin] = tie_points[0..6] else { + return Err(TiffError::FormatError(TiffFormatError::InvalidTag)); + }; + + // Create affine transformation matrix + let transform = AffineTransform::new( + x_scale, x_rotation, x_origin, y_rotation, -y_scale, y_origin, + ); + + Ok(transform) + } } /// Synchronously read a GeoTIFF file into an [`ndarray::Array`] @@ -42,12 +93,15 @@ pub fn read_geotiff(stream: R) -> TiffResult> { #[cfg(test)] mod tests { - use std::io::{Seek, SeekFrom}; + use std::io::{Cursor, Seek, SeekFrom}; + use geo::AffineTransform; + use object_store::parse_url; use tempfile::tempfile; use tiff::encoder::{colortype, TiffEncoder}; + use url::Url; - use crate::io::geotiff::read_geotiff; + use crate::io::geotiff::{read_geotiff, CogReader}; #[test] fn test_read_geotiff() { @@ -76,4 +130,24 @@ mod tests { assert_eq!(arr.ncols(), 20); // x-axis assert_eq!(arr.mean(), Some(14.0)); } + + #[tokio::test] + async fn test_cogreader_transform() { + let cog_url: &str = + "https://github.com/cogeotiff/rio-tiler/raw/6.4.0/tests/fixtures/cog_nodata_nan.tif"; + let tif_url = Url::parse(cog_url).unwrap(); + let (store, location) = parse_url(&tif_url).unwrap(); + + let result = store.get(&location).await.unwrap(); + let bytes = result.bytes().await.unwrap(); + let stream = Cursor::new(bytes); + + let mut reader = CogReader::new(stream).unwrap(); + let transform = reader.transform().unwrap(); + + assert_eq!( + transform, + AffineTransform::new(200.0, 0.0, 499980.0, 0.0, -200.0, 5300040.0) + ); + } }