Una mesh a 30 metri, gratuita,

Lo scopo di questo articolo è fornire indicazioni su una fonte di dati satellitari, poco nota e gratuita, utili per creare un mesh terrain fino a 30 metri. Mi auguro che possa essere d’aiuto a coloro che, con ben più scarsi mezzi, ci hanno già dato ottimi mesh terrain dell’Italia.

Ciascuno di noi vorrebbe ritrovare nei suoi simulatori, almeno per i voli VFR, il territorio che meglio conosce e questo reso con il maggior dettaglio possibile.

Compreso, come mi è capitato di leggere, lo stato della flora nel vicinato di casa. Ovviamente ciò non è possibile!

Anche solo per quanto riguarda l’orografia, la rappresentazione, incluse le migliori mesh in circolazione, è tutt’altro che ottimale.

Provengo dalla provincia di Verona, quella parte immersa nelle morene glaciali compresa tra il lago di Garda e la valle dell’Adige, e non sono mai rimasto soddisfatto della rappresentazione della zona nelle varie mesh che mi sono capitate sotto mano.

Non una con la rocca di Garda svettante come dovrebbe e separata dall’antirocca, per non parlare della penisola di Sirmione, ridotta sempre ad una lingua di sabbia. E che dire della chiusa di Ceraino, che controllava la via del Tirolo (o della pianura Padana, a seconda del verso di percorrenza della valle dell’Adige)?.

Ho deciso quindi di "fabbricarmi" la mia mesh. Purtroppo gli ingradienti sono dei punti quota di posizione nota sul territorio, e ne servono parecchi per una mesh a 30m.

