Shapefiles : avalanches et bouquetins

  • Pascal Thormeier

Les shapefiles, ou fichiers de forme, sont un excellent moyen de stocker des donnĂ©es spatiales et gĂ©ographiques, comme des biens immobiliers, des limites communales, des riviĂšres, des lacs, des rues, etc. ainsi que des mĂ©tadonnĂ©es comme des noms, des identifiants, l’annĂ©e d’établissement, etc.

En savoir plus sur le service Open data de notre agence digitale.

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: '© 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.


Qu’en penses-tu?