Etape 10 – Modèle numérique du terrain

Un modèle numérique du terrain, abrégé MNT (digital elevation model, abrégé DEM en anglais) est une description numérique de la topographie d'une zone de la Terre.

Les MNTs utilisés dans le cadre de ce projet sont discrets, dans le sens où ils donnent l'altitude d'un nombre fini de points répartis à la surface de la Terre. La figure ci-dessous présente une topographie — en bleu — et dix points dont l'altitude — en rouge clair — pourrait constituer un MNT la représentant.

Le format HGT est un format de fichier binaire simple représentant un modèle numérique du terrain.

Un fichier HGT couvre toujours une zone de 1° de côté en longitude et en latitude. Le nom du fichier donne les coordonnées du coin sud-ouest de cette zone, dans le système WGS 84. Par exemple, un fichier HGT nommé N46E007.hgt décrit la zone dont le coin sud-ouest a les coordonnées (+7°, +46°). Son coin nord-est a donc les coordonnées (+8°, +47°). L'image ci-dessous montre la zone couverte par ce fichier.

Un fichier HGT contient les altitudes d'un certain nombre de points séparés les uns des autres par une distance angulaire constante et identique en longitude et en latitude, appelée sa résolution (angulaire) et notée \(\delta\) dans la figure ci-dessus.

Par exemple, un fichier HGT ayant une résolution de 1" (une seconde d'arc) donne les altitudes pour des points distants d'une seconde d'arc les uns des autres, aussi bien en longitude qu'en latitude. Sachant que chaque degré est composé de 60 minutes d'arc, composées chacune de 60 secondes, un tel fichier contient l'altitude de \((60\cdot 60+1)^2 = 12\,967\,201\) points. Le \(+1\) dans cette formule provient du fait qu'il y a chevauchement entre les fichiers HGT voisins, chacun partageant une ligne avec son voisin du dessus, une avec son voisin du dessous, une colonne avec son voisin de gauche et une avec son voisin de droite.

Les altitudes sont données en mètres, sous forme de valeurs entières de 16 bits. Chacune de ces valeurs est représentée par deux octets successifs dans le fichier, le premier contenant les 8 bits de poids fort, le second les 8 bits de poids faible. Par exemple, l'altitude 429 m est représentée par la séquence d'octets (01, AD) (en hexadécimal), étant donné que 429 en hexadécimal est 1AD.

La première altitude figurant dans un fichier HGT est celle de son coin haut-gauche (!). La seconde est celle de son voisin de droite, et ainsi de suite jusqu'à la fin de la première ligne. Viennent ensuite les altitudes de la seconde ligne, de gauche à droite, et ainsi de suite. Le point dont l'altitude figure en premier dans le fichier HGT de la figure ci-dessus est donc le point mis en évidence en rouge. Attention : ce point n'est pas celui qui donne son nom au fichier !

Comme nous le verrons à l'étape suivante, le dessin d'un relief ombré nécessite de connaître le vecteur normal à la surface de la Terre pour chaque pixel de l'image de la carte. Bien entendu, celui-ci se détermine en fonction d'un modèle numérique du terrain.

Pour calculer ce vecteur normal, nous utiliserons une méthode qui n'exploite pas au maximum les données du MNT mais qui a l'avantage d'être simple et suffisante pour nos besoins. Elle consiste à faire la moyenne entre deux vecteurs normaux correspondant aux quatre points du MNT voisins de celui qui nous intéresse.

La figure ci-dessous illustre comment cette technique permet de déterminer le vecteur, désigné par \(\vec{n}\), normal à la surface de la Terre au point \(P\). Au moyen des quatres points voisins dont le MNT donne l'altitude, on calcule deux vecteurs normaux \(\vec{n_1} = \vec{a}\times\vec{b}\) et \(\vec{n_2}=\vec{c}\times\vec{d}\), puis on utilise leur moyenne comme approximation du vecteur \(\vec{n}\).

Pour simplifier les calculs, on fait l'hypothèse que la distance séparant les points du MNT, désignée par \(s\) dans la figure ci-dessus, est la même dans les deux directions. Pour un fichier HGT cela est incorrect sauf à l'équateur mais cette approximation suffit sous nos latitudes. La formule suivante permet donc d'obtenir \(s\), en mètres, en fonction de la résolution angulaire \(\delta\) du fichier HGT (en radians) :

où \(R = 6\,378\,137\) est le rayon, en mètres, de la Terre supposée sphérique.

Les vecteurs \(\vec{a}\) à \(\vec{d}\) de la figure 3 sont donnés par :

où

et \(z_{i,j}\) est l'altitude, en mètres, donnée par le MNT pour le point d'index \((i, j)\).

Les deux vecteurs normaux \(\vec{n_1}\) et \(\vec{n_2}\) sont quant à deux donnés par :

et leur moyenne vaut :

Pour commencer, définissez une interface nommée p.ex. Earth et dont le seul but est de contenir des constantes liées à la Terre. La seule qui nous intéresse ici est son rayon en mètres, stocké dans un champ final et statique de l'interface nommé p.ex. RADIUS et valant 6'378'137.

Un modèle numérique du terrain est représenté par une interface, nommée p.ex. DigitalElevationModel. Cette interface étend l'interface AutoCloseable étant donné qu'il est probable qu'un MNT soit stocké dans un fichier et nécessite donc d'être fermé après utilisation.

En plus de la méthode close héritée de AutoCloseable, l'interface représentant un MNT est très simple et ne fournit qu'une seule méthode. Celle-ci, nommée p.ex. normalAt, prend un point en coordonnées WGS 84 en argument et retourne le vecteur normal à la Terre en ce point, déterminé au moyen de la technique approximative présentée ci-dessus. Cette méthode lève une exception si le point pour lequel la normale est demandé ne fait pas partie de la zone couverte par le MNT.

Bien entendu, le point en coordonnée WGS 84 a le type PointGeo défini à l'étape 1 tandis que le vecteur normal est un vecteur tridimensionnel du type défini à l'étape 9.

Un MNT stocké dans un fichier au format HGT est représenté par une classe, nommée p.ex. HGTDigitalElevationModel et implémentant bien entendu l'interface des MNT décrite ci-dessus.

Le constructeur de cette classe prend un unique argument, de type File, qui désigne le fichier HGT contenant le modèle numérique du terrain. Le nom de ce fichier, qui s'obtient au moyen de la méthode getName, est structuré de la manière suivante :

• la première lettre donne le signe de la latitude du coin sud-ouest et est soit N (pour nord, donc +) soit S (pour sud, donc –),
• les deux prochaines lettres donnent la valeur absolue de la latitude du coin sud-ouest,
• la quatrième lettre donne le signe de la longitude du coin sud-ouest et est soit E (pour est, donc +) soit W (pour ouest, donc –),
• les trois prochaines lettres donnent la valeur absolue de la longitude du coin sud-ouest,
• les quatre dernières lettres constituent l'extension du nom de fichier et valent toujours .hgt.

Ces conventions de nommage permettent de déterminer automatiquement la zone couverte par un fichier HGT. Sa résolution peut quant à elle être déterminée grâce à sa taille en octets, qui s'obtient avec la méthode length.

Bien entendu, le constructeur doit lever une exception si le nom du fichier HGT n'obéit pas aux conventions de nommage susmentionnées, ou si sa taille en octets divisée par deux n'a pas une racine carrée entière.

Pour lire les données d'un fichier HGT, deux solutions existent :

1. la totalité des données (c-à-d des altitudes) du fichier peuvent être lues en mémoire et stockées p.ex. dans un tableau bidimensionnel de type short[][], d'où elles sont ensuite lues,
2. le fichier peut être mappé en mémoire et les données directement lues depuis là.

La seconde solution est plus courte à mettre en œuvre et probablement plus efficace que la première, mais nécessite l'utilisation du concept de fichier mappé en mémoire — et des classes de la bibliothèque Java correspondantes — qui n'a pas été vu au cours.

Les personnes désirant adopter la première solution peuvent le faire simplement au moyen des concepts vus au cours et de la classe DataInputStream, dont la méthode readShort permet directement de lire une valeur de 16 bits du flot sous-jacent.

Les personnes désirant adopter la seconde solution pourront s'aider de la rapide introduction au concept de fichier mappé donnée ci-dessous, ainsi que de l'exemple de code qui suit.

Le « mappage » d'un fichier en mémoire est une opération permettant — conceptuellement en tout cas — de charger la totalité de ses octets en mémoire afin d'y accéder comme aux éléments d'un tableau, c'est-à-dire dans un ordre quelconque. Le mappage offre donc un accès aléatoire au contenu d'un fichier, contrairement aux flots d'entrée-sortie qui y offrent un accès séquentiel.

Comparé à un véritable tableau contenant les octets d'un fichier, un fichier mappé en mémoire (memory mapped file) possède toutefois deux avantages importants :

