Créer un fichier raster depuis cadastre.gouv.fr🔗

Posted by Médéric Ribreux 🗓 In blog/ OpenStreetMap/

#osm

Introduction

La Direction Générale des Finances Publique fait vraiment un travail formidable sur le domaine de l'Opendata et de l'information géographique. En effet, elle produit, dans un silence relatif si l'on compare à ce qui se passe sur les différentes plateformes Opendata des collectivités locales, un outil essentiel à la connaissance de notre territoire. Cet outil, c'est le Plan Cadastral. Avec la modernité, ce Plan Cadastral est maintenant numérique. Pour aller plus loin, il est même complètement disponible et gratuitement sur le site http://www.cadastre.gouv.fr que je vous recommande de consulter. Depuis quelques années maintenant, le projet OpenStreetMap a obtenu l'autorisation écrite d'utiliser les ressources de cadastre.gouv.fr pour générer du contenu dans OpenStreetMap.

Car il faut bien reconnaître que le cadastre est très intéressant pour cet usage. En effet, c'est un document qui est globalement assez complet, assez précis. Il permet notamment d'importer l'intégralité des bâtiments cadastrés, les emplacements des rues ainsi que les noms des lieux-dits et, pour finir, l'ensemble des surfaces en eaux. C'est également un document dont la couverture est quasi-totale et égale sur le territoire (à de rares exceptions près). Il permet donc de régler un des problèmes majeurs d'OpenStreetMap, à savoir, la faible couverture dans les communes rurales ou isolées qui ne bénéficient pas d'un grand nombre de "mappeurs" OpenStreetMap. Grâce au cadastre, nul besoin de se déplacer pour importer bâti et routes, voire adresses postales.

Si je devais vous montrer l'intérêt du cadastre dans OpenStreetMap, je prendrais un seul exemple sous forme d'images:

À gauche, le cœur d'une ville bien connue: San Francisco. À droite, un petit village perdu du Maine-et-Loire. La différence est saisissante: d'un côté on a bien les rues, de l'autre côté, il y a tout: les maisons, les rues, les lieux-dits, les cours d'eau, les étangs, le tout avec le nommage qui va bien. Il suffit d'un peu de motivation et surtout des bons outils pour réaliser ce travail d'import.

Sur ce point technique, certes il existe l'interrogation du site www.cadastre.gouv.fr qui présente un outil Web pour afficher le cadastre à un niveau de détail avancé. Cette application Web repose sur l'interrogation d'un serveur WMS accessible après connexion (cookie nécessaire). Néanmoins il existe également un plugin dédié dans JOSM, nomme cadastre-fr. Ce dernier agit en téléchargeant des dalles rasters depuis le serveur WMS du cadastre.

Néanmoins, la solution offerte par JOSM pour télécharger le cadastre comporte quelques inconvénients:

Toutes ces raisons m'ont amené à me pencher sur une solution permettant de récupérer une dalle raster du cadastre indépendamment de JOSM pour permettre une visualisation instantanée (ou rapide du moins) dans cet outil pour pouvoir facilement importer des éléments du cadastre. J'ai donc essayé de mettre au point un outil permettant de récupérer ces données directement depuis le site cadastre.gouv.fr.

Utilité et couverture du besoin

L'objectif principal est de créer une dalle raster depuis le serveur WMS de cadastre. Cela permet de disposer d'un document unique manipulable par n'importe quel logiciel de SIG, après récupération.

Voici un bref cahier des charges pour notre outil:

Méthode

Après avoir passé en revue nos besoins, voici une petite implémentation de cet outil. J'ai commencé à le réaliser en Python3 avant de m'apercevoir que je pouvais tout faire bien plus rapidement avec un simple script Shell (sous Bash). En effet, une requête WMS n'est rien d'autre qu'une requête HTTP GET bien construite qui retourne un fichier image. Des outils comme wget ou curl savent justement bien faire ce genre de travail.

La complexité vient du fait que le serveur WMS du cadastre requiert un cookie pour donner une réponse correcte. Il est donc indispensable de suivre un minimum la procédure de connexion au site cadastre.gouv.fr et son formulaire de recherche. Requêter ce dernier permet de récupérer le cookie qui nous intéresse ainsi que de vérifier que la commune demandée existe. Une fois munis du cookie et de l'existence de la commune, il faut au préalable récupérer la page HTML de navigation WMS, car cette dernière contient les coordonnées de l'étendue de la commune (la bounding-box). Une fois qu'on dispose de cette étendue, on peut commencer le téléchargement WMS.

