Modèles numériques du terrain
Etape 3

1 Introduction

Le but de cette troisième étape est d'écrire des classes et interfaces représentant ce que l'on nomme des modèles numériques du terrain, souvent abrégé MNT — ou DEM en anglais, pour digital elevation model.

Un modèle numérique du terrain permet de connaître l'altitude d'un point à la surface de la Terre. Il s'agit donc d'une fonction à deux variables — les coordonnées d'un point — dont la valeur est l'altitude en ce point.

Bien entendu, la surface de la Terre est une fonction continue mais, comme souvent en informatique, on la transforme en fonction discrète par échantillonnage, afin de lui donner une représentation finie. Cette fonction discrète est ensuite interpolée pour obtenir à nouveau une fonction continue. Dès lors, nous distinguerons deux types de MNT : les MNT discrets et les MNT continus.

Dans le cadre de ce projet, les MNT discrets proviendront de fichiers — qui vous seront fournis à l'étape suivante — tandis que les MNT continus seront obtenus par interpolation bilinéaire des MNT discrets.

1.1 Modèle discret du terrain

Un MNT discret donne l'altitude de la surface terrestre pour un nombre fini de points régulièrement espacés.

Les altitudes données par un MNT sont appelées échantillons (samples en anglais). Pour faciliter les choses, ces échantillons sont indexés par des coordonnées entières qui diffèrent des coordonnées géographiques du point correspondant. Mathématiquement, un MNT discret est donc une fonction \(f: \Bbb{Z}\times\Bbb{Z}\rightarrow\Bbb{R}\).

Les échantillons des MNT discrets utilisés dans ce projet sont espacés d'une seconde d'arc, soit 1°/3600, en longitude comme en latitude. Notez que même si cet espacement est régulier en degrés sur toute la surface terrestre, il varie en mètres car la longueur d'un parallèle — un cercle de latitude constante — varie en fonction de la latitude. Ainsi, à l'équateur, le parallèle a une longueur d'environ 40000 km, tandis qu'aux pôles il a une longueur nulle.

Par exemple, les Alpes étant à une latitude proche de 45°, une seconde d'arc correspond environ à 30 m dans la direction nord-sud et 20 m dans la direction est-ouest :

\begin{align*} d_{\text{NS}} &= \frac{2\pi\cdot 6\,371\,000}{360\cdot 3\,600} \approx 30.89\\ d_{\text{EO}} &= \frac{2\pi\cdot 6\,371\,000}{360\cdot 3\,600}\,\cos 45° \approx 21.84 \end{align*}

Ces distances sont celles séparant, sous nos latitudes, deux échantillons des MNT que nous utiliserons. Elles donnent une idée de la qualité des panoramas qu'il nous sera possible d'obtenir.

1.1.1 Indexation

Les échantillons d'un MNT discret sont indexés par une paire d'entiers dont la première composante représente la longitude, et la seconde la latitude. De plus, l'échantillon d'index (0,0) correspond au point de coordonnées (0°,0°). En conséquence, les composantes d'un index peuvent être négatives. C'est le cas p.ex. pour tous les points se trouvant dans l'hemisphère sud ou à l'ouest du méridien de référence.

Etant donné que deux échantillons sont toujours séparés d'une seconde d'arc dans chacune des dimensions, les coordonnées \((\lambda, \varphi)\) du point correspondant à l'échantillon d'index \((x,y)\) sont simplement données par :

\[ (\lambda, \varphi) = \left(\frac{x}{3\,600}, \frac{y}{3\,600}\right) \]

Attention : ces coordonnées sont en degrés !

Par exemple, l'échantillon correspondant au point de coordonnées (6°34′04″, 46°31′07″) a l'index :

\begin{align*} (x, y) &= ([[6\cdot 60] + 34]\cdot 60 + 4, [[46\cdot 60] + 31]\cdot 60 + 7)\\ &= (23\,644, 167\,467) \end{align*}

1.1.2 Etendue

Un MNT discret donné ne couvre généralement qu'une partie de la surface de la Terre, et non pas sa totalité. Dès lors, il est utile de pouvoir connaître ce que nous appellerons l'étendue d'un MNT (extent en anglais), c'est-à-dire les index de tous les échantillons qu'il contient.

Notez que, en termes mathématiques, l'étendue n'est rien d'autre que le domaine de définition de la fonction \(f\) que constitue le MNT, que nous noterons \(D_f\).

