How the System Call That Ties Unix Together Came About – The New Stack

Pipe: How the System Call That Ties Unix Together Came About – The New Stack

It’s an everyday command in the life of developers, sysadmins, and Unix lovers everywhere. So it’s remarkable to remember that Unix’s pipe command was implemented in a single day, representing not only a great moment in computing history, but also a uniquely important moment for its profound impact on the culture of Unix.

And it changed the way we program ever since.

People were asking how I was extracting the digital terrain model (tree/building elevation) from the NYS GIS LiDAR point clouds (LAS) files

People were asking how I was extracting the digital terrain model (tree/building elevation) from the NYS GIS LiDAR point clouds (LAS) files. While the digital surface models (bare earth) models are widely available in raster digital elevation model GeoTIFF, these files are bare earth and don’t include vegetation or buildings. But you can easily create raster digital elevation model GeoTIFF from the LIDAR point clouds with the PDAL – Point Data Abstraction Library.

Install PDAL

It’s easy if you are running Ubuntu Linux — it’s in standard distribution with 20.04 LTS or later. Most modern distributions of Linux include it in their repositories.

sudo apt install pdal

Download LiDAR Point Clouds

You can download the individual cloud files from NYS GIS FTP site. You can use the LiDAR Shapefile indexes on that website to figure out which file you need. Be aware you may have to download and store a significant amount of data — each roughly 556 acre tile (less then a square mile) is between 500 and 1,000 MB. So you will want to do this somewhere you have a lot of free unlimited data service and storage on your hard drive. And patience while it downloads.

Extract the Digital Terrain Model

You should create a text file called pdal_dtm.json or something similar. The text file should contain:

{ 
	"pipeline": [
		"18TWN220520.las",
		{ 
			"type": "filters.range",
			"limits": "Classification[1:1]"
		},
		{
			"filename":"dsm.tif",
			"gdaldriver":"GTiff",
			"resolution": 0.5,
			"output_type": "max",
			"type":"writers.gdal"
		}
	]
}

This creates the pdal “pipeline” configuration file for the conversion. An explanation of important lines you will need to change:

“18TWN220520.las” – Input point cloud LAS file
“type”: “filters.range” – Tells PDAL to filter points based on request
“limits”: “Classification[1:1]” – Layer to extract *
“filename”:”dsm.tif” – Export file name
“gdaldriver”:”GTiff” – Use GeoTIFF export (same as standard DEM files)
“resolution”: 0.5 – Resolution in meters for export, typically 0.5m or 1m depending on survey
“output_type”: “max” – Highest point reflected back within that point, to get digital surface elevation. You can also use min/max/mean/count/stdev/idw.
“type”:”writers.gdal” – Use the GDAL library for create GeoTIFF

* Layers in NYS GIS Point Clouds – “These point clouds will have at a minimum 2 classifications; Class 1 Unclassified, and Class 2 Bare Earth.” Layer 1 is the digital surface model, it contains building heights and tree/crop cover. Layer 2 “Bare Earth” is the same as what you get from the LiDAR Digital Elevation GeoTIFFs that are widely available on NYSGIS website. Be aware that the Class 2 points may have no-data areas, where no ground elevation was detected due to tree cover or other things blocking ground. The DEM files you download from NYSGIS have these areas filled in with interpolation from surrounding areas.

Once your configuration file is done, run the command:

pdal pipeline pdal_dtm.json

As the Point Cloud files are internally indexed, the export should take only a matter of seconds on a modern desktop computer. The exported DEM/GeoTIFF file will be between 15 to 25 MB, which is much smaller then the point cloud. Then you can load in Quantum GIS or your favorite non-free software GIS client like any other Raster DEM file.

To calculate the building or tree height, just subtract the Class 1 raster from Class 2 raster, using the Raster Calculator in QGIS or your favorite GIS program.Β  Or better yet, subtract Class 1 points from the already processed DEM files on NYSGIS website, as they have the missing bare earth data filled in.