Dopo qualche girovagare ho trovato la soluzione sul sito (bellissimo) di John Childs (http://www.terrainmap.com), al quale vivamente rimando chi sia interessato a "giocare" con DEM e DTM.

Avendo subito scartato la soluzione di digitalizzare le curve di livello delle mappe della zona, non restava che cercare qualche altra fonte di DEM. Apparentemente i DEM disponibili liberamente per il territorio Italiano (di pagare non se ne parla) sono quelli della serie GTOPO30, o simili, con un punto quota a circa 500m alla nostra latitudine.

Apparentemente, però.

Sul sito di Childs è apparso, poco tempo fa, un articolo su una serie di dati ottenuti da un satellite dal nome gentile: ASTER.

Anche se creato per scopi diversi, dai sensori di questo satellite possono essere ricavati DEM a 30m per tutto il modo, e gratuiti! La terra promessa era di fronte a me, o quasi. Come ho avuto modo di scoprire in seguito.

I DEM a 30 metri ottenuti dai dati di ASTER hanno diversi lati negativi, il primo è sicuramente metterci le mani sopra.

Ancorché gratuiti, la procedura per ottenerli è quantomeno lunga, ed anche non poco complicata. Il secondo è la dimensione: 24Mb per un quadrato di circa 75 km di lato non sono pochi. Soprattutto se scaricato via telefono. Il terzo aspetto negativo, e questo e’ stato davvero una sorpresa, sono state le quote dei punti della mesh. A differenza dall’errore suggerito dai creatori (circa 30m), punti quota a 2200m erano finiti a 2600 ed oltre.

Si inizia registrandosi come utente presso il centro EOS Data Gateway che distribuisce i dati ottenuti da ASTER, si deve poi fare la login come l’utente che si è registrato e, nella pagina di ricerca dei dati selezionare come Data Set il set ASTER DEM PRODUCT V002. Si seleziona come fonte dei dati il Data Center EDC-ECS, e più avanti nella pagina si setta la latitudine e longitudine della zona di cui si è interessati ad avere il DEM. Ogni file di dati che include tale zona viene denominato "granulo" (data granule). Si avvia la ricerca con l’apposito bottone Start posto alla fine della pagina e si attende.

Se non viene trovato alcun granulo, allora per la zona in esame non è ancora stato creato alcun DEM e si deve ordinarne uno. Tranquilli, l’operazione è gratuita. Per ordinare il DEM bisogna sapere da quali dati deve essere creato, a questo fine si ripete l’operazione di ricerca ma selezionando, questa volta ASTER L1A RECONSTRUCTED UNPROCESSED INSTRUMENT DATA V002.

Se viene segnalato che sono stati trovati dei granuli questi vanno esaminati.

Controllare che il granulo abbia attributi (longitudine e latitudine) tali da comprendere la zona di interesse e che non abbia una percentuale di copertura nuvolosa troppo alta (cloud cover). Poiché i DEM sono ottenuti da queste immagini, maggiore la percentuale maggiori saranno le aree prive di dati di quota. Un valore limite è del 30%. E’ possibile vedere l’immagine L1A dalla quale verrà creato il DEM cliccando sul bottone "View Image". Quando si è deciso cosa "ordinare" cliccare sul bottone "Add to Cart" e premere poi il bottone "Order Options". Scegliere come media "FtpPull" che significa che, quando pronto, scaricherete il DEM via FTP anonimo. Cliccare su "OK, accept my choice and return to shopping cart" e all’apparire della schermata successiva su 'Go to Step 2: Order Form'.

Controllare i dati sulla scheda d’ordine e cliccare su 'Go to Step 3: Review Order Summary’ e poi su 'Go to Step 4: Submit Order'.

Se tutto è andato bene si vedrà una schermata 'Order Submitted', il DEM verrà creato e si avrà conferma del ricevimento dell’ordine di creazione via email.

Una seconda email ci avvertirà che il DEM è pronto. In quest’ultima ci sarà un indirizzo su cui fare un FTP anonimo e dal quale si potrà scaricare l’ASTER DEM nel formato HDF-EOS (non GeoTiff), assieme ad un file ASCII, avente suffisso .hdf.met. Una versione GeoTif del DEM verrà automaticamente messa nella directory pub/asterdem/relative su un diverso sito con indirizzo IP 152.61.128.25. Da questo può essere scaricata con un FTP anonimo cioè usando 'anonymous' come user id e il vostro indirizzo email come password. Ricordarsi di trasferire i dati in modalità binary (comando bin nell’ftp) e ricordarsi che il file GeoTif si chiamerà come il quello hdf indicato nella seconda email –estensione a parte.

Come già vi ho detto, una volta trovate le immagini LA1 della zona di interesse è bene guardarle e richiedere la produzione del DEM solo se non sono presenti troppe nubi.

Dato che il DEM viene prodotto a partire da queste immagini, dove ci sono nubi ci saranno delle zone prive di dati di quota. Tutte le zone con dati senza valore sono codificate come aventi una quota significativamente negativa, diversa da DEM a DEM (nel DEM utilizzato in questo esempio –250m). Questo porta ad avere dei DEM a macchia di leopardo e a dover integrare i dati mancanti pena antiestetiche cave sparse qua e là sul territorio. Dovendo integrare i dati mancanti, ho deciso di utilizzare un DEM dell’Europa della serie Gtopo30: il W020N90. Quest’ultimo, che utilizza un Datum WGS84, ha le quote coincidenti con quelle dello scenario originale di FS2002.

Per passare dalle quote date dall’Atlante Touring (basato su carte IGM che utilizzano un Datum diverso – mi pare IGS50), alle quote nel riferimento WGS84 bisogna aggiungere 16m. Il lago di Garda ha, infatti, una quota di 65m nell’atlante Touring ed una quota di 81m nel DEM Gtopo30.

Purtroppo la precisione di quota dei DEM ottenuti da dati ASTER è molto variabile.

Dovrebbe essere attorno ai 7-20m, ma in quello dell’esempio la vetta del monte Baldo era finita a 2600m (in realtà è circa 2200) mentre la vicina Garda sprofondava a –1m, fronte lago (+69m secondo l’Atlante Automobilistico del Touring Club Italiano). Riva, all’estremo Nord del lago, era a +83m (quindi ad una quota ragionevole rispetto ai +81m del lago) e Sirmione, al capo opposto, era finita a –54m. Per un secondo DEM, relativo alla zona dal basso lago al Po le quote si sono dimostrate, invece, perfettamente coincidenti con quelle dell’Atlante Touring. Quindi corrette l’una relativamente alle altre, ma con una scarto di –15m rispetto a quelle che sarebbero dovute essere se riferite allo zero WGS84.

Peccato, perché la precisione ci sarebbe. Nella prima immagine seguente ho riportato la zona attorno a S. Benedetto Po (il quadrato D1 della Tavola 47 dell’Atlante Touring, per chi si volesse prendere la briga di controllare). Cosa pensate che sia la striscia quasi verticale piu’ chiara evidenziata con le frecce? Ebbene si, è il rilevato dell’autostrada che scavalca il Po!

Lascio alla vostra curiosita’ scoprire a cosa corrispondano gli altri rilevati evidenziati nella seconda figura, che potrete ritrovare nella Tavola 31 del sopra citato Atlante.

In buona sostanza, le operazioni da fare per poter utilizzare il DEM sono al più quelle di una scalatura, una distorsione delle quote per portare in piano i punti a quota costante (linea di costa del lago), un riempimento dei buchi ed una allineamento con un DEM di maggiore superficie che permetta di ritagliare il rettangolo di dati da dare in pasto al Terrain SDK di FS2002 (o FS2000).

Il programma RESAMPLE del Terrain SDK accetta DEM in formato binario, con le quote espresse da interi a 16bit con segno e disposte nel file sequenzialmente, riga per riga dall’alto (Nord) verso il basso (Sud). Nel file dell’esempio, che ha 2531 valori per riga e 2433 righe, la quota dell’angolo NW del DEM corrisponde al primo valore nel file, quella dell’angolo NE al 2531- esimo valore, quello dell’angolo SW al (2531*2431+1) esimo valore e cosi’ via. Un file di questo tipo ha molti nomi, uno per ogni programma che manipola DEM che esiste in giro. Io lo chiamerò file BSQ.

Dato che i DEM ASTER vengono forniti anche in formato GeoTiff, utilizzeremo quest’ultimo. Anche perché i dati in formato hdf restano per soli pochi giorni sul server e se si è solo un poco distratti... I file GeoTiff sono invece sempre accessibili con un FTP anonimo.

Scaricato il file GeoTiff, si tratta di trasformare il file GeoTiff in un file BSQ, per questo viene utile il programma MICRODEM (che si può scaricare gratuitamente dalla rete). Io ho usato la versione 5, con la quale ho ottenuto le figure da 1 a 5.

Si apre il GeoTiff da MICRODEM, forzandolo ad aprire un DEM e dandogli a mano il nome del file GeoTiff. Si salva il GeoTiff come DEM nel formato proprietario di MICRODEM e si utilizza il programma MDEM2BSQ.EXE di John Childs (la URL ve l’ho già data) che si macina il DEM MICRODEM e sforna un DEM BSQ. Ed ecco spiegato da dove ho preso il nome del file dati da dare in pasto al TerrainSDK!

A questo punto il più è fatto, a parte leggersi la documentazione del TerrainSDK ovviamente. Ma a voi che mi avete seguito fin qui darò direttamente il file .INF necessario a RESAMPLE.EXE.

Come potete vedere dalla figura 3, una buona parte del DEM è priva di dati. Per utilizzare comunque tutti i dati presenti, ho espanso il DEM della figura 4, interpolando, in modo da portarlo ad avere gli stessi punti del DEM di figura precedente. Poi glielo ho sovrapposto distorcendolo in modo da far tornare il più possibile alcuni punti notevoli (vette e valli) comuni ad entrambi i DEM. Come alternativa, ho provato anche l’operazione opposta, cioè riempire le zone senza dati attorno al DEM Aster distorcendo il DEM Gtopo30. Nel primo caso mi sembra si abbia il migliore inserimento nel terreno di serie, nel secondo si ha una più precisa resa della posizione dei dettagli topografici. Ad esempio la fine del lago verso NE.

I DEM ASTER non sono esattamente prefetti, hanno ancora molti artefatti, principalmente picchi. In media un venti o trenta per DEM. Nel mio esempio quelli presenti non sono stati corretti, quindi occhio a non sbatterci dentro! Un ultimo appunto prima del volo; la posizione del lago di Garda (o meglio del DEM) non è ancora perfetta. Confido che mi scuserete ed aspetto che qualcuno, là fuori, la corregga.

I due DEM ASTER di figura 3 e figura 5 sono, rispettivamente: ASTER_DEM20020319111246.tif e ASTER_DEM20020401130738.tif.

Il file INF che ho utilizzato, e che vi ho promesso, ha il seguente aspetto:

[Destination]

; LOD = (Auto, 0..n)

; LOD = Auto tells the resampler to pick the closest LOD to the source

; cell size. Specifying a number will force an LOD. This can

; be used for undersampling of oversampling.

; UseSourceDimension = 1 means to retain le dimension specified in the

; [Source] section oterhwise specify the bounding box to be resapled wiht

; the NothLat, SouthLat, EastLong and WestLong keywords

;

; Aster limits: Lat: 46.0762, 46.0612, 45.4048, 45.4194

; Lon: 10.2586, 11.2400, 11.2139, 10.2439

; dlat=46.0762-45.4194=0.6568

; dlat=46.0612-45.4048=0.6564

; dlon=10.2586-11.2400=0.9814

; dlon=11.2139-10.2439=0.9700

: CellXdimensionDeg = 3.87751876729E-4

; CellYdimensionDeg = 2.69954788327E-4

; CellXdimensionDeg = 3.83247728171E-4 ; .9700/2531

; CellYdimensionDeg = 2.69790382244E-4 ; .6564/2433

; CellXdimensionDeg = 3.83399209486E-4 ; .9700/2530

; CellYdimensionDeg = 2.69901315789E-4 ; .6564/2432

; cells: grid spacind 30 m.

;

; max row 2530

; max col 2432

;

LOD = 10

DestDir = "tmf"

DestBaseFileName = "garda_aster"

UseSourceDimensions = 1

;NorthLat =

;SouthLat =

;EasLong =

;WestLong =

;

;FractionBits = 0

;BaseValue = 0

[Source]

; Type = (ElevS16LSB, ElevS16MSB)

; ElevS16LSB is the normal unless your data is in big endian format;

; ElevS16LSB means 16 bit signed elevation data with the least significant

; byte first, ElevS16MSB as before but most significant first

; Lat and Lon are the latidute and longitude of the NW corner of the bounding area

; BaseValue is the height, in meters, corresponding to a data value of zero

Type = ElevS16LSB

SourceDir = "aster"

SourceFile = "garda_aster.bsq"

CellType = S16LSB

NullCellValue = -250

Lat = 46.0762

Lon = 10.2586

NumOfCellsPerLine = 2531

NumOfLines = 2433

CellXdimensionDeg = 3.83247728100E-4 ; circa .9700/2531

CellYdimensionDeg = 2.69790382244E-4 ; .6564/2433

ScaleinMeters = 1.0

;FractionBits = 0

BaseValue = 0.

Per chi fosse interessato a dare una veste fotorealistica, anche se a falsi colori, alla mesh appena creata, consiglio il seguente articolo: ASTER L1A/L1B False Color Composites.

Con questo concludo, sperando di essere stato sufficientemente descrittivo su come creare le mesh partendo dai dati ASTER e sperando, anche, di aver suscitato il vostro interesse. Happy Flying to everybody!

Luca Martinelli