Shapefiles : avalanches et bouquetins · Blog

Chez Liip, nous venons de mettre au point une application destinée aux étudiant·e·s en architecture. Celle-ci utilise les shapefiles comme données de base pour fournir des informations sur les terrains constructibles dans la ville de Zurich. Ils apparaissent sur une carte et peuvent être consultés par l’utilisateur·trice. Un projet similaire a été mené par une équipe d’étudiant·e·s, dont je faisais partie, pour la filière FHNW HT: une application web ludique qui permet aux enfants de l’école primaire d’apprendre à connaître leur commune.

Non seulement les shapefiles contribuent à améliorer l’expérience utilisateur·trice, mais ils contiennent également des données qui, lorsqu’elles sont mises en contexte, peuvent fournir des informations intéressantes. Regardons ce que je peux trouver en fouillant dans quelques shapefiles et en créant une petite application pour les explorer.
Dans ce billet de blog, j’adopte une approche pratique de l’utilisation des shapefiles et je vous montre comment vous pouvez les utiliser.

La Convention alpine

Pour utiliser des shapefiles, il faut d’abord comprendre de quoi il s’agit exactement et pouvoir s’en procurer.

Les shapefiles sont des fichiers binaires qui contiennent des données géographiques. Ils se présentent généralement par groupes d’au moins 3 fichiers différents qui portent tous le même nom. Ceux-ci se différencient seulement par leur extension: .shp, .shx et .dbf. Cependant, un ensemble complet de données peut contenir beaucoup plus de fichiers. Ce format de données a été mis au point par Esri et introduit au début des années 1990. Les shapefiles contiennent des « caractéristiques », c’est-à-dire des formes géométriques avec leurs métadonnées. Une caractéristique peut être un point (coordonnées X/Y uniques) et un texte correspondant à ce point. Ou un polygone, constitué de plusieurs points X/Y, d’une ligne simple, d’une ligne multiple, etc.

Pour obtenir des shapefiles, je consulte le site OpenData.swiss. OpenData propose 102 ensembles de données de shapefiles différents ( novembre 2017). Il me suffit de les télécharger pour pouvoir les utiliser. Pour cet exemple, j’ai choisi un shapefile plutôt petit : Alpine convention, qui comprend les périmètres de la Convention alpine en Suisse.

Ces fichiers étant binaires et ne pouvant pas être consultés dans un éditeur de texte, j’ai besoin d’un outil pour consulter le fichier que je viens de télécharger.

Présentation de QGIS, un système d’information géographique gratuit et open source

QGIS est une solution complète pour tous les types de SIG (Systèmes d’information géographique). Il peut être téléchargé et installé gratuitement, et dispose de nombreux plugins et d’une communauté active. Je vais l’utiliser pour regarder le shapefile que je viens de télécharger.
Pour cet exemple, j’ai installé QGIS avec le plugin OpenLayers. Voici à quoi il ressemble, une fois lancé:

Pour obtenir des références sur l’endroit où je me trouve réellement et où vont mes données, j’ajoute OpenStreetMap comme calque. Je fais également un panoramique et un zoom sur la carte pour afficher la Suisse en gros plan.

L’étape suivante est d’ouvrir réellement le shapefile afin que QGIS puisse l’afficher sur la carte. Pour cela, il me suffit de double-cliquer sur le fichier .shp dans le navigateur à gauche.

Cette vue me donne déjà quelques informations sur le shapefile que je traite, sur la zone qu’il couvre et sur le nombre de caractéristiques qu’il contient. Le shapefile Alpine convention semble ne comporter qu’une seule caractéristique, un grand polygone, ce qui est suffisant, étant donné qu’il ne montre que les périmètres de la Convention alpine en Suisse.
Pour afficher les zones et/ou les détails qu’il couvre, je peux modifier le style du calque, en le rendant transparent. Je fais également un zoom pour voir si la forme est précise. Ses bords doivent recouvrir exactement les frontières de la Suisse.

