diff --git a/.github/workflows/maturin.yml b/.github/workflows/maturin.yml new file mode 100644 index 0000000..5ef1f58 --- /dev/null +++ b/.github/workflows/maturin.yml @@ -0,0 +1,120 @@ +# This file is autogenerated by maturin v1.3.2 +# To update, run +# +# maturin generate-ci github --manifest-path=japan-geoid-py/Cargo.toml +# +name: CI + +on: + push: + branches: + - main + - master + tags: + - '*' + pull_request: + workflow_dispatch: + +permissions: + contents: read + +jobs: + linux: + runs-on: ubuntu-latest + strategy: + matrix: + target: [x86_64, x86, aarch64, armv7, s390x, ppc64le] + steps: + - uses: actions/checkout@v3 + - uses: actions/setup-python@v4 + with: + python-version: '3.10' + - name: Build wheels + uses: PyO3/maturin-action@v1 + with: + target: ${{ matrix.target }} + args: --release --out dist --find-interpreter --manifest-path japan-geoid-py/Cargo.toml + sccache: 'true' + manylinux: auto + - name: Upload wheels + uses: actions/upload-artifact@v3 + with: + name: wheels + path: dist + + windows: + runs-on: windows-latest + strategy: + matrix: + target: [x64, x86] + steps: + - uses: actions/checkout@v3 + - uses: actions/setup-python@v4 + with: + python-version: '3.10' + architecture: ${{ matrix.target }} + - name: Build wheels + uses: PyO3/maturin-action@v1 + with: + target: ${{ matrix.target }} + args: --release --out dist --find-interpreter --manifest-path japan-geoid-py/Cargo.toml + sccache: 'true' + - name: Upload wheels + uses: actions/upload-artifact@v3 + with: + name: wheels + path: dist + + macos: + runs-on: macos-latest + strategy: + matrix: + target: [x86_64, aarch64] + steps: + - uses: actions/checkout@v3 + - uses: actions/setup-python@v4 + with: + python-version: '3.10' + - name: Build wheels + uses: PyO3/maturin-action@v1 + with: + target: ${{ matrix.target }} + args: --release --out dist --find-interpreter --manifest-path japan-geoid-py/Cargo.toml + sccache: 'true' + - name: Upload wheels + uses: actions/upload-artifact@v3 + with: + name: wheels + path: dist + + sdist: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - name: Build sdist + uses: PyO3/maturin-action@v1 + with: + command: sdist + args: --out dist --manifest-path japan-geoid-py/Cargo.toml + - name: Upload sdist + uses: actions/upload-artifact@v3 + with: + name: wheels + path: dist + + release: + name: Release + runs-on: ubuntu-latest + if: "startsWith(github.ref, 'refs/tags/')" + needs: [linux, windows, macos, sdist] + permissions: + id-token: write + steps: + - uses: actions/download-artifact@v3 + with: + name: wheels + - name: Publish to PyPI + uses: PyO3/maturin-action@v1 + with: + command: upload + args: --non-interactive --skip-existing * diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6775140 --- /dev/null +++ b/.gitignore @@ -0,0 +1,73 @@ +/target +gsigeo2011_ver2_2.* + +# Byte-compiled / optimized / DLL files +__pycache__/ +.pytest_cache/ +*.py[cod] + +# C extensions +*.so + +# Distribution / packaging +.Python +.venv/ +env/ +bin/ +build/ +develop-eggs/ +dist/ +eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +include/ +man/ +venv/ +*.egg-info/ +.installed.cfg +*.egg + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt +pip-selfcheck.json + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.cache +nosetests.xml +coverage.xml + +# Translations +*.mo + +# Mr Developer +.mr.developer.cfg +.project +.pydevproject + +# Rope +.ropeproject + +# Django stuff: +*.log +*.pot + +.DS_Store + +# Sphinx documentation +docs/_build/ + +# PyCharm +.idea/ + +# VSCode +.vscode/ + +# Pyenv +.python-version \ No newline at end of file diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..e2903a0 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,401 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "adler" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f26201604c87b1e01bd3d98f8d5d9a8fcbb815e8cedb41ffccbeb4bf593a35fe" + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "crc32fast" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b540bd8bc810d3885c6ea91e2018302f68baba2129ab3e88f32389ee9370880d" +dependencies = [ + "cfg-if", +] + +[[package]] +name = "flate2" +version = "1.0.27" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c6c98ee8095e9d1dcbf2fcc6d95acccb90d1c81db1e44725c6a984b1dbdfb010" +dependencies = [ + "crc32fast", + "miniz_oxide", +] + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "indoc" +version = "2.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e186cfbae8084e513daff4240b4797e342f988cecda4fb6c939150f96315fd8" + +[[package]] +name = "japan-geoid" +version = "0.1.0" +dependencies = [ + "flate2", +] + +[[package]] +name = "japan-geoid-py" +version = "0.0.2" +dependencies = [ + "japan-geoid", + "numpy", + "pyo3", +] + +[[package]] +name = "libc" +version = "0.2.150" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "89d92a4743f9a61002fae18374ed11e7973f530cb3a3255fb354818118b2203c" + +[[package]] +name = "lock_api" +version = "0.4.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c168f8615b12bc01f9c17e2eb0cc07dcae1940121185446edc3744920e8ef45" +dependencies = [ + "autocfg", + "scopeguard", +] + +[[package]] +name = "matrixmultiply" +version = "0.3.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "memoffset" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "miniz_oxide" +version = "0.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e7810e0be55b428ada41041c41f32c9f1a42817901b4ccf45fa3d4b6561e74c7" +dependencies = [ + "adler", +] + +[[package]] +name = "ndarray" +version = "0.15.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" +dependencies = [ + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "rawpointer", +] + +[[package]] +name = "num-complex" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ba157ca0885411de85d6ca030ba7e2a83a28636056c7c699b07c8b6f7383214" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "39e3200413f237f41ab11ad6d161bc7239c84dcb631773ccd7de3dfe4b5c267c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "numpy" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef41cbb417ea83b30525259e30ccef6af39b31c240bda578889494c5392d331" +dependencies = [ + "libc", + "ndarray", + "num-complex", + "num-integer", + "num-traits", + "pyo3", + "rustc-hash", +] + +[[package]] +name = "once_cell" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" + +[[package]] +name = "parking_lot" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3742b2c103b9f06bc9fff0a37ff4912935851bee6d36f3c02bcc755bcfec228f" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4c42a9226546d68acdd9c0a280d17ce19bfe27a46bf68784e4066115788d008e" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + +[[package]] +name = "proc-macro2" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "134c189feb4956b20f6f547d2cf727d4c0fe06722b20a0eec87ed445a97f92da" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "04e8453b658fe480c3e70c8ed4e3d3ec33eb74988bd186561b0cc66b85c3bc4b" +dependencies = [ + "cfg-if", + "indoc", + "libc", + "memoffset", + "parking_lot", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", + "unindent", +] + +[[package]] +name = "pyo3-build-config" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a96fe70b176a89cff78f2fa7b3c930081e163d5379b4dcdf993e3ae29ca662e5" +dependencies = [ + "once_cell", + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "214929900fd25e6604661ed9cf349727c8920d47deff196c4e28165a6ef2a96b" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dac53072f717aa1bfa4db832b39de8c875b7c7af4f4a6fe93cdbf9264cf8383b" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7774b5a8282bd4f25f803b1f0d945120be959a36c72e08e7cd031c792fdfd424" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "quote" +version = "1.0.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "redox_syscall" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4722d768eff46b75989dd134e5c353f0d6296e5aaa3132e776cbdb56be7731aa" +dependencies = [ + "bitflags", +] + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "smallvec" +version = "1.11.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4dccd0940a2dcdf68d092b8cbab7dc0ad8fa938bf95787e1b916b0e3d0e8e970" + +[[package]] +name = "syn" +version = "2.0.39" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23e78b90f2fcf45d3e842032ce32e3f2d1545ba6636271dcbf24fa306d87be7a" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.12.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "14c39fd04924ca3a864207c66fc2cd7d22d7c016007f9ce846cbb9326331930a" + +[[package]] +name = "unicode-ident" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" + +[[package]] +name = "unindent" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" + +[[package]] +name = "windows-targets" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a2fa6e2155d7247be68c096456083145c183cbbbc2764150dda45a87197940c" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b38e32f0abccf9987a4e3079dfb67dcd799fb61361e53e2882c3cbaf0d905d8" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc35310971f3b2dbbf3f0690a219f40e2d9afcf64f9ab7cc1be722937c26b4bc" + +[[package]] +name = "windows_i686_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a75915e7def60c94dcef72200b9a8e58e5091744960da64ec734a6c6e9b3743e" + +[[package]] +name = "windows_i686_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f55c233f70c4b27f66c523580f78f1004e8b5a8b659e05a4eb49d4166cca406" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53d40abd2583d23e4718fddf1ebec84dbff8381c07cae67ff7768bbf19c6718e" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b7b52767868a23d5bab768e390dc5f5c55825b6d30b86c844ff2dc7414044cc" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed94fce61571a4006852b7389a063ab983c02eb1bb37b47f8272ce92d06d9538" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..5d4b63a --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,16 @@ +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[package] +name = "japan-geoid" +version = "0.1.0" +edition = "2021" + +[dev-dependencies] +flate2 = "1.0.27" + +[workspace] +resolver = "2" +members = ["./japan-geoid-py"] + +[profile.dev.package."*"] +opt-level = 1 \ No newline at end of file diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..1205cac --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,19 @@ +Copyright (c) 2023 MIERUNE Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..04dbfef --- /dev/null +++ b/README.md @@ -0,0 +1,7 @@ +# japan-geoid + +Rust and Python library to calculate geoid heights in Japan using [GSI's geoid model](https://fgd.gsi.go.jp/download/geoid.php). + +[国土地理院のジオイドモデル](https://fgd.gsi.go.jp/download/geoid.php)を用いて日本のジオイド高を計算する Rust および Python 用ライブラリ。 + +Work in progress. diff --git a/examples/convert_asc_to_bin.rs b/examples/convert_asc_to_bin.rs new file mode 100644 index 0000000..dc8fd18 --- /dev/null +++ b/examples/convert_asc_to_bin.rs @@ -0,0 +1,35 @@ +use flate2::{read::GzDecoder, write::GzEncoder}; +use std::{ + fs::File, + io::{BufReader, BufWriter}, +}; + +use japan_geoid::{Geoid, MemoryGrid}; + +fn main() -> std::io::Result<()> { + let (lng, lat) = (138.2839817085188, 37.12378643088312); + + // Load from the original ascii format. + let file = File::open("./gsigeo2011_ver2_2.asc")?; + let mut reader = BufReader::new(file); + let geoid = MemoryGrid::from_ascii_reader(&mut reader)?; + + let z = geoid.get_height(lng, lat); + println!("Input: (lng: {}, lat: {}) -> Geoid height: {}", lng, lat, z); + + // Dump as the binary format. + let file = File::create("./gsigeo2011_ver2_2.bin.gz")?; + let mut writer = GzEncoder::new(BufWriter::new(file), flate2::Compression::fast()); + geoid.to_binary_writer(&mut writer)?; + writer.finish()?; + + // Load from the binary. + let mut file = File::open("./gsigeo2011_ver2_2.bin.gz")?; + let mut reader = GzDecoder::new(&mut file); + let geoid = MemoryGrid::from_binary_reader(&mut reader)?; + + let z = geoid.get_height(lng, lat); + println!("Input: (lng: {}, lat: {}) -> Geoid height: {}", lng, lat, z); + + Ok(()) +} diff --git a/japan-geoid-py/.gitignore b/japan-geoid-py/.gitignore new file mode 100644 index 0000000..af3ca5e --- /dev/null +++ b/japan-geoid-py/.gitignore @@ -0,0 +1,72 @@ +/target + +# Byte-compiled / optimized / DLL files +__pycache__/ +.pytest_cache/ +*.py[cod] + +# C extensions +*.so + +# Distribution / packaging +.Python +.venv/ +env/ +bin/ +build/ +develop-eggs/ +dist/ +eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +include/ +man/ +venv/ +*.egg-info/ +.installed.cfg +*.egg + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt +pip-selfcheck.json + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.cache +nosetests.xml +coverage.xml + +# Translations +*.mo + +# Mr Developer +.mr.developer.cfg +.project +.pydevproject + +# Rope +.ropeproject + +# Django stuff: +*.log +*.pot + +.DS_Store + +# Sphinx documentation +docs/_build/ + +# PyCharm +.idea/ + +# VSCode +.vscode/ + +# Pyenv +.python-version \ No newline at end of file diff --git a/japan-geoid-py/Cargo.lock b/japan-geoid-py/Cargo.lock new file mode 100644 index 0000000..5ba934b --- /dev/null +++ b/japan-geoid-py/Cargo.lock @@ -0,0 +1,359 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "indoc" +version = "2.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e186cfbae8084e513daff4240b4797e342f988cecda4fb6c939150f96315fd8" + +[[package]] +name = "japan-geoid-py" +version = "0.1.0" +dependencies = [ + "numpy", + "pyo3", +] + +[[package]] +name = "libc" +version = "0.2.150" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "89d92a4743f9a61002fae18374ed11e7973f530cb3a3255fb354818118b2203c" + +[[package]] +name = "lock_api" +version = "0.4.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c168f8615b12bc01f9c17e2eb0cc07dcae1940121185446edc3744920e8ef45" +dependencies = [ + "autocfg", + "scopeguard", +] + +[[package]] +name = "matrixmultiply" +version = "0.3.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "memoffset" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "ndarray" +version = "0.15.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" +dependencies = [ + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "rawpointer", +] + +[[package]] +name = "num-complex" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ba157ca0885411de85d6ca030ba7e2a83a28636056c7c699b07c8b6f7383214" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "39e3200413f237f41ab11ad6d161bc7239c84dcb631773ccd7de3dfe4b5c267c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "numpy" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef41cbb417ea83b30525259e30ccef6af39b31c240bda578889494c5392d331" +dependencies = [ + "libc", + "ndarray", + "num-complex", + "num-integer", + "num-traits", + "pyo3", + "rustc-hash", +] + +[[package]] +name = "once_cell" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" + +[[package]] +name = "parking_lot" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3742b2c103b9f06bc9fff0a37ff4912935851bee6d36f3c02bcc755bcfec228f" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4c42a9226546d68acdd9c0a280d17ce19bfe27a46bf68784e4066115788d008e" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + +[[package]] +name = "proc-macro2" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "134c189feb4956b20f6f547d2cf727d4c0fe06722b20a0eec87ed445a97f92da" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "04e8453b658fe480c3e70c8ed4e3d3ec33eb74988bd186561b0cc66b85c3bc4b" +dependencies = [ + "cfg-if", + "indoc", + "libc", + "memoffset", + "parking_lot", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", + "unindent", +] + +[[package]] +name = "pyo3-build-config" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a96fe70b176a89cff78f2fa7b3c930081e163d5379b4dcdf993e3ae29ca662e5" +dependencies = [ + "once_cell", + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "214929900fd25e6604661ed9cf349727c8920d47deff196c4e28165a6ef2a96b" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dac53072f717aa1bfa4db832b39de8c875b7c7af4f4a6fe93cdbf9264cf8383b" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7774b5a8282bd4f25f803b1f0d945120be959a36c72e08e7cd031c792fdfd424" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "quote" +version = "1.0.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "redox_syscall" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4722d768eff46b75989dd134e5c353f0d6296e5aaa3132e776cbdb56be7731aa" +dependencies = [ + "bitflags", +] + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "smallvec" +version = "1.11.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4dccd0940a2dcdf68d092b8cbab7dc0ad8fa938bf95787e1b916b0e3d0e8e970" + +[[package]] +name = "syn" +version = "2.0.39" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23e78b90f2fcf45d3e842032ce32e3f2d1545ba6636271dcbf24fa306d87be7a" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.12.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "14c39fd04924ca3a864207c66fc2cd7d22d7c016007f9ce846cbb9326331930a" + +[[package]] +name = "unicode-ident" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" + +[[package]] +name = "unindent" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" + +[[package]] +name = "windows-targets" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a2fa6e2155d7247be68c096456083145c183cbbbc2764150dda45a87197940c" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b38e32f0abccf9987a4e3079dfb67dcd799fb61361e53e2882c3cbaf0d905d8" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc35310971f3b2dbbf3f0690a219f40e2d9afcf64f9ab7cc1be722937c26b4bc" + +[[package]] +name = "windows_i686_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a75915e7def60c94dcef72200b9a8e58e5091744960da64ec734a6c6e9b3743e" + +[[package]] +name = "windows_i686_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f55c233f70c4b27f66c523580f78f1004e8b5a8b659e05a4eb49d4166cca406" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53d40abd2583d23e4718fddf1ebec84dbff8381c07cae67ff7768bbf19c6718e" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b7b52767868a23d5bab768e390dc5f5c55825b6d30b86c844ff2dc7414044cc" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed94fce61571a4006852b7389a063ab983c02eb1bb37b47f8272ce92d06d9538" diff --git a/japan-geoid-py/Cargo.toml b/japan-geoid-py/Cargo.toml new file mode 100644 index 0000000..be7ed07 --- /dev/null +++ b/japan-geoid-py/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "japan-geoid-py" +version = "0.0.2" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html +[lib] +name = "japan_geoid" +crate-type = ["cdylib"] + +[dependencies] +pyo3 = "0.20.0" +numpy = "0.20.0" +japan-geoid = { path = "../" } diff --git a/japan-geoid-py/LICENSE.txt b/japan-geoid-py/LICENSE.txt new file mode 100644 index 0000000..1205cac --- /dev/null +++ b/japan-geoid-py/LICENSE.txt @@ -0,0 +1,19 @@ +Copyright (c) 2023 MIERUNE Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/japan-geoid-py/README.md b/japan-geoid-py/README.md new file mode 100644 index 0000000..94826b3 --- /dev/null +++ b/japan-geoid-py/README.md @@ -0,0 +1,11 @@ +# japan-geoid-py + +Python library to calculate geoid heights in Japan using [GSI's geoid model](https://fgd.gsi.go.jp/download/geoid.php). + +[国土地理院のジオイドモデル](https://fgd.gsi.go.jp/download/geoid.php)を用いて日本のジオイド高を計算する Python 用ライブラリ。 + +``` +pip install japan-geoid +``` + +Work in progress. diff --git a/japan-geoid-py/examples/convert_asc_to_bin.py b/japan-geoid-py/examples/convert_asc_to_bin.py new file mode 100644 index 0000000..5e839d8 --- /dev/null +++ b/japan-geoid-py/examples/convert_asc_to_bin.py @@ -0,0 +1,23 @@ +""" +$ python3 convert_to_bin.py ./gsigeo2011_ver2_2.asc +""" + +import gzip +import sys +from pathlib import Path + +from japan_geoid import GsiGeoid + +if __name__ == "__main__": + path = Path(sys.argv[1]) + + with open(path, "r") as f: + geoid = GsiGeoid.from_ascii(f.read()) + + with gzip.open(path.with_suffix(".bin.gz"), "wb") as dest_f: + dest_f.write(geoid.to_binary()) + + (lng, lat) = (138.2839817085188, 37.12378643088312) + + height = geoid.get_height(lng, lat) + print(f"{lng=} {lat=} {height=}") diff --git a/japan-geoid-py/examples/geog2geoc.py b/japan-geoid-py/examples/geog2geoc.py new file mode 100644 index 0000000..bfe7b0f --- /dev/null +++ b/japan-geoid-py/examples/geog2geoc.py @@ -0,0 +1,66 @@ +""" +JGD2011 3D (EPSG:6697) -> WGS 84 Geograhic 3D (EPSG:4979) -> WGS 84 Geocentric 3D (EPSG:4978) + +$ echo 34.290 135.630 0 | python3 geog2geoc.py +""" + +import gzip +import math +import sys + +import pyproj +from japan_geoid import GsiGeoid + +# WGS 84 楕円体のパラメータ +A = 6378137.0 # 長半径 +INV_F = 298.257223563 # 偏平率の逆数 + +# 楕円体パラメータから導ける値 +_F = 1 / INV_F # 偏平率 +# _B = A * (1 - _F) # 短半径 +_E_SQ = _F * (2 - _F) # 離心率の2乗 +# _E = _E_SQ**0.5 # 離心率 + + +def geog2geoc(lat, lng, ellipsoidal_heght): + (lat_deg, lng_deg, height) = (lat, lng, ellipsoidal_heght) + (lat_rad, lng_rad) = (math.radians(lat_deg), math.radians(lng_deg)) + + tan_psi = (1 - _E_SQ) * math.tan(lat_rad) + z = A * (1 / (1 / (tan_psi * tan_psi) + 1 / ((1 - _F) ** 2))) ** 0.5 + r = A * (1 / (1 + (tan_psi * tan_psi) / ((1 - _F) ** 2))) ** 0.5 + + x = r * math.cos(lng_rad) + y = r * math.sin(lng_rad) + dhz = math.sin(lat_rad) + dhx = math.cos(lat_rad) * math.cos(lng_rad) + dhy = math.cos(lat_rad) * math.sin(lng_rad) + return (x + dhx * height, y + dhy * height, z + dhz * height) + + +if __name__ == "__main__": + with gzip.open("../gsigeo2011_ver2_2.bin.gz", "rb") as f: + geoid = GsiGeoid.from_binary(f.read()) + + proj_geog_tr = pyproj.Transformer.from_crs("EPSG:6697", "EPSG:4978") + proj_geoc_tr = pyproj.Transformer.from_crs("EPSG:6697", "EPSG:4979") + + for line in sys.stdin: + (lat, lng, elevation) = [float(v) for v in line.strip().split()] + + # PROJ による変換 (JGD2011 -> WGS84 Geographic 3D) + proj_lat, proj_lng, proj_height = proj_geoc_tr.transform(lat, lng, elevation) + # PROJ による変換 (JGD2011 -> WGS84 Geocentric 3D) + proj_x, proj_y, proj_z = proj_geog_tr.transform(lat, lng, elevation) + + # 自前で変換 + geoid_height = geoid.get_height(lng, lat) + height = geoid_height + elevation # ellipsoidal height + x, y, z = geog2geoc(lat, lng, height) + + print(f"{geoid_height=}") + print(f"{height=} {proj_height=} diff[m]={height - proj_height}") + print(f"OURS: {x=} {y=} {z=}") + print(f"PROJ: {proj_x=} {proj_y=} {proj_z=}") + geoc_diff = ((proj_x - x) ** 2 + (proj_y - y) ** 2 + (proj_z - z) ** 2) ** 0.5 + print(f"distance[m]={geoc_diff}") diff --git a/japan-geoid-py/examples/pure_python_gsi_geoid.py b/japan-geoid-py/examples/pure_python_gsi_geoid.py new file mode 100644 index 0000000..8d26786 --- /dev/null +++ b/japan-geoid-py/examples/pure_python_gsi_geoid.py @@ -0,0 +1,128 @@ +"""Pure-Python implementation""" + +import gzip +from array import array +from math import nan +from struct import pack, unpack + + +class GsiGeoid: + def __init__(self, f): + self._load_bin(f) + + @staticmethod + def convert_asc_to_bin(src_f, dest_f): + # Header + h0, h1, h2, h3, h4, h5, h6, h7 = src_f.readline().strip().split() + _lat_min = float(h0) + _lng_min = float(h1) + assert h2 == "0.016667" + assert h3 == "0.025000" + _lat_denom = 60 + _lng_denom = 40 + _n_lat = int(h4) + _n_lng = int(h5) + assert _lat_min == 20 + assert _lng_min == 120 + assert _n_lat == 1801 + assert _n_lng == 1201 + ikind = int(h6) + assert ikind == 1 + version = h7.strip() + assert version == "ver2.2" + + dest_f.write( + pack( + " float: + """Bilinear interpolation""" + + if x == 0: + if y == 0: + return v00 + else: + return v00 * (1 - y) + v10 * (1 - x) * y + elif y == 0: + return v00 * (1 - x) + v01 * x + else: + return ( + v00 * (1 - x) * (1 - y) + + v01 * x * (1 - y) + + v10 * (1 - x) * y + + v11 * x * y + ) + + def get_height(self, lng: float, lat: float) -> float: + x, x_residual = divmod((lng - self._lng_min) * self._lng_denom, 1) + y, y_residual = divmod((lat - self._lat_min) * self._lat_denom, 1) + x = int(x) + y = int(y) + if x < 0 or self._n_lng <= x or y < 0 or self._n_lat <= y: + return nan + + n_lng = self._n_lng + n_lat = self._n_lat + return self._bilinear( + x_residual, + y_residual, + self._points[n_lng * y + x], + self._points[n_lng * y + (x + 1)] if x < n_lng - 1 else nan, + self._points[n_lng * (y + 1) + x] if y < n_lat - 1 else nan, + ( + self._points[n_lng * (y + 1) + (x + 1)] + if y < n_lat - 1 and x < n_lng - 1 + else nan + ), + ) + + +if __name__ == "__main__": + # 初回のみ必要（元データをバイナリに変換） + with open("./gsigeo2011_ver2_2.asc") as src_f: + with gzip.open("./gsigeo2011_ver2_2.bin.gz", "wb") as dest_f: + GsiGeoid.convert_asc_to_bin(src_f, dest_f) + + # 2回目以降はこれだけでいい + with gzip.open("./gsigeo2011_ver2_2.bin.gz", "rb") as f: + geoid = GsiGeoid(f) + + z = geoid.get_height(138.2839817085188, 37.12378643088312) + print(z) + z = geoid.get_height(141.36199967724426, 43.06539278249951) + print(z) diff --git a/japan-geoid-py/japan_geoid.pyi b/japan-geoid-py/japan_geoid.pyi new file mode 100644 index 0000000..c96adb9 --- /dev/null +++ b/japan-geoid-py/japan_geoid.pyi @@ -0,0 +1,24 @@ +from typing import Self + +from numpy import ndarray + +class GsiGeoid: + """GSIGEO2011 geoid model for Japan.""" + + @classmethod + def from_ascii(cls, content: str) -> Self: + """Load the geoid model from the original ascii format.""" + return japan_geoid.GsiGeoid.from_ascii(content) + + @classmethod + def from_binary(cls, content: bytes) -> Self: + """Load the geoid model from the efficient binary format.""" + + def to_binary(self) -> bytes: + """Serialize the geoid model in the efficient binary format.""" + + def get_height(self, lng: float, lat: float) -> float: + """Get the geoid height at a specified point.""" + + def get_heights(self, lng: ndarray, lat: ndarray) -> float: + """Get the geoid height at each specified point.""" diff --git a/japan-geoid-py/pyproject.toml b/japan-geoid-py/pyproject.toml new file mode 100644 index 0000000..a178e94 --- /dev/null +++ b/japan-geoid-py/pyproject.toml @@ -0,0 +1,25 @@ +[build-system] +requires = ["maturin>=1.2,<2.0"] +build-backend = "maturin" + +[project] +name = "japan-geoid" +description = "Calculate geoid heights in Japan using GSI's geoid model." +readme = "README.md" +authors = [ + { name = "MIERUNE Inc.", email = "info@mierune.co.jp" }, +] +license = { file = "LICENSE.txt" } +requires-python = ">=3.8" +classifiers = [ + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: Implementation :: CPython", + "Programming Language :: Python :: Implementation :: PyPy", + "Programming Language :: Rust", +] + +[project.urls] +repository = "https://github.com/MIERUNE/japan-geoid/japan-geoid-py" + +[tool.maturin] +features = ["pyo3/extension-module"] diff --git a/japan-geoid-py/src/lib.rs b/japan-geoid-py/src/lib.rs new file mode 100644 index 0000000..bbbc9f2 --- /dev/null +++ b/japan-geoid-py/src/lib.rs @@ -0,0 +1,70 @@ +use ndarray::{Array, Zip}; +use numpy::*; +use pyo3::prelude::*; +use pyo3::types::PyType; +use std::borrow::Cow; + +use ::japan_geoid::*; + +#[pyclass] +struct GsiGeoid { + geoid: MemoryGrid<'static>, +} + +#[pymethods] +impl GsiGeoid { + /// Load the geoid model from the original ascii format. + #[classmethod] + fn from_ascii(_cls: &PyType, content: &str) -> PyResult { + let mut reader = std::io::Cursor::new(content); + let geoid = MemoryGrid::from_ascii_reader(&mut reader)?; + Ok(GsiGeoid { geoid }) + } + + /// Load the geoid model from the efficient binary format. + #[classmethod] + fn from_binary(_cls: &PyType, content: &[u8]) -> PyResult { + let mut reader = std::io::Cursor::new(content); + let geoid = MemoryGrid::from_binary_reader(&mut reader)?; + Ok(GsiGeoid { geoid }) + } + + /// Serialize the geoid model in the efficient binary format. + fn to_binary(&self) -> PyResult> { + let mut buf = Vec::new(); + self.geoid.to_binary_writer(&mut buf)?; + Ok(buf.into()) + } + + /// Get the geoid height at a specified point. + fn get_height(&self, lng: f64, lat: f64) -> f64 { + self.geoid.get_height(lng, lat) + } + + /// Get the geoid height at each specified point. + fn get_heights<'py>( + &self, + py: Python<'py>, + lng: PyReadonlyArrayDyn<'py, f64>, + lat: PyReadonlyArrayDyn<'py, f64>, + ) -> PyResult<&'py PyArrayDyn> { + if lng.shape() != lat.shape() { + return Err(PyErr::new::( + "lng and lat must have the same shape", + )); + } + let mut c = Array::zeros(lng.shape()); + Zip::from(&mut c) + .and(lng.as_array()) + .and(lat.as_array()) + .for_each(|c, &a, &b| *c = self.geoid.get_height(a, b)); + Ok(c.into_pyarray(py)) + } +} + +/// A Python module implemented in Rust. +#[pymodule] +fn japan_geoid(_py: Python, m: &PyModule) -> PyResult<()> { + m.add_class::()?; + Ok(()) +} diff --git a/japan-geoid-py/tests/__init__.py b/japan-geoid-py/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/japan-geoid-py/tests/test_simple.py b/japan-geoid-py/tests/test_simple.py new file mode 100644 index 0000000..73fc60c --- /dev/null +++ b/japan-geoid-py/tests/test_simple.py @@ -0,0 +1,28 @@ +import gzip + +import japan_geoid +import numpy as np +from pytest import approx + + +def test_test(): + with gzip.open("../../gsigeo2011_ver2_2.bin.gz", "rb") as f: + geoid = japan_geoid.GsiGeoid(f.read()) + + assert geoid.get_height(138.2839817085188, 37.12378643088312) == approx( + 39.47387210863509, 1e-12 + ) + assert geoid.get_height(141.36199967724426, 43.06539278249951) == approx( + 31.900711649958033, 1e-12 + ) + + print( + geoid.get_heights( + np.array([138.2839817085188, 141.36199967724426]), + np.array([37.12378643088312, 43.06539278249951]), + ) + ) + + +if __name__ == "__main__": + test_test() diff --git a/src/gsi.rs b/src/gsi.rs new file mode 100644 index 0000000..1b1a356 --- /dev/null +++ b/src/gsi.rs @@ -0,0 +1,227 @@ +use std::borrow::Cow; +use std::io::{self, BufRead, Read, Write}; + +#[derive(Debug)] +pub struct GsiGridInfo { + /// Number of grid points along X-axis + x_num: u32, + /// Number of grid points along Y-axis + y_num: u32, + /// Denominator of grid interval along X-axis + x_denom: u32, + /// Denominator of grid interval along Y-axis + y_denom: u32, + /// Minimum value of X-axis + x_min: f32, + /// Minimum value of Y-axis + y_min: f32, + /// ikind (not used) + ikind: u16, + /// Version + version: String, +} + +#[derive(Debug)] +pub struct MemoryGrid<'a> { + pub grid_info: GsiGridInfo, + points: Cow<'a, [i32]>, +} + +/// Bilinear interpolation +fn bilinear(x: f64, y: f64, v00: f64, v01: f64, v10: f64, v11: f64) -> f64 { + if x == 0.0 && y == 0.0 { + v00 + } else if x == 0.0 { + v00 * (1.0 - y) + v10 * y + } else if y == 0.0 { + v00 * (1.0 - x) + v01 * x + } else { + v00 * (1.0 - x) * (1.0 - y) + v01 * x * (1.0 - y) + v10 * (1.0 - x) * y + v11 * x * y + } +} + +pub trait Geoid { + fn get_height(&self, lng: f64, lat: f64) -> f64; +} + +pub trait Grid { + fn grid_info(&self) -> &GsiGridInfo; + fn lookup_grid_points(&self, ix: u32, iy: u32) -> f64; + + fn get_interpolated_value(&self, x: f64, y: f64) -> f64 { + use std::f64::NAN; + let grid = self.grid_info(); + let grid_x = (x - grid.x_min as f64) * (grid.x_denom as f64); + let grid_y = (y - grid.y_min as f64) * (grid.y_denom as f64); + if grid_x < 0.0 || grid_y < 0.0 { + return NAN; + } + + let ix = grid_x.floor() as u32; + let iy = grid_y.floor() as u32; + let x_residual = grid_x % ix as f64; + let y_residual = grid_y % iy as f64; + + if ix >= grid.x_num || iy >= grid.y_num { + NAN + } else { + let lookup_or_nan = |x, y, cond: bool| { + if cond { + self.lookup_grid_points(x, y) + } else { + NAN + } + }; + + bilinear( + x_residual, + y_residual, + self.lookup_grid_points(ix, iy), + lookup_or_nan(ix + 1, iy, ix < grid.x_num - 1), + lookup_or_nan(ix, iy + 1, iy < grid.y_num - 1), + lookup_or_nan(ix + 1, iy + 1, ix < grid.x_num - 1 && iy < grid.y_num - 1), + ) + } + } +} + +impl<'a> Grid for MemoryGrid<'a> { + fn grid_info(&self) -> &GsiGridInfo { + &self.grid_info + } + fn lookup_grid_points(&self, ix: u32, iy: u32) -> f64 { + match self.points[(self.grid_info.x_num * iy + ix) as usize] { + 9990000 => f64::NAN, + v => v as f64 / 10000.0, + } + } +} + +impl<'a> Geoid for MemoryGrid<'a> { + fn get_height(&self, lng: f64, lat: f64) -> f64 { + self.get_interpolated_value(lng, lat) + } +} + +impl<'a> MemoryGrid<'a> { + pub fn from_binary_reader(reader: &mut R) -> io::Result { + // Read header + let mut buf = [0; 28]; + reader.read_exact(&mut buf)?; + let grid_info = GsiGridInfo { + x_num: u16::from_le_bytes(buf[0..2].try_into().unwrap()) as u32, + y_num: u16::from_le_bytes(buf[2..4].try_into().unwrap()) as u32, + x_denom: u16::from_le_bytes(buf[4..6].try_into().unwrap()) as u32, + y_denom: u16::from_le_bytes(buf[6..8].try_into().unwrap()) as u32, + x_min: f32::from_le_bytes(buf[8..12].try_into().unwrap()), + y_min: f32::from_le_bytes(buf[12..16].try_into().unwrap()), + ikind: u16::from_le_bytes(buf[16..18].try_into().unwrap()), + version: String::from_utf8_lossy(&buf[18..28]).into(), + }; + + // Read grid point values + let mut points = Vec::with_capacity((grid_info.x_num * grid_info.y_num) as usize); + while points.len() < (grid_info.y_num * grid_info.x_num) as usize { + let mut buf = [0; 4]; + reader.read_exact(&mut buf)?; + points.push(i32::from_le_bytes(buf)); + } + + Ok(MemoryGrid { + grid_info, + points: points.into(), + }) + } + + pub fn to_binary_writer(&self, writer: &mut W) -> io::Result<()> { + // Write header + writer.write_all(&(self.grid_info.x_num as u16).to_le_bytes())?; + writer.write_all(&(self.grid_info.y_num as u16).to_le_bytes())?; + writer.write_all(&(self.grid_info.x_denom as u16).to_le_bytes())?; + writer.write_all(&(self.grid_info.y_denom as u16).to_le_bytes())?; + writer.write_all(&self.grid_info.x_min.to_le_bytes())?; + writer.write_all(&self.grid_info.y_min.to_le_bytes())?; + writer.write_all(&(self.grid_info.ikind).to_le_bytes())?; + if self.grid_info.version.len() > 10 { + panic!("version string must be less than 10 characters"); + } + writer.write_all(self.grid_info.version.as_bytes())?; + for _ in 0..10 - self.grid_info.version.len() { + writer.write_all(&[0])?; + } + + // Write grid point values + for p in self.points.iter() { + writer.write_all(&p.to_le_bytes())?; + } + Ok(()) + } + + pub fn from_ascii_reader(reader: &mut R) -> io::Result { + use io::{Error, ErrorKind::InvalidData}; + let mut reader = io::BufReader::new(reader); + let mut line = String::new(); + reader.read_line(&mut line)?; + + let c: Vec<&str> = line.split_whitespace().collect(); + if c.len() != 8 { + return Err(Error::new(InvalidData, "header line must have 8 values")); + } + if c[2] != "0.016667" { + return Err(Error::new( + InvalidData, + "latitude interval must be 0.016667", + )); + } + if c[3] != "0.025000" { + return Err(Error::new( + InvalidData, + "longitude interval must be 0.025000", + )); + } + + let grid_info = GsiGridInfo { + x_num: c[5] + .parse() + .map_err(|_| Error::new(InvalidData, "cannot parse header"))?, + y_num: c[4] + .parse() + .map_err(|_| Error::new(InvalidData, "cannot parse header"))?, + x_denom: 40, + y_denom: 60, + x_min: c[1] + .parse() + .map_err(|_| Error::new(InvalidData, "cannot parse header"))?, + y_min: c[0] + .parse() + .map_err(|_| Error::new(InvalidData, "cannot parse header"))?, + ikind: c[6] + .parse() + .map_err(|_| Error::new(InvalidData, "cannot parse header"))?, + version: c[7].to_string(), + }; + + let mut points = Vec::with_capacity((grid_info.x_num * grid_info.y_num) as usize); + for line_or_err in reader.lines() { + match line_or_err { + Ok(line) => { + for s in line.split_ascii_whitespace() { + let s = s.replace('.', ""); + let Ok(n) = s.parse::() else { + return Err(Error::new(InvalidData, "Invalid data")); + }; + points.push(n); + } + } + Err(e) => { + return Err(e); + } + } + } + + Ok(MemoryGrid { + grid_info, + points: points.into(), + }) + } +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..c4f389c --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,3 @@ +mod gsi; + +pub use gsi::*; diff --git a/tests/geoid_encrypted_for_ci.zip b/tests/geoid_encrypted_for_ci.zip new file mode 100644 index 0000000..f6331d5 Binary files /dev/null and b/tests/geoid_encrypted_for_ci.zip differ