Mais au préalable, il faut requêter intelligemment le serveur WMS. En effet, si on souhaite utiliser le cadastre, c'est pour avoir un niveau de détails assez fin. En conséquence, l'étendue de chaque requête WMS sera donc faible. Ainsi, il nous faudra découper notre étendue globale en autant de dalles WMS que nécessaires. Pour avoir un bon niveau de détails, j'ai utilisé le même pas que celui du plugin cadastre-fr de JOSM, à savoir un carré de 100m de côté. La tâche suivante consiste donc à faire une boucle qui permet de requêter toutes ces dalles WMS qui donnent pour chacune d'entre elles une image PNG. J'utilise le format PNG car l'algorithme de compression est non destructif ce qui permet d'avoir des images plus nettes, sans artéfacts de compression (comme le Jpeg). Pour une commune, le nombre de fichiers PNG peut atteindre les 8000 ! Cela donnera donc 8000 requêtes WMS ! Pour chaque fichier PNG, on va également générer un fichier pngw, un "worldfile" qui permet de préciser les coordonnées de l'emprise de chaque fichier PNG. C'est indispensable pour pouvoir effectuer un assemblage de ces nombreuses dalles rasters.

Une fois tous les fichiers récupérés, on peut utiliser les outils courants de GDAL pour assembler les fichiers en une seule dalle en deux étapes:

À la fin, on peut ajouter des pyramides (images pré-calculées avec une définition moindre, directement insérées dans le fichier de dalle) pour accélérer la visualisation avec les outils courants de SIG.

Une fois terminé, il faut bien entendu faire du ménage pour ne garder que l'unique dalle finale.

Pour fonctionner le script s'appuiera sur les outils suivants:

sous Debian, vous aurez donc besoin d'installer le paquet gdal-bin ainsi que curl.

Le Code

Après les concepts, voici le code. Je le place en lecture directe pour faciliter sa compréhension. Un simple copier coller et il est à vous !

#!/bin/bash

# Script de téléchargement de dalles raster de communes depuis cadastre.gouv.fr

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

# TODO:
# Gestion de la projection

# Vérification de la présence des outils indispensables: sed, gdalbuildvrt, gdal_translate, curl
CMDS="curl sed gdalbuildvrt gdal_translate gdaladdo"
for i in $CMDS
do
command -v $i >/dev/null && continue || { echo "Vous devez installer $i pour que ce script fonctionne…"; exit 1; }
done

if [ ! -n "$1" ]
then
  echo "Usage: $(basename $0) DEP COMMUNE"
  exit 1
else
  if [ "${#1}" -ne "3" ]
  then
echo "Le département doit être écrit sur 3 caractères."
echo "Ex: 049 pour le Maine-et-Loire."
exit 1
  fi
  DEP=$1
  COM=$(echo $2 | tr a-z A-Z)
fi

echo "Téléchargement de la commune $COM du département $departement…"
cook=$(mktemp cookie.XXXXXX)
curl_options="-s -c $cook"
rm $cook

# Gestion du proxy
if [ -n "$http_proxy" ]
then
  echo "Utilisation du proxy: $http_proxy"
  curl_options="$curl_options --proxy-anyauth -x $http_proxy"
fi

# Etape 1: Connexion et récupération du cookie de session
curl $curl_options "http://www.cadastre.gouv.fr/scpc/rechercherPlan.do" >> /dev/null
curl_options="-b $cook $curl_options"