Parfait. Maintenant que je vois la forme de la caractéristique, je vais regarder les métadonnées contenues dans le shapefile. Elles peuvent être consultées lorsque l’on ouvre le tableau des attributs du shapefile dans le navigateur situé en bas à gauche.

Ce fichier ne contient pas non plus beaucoup de métadonnées, mais je peux maintenant voir qu’il contient en fait deux caractéristiques. Il pourrait y avoir des shapefiles plus intéressants à utiliser, mais je vais m’en tenir au sujet alpin.

Des avalanches et des bouquetins

En fouillant un peu plus dans le répertoire de shapefiles d’OpenData, j’ai trouvé deux autres shapefiles qui pourraient fournir des données intéressantes: Distribution of ibex colonies (la répartition des colonies de bouquetins) et State of natural hazard mapping (Communes): Avalanches (la carte des zones avalancheuses).

Même si QGIS est un outil formidable pour examiner les données disponibles, je veux créer mon propre explorateur, peut-être même à partir de plusieurs shapefiles et selon ma propre conception. Pour cela, il existe de multiples bibliothèques pour tous types de langages. Pour le développement web, les trois bibliothèques les plus intéressantes seraient: PHP, Python et JavaScript.

Pour mon exemple, je vais créer une petite application front-end en JavaScript pour afficher mes trois shapefiles : la Convention alpine, les colonies de bouquetins et la cartographie des terrains avalancheux. Pour cela, je vais utiliser un package JavaScript tout simplement appelé «shapefile» ainsi que leafletjs et afficher les polygones à partir des caractéristiques. Mais chaque chose en son temps.

Lecture des caractéristiques

Avis de non-responsabilité : les exemples de code suivants utilisent ES6 et impliquent une configuration webpack/babel. Ils ne peuvent pas être copiés/collés, mais ils montrent comment travailler avec les outils disponibles.
Tout d’abord, j’essaie de charger le shapefile de la convention alpine et de l’enregistrer dans la console. Le package shapefile que j’ai mentionné est accompagné de deux exemples dans le README. Donc, pour simplifier, je me contente de cette approche :

import { open } from 'shapefile'

open('assets/shapefiles/alpine-convention/shapefile.shp').then(source => {
  source.read().then(function apply (result) { // Read feature from the shapefile
    if (result.done) { // If there's no features left, abort
      return
    }

    console.log(result)

    return source.read().then(apply) // Iterate over result
  })
})

Cela fonctionne à merveille ! On obtient les deux caractéristiques du shapefile dans la console. Ce que je vois déjà ici, c’est que les propriétés de chaque caractéristique, normalement stockées dans des fichiers séparés, sont déjà incluses dans la caractéristique. C’est le résultat de la bibliothèque qui récupère les deux fichiers et les traite ensemble:

Maintenant, je regarde les coordonnées qui sont liées à la première caractéristique. La première chose qui se produit est l’affichage de deux groupes différents, ce qui implique deux ensembles de coordonnées distincts, mais je vais ignorer le second pour le moment.

Et c’est là le premier piège. Ces coordonnées ne ressemblent en rien au WGS84 (latitude/longitude), mais plutôt à quelque chose d’autre. Les coordonnées en WGS84 devraient être un nombre à virgule flottante, commençant idéalement par 47 et 8 pour la Suisse. Donc, ce shapefile m’a confronté à un système de référence de coordonnées différent (ou CRS). Afin d’afficher correctement la forme en haut d’une carte, ou de l’utiliser avec d’autres données, je vais d’abord convertir ces coordonnées en WGS84, mais pour ce faire, je dois déterminer quel CRS ce shapefile utilise. Puisque ce shapefile provient de l’Office fédéral du développement territorial ARE, le CRS utilisé est très probablement CH1903, alias le système de coordonnées géographiques suisse.

Donc, pour faire la conversion, j’ai besoin de faire des maths. En cherchant un peu sur le web, j’ai trouvé une solution JavaScript pour calculer les allers-retours entre CH1903 et WGS84. Mais je n’ai besoin que de certaines parties, alors je copie et modifie un peu le code :

