Fusion, redimensionnement et masque de rasters dans QGIS et publication du GeoTiff dans GeoServer

de | | , ,

L’objectif de ce tutoriel est de masquer des rasters fusionnés en un fichier selon des limites territoriales définies puis de diffuser la donnée finale via un serveur cartographique.
Après le téléchargement des données SIG, on utilise des fonctions basiques pour traiter les rasters dans le logiciel libre QGIS. La donnée en sortie de format GeoTiff est ensuite stockée et stylisée dans GeoServer.

1. Rasters utilisés

Pour illustrer cet article, on utilise des données issues d’un modèle d’élévation numérique (DEM : digital elevation model) construit par l’ U.S. Geological Survey’s Center for Earth Resources Observation and Science (EROS) : le Global 30 Arc-Second Elevation (GTOPO30). Le DEM GTOPO30, construit avec des données topographiques vectorielles et satellites, a une résolution spatiale d’environ 1 km (par pixel). Le produit est divisé en 33 dalles dont 27 couvrent chacune 50° en latitude et 40° en longitude sur la zone d’étendue [-180°W ; -60°S ; 180°E ; 90°N] et 6 dalles d’une résolution de 30° en latitude et 60° en longitude sur l’étendue de l’Antarctique. Certes, la précision verticale de GTOPO30 n’est qu’en moyenne d’environ 30 m mais son principal avantage est d’avoir une couverture géographique mondiale comme son nom l’indique. A une échelle fine, il existe bien entendu d’autres données spatiales plus précises pour estimer l’altitude (par exemple les données ASTER, LIDAR, etc).

Ici, on ne cherche pas à construire un jeu de données exhaustifs et on utilise seulement quatre dalles : gt30e020n40,  gt30e020n90, gt30w020n40 et gt30w020n90. Les rasters GTOPO30 sont téléchargés sur le site de l’USGS (figure 1).

Figure 1 : Téléchargement des données GTOPO30 sur le site de l’USGS.

2. Traitements des rasters dans QGIS

A l’aide du logiciel QGIS et de la bibliothèque GDAL implémentée, trois traitements basiques des rasters sont mis en œuvre dans ce paragraphe : la fusion (2.1), le redimensionnement (2.2) et le masque de la donnée fusionnée (2.3). Je préfère considérer le terme de masque plutôt que de découpage utilisé dans QGIS car ce dernier terme porte à confusion avec le redimensionnement du raster. En effet, le redimensionnement est défini par la diminution ou l’augmentation du nombre de lignes et de colonnes du raster en fonction de limites géographiques prédéfinies sans modifier les valeurs et la taille des pixels.  On pourrait donc parler de découpage du raster selon une fenêtre géographique. Tandis que le masque du raster consiste à modifier des valeurs de quelques ou tous les pixels du raster en fonction d’un référentiel géographique ou de certaines valeurs du pixel.

On note que les étapes de redimensionnement et de masque de la donnée peuvent être inversées mais il est préférable de traiter une donnée raster avec des dimensions minimales afin d’abaisser l’utilisation de ressources informatiques.

2.1. Fusion des données

Le premier traitement consiste à assembler les dalles afin d’obtenir un fichier unique en sortie. Après le chargement des rasters dans QGIS, dans le menu Raster > Divers, on choisit l’onglet Fusionner… comme le montre la figure 2 ci-dessous :

Figure 2: Fusion de quatre rasters du DEM GTOPO30 dans QGIS.

Dans la fenêtre Fusionner, on sélectionne les fichiers rasters en input – ici situés dans un même répertoire – puis on saisit un nom en output avec l’extension désirée – ici, on choisit GeoTiff, format de fichier largement reconnu -. Les dalles du GTOPO30 ne se superposent pas et ne sont pas issus de capteurs différents. Il n’est donc pas nécessaire de sélectionner une bande séparée pour chaque fichier ou utiliser l’emprise intersectée. Pour le moment, le style du raster n’est pas considéré et on peut finalement lancer la procédure de fusion. Le raster en sortie de 9600 colonnes et de 12000 lignes s’étend sur la zone géographique [-20°W ; -10°S ; 60°E ; 90°N].

2.2. Redimensionnement du raster

Le raster fusionné de format GeoTiff couvre une surface considérable et pour des objectifs divers, on décide de cibler un espace spécifique. Il faut donc diminuer les dimensions du raster, c’est-à-dire qu’il faut déterminer le nombre de colonnes (X) et de lignes (Y) en fonction d’une fenêtre géographique définie.  Le nombre de colonnes (X) correspond à la différence entre la longitude maximale (Xmax) et la longitude minimale (Xmin) et le nombre de lignes (Y) est la différence entre la latitude maximale (Ymax) et la latitude minimale (Ymin). Pour le redimensionnement du raster, il est conseillé de garder la même taille de pixel avec des colonnes et des lignes entières. Par rapport à la donnée de base b, on a donc :