# Etape 2: Vérifier que la commune existe dans le département
data="numeroVoie=&indiceRepetition=&nomVoie=&lieuDit=&ville=${COM}&codePostal=&codeDepartement=${DEP}&nbResultatParPage=10&x=45&y=11"
code_com=$(curl --data "$data" $curl_options http://www.cadastre.gouv.fr/scpc/rechercherPlan.do | grep -m 1 -o "afficherCarteCommune.do?c=[^']*" | head -n 1 | cut -d "=" -f 2)

# Récupération du code de la commune
echo "Code commune: $code_com"
#if not alf:
#    print("Desolé mais la commune {0} est introuvable dans le département {1}".format(args.COMMUNE, args.DEP))
#    exit(1)

#print("{0} a pour code commune: {1}".format(args.COMMUNE, code_commune))

# Etape 3: Connexion à la carte en ligne
# et récupération de la boundingbox
url="http://www.cadastre.gouv.fr/scpc/afficherCarteCommune.do?c=${code_com}&dontSaveLastForward&keepVolatileSession="
coords=$(curl $curl_options "$url" | grep -A 4 "new GeoBox(" | tail -n 4 | sed -e 's/,//g' -e 's/)//g' -e 's/\t//g' | tr '\n' ' ')
echo $coords | read -d " " -a bbox
OLDFS=$FS
FS=
read -a bbox <<EOF
  $coords
EOF
FS=$OLDFS

# Calcul des coordonnées de bounding box à prendre en compte
xmin=$(expr ${bbox[0]%.*} / 100)
xmin=$(expr $xmin \* 100);
xmax=$(expr ${bbox[2]%.*} / 100)
xmax=$(expr $xmax + 1)
xmax=$(expr $xmax \* 100)
ymin=$(expr ${bbox[1]%.*} / 100)
ymin=$(expr $ymin \* 100);
ymax=$(expr ${bbox[3]%.*} / 100)
ymax=$(expr $ymax + 1)
ymax=$(expr $ymax \* 100)

stat_x=$(expr $xmax - $xmin)
stat_x=$(expr $stat_x / 100)
stat_y=$(expr $ymax - $ymin)
stat_y=$(expr $stat_y / 100)
nb_dalles=$(expr $stat_x \* $stat_y)

echo "Téléchargement de $nb_dalles fichiers de dalles…"
# Etape 4: Boucle de récupération des fichiers au format PNG
x="$xmin"
i=0
while [ "$x" -lt "$xmax" ]
do
  y="$ymin"
  # Définition de la bbox à télécharger:
  a="${x}.0"
  c=$(expr $x + 100)
  c="$c.0"
  while [ "$y" -lt "$ymax" ]
  do
b="$y.0"
d=$(expr $y + 100)
d="$d.0"
this_bbox="${a},${b},${c},${d}"
# récupération du fichier PNG
url="http://www.cadastre.gouv.fr/scpc/wms?version=1.1&request=GetMap&layers=CDIF:LS3,CDIF:LS2,CDIF:LS1,CDIF:PARCELLE,CDIF:NUMERO,CDIF:PT3,CDIF:PT2,CDIF:PT1,CDIF:LIEUDIT,CDIF:SUBSECTION,CDIF:SECTION,CDIF:COMMUNE&format=image/png&bbox=${this_bbox}&width=1000&height=800&exception=application/vnd.ogc.se_inimage&styles=LS3_90,LS2_90,LS1_90,PARCELLE_90,NUMERO_90,PT3_90,PT2_90,PT1_90,LIEUDIT_90,SUBSECTION_90,SECTION_90,COMMUNE_90"
filename="$i.png"
echo "Ecriture de $filename …"
curl -o $filename $curl_options "$url"

# Gestion du world file:
pngw="${filename}w"
echo -e "0.10\n" > $pngw
echo -e "0\n" >> $pngw 
echo -e "0\n" >> $pngw
echo -e "-0.125\n" >> $pngw
echo -e "$a\n" >> $pngw
echo -e "$d\n" >> $pngw

y=$(expr $y + 100)
i=$(expr $i + 1)
  done
  x=$(expr $x + 100)
done

# Fabrication du fichier VRT
gdalbuildvrt ${COM}.vrt *.png

# Génération au format GeoTiff
gdal_translate -of GTiff -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=IF_NEEDED -co ZLEVEL=9 ${COM}.vrt ${COM}.tiff

# Ajout des pyramides au GéoTiff
gdaladdo -r average ${COM}.tiff 2 4 8 16

# Nettoyage:
rm $cook *.png *.pngw

exit 0

Utilisation et limites

L'utilisation de l'outil est très simple, la commande suivante:

$ ./wms_cadastre.sh 049 ANGRIE

permet de télécharger un raster cadastre dans le répertoire courant de la commune ANGRIE située dans le département 049.

Le téléchargement va prendre pas mal de temps puisqu'il faudra réaliser de nombreuses requêtes WMS. L'assemblage est, en comparaison, assez rapide avec les options de création. Il faudra toutefois un peu plus de temps pour générer les pyramides.

Cet outil est loin d'être parfait (mais vous pouvez l'améliorer). Déjà, il est intolérant aux pannes serveurs ou de connexion réseau: si votre requête plante au milieu d'un téléchargement, il vous manquera vraisemblablement un fichier ce qui occasionnera un trou dans la dalle. Mais dans la pratique, il est tout à fait possible de télécharger 8000 dalles sans incidents.

Un deuxième inconvénient majeur, se situe au niveau de la transparence. Lorsqu'on télécharge une commune donnée, on télécharge tout ce qui se trouve dans une emprise qui est rectangulaire. Les espaces situés en dehors des limites communales retournent du blanc. Ainsi, on télécharge une grosse dalle par commune et il est impossible de l'assembler avec celle d'une autre commune qui serait limitrophe. Il est donc difficile (mais pas impossible si l'on joue avec les options de transparence de GDAL) d'obtenir une dalle raster d'un groupe de communes ou d'un département.

