First, we need to convert our LAS file to a format that can be used by PDAL. We created a Python script, which reads from a directory of LAS files and generates its corresponding JSON. With this script, we can automate the generation if we have a large directory of files. Also, we chose Python for its simplicity and because you can execute the script regardless of the operating system you are using. To execute the script, run the following in the console (for Windows users, make sure you have the Python interpreter included in the PATH variable):
$ python insert_files.py -f <lasfiles_path>
This script will read each LAS file, and will store within a folder called pipelines all the metadata related to the LAS file that will be inserted into the database.
Now, using PDAL, we execute a for loop to insert LAS files into Postgres:
$ for file in `ls pipelines/*.json`;
do
pdal pipeline $file;
done
This point cloud data is split into three different tables. If we want to merge them, we need to execute the following SQL command:
DROP TABLE IF EXISTS chp07.lidar;
CREATE TABLE chp07.lidar AS WITH patches AS
(
SELECT
pa
FROM "chp07"."N2210595"
UNION ALL
SELECT
pa
FROM "chp07"."N2215595"
UNION ALL
SELECT
pa
FROM "chp07"."N2220595"
)
SELECT
2 AS id,
PC_Union(pa) AS pa
FROM patches;
The postgres-pointcloud extension uses two main point cloud objects as variables: the PcPoint object, which is a point that can have many dimensions, but a minimum of X and Y values that are placed in a space; and the PcPatch object,which is a collection of multiple PcPoints that are close together. According to the documentation of the plugin, it becomes inefficient to store large amounts of points as individual records in a table.
Now that we have all of our data into our database within a single table, if we want to visualize our point cloud data, we need to create a spatial table to be understood by our layer viewer; for instance, QGIS. The point cloud plugin for Postgres has PostGIS integration, so we can transform our PcPatch and PcPoint objects into geometries and use PostGIS functions for analyzing the data:
CREATE TABLE chp07.lidar_patches AS WITH pts AS
(
SELECT
PC_Explode(pa) AS pt
FROM chp07.lidar
)
SELECT
pt::geometry AS the_geom
FROM pts;
ALTER TABLE chp07.lidar_patches ADD COLUMN gid serial;
ALTER TABLE chp07.lidar_patches ADD PRIMARY KEY (gid);
This SQL script performs an inner query, which initially returns a set of PcPoints from the PcPatch using the PC_Explode function. Then, for each point returned, we cast from PcPoint object to a PostGIS geometry object. Finally, we create the gid column and add it to the table as a primary key.
Now, we can view our data using our favorite desktop GIS, as shown in the following image:
