Traitement de données LIDAR acquises en sous bois, extraction par méthodes tomographiques de données individuelles d'arbre, localisation, diamètre, volumétrie.

Le jeu de données initial,  les étapes de développement (Mise à jour 23/11/2019)

Les travaux relatifs à ce projet sont présentés dans l'article à lire ICI.

Un jeu d'écrans (une centaine), collecte de manière peu structurée plus d’informations sur la procédure de travail ICI

Jouons avec les mots, dans les travaux décrits ici, la topographie a son importance car on s'intéresse à des reliefs et des localisations. Mais la tomographie est tout aussi importante, car les traitements de nuage de points LIDAR, relevant de la 3D par nature, se font ici totalement en 2D sous forme d'images avec les traitements classique au processus d'imagerie, selon des coupes verticales, horizontales, ou en respect de tout autre plan défini par deux vecteurs.

Ce qui suit est pour mémoire et donne les sources et étapes de développement pour aboutir à une extraction automatique des attributs d’arbre en forêts. Hauteur, diamètre, géolocalisation auxquels peuvent être ajoutés les attributs les plus divers qui soient, sachant qu’avec cette approche tout arbre peut être identifié individuellement et suivi sur tout son cycle de vie. Le tout est évidemment gérable en SIG de préférence, et en base de données géométrique (PostGres/postgis, SpatiaLite).

L’essentiel des travaux de R&D et développement porte sur des relevés en sous bois effectués par la société EXAMETRICS, en utilsant un LIDAR terrestre GEOSLAM Z-REVO. Porté par un marcheur, un relevé de 10 ha se fait 5 jours, la durée est liée tant à la topographie du terrain qu'à son niveau d'embroussaillement. Ce premier relevé de caractère expérimental au sein de la RNN la Massane a fait l’objet de l’étude qui suit. Ceci a permis de mieux cerner des protocoles de travail à mettre en oeuvre pour l’avenir, afin de s’assurer de la meilleure qualité de données qui soit pour mettre en oeuvre les outils et algorithmes issus de ces travaux et  développements.

Bien qu’à vérifier plus profondément, nous pensons que nous proposons une approche qui devrait permettre une exhaustivité dendrométrique et cela à des conditions économiques relativement attractives avec une procédure de travail qui, quoique rigoureuse, est assez simple à mettre en oeuvre. Le lidar GEOSLAM de performances modestes somme toute mais suffisantes, fournit, en resserrant davantage le protocole d’acquisition,  un flot de données permettant de maîtriser les modèles numériques de terrain (MNT) et la géométrie des arbres juqu’aux premiers branchages car nous avons limité nos algorithme à ces premiers objectifs. Maintenant la voie est ouverte pour aller plus loin.

Nos algorithmes, relativement simples,  certainement à enrichir, s’appuient sur les fondamentaux de la géométrie 3D, vecteurs, plans,  n’allant pas au-delà des matrices homogènes de rotation voire des quaternions, Les méthodes statistiques se limitent à moyenne et écart-type et nous utilisons bien sûr quelques filtres de lissage ou de voisinage de type KNN.

Nous nous sommes également appuyés également sur des données (nuages LIDAR et Orthophotos) fournies à la RNN La Massane par la société SINTEGRA suite à relevés aériens, à titre de contrôle. Au-delà nous n’en avons pas fait grand usage, la résolution étant insuffisante.