1. les octets du fichier ne sont chargés en mémoire que lorsqu'ils sont effectivement lus (ou écrits) par le programme,
2. les écritures faites en mémoire sont automatiquement répercutées sur le contenu du fichier stocké sur disque.

Le premier de ces deux avantages peut être intéressant pour ce projet, étant donné que les fichiers HGT sont gros mais que seul un petit sous-ensemble de leurs valeurs est généralement nécessaire pour le dessin d'une carte donnée.

Pour mapper un fichier en mémoire en Java, il faut tout d'abord obtenir un canal (de type FileChannel) pour ce fichier. Un moyen simple d'y parvenir consiste à d'abord obtenir un flot d'entrée sur le fichier, de type FileInputStream, puis d'utiliser la méthode getChannel de cette classe.

Une fois le canal obtenu, on peut utiliser sa méthode map pour mapper le fichier correspondant en mémoire. (Au passage, cette méthode map n'a absolument rien à voir avec la méthode du même nom de l'interface Stream.) Cette méthode retourne un objet de type ByteBuffer — en réalité de son sous-type MappedByteBuffer, mais c'est sans importance — qui représente (en gros) un tableau d'octets, ceux du fichier mappé.

Un tableau d'octet de type ByteBuffer peut être vu comme un tableau d'entiers de 16 bits de type ShortBuffer au moyen de sa méthode asShortBuffer, très utile ici étant donné que les fichiers HGT sont composés de tels entiers. Une fois cette vue obtenue, on peut l'utiliser pour lire un entier de 16 bits dont on connaît l'index au moyen de la méthode get.

Le programme d'exemple suivant mappe en mémoire le fichier HGT nommé N46E007.hgt et en affiche les onze premières altitudes, lues sous la forme de valeurs de type short. Etant donné ce qui a été dit ci-dessus, ces valeurs sont les altitudes des points aux coordonnées :

1. (7° 0' 0'', 47° 0' 0''),
2. (7° 0' 1'', 47° 0' 0''),
3. …,
11. (7° 0' 10'', 47° 0' 0'').

Tous ces points se trouvant dans le lac de Neuchâtel, et ce dernier étant à une altitude de 429 m, cette valeur est imprimée onze fois par le programme.

import java.io.File;
import java.io.IOException;
import java.io.FileInputStream;
import java.nio.ShortBuffer;
import java.nio.channels.FileChannel.MapMode;

public class MemMapExample {
    public static void main(String[] args)
            throws IOException {
        File file = new File("N46E007.hgt");
        long length = file.length();
        try (FileInputStream stream =
                 new FileInputStream(file)) {
            ShortBuffer buffer = stream.getChannel()
                .map(MapMode.READ_ONLY, 0, length)
                .asShortBuffer();

            for (int i = 0; i <= 10; ++i)
                System.out.println(buffer.get(i));
        }
    }
}

Pour tester votre code de lecture des fichiers HGT, nous vous fournissons une archive Zip contenant les fichiers HGT N46E006.hgt et N46E007.hgt qui, à eux deux, couvrent les zones correspondant aux données OSM fournies à l'étape 5 pour les villes de Lausanne, Berne et Interlaken. Ces fichiers proviennent du site Viewfinder Panoramas de Jonathan de Ferranti, sur lequel les personnes intéressées en trouveront d'autres.

Attention, cette archive est d'une taille assez conséquente (23 Mo) et vous êtes donc priés de la stocker une fois téléchargée.

Pour voir si les vecteur normaux que vous calculez sont valides, vous pouvez soit attendre l'étape suivante, soit écrire un petit programme principal produisant une image en niveaux de gris représentant l'une ou l'autre des composantes du vecteur normal, ou une combinaison d'entre-elles. Par exemple, on peut utiliser la composante \(y\) du vecteur normal normé, \(n_y\), pour attribuer à chaque pixel la couleur \((g, g, g)\) où \(g\) est donnée par : \(g = \tfrac{1}{2}(n_y + 1)\). Avec cette formule, les pentes verticales orientées plein nord sont blanches (car \(n_y = 1\)), celles orientées plein sud sont noires (car \(n_y = -1\)), et les autres adoptent différentes teintes de gris.

En utilisant cette technique pour dessiner une image de 800 pixels de côté représentant la zone dont le point bas-gauche a la position (7.2°, 46.2°) et le point haut-droite a la position (7.8°, 46.8°), on obtient l'image ci-dessous. Le lac de Thoune est bien visible en haut à droite, de même que la vallée du Rhone dans la partie inférieure.