S'abonner à un flux RSS
Classement en iso-valeurs des données LIDAR
.jpg)
Sommaire |
Contexte
Cette page est un sous-thème de la page sur l'Utilisation des données LIDAR pour la directive inondation. Le lecteur est invité à s'y reporter pour comprendre le contexte général sur l'exploitation des données LIDAR.
Un préalable à l'utilisation des fonctions définies est la réalisation des Préalables pour l'utilisation de Qgis-GRASS sur le LIDAR.
Cette page a pour objet de fournir rapidement une méthode pour réaliser des polygones d'iso-valeurs à partir de classement de raster (et non d'interpolations lourdes en temps de calcul) comme:
- la création de courbes en NGF à équidistance
- les zones basses littorales définies sur des niveaux marins extrêmes (méthode dans l'étude « Vulnérabilité du territoire national aux risques littoraux » réalisée par le CETMEF, le CETE Méditerranée et le CETE Ouest
- la création de classes de hauteurs d'eau.
Le résultat de ce type de traitement est proposé sur les cas suivants (Images en cours de réalisation):
- Classement tous les 10m sur le MNT IGN 1m du département du Var levé après Juin 2010
- Classement tous les 10m sur le MNT IGN 1m de Mayotte
- Classement des zones littorales du Golfe du Lion
- Classement des zones levées suite à la tempête Xynthia
- Classement sur les surfaces drainées de l'opération d'Exzeco sur un MNT créé au pas de 5m sur le département du Gard
Principes
La méthode retenue est fournie en annexe 5 du rapport cité ci-dessus pour les DOM-COM.UNIQ3832793f3a0ac951-anyweb-00000000-QINU
Au lieu d'interpoler des courbes d'iso-valeurs, un classement est réalisé. Le classement est plus pertinent car moins consommateur de ressource informatique et l'écart géométrique entre les deux résultats n'excède pas la largeur de la maille de la grille (plus exactement la largeur de maille multipliée par racine carrée de 2).
Les algorithmes classiques sont les suivants:
- Classement du fichier raster (MNT par exemple) en iso-classes (de 0 à 1m NGF dans la classe 1, de 1 à 2 dans la classe 2, de 2 à 5 dans la classe 3)
- Conversion de ce raster classé en polygone sans lissage
- Ajout des métadonnées (producteur, maitre d’ouvrage, date, minimum et maximum de chaque classe), et les informations de couleurs
Ce travail fait sur une dalle peut être automatisé sur toutes les dalles, la jonction géométrique est parfaite.
Le résultat final se compose de:
- Des tables sous forme vecteur sur la même emprise que chaque dalle qui peuvent être ouvertes à partir d'une table d'assemblage, en utilisant des fonctions de type hotlink d'outil SIG, ces données permettent une grande facilité de lecture et une gestion simple et légère
- une table fusionnée de tous ces vecteurs qui peut être lourde en fonction de la taille de la zone traitée (elle le sera beaucoup moins que les dalles LIDAR initiales)
En fonction de la qualité de rendu que l'on souhaite, il peut être utile de simplifier le résultat final. Quelques points sont fournis sur cette probélamtique en fin de page.
Méthode directe dans Qgis-GRASS
Classement du raster GRASS
Pour cela, on peut utiliser l'outil r.mapcalc. Une fenêtre s'ouvre ou la ligne de commande est celle-ci: r.mapcalc "IGN_MNT_1m_0981_6271GRASS1000 = ((IGN_MNT_1m_0981_6271GRASS@Class*1000))"
L'outil r.reclass demande ensuite un fichier de relation, à créer au préalable, les classes voulues étant tous les 10 mètres, donc pour notre raster multiplié par 1000 cela donne:
0 thru 10000 = 1
10000 thru 20000 = 2
20000 thru 30000 = 3
30000 thru 40000 = 4
40000 thru 50000 = 5
...
Le raster en sortie aura le nombre de classes choisies en sortie, on peut faire des intervalles variables...Pour l'automatisation, la jointure (...), il semble plus intéressant de ne pas utiliser la fonction NULL afin d'avoir le vecteur sur l'ensemble de l'emprise du Raster, certains "trous" dans la géométrie étant parfois mal accepté dans certains outils SIG.
L'onglet manuel permet de voir les différents formats ou règles à constituer.
ATTENTION LE CLASSEMENT SE FAIT PAR VALEURS ENTIERES, IL FAUDRA TOUT D'ABORD MULTIPLIER PAR UN GRAND NOMBRE.
Conversion du raster classé en vecteur
Le raster classé peut ensuite être converti en vecteur avec la routine r.to.vect ATTENTION, IL NE FAUT PAS GARDER LE RESULTAT DE L'INTERFACE GRASS CAR L'OPTION PAR DEFAUT LISSE LES BORDS CE QUI AUGMENTE LE POIDS DU FICHIER FINAL ET N'EST PAS EN ACCORD AVEC LE RASTER INTIAL Pour faire ce travail, il est conseillé de passer dans l'interface GRASS r.to.vect.area Quand le travail est terminé, dans l'onglet « Rendu », copier la ligne de commande sans le « flag » -s à la fin (voir sa signification dans l'onglet manuel) et changer une lettre du nom (il ne veut pas écraser). Ouvrir l'interface SHELL Console GRASS, coller (clic droit de la souris et lancer Un autre moyen de taper r.to.vect dans la console , bienvenue dans le monde GRASS...
Jointure avec une table de méta-données
UNIQ3832793f3a0ac951-anyweb-00000001-QINU Dans Qgis, il suffit d'ouvrir par exemple un fichier csv dans lequel une colonne servira de jointure, d'aller dans les propritété du vecteur classé , l'onglet jointure. Choisir ajouter afin de faire la liaison.
Pour faciliter la jointure et le formatage des colonnes, il est utile d'y associer un fichier csvt qui permet la déclaration des types de champs. Un exemple de ces fichiers et fourni dans la partie automatisation ci-dessous.
Dans GRASS, la procédure est plus complexe, elle est fournie dans le fichier pdf joint. Cette procédure est utilisée dans la méthode d'automatisation.
Retour d'expérience
Il ne reste plus qu'à exporter au format vecteur qui vous va bien dans l'arborescence des outils GRASS et dans le cas de submersions marines croiser avec des zones de niveaux marins extrêmes dans votre outil SIG, si cette méthode convient.
Un test a été fait avec 32 dalles raster de 7,5Mo chacune, le résultat est un seul fichier Mapinfo de moins de 1Mo.
Le temps de travail et calcul était de 1/2 heure.
Automatisation dans Qgis-GRASS
Avant l'automatisation, il faut faire un nouveau jeu de donnée GRASS (le faire sur la première dalle qui vous sert pour copier-coller vos lignes de commande), ne pas oublier de mettre la taille de votre maille dans la définition de la région.
Pour automatiser cette procédure à de nombreuses dalles, le principe est de prendre les fonctions utilisées précédement ou faite sur votre première dalle avec des copiers-collers:
- Importer une dalle
r.in.gdal -o input=D:/Topo/Dpt83-Juin2010/MNT/IGN_MNT_1m_0961_6270.asc output=IGN_MNT_1m_0961_6270GRASS
- Réduire la région sur votre dalle, cette étape permet de définir l'emprise de GRASS et améliore de manière spectaculaire les perfrormances en temps de calcul
g.region rast=IGN_MNT_1m_0961_6270GRASS@Class
- Multiplier votre raster par un grand nombre comme 1000 car le classement ne se fait que sur des entiers
r.mapcalc "IGN_MNT_1m_0961_6270GRASS1000 = ((IGN_MNT_1m_0961_6270GRASS@Class*1000))"
- Classer ce raster (fichier de reclassement de 0 à 1000 par par de 10 mètres qui est fonction de votre multiplication précédente)
r.reclass input=IGN_MNT_1m_0961_6270GRASS1000@Class rules=D:/Topo/Dpt83-Juin2010/MNT/fichierclassement.txt output=IGN_MNT_1m_0961_6270GRASSclasse
- Le convertir en polygone
r.to.vect input=IGN_MNT_1m_0961_6270GRASSclasse@Class output=IGN_MNT_1m_0961_6270GRASSclassevecteur feature=area
- L'exporter dans un format vecteur comme le format shape
v.out.ogr input=IGN_MNT_1m_0961_6270GRASSclassevecteur@Class type=area layer=1 format=ESRI_Shapefile dsn=D:/Topo/Dpt83-Juin2010/MNT/IGN_MNT_1m_0961_6270vect.shp olayer=default -e -c
- Supprimer tous les fichiers créés dans GRASS
g.remove vect=IGN_MNT_1m_0961_6270GRASSclassevecteur@Class
g.remove rast=IGN_MNT_1m_0961_6270GRASSclasse@Class
g.remove rast=IGN_MNT_1m_0961_6270GRASS1000@Class
g.remove rast=IGN_MNT_1m_0961_6270GRASS@Class
Ensuite, si vous avez 10 lignes, vous copiez 10 fois votre liste de dalles à la suite dans un tableur, vous triez par le nom et vous modifiez dans les lignes de commandes uniquement les noms des dalles.
Vous avez une liste de la forme suivante:
r.in.gdal -o input=D:/Topo/Dpt83-Juin2010/MNT/IGN_MNT_1m_0961_6270.asc output=IGN_MNT_1m_0961_6270GRASS
g.region rast=IGN_MNT_1m_0961_6270GRASS@Class
r.mapcalc "IGN_MNT_1m_0961_6270GRASS1000 = ((IGN_MNT_1m_0961_6270GRASS@Class*1000))"
r.reclass input=IGN_MNT_1m_0961_6270GRASS1000@Class rules=D:/Topo/Dpt83-Juin2010/MNT/fichierclassement.txt output=IGN_MNT_1m_0961_6270GRASSclasse
r.to.vect input=IGN_MNT_1m_0961_6270GRASSclasse@Class output=IGN_MNT_1m_0961_6270GRASSclassevecteur feature=area
v.out.ogr input=IGN_MNT_1m_0961_6270GRASSclassevecteur@Class type=area layer=1 format=ESRI_Shapefile dsn=D:/Topo/Dpt83-Juin2010/MNT/IGN_MNT_1m_0961_6270vect.shp olayer=default -e -c
g.remove vect=IGN_MNT_1m_0961_6270GRASSclassevecteur@Class
g.remove rast=IGN_MNT_1m_0961_6270GRASSclasse@Class
g.remove rast=IGN_MNT_1m_0961_6270GRASS1000@Class
g.remove rast=IGN_MNT_1m_0961_6270GRASS@Class
r.in.gdal -o input=D:/Topo/Dpt83-Juin2010/MNT/IGN_MNT_1m_0962_6269.asc output=IGN_MNT_1m_0962_6269GRASS
g.region rast=IGN_MNT_1m_0962_6269GRASS@Class
r.mapcalc "IGN_MNT_1m_0962_6269GRASS1000 = ((IGN_MNT_1m_0962_6269GRASS@Class*1000))"
r.reclass input=IGN_MNT_1m_0962_6269GRASS1000@Class rules=D:/Topo/Dpt83-Juin2010/MNT/fichierclassement.txt output=IGN_MNT_1m_0962_6269GRASSclasse
r.to.vect input=IGN_MNT_1m_0962_6269GRASSclasse@Class output=IGN_MNT_1m_0962_6269GRASSclassevecteur feature=area
v.out.ogr input=IGN_MNT_1m_0962_6269GRASSclassevecteur@Class type=area layer=1 format=ESRI_Shapefile dsn=D:/Topo/Dpt83-Juin2010/MNT/IGN_MNT_1m_0962_6269vect.shp olayer=default -e -c
g.remove vect=IGN_MNT_1m_0962_6269GRASSclassevecteur@Class
g.remove rast=IGN_MNT_1m_0962_6269GRASSclasse@Class
g.remove rast=IGN_MNT_1m_0962_6269GRASS1000@Class
g.remove rast=IGN_MNT_1m_0962_6269GRASS@Class
etc...
L'opération ou plutôt les 10 lignes de commandes par dalles prennent 6 à 7 secondes.
La jointure de méta-données dans chaque fichier vecteur résultat de chaque dalle est en cours de recherche.
A ce stade (sans la jointure), le travail sur 701 dalles prend 1h23.
Vous pouvez ensuite assembler tous vos fichiers vecteurs dans Qgis (Menu Vecteur/Outil de gestions de données/Fusionner les shapefile en un seul).
Cette phase prend moins d'une minute, le fichier résultat fait 330 Mo comparé aux 4.71 Go initiaux.
Vous pouvez réaliser une jointure avec des méta-données à partir de ce fichier fusionné.
Automatisation avec jointure dalle à dalle
L'automatisation avec une jointure dalle à dalle est possible mais plus complexe, il est utile de tester sur quelques dalles la procédure précédente. Le fichier pdf visible dans wikhydro montre la procédure, le lien permet d'accéder au tableau. Si ce lien ne fonctionne plus, merci de contacter un des auteurs de cette page.
Finalisation des résultats
Lissage
Les quelques points suivants donnent notre avis sur le problème de "pixellisation" possible en fin de traitement et des options pour le gérer:
- De manière générale, il nous semble inutile de dégrader l'information initiale par exemple en ré-échantillonant le raster au pas de 2m alors qu'il était initialement au pas de 1m pour "lisser" les résultats sauf si le travail nécessite du temps réel.
- Utiliser des algorithmes pour lisser les bords de polygones ne nous semblent pas utile
- Des filtres sur le raster avec la fonction r.neighbors sont possibles mais sont délicats à utiliser directement pour un travail dalle à dalle car les limites seront mal traitées. De plus, cette technique sur la donnée brute peut effacer les éléments parfois structurants mais fins comme des digues, remblais, présence de chenaux d'écoulements... Ils sont intégrés dans l'outil DICARTO après l'option de classement avec l'option mode qui permet de légèrement simplifier les géométries.
- L'option la plus intéressante est l'utilisation de la fonction v.clean avec la sous-option de "Supprimer de petites zones, la limite la plus longue avec la zone adjacente est supprimée. Elle est par contre couteuse en temps.
En fonction des besoins de qualité de restitution liés à l'échelle de rendu cartographique par exemple, différents seuils de l'option précédente pourront être utilisés. L'utilisation de plusieurs seuils à la suite en allant du plus petit au plus grand peut être réalisé.
L'exemple fourni montre sur un croisement topographie - niveaux marins sur une grille homogène de 3m le fonctionnement du processus. Le nettoyage progressif du phénomènes de pixelisation se fait en supprimant les petits objets qui ont la plus longue zone adjacente. Une amélioration sécuritaire à trouver pour la cartographie de hateur d'eau serait la plus forte valeur jointive.
Le codage ci-dessous est simple:UNIQ3832793f3a0ac951-anyweb-00000002-QINU
v.clean input="PSL9m@PSL" output="PSL18m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=18
v.clean input="PSL18m@PSL" output="PSL27m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=27
v.clean input="PSL27m@PSL" output="PSL36m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=36
v.clean input="PSL36m@PSL" output="PSL45m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=45
v.clean input="PSL45m@PSL" output="PSL54m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=54
v.clean input="PSL54m@PSL" output="PSL63m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=63
v.clean input="PSL63m@PSL" output="PSL72m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=72
v.clean input="PSL72m@PSL" output="PSL81m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=81
v.clean input="PSL81m@PSL" output="PSL90m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=90
v.clean input="PSL90m@PSL" output="PSL99m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=99
v.clean input="PSL99m@PSL" output="PSL108m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=108
v.clean input="PSL108m@PSL" output="PSL117m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=117
v.clean input="PSL117m@PSL" output="PSL126m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=126
v.clean input="PSL126m@PSL" output="PSL135m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=135
v.clean input="PSL135m@PSL" output="PSL144m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=144
v.clean input="PSL144m@PSL" output="PSL153m" type="point,line,boundary,centroid,area" tool="rmarea" thresh=153
Le temps est assez long, surtout pour la première étape (32000 objets) et devient très rapide ensuite pour obtenir 2000 objets pour une surface de suppression finale de 153m². La disparition de sorte de canaux pourrait être minimisée ou retardée avec l'utilisation d'un LIDAR avec un pas plus précis (1m au lieu de 3m).
Fourniture finalisée
Le travail a été réalisé en fichier shape. Pour le diffuser en interne au ministère, il est encore utile de le convertir au format Mapinfo dans la bonne projection. Un utilitaire mbx permet ensuite de "colorier" les différents objets avec les valeurs RVB fournies dans la table. Ce fichier est disponible sur demande.
La "coloriage" automatique est intégré avec les fichier qml dans les outils DICARTO.
Liens externes
- http://geoinfo.metier.i2/qgis-r625.html, site intranet MEDDE sur l'outil Qgis avec version ministère, Qgis - Documentation
Qgis - Contributions des services
Qgis - Groupe de travail national
Qgis - Plugin
Qgis - Fiche pratique
Qgis - Faq
Qgis - Le coin des développeurs
Packages d’installation
GRASS
- https://github.com/qgis, les sources du logiciel
- http://forum.qgis.org/index.php, le forum Qgis
- http://grass.osgeo.org/community/support.php, le forum GRASS
- http://grass.osgeo.org/wiki/GRASS-QGIS_relevant_module_list qui permet de connaitre les modules grass présent dans les outils GRASS à travers Qgis et ceux non disponibles mais parfois utiles
- http://www.linux-nantes.org/wws/info/grass-fr, liste de discussion GRASS française, que nous remercions pour sa collaboration
- http://fr.wikipedia.org/wiki/Wikip%C3%A9dia:Atelier_graphique/Didacticiels_cartographiques/Didacticiel_de_cr%C3%A9ation_de_carte_avec_QGIS, aide à la création de carte sous Qgis
Les intervenants
Maîtrise d'Ouvrage: MEDDM/DGPR
- Responsable de l'étude: Sabine Baillarguet
Réalisation: Centre d'Études Techniques de l'Équipement Méditerranée
- Responsable de l'étude: Frédéric Pons
- Intervenants: Mathieu Alquier, Isabelle Roux, Christine Alpin, Yves Bonin, Bertrand Bouteille, Christophe Esposito, Céline Trmal , Christophe Laroche
- Stagiaire: Arnaud Alvarez (TFE ESGT), Adeline Agostini (Stage Master II GERINAT)
Partenariats:
- CPII DO Ouest: René Jonchere
- CETMEF: Bruno Bader
- CETE Sud-Ouest: Yves Nedelec
- DGPR/SRNH/SCHAPI: Christophe Astier, Jean-Luc Souldadié
- IGN: Alain Duperet, Sébastien Saur
Article rédigé par Appui CETE/CETMEF pour la Directive Inondation piloté par le PCI Inondations et Aléas Côtiers
Note : d'autres personnes peuvent avoir contribué au contenu de cet article, [Consultez l'historique]. |
