Importation rapide d'un fichier las dans une table de Postgis et analyse de densité.

Je viens d'être confronté à un besoin d'estimation rapide de la densité de points par m2 d'un fihier "las" provenant d'un vol Lidar. Il fallait pouvoir répondre à la question combien de zones de 1m2 contiennnent des densités de points inférieur à un certain seuil.
Wikipedia l'explique bien mieux mais en résumé pour les pressés : La télédétection par laser ou LIDAR « light detection and ranging » est une technique qui permet en substance de "mesurer" depuis un avion ou un hélicoptère la position d'objets aux sol et d'obtenir un nuage de points x,y,z  qui donne un modèle numérique de surface ou MNS.
En fonction du plan vol, de la météo des conditions topologiques en présence (grand bâtiments) on peut avoir des "zones d'ombres" ou on a peu ou pas de  données, il faut néanmoins pouvoir quantifier et visualiser ces zones.
Pour les curieux je recommande l'excellente présentation de M. Marc Riedo responsable du SITN sur les applications concrètes du Lidar dans le canton de Neuchâtel.


J'étais curieux de voir les informations  que je pouvais obtenir rapidement  depuis Linux avec un échantillon de données au format  ASPRS LAS format fourni par notre prestataire de service.

M. Riedo m'avait d'ailleurs déjà rendu attentif à l'existence de la librairie opensource liblas . Vous pouvez l'installer sur votre distribution Linux via votre gestionnaire de paquets ou en ligne de commande sur une debian et autre dérivés (*ubuntu) en tappant la commande suivante : 
sudo apt-get install  liblas-bin liblas-dev python-liblas liblas1
Dans le paquet liblas-bin il y a toute une série d'utilitaires dont "lasinfo" dont vous pouvez tout de suite tirer parti pour avoir toute une série d'informations synthétique sur votre fichier las min,max des différentes valeurs, nbre points par classe etc.:
lasinfo swissphoto_1164212_Auschnitt.las
...
---------------------------------------------------------
  Point Inspection Summary
---------------------------------------------------------
  Header Point Count: 6862236
  Actual Point Count: 6862236

  Minimum and Maximum Attributes (min,max)
---------------------------------------------------------
  Min X, Y, Z:         561535.84, 204627.33, 426.14
  Max X, Y, Z:         562535.83, 205627.32, 595.77
  Bounding Box:        561535.84, 204627.33, 562535.83, 205627.32
  Time:            0.000000, 318694.920374
  Return Number:    1, 4
  Return Count:        1, 4
  Flightline Edge:    0, 1
  Intensity:        0, 5100
  Scan Direction Flag:    0, 1
  Scan Angle Rank:    -20, 19
  Classification:    1, 10
  Point Source Id:    71, 25524
  User Data:        0, 0
  Minimum Color (RGB):    0 0 0
  Maximum Color (RGB):    0 0 0

  Number of Points by Return
---------------------------------------------------------
    (1) 6791699    (2) 59895    (3) 5697    (4) 4945

  Number of Returns by Pulse
---------------------------------------------------------
    (1) 6731951    (2) 108308    (3) 2527    (4) 19450

  Point Classifications
---------------------------------------------------------
    2777140 Unclassified (1)
    534520 Ground (2)
    1292256 Low Vegetation (3)
    2103334 Medium Vegetation (4)
    6844 High Vegetation (5)
    20293 Building (6)
    1033 Low Point (noise) (7)
    122078 Model Key-point (mass point) (8)
    4738 Reserved for ASPRS Definition (10)
  -------------------------------------------------------
      0 withheld
      0 keypoint
      0 synthetic
  -------------------------------------------------------
Donc ici nous avons 6.8 millions de points dont il s'agirait d'analyser par exemple la densité par m2 (sans utiliser un fme...) avec des outils opensources

QGIS avec les points du LAS après import dans postgis

J'ai choisi d'importer un certain nombre d'attributs du format las dans une table postgis pour pouvoir ensuite triturer ces informations via le couple QGIS & Postgis pour visualiser les zones de 1m2 qui ne contiennent pas de points (carrés jaune) ou trés peu de points (entre 1 et 2 pts/m2 carré vert ci-dessous) tout en affichant une couleur différente  pour les point_source id.

Le temps d'affichage est excellent grâce à l'index spatial Gist de postgis
Bien entendu ceci est vrai tant que l'on essaye pas d'afficher un zoom sur l'étendue de tout la zone car avec tous les points ça sera un peu long et pas vraiment très utile(6.8 millions de features afficher c'est rarement une bonne
idée...)
Voici comment  j'ai procédé dans ce cas :

1) Génération d'un fichier asci de type csv (ou séparé par ce que vous voulez) à partir du fichier las: on peut le faire avec l'utilitaire las2txt
par exemple avec une commande comme :
las2txt --parse xyztrainrcupedRGBM -i swissphoto_1164212_Auschnitt.las -o swissphoto_1164212_Auschnitt_zone_test.txt  --minx 561985 --miny 205284 --maxx 561988  --maxy 205287
on peut instantanément  extraire tous les points avec certains attributs sur une zone d’intérêt uniquement 
head -n 3 swissphoto_1164212_Auschnitt_zone_test.txt 
561985.09,205285.86,479.36,224702.63878056,1,-7,37,1,1,1,0,25523,0,0,0,0,0,0
561985.29,205285.38,479.33,224702.63878856,1,-7,31,1,1,1,0,25523,0,0,0,0,0,1
561985.60,205285.29,479.21,224702.64446056,1,-7,30,1,1,1,0,25523,0,1,0,0,0,2
...
Seuleument dans mon cas je voulais aussi importer l'attribut  gps-time en tant  que datetime dans postgis  or dans le cas de las2txt il s'affiche comme un nombre float double,  il y a probablement moyen de convertir mais j'ai décidé pour am culture de tester l'interface python de la librairie liblas qui est vraiment complète et bien documentée. Ce qui fait qu'avec ces lignes de python ci-après, j'ai pu générer mon fichier texte avec les informations que je voulais et le champ gps-time au format datetime




#!/usr/bin/python
from liblas import file
import os.path
lasfile2import = sys.argv[1];
if (os.path.isfile(lasfile2import)):
    f = file.File('swissphoto_1164212_Auschnitt.las', mode = 'r')
    sep = ";"
    i=1
    for p in f:  
        print str(i) + sep + '1' + sep + str(p.x) +  sep +  str(p.y) +  sep + str(p.z) +  sep + str(p.get_classification()) +  sep + str(p.get_intensity()) +  sep + str(p.get_scan_angle_rank()) + sep + str(p.get_point_source_id()) + sep + str(p.time)
        i+=1
 
    f.close()
else:
    print "Le premier parametre n'est pas un fichier : " + lasfile2import

ce fichier python permet de générer un fichier de données texte séparé par un point-virgule.
python cmdReadLasFile.py swissphoto_1164212_Auschnitt.las > swissphoto_1164212_Auschnitt.csv

SAGA est un logiciel de GIS Libre

Ah au fait si comme moi vous avez installé liblas manuellement en compilant le source de la dernière version (1.7.0 ces jours) et bien vous aurez quelques fonctionnalités intéressantes en plus par rapport à la version disponible en standard (liblas1 1.2.1-1) sur ma distribution Ubuntu 10.04.4 LTS mais cette installation manuelle d'une liblas plus récente  risque de nuire à la bonne marche de l'importation de fichier las dans l'excellent produit saga-gis. En résumé si vous voulez compiler la version stable 2.0.8 de saga-gis il faut pour l'instant prendre garde à utiliser la version 1.2.1 de liblas et ne pas déployer une version plus récente dans les chemins par défaut accessibles pour les librairies (/etc/ld.so.conf)
D'ailleurs pour générer des images de densité de points saga-gis est excellent comme expliqué sur ce blog : http://dominoc925.blogspot.com/2010/08/create-lidar-point-density-map-with.html

Bon revenons à nos moutons... Maintenant que j'ai un fichier texte avec les valeurs qui m'intéressent de mon fichier las, je vais pouvoir l'importer dans une table postgresql avec la commande "COPY"
au préalable il faut bien entendu avoir créé cette table dans postgres voici ma définition de table qu'il faudra ajuster en fonction des attributs du fichier las que vous avez décidé d'importer

CREATE TABLE lidar_point
(
  id integer NOT NULL DEFAULT 0,
  iddatafile integer NOT NULL DEFAULT 1,
  x double precision NOT NULL DEFAULT 0.0,
  y double precision NOT NULL DEFAULT 0.0,
  z double precision NOT NULL DEFAULT 0.0,
  classification integer NOT NULL DEFAULT 0,
  intensity integer DEFAULT 0,
  scan_angle_rank integer DEFAULT 0,
  point_source_id integer DEFAULT 0,
  gps_time timestamp without time zone,
 CONSTRAINT pk_lidar_point PRIMARY KEY (id)
)
WITH (
  OIDS=FALSE
);

Ensuite il faut procéder à l'insertion et pour des questions de performances il faut impérativement utiliser cette variante de chargement de données ("bulk" upload) au lieu de faire des INSERT successif.

COPY lidar_point FROM 'swissphoto_1164212_Auschnitt.csv'
     USING DELIMITERS ';' CSV;

vous pouvez stocker cette commande sql dans un fichier et l'executer avec psql (avec -f monfichier.sql) ou alors la lancer directement :
psql -c "COPY lidar_point FROM 'swissphoto_1164212_Auschnitt.csv' USING DELIMITERS ';' CSV;" dblidartest
une fois que vos données ont été chargées vous pouvez rajouter une géométrie pour le point en 2d

SELECT AddGeometryColumn('','lidar_point','the_geom',21781,'POINT',2);

puis vous pouvez enchaîner en mettant a jour ce nouveau champ géométrique, ici dans le système de projection Suisse

UPDATE lidar_point SET  the_geom=ST_GeomFromText('POINT(' || x || ' ' || y || ')',21781);

et finalement vous pouvez lancer la création de l'index spatial gist sur cette table

CREATE INDEX idxgist_lidar_point2d on lidar_point USING GIST (the_geom);

Selon le nombre de point que vous avez insérer il faut maintenant patienter un peu... Non mais sérieusement vous pouvez envisager une autre activité en attendant que cela se fasse...
Une fois que c'est terminé vous pouvez vérifier l'emprise de tous les points avec un

SELECT ST_Extent(the_geom) FROM lidar_point;

qui vous retournera qqchose du style : BOX(561535.84 204627.33,562535.83 205627.32)
Ces informations vous permettrons le cas échéant de générer votre grille, dans mon cas c'est déjà fait et je n'ai plus qu'a lancer un :

UPDATE lidar_grid1m g SET num_points = (SELECT COUNT(*) FROM lidar_point p WHERE ST_Within(p.the_geom,g.the_geom)) WHERE g.the_geom && ST_SetSRID('BOX3D(561535 204627,562536 205628)'::box3d,21781);  

pour alimenter un attribut de ma table grille avec le nombre de points du lidar trouvé pour chaque carré de ma grille, je limite ici l'update de ma grille à la bounding box correspondant aux points de mon fichier las, cette clause where étant surtout utile dans le cas ou la table grille (ici lidar_grid1m) est beaucoup plus grande que la zone de points lidar.