MobilApp/lib/controls/geoid_grid.dart

68 lines
2.4 KiB
Dart

import 'dart:typed_data';
import 'package:flutter/services.dart';
class GeoidGrid {
final double southLat, westLon, latStep, lonStep;
final int nRows, nCols;
final Float32List grid; // big-endian, south-to-north
GeoidGrid._(this.southLat, this.westLon, this.latStep, this.lonStep,
this.nRows, this.nCols, this.grid);
static Future<GeoidGrid> load(String assetPath) async {
final bytes = (await rootBundle.load(assetPath)).buffer.asByteData();
final southLat = bytes.getFloat64(0, Endian.big);
final westLon = bytes.getFloat64(8, Endian.big);
final latStep = bytes.getFloat64(16, Endian.big);
final lonStep = bytes.getFloat64(24, Endian.big);
final nRows = bytes.getInt32(32, Endian.big);
final nCols = bytes.getInt32(36, Endian.big);
final grid = Float32List(nRows * nCols);
for (int i = 0; i < grid.length; i++) {
grid[i] = bytes.getFloat32(40 + i * 4, Endian.big);
}
return GeoidGrid._(southLat, westLon, latStep, lonStep, nRows, nCols, grid);
}
/// Returns EHT2014 geoid undulation [m] at (lat, lon) via bilinear interpolation.
/// Returns null if outside grid or in no-data area.
double? getUndulation(double lat, double lon) {
final rowF = (lat - southLat) / latStep;
final colF = (lon - westLon) / lonStep;
final r0 = rowF.floor();
final c0 = colF.floor();
if (r0 < 0 || r0 >= nRows - 1 || c0 < 0 || c0 >= nCols - 1) return null;
final dr = rowF - r0;
final dc = colF - c0;
double v(int r, int c) => grid[r * nCols + c].toDouble();
final noData = -80.0;
final v00 = v(r0, c0);
if (v00 < noData) return null;
final v10 = v(r0 + 1, c0);
if (v10 < noData) return null;
final v01 = v(r0, c0 + 1);
if (v01 < noData) return null;
final v11 = v(r0 + 1, c0 + 1);
if (v11 < noData) return null;
return v00 * (1 - dr) * (1 - dc) +
v10 * dr * (1 - dc) +
v01 * (1 - dr) * dc +
v11 * dr * dc;
}
/// Orthometric height in the Hungarian EHT2014 system.
/// [nmea_ortho] = GGA field 9 (altitudeAboveMeanSeaLevel)
/// [nmea_sep] = GGA field 11 (geoidSeparation from receiver)
double? toEovHeight(
double lat, double lon, double nmeaOrtho, double nmeaSep) {
final nEht2014 = getUndulation(lat, lon);
if (nEht2014 == null) return null;
final ellipsoidal = nmeaOrtho + nmeaSep; // h = H + N_receiver
return ellipsoidal - nEht2014; // H_EOV = h - N_EHT2014
}
}