Caminando hacia un DTM con Software Libre 1

Hola:

Esto no es un tutorial, son mis apuntes. Si alguien ve un error le agradecería me lo hiciera saber.

Aunque es sobre software libre, lastools no es completamente libre, pero eso si, lastools es bastante chido.

cmd_lasview_tiles_ground_rgb_filtered

¿Qué quiero? –> Quiero obtener un DTM (Modelo Digital de Terreno) a partir de un DSM (Modelo Digital de Superficie).
Tengo un DSM generado con Pix4d.
Y ya tengo el DTM que arrojó Pix4, pero en algunos lugares la elevación del DTM quedó hasta por 4 metros por debajo de la elevación de los Puntos de Control Terreste. Pero fuera de eso se ve muy bien este DTM de Pix4d.
Quise probar si había otra herramienta que no “rasarura las elevaciones en la cúspide del cerro. Y por eso usé Lastools, que es bueno, rápido, elegante, chingón!!!
Los resultados con lastools fueron mejores que con Pix4d en unas partes pero en otras Pix4d lo hace mejor. Estoy seguro que mucho de todo esto se debe a que la fuente de los datos una nube de puntos resultado de visión computacional y no un LIDAR, y a partír de ahí el algoritmo que usan para encontrar e interpolar el suelo desnudo.

También podrían limpiarse a mano árboles y construcciones. El Pix4d lo permite, claro que hay que pagar el costo de la licencia. Y en lastools habría que reclasificar los puntos que no son considerados suelo desnudo, pero para que un cambio de clasificación surta efecto hay que pagar licencia, y ésta oscila de 1,000 a 5,000 Euros por año. Cabe decir que aquí en México, para las chambas que tiene un despacho de Topografía pequeño el costo de una licencia de este monto está fuera de contexto.
Y esta es la razón principal y el comparar los resutados entre dos paquetes. Esto tiene la ventaja de que se aprende mucho en el camino, y la enorme desventaja de poder frustrarse por lo lento que se avanza.

Manos a la obra:

Leí los siguientes enlaces para entender cómo hacerlo:
1.- https://rapidlasso.com/2014/06/13/dtm-from-dense-matching-points/
2.- http://canalgeo.com/obtener-dtm-de-archivo-lidar-fotogrametria/
3.- https://rapidlasso.com/2016/11/22/creating-a-better-dtm-from-photogrammetic-points-by-avoiding-shadows/
4.- https://rapidlasso.com/2014/08/25/processing-large-dense-matching-point-clouds/
5.- https://rapidlasso.com/2015/09/21/creating-dtms-from-dense-matched-points-of-uav-imagery-from-senseflys-ebee/

Qué haré:
1.- Copiar el DSM a un lugar cercano a los binarios de lastools.
Tengo 2 archivos:
1) terreno_dsm.laz
2) terreno_dsm_200cm.laz (que es un muestreo a cada 2 metros, una retícula) NOTA: No tengo claro si este muestreo contiene los puntos de elevación más baja, lo cual sería de mucha ayuda para aplicar PDAL.
2.- Voy a obtener información de ellos usando los comandos, a través del cmd: de Windows:
C:\LAStools\bin>lasinfo.exe ..\dsm_to_dtm\terreno_dsm.laz
lasinfo (161114) report for ..\dsm_to_dtm\terreno_dsm.laz
reporting all LAS header entries:



number of first returns: 85608038
number of intermediate returns: 0
number of last returns: 85608038
number of single returns: 85608038
WARNING: there are 85608038 points with return number 0
WARNING: there are 85608038 points with a number of returns of given pulse of 0
histogram of classification of points:
85608038 never classified (0)

Tengo 85608038 en un solo retorno, con la misma intensidad, no es el resultado de un levantamineto con LIDAR, sino de un levantamiento con drone con cámara digitial pues.
Ahora veré la información del otro archivo:
C:\LAStools\bin>lasinfo.exe ..\dsm_to_dtm\terreno_dsm_200cm.laz
lasinfo (161114) report for ..\dsm_to_dtm\terreno_dsm_200cm.laz