X = Xb / (Xbmax  – Xbmin ) x (Xmax – Xmin);

Y = Yb / (Ybmax  – Ybmin ) x (Ymax – Ymin);

Dans notre exemple, on choisit l’étendue géographique [-20°W ; 30°N ; 20°E ; 70°N] qui correspond approximativement à celle des pays délimités en rouge sur la figure 3, soit :

X = 9600 / (60-(-20)) x (20-(-20));

X = 4800 colonnes;

Y = 12000 / (90-(-10)) x (70-30);

Y = 4800 lignes;

Dans QGIS, on ouvre la Calculatrice Raster à partir du menu Raster puis on définit l’emprise géographique du nouveau raster. Dans l’expression, on ne modifie pas les valeurs du raster de base, on saisit donc le nom du fichier que l’on multiplie par un.

Figure 3 : Redimensionnement du raster à l’aide de la Calculatrice Raster.

Vous l’aurez compris ou vous le savez déjà, de nombreuses opérations de masque avec l’altération des valeurs de pixels sont possibles avec l’outil de Calculatrice Raster. Toutefois, on utilise un autre outil de masque dans le logiciel.

2.3. Masque du raster fusionné et redimensionné

Dans QGIS, dans le menu Raster > Extraction > Découper, deux méthodes de masque sont possibles :

  • la première se base sur des limites définies par nos soins ou à l’aide d’un rectangle dessiné sur le canvas;
  • l’autre découpage suit les contours du vecteur référentiel. Pour cette méthode que l’on choisit dans ce tutoriel, il faut faire attention que le vecteur soit bien dessiné pour ne pas « planter » le traitement.

Dans la fenêtre, on choisit la donnée raster à masquer puis on saisit le nom du fichier de sortie. Il faut bien comprendre que le raster créé aura les mêmes dimensions avec des nombres de colonnes et de lignes identiques au fichier d’entrée. Ici, on utilise un vecteur représentant quelques pays délimités en rouge sur la figure 4 comme couche de masquage. C’est-à-dire que seuls les pixels du raster en entrée inclus dans le périmètre du vecteur référentiel sont pris en compte. On affecte une valeur -9999 aux pixels externes aux contours des pays considérés et on sélectionne « Découper l’emprise du jeu de données cible selon l’emprise du trait de coupe ». La résolution de la taille des pixels est conservée et on peut enfin lancer le traitement.

Figure 4 : Masque du raster à partir des contours d’un vecteur dans QGIS.

La figure 5 montre le résultat du masque du DEM GTOPO30 préalablement fusionné et redimensionné. On lui affecte un style de type « Pseudo-couleur à bande unique » avec une échelle d’altitude définie pour l’exemple. Une couleur est affectée à chaque intervalle de valeurs et pour un meilleur rendu visuel, on choisit l’interpolation linéaire qui étale les couleurs selon les valeurs réelles du pixel. Une interpolation de type discret implique que tous les pixels d’un intervalle se voient attribuer une même couleur; ce qui a pour effet un rendu moins lisse.

Figure 5 : Création d’un style par interpolation linéaire pour le raster dans QGIS.

Enfin, le style est enregistré dans un fichier de format QML détaillé dans la paragraphe 3.2. Il est à présent temps de diffuser la donnée grâce au serveur cartographique GeoServer.

3. Publication du GeoTiff dans GeoServer

Le raster créé dans la section 2 est publié puis stylisé.

3.1. Création de l’espace de travail et de l’entrepôt

Dans GeoServer, il existe un type d’entrepôt spécifique pour les données GTOPO30 avec l’extension DEM moins volumineuse que le GeoTiff. Toutefois, cette dernière extension est le format de fichier utilisé dans notre tutoriel qui est aussi une source de donnée image gérée par le serveur cartographique. Comme pour le tutoriel sur la création d’un fond cartographique tuilé, on créé un espace de travail et un entrepôt pour une source de données images de type GeoTiff (figure 6). Ce type d’entrepôt prend en compte qu’un seul raster, on verra prochainement qu’il est possible de charger plusieurs rasters dans d’autres types d’entrepôts comme la source ImageMosaic. Dans cet exemple, la donnée est dans le répertoire de données par défaut de GeoServer data_dir, en production il est conseillé de stocker les données SIG dans un répertoire externe à GeoServer.

Figure 6 : Création de la source de donnée GeoTiff dans GeoServer.

Ensuite, on peut publier la couche SIG en spécifiant au minimum l’emprise géographique et son nom.

3.2. Style du raster dans GeoServer

Dans ce paragraphe, on cherche à attribuer le style exporté depuis QGIS au raster. Toutefois, même si certaines améliorations de compatibilité de styles entre QGIS et GeoServer ont été réalisées, le fichier de format QML décrit ci-dessous n’est pas pris en charge par GeoServer.

