Per passare da quote ellissoidiche GPS a quote sul geoide ci sono varie opzioni in base al modello scelto
 |
| ITG2009 |
Se si usa EGM 2008 si ha un dato ogni 2.5 minuti di arco con una accuratezza dichiarata di 10 cm
Per effettuare la correzione di deve scaricare il file egm2008-2.5.pgm
https://sourceforge.net/projects/geographiclib/files/geoids-distrib/
e lanciare lo script indicando le coordinate GPS
from pygeodesy.geoids import GeoidPGM
geoid = GeoidPGM(r'.\geoids\egm2008-2_5.pgm', kind=3)
lat, lon = 43.7332442459068, 10.2760008626873
#lat, lon = 43.7337747, 10.2762272 # San Rossore
h_ell = 51.152
N = geoid.height(lat, lon)
H_ortho = h_ell - N
print(f"Differenza ellissoide ortometrica = {N:.3f} m")
print(f"H ortometrica = {H_ortho:.3f} m")
altrimenti si puo' usare ITG 2009 specifico per Italia e fornito da Polimi con spaziatura 1.5 minuti primi
https://www.isgeoid.polimi.it/Geoid/Europe/Italy/ITG2009_g.html
"""
Parser per file .isg (International Service for the Geoid, Politecnico di Milano)
Formato griglia regolare ASCII: header + valori di ondulazione N riga per riga,
ordinamento default: nord -> sud, ogni riga ovest -> est.
Uso:
geoid = IsgGeoid("ITG2009_20170606_isg.txt")
N = geoid.height(lat=43.6833, lon=10.2833)
H_ortho = h_ellissoidica - N
"""
import numpy as np
from scipy.interpolate import RegularGridInterpolator
import os
class IsgGeoid:
def __init__(self, path):
with open(path, encoding="utf-8", errors="replace") as f:
lines = f.readlines()
header = {}
end_idx = None
for i, line in enumerate(lines):
if "end_of_head" in line:
end_idx = i
break
if "=" in line and "begin_of_head" not in line:
key, val = line.split("=", 1)
header[key.strip()] = val.strip()
if end_idx is None:
raise ValueError("Non trovo 'end_of_head' nel file .isg")
self.lat_min = float(header["lat min"])
self.lat_max = float(header["lat max"])
self.lon_min = float(header["lon min"])
self.lon_max = float(header["lon max"])
self.dlat = float(header["delta lat"])
self.dlon = float(header["delta lon"])
self.nrows = int(header["nrows"])
self.ncols = int(header["ncols"])
self.nodata = float(header.get("nodata", -9999.0))
data_lines = lines[end_idx + 1:]
vals = []
for line in data_lines:
line = line.strip()
if line:
vals.extend(float(x) for x in line.split())
expected = self.nrows * self.ncols
if len(vals) != expected:
raise ValueError(
f"Numero di valori letti ({len(vals)}) diverso da nrows*ncols "
f"({expected}). File troncato o formato inatteso."
)
grid = np.array(vals, dtype=np.float64).reshape(self.nrows, self.ncols)
grid[grid == self.nodata] = np.nan
# riga 0 = lat_max, riga -1 = lat_min (ordinamento default ISG: nord->sud)
# capovolgo per avere lat crescente, richiesto da RegularGridInterpolator
grid = grid[::-1, :]
lat_arr = np.linspace(self.lat_min, self.lat_max, self.nrows)
lon_arr = np.linspace(self.lon_min, self.lon_max, self.ncols)
self._interp = RegularGridInterpolator(
(lat_arr, lon_arr), grid, method="linear", bounds_error=True
)
def height(self, lat, lon):
"""Ritorna l'ondulazione del geoide N (metri) al punto lat/lon (gradi decimali)."""
val = self._interp([[lat, lon]])[0]
if np.isnan(val):
raise ValueError(f"Punto ({lat}, {lon}) cade in area nodata della griglia")
return float(val)
def height_array(self, lats, lons):
"""Versione vettorizzata: lats, lons array numpy della stessa lunghezza."""
pts = np.column_stack([lats, lons])
vals = self._interp(pts)
if np.any(np.isnan(vals)):
raise ValueError("Alcuni punti cadono in area nodata della griglia")
return vals
if __name__ == "__main__":
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
path = os.path.join(SCRIPT_DIR, "ITG2009_20170606.isg.txt")
geoid = IsgGeoid(path)
punti = {
"SanRossore": (43.7332442459068, 10.2760008626873),
}
# esempio conversione completa
h_ell = 51.152 # quota ellissoidica misurata dal drone (esempio)
lat, lon = punti["SanRossore"]
N = geoid.height(lat, lon)
print(N)
H_ortho = h_ell - N
print(f"\nEsempio: h_ell={h_ell} m -> H_ortho={H_ortho:.3f} m (N={N:.3f} m)")
modelli piu' grezzi ASTER GDEM v3 con passo di campionamento ad 1 grado
per confronto
Punto: lat=43.7332442459068, lon=10.2760008626873
ASTER GDEM V3 (Open Topo Data) : 7.00 m (riferito a EGM96, ~quota ortometrica/DSM)
Ondulazione geoide EGM2008 (N) : 46.5572 m
Ondulazione geoide EGM96 (N) : 46.9583 m
Ondulazione geoide EGM84 (N) : 47.5683 m
Quota ASTER ellissoidica WGS84 (via N_EGM96) : 53.96 m
Quota ASTER ellissoidica WGS84 (via N_EGM2008, non corretto) : 53.56 m
Differenza N(EGM96) - N(EGM2008) nel punto: +0.4011 m
Nessun commento:
Posta un commento