Comme il s'agit de recaler des données issues de relevés sur le terrain, nous en donnons ci-après les diverses sources et approches de traitement. Nos travaux s’appuient exclusivement sur des logiciels libres, et nous remercions toutes les personnes ayant oeuvré à leur développement, citons CloudCompare, Point Cloud Library (PCL), Point Data Abstract Library (PDAL), LilBlas pour les principales. Les travaux présentés ci-dessous sont codés selon les besoins en C++, C_CUDA, Python.

  1. Un ensemble de données qui est une extraction de MapInfo fournit la cartographie issu d'identification manuelle effectuée sur le terrain, =>    /media/mbariou/LIBPROG/PROJ_RNN//MapInfo_Data, reçu en décembre 2018,
  2. des rapports d'études, 55, 62, 72, 80 et rapport lidar Sintegra, =>  /media/mbariou/LIBPROG/PROJ_RNN/LaMassane, reçu en décembre 2018,
  3. données récentes, extraction faites par EXAMETRICS + arbres sur pied LaMassane => /media/mbariou/LIBPROG/PROJ_RNN/DonneesRecentes/ reçu en mars 2018,
  4. Documentation projet => /media/mbariou/LIBPROG/PROJ_RNN/DocConstructionProj/,
  5.  Arbres sur pied, extraction MapInfo  => /media/mbariou/LIBPROG/PROJ_RNN/Arbres_Sur_Pied
  6. A la racine divers documents à reclasser => /media/mbariou/LIBPROG/PROJ_RNN/
  7. Nuages de points :
    1. nuage EXAMETRICS de la cloture et du refuge => ~/PointCloud/ClotureRefuge
    2. ensemble des relevés LIDAR Sintegra ~/PointCloud/laMassaneCanopée
    3. ensemble des relevés EXAMETRICS en sous-bois => ~/PointCloud/LaMassaneSol
    4. fichiers test et tuto SpatiaLITE => ~/PointCloud/LaMassaneSol
    5. OrthoPhoto Sintegra => /home/mbariou/WKS_QGIS/QGIS_ORG/Projets_en_Cours/06-Orthophoto_10cm
    6. Notes d'extraction de données  =>  /home/mbariou/WKS_QGIS/QGIS_ORG/Projets_en_Cours/Lamassane
  8. Informations WEB => brasnah.fr/articles/Extraction de données la Massane
  9. Transformarion 3D/Helmert => Transf3DLaMassane.cu .../WKSCUDA/C_CUDA/CodeSamples/chapter02/ ( nvcc -O2 -arch=sm_35 -std=c++11 Transf3DLaMassane.cu -o Transf3DLaMassane)
  10. Exraction de donnée par Zone -> Extract_Tile_from_las dans /PointCloud/Cloture_Refuge
  11. Création de MNT cellulaire => CreateStructCellMnt dans /PointCloud/Cloture_Refuge => Aec table Hash, essai abandonné
  12. Moyennage de MNT => CreateMnt dans /PointCloud/Cloture_Refuge (Moyennage par cellule à taille paramétrable et évaluation de l'écart-type) organisation en conteneur Map et mise en fichier
  13. Test de conteneur MAP =>mapCellPrj  .../WKSCPLUSPLUS/ProgCode/Project/mapCellPrj Créationde Map et sérialisation pour mise en fichier, texte ou binaire
  14. Test de précision de calcul de cacul pour les recalages de nuage Amer/Amer homologue, Triangle sur triangle :
    1. /home/mbariou/WKS_PYTHON/SimuQuaternion/transformations.py (addition quaternion) + couture nuage ,
    2. /home/mbariou/WKS_PYTHON/SimuQuaternion/TestQuaternionTriangle.py Triangles homologue, deux nuages dans un même référentiel
    3. /home/mbariou/WKS_PYTHON/SimuQuaternion/RebuiltTriangleForMatTrsf.py bibliothèque de travail, construction des triangles plan en nuage de points (méthode barycentrique), évaluation de la matrice de transformation, contrôle de la précison, métrique entre deux nuages
    4. /home/mbariou/WKS_PYTHON/SimuQuaternion/RebuiltTriangleForMatTrsf.py mise en oeuvre des transformations et contrôle, extraction de la matrice de transformation pour reversement dans CUDA
  15. CUDA ~/WKSCUDA/C_CUDA/CodeSamples/Chapter02/LidarTransfTxT.cu Code de transformation de nuage de points, précision insuffisante mais rapide:
    1. nvcc -O2 -arch=sm_35 -std=c++11 ~/WKSCUDA/C_CUDA/CodeSamples/chapter02/LidarTransfTxT.cu -o ~/WKSCUDA/C_CUDA/CodeSamples/chapter02/LidarTransfTxT
    2. ~/WKSCUDA/C_CUDA/CodeSamples/chapter02/LidarTransfTxT
    3. ~/WKSCUDA/C_CUDA/CodeSamples/chapter02/LidarTransfTxTDouble.cu prise en compte des erreurs de calcul, passage en double (envisager la translation en calcul entier grand nombre)
  16. Voir les fichiers JSON dans ~/PointCloud/LaMassaneSol notamment la création des cellues numériques de 4 m2 BigFile & cellsOfBig
  17. Finir le reccensement des fichiers *.py, *.cpp, *.cu et rappeler leur fonctionnalité et partir sur PyQgis pour l'intégration
  18. Carroyer à l'identique le nuage de point en sous bois et le modèle numérique de terrain en envisageant sur ce dernier des cellules de 2 m2 (2 x 2) et 25 m2 (5 x 5)
  19. Creation de MinMaxTiling carroyage en place, en une seule passe  du nuage de point sur la base de définition de nombre de cellule (M x N) orthogonales au référentiel X,Y, accès aléatoire
  20. Dérivé de MinMaXTiling => MinMaxVoxel découpage en mode non récursif, mais en une seule passe  et en place du nuage de point sur la base d'une définition d'une matrice 3D (M x N xP), accès aléatoire 
  21. extract-Tile_From_Las assure extraction (parallélipipède),  d’un nuage de points,  d’une section horizontale ou verticale, ce qui permet de coduire des traitements de type tomographique sur des tranches de nuages de points totalement paramétrables.
  22. createMNT_Optimal_CM, partant d’un fichier *.las de mesure, permet d’extraire le MNT correspondant, selon les besoins la résolution permet de traiter des cellules de 1 cm2 (1 x1). deux filtres un premier de lissage pondéré et un second de voisinage permettent de retirer les cellules aberrantes.
  23. substract MNT sur la base du fichier MNT, la soustration de ce dernier permet de mettre le fichier dont il est extrait au relief Zéro, ceci est particulièrment intéressant pour le traitemnt de nuage de points sur des reliefs très accidentés (par ex. forêts en montagne).
  24. WG84_2_RFG93 convertit WGS84 (srid4326)  vers le CRS RFG93/Lambert 93 (srid2154),  calcule (rotation autour de 2 vecteurs avec point d’ancrage ou quaternions) la matrice de transformation pour caler un nuage de points source (nouvelle acquisition) sur un nuage cible de référence (calé nettoyé et traité) la méthode utilise la superposition de deux triangles semblabes le premier est situé dans le nuage cible (issu de relevé d’amers lors de l’acquistion des données sur le terrain), dans le nuage source il doit être aisé d’identifier le triangle homologue. Cette matrice (coordonnées homogènes) sera ensuite utilisée pour toutes transformations des données issues du fichier source. Les conditions de mesures changeant à chaque nouvelle acquisition de données (du moins entre deux recalages de la centrale inertielle)  cette procédure doit être reconduite.
  25. createPNG-X-Y-Zmoy, sur la base d’une extraction (parallélipipède)  d’un nuage d’une section horizontale ou verticale, crée par projection sur un plan  une image PNG normalisée en terme de géolocalisation en la calant sur la géolocalisation du nuage de points extrait,
  26. FilterCircleStandard Recherche dans une image PNG normalisée (reférentiel de localiation point/image calé sur son image PNG de réfrence /  attributs à introduire dan les données EXIF), fournit en sortie la liste des arbres détectés (pic de corréléation) pour un filtre de diamètre donné,
  27. neigbourHoodCluster Reprend les données en qualifiant définitivement les pics de corrélation en réduisant le taux de fausse alarme, par la méthode du pivot d'accumulation.
  28. benchMarkFilter les section horizontales des tronc d'arbre se présentent comme des enveloppes convexes, non nécessairement fermées et à densité de points variable, des précautions sont à prendre pour le paramétrage des filtres de corrélations et de détection des arbre. Cet outil permet d'ajuster ces paramères en appliquant également des traitements morphologiques  sur des panels de sections d'arbres collectée (banque d'image pour test), afin de normaliser leur densité. Il intègre également les fonctionnalités de neigbourHoodCluster, fournissant ainsi les coordonnées image (x, y en pixel) des enveloppes convexes (centre barycentrique) détectées pour un filtre données.

 