Pour faciliter les choses, nous ne manipulerons dans ce projet que des MNT discrets dont l'étendue est représentable par un intervalle bidimensionnel. Ainsi, un MNT discret couvrant la zone comprise entre 6° et 7° de longitude et 45° et 47° de latitude aura pour étendue l'intervalle bidimensionnel :

\[\newcommand{\iint}[2]{[#1\,{.\!.}\,#2]} \iint{6\cdot 3\,600}{7\cdot 3\,600}\times\iint{45\cdot 3\,600}{47\cdot 3\,600}\\ \ \ = \iint{21\,600}{25\,200}\times\iint{162\,000}{169\,200} \]

1.1.3 Extension

Dans certaines situations, le fait qu'un MNT discret ne couvre qu'une partie de la surface de la Terre est peu pratique, et il serait plus simple d'avoir un MNT la couvrant en totalité.

Pour cette raison, nous définirons l'extension d'un MNT discret donné comme un MNT discret couvrant la totalité de la surface de la Terre, égal au MNT qu'il étend partout où celui-ci est défini, et nul ailleurs.

Mathématiquement, l'extension d'un MNT \(f\), notée \(\uparrow f\), est définie ainsi :

\[ (\uparrow f)(x) = \begin{cases} f(x) & \text{si } x\in D_f\\ 0 & \text{sinon} \end{cases} \]

1.1.4 Union

Pour peu que leurs étendues soient unionables, il est tout à fait possible de combiner deux MNT discrets pour en obtenir un nouveau dont l'étendue est l'union de celle des MNT combinés.

Mathématiquement, cette union de deux MNT discrets \(f_1\) et \(f_2\), notée par \(\oplus\), est définie ainsi :

\[ (f_1\oplus f_2)(x) = \begin{cases} f_1(x) & \text{si } x\in D_{f_1}\\ f_2(x) & \text{si } x\in D_{f_2} \end{cases} \]

Cette opération est très utile pour obtenir des MNT couvrant une zone importante — p.ex. la totalité des Alpes — à partir de MNT plus petits. Comme nous le verrons à l'étape suivante, les MNT discrets que nous utiliserons proviennent de fichiers couvrant chacun une zone de 1° de côté, en latitude et en longitude, et un seul d'entre eux ne suffit que rarement au dessin d'un panorama.

1.2 Modèle continu du terrain

Un MNT continu donne l'altitude de la surface terrestre pour un point quelconque, dont la position est exprimée en coordonnées sphériques. Mathématiquement, un MNT continu est donc une fonction \(f: \Bbb{R}\times\Bbb{R}\rightarrow\Bbb{R}\).

Tous les MNT continus manipulés dans le cadre de ce projet sont obtenus par interpolation bilinéaire de l'extension d'un MNT discret. Il en découle qu'un MNT continu couvre toujours la totalité de la surface terrestre, et aucune étendue ne lui est donc associée.

En plus de l'altitude, obtenue par simple interpolation, le dessin de panorama requiert la connaissance de la pente du terrain en un point donné. Son calcul n'est pas trivial et mérite d'être décrit en détail.

1.2.1 Pente

La pente est l'angle formé entre un vecteur normal au terrain et un vecteur vertical.

En admettant que l'on connaisse un vecteur normal au terrain \(\vec{n}\), on peut calculer la pente \(\theta\) au moyen du produit scalaire, sachant que :

\[ \vec{n}\cdot\vec{v} = \|\vec{n}\|\,\|\vec{v}\|\,\cos(\theta) \]

où \(\vec{v}\) est un vecteur unitaire vertical, \(\vec{v} = (0\ 0\ 1)^T\) et donc \(\|\vec{v}\| = 1\). On en déduit la formule suivante pour déterminer la pente :

\[\theta = \arccos\left(\frac{n_z}{\|\vec{n}\|}\right)\]

où \(n_z\) est la composante verticale de \(\vec{n}\).

Il reste donc à déterminer le vecteur \(\vec{n}\). Pour un point du MNT discret d'index \((x,y)\), on peut le calculer au moyen du produit vectoriel de deux vecteurs \(\vec{a}\) et \(\vec{b}\) définis à partir de deux points voisins du MNT, comme l'illustre la figure ci-dessous.

Sorry, your browser does not support SVG.

Figure 1 : Vecteur normal au terrain

En effet, \(\vec{n} = \vec{a}\wedge\vec{b}\). En faisant l'hypothèse que la distance séparant deux échantillons est la même dans les deux directions — ce qui est faux, comme nous l'avons vu plus haut, mais suffit à nos besoins — et vaut \(d\), les vecteurs \(\vec{a}\) et \(\vec{b}\) sont donnés par :

\[ \vec{a} = \begin{pmatrix} d\\ 0\\ \Delta z_a \end{pmatrix}, \ \ \vec{b} = \begin{pmatrix} 0\\ d\\ \Delta z_b \end{pmatrix} \]

où \(\Delta z_a\) est la différence d'altitude entre le point d'index \((x + 1, y)\) et celui d'index \((x, y)\), et \(\Delta z_b\) est la différence d'altitude entre le point d'index \((x, y + 1)\) et celui d'index \((x,y)\).

En substituant \(\vec{n}\) par \(\vec{a}\wedge\vec{b}\) dans la formule de calcul de pente présentée plus haut puis en simplifiant, on obtient finalement la formule suivante : \[ \theta = \arccos\,\left(\frac{d}{\sqrt{\Delta z_a^2 + \Delta z_b^2 + d^2}}\right) \]

qui donne la pente \(\theta\) du point du MNT discret d'index \((x, y)\). Pour déterminer la pente d'un point quelconque du MNT continu, on utilise simplement cette formule quatre fois pour déterminer la pente des quatres points voisins du MNT discret, puis on effectue une interpolation bilinéaire entre les valeurs obtenues.

2 Mise en œuvre Java

La totalité des classes de cette étape se trouve dans le paquetage ch.epfl.alpano.dem, réservé au code lié aux modèles numériques du terrain.

2.1 Interface DiscreteElevationModel

L'interface DiscreteElevationModel représente un MNT discret. Pour des raisons qui deviendront claires plus tard, elle étend l'interface AutoCloseable du paquetage java.lang. Elle offre les attributs statiques suivants :

  • int SAMPLES_PER_DEGREE, qui contient le nombre d'échantillons par degré d'un MNT discret, soit 3600,
  • double SAMPLES_PER_RADIAN, qui contient le nombre d'échantillons par radian d'un MNT discret.

En plus de ces attributs, elle offre une unique méthode statique :

  • double sampleIndex(double angle), qui retourne l'index correspondant à l'angle donné, exprimé en radians ; cet angle peut représenter soit une longitude, soit une latitude, étant donné que la distance angulaire entre deux échantillons est égale dans les deux cas.

De plus, l'interface DiscreteElevationModel possède deux méthodes abstraites :

  • Interval2D extent(), qui retourne l'étendue du MNT,
  • double elevationSample(int x, int y), qui retourne l'échantillon d'altitude à l'index donné, en mètres, ou lève l'exception IllegalArgumentException si l'index ne fait pas partie de l'étendue du MNT.

Finalement, l'interface DiscreteElevationModel offre une méthode par défaut dont le but est de faciliter la mise en commun de deux modèles discrets :

  • default DiscreteElevationModel union(DiscreteElevationModel that), qui retourne un MNT discret représentant l'union du récepteur et de l'argument that, ou lève l'exception IllegalArgumentException si leurs étendues ne sont pas unionables.

Notez que la méthode union ne fait rien d'autre que créer et retourner une instance de la classe CompositeDiscreteElevationModel, décrite ci-après.

2.2 Classe CompositeDiscreteElevationModel

La classe CompositeDiscreteElevationModel, visible uniquement dans son paquetage et immuable (donc finale), représente l'union de deux modèles du terrain discrets. Logiquement, cette classe implémente l'interface DiscreteElevationModel. Son constructeur prend en arguments les deux modèles du terrain à composer :

  • CompositeDiscreteElevationModel(DiscreteElevationModel dem1, DiscreteElevationModel dem2), qui retourne un MNT discret représentant l'union des MNT donnés, ou lève l'exception NullPointerException si l'un de ces MNT est nul.

En dehors de ce constructeur, cette classe ne contient rien d'autre que des mises en œuvre des méthodes abstraites des interfaces qu'elle implémente, à savoir DiscreteElevationModel et AutoCloseable.

2.2.1 Interface AutoCloseable

L'interface AutoCloseable ne contient qu'une seule méthode, close. Le but de cette méthode est de libérer les resources associées à l'objet auquel on l'applique. Par exemple, lorsqu'un tel objet représente un fichier ouvert, la méthode close le ferme et libère ainsi la mémoire qui lui est réservée par le système d'exploitation. Cette notion de resource sera examinée plus en détails dans le cours sur les entrées-sorties.

Dans le cas qui nous intéresse ici, à savoir la méthode close de la classe CompositeDiscreteElevationModel, la seule chose à faire est d'appeler les méthodes close de chacun des deux modèles du terrain qu'elle compose. Cela est logique, car pour libérer les resources associées à une composition de deux objets, il faut libérer les resources associées à chacun d'eux.

2.2.2 Vérification de nullité

Le constructeur de CompositeDiscreteElevationModel doit lever l'exception NullPointerException si l'un de ses deux arguments est nul.

Dans la plupart des cas, il n'est pas nécessaire de lever cette exception explicitement, car elle l'est automatiquement lorsqu'on essaie d'accéder à un membre d'un objet nul. Toutefois, dans le cas de ce constructeur, cette exception ne serait pas levée automatiquement, car la seule chose que le constructeur fait est de stocker les MNT dans des champs.

Etant donné qu'il est toujours souhaitable de signaler les erreurs de programmation aussi rapidement que possible, car cela facilite leur diagnostic, le constructeur de CompositeDiscreteElevationModel doit vérifier explicitement que ses arguments ne sont pas nuls avant de les stocker. Cela peut se faire très facilement au moyen de la méthode requireNonNull de la classe Objects (avec un s, qui provient du paquetage java.util), qui retourne son argument si celui-ci n'est pas nul, et lève NullPointerException sinon. On peut donc l'utiliser ainsi dans un constructeur : 

import static java.util.Objects.requireNonNull;

public final class MyClass {
  private final Object someObject;

  public MyClass(Object someObject) {
    this.someObject = requireNonNull(someObject);
  }

  // … reste de la classe
}

Notez bien que cette méthode n'est destinée à être utilisée que dans des cas similaires ! Il ne faut pas l'utiliser à tout bout de champ dans le code, car dans la quasi-totalité des cas, elle est redondante. Par exemple, si le constructeur de MyClass ci-dessus accédait à un membre de someObject, l'exception NullPointerException serait levée à ce moment-là, et l'appel à requireNonNull pourrait être supprimé.

2.3 Classe ContinuousElevationModel

La classe ContinuousElevationModel, publique et immuable (donc finale), représente un MNT continu, obtenu par interpolation d'un MNT discret. Elle possède le constructeur suivant :

  • ContinuousElevationModel(DiscreteElevationModel dem), qui construit un MNT continu basé sur le MNT discret passé en argument, ou lève l'exception NullPointerException si celui-ci est nul.

En plus de ce constructeur, cette classe offre les méthodes publiques suivantes :

  • double elevationAt(GeoPoint p), qui retourne l'altitude au point donné, en mètres, obtenue par interpolation bilinéaire de l'extension du MNT discret passé au constructeur,
  • double slopeAt(GeoPoint p), qui retourne la pente du terrain au point donné, en radians.

2.3.1 Méthodes privées

Vous vous faciliterez beaucoup la tâche en définissant quelques méthodes auxiliaires privées pour cette classe. Vous êtes libres de choisir lesquelles définir, mais nous vous suggérons d'en définir au moins deux :

  • une première retournant l'altitude d'un point de l'extension du DEM discret passé au constructeur,
  • une seconde retournant la pente d'un point de l'extension du DEM discret passé au constructeur.

2.3.2 Distance entre échantillons

Pour la distance entre deux échantillons, notée \(d\) dans les formules ci-dessus, utilisez la distance nord-sud, notée \(d_{NS}\) plus haut et valant approximativement 30.89 m.

N'entrez toutefois pas directement cette valeur dans vos formules, mais calculez-la au moyen de la constante SAMPLES_PER_RADIAN et de la méthode toMeters de l'interface Distance. Stockez le résultat de ce calcul dans un champ privé statique de la classe, pour éviter de le faire plusieurs fois.

2.4 Tests

Comme pour l'étape précédente, nous ne vous fournissons plus de tests mais un fichier de vérification de noms contenu dans une archive Zip à importer dans votre projet.

Pour vérifier que votre code fonctionne, nous vous conseillons bien entendu d'écrire des tests unitaires, mais aussi de générer des images basées sur des MNT fictifs. Pour vous aider, nous vous fournissons tout d'abord la classe WavyDEM ci-dessous, qui modélise un MNT discret obtenu par échantillonnage d'une fonction ayant la même forme que celle-ci :

\[f(x, y) = \sin x\,\cos y\]

dont le graphe tri-dimensionnel est :

sinx_cosy.png

Figure 2 : Graphe de \(\sin x\,\cos y\)

La classe WavyDEM ajuste cette fonction pour que les altitudes soient comprises entre 0 et 1000 m, et pour que sa période soit de 100 échantillons. 

final class WavyDEM implements DiscreteElevationModel {
  private final static double PERIOD = 100, HEIGHT = 1000;
  private final Interval2D extent;

  public WavyDEM(Interval2D extent) {
    this.extent = extent;
  }

  @Override
  public void close() throws Exception { }

  @Override
  public Interval2D extent() { return extent; }

  @Override
  public double elevationSample(int x, int y) {
    double x1 = PI * 2d * x / PERIOD;
    double y1 = PI * 2d * y / PERIOD;
    return (1 + sin(x1) * cos(y1)) / 2d * HEIGHT;
  }
}

Pour générer des images à partir d'un MNT continu basé sur WavyDEM, vous pouvez utiliser la classe DrawDEM ci-dessous. Cette classe construit deux MNT de type WavyDEM, en fait l'union puis construit un MNT continu de l'extension de cette union. Ensuite, elle génère deux images :

  1. la première, stockée dans le fichier elevation.png, représente l'altitude du MNT au moyen de niveaux de gris, où noir correspond à 0 m et blanc à 1000 m,
  2. la seconde, stockée dans le fichier slope.png, représente la pente du MNT au moyen de niveaux de gris, où noir correspond à 0° et blanc à 90°.

La zone représentée sur ces deux images correspond, pour le MNT discret, à l'intervalle bidimensionnel [0..100]×[0..100]. Comme les images font 300 pixels de large, l'interpolation entre en jeu. 

import static java.awt.image.BufferedImage.TYPE_INT_RGB;
import static java.lang.Math.*;

import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;

import javax.imageio.ImageIO;

public final class DrawDEM {
  @SuppressWarnings("resource")
  public static void main(String[] args)
    throws IOException {

    DiscreteElevationModel dDEM1 =
      new WavyDEM(new Interval2D(new Interval1D(0, 50),
                                 new Interval1D(0, 100)));
    DiscreteElevationModel dDEM2 =
      new WavyDEM(new Interval2D(new Interval1D(50, 100),
                                 new Interval1D(0, 100)));
    DiscreteElevationModel dDEM =
      dDEM1.union(dDEM2);
    ContinuousElevationModel cDEM =
      new ContinuousElevationModel(dDEM);

    int size = 300;
    double scale = (100d / 3600d) / (size - 1);
    BufferedImage elI =
      new BufferedImage(size, size, TYPE_INT_RGB);
    BufferedImage slI =
      new BufferedImage(size, size, TYPE_INT_RGB);
    for (int x = 0; x < size; ++x) {
      for (int y = 0; y < size; ++y) {
        GeoPoint p = new GeoPoint(toRadians(x * scale),
                                  toRadians(y * scale));
        double el = cDEM.elevationAt(p);
        elI.setRGB(x, y, gray(el / 1000d));

        double sl = cDEM.slopeAt(p);
        slI.setRGB(x, y, gray(sl / (PI / 2d)));
      }
    }

    ImageIO.write(elI, "png", new File("elevation.png"));
    ImageIO.write(slI, "png", new File("slope.png"));
  }

  private static int gray(double v) {
    double clampedV = max(0, min(v, 1));
    int gray = (int) (255.9999 * clampedV);
    return (gray << 16) | (gray << 8) | gray;
  }
}

Après exécution de ce programme, les fichiers susmentionnés se trouveront dans le dossier de votre projet, et en les visualisant vous devriez constater qu'ils ont l'apparence suivante :

elevation.png

Figure 3 : Dessin de l'altitude (elevation.png)

slope.png

Figure 4 : Dessin de la pente (slope.png)

3 Résumé

Pour cette étape, vous devez :

  • écrire :

    • l'interface DiscreteElevationModel,
    • la classe CompositeDiscreteElevationModel, et
    • la classe ContinuousElevationModel

    en fonction des indications données plus haut,

  • documenter la totalité des entités publiques que vous avez définies,
  • rendre votre code au plus tard le 10 mars 2017 à 16h30, via le système de rendu.

N'attendez surtout pas le dernier moment pour effectuer votre rendu, car vous n'êtes pas à l'abri d'imprévus et aucun retard, aussi insignifiant soit-il, ne sera toléré !