¿Qué es laspy?
Laspy es una librería de Python que nos permite trabajar con archivos LiDAR, leyéndolos, modificándolos y creando nuevos archivos.
Esta librería funciona con versiones de Python 3.6 o superiores. En la web que recoge la mayoría de los paquetes y librerías de Python, podemos encontrar mas información sobre laspy y algunos ejemplos.
Preparando un entorno de trabajo
Vamos a ver algunos ejemplos sobre esta librería. Para realizar este ejercicio por tu cuenta solo necesitarás:
- Una instalación de Python 3.6 o superior
- Un entorno de trabajo de Python o «environment» donde instalaremos laspy
- Datos LiDAR
En este caso, vamos a usar la versión de Python 3.7.8 que puedes descargar desde la página oficial de Python.
Lo primero que vamos a hacer es crear nuestro «Python Environment» así, podemos inslatar la librería de laspy y evitar posibles conflictos entre librerías.
Una vez tenemos nuestro entorno de trabajo listo, es hora de activarlo. Para ello ejecuta el siguiente comando que básicamente activa nuestro «environment»:
env_laspy\Scripts\activate
Activando el espacio de trabajo nos permitirá instalar las librerías que queramos en él. Instala el paquete de «laspy», por defecto se instalará la versión mas reciente. Si no tenías «numpy» instalado en tu PC este también se instalará.
Ahora ya tenemos todo listo. Es el momento de crear un script de Python y leer el primer archivo LiDAR.
Leer un archivo LiDAR con LASPY
Para la realización de este ejemplo, vamos a usar Pycharm, pero podrías utilizar cualquier otra IDE para trabajar con Python, incluso el bloc de notas. Para descargar la versión gratuita de Pycharm pincha aquí.
Abre un proyecto de Pycharm en la carpeta donde tienes tu «environment» de Python, de este modo, Pycharm encontrar automáticamente el intérprete de Python que necesita. De forma alternativa, puedes manualmente establecer el intérprete tu mismo desde «settings»:
Ahora, crea un archivo de python desde la ventana de proyectos y llamalo «test_laspy»:
En primer lugar, vamos a comenzar importando laspy y numpy. Usaremos ambas librerías para nuestro propósito:
import laspy as lp
import numpy as np
El siguiente paso será añadir nuestros datos de entrada, es decir, la nube de puntos que usaremos como test:
input_path = "C:\Javi\Movidas\laspy\PNOA_1.las"
Ahora ya tenemos definida la variable con el archivo de entrada. Vamos a leer el archivo LiDAR y obtener el número de puntos que hay según la información de cabecera («header information»)
with lp.open(input_path) as fh:
print('Numero de puntos en la cabecera:', fh.header.point_count)
Con la función «print», obtendremos en pantalla el número de puntos o retornos que contiene nuestro archivo LiDAR, para la nube de puntos que estoy usando del PNOA, he obtenido algo más de 3 millones de puntos: 3.221.429.
Ahora, en lugar de leer el número de retornos desde la información de cabecera, vamos a leer el número total de puntos que tiene nuestro LiDAR. Esto es una práctica habitual y puede ayudarnos a identificar archivos corruptos cuando hay discrepancias entre ambos valores. Para ello, ejecuta el siguiente código:
with lp.open(input_path) as fh:
print('Numero de puntos en la cabecera:', fh.header.point_count)
las = fh.read()
print('Points from data:', len(las.points))
Ahora tienes dos funciones «print» y ambas te deberían dar el mismo resultado en tus datos, que en este caso es 3221429. Uno de los valores se consigue usando «point_count» y el otro con la función «len()».
Filtrar una nube de puntos por clase y retorno
Cuando tenemos nubes de puntos clasificadas, es interesante conocer formas de filtrar el LiDAR para obtener información sobre cuantos puntos tenemos en cada clase, así como la distribución del número de retornos. En cualquier nube de puntos clasificada, deberíamos esperar que el suelo tuviera valores muy altos para el primer (y único retorno). Vamos a ver si esto se cumple en nuestra nube de puntos:
with lp.open(input_path) as fh:
las = fh.read()
ground_pts = las.classification == 2
bins, counts = np.unique(las.return_number[ground_pts], return_counts=True)
print('Ground Point Return Number distribution:')
for r,c in zip(bins,counts):
print('{}:{}'.format(r, c))
Como resultado obtenemos:
Trabajar con numpy arrays para obtener informacion de la nube de puntos
En este último paso vamos a crear un «numpy array» con los valores en Z de nuestra nube de puntos y vamos a calcular los valores máximos, mínimos y medios de dicha nube de puntos. Para ello, utiliza el siguiente código:
import laspy as lp
import numpy as np
input_path = "C:\Javi\Movidas\laspy\PNOA_1.las"
point_cloud = lp.read(input_path)
pointsZ = np.vstack((point_cloud.z)).transpose()
print('\nAltitud maxima:', np.max(pointsZ))
print('\nAlitud media:', np.mean(pointsZ))
print('\nAltitud minima:', np.min(pointsZ))
Paso por paso, importamos las librerías de laspy y numpy, después leemos con laspy nuestra nube de puntos. A continuación usamos «vstak» para convertir los datos en Z de nuestra nube de puntos en un «numpy array».
El último paso consiste en usar las funciones max, mean y min para leer los valores máximos, mínimos y medios en Z. También podrías obtener información sobre quartiles, desviación estándar, mediana etc.
Ahora que sabemos los valores mínimos y máximo de alturas, ¿no sería interesante conocer la distribución de los mismos?
Podemos conseguir esto creando un histograma con 5 grupos iguales (bins), así averiguaremos como se distribuyen los datos desde 4 metros hasta 963 metros de altura. Usaremos nuevamente numpy para conseguir nuestro objetivo:
input_path = "C:\Javi\Movidas\laspy\PNOA_1.las"
point_cloud = lp.read(input_path)
pointsZ = np.vstack((point_cloud.z)).transpose()
print('\nAltitud maxima:', np.max(pointsZ))
print('\nAlitud media:', np.mean(pointsZ))
print('\nAltitud minima:', np.min(pointsZ))
histo = (np.histogram(pointsZ, bins=[4, 200, 400, 600, 800, 963]))
La función histograma de numpy, crea la distribución de puntos acorde con el tamaño de bin que nosotros deseemos.
Si nos fijamos en los resultados, vemos que la inmensa mayoría de nuestros valores se encuentran entre 4 y 200 metros. Lo cual ya nos indica que valores superiores a 200 metros serán probablemente «outliers» o valores atípicos.
ConclusiÓn
Como hemos podido comprobar, las librerías de laspy en combinación con numpy, nos dan mucha versatilidad para obtener información importante sobre nubes de puntos lidar. En posteriores posts veremos como crear nuevas nubes de puntos, filtrar por RGB, eliminar valores atípicos, etc.
El código que hemos usado quedaría del siguiente modo:
import laspy as lp
import numpy as np
## open file
input_path = "C:\Javi\Movidas\laspy\PNOA_1.las"
point_cloud = lp.read(input_path) # read file
pointsZ = np.vstack((point_cloud.z)).transpose() # convert Z values into np array
print('\nAltitud maxima:', np.max(pointsZ)) # print some metrics
print('\nAlitud media:', np.mean(pointsZ))
print('\nAltitud minima:', np.min(pointsZ))
histo = (np.histogram(pointsZ, bins=[4, 200, 400, 600, 800, 963])) ## find out Z distribution
print('\nHistograma: ', histo)
with lp.open(input_path) as fh:
las = fh.read()
print('Points from data:', len(las.points)) ## read number of points from file
with lp.open(input_path) as fh:
las = fh.read()
ground_pts = las.classification == 2
bins, counts = np.unique(las.return_number[ground_pts], return_counts=True)
for r,c in zip(bins,counts):
print('{}:{}'.format(r, c)) ## find out ground point distribution
gistonic 