// Inspired by https://raw.githubusercontent.com/ValentinMinder/Swisstopo-WGS84-LV03/master/scripts/js/wgs84_ch1903.js

/**
 * Converts CH1903(+) to Latitude
 * @param y
 * @param x
 * @return {number}
 * @constructor
 */
const CHtoWGSlat = (y, x) => {
  // Converts military to civil and to unit = 1000km
  // Auxiliary values (% Bern)
  const yAux = (y - 600000) / 1000000
  const xAux = (x - 200000) / 1000000

  // Process lat
  const lat = 16.9023892 +
    (3.238272 * xAux) -
    (0.270978 * Math.pow(yAux, 2)) -
    (0.002528 * Math.pow(xAux, 2)) -
    (0.0447 * Math.pow(yAux, 2) * xAux) -
    (0.0140 * Math.pow(xAux, 3))

  // Unit 10000" to 1" and converts seconds to degrees (dec)
  return lat * 100 / 36
}

/**
 * Converts CH1903(+) to Longitude
 * @param y
 * @param x
 * @return {number}
 * @constructor
 */
const CHtoWGSlng = (y, x) => {
  // Auxiliary values (% Bern)
  const yAux = (y - 600000) / 1000000
  const xAux = (x - 200000) / 1000000

  // Process lng
  const lng = 2.6779094 +
    (4.728982 * yAux) +
    (0.791484 * yAux * xAux) +
    (0.1306 * yAux * Math.pow(xAux, 2)) -
    (0.0436 * Math.pow(yAux, 3))

  // Unit 10000" to 1 " and converts seconds to degrees (dec)
  return lng * 100 / 36
}

/**
 * Convert CH1903(+) to WGS84 (Latitude/Longitude)
 * @param y
 * @param x
 */
export default (y, x) => [
  CHtoWGSlat(y, x),
  CHtoWGSlng(y, x)
]

Je peux maintenant l’utiliser dans mon application principale:

import { open } from 'shapefile'
import ch1903ToWgs from './js/ch1903ToWgs'

open('assets/shapefiles/alpine-convention/shapefile.shp').then(source => {
  source.read().then(function apply (result) { // Read feature from the shapefile
    if (result.done) { // If there's no features left, abort
      return
    }

    // Convert CH1903 to WGS84
    const coords = result.value.geometry.coordinates[0].map(coordPair => {
      return ch1903ToWgs(coordPair[0], coordPair[1])
    })

    console.log(coords)

    return source.read().then(apply) // Iterate over result
  })
})

Ce qui donne les coordonnées ajustées suivantes:

C’est beaucoup mieux. Maintenant, je vais les superposer avec une carte pour obtenir quelque chose que je peux vraiment regarder.

Rendre les choses visibles

Pour cela, j’ajoute leaflet à l’application avec les tuiles Positron (lite) de CartoDB. Le thème Positron (lite) est gris clair/blanc, ce qui offre un bon contraste avec les caractéristiques que je vais afficher par-dessus. OpenStreetMap est génial mais trop coloré, ce qui fait que les polygones ne sont pas suffisamment visibles.

import leaflet from 'leaflet'
import { open } from 'shapefile'
import ch1903ToWgs from './js/ch1903ToWgs'

/**
 * Add the map
 */
const map = new leaflet.Map(document.querySelector('#map'), {
  zoomControl: false, // Disable default one to re-add custom one
}).setView([46.8182, 8.2275], 9) // Show Switzerland by default

// Move zoom control to bottom right corner
map.addControl(leaflet.control.zoom({position: 'bottomright'}))

// Add the tiles
const tileLayer = new leaflet.TileLayer('//cartodb-basemaps-{s}.global.ssl.fastly.net/light_all/{z}/{x}/{y}.png', {
  minZoom: 9,
  maxZoom: 20,
  attribution: '&copy; CartoDB basemaps'
}).addTo(map)