Enfin, le script ne gère pas la projection. Pour ma part, comme je compte travailler sur un secteur limité, la projection sera toujours identique (RGF CC47 soit le code EPSG:3947). Vous pouvez toujours modifier le code EPSG dans la commande gdal_translate si vous souhaitez travailler sur une autre zone.

Intégration dans JOSM

Nous disposons maintenant d'une dalle raster du cadastre d'une commune. Il faut pouvoir la visualiser dans JOSM. J'ai testé plusieurs méthodes pour y parvenir. La première consistait à utiliser un plugin nommé ImportImagePlugin. Son but est simplement de charger une dalle raster géo-référencée dans JOSM. Mais ce plugin est complètement bugué et je n'ai jamais réussi à le faire fonctionner.

Une autre technique consiste à fabriquer un ensemble de tuiles incorporables dans JOSM grâce à sa gestion des couches TMS. Gdal fournit un utilitaire nommé gdal2tiles.py pour faire ce travail de découpage. Mais encore une fois, je n'ai pas réussi à faire correspondre les tuiles à celles requêtées par JOSM. Il y avait toujours un décalage sur les coordonnées y. Impossible de corriger ça, même à la main. De plus, gdal2tiles n'est pas très performant il lui faut quelques heures pour générer un tas de tuiles proposant un bon niveau de détails.

J'ai finalement mis la main sur la bonne méthode qui consiste à utiliser le protocole WMS. Le principe est assez simple: on installe un serveur WMS sur sa machine (en localhost) et on lui demande de servir la couche raster du cadastre téléchargée. Bien entendu, il faut un serveur WMS. Le plus simple que j'ai trouvé est celui qui est installé avec Qgis. Dans Debian, il suffit d'installer le paquet qgis-mapserver. Ce dernier installe un programme accessible en CGI (vous aurez besoin d'un serveur HTTP. Apache2 sera votre ami). Si vous avez une configuration par défaut d'Apache, le serveur WMS est simplement accessible via l'url http://localhost/cgi-bin/qgis_mapserv.fcgi.

Le serveur WMS de Qgis est finalement très simple à configurer. Vous aurez besoin de créer un fichier de projet (.qgs) contenant votre couche raster, de renommer cette couche en CADASTRE dans QGis et d'enregistrer le fichier de projet en ayant soin d'indiquer la bonne projection qui sera celle que vous utiliserez dans JOSM. Je vous invite à utiliser le RGF CC de la commune que vous comptez importer. Bien entendu, le fichier de projet doit être situé dans un répertoire accessible en lecture à l'utilisateur Apache (www-data sous Debian).

Pour configurer la couche WMS dans JOSM, vous aurez besoin de faire un tour dans les Préférences (Menu Modifier → Préférences) et d'ouvrir l'onglet "WMS/TMS". Ajoutez l'URL vers la bonne couche WMS. En voici un exemple:

wms:http://localhost/cgi-bin/qgis_mapserv.fcgi?MAP=/opt/OSM/WMS.qgs&FORMAT=PNG&VERSION=1.1.1&SERVICE=WMS&REQUEST=GetMap&LAYERS=CADASTRE&STYLES=&SRS={proj}&WIDTH={width}&HEIGHT={height}&BBOX={bbox}

Décomposons le contenu de cette URL pour mieux la comprendre:

Une fois la couche WMS configurée, vous pouvez l'activer dans JOSM et commencer à importer votre cadastre.

En termes de format, j'ai remarqué que les données les plus rapides sont celles qui sont issues d'une source en GéoTiff avec pyramides internes. A l'utilisation, sur une machine moyenne, l'affichage de la couche WMS est assez rapide et dans tous les cas, bien plus qu'en passant par le plugin cadastre-fr.

Il vous faudra modifier votre fichier de projet Qgis (.qgs) lorsque vous souhaiterez travailler sur une autre commune.

Conclusion

Avec le script exposé dans cet article, il est tout à fait possible de récupérer un fichier raster de cadastre et de l'afficher dans JOSM. Le temps de capture de la dalle raster peut sembler un peu long, mais il faut le relativiser par rapport au travail d'import manuel du contenu du cadastre dans OpenStreetMap. En effet, il faut compter à peu près 1 heure de travail pour intégrer les routes, les rues et les lieux-dits d'une petite commune. Pendant ce temps, rien n'empêche de télécharger une autre commune en tâche de fond !