{"id":201,"date":"2018-11-12T16:30:33","date_gmt":"2018-11-12T15:30:33","guid":{"rendered":"https:\/\/sigeo.cerege.fr\/?p=201"},"modified":"2018-11-13T14:14:22","modified_gmt":"2018-11-13T13:14:22","slug":"lidar-filtrage-du-sol-et-creation-de-mnt","status":"publish","type":"post","link":"https:\/\/sigeo.cerege.fr\/?p=201","title":{"rendered":"Lidar : Filtrage du sol et Cr\u00e9ation de MNT"},"content":{"rendered":"

Workflow pour filtrer les points du sol dans des donn\u00e9es Lidar et cr\u00e9er un MNT<\/h1>\n

Ce workflow est un exemple de cas d’\u00e9tude que nous avons eu \u00e0 traiter, il est \u00e0 adapter en fonction des donn\u00e9es utilis\u00e9es<\/p>\n

Nous allons utiliser LibLAS, PDAL et Points2Grid<\/p>\n

Donn\u00e9es en entr\u00e9e<\/h2>\n

Nous avons des nuages de points Lidar au format LAS, avec des coordonn\u00e9es en WGS84 pour le X et Y, et en m\u00e8tres pour le Z<\/p>\n

Affichage des m\u00e9tadonn\u00e9es (avec libLAS):<\/p>\n

lasinfo input.las<\/code><\/p>\n

Affichage des m\u00e9tadonn\u00e9es (avecPDAL):<\/p>\n

$ pdal info input.las<\/code><\/p>\n

il faut v\u00e9rifier si le syst\u00e8me de coordonn\u00e9es est d\u00e9finit, s’il y a un facteur d’\u00e9chelle selon X, Y ou Z, et toutes autres informations<\/p>\n

Supposons que le syst\u00e8me de coordonn\u00e9es n’est pas d\u00e9finit ou mal d\u00e9finit dans les m\u00e9tadonn\u00e9es du LAS, cela va causer des probl\u00e8mes dans la suite du traitement<\/p>\n

D\u00e9finition du syst\u00e8me de coordonn\u00e9es<\/h2>\n

Cr\u00e9er un fichier JSON pour le traitement PDAL, nomm\u00e9 par exemple “assign_scr.json” :<\/p>\n

{<\/span>\r\n    \"pipeline\"<\/span>:<\/span> [<\/span>\r\n        {<\/span>\r\n            \"type\"<\/span> :<\/span> \"readers.las\"<\/span>,<\/span>\r\n            \"filename\"<\/span> :<\/span> \"input.las\"<\/span>\r\n        },<\/span>\r\n        {<\/span>\r\n            \"type\"<\/span> :<\/span> \"writers.las\"<\/span>,<\/span>\r\n            \"a_srs\"<\/span>:<\/span> \"EPSG:4326\"<\/span>,<\/span>\r\n<\/span>            \"filename\"<\/span> :<\/span> \"output.las\"<\/span>\r\n        }<\/span>\r\n    ]<\/span>\r\n}\r\n<\/span><\/pre>\n
Puis lancer la commande:<\/p>\n
pdal pipeline assign_scr.json<\/code><\/p>\n
Maintenant la d\u00e9finition est correcte, mais il est pr\u00e9f\u00e9rable pour les interpolations et les traitements sur les MNT d’\u00eatre en projection avec les m\u00eames unit\u00e9s en horizontal qu’en vertical, nous allons donc reprojeter en UTM<\/p>\n
Reprojection<\/h2>\n
cr\u00e9er un pipeline JSON comme suit, nomm\u00e9 par exemple “reproject.json” :<\/p>\n
{<\/span>\r\n  \"pipeline\"<\/span>:[<\/span>\r\n    {<\/span> \r\n        \"type\"<\/span>:<\/span>\"writers.las\"<\/span>,<\/span> \r\n        \"filename\"<\/span>:<\/span>\"input.las\"<\/span> \r\n    }<\/span>\r\n    {<\/span>\r\n        \"type\"<\/span>:<\/span>\"filters.reprojection<\/span>\"<\/span>,<\/span>\r\n        \"out_srs\"<\/span>:<\/span>\"EPSG:32633\"<\/span>\r\n    },<\/span>\r\n    {<\/span>\r\n      \"type\"<\/span>:<\/span>\"writers.las\"<\/span>,<\/span>\r\n      \"filename\"<\/span>:<\/span>\"output.las\"<\/span>\r\n    }<\/span>\r\n  ]<\/span>\r\n}<\/span><\/pre>\n
Puis lancer la commande:<\/p>\n
pdal pipeline reproject.json<\/code><\/p>\n
Il faut ensuite faire une classification sur les points<\/p>\n
Classification du sol<\/h2>\n
cr\u00e9er un pipeline JSON comme suit, nomm\u00e9 par exemple “ground_classif.json” :<\/p>\n
{\r\n  \"pipeline\":[\r\n    {\r\n         \"type\":\"filters.assign\",\r\n         \"assignment\":\"Classification[:]=0\"\r\n    },\r\n    {\r\n         \"type\":\"filters.elm\"\r\n    },\r\n   {\r\n         \"type\":\"filters.outlier\",\r\n         \"mean_k\":8,\r\n         \"multiplier\":3 \r\n   },\r\n   {\r\n         \"type\":\"filters.smrf\",\r\n         \"ignore\":\"Classification[7:7]\",\r\n         \"slope\":0.2,\r\n         \"window\":16,\r\n         \"threshold\":0.45,\r\n         \"scalar\":1.2\r\n    },\r\n    {\r\n         \"type\":\"filters.range\",\r\n         \"limits\":\"Classification[2:2]\"\r\n     }\r\n  ]\r\n}<\/pre>\n
Ce pipeline est bien d\u00e9crit dans les aides de PDAL, les \u00e9tapes en sont les suivantes :<\/p>\n
    \n
  1. r\u00e9assigner tous les points \u00e0 la classe 0 car on suppose que la classification d’origine n’est pas bonne<\/li>\n
  2. utilisation de l’algo ELM pour assigner les points consid\u00e9r\u00e9s comme du bruit \u00e0 la classe 7<\/li>\n
  3. faire un filtrage statistique pour les points abb\u00e9rants<\/li>\n
  4. utiliser le filtre SMRF sur les points restants pour classer les points<\/li>\n
  5. ne r\u00e9cup\u00e9rer que les points class\u00e9s comme ground<\/li>\n<\/ol>\n
    Puis lancer la commande:<\/p>\n
    pdal pipeline ground_classif.json<\/code><\/p>\n
    on r\u00e9cup\u00e8re donc un fichier las ne contenant que les points class\u00e9s comme sol (classe 2)<\/p>\n
    il est maintenant possible d’interpoler les points pour cr\u00e9er un MNT<\/p>\n
    G\u00e9n\u00e9ration du MNT<\/h2>\n
    l’interpolation peut se faire avec PDAL ou avec Points2Grid (outil utilis\u00e9 par OpenTopography pour g\u00e9n\u00e9rer des MNT). L’impl\u00e9mentation de PDAL est bas\u00e9e sur celle de Points2Grid, donc au final c’est \u00e9quivalent je pense.<\/p>\n
    La version de points2grid que nous utilisons est celle de d\u00e9veloppement ici<\/a><\/p>\n
    points2grid -i ground.las -o dtm.tif -r 1 --output_format grid --input_format las --interpolation_mode incore --idw --resolution 1 --fill --fill_window_size 5<\/code><\/p>\n
    pour des explications sur les param\u00e8tres taper : points2grid -h<\/code><\/p>\n
    R\u00e9f\u00e9rences:<\/h2>\n