feat: add ISG format support and geoid2024 model

Author: rot1024

.gitignore

@@ -1 +1,3 @@
 *.asc
+*.isg
+/conv

README.md

@@ -2,13 +2,40 @@
 
 日本のジオイド高を計算するための Go ライブラリ / Go library for calculating geoid height in Japan
 
-- 本ライブラリは国土地理院が提供するものではありません。
-- 本ライブラリは [ciscorn/japan-geoid](https://github.com/ciscorn/japan-geoid/tree/main) （MIT License）を参考に実装しています。
+- **外部依存ゼロ** - 標準ライブラリのみで動作します
+- 本ライブラリは国土地理院が提供するものではありません
+- 本ライブラリは [ciscorn/japan-geoid](https://github.com/ciscorn/japan-geoid/tree/main) （MIT License）を参考に実装しています
 
 ## 対応ジオイドモデル
 
-- 日本のジオイド 2011（Ver.2.2）（[出典](https://fgd.gsi.go.jp/download/geoid.php)）: `gsigeoid2011`
-- 日本のジオイド 2025（[出典](https://www.gsi.go.jp/buturisokuchi/grageo_reference.html)）: 正式公開後対応予定
+| パッケージ | 関数 | 内容 | サイズ |
+|-----------|------|------|--------|
+| `gsigeoid2011` | `Load()` | 日本のジオイド2011（Ver.2.2） | 約300KB |
+| `gsigeoid2024` | `Load()` | ジオイド2024＋基準面補正パラメータ（合成版、推奨） | 約400KB |
+| `gsigeoid2024` | `LoadHrefconv()` | 基準面補正パラメータのみ | 約15KB |
+| `gsigeoid2024geoid` | `Load()` | ジオイド2024のみ | 約3.7MB |
+
+### 日本のジオイド2011（Ver.2.2）
+
+- 出典: https://fgd.gsi.go.jp/download/geoid.php
+
+### ジオイド2024（日本とその周辺）＋基準面補正パラメータ
+
+- 出典: https://service.gsi.go.jp/kiban/app/geoid/
+
+ジオイド・モデル「ジオイド2024日本とその周辺」は、世界測地系（日本測地系2024）における座標値（緯度と経度）で示された任意の位置でのジオイド高を表したデータです。
+
+また、「基準面補正パラメータ」は、東京湾平均海面と離島独自の平均海面の差を基準面補正量と定義し、前述した任意の位置での基準面補正量を表したデータです。
+
+一部の離島において衛星測位によって標高を求める際には、ジオイド高と基準面補正量を使用する必要があります。
+
+ジオイドのみ（`gsigeoid2024geoid`）はサイズが大きいため、別パッケージとして分離しています。必要な場合のみインポートしてください。
+
+※ ジオイドのみのデータが大きい理由：合成版は海域がnodata扱いですが、ジオイドのみのデータは海域にも値が含まれており、圧縮効率が低くなっています。
+
+## 使い方
+
+### gsigeoid2011
 
 ## 使い方
 
@@ -20,26 +47,116 @@ import (
 	"github.com/eukarya-inc/japan-geoid-go/gsigeoid2011"
 )
 
-func Example() {
+func main() {
 	g, err := gsigeoid2011.Load()
 	if err != nil {
 		panic(err)
 	}
 
 	lng, lat := 138.2839817085188, 37.12378643088312
 	height := g.GetHeight(lng, lat)
-	fmt.Printf("Input: (lng: %.6f, lat: %.6f) -> Geoid height: %.6f\n", lng, lat, height)
+	fmt.Printf("Geoid height: %.6f\n", height)
+	// Output: Geoid height: 39.473871
+}
+```
+
+### gsigeoid2024
+
+```go
+package main
 
-	lng, lat = 10, 10
-	height = g.GetHeight(lng, lat)
-	fmt.Printf("Input: (lng: %.6f, lat: %.6f) -> Geoid height: %.6f\n", lng, lat, height)
+import (
+	"fmt"
+	"github.com/eukarya-inc/japan-geoid-go/gsigeoid2024"
+)
+
+func main() {
+	// ジオイド2024＋基準面補正パラメータ（合成版）
+	g, err := gsigeoid2024.Load()
+	if err != nil {
+		panic(err)
+	}
+
+	lng, lat := 138.2839817085188, 37.12378643088312
+	height := g.GetHeight(lng, lat)
+	fmt.Printf("Geoid height: %.6f\n", height)
+	// Output: Geoid height: 39.596702
 
-	// Output:
-	// Input: (lng: 138.283982, lat: 37.123786) -> Geoid height: 39.473871
-	// Input: (lng: 10.000000, lat: 10.000000) -> Geoid height: NaN
+	// 基準面補正パラメータのみ（離島で使用）
+	gHref, _ := gsigeoid2024.LoadHrefconv()
+	// 小笠原諸島（父島）付近
+	fmt.Printf("Hrefconv (Ogasawara): %.6f\n", gHref.GetHeight(142.19, 27.09))
+	// Output: Hrefconv (Ogasawara): 0.642000
 }
 ```
 
+### gsigeoid2024geoid（ジオイドのみ、約3.7MB）
+
+```go
+package main
+
+import (
+	"fmt"
+	"github.com/eukarya-inc/japan-geoid-go/gsigeoid2024geoid"
+)
+
+func main() {
+	g, err := gsigeoid2024geoid.Load()
+	if err != nil {
+		panic(err)
+	}
+
+	lng, lat := 138.2839817085188, 37.12378643088312
+	height := g.GetHeight(lng, lat)
+	fmt.Printf("Geoid height: %.6f\n", height)
+	// Output: Geoid height: 39.596702
+}
+```
+
+### ISGパーサー
+
+ISG形式のファイルを直接パースする場合は、`isg` パッケージを使用できます。
+
+```go
+package main
+
+import (
+	"fmt"
+	"os"
+
+	"github.com/eukarya-inc/japan-geoid-go/isg"
+)
+
+func main() {
+	f, _ := os.Open("geoid.isg")
+	defer f.Close()
+
+	grid, err := isg.Parse(f)
+	if err != nil {
+		panic(err)
+	}
+
+	fmt.Printf("Model: %s\n", grid.Header.ModelName)
+	fmt.Printf("Rows: %d, Cols: %d\n", grid.Header.NRows, grid.Header.NCols)
+}
+```
+
+## 補間アルゴリズム
+
+本ライブラリでは、グリッド点間のジオイド高を計算する際に **双線形補間（Bilinear Interpolation）** を使用しています。
+
+指定された座標を囲む4つのグリッド点の値から、以下の式で補間値を計算します：
+
+```
+h = h00 * (1-x) * (1-y) + h01 * x * (1-y) + h10 * (1-x) * y + h11 * x * y
+```
+
+ここで：
+- `h00`, `h01`, `h10`, `h11`: 4つのグリッド点のジオイド高
+- `x`, `y`: グリッド内での相対位置（0.0〜1.0）
+
+グリッド点がデータ範囲外または nodata の場合は `NaN` を返します。
+
 ## License
 
-[MIT Lisence](LICENSE)
+[MIT License](LICENSE)

cmd/conv.go

@@ -2,6 +2,7 @@ package main
 
 import (
 	"compress/gzip"
+	"fmt"
 	"io"
 	"net/http"
 	"os"
@@ -13,33 +14,48 @@ import (
 
 func main() {
 	if len(os.Args) < 2 {
-		println("Usage: conv <file>")
+		println("Usage: conv <file.asc|file.isg>")
 		os.Exit(1)
 	}
 
 	p := os.Args[1]
 	name := path.Base(p)
 
-	f, err := loadAsc(p)
+	f, err := loadFile(p)
 	if err != nil {
 		println(err.Error())
 		os.Exit(1)
 	}
-
 	defer f.Close()
 
-	grid, err := japangeoid.FromAsc(f)
-	if err != nil {
-		println(err.Error())
+	// 拡張子に応じてパーサーを選択
+	var grid *japangeoid.MemoryGrid
+	switch {
+	case strings.HasSuffix(strings.ToLower(name), ".isg"):
+		grid, err = japangeoid.FromIsg(f)
+		if err != nil {
+			println(err.Error())
+			os.Exit(1)
+		}
+		name = strings.TrimSuffix(name, ".isg")
+	case strings.HasSuffix(strings.ToLower(name), ".asc"):
+		grid, err = japangeoid.FromAsc(f)
+		if err != nil {
+			println(err.Error())
+			os.Exit(1)
+		}
+		name = strings.TrimSuffix(name, ".asc")
+	default:
+		println("Unknown file extension. Supported: .asc, .isg")
 		os.Exit(1)
 	}
 
-	f2, err := os.Create(strings.TrimSuffix(name, ".asc") + ".bin.gz")
+	outName := name + ".bin.gz"
+	f2, err := os.Create(outName)
 	if err != nil {
 		println(err.Error())
 		os.Exit(1)
 	}
-
 	defer f2.Close()
 
 	w := gzip.NewWriter(f2)
@@ -50,10 +66,11 @@ func main() {
 		os.Exit(1)
 	}
 
-	println("Done")
+	fmt.Printf("Done: %s (XNum=%d, YNum=%d, XMin=%.2f, YMin=%.2f)\n",
+		outName, grid.Info.XNum, grid.Info.YNum, grid.Info.XMin, grid.Info.YMin)
 }
 
-func loadAsc(url string) (io.ReadCloser, error) {
+func loadFile(url string) (io.ReadCloser, error) {
 	if strings.HasPrefix(url, "http://") || strings.HasPrefix(url, "https://") {
 		req, err := http.NewRequest("GET", url, nil)
 		if err != nil {

gsigeoid2024/Hrefconv2024.bin.gz

@@ -0,0 +1,38 @@
+package gsigeoid2024
+
+import (
+	"bytes"
+	"compress/gzip"
+	_ "embed"
+	"fmt"
+
+	japangeoid "github.com/eukarya-inc/japan-geoid-go"
+)
+
+//go:embed JPGEO2024+Hrefconv2024.bin.gz
+var combinedRaw []byte
+
+//go:embed Hrefconv2024.bin.gz
+var hrefconvRaw []byte
+
+// Load はジオイド2024と基準面補正パラメータを合成したデータを読み込みます。
+// 一部の離島において衛星測位によって標高を求める際には、この合成データを使用してください。
+func Load() (*japangeoid.MemoryGrid, error) {
+	return load(combinedRaw)
+}
+
+// LoadHrefconv は基準面補正パラメータのみを読み込みます。
+// 東京湾平均海面と離島独自の平均海面の差（基準面補正量）を表したデータです。
+func LoadHrefconv() (*japangeoid.MemoryGrid, error) {
+	return load(hrefconvRaw)
+}
+
+func load(data []byte) (*japangeoid.MemoryGrid, error) {
+	b := bytes.NewReader(data)
+	g, err := gzip.NewReader(b)
+	if err != nil {
+		return nil, fmt.Errorf("failed to create gzip reader: %w", err)
+	}
+
+	return japangeoid.FromBinary(g)
+}

gsigeoid2024/JPGEO2024+Hrefconv2024.bin.gz

@@ -0,0 +1,92 @@
+package gsigeoid2024
+
+import (
+	"fmt"
+	"math"
+	"testing"
+)
+
+func TestLoad(t *testing.T) {
+	g, err := Load()
+	if err != nil {
+		t.Fatalf("failed to load: %v", err)
+	}
+
+	// グリッド情報の確認
+	if g.Info.XNum != 1601 {
+		t.Errorf("XNum: got %d, want 1601", g.Info.XNum)
+	}
+	if g.Info.YNum != 2101 {
+		t.Errorf("YNum: got %d, want 2101", g.Info.YNum)
+	}
+
+	// 日本国内の座標でジオイド高を取得できることを確認
+	lng, lat := 138.2839817085188, 37.12378643088312
+	height := g.GetHeight(lng, lat)
+	if math.IsNaN(height) {
+		t.Errorf("GetHeight(%f, %f): got NaN, want valid height", lng, lat)
+	}
+	t.Logf("Load: GetHeight(%f, %f) = %f", lng, lat, height)
+
+	// 範囲外の座標はNaNを返すことを確認
+	lng, lat = 10, 10
+	height = g.GetHeight(lng, lat)
+	if !math.IsNaN(height) {
+		t.Errorf("GetHeight(%f, %f): got %f, want NaN", lng, lat, height)
+	}
+}
+
+func TestLoadHrefconv(t *testing.T) {
+	g, err := LoadHrefconv()
+	if err != nil {
+		t.Fatalf("failed to load: %v", err)
+	}
+
+	// グリッド情報の確認
+	if g.Info.XNum != 1601 {
+		t.Errorf("XNum: got %d, want 1601", g.Info.XNum)
+	}
+	if g.Info.YNum != 2101 {
+		t.Errorf("YNum: got %d, want 2101", g.Info.YNum)
+	}
+
+	// 日本国内の座標で値を取得できることを確認
+	lng, lat := 138.2839817085188, 37.12378643088312
+	height := g.GetHeight(lng, lat)
+	// Hrefconvは本土ではNaNになる可能性があるのでログ出力のみ
+	t.Logf("LoadHrefconv: GetHeight(%f, %f) = %f", lng, lat, height)
+}
+
+func Example() {
+	g, err := Load()
+	if err != nil {
+		panic(err)
+	}
+
+	lng, lat := 138.2839817085188, 37.12378643088312
+	height := g.GetHeight(lng, lat)
+	fmt.Printf("Input: (lng: %.6f, lat: %.6f) -> Geoid height: %.6f\n", lng, lat, height)
+
+	lng, lat = 10, 10
+	height = g.GetHeight(lng, lat)
+	fmt.Printf("Input: (lng: %.6f, lat: %.6f) -> Geoid height: %.6f\n", lng, lat, height)
+
+	// Output:
+	// Input: (lng: 138.283982, lat: 37.123786) -> Geoid height: 39.596702
+	// Input: (lng: 10.000000, lat: 10.000000) -> Geoid height: NaN
+}
+
+func ExampleLoadHrefconv() {
+	g, err := LoadHrefconv()
+	if err != nil {
+		panic(err)
+	}
+
+	// 小笠原諸島（父島）付近 - 基準面補正量がある
+	lng, lat := 142.19, 27.09
+	height := g.GetHeight(lng, lat)
+	fmt.Printf("Hrefconv: %.6f\n", height)
+
+	// Output:
+	// Hrefconv: 0.642000
+}