map.addLayer(tileLayer)

/**
 * Process the shapefile
 */
open('assets/shapefiles/alpine-convention/shapefile.shp').then(source => {
  source.read().then(function apply (result) { // Read feature from the shapefile
    if (result.done) { // If there's no features left, abort
      return
    }

    // Convert CH1903 to WGS84
    const coords = result.value.geometry.coordinates[0].map(coordPair => {
      return ch1903ToWgs(coordPair[0], coordPair[1])
    })

    console.log(coords)

    return source.read().then(apply) // Iterate over result
  })
})

J’obtiens déjà une belle carte avec laquelle je peux travailler:

L’étape suivante consiste à connecter la carte au chargement du shapefile. Pour cela, je crée des polygones leaflet à partir des coordonnées que j’ai précédemment transformées en WGS84:

// ...

// Convert CH1903 to WGS84
const coords = result.value.geometry.coordinates[0].map(coordPair => {
  return ch1903ToWgs(coordPair[0], coordPair[1])
})

const leafletPolygon = new leaflet.Polygon(coords, {
  weight: 0.5,
  color: '#757575',
  fillOpacity: 0.3,
  opcacity: 0.3,
})

leafletPolygon.addTo(map)

// ...

Et je regarde le résultat:

Sympa!

Le résultat final

Avec un peu plus d’interface utilisateur·trice et l’ajout d’un bouton à bascule pour les multiples shapefiles préalablement téléchargés, je peux créer une petite application qui me montre les zones à risque d’avalanche et les territoires où vivent les bouquetins en Suisse :

Nous pouvons voir les deux shapefiles en action : les polygones vert clair, vert, jaune et rouge sont des communes présentant un certain risque d’avalanche (vert clair = faible, vert = moyen-faible, jaune = moyen, rouge = élevé), les polygones plus foncés en haut représentent les colonies de bouquetins.

Cette carte offre déjà un bon aperçu du comportement des bouquetins : apparemment, ils vivent seulement dans des zones à risque d’avalanche faible à moyen-faible. Ils ont donc tendance à éviter les zones à risque moyen et élevé. Étant donné que les bouquetins préfèrent les terrains alpins et que ces terrains présentent généralement des risques d’avalanches, cette affirmation est plausible. Je dispose maintenant de données visibles qui la confirment.

L’application terminée est disponible ici, et son code source se trouve dans ce répertoire.

Les pièges

Bien que les shapefiles soient un moyen formidable de traiter les données SIG, il y a quelques pièges à éviter. L’un d’entre eux, que j’ai déjà mentionné, est un CRS inattendu. Dans la plupart des cas, il peut être identifié au point d’utilisation, mais quand on consulte un shapefile, il est recommandé de vérifier quel CRS il utilise. Le deuxième piège à éviter est la taille des shapefiles. Lors de l’utilisation de certains shapefiles avec JavaScript directement dans le navigateur, ce dernier peut « planter » car il tente de traiter la quantité astronomique de polygones. Il existe cependant des solutions : on peut soit simplifier le shapefile en supprimant les polygones inutiles, soit le prétraiter et stocker ses polygones dans une sorte de base de données et n’interroger que les polygones actuellement visibles.

Éléments à retenir

Personnellement, j’adore travailler avec les shapefiles. Les domaines d’utilisation sont nombreux et l’exploration des données qu’ils contiennent est vraiment passionnante. Même s’il peut y avoir un ou deux pièges, en général, c’est plaisant, car on peut obtenir quelque chose d’abouti sans trop d’efforts. La communauté OpenData utilise beaucoup les shapefiles, c’est donc une norme bien établie. Il existe aussi un grand nombre de bibliothèques et d’outils qui rendent le travail avec les shapefiles vraiment génial.

Tu as une question, un commentaire ou cet article t’as simplement inspiré? Mentionne-nous ou partage l’article sur Mastodon, Twitter ou LinkedIn.