Modèles HGT et profils altimétriques
Etape 4

Nous mettons à votre disposition une archive Zip contenant 18 fichiers HGT couvrant la zone allant de 6° à 12° de longitude, et de 45° à 48° de latitude. Cette archive est à importer dans votre projet comme d'habitude. La zone couverte par les fichiers qu'elle contient est illustrée ci-dessous.

Tous ces fichiers HGT proviennent du site Viewfinder Panoramas de Jonathan de Ferranti. Les personnes intéressées à en télécharger d'autres les trouveront sur la page de téléchargement.

La classe File du paquetage java.io représente un fichier. Elle possède un certain nombre de méthodes permettant d'obtenir des informations sur le fichier, dont deux sont utiles à cette étape :

1. getName, qui retourne le nom du fichier sous forme de chaîne de caractères,
2. length, qui retourne la longueur du fichier, en octets.

La classe File possède plusieurs constructeurs, mais le seul utile à cette étape est celui prenant le nom du fichier en argument, sous la forme d'une chaîne de caractères.

Notez qu'une instance de la classe File peut être passée à l'un des constructeurs de la classe FileInputStream pour créer un flot d'entrée d'octets pour ce fichier, comme l'illustre l'exemple de la section 2.1.4 ci-dessous.

Pour extraire les différentes parties du nom du fichier HGT, vous pouvez utiliser les méthode charAt et substring de la classe String.

Pour convertir les sous-chaînes représentant les longitudes et latitudes du coin sud-ouest du fichier en entiers, utilisez la méthode parseInt de la classe Integer.

Le « mappage » d'un fichier en mémoire est une opération permettant — conceptuellement en tout cas — de charger d'un coup la totalité de son contenu en mémoire afin d'y accéder comme à 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 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. le contenu du fichier n'est lu que lorsque le programme l'utilise effectivement,
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'un panorama donné.

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. Cette méthode retourne un objet de type ByteBuffer 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°00′00″, 47°00′00″),
2. (7°00′01″, 47°00′00″),
3. …,
11. (7°00′10″, 47°00′00″).

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 f = new File("N46E007.hgt");
    long l = f.length();
    try (FileInputStream s = new FileInputStream(f)) {
      ShortBuffer b = s.getChannel()
        .map(MapMode.READ_ONLY, 0, l)
        .asShortBuffer();

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

Les méthodes elevationAt, slopeAt et positionAt doivent toutes trois calculer les coordonnées (longitude et latitude) d'un point à une position donnée du profil. Pour ce faire, elles pourraient bien entendu appliquer simplement les formules données à la section 1.3.

Malheureusement, ces formules sont assez coûteuses à calculer, car elles utilisent un grand nombre de fonction trigonométriques. Ce coût ne serait pas perceptible par l'utilisateur final si les méthodes susmentionnées n'étaient utilisées que rarement par le programme. Malheureusement, comme nous le verrons plus tard, tel n'est pas le cas. Par exemple, le dessin du panorama donné en exemple dans le document d'introduction au projet nécessite plus de 20 millions d'appels à la méthode positionAt.

Dès lors, il peut être utile d'optimiser ces méthodes. Une manière simple de le faire consiste à utiliser les formules de la section 1.3 pour calculer une fois pour toutes, à la construction d'un profil, la position d'un certain nombre de points régulièrement espacés (p.ex. tous les 4 km environ) le long de la trajectoire du grand cercle correspondant, puis d'interpoler linéairement entre ces positions par la suite.

Par exemple, imaginons qu'une instance de ElevationProfile soit construite ainsi :

ContinuousElevationModel dem = …;
GeoPoint origin =
  new GeoPoint(toRadians(6), toRadians(46));
double azimuth = toRadians(45);
double length = 100_000;

ElevationProfile profile =
  new ElevationProfile(dem, origin, azimuth, length);

Le constructeur de ElevationProfile calcule, au moyen des formules de la section 1.3, la position (longitude et latitude) de points régulièrement espacés de 4096 m les uns des autres. Etant donné que la longueur du profil est de 100 km, il y en a 26. La table ci-dessous donne la position des quatre premiers et du dernier point ainsi calculés :

Une fois ces valeurs calculées et stockées, la méthode positionAt — qui est aussi appelée par elevationAt et slopeAt — les utilise pour déterminer, par interpolation linéaire, la position des points intermédiaires. Par exemple, pour l'appel suivant :

double e = profile.elevationAt(10240);

le résultat est calculé par interpolation linéaire entre les positions des points d'index 2 et 3, étant donné que 10240 se trouve entre 8192 et 12288. Comme 10240 se trouve exactement au milieu de 8192 et 12288, les coordonnées sont déterminées ainsi par l'interpolation :

Comme on peut le voir en comparant ces valeurs à celles obtenues au moyen des formules exactes à la section 1.3, l'erreur introduite par l'interpolation est très faible, de l'ordre du mètre.

Notez que l'espacement des points choisi, 4096 m, est un peu arbitraire mais convient bien pour nos latitudes, dans le sens où il permet de beaucoup accélérer les calculs sans que l'erreur due à l'interpolation ne soit visible sur les images finales. De plus, 4096 est une puissance de deux, ce qui permet d'optimiser encore la méthode positionAt en remplaçant une division par un appel à la méthode scalb, passablement plus rapide. Vous êtes libres de mettre en œuvre cette optimisation si vous le désirez — et si vous la comprenez.

Pour vérifier que les données des fichiers HGT sont bien chargées, une technique consiste à générer une image se basant sur les altitudes des points et de voir si elle correspond à la réalité. Une telle image est similaire à celle produite par le programme DrawDEM de l'étape précédente.

Le programme ci-dessous fait exactement cela, en générant une image carrée de 300 pixels de côté représentant une zone d'un demi degré dont le coin sud-ouest a les coordonnées (6.25°, 46.25°) et le coin nord-est les coordonnées (6.75°, 46.75°). L'extrait de carte ci-dessous montre cette zone.

Notez que la raison pour laquelle cette zone apparaît rectangulaire ici mais carré plus bas est que cette carte utilise la projection de Mercator alors que le programme plus bas utilise une projection dite « plate carrée ».

Pour alléger la présentation du programme, les énoncés d'importation (import) ont été omis mais sont similaires à ceux du programme d'exemple de l'étape précédente. Pour la même raison, la méthode gray, en tout point identique à celle de l'étape précédente, a également été omise.

final class DrawHgtDEM {
  final static File HGT_FILE = new File("N46E006.hgt");
  final static double ORIGIN_LON = toRadians(6.25);
  final static double ORIGIN_LAT = toRadians(46.25);
  final static double WIDTH = toRadians(0.5);
  final static int IMAGE_SIZE = 300;
  final static double MIN_ELEVATION = 200;
  final static double MAX_ELEVATION = 1_500;

  public static void main(String[] as) throws Exception {
    DiscreteElevationModel dDEM =
      new HgtDiscreteElevationModel(HGT_FILE);
    ContinuousElevationModel cDEM =
      new ContinuousElevationModel(dDEM);

    double step = WIDTH / (IMAGE_SIZE - 1);
    BufferedImage i = new BufferedImage(IMAGE_SIZE,
                                        IMAGE_SIZE,
                                        TYPE_INT_RGB);
    for (int x = 0; x < IMAGE_SIZE; ++x) {
      double lon = ORIGIN_LON + x * step;
      for (int y = 0; y < IMAGE_SIZE; ++y) {
        double lat = ORIGIN_LAT + y * step;
        GeoPoint p = new GeoPoint(lon, lat);
        double el =
          (cDEM.elevationAt(p) - MIN_ELEVATION)
          / (MAX_ELEVATION - MIN_ELEVATION);
        i.setRGB(x, IMAGE_SIZE - 1 - y, gray(el));
      }
    }
    dDEM.close();

    ImageIO.write(i, "png", new File("dem.png"));
  }
}

En exécutant ce programme, on obtient l'image ci-dessous, sur laquelle le Léman est bien visible dans la moitié inférieure. On devine aussi le lac de Joux en haut à gauche, et les vallées creusées par les principaux cours d'eau apparaissent clairement.

Pour vérifier le calcul des profils altimétriques, il est aussi possible de les dessiner.

Le programme ci-dessous dessine le même profil que celui de la figure 3 — qui part des Pierres du Niton et va jusqu'à la gare de la Chaux-de-Fonds — mais en utilisant les données des fichiers HGT fournis. Une fois encore, les énoncés d'importation ont été omis pour alléger la présentation.

final class DrawElevationProfile {
  final static File HGT_FILE = new File("N46E006.hgt");
  final static double MAX_ELEVATION = 1_500;
  final static int LENGTH = 111_000;
  final static double AZIMUTH = toRadians(27.97);
  final static double LONGITUDE = toRadians(6.15432);
  final static double LATITUDE = toRadians(46.20562);
  final static int WIDTH = 800, HEIGHT = 100;

  public static void main(String[] as) throws Exception {
    DiscreteElevationModel dDEM =
      new HgtDiscreteElevationModel(HGT_FILE);
    ContinuousElevationModel cDEM =
      new ContinuousElevationModel(dDEM);
    GeoPoint o =
      new GeoPoint(LONGITUDE, LATITUDE);
    ElevationProfile p =
      new ElevationProfile(cDEM, o, AZIMUTH, LENGTH);

    int BLACK = 0x00_00_00, WHITE = 0xFF_FF_FF;

    BufferedImage i =
      new BufferedImage(WIDTH, HEIGHT, TYPE_INT_RGB);
    for (int x = 0; x < WIDTH; ++x) {
      double pX = x * (double) LENGTH / (WIDTH - 1);
      double pY = p.elevationAt(pX);
      int yL = (int)((pY / MAX_ELEVATION) * (HEIGHT - 1));
      for (int y = 0; y < HEIGHT; ++y) {
        int color = y < yL ? BLACK : WHITE;
        i.setRGB(x, HEIGHT - 1 - y, color);
      }
    }
    dDEM.close();

    ImageIO.write(i, "png", new File("profile.png"));
  }
}

En l'exécutant, on obtient l'image ci-dessous, qui a été redimensionnée ici pour des questions de présentation mais est également disponible en taille réelle.

En comparant ce profil avec celui obtenu sur map.geo.admin.ch, présenté à la figure 3 plus haut, on constate que les deux sont généralement assez similaires, avec toutefois deux différences importantes :

• le nôtre s'arrête brusquement à environ 100 km, car le fichier HGT utilisé ne couvre pas la totalité de la zone couverte par le profil,
• leur profil est nettement plus lisse que le nôtre, car :
  1. leurs MNT sont plus précis que les nôtres, et
  2. ils appliquent probablement un filtre de lissage pour rendre leur profil plus doux.

Néanmoins, les ressemblances sont assez importantes pour nous donner confiance en notre programme.