diff --git a/Cargo.lock b/Cargo.lock index 8822144..0bf02e8 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,18 +8,109 @@ version = "2.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "512761e0bb2578dd7380c6baaa0f4ce03e84f95e960231d1dec8bf4d7d6e2627" +[[package]] +name = "aho-corasick" +version = "1.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916" +dependencies = [ + "memchr", +] + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + [[package]] name = "autocfg" version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" +[[package]] +name = "bindgen" +version = "0.66.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2b84e06fc203107bfbad243f4aba2af864eb7db3b1cf46ea0a023b0b433d2a7" +dependencies = [ + "bitflags", + "cexpr", + "clang-sys", + "lazy_static", + "lazycell", + "log", + "peeking_take_while", + "prettyplease", + "proc-macro2", + "quote", + "regex", + "rustc-hash", + "shlex", + "syn", + "which", +] + +[[package]] +name = "bitflags" +version = "2.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de" + +[[package]] +name = "bytemuck" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94bbb0ad554ad961ddc5da507a12a29b14e4ae5bda06b19f575a3e6079d2e2ae" + +[[package]] +name = "cc" +version = "1.1.28" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2e80e3b6a3ab07840e1cae9b0666a63970dc28e8ed5ffbcdacbfc760c281bfc1" +dependencies = [ + "shlex", +] + +[[package]] +name = "cexpr" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6fac387a98bb7c37292057cffc56d62ecb629900026402633ae9160df93a8766" +dependencies = [ + "nom", +] + [[package]] name = "cfg-if" version = "1.0.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" +[[package]] +name = "clang-sys" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b023947811758c97c59bf9d1c188fd619ad4718dcaa767947df1cadb14f39f4" +dependencies = [ + "glob", + "libc", + "libloading", +] + +[[package]] +name = "cmake" +version = "0.1.51" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb1e43aa7fd152b1f968787f7dbcdeb306d1867ff373c69955211876c053f91a" +dependencies = [ + "cc", +] + [[package]] name = "crc32fast" version = "1.4.2" @@ -29,12 +120,49 @@ dependencies = [ "cfg-if", ] +[[package]] +name = "delaunator" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ab46e386c7a38300a0d93b0f3e484bc2ee0aded66c47b14762ec9ab383934fa" +dependencies = [ + "robust", +] + +[[package]] +name = "either" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" + [[package]] name = "equivalent" version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5" +[[package]] +name = "errno" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "534c5cf6194dfab3db3242765c03bbe257cf92f22b38f6bc0c58d59108a820ba" +dependencies = [ + "libc", + "windows-sys 0.52.0", +] + +[[package]] +name = "filetime" +version = "0.2.25" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "35c0522e981e68cbfa8c3f978441a5f34b30b96e146b33cd3359176b50fe8586" +dependencies = [ + "cfg-if", + "libc", + "libredox", + "windows-sys 0.59.0", +] + [[package]] name = "flate2" version = "1.0.33" @@ -45,21 +173,70 @@ 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 = "geo-index" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "818dd7464a1edadbe7932a09b8e3672216597557100f23d1315f351e46c2c20e" +dependencies = [ + "bytemuck", + "float_next_after", + "num-traits", + "thiserror", + "tinyvec", +] + +[[package]] +name = "geo-types" +version = "0.7.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ff16065e5720f376fbced200a5ae0f47ace85fd70b7e54269790281353b6d61" +dependencies = [ + "approx", + "num-traits", + "serde", +] + [[package]] name = "geotiff" version = "0.0.2" dependencies = [ + "delaunator", + "geo-index", + "geo-types", "num-traits", "num_enum", + "proj", "tiff", ] +[[package]] +name = "glob" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" + [[package]] name = "hashbrown" version = "0.14.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" +[[package]] +name = "home" +version = "0.5.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3d1354bf6b7235cb4a0576c2619fd4ed18183f689b12b006a0ee7329eeff9a5" +dependencies = [ + "windows-sys 0.52.0", +] + [[package]] name = "indexmap" version = "2.5.0" @@ -76,12 +253,75 @@ version = "0.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f5d4a7da358eff58addd2877a45865158f0d78c911d43a5784ceb7bbf52833b0" +[[package]] +name = "lazy_static" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" + +[[package]] +name = "lazycell" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "830d08ce1d1d941e6b30645f1a0eb5643013d835ce3779a5fc208261dbe10f55" + +[[package]] +name = "libc" +version = "0.2.159" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "561d97a539a36e26a9a5fad1ea11a3039a67714694aaa379433e580854bc3dc5" + +[[package]] +name = "libloading" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4979f22fdb869068da03c9f7528f8297c6fd2606bc3a4affe42e6a823fdb8da4" +dependencies = [ + "cfg-if", + "windows-targets", +] + +[[package]] +name = "libm" +version = "0.2.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" + +[[package]] +name = "libredox" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c0ff37bd590ca25063e35af745c343cb7a0271906fb7b37e4813e8f79f00268d" +dependencies = [ + "bitflags", + "libc", + "redox_syscall", +] + +[[package]] +name = "linux-raw-sys" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78b3ae25bc7c8c38cec158d1f2757ee79e9b3740fbc7ccf0e59e4b08d793fa89" + +[[package]] +name = "log" +version = "0.4.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a7a70ba024b9dc04c27ea2f0c0548feb474ec5c54bba33a7f72f873a39d07b24" + [[package]] name = "memchr" version = "2.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" +[[package]] +name = "minimal-lexical" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68354c5c6bd36d73ff3feceb05efa59b6acb7626617f4962be322a825e61f79a" + [[package]] name = "miniz_oxide" version = "0.8.0" @@ -91,6 +331,16 @@ dependencies = [ "adler2", ] +[[package]] +name = "nom" +version = "7.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d273983c5a657a70a3e8f2a01329822f3b8c8172b73826411a55751e404a0a4a" +dependencies = [ + "memchr", + "minimal-lexical", +] + [[package]] name = "num-traits" version = "0.2.19" @@ -98,6 +348,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", + "libm", ] [[package]] @@ -121,6 +372,34 @@ dependencies = [ "syn", ] +[[package]] +name = "once_cell" +version = "1.20.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1261fe7e33c73b354eab43b1273a57c8f967d0391e80353e51f764ac02cf6775" + +[[package]] +name = "peeking_take_while" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19b17cddbe7ec3f8bc800887bab5e717348c95ea2ca0b1bf0837fb964dc67099" + +[[package]] +name = "pkg-config" +version = "0.3.31" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "953ec861398dccce10c670dfeaf3ec4911ca479e9c02154b3a215178c5f566f2" + +[[package]] +name = "prettyplease" +version = "0.2.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "479cf940fbbb3426c32c5d5176f62ad57549a0bb84773423ba8be9d089f5faba" +dependencies = [ + "proc-macro2", + "syn", +] + [[package]] name = "proc-macro-crate" version = "3.2.0" @@ -139,6 +418,32 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "proj" +version = "0.27.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ad1830ad8966eba22c76e78440458f07bd812bef5c3efdf335dec55cd1085ab" +dependencies = [ + "geo-types", + "libc", + "num-traits", + "proj-sys", + "thiserror", +] + +[[package]] +name = "proj-sys" +version = "0.23.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "601bf4fa1e17fde1a56d303f7bed5c65969cf1822c6baf5d6c2c12c593638fec" +dependencies = [ + "bindgen", + "cmake", + "flate2", + "pkg-config", + "tar", +] + [[package]] name = "quote" version = "1.0.37" @@ -148,6 +453,95 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "redox_syscall" +version = "0.5.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b6dfecf2c74bce2466cabf93f6664d6998a69eb21e39f4207930065b27b771f" +dependencies = [ + "bitflags", +] + +[[package]] +name = "regex" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "38200e5ee88914975b69f657f0801b6f6dccafd44fd9326302a4aaeecfacb1d8" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "368758f23274712b504848e9d5a6f010445cc8b87a7cdb4d7cbee666c1288da3" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" + +[[package]] +name = "robust" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5864e7ef1a6b7bcf1d6ca3f655e65e724ed3b52546a0d0a663c991522f552ea" + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustix" +version = "0.38.37" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8acb788b847c24f28525660c4d7758620a7210875711f79e7f663cc152726811" +dependencies = [ + "bitflags", + "errno", + "libc", + "linux-raw-sys", + "windows-sys 0.52.0", +] + +[[package]] +name = "serde" +version = "1.0.210" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8e3592472072e6e22e0a54d5904d9febf8508f65fb8552499a1abc7d1078c3a" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.210" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "243902eda00fad750862fc144cea25caca5e20d615af0a81bee94ca738f1df1f" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + [[package]] name = "syn" version = "2.0.77" @@ -159,6 +553,37 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "tar" +version = "0.4.42" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ff6c40d3aedb5e06b57c6f669ad17ab063dd1e63d977c6a88e7f4dfa4f04020" +dependencies = [ + "filetime", + "libc", + "xattr", +] + +[[package]] +name = "thiserror" +version = "1.0.64" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d50af8abc119fb8bb6dbabcfa89656f46f84aa0ac7688088608076ad2b459a84" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.64" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08904e7672f5eb876eaaf87e0ce17857500934f4981c4a0ab2b4aa98baac7fc3" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + [[package]] name = "tiff" version = "0.9.1" @@ -170,6 +595,21 @@ dependencies = [ "weezl", ] +[[package]] +name = "tinyvec" +version = "1.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "445e881f4f6d382d5f27c034e25eb92edd7c784ceab92a0937db7f2e9471b938" +dependencies = [ + "tinyvec_macros", +] + +[[package]] +name = "tinyvec_macros" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f3ccbac311fea05f86f61904b462b55fb3df8837a366dfc601a0161d0532f20" + [[package]] name = "toml_datetime" version = "0.6.8" @@ -199,6 +639,100 @@ version = "0.1.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "53a85b86a771b1c87058196170769dd264f66c0782acf1ae6cc51bfd64b39082" +[[package]] +name = "which" +version = "4.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "87ba24419a2078cd2b0f2ede2691b6c66d8e47836da3b6db8265ebad47afbfc7" +dependencies = [ + "either", + "home", + "once_cell", + "rustix", +] + +[[package]] +name = "windows-sys" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "282be5f36a8ce781fad8c8ae18fa3f9beff57ec1b52cb3de0789201425d9a33d" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-targets" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_gnullvm", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469" + +[[package]] +name = "windows_i686_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b" + +[[package]] +name = "windows_i686_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66" + +[[package]] +name = "windows_i686_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" + [[package]] name = "winnow" version = "0.6.18" @@ -207,3 +741,14 @@ checksum = "68a9bda4691f099d435ad181000724da8e5899daa10713c2d432552b9ccd3a6f" dependencies = [ "memchr", ] + +[[package]] +name = "xattr" +version = "1.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8da84f1a25939b27f6820d92aed108f83ff920fdf11a7b19366c27c4cda81d4f" +dependencies = [ + "libc", + "linux-raw-sys", + "rustix", +] diff --git a/Cargo.toml b/Cargo.toml index 70fa012..866af36 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,10 +3,22 @@ name = "geotiff" description = "A GeoTIFF library for Rust" version = "0.0.2" edition = "2021" -authors = ["Dominik Bucher "] +authors = [ + "Dominik Bucher ", + "Gunnar Schulze " +] repository = "https://github.com/georust/geotiff" [dependencies] +delaunator = { version = "1.0", optional = true } +geo-index = { version = "0.1", optional = true } +geo-types = { version = "0.7" } num_enum = "0.7" num-traits = "0.2" tiff = "0.9" + +[dev-dependencies] +proj = "0.27" + +[features] +tie-points = ["dep:delaunator", "dep:geo-index"] diff --git a/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif b/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif new file mode 100644 index 0000000..b18b0b0 Binary files /dev/null and b/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif differ diff --git a/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif b/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif new file mode 100644 index 0000000..2bb7fa3 Binary files /dev/null and b/resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif differ diff --git a/resources/austrian_capitals_model_tie_points_pixel_is_area.tif b/resources/austrian_capitals_model_tie_points_pixel_is_area.tif new file mode 100644 index 0000000..3875c2c Binary files /dev/null and b/resources/austrian_capitals_model_tie_points_pixel_is_area.tif differ diff --git a/resources/austrian_capitals_model_tie_points_pixel_is_point.tif b/resources/austrian_capitals_model_tie_points_pixel_is_point.tif new file mode 100644 index 0000000..b8b7dc2 Binary files /dev/null and b/resources/austrian_capitals_model_tie_points_pixel_is_point.tif differ diff --git a/resources/austrian_capitals_model_transformation_pixel_is_area.tif b/resources/austrian_capitals_model_transformation_pixel_is_area.tif new file mode 100644 index 0000000..9a976dc Binary files /dev/null and b/resources/austrian_capitals_model_transformation_pixel_is_area.tif differ diff --git a/resources/austrian_capitals_model_transformation_pixel_is_point.tif b/resources/austrian_capitals_model_transformation_pixel_is_point.tif new file mode 100644 index 0000000..79004ba Binary files /dev/null and b/resources/austrian_capitals_model_transformation_pixel_is_point.tif differ diff --git a/src/coordinate_transform.rs b/src/coordinate_transform.rs new file mode 100644 index 0000000..a6b0c20 --- /dev/null +++ b/src/coordinate_transform.rs @@ -0,0 +1,179 @@ +#[cfg(feature = "tie-points")] +use std::rc::Rc; + +#[cfg(feature = "tie-points")] +use geo_index::rtree::OwnedRTree; +use geo_types::Coord; +use tiff::{TiffError, TiffFormatError, TiffResult}; + +#[cfg(feature = "tie-points")] +use crate::coordinate_transform::tie_points::Face; + +mod affine_transform; +mod tie_point_and_pixel_scale; +#[cfg(feature = "tie-points")] +mod tie_points; + +const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag"; +const MODEL_PIXEL_SCALE_TAG: &str = "ModelPixelScaleTag"; +const MODEL_TRANSFORMATION_TAG: &str = "ModelTransformationTag"; + +/// Defines the transformation between raster space and model space. +/// +/// Ref: https://docs.ogc.org/is/19-008r4/19-008r4.html#_raster_to_model_coordinate_transformation_requirements +#[derive(Debug)] +pub enum CoordinateTransform { + AffineTransform { + transform: [f64; 6], + inverse_transform: [f64; 6], + }, + TiePointAndPixelScale { + raster_point: Coord, + model_point: Coord, + pixel_scale: Coord, + }, + #[cfg(feature = "tie-points")] + TiePoints { + raster_mesh: Rc>, + raster_index: OwnedRTree, + model_mesh: Rc>, + model_index: OwnedRTree, + }, +} + +impl CoordinateTransform { + pub(super) fn from_tag_data( + pixel_scale_data: Option>, + model_tie_points_data: Option>, + model_transformation_data: Option>, + ) -> TiffResult { + let pixel_scale = pixel_scale_data + .map(|data| { + <[f64; 3]>::try_from(data).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Number values in {MODEL_PIXEL_SCALE_TAG} must be equal to 3" + ))) + }) + }) + .transpose()?; + let tie_points = model_tie_points_data + .map(|data| { + let len = data.len(); + if len == 0 { + return Err(TiffError::FormatError(TiffFormatError::Format(format!( + "Number of values in {MODEL_TIE_POINT_TAG} must be greater than 0" + )))); + } + + if len % 6 != 0 { + return Err(TiffError::FormatError(TiffFormatError::Format(format!( + "Number of values in {MODEL_TIE_POINT_TAG} must be divisible by 6" + )))); + } + + Ok(data) + }) + .transpose()?; + let transformation_matrix = model_transformation_data + .map(|data| { + <[f64; 16]>::try_from(data).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Number of values in {MODEL_TRANSFORMATION_TAG} must be equal to 16" + ))) + }) + }) + .transpose()?; + + if let Some(transformation_matrix) = transformation_matrix { + if pixel_scale.is_some() { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_PIXEL_SCALE_TAG} must not be specified when {MODEL_TRANSFORMATION_TAG} is present"), + ))); + } + if tie_points.is_some() { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_TIE_POINT_TAG} must not be specified when {MODEL_TRANSFORMATION_TAG} is present"), + ))); + } + + Self::from_transformation_matrix(transformation_matrix) + } else { + let Some(tie_points) = tie_points else { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_TIE_POINT_TAG} must be present when {MODEL_TRANSFORMATION_TAG} is missing"), + ))); + }; + + if tie_points.len() == 6 { + let Some(pixel_scale) = pixel_scale else { + return Err(TiffError::FormatError(TiffFormatError::Format( + format!("{MODEL_PIXEL_SCALE_TAG} must be specified when {MODEL_TIE_POINT_TAG} contains 6 values"), + ))); + }; + + Self::from_tie_point_and_pixel_scale(&tie_points, &pixel_scale) + } else { + #[cfg(feature = "tie-points")] + { + Self::from_tie_points(&tie_points) + } + #[cfg(not(feature = "tie-points"))] + { + Err(TiffError::FormatError(TiffFormatError::Format( + "Transformation by tie points is not supported".into(), + ))) + } + } + } + } + + pub fn transform_to_model(&self, coord: &Coord) -> Coord { + match self { + CoordinateTransform::AffineTransform { transform, .. } => { + Self::transform_by_affine_transform(transform, coord) + } + CoordinateTransform::TiePointAndPixelScale { + raster_point, + model_point, + pixel_scale, + } => Self::transform_to_model_by_tie_point_and_pixel_scale( + raster_point, + model_point, + pixel_scale, + coord, + ), + #[cfg(feature = "tie-points")] + CoordinateTransform::TiePoints { + raster_index, + raster_mesh, + model_mesh, + .. + } => Self::transform_by_tie_points(raster_index, raster_mesh, model_mesh, coord), + } + } + + pub(super) fn transform_to_raster(&self, coord: &Coord) -> Coord { + match self { + CoordinateTransform::AffineTransform { + inverse_transform, .. + } => Self::transform_by_affine_transform(inverse_transform, coord), + CoordinateTransform::TiePointAndPixelScale { + raster_point, + model_point, + pixel_scale, + } => Self::transform_to_raster_by_tie_point_and_pixel_scale( + raster_point, + model_point, + pixel_scale, + coord, + ), + #[cfg(feature = "tie-points")] + CoordinateTransform::TiePoints { + model_index, + model_mesh, + raster_mesh, + .. + } => Self::transform_by_tie_points(model_index, model_mesh, raster_mesh, coord), + } + } +} diff --git a/src/coordinate_transform/affine_transform.rs b/src/coordinate_transform/affine_transform.rs new file mode 100644 index 0000000..c97c3ba --- /dev/null +++ b/src/coordinate_transform/affine_transform.rs @@ -0,0 +1,45 @@ +use geo_types::Coord; +use tiff::{TiffError, TiffFormatError, TiffResult}; + +use crate::coordinate_transform::CoordinateTransform; + +impl CoordinateTransform { + pub fn from_transformation_matrix(transformation_matrix: [f64; 16]) -> TiffResult { + let transform = [ + transformation_matrix[0], + transformation_matrix[1], + transformation_matrix[3], + transformation_matrix[4], + transformation_matrix[5], + transformation_matrix[7], + ]; + + let det = transform[0] * transform[4] - transform[1] * transform[3]; + if det.abs() < 0.000000000000001 { + return Err(TiffError::FormatError(TiffFormatError::Format( + "Provided transformation matrix is not invertible".into(), + ))); + } + + let inverse_transform = [ + transform[4] / det, + -transform[1] / det, + (transform[1] * transform[5] - transform[2] * transform[4]) / det, + -transform[3] / det, + transform[0] / det, + (-transform[0] * transform[5] + transform[2] * transform[3]) / det, + ]; + + Ok(CoordinateTransform::AffineTransform { + transform, + inverse_transform, + }) + } + + pub(super) fn transform_by_affine_transform(transform: &[f64; 6], coord: &Coord) -> Coord { + Coord { + x: coord.x * transform[0] + coord.y * transform[1] + transform[2], + y: coord.x * transform[3] + coord.y * transform[4] + transform[5], + } + } +} diff --git a/src/coordinate_transform/tie_point_and_pixel_scale.rs b/src/coordinate_transform/tie_point_and_pixel_scale.rs new file mode 100644 index 0000000..6075752 --- /dev/null +++ b/src/coordinate_transform/tie_point_and_pixel_scale.rs @@ -0,0 +1,50 @@ +use geo_types::Coord; +use tiff::TiffResult; + +use crate::coordinate_transform::CoordinateTransform; + +impl CoordinateTransform { + pub(super) fn from_tie_point_and_pixel_scale( + tie_points: &[f64], + pixel_scale: &[f64], + ) -> TiffResult { + Ok(CoordinateTransform::TiePointAndPixelScale { + raster_point: Coord { + x: tie_points[0], + y: tie_points[1], + }, + model_point: Coord { + x: tie_points[3], + y: tie_points[4], + }, + pixel_scale: Coord { + x: pixel_scale[0], + y: pixel_scale[1], + }, + }) + } + + pub(super) fn transform_to_model_by_tie_point_and_pixel_scale( + raster_point: &Coord, + model_point: &Coord, + pixel_scale: &Coord, + coord: &Coord, + ) -> Coord { + Coord { + x: (coord.x - raster_point.x) * pixel_scale.x + model_point.x, + y: (coord.y - raster_point.y) * -pixel_scale.y + model_point.y, + } + } + + pub(super) fn transform_to_raster_by_tie_point_and_pixel_scale( + raster_point: &Coord, + model_point: &Coord, + pixel_scale: &Coord, + coord: &Coord, + ) -> Coord { + Coord { + x: (coord.x - model_point.x) / pixel_scale.x + raster_point.x, + y: (coord.y - model_point.y) / -pixel_scale.y + raster_point.y, + } + } +} diff --git a/src/coordinate_transform/tie_points.rs b/src/coordinate_transform/tie_points.rs new file mode 100644 index 0000000..4498170 --- /dev/null +++ b/src/coordinate_transform/tie_points.rs @@ -0,0 +1,347 @@ +use std::array; +use std::rc::Rc; + +use delaunator::{Point, Triangulation}; +use geo_index::rtree::sort::STRSort; +use geo_index::rtree::{OwnedRTree, RTreeBuilder, RTreeIndex}; +use geo_types::Coord; +use tiff::TiffResult; + +use crate::coordinate_transform::CoordinateTransform; + +impl CoordinateTransform { + pub(super) fn from_tie_points(tie_points: &[f64]) -> TiffResult { + let capacity = tie_points.iter().len() / 6; + let mut raster_points = Vec::with_capacity(capacity); + let mut model_points = Vec::with_capacity(capacity); + + for chunk in tie_points.chunks(6) { + raster_points.push(Point { + x: chunk[0], + y: chunk[1], + }); + model_points.push(Point { + x: chunk[3], + y: chunk[4], + }); + } + + let triangulation = delaunator::triangulate(&raster_points); + let raster_mesh = Rc::new(Self::build_faces(raster_points, &triangulation)); + let model_mesh = Rc::new(Self::build_faces(model_points, &triangulation)); + let raster_index = Self::build_index(&raster_mesh); + let model_index = Self::build_index(&model_mesh); + + Ok(Self::TiePoints { + raster_mesh: raster_mesh.clone(), + raster_index, + model_mesh: model_mesh.clone(), + model_index, + }) + } + + fn build_faces(points: Vec, triangulation: &Triangulation) -> Vec { + let Triangulation { + triangles, hull, .. + } = triangulation; + + let len = hull.len(); + let mut angle_bisectors = vec![None; points.len()]; + for i in 0..len { + let pi = hull[i]; + let ci = hull[(i + 1) % len]; + let ni = hull[(i + 2) % len]; + + let prev = &points[pi]; + let curr = &points[ci]; + let next = &points[ni]; + + let prev_curr = Coord { + x: curr.x - prev.x, + y: curr.y - prev.y, + } + .normalize(); + let next_curr = Coord { + x: curr.x - next.x, + y: curr.y - next.y, + } + .normalize(); + let direction = Coord { + x: prev_curr.x + next_curr.x, + y: prev_curr.y + next_curr.y, + } + .normalize(); + + angle_bisectors[ci] = Some(direction); + } + triangles + .chunks(3) + .map(|chunk| { + let i1 = chunk[0]; + let i2 = chunk[1]; + let i3 = chunk[2]; + + let b12 = hull.as_slice().contains_sequence(&chunk[0..2]); + let b23 = hull.as_slice().contains_sequence(&chunk[1..3]); + let b31 = hull.as_slice().contains_sequence(&[i3, i1]); + + let c1 = Coord { + x: points[i1].x, + y: points[i1].y, + }; + let c2 = Coord { + x: points[i2].x, + y: points[i2].y, + }; + let c3 = Coord { + x: points[i3].x, + y: points[i3].y, + }; + + let boundary = if b12 { + if b23 { + if b31 { + // Open + None + } else { + // Closed at edge 3-1 + Some(Boundary::Open { + coords: vec![c3, c1], + from_direction: angle_bisectors[i3].unwrap(), + to_direction: angle_bisectors[i1].unwrap(), + }) + } + } else if b31 { + // Closed at edge 2-3 + Some(Boundary::Open { + coords: vec![c2, c3], + from_direction: angle_bisectors[i2].unwrap(), + to_direction: angle_bisectors[i3].unwrap(), + }) + } else { + // Closed at edges 2-3 and 3-1 + Some(Boundary::Open { + coords: vec![c2, c3, c1], + from_direction: angle_bisectors[i2].unwrap(), + to_direction: angle_bisectors[i1].unwrap(), + }) + } + } else if b23 { + if b31 { + // Closed at edge 1-2 + Some(Boundary::Open { + coords: vec![c1, c2], + from_direction: angle_bisectors[i1].unwrap(), + to_direction: angle_bisectors[i2].unwrap(), + }) + } else { + // Closed at edges 1-2 and 3-1 + Some(Boundary::Open { + coords: vec![c3, c1, c2], + from_direction: angle_bisectors[i3].unwrap(), + to_direction: angle_bisectors[i2].unwrap(), + }) + } + } else if b31 { + // Closed at edges 1-2 and 2-3 + Some(Boundary::Open { + coords: vec![c1, c2, c3], + from_direction: angle_bisectors[i1].unwrap(), + to_direction: angle_bisectors[i3].unwrap(), + }) + } else { + // Closed + Some(Boundary::Closed { + coords: vec![c1, c2, c3, c1], + }) + }; + + Face { + boundary, + support_points: array::from_fn(|i| { + let point = &points[chunk[i]]; + Coord { + x: point.x, + y: point.y, + } + }), + } + }) + .collect() + } + + fn build_index(mesh: &[Face]) -> OwnedRTree { + let mut builder = RTreeBuilder::new(mesh.len()); + for face in mesh.iter() { + let (min_x, min_y, max_x, max_y) = face.compute_envelope(); + builder.add(min_x, min_y, max_x, max_y); + } + builder.finish::() + } + + pub(super) fn transform_by_tie_points( + source_index: &OwnedRTree, + source_mesh: &Rc>, + target_mesh: &Rc>, + coord: &Coord, + ) -> Coord { + let index = source_index + .search(coord.x, coord.y, coord.x, coord.y) + .into_iter() + .find(|face_index| source_mesh[*face_index].contains(coord)) + .unwrap(); + let uv = source_mesh[index].locate(coord); + target_mesh[index].interpolate(uv) + } +} + +#[derive(Debug)] +pub struct Face { + boundary: Option, + support_points: [Coord; 3], +} + +impl Face { + fn contains(&self, coord: &Coord) -> bool { + let Some(boundary) = &self.boundary else { + return true; + }; + + let check = |c1: &Coord, c2: &Coord| -> bool { + ((c2.x - c1.x) * (coord.y - c1.y) - (c2.y - c1.y) * (coord.x - c1.x)).is_sign_positive() + }; + + match boundary { + Boundary::Open { + coords, + from_direction, + to_direction, + } => { + check(&(coords[0] + *from_direction), &coords[1]) + && check(&coords[1], &(coords[1] + *to_direction)) + && coords.windows(2).all(|w| check(&w[0], &w[1])) + } + Boundary::Closed { coords } => { + let len = self.support_points.len(); + (0..=len).all(|i| check(&coords[i], &coords[(i + 1) % len])) + } + } + } + + fn locate(&self, coord: &Coord) -> [f64; 2] { + let [a, b, c] = &self.support_points; + let d = c.x * (a.y - b.y) - b.x * (a.y - c.y) + a.x * (b.y - c.y); + [ + -(coord.x * (a.y - c.y) - c.x * (a.y - coord.y) + a.x * (c.y - coord.y)) / d, + (coord.x * (a.y - b.y) - b.x * (a.y - coord.y) + a.x * (b.y - coord.y)) / d, + ] + } + + fn interpolate(&self, params: [f64; 2]) -> Coord { + let [u, v] = params; + let [a, b, c] = &self.support_points; + Coord { + x: -u * a.x - v * a.x + a.x + u * b.x + v * c.x, + y: -u * a.y - v * a.y + a.y + u * b.y + v * c.y, + } + } + + fn compute_envelope(&self) -> (f64, f64, f64, f64) { + let Some(boundary) = &self.boundary else { + return (f64::MIN, f64::MIN, f64::MAX, f64::MAX); + }; + + let mut min_x = f64::MAX; + let mut min_y = f64::MAX; + let mut max_x = f64::MIN; + let mut max_y = f64::MIN; + + let mut update_envelope = |coords: &[Coord]| { + for c in coords { + min_x = min_x.min(c.x); + min_y = min_y.min(c.y); + max_x = max_x.max(c.x); + max_y = max_y.max(c.y); + } + }; + + match boundary { + Boundary::Open { + coords, + from_direction, + to_direction, + } => { + update_envelope(coords); + + for direction in [from_direction, to_direction] { + // Compute right-hand side normal vector + let nx = direction.y; + let ny = -direction.x; + + if nx.is_sign_positive() { + max_x = f64::MAX; + } else { + min_x = f64::MIN; + } + + if ny.is_sign_positive() { + max_y = f64::MAX; + } else { + min_y = f64::MIN; + } + } + } + Boundary::Closed { coords } => update_envelope(coords), + } + + (min_x, min_y, max_x, max_y) + } +} + +#[derive(Debug)] +enum Boundary { + Open { + coords: Vec, + from_direction: Coord, + to_direction: Coord, + }, + Closed { + coords: Vec, + }, +} + +pub(crate) trait CoordExt { + fn normalize(self) -> Self; +} + +impl CoordExt for Coord { + fn normalize(mut self) -> Self { + let len = (self.x.powi(2) + self.y.powi(2)).sqrt(); + self.x /= len; + self.y /= len; + self + } +} + +trait HullExt +where + T: PartialEq, +{ + fn contains_sequence(&self, seq: &[T]) -> bool; +} + +impl HullExt for &[T] +where + T: PartialEq, +{ + fn contains_sequence(&self, seq: &[T]) -> bool { + let len = self.len(); + for i in 0..len { + let (a, b) = seq.split_at(seq.len().min(len - i)); + if self[i..].starts_with(a) && self.starts_with(b) { + return true; + } + } + false + } +} diff --git a/src/decoder_ext.rs b/src/decoder_ext.rs index 5ea1a83..adfd560 100644 --- a/src/decoder_ext.rs +++ b/src/decoder_ext.rs @@ -4,13 +4,44 @@ use tiff::decoder::Decoder; use tiff::tags::Tag; use tiff::TiffResult; +use crate::coordinate_transform::CoordinateTransform; use crate::geo_key_directory::GeoKeyDirectory; pub(super) trait DecoderExt { + fn coordinate_transform(&mut self) -> TiffResult>; + fn geo_key_directory(&mut self) -> TiffResult; } impl DecoderExt for Decoder { + fn coordinate_transform(&mut self) -> TiffResult> { + let pixel_scale_data = self + .find_tag(Tag::ModelPixelScaleTag)? + .map(|value| value.into_f64_vec()) + .transpose()?; + let tie_points_data = self + .find_tag(Tag::ModelTiepointTag)? + .map(|value| value.into_f64_vec()) + .transpose()?; + let model_transformation_data = self + .find_tag(Tag::ModelTransformationTag)? + .map(|value| value.into_f64_vec()) + .transpose()?; + + if pixel_scale_data.is_none() + && tie_points_data.is_none() + && model_transformation_data.is_none() + { + return Ok(None); + } + + Ok(Some(CoordinateTransform::from_tag_data( + pixel_scale_data, + tie_points_data, + model_transformation_data, + )?)) + } + fn geo_key_directory(&mut self) -> TiffResult { let Some(directory_data) = self .find_tag(Tag::GeoKeyDirectoryTag)? diff --git a/src/geo_key_directory.rs b/src/geo_key_directory.rs index c74b7df..f0f890e 100644 --- a/src/geo_key_directory.rs +++ b/src/geo_key_directory.rs @@ -12,7 +12,7 @@ pub struct GeoKeyDirectory { pub key_revision: u16, pub minor_revision: u16, pub model_type: Option, - pub raster_type: Option, + pub raster_type: Option, pub citation: Option, pub geographic_type: Option, pub geog_citation: Option, @@ -99,8 +99,14 @@ impl GeoKeyDirectory { Self::get_short(key_tag, location_tag, *count, *value_or_offset)?.into() } GeoKeyDirectoryTag::RasterType => { + let raster_type = + Self::get_short(key_tag, location_tag, *count, *value_or_offset)?; directory.raster_type = - Self::get_short(key_tag, location_tag, *count, *value_or_offset)?.into() + Some(RasterType::try_from(raster_type).map_err(|_| { + TiffError::FormatError(TiffFormatError::Format(format!( + "Unknown raster type: {raster_type}" + ))) + })?) } GeoKeyDirectoryTag::Citation => { directory.citation = Self::get_string( @@ -657,3 +663,15 @@ enum GeoKeyDirectoryTag { VerticalDatum = 4098, VerticalUnits = 4099, } + +/// The raster type establishes the raster space used. +/// +/// Ref: https://docs.ogc.org/is/19-008r4/19-008r4.html#_requirements_class_gtrastertypegeokey +#[derive(Debug, Clone, Copy, PartialEq, TryFromPrimitive, IntoPrimitive)] +#[repr(u16)] +pub enum RasterType { + Undefined = 0, + RasterPixelIsArea = 1, + RasterPixelIsPoint = 2, + UserDefined = 32767, +} diff --git a/src/lib.rs b/src/lib.rs index 6ac7137..f69276d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,6 +2,7 @@ use std::any::type_name; use std::io::{Read, Seek}; +use geo_types::{Coord, Rect}; use num_traits::FromPrimitive; use tiff::decoder::{Decoder, DecodingResult}; use tiff::tags::Tag; @@ -9,9 +10,11 @@ use tiff::TiffResult; pub use crate::geo_key_directory::*; +use crate::coordinate_transform::*; use crate::decoder_ext::*; use crate::raster_data::*; +mod coordinate_transform; mod decoder_ext; mod geo_key_directory; mod raster_data; @@ -39,14 +42,17 @@ pub struct GeoTiff { pub raster_width: usize, pub raster_height: usize, pub num_samples: usize, + coordinate_transform: Option, raster_data: RasterData, } impl GeoTiff { + /// Reads a GeoTIFF from the given source. pub fn read(reader: R) -> TiffResult { let mut decoder = Decoder::new(reader)?; let geo_key_directory = decoder.geo_key_directory()?; + let coordinate_transform = decoder.coordinate_transform()?; let (raster_width, raster_height) = decoder .dimensions() @@ -74,27 +80,43 @@ impl GeoTiff { raster_width, raster_height, num_samples, + coordinate_transform, raster_data, }) } - pub fn get_value_at(&self, x: usize, y: usize, sample: usize) -> T { - let GeoTiff { - raster_width, - num_samples, - raster_data, - .. - } = self; + /// Returns the extent of the image in model space. + pub fn model_extent(&self) -> Rect { + let offset = self.raster_offset(); + let lower = Coord { + x: offset, + y: offset, + }; + let upper = Coord { + x: self.raster_width as f64 + offset, + y: self.raster_height as f64 + offset, + }; - if &sample >= num_samples { - panic!( - "sample out of bounds: the number of samples is {} but the sample is {}", - num_samples, sample + if let Some(coordinate_transform) = &self.coordinate_transform { + Rect::new( + coordinate_transform.transform_to_model(&lower), + coordinate_transform.transform_to_model(&upper), ) + } else { + Rect::new(lower, upper) } + } - let index = (y * raster_width + x) * num_samples + sample; - match raster_data { + /// Returns the value at the given location for the specified sample. + /// The coordinates are in model space. + pub fn get_value_at( + &self, + coord: &Coord, + sample: usize, + ) -> Option { + let index = self.compute_index(coord, sample)?; + + Some(match &self.raster_data { RasterData::U8(data) => unwrap_primitive_type!(T::from_u8(data[index]), u8, T), RasterData::U16(data) => unwrap_primitive_type!(T::from_u16(data[index]), u16, T), RasterData::U32(data) => unwrap_primitive_type!(T::from_u32(data[index]), u32, T), @@ -105,6 +127,50 @@ impl GeoTiff { RasterData::I16(data) => unwrap_primitive_type!(T::from_i16(data[index]), i16, T), RasterData::I32(data) => unwrap_primitive_type!(T::from_i32(data[index]), i32, T), RasterData::I64(data) => unwrap_primitive_type!(T::from_i64(data[index]), i64, T), + }) + } + + fn compute_index(&self, coord: &Coord, sample: usize) -> Option { + let GeoTiff { + raster_width, + raster_height, + num_samples, + coordinate_transform, + .. + } = self; + + if &sample >= num_samples { + panic!( + "sample out of bounds: the number of samples is {} but the sample is {}", + num_samples, sample + ) + } + + let mut coord = match coordinate_transform { + None => *coord, + Some(transform) => transform.transform_to_raster(coord), + }; + + // See https://docs.ogc.org/is/19-008r4/19-008r4.html#_raster_space for reference + let raster_offset = self.raster_offset(); + coord.x -= raster_offset; + coord.y -= raster_offset; + + if coord.x < 0.0 + || coord.x >= *raster_width as f64 + || coord.y < 0.0 + || coord.y >= *raster_height as f64 + { + return None; + } + + Some((coord.y as usize * raster_width + coord.x as usize) * num_samples + sample) + } + + fn raster_offset(&self) -> f64 { + match self.geo_key_directory.raster_type { + Some(RasterType::RasterPixelIsPoint) => -0.5, + _ => 0.0, } } } diff --git a/tests/common/mod.rs b/tests/common/mod.rs new file mode 100644 index 0000000..f6774ee --- /dev/null +++ b/tests/common/mod.rs @@ -0,0 +1,8 @@ +use std::fs::File; +use std::path::Path; + +use geotiff::GeoTiff; + +pub fn read_geotiff>(path: P) -> GeoTiff { + GeoTiff::read(File::open(path).expect("File I/O error")).expect("File I/O error") +} diff --git a/tests/integration.rs b/tests/integration.rs index 2a2181f..ba06fae 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -1,11 +1,8 @@ -use std::fs::File; -use std::path::Path; +use common::read_geotiff; +use geo_types::{Coord, Rect}; +use geotiff::{GeoKeyDirectory, RasterType}; -use geotiff::{GeoKeyDirectory, GeoTiff}; - -fn read_geotiff>(path: P) -> GeoTiff { - GeoTiff::read(File::open(path).expect("File I/O error")).expect("File I/O error") -} +mod common; #[test] fn test_load_marbles() { @@ -15,9 +12,28 @@ fn test_load_marbles() { assert_eq!(geotiff.raster_width, 1419); assert_eq!(geotiff.raster_height, 1001); assert_eq!(geotiff.num_samples, 3); - assert_eq!(geotiff.get_value_at::(761, 599, 0), 147); - assert_eq!(geotiff.get_value_at::(761, 599, 1), 128); - assert_eq!(geotiff.get_value_at::(761, 599, 2), 165); + assert_eq!( + geotiff.model_extent(), + Rect::new( + Coord { x: 0.0, y: 0.0 }, + Coord { + x: 1419.0, + y: 1001.0 + } + ) + ); + assert_eq!( + geotiff.get_value_at::(&Coord { x: 761.0, y: 599.0 }, 0), + Some(147) + ); + assert_eq!( + geotiff.get_value_at::(&Coord { x: 761.0, y: 599.0 }, 1), + Some(128) + ); + assert_eq!( + geotiff.get_value_at::(&Coord { x: 761.0, y: 599.0 }, 2), + Some(165) + ); } #[test] @@ -28,9 +44,49 @@ fn test_load_zh_dem_25() { assert_eq!(geotiff.raster_width, 399); assert_eq!(geotiff.raster_height, 366); assert_eq!(geotiff.num_samples, 1); - assert_eq!(geotiff.get_value_at::(0, 0, 0), 551); - assert_eq!(geotiff.get_value_at::(67, 45, 0), 530); - assert_eq!(geotiff.get_value_at::(325, 142, 0), 587); + assert_eq!( + geotiff.model_extent(), + Rect::new( + Coord { + x: 677562.5, + y: 243862.5 + }, + Coord { + x: 687537.5, + y: 253012.5 + } + ) + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: 677575.0, + y: 253000.0 + }, + 0 + ), + Some(551) + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: 679250.0, + y: 251875.0 + }, + 0 + ), + Some(530) + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: 685700.0, + y: 249450.0 + }, + 0 + ), + Some(587) + ); assert_eq!( geotiff.geo_key_directory, @@ -56,7 +112,7 @@ fn test_load_merc() { key_revision: 1, minor_revision: 2, model_type: Some(1), - raster_type: Some(1), + raster_type: Some(RasterType::RasterPixelIsArea), geog_geodetic_datum: Some(6267), geog_ellipsoid: Some(7008), projected_type: Some(32767), @@ -73,4 +129,18 @@ fn test_load_merc() { ..Default::default() } ); + + assert_eq!( + geotiff.model_extent(), + Rect::new( + Coord { + x: 1871032.9538880002, + y: 662408.6726400064 + }, + Coord { + x: 1901982.949391994, + y: 693358.6681440001 + } + ) + ); } diff --git a/tests/transform.rs b/tests/transform.rs new file mode 100644 index 0000000..2634933 --- /dev/null +++ b/tests/transform.rs @@ -0,0 +1,215 @@ +use std::path::Path; + +use common::read_geotiff; +use geo_types::{Coord, Rect}; +use geotiff::RasterType; +use proj::Proj; + +mod common; + +trait RectExt { + fn rounded(self) -> Self; +} + +impl RectExt for Rect { + fn rounded(mut self) -> Self { + let factor = 10f64.powi(8); + let mut min = self.min(); + let mut max = self.max(); + min.x = (min.x * factor).round() / factor; + min.y = (min.y * factor).round() / factor; + max.x = (max.x * factor).round() / factor; + max.y = (max.y * factor).round() / factor; + self.set_min(min); + self.set_max(max); + self + } +} + +const BREGENZ: Coord = Coord { + x: 9.74926, + y: 47.50315, +}; +const EISENSTADT: Coord = Coord { + x: 15.43301, + y: 47.06298, +}; +const GRAZ: Coord = Coord { + x: 15.43301, + y: 47.06298, +}; +const INNSBRUCK: Coord = Coord { + x: 11.39960, + y: 47.26239, +}; +const KLAGENFURT: Coord = Coord { + x: 14.31528, + y: 46.62366, +}; +const LINZ: Coord = Coord { + x: 14.30571, + y: 48.27532, +}; +const SALZBURG: Coord = Coord { + x: 13.05345, + y: 47.80763, +}; +const SANKT_POELTEN: Coord = Coord { + x: 15.62291, + y: 48.20440, +}; +const VIENNA: Coord = Coord { + x: 16.37499, + y: 48.22158, +}; + +const WHITE: u8 = 255; +const BLACK: u8 = 0; + +#[test] +fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_tie_point_and_pixel_scale_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_tie_point_and_pixel_scale_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +#[test] +fn test_transform_by_model_transformation_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_transformation_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[test] +fn test_transform_by_model_transformation_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_transformation_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +#[cfg(feature = "tie-points")] +#[test] +fn test_transform_by_model_tie_points_pixel_is_area() { + test_transform( + "resources/austrian_capitals_model_tie_points_pixel_is_area.tif", + RasterType::RasterPixelIsArea, + ) +} + +#[cfg(feature = "tie-points")] +#[test] +fn test_transform_by_model_tie_points_pixel_is_point() { + test_transform( + "resources/austrian_capitals_model_tie_points_pixel_is_point.tif", + RasterType::RasterPixelIsPoint, + ) +} + +fn test_transform>(path: P, raster_type: RasterType) { + let geotiff = read_geotiff(path); + + let proj = Proj::new_known_crs( + "EPSG:4326", + &format!("EPSG:{}", geotiff.geo_key_directory.projected_type.unwrap()), + None, + ) + .unwrap(); + + let mut bregenz = proj.project(BREGENZ, false).unwrap(); + let mut eisenstadt = proj.project(EISENSTADT, false).unwrap(); + let mut graz = proj.project(GRAZ, false).unwrap(); + let mut innsbruck = proj.project(INNSBRUCK, false).unwrap(); + let mut klagenfurt = proj.project(KLAGENFURT, false).unwrap(); + let mut linz = proj.project(LINZ, false).unwrap(); + let mut salzburg = proj.project(SALZBURG, false).unwrap(); + let mut sankt_poelten = proj.project(SANKT_POELTEN, false).unwrap(); + let mut vienna = proj.project(VIENNA, false).unwrap(); + let mut expected_model_extent = Rect::new( + Coord { + x: 4302000.0, + y: 2621000.0, + }, + Coord { + x: 4809000.0, + y: 2811000.0, + }, + ); + + // If raster type is RasterPixelIsPoint, shift coordinates accordingly + if raster_type == RasterType::RasterPixelIsPoint { + let mut min = expected_model_extent.min(); + let mut max = expected_model_extent.max(); + + let mut coords = [ + &mut bregenz, + &mut eisenstadt, + &mut graz, + &mut innsbruck, + &mut klagenfurt, + &mut linz, + &mut salzburg, + &mut sankt_poelten, + &mut vienna, + &mut min, + &mut max, + ]; + + for coord in coords.iter_mut() { + coord.x -= 500.0; + coord.y += 500.0; + } + + expected_model_extent.set_min(min); + expected_model_extent.set_max(max); + } + + assert_eq!(geotiff.model_extent().rounded(), expected_model_extent); + + // Non-capital location should return background color + assert_eq!( + geotiff.get_value_at::(&expected_model_extent.center(), 0), + Some(WHITE) + ); + + // A location outside the model extent should return None + + let min = expected_model_extent.min(); + let max = expected_model_extent.max(); + + assert_eq!( + geotiff.get_value_at::(&Coord { x: min.x, y: min.y }, 0), + None + ); + assert_eq!( + geotiff.get_value_at::( + &Coord { + x: max.x + 1.0, + y: max.y + 1.0, + }, + 0 + ), + None + ); + + // Each capital location should return Some(BLACK) + assert_eq!(geotiff.get_value_at::(&bregenz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&eisenstadt, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&graz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&innsbruck, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&klagenfurt, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&linz, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&salzburg, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&sankt_poelten, 0), Some(BLACK)); + assert_eq!(geotiff.get_value_at::(&vienna, 0), Some(BLACK)); +}