Construction SIG et Raster

  1. Raster des cellules Sintegra => /home/mbariou/PointCloud/LaMassaneCanopee/lidarLamassane.pdf, l’images est à l’écelle 1/10000 et les coordonnées latitudes longitudes sont données sur le pourtour de celle-ci afin d’en assurer la calibration lors de l’intégration SIG avec GeoReferencer Raster
  2. Fichier de points précalibrés du raster des cellules Sintegra est ici => /home/mbariou/PointCloud/LaMassaneCanopee/lidarLamassane.pdf.points
  3. Des fichiers utilitaires et/ ou intermédiaires dans le répertoire =>  /home/mbariou/PointCloud/LaMassaneCanopee/

Caractérisation du jeu de données

Relevé des empreintes au sol d'un jeu de données *.las => Obtention d'un fichier géométrique JSON, multipolygonale et duquel il faut extraire les polygônes principaux.

cde en ligne =>  pdal info --boundary Massane2.las > boundary_massane.json

Changement de format d'un jeu de données, conversion d'un fichier texte vers un format *.las

Urilisation de la commande translate appelant la fichier *.json qui suit.

{
  "pipeline":[
  {
   "type":"readers.text",
   "separator":" ",
   "filename":"12616_Massane_MNT50cm_L93-IGN69_000001.xyz"
  },
   {
   "type":"writers.las",
   "filename":"outputfile.las"
  }
  ]
}

