|
| 1 | +package japangeoid |
| 2 | + |
| 3 | +import "math" |
| 4 | + |
| 5 | +// Grid はジオイドデータを表すインターフェースです。 |
| 6 | +type Grid interface { |
| 7 | + // GetHeight は指定した経度(lng)・緯度(lat)におけるジオイドの高さを返します。 |
| 8 | + GetHeight(lng, lat float64) float64 |
| 9 | +} |
| 10 | + |
| 11 | +// MemoryGrid はジオイドデータを保持する構造体です。 |
| 12 | +type MemoryGrid struct { |
| 13 | + Info GridInfo |
| 14 | + Points []int32 // グリッド点の高さ: [y_num * x_num] |
| 15 | +} |
| 16 | + |
| 17 | +// GetHeight は指定した経度(lng)・緯度(lat)におけるジオイドの高さを返します。無効な座標の場合は NaN を返します。 |
| 18 | +func (g *MemoryGrid) GetHeight(lng, lat float64) float64 { |
| 19 | + x, y, info := lng, lat, g.Info |
| 20 | + |
| 21 | + // 1. convert lng, lat to grid index |
| 22 | + gridX := (x - float64(info.XMin)) * float64(info.XDenom) |
| 23 | + gridY := (y - float64(info.YMin)) * float64(info.YDenom) |
| 24 | + if gridX < 0.0 || gridY < 0.0 { |
| 25 | + return math.NaN() |
| 26 | + } |
| 27 | + |
| 28 | + ix := int(math.Floor(gridX)) |
| 29 | + iy := int(math.Floor(gridY)) |
| 30 | + if ix < 0 || iy < 0 || ix >= int(info.XNum) || iy >= int(info.YNum) { |
| 31 | + return math.NaN() |
| 32 | + } |
| 33 | + |
| 34 | + // 2. get values of 4 points around |
| 35 | + v00 := g.lookupGridPoints(ix, iy) |
| 36 | + |
| 37 | + var v01 float64 |
| 38 | + if ix < int(info.XNum)-1 { |
| 39 | + v01 = g.lookupGridPoints(ix+1, iy) |
| 40 | + } else { |
| 41 | + v01 = math.NaN() |
| 42 | + } |
| 43 | + |
| 44 | + var v10 float64 |
| 45 | + if iy < int(info.YNum)-1 { |
| 46 | + v10 = g.lookupGridPoints(ix, iy+1) |
| 47 | + } else { |
| 48 | + v10 = math.NaN() |
| 49 | + } |
| 50 | + |
| 51 | + var v11 float64 |
| 52 | + if ix < int(info.XNum)-1 && iy < int(info.YNum)-1 { |
| 53 | + v11 = g.lookupGridPoints(ix+1, iy+1) |
| 54 | + } else { |
| 55 | + v11 = math.NaN() |
| 56 | + } |
| 57 | + |
| 58 | + // 3. apply bilinear interpolation |
| 59 | + xFrac := gridX - float64(ix) |
| 60 | + yFrac := gridY - float64(iy) |
| 61 | + return bilinear(xFrac, yFrac, v00, v01, v10, v11) |
| 62 | +} |
| 63 | + |
| 64 | +func (g *MemoryGrid) lookupGridPoints(ix, iy int) float64 { |
| 65 | + pos := iy*int(g.Info.XNum) + ix |
| 66 | + v := g.Points[pos] |
| 67 | + |
| 68 | + // 9990000 is NaN |
| 69 | + if v == 9990000 { |
| 70 | + return math.NaN() |
| 71 | + } |
| 72 | + |
| 73 | + return float64(v) / 10000.0 |
| 74 | +} |
0 commit comments