[Do Not Merge] GSIGEO2024

.gitignore

@@ -5,6 +5,7 @@ Cargo.lock
 wasm-pack.log
 gsigeo2011_ver2_2.*
 !gsigeo2011_ver2_2.bin.*
+GSIGEO*.isg
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

examples/convert_asc_to_bin_lz4.rs

@@ -1,14 +1,23 @@
+//! Convert .asc to .bin.lz4
+//!
+//! Usage:
+//! cargo run --example convert_asc_to_bin_lz4 -- input.asc
+
 use std::fs::{self, File};
 use std::io::{BufReader, Write};
+use std::path::Path;
 
 use japan_geoid::gsi::MemoryGrid;
 use japan_geoid::Geoid;
 
 fn main() -> std::io::Result<()> {
+    let argv = std::env::args().collect::<Vec<_>>();
+    let path = Path::new(&argv[1]);
+
     let (lng, lat) = (138.2839817085188, 37.12378643088312);
 
     // Load from the original ascii format.
-    let mut reader = BufReader::new(File::open("./gsigeo2011_ver2_2.asc")?);
+    let mut reader = BufReader::new(File::open(path)?);
     let geoid = MemoryGrid::from_ascii_reader(&mut reader)?;
 
     // Calculate the geoid height.
@@ -21,12 +30,11 @@ fn main() -> std::io::Result<()> {
     // Dump as the efficient binary format.
     let mut buf = Vec::new();
     geoid.to_binary_writer(&mut buf)?;
-    File::create("./gsigeo2011_ver2_2.bin.lz4")?
-        .write_all(&lz4_flex::compress_prepend_size(&buf))?;
+    let out_path = path.with_extension("bin.lz4");
+    File::create(&out_path)?.write_all(&lz4_flex::compress_prepend_size(&buf))?;
 
     // Load the binary model.
-    let decompressed =
-        &lz4_flex::decompress_size_prepended(&fs::read("./gsigeo2011_ver2_2.bin.lz4")?).unwrap();
+    let decompressed = &lz4_flex::decompress_size_prepended(&fs::read(&out_path)?).unwrap();
     let geoid = MemoryGrid::from_binary_reader(&mut std::io::Cursor::new(decompressed))?;
 
     let height = geoid.get_height(lng, lat);

isg2asc.py

@@ -0,0 +1,36 @@
+"""
+usage: isg2asc GSIGEO2024beta.isg
+"""
+
+import sys
+from pathlib import Path
+
+data_points = []
+
+filename = Path(sys.argv[1])
+
+with open(filename, "r") as f:
+    for line in f:
+        if line.strip().startswith("end_of_head ========"):
+            break
+
+    for line in f:
+        data_points.extend(v for v in line.split())
+
+NX = 1601
+NY = 2101
+
+assert len(data_points) == NX * NY
+
+with open(filename.with_suffix(".asc"), "w") as f:
+    f.write("15.00000 120.00000 0.016667 0.025000 2101 1601 1 ver-beta")
+    f.write("\n")
+    f.write("\n")
+    for i in range(NY):
+        for j in range(NX):
+            v = data_points[(NY - i - 1) * NX + j]
+            if v == "-9999.0000":
+                v = "999.0000"
+            f.write(v)
+            f.write(" ")
+        f.write("\n")

src/gsi.rs

@@ -19,8 +19,8 @@ pub trait Grid {
 
         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;
+        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
@@ -278,6 +278,24 @@ pub fn load_embedded_gsigeo2011() -> MemoryGrid<'static> {
     .unwrap()
 }
 
+/// Loads the embedded GSIGEO2024 Japan geoid model.
+///
+/// ```
+/// use japan_geoid::gsi::load_embedded_gsigeo2024;
+/// use japan_geoid::Geoid;
+///
+/// let geoid = load_embedded_gsigeo2024();
+/// let height = geoid.get_height(138.2839817085188, 37.12378643088312);
+/// assert!((height - 39.61005876363226).abs() < 1e-6)
+/// ```
+pub fn load_embedded_gsigeo2024() -> MemoryGrid<'static> {
+    const EMBEDDED_MODEL: &[u8] = include_bytes!("gsigeo2024_beta.bin.lz4");
+    MemoryGrid::from_binary_reader(&mut std::io::Cursor::new(
+        lz4_flex::decompress_size_prepended(EMBEDDED_MODEL).unwrap(),
+    ))
+    .unwrap()
+}
+
 #[cfg(test)]
 mod tests {
     use core::panic;
@@ -287,7 +305,7 @@ mod tests {
     use super::*;
 
     #[test]
-    fn embedded() {
+    fn embedded2011() {
         let geoid = load_embedded_gsigeo2011();
         let _ = format!("{:?}", geoid);
 
@@ -334,6 +352,62 @@ mod tests {
         assert_eq!(info.y_min, 20.0);
     }
 
+    #[test]
+    fn embedded2024() {
+        let geoid = load_embedded_gsigeo2024();
+        let _ = format!("{:?}", geoid);
+
+        let height = geoid.get_height(138.2839817085188, 37.12378643088312);
+        assert!((height - 39.61005876363226).abs() < 1e-6);
+
+        // compare with the results of GSI's 'geoidcalc' implementation
+        let height = geoid.get_height(140.085365000, 36.104394000);
+        assert!((height - 40.3059).abs() < 1e-4);
+        let height = geoid.get_height(139.615526456, 35.160410123);
+        assert!((height - 36.7568).abs() < 1e-4);
+        let height = geoid.get_height(138.215695342, 36.832842854);
+        assert!((height - 41.6041).abs() < 1e-4);
+        let height = geoid.get_height(130., 30.);
+        assert!((height - 30.5669).abs() < 1e-4);
+
+        let height = geoid.get_height(120.0, 15.0);
+        assert!(!f64::is_nan(height));
+        let height = geoid.get_height(120.0, 50.0);
+        assert!(!f64::is_nan(height));
+        let height = geoid.get_height(160.0, 15.0);
+        assert!(!f64::is_nan(height));
+        let height = geoid.get_height(160.0, 50.0);
+        assert!(!f64::is_nan(height));
+
+        let height = geoid.get_height(130.0, 15.0);
+        assert!(!f64::is_nan(height));
+        let height = geoid.get_height(130.0, 50.0);
+        assert!(!f64::is_nan(height));
+        let height = geoid.get_height(120.0, 20.0);
+        assert!(!f64::is_nan(height));
+        let height = geoid.get_height(160.0, 20.0);
+        assert!(!f64::is_nan(height));
+
+        let height = geoid.get_height(130.0, 14.99);
+        assert!(f64::is_nan(height));
+        let height = geoid.get_height(130.0, 50.01);
+        assert!(f64::is_nan(height));
+        let height = geoid.get_height(119.99, 20.0);
+        assert!(f64::is_nan(height));
+        let height = geoid.get_height(160.01, 20.0);
+        assert!(f64::is_nan(height));
+
+        let info = geoid.grid_info();
+        let _ = format!("{:?}", info);
+        assert_eq!(info.x_num, 1601);
+        assert_eq!(info.y_num, 2101);
+        assert_eq!(info.version, "ver-beta\0\0");
+        assert_eq!(info.x_denom, 40);
+        assert_eq!(info.y_denom, 60);
+        assert_eq!(info.x_min, 120.0);
+        assert_eq!(info.y_min, 15.0);
+    }
+
     #[test]
     fn ascii_to_binary() {
         // from ascii