intensity 0 0
return_number 0 0
number_of_returns 0 0
edge_of_flight_line 0 0
scan_direction_flag 0 0
classification 0 0
scan_angle_rank 0 0
user_data 0 0
point_source_ID 0 0
gps_time 0.000000 0.000000
Color R 4096 65280
G 4352 65280
B 4864 65280
number of first returns: 166874
number of intermediate returns: 0
number of last returns: 166874
number of single returns: 166874
WARNING: there are 166874 points with return number 0
WARNING: there are 166874 points with a number of returns of given pulse of 0
histogram of classification of points:
166874 never classified (0)

Tengo mucho menos puntos, podría usar este archivo y obtener a partir de él el MDT, pero no estoy seguro que el muestreo de su retícula lo haya hecho a partir del punto de elevación más bajo encontrado, que podría acercarse a la elevación del suelo desnudo.
Por otro lado, también hay que notar que en ambos archivos no está considerando el tiempo gps (gps_time 0.000000 0.000000), y me parece que es porque no estoy usando el lastools con licencia.

3.- Partir el archivo en pequeñas teselas de 50 x 50 metros, con un buffer de 10 metros, con la opción “flag_as_withheld”, entiendo que para que sea fácil posteriormente quitar el bufer y renombrar la salida en un directorio llamado “tiles_raw” con el sufijo “_tiled” para saber que cosa tienen esos archivos.

C:\LAStools\bin>mkdir ..\dsm_to_dtm\tiles_raw