<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="2.18.9" minimumScale="inf" maximumScale="1e+08" hasScaleBasedVisibilityFlag="0">
  <pipe>
    <rasterrenderer opacity="1" alphaBand="0" classificationMax="1380" classificationMinMaxOrigin="CumulativeCutFullExtentEstimated" band="1" classificationMin="2" type="singlebandpseudocolor">
      <rasterTransparency/>
      <rastershader>
        <colorrampshader colorRampType="INTERPOLATED" clip="0">
          <item alpha="255" value="-10" label="-9999 - -10" color="#fffffd"/>
          <item alpha="255" value="0" label="-10 - 0" color="#c5ffea"/>
          <item alpha="255" value="50" label="0 - 50" color="#ffe98c"/>
          <item alpha="255" value="100" label="50- 100" color="#ffd265"/>
          <item alpha="255" value="150" label="100 - 150" color="#feb751"/>
          <item alpha="255" value="300" label="150 -300" color="#fe9b43"/>
          <item alpha="255" value="600" label="300 - 600" color="#fb7b35"/>
          <item alpha="255" value="1200" label="600 - 1200" color="#f55629"/>
          <item alpha="255" value="1800" label="1200 - 1800" color="#eb3420"/>
          <item alpha="255" value="2400" label="1800 - 2400" color="#d41a23"/>
          <item alpha="255" value="3200" label="2400 - 3200" color="#970000"/>
          <item alpha="255" value="5000" label="3200 - 5000" color="#51000e"/>
        </colorrampshader>
      </rastershader>
    </rasterrenderer>
    <brightnesscontrast brightness="0" contrast="0"/>
    <huesaturation colorizeGreen="128" colorizeOn="0" colorizeRed="255" colorizeBlue="128" grayscaleMode="0" saturation="0" colorizeStrength="100"/>
    <rasterresampler maxOversampling="2"/>
  </pipe>
  <blendMode>0</blendMode>
</qgis>

Il faut donc le modifier et l’adapter en format SLD. La documentation de GeoServer sur les styles des rasters ainsi que des exemples sont d’une aide précieuse.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
  <NamedLayer>
    <Name>Gtopo30</Name>
    <UserStyle>
      <Name>Gtopo30</Name>
      <Title>Gtopo30 distribution</Title>
      <FeatureTypeStyle>
        <Rule>
          <RasterSymbolizer>
            <Opacity>1.0</Opacity>
            <ColorMap>
              <ColorMapEntry color="#ffffff" quantity="-10" label="-9999 - -10"/>
              <ColorMapEntry color="#c5ffea" quantity="0" label="-10 - 0"/>
              <ColorMapEntry color="#ffe98c" quantity="50" label="0 - 50"/>
              <ColorMapEntry color="#ffd265" quantity="100" label="50 - 100"/>
              <ColorMapEntry color="#feb751" quantity="150" label="100 - 150"/>
              <ColorMapEntry color="#fe9b43" quantity="300" label="150 - 300"/>
              <ColorMapEntry color="#fb7b35" quantity="600" label="300 - 600"/>
              <ColorMapEntry color="#f55629" quantity="1200" label="600 - 1200"/>
              <ColorMapEntry color="#eb3420" quantity="1800" label="1200 - 1800"/>
              <ColorMapEntry color="#d41a23" quantity="2400" label="1800 - 2400"/>
              <ColorMapEntry color="#970000" quantity="3200" label="2400 - 3200"/>
              <ColorMapEntry color="#51000e" quantity="5000" label="3200 - 5000"/>
            </ColorMap>
          </RasterSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

On note que les structures des fichiers sont relativement proches. Pour le fichier SLD, on conserve l’interpolation linéaire des couleurs en lien avec les classes de valeurs. Le style est validé puis affecté au raster que l’on peut prévisualiser comme le montre la figure 7 ci-après :

Figure 7 : Style et prévisualisation du raster représentant l’élévation dans GeoServer.

La couche SIG est désormais diffusable via les protocoles WM(T)S.

En conclusion, ce tutoriel nous a permis d’introduire trois traitements basiques réalisés sur des rasters issus du modèle d’élévation numérique GTOPO30. La fusion, le redimensionnement puis le masque des rasters utilisés sont des opérations relativement simples dans QGIS. Pour la manipulation de nombreux fichiers, il est préférable d’utiliser directement GDAL. Le raster créé de format GeoTiff a été publié dans GeoServer puis stylisé en adaptant le fichier QML exporté depuis QGIS. En perspective, les données d’élévation issues du fichier peuvent être utilisées en webmapping.

Partager l'article
Taggé , , , .Mettre en favori le Permaliens.

A propos Florian Delahaye

Passionné de Géomatique

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *

  • ☰ Tutoriels SIG les plus populaires