Le fichier texte à transformer, 12616_Massane_MNT50cm_L93-IGN69_000001.xyz,  se présente comme suit :

X Y Z
702049.750 6155963.250 633.17
702049.750 6155962.750 633.35
702049.750 6155962.250 633.54
702049.750 6155961.750 633.72
702049.750 6155961.250 633.95
702049.750 6155960.750 634.15
702049.750 6155960.250 634.33
...

La ligne d'en-tête, avec séparateur (un espace comme indiqué dans le fichier de paramétrage *.json, translatetxt2las.json, donne les trois coordonnées présentes dans le fichier. Le fichier de sortie est un fichier *.las

La commande en ligne correspondante est :

pdal pipeline translatetxt2las.json

 A l'issue de dette commande nous disposonsons d''un fichier *.las, dont nous extrayons l'empreinte au sol

pdal info --boundary outputfile.las > boundary_massane-01.json

 

Données de test C++ Entrée de données map<string, struct> en fichier lecture et écriture

Rassemblés dans WKSCPLUPLUS/CodeProg/chp21 une série de miniprojets de test de nuage de points.

Un ostream convertit tout type d'objet en flux de bytes caractères (character) pour écriture en fichier

Un istream convertit un flux de bytes caractères (character) en tout type d'objet  pour lecture en fichier

un iostream assre les deux conversions précédentes selon les situations d'écriture ou de lecturetranslate

 

Construction de modèle numérique cellulaire de terrain .

Dans le cadre de travaux de foresterie, des modèles numériques de terrain issus de relevés laser aériens, collectent de l'ordre de 10 à 25 pts par m2 de terrain, après extraction sur les semis brut des nuages de points, des éléments classifiés comme dernier écho d'un tir laser aérien, le dernier écho donc le plus éloigné, l'hypothèse étant faite qu'il s'agit de l'écho du sol.

Par ailleurs, les résultats sont fournis sur des carroyages de 25 000 m2 (500m x 500m). En relief accidenté ou à variation notable de relevé, l'analyse d'un nuage de point permettant les extractions d'arbres et l'évaluation de leurs attributs dimensionnel est facilité par la remise de tous ces arbres à l'altitude 0. Ceci est obtenu en réduisant la coordonnées Z (altitude) de chaque point lidar, de relevé en sous-bois,  de la valeur de l'altitude du sol fourni par la cellule moyennée de MNT correspondant à la géolocalisation de la zone traitée un instant donné. Afin de limiter le volume de calcul, on s'appuie sur une altitude moyenne évaluée sur une cellule de 5m x 5m (25 m2) , cellule dans laquelle sont les points courants du nuage analysé. Il peut être envisagé de retenir des cellules de tailles différentes, cest une affaire de compromis entre le temps de calcul souhaitable, l'écart-type toléré par chaque cellule moyennée.

Dans les techniques de carroyages que nous proposons, dans le réfrentiel 2D XOY, chaque grille élémentaire est identifiée par la concaténation Xmin.Ymin (modulo le pas de grille retenu), en 3D pour un voxel l'identification se fait par la concaténation  Xmin.Ymin.Zmin (modulo le pas de grille retenu). Ainsi le quotient d'une opération modulo, permet une classification rapide d'un point d'écho soit dans sa gille soit dans son voxel de localisation.

Ceci est obtenu au prix de la gestion, des grilles, des voxels, des cellules moyennées de terrain, dans des conteneurs <map> non ordonnés et à accès aléatoire (ce qu'ils sont par définition).

 

Dans les différentes analyses que nous conduisons, nous travallons sur des bandes de terrains et nos méthodes de rangement des données permettent de manipuler une seule fois chque point que nous extrayons du nuage principal.

 

 

 

 

 

 

 

 

 

 

Nous utilisons des cookies. Ils sont activés sur ce sites/, êtes vous d'accord ? Des infos sur les Cookies => ICI informations privées.

  J'accepte les cookies de ce site.
EU Cookie Directive plugin by www.channeldigital.co.uk