C:\LAStools\bin>lastile.exe -i ..\dsm_to_dtm\terreno_dsm.laz -tile_size 50 -buffer 10 -flag_as_withheld -epsg 32614 -o ..\dsm_to_dtm\tiles_raw\terreno_dsm_tiled.laz
Please note that LAStools is not "free" (see http://lastools.org/LICENSE.txt)
contact 'martin.isenburg@rapidlasso.com' to clarify licensing terms if needed.

Después de un ratito obtengo varios archivos “.laz” que son las teselas de todo el “DMS”. En este caso fueron 316 archivos. En cuestión de un par de minutos. Por eso sigo pensando que lastools es una chulada.

4.- Reordenaré los puntos a través de un rellenado espacial con curvas de Hibelert. Utilizando los 4 núcleos del procesador para darle poder!!! Le decimos al comando que use el sufijo “_sorted” y que los archivos de salida tengan el formato “.laz”.

C:\LAStools\bin>lassort.exe -i ..\dsm_to_dtm\tiles_raw\*.laz -odir ..\dsm_to_dtm\tiles_sorted\ -olaz -odix _sorted -cores 4

5.- Luego uso lasthin. Según el enlace (3), de arriba, con lasthin disminuimos el número de puntos a considerar como puntos sobre la tierra. En este caso yo obtendré los puntos de menor elevación en una retícula de 1 metroX1metro y los clasifico con el código “9”.

C:\LAStools\bin>lasthin -i ..\dsm_to_dtm\tiles_sorted\*.laz -step 1.0 -lowest -classify_as 9 -odir ..\dsm_to_dtm\tiled_thinned\ -odix _thinned -olaz -cores 4

Según el post (3) debe haber una reconsideración: hay puntos que están en la sombra cuyo valoren RGB es cercano a cero, y considerando que la imagen es de 16 bits. Entonces hay que usar las2las para reclasificar todos los puntos cuyo color en RBG están cercanos al cero y asignarles el código 7. Para el red entre 0 y 30, green 0 y 35, blue 0 y 40. Pero como es una imagen de 16 bits, hay que multiplicar estos números por 256.

C:\LAStools\bin>las2las.exe -i ..\dsm_to_dtm\tiles_thinned\*.laz -keep_RGB_red 0 7680 -keep_RGB_green 0 8960 -keep_RGB_blue 0 10240 -filtered_transform -set_classification 7 -odir ..\dsm_to_dtm\tiles_rgb_filtered\ -odix _rgb_filtered -olaz -cores 4

6.- Ahora hay que correr lasground para obtener el terreno pero ignorando los puntos clasificados como 0 y como 7. Debo tener en cuenta que hay puntos clasificados con el código 9, producto de comando lasthin que use anteriormente. En el enlace (1) el las sentencias de la línea de comandos muestra que el lasgroud debe correrese para los archivos que están en \tiles_thinned\ lo cual no me parece que esté bien, porque lo que se acaba de hacer es clasificar los puntos con determinado RGB y ponerlos en \tiles_rgb_filtered\, así que yo lo haré desde ahí:

C:\LAStools\bin>lasground.exe -i ..\dsm_to_dtm\tiles_rgb_filtered\*.laz -step 20 -bulge 0.5 -spike 0.5 -fine -ignore_class 0 7 -odir ..\dsm_to_dtm\tiles_ground_rgb_filtered\ -olaz -odix _ground -cores 4

¿Cómo se ve?
lasview_tiles_ground_rgb_filtered

7.- Finalmente hay que unir todas las teselas o tiles en un gran archivo (merge), para eso se debe correr el blast2dem y/o lasgrid. Lasgrid genera un raster de nxn metros a partir de un archivo “laz”. Y blas2dem crea un TIN (Red de Tríangulos Irregulares) y luego genera un raster a partir de él que resulta en un DEM (Modelo Digital de Elevaciones). Más parecido a lo que nosotros los topos hacemos con CivilCAD de ArqCOM. ¿No?

Primero corro el lasgrid, y obvio, como no tengo licencia el resultado no es el esperado, tal como lo indica en la línea de comandos:
lasgrid.exe -i ..\dsm_to_dtm\tiles_ground_rgb_filtered\*.laz -merged -drop_withheld -step 0.5 -average -o ..\dsm_to_dtm\terreno_dsm.bil

cmd_lasgrid_tiles_ground_rgb_filtered
cmd_lasgrid_tiles_ground_rgb_filtered

Y efectivamente le inserta una línea negra y por si fuera poco… digaonal. ¡Qué cosas!
lasview_lasgrid_tiles_ground_rgb_filtered

Ok. Ahora blast2dem, en realidad no entiendo porqué tiene el parámetro -keep_class 2 , me parece que debiera ser -keep_class_ 9. Porque 9 es el suelo desnudo, o sea, ground. Probaré ambas opciones. Esto lo sugiere igualmente el post (3).
C:\LAStools\bin>blast2dem.exe -i ..\dsm_to_dtm\tiles_ground_rgb_filtered\*.laz -merged -drop_withheld -keep_class 2 -step 0.5 -o ..\dsm_to_dtm\terreno_dtm.bil
Please note that LAStools is not "free" (see http://lastools.org/LICENSE.txt)
contact 'martin.isenburg@rapidlasso.com' to clarify licensing terms if needed.
WARNING: first file is a buffered tile. maybe remove buffers first?
WARNING: first file is a buffered tile. maybe remove buffers first?

C:\LAStools\bin>blast2dem.exe -i ..\dsm_to_dtm\tiles_ground_rgb_filtered\*.laz -merged -drop_withheld -keep_class 9 -step 0.5 -o ..\dsm_to_dtm\terreno_dtm_9.bil
Please note that LAStools is not "free" (see http://lastools.org/LICENSE.txt)
contact 'martin.isenburg@rapidlasso.com' to clarify licensing terms if needed.
WARNING: first file is a buffered tile. maybe remove buffers first?
WARNING: first file is a buffered tile. maybe remove buffers first?
ERROR: no input points
ERROR: wrong reader (data is -1 but reader is SPB 7)
ERROR: wrong SMreader (need -1 but this is SMreader_smb 0)
check ... maybe the file (or the stream) is empty

Ok, ahora leyendo: “https://rapidlasso.com/lastools/lasground/”, me doy cuenta que lasground clasifica los puntos “ground” como class = 2 y los “non-ground” como class = 1.

Y MUCHO OJO, TENER PRESENTE: que hay que probar lo siguientes, según el README de lasground (http://www.cs.unc.edu/~isenburg/lastools/download/lasground_README.txt):

The tool also produces excellent results for town or cities
  but buildings larger than the step size can be problematic.
  The default step size is 5 meters, which is good for forest
  or mountains. For towns or flat terrains '-town' the step
  size is increased to 10 meters. For cities or warehouses
  '-city' the step size is increased to 25 meters. For very
  large cities use '-metro' and the step size is increased
  to 50 meters You can also set it directly with '-step 35'.

Así que yo utilicé en un principio -step 20, ahora voy a usar -step 5. Y también leo que para puntos obtenidos con fotogrametría por Pix4d y Agisoft conviene usar el parámetro -bulge que según mi la inglés debería tener un valor de la décima parte del step cuando éste vale 5. Así que serí algo así como -bulge 0.5. En resulmen quedaría así:

lasground.exe -i ..\dsm_to_dtm\tiles_rgb_filtered\*.laz -step 5 -bulge 0.5 -spike 0.5 -fine -ignore_class 0 7 -odir ..\dsm_to_dtm\tiles_ground_rgb_filtered_step5\ -olaz -odix _ground_step5 -cores 4
Si compraro los resultados, la imagen de la izquierda en con el -step 5 y la derecha con -step 20. Algo cambia, pero creo que podría apreciarse mejor teniendo la ortofoto de fondo y curvas de nivel de ambos caso.

lasgrond_step5_izq_vs_step_20_der
lasgrond_step5_izq_vs_step_20_der

Ahora faltaría hacer un blast2dem a los archivos resultado del lasground con -step 5 y luego mostrar todo traslapado.

C:\LAStools\bin>blast2dem.exe -i ..\dsm_to_dtm\tiles_ground_rgb_filtered_step5\*.laz -merged -drop_withheld -keep_class 2 -step 0.5 -o ..\dsm_to_dtm\terreno_dtm_step5.bil
Please note that LAStools is not "free" (see http://lastools.org/LICENSE.txt)
contact 'martin.isenburg@rapidlasso.com' to clarify licensing terms if needed.
WARNING: first file is a buffered tile. maybe remove buffers first?
WARNING: first file is a buffered tile. maybe remove buffers first?

Sólo falta pasar las curvas que decida que están mejor utilizando AutoCAD Map, claro que sería más elegante con gdal. Pero no me sale esto de pasar con gdal de SHP a DXF. Así que esto luego.
También se me ocurre que pueden quitarse las construcciones haciendo una especie de CROP o en términos de SIG un geoproceso llamado CLIP o RECORTE, y luego una interpolación. Y al final inciar a todo este rollo con lastools teniendo resuelto quitar las construcciones desde un inicio. Se puede hacer con postgis, creo que también con el lastools y seguro con qgis o grass y con Arcmap.

¿Es mucho rollo y muy tardado? Bueno, teniendo todo claro y en un script sería mucho más rápido que dejar al Pix4d trabajando por horas y consumiéndose todos los recursos de la máquina. ¿No? Peeero, al final no salió tan mal con Pix4d.

COMPARACION DE RESULTADOS:

En la siguientes imagenes se ven las curvas de nivel en rojo con el -step 5 y en azul con -step 20 y en amarillo las de Pix4d.
Aquí cuado todas las curvas se ven igual de bien:

comparacion_curvas_pix3d_vs_lastools_step5y20-1
comparacion_curvas_pix3d_vs_lastools_step5y20-1

Aquí donde la interpolación de Pix4d parece que rasura la cúspide del monte:

comparacion_curvas_pix3d_vs_lastools_step5y20-2
comparacion_curvas_pix3d_vs_lastools_step5y20-2

Y aquí cuando Pix4d parece que interpola mejor que lastools:

comparacion_curvas_pix3d_vs_lastools_step5y20-3
comparacion_curvas_pix3d_vs_lastools_step5y20-3

Lo que concluyo son dos cosas:
1.- Hay que usar Filtros Morfológicos Progresivos, no se que algoritmo usen Pix4d y lastools.
2.- Hay que usar LIDAR. Esos árboles tupidos no dejan ver el suelo desnudo para nada.

Nota: Recuerdo que en un texto que habla de usan kinetic de microsoft como sensor lidar, me parece que con una resolución espacial de unos 35 cm. Aunque no sé el número de retornos. Aquí se puede ver: http://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XLI-B7/945/2016/