Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Taking Slices from LiDAR data: Part IV

Posted by smathermather on February 5, 2016

Trans-Alaska Pipeline System Luca Galuzzi 2005

In PDAL, a pipeline file can be used to do a variety of operations. Within the following context, I think of a pipeline file like an ad hoc preferences file, where I can use an external command to iterate through the things I want to change, while holding constant everything else in the pipeline file.

In my use case for this vignette, I’ll use the pipeline file to hold my database preferences for putting the point clouds into my PostGIS database. For the record, I’m using vpicavet’s docker pggis as the starting place for installing PostGIS with the pgpointcloud extension. I have adapted the following pipeline file from the PDAL writers.pgpointcloud example.

<?xml version="1.0" encoding="utf-8"?>
<Pipeline version="1.0">
  <Writer type="writers.pgpointcloud">
    <Option name="connection">
      host='localhost' dbname=‘lidar’ user='user' password=‘password’
    <Option name="table">54001640PAN_heightasz_0-1.5</Option>
    <Option name="compression">dimensional</Option>
<!--    <Option name="srid">3270</Option> -->
    <Filter type="filters.chipper">
      <Option name="capacity">400</Option>
      <Reader type="readers.las">
          <Option name="filename">54001640PAN_heightasz_0-1.5.las</Option>
<!--          <Option name="spatialreference">EPSG:3270</Option> -->

Some things to note. I have commented out the SRID and readers.las.spatialreference in the XML above. We’ll rely on PDAL to discover this, and use the default output of epsg:4326 to store our point clouds for the time being.

Our wrapper script for the pipeline file is very simple. We will use the wrapper script to specify the table and the input file name.


TABLE=`basename $1 .bpf`

#echo $TABLE
pdal pipeline -i pipeline.xml --writers.pgpointcloud.table=$TABLE --readers.bpf.filename=$INPUT

Now to use, we’ll use GNU Parallel, as much for it’s XArgs like functionality as scalability:

ls *.bpf | parallel -j6 ./ {}

Now we can see what tables got loaded:

psql -d lidar
psql (9.5.0)
Type "help" for help.

lidar=# \dt
                    List of relations
 Schema |             Name              | Type  |  Owner  
 public | 54001640PAN_heightasz_0-1.5   | table | user
 public | 54001640PAN_heightasz_1.5-6   | table | user
 public | 54001640PAN_heightasz_100-200 | table | user
 public | 54001640PAN_heightasz_30-45   | table | user
 public | 54001640PAN_heightasz_45-55   | table | user
 public | 54001640PAN_heightasz_55-65   | table | user
 public | 54001640PAN_heightasz_6-30    | table | user
 public | 54001640PAN_heightasz_65-75   | table | user
 public | 54001640PAN_heightasz_75-85   | table | user
 public | 54001640PAN_heightasz_85-100  | table | user
 public | pointcloud_formats            | table | user
 public | spatial_ref_sys               | table | user
(12 rows)

W00t! We’ve got point clouds in our database! Next, we will visualize the data, and extract some analyses from it.

Posted in 3D, LiDAR, Other, PDAL, pointcloud | Tagged: , , | Leave a Comment »

Taking Slices from LiDAR data: Part III

Posted by smathermather on February 5, 2016

Borrowed from:

Continuing my series on slicing LiDAR data in order to analyze a forest, one of the objectives of the current project is to understand the habitats that particular species of birds prefer. This will be accomplished using field info from breeding bird surveys combined with LiDAR data of forest structure to help predict what habitats are necessary for particular species of breeding birds.

There are a number of studies doing just this for a variety of regions which I will detail later, but suffice it to say, structure of vegetation matters a lot to birds, so using LiDAR to map out structure can be an important predictive tool for mapping bird habitat. Being able to do that at scale across entire ecoregions—well, that’s just an exciting prospect.
Let’s get to some coding. I would like to take a few slices through our forest based on functional height groups:

Forest Canopy >15 meters
Sub-Canopy 10-15 meters
Tall Shrub 2-10 meters
Short Shrub 0.5-2 meters

(For the record, we don’t have to get these perfectly right. Sub-canopy, for example could be taller or shorter depending on how tall a forest it is in. This is just a convenience for dividing and reviewing the data for the time being. Also, we’ll split a little finer than the above numbers just for giggles.)

We’ll start by just echoing our pdal translate to make sure we like what we are getting for output:

for START in 0:0.5 0.5:1 1:2 2:5 5:10 10:15 15:20 20:35 35:50 50:75 75:100 100:200
  nameend=`echo $START | sed s/:/-/g`
  echo pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_$nameend.bpf -f range --filters.range.limits="Height[$START)"

Thus we get the following output:

pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_0-0.5.bpf -f range --filters.range.limits=Height[0:0.5)
pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_0.5-1.bpf -f range --filters.range.limits=Height[0.5:1)
pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_1-2.bpf -f range --filters.range.limits=Height[1:2)
pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_2-5.bpf -f range --filters.range.limits=Height[2:5)
... etc.

Let’s remove our echo statement so this actually runs:

for START in 0:0.5 0.5:1 1:2 2:5 5:10 10:15 15:20 20:35 35:50 50:75 75:100 100:200
  nameend=`echo $START | sed s/:/-/g`
  pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_$nameend.bpf -f range --filters.range.limits="Height[$START)"

We should generalize this a bit too. Let’s make this a script to which we can pass our filenames and ranges:

namer=`basename $1 .bpf`
for START in $2
  nameend=`echo $START | sed s/:/-/g`
  pdal translate $namer.bpf $namer"_"$nameend".bpf" -f range --filters.range.limits="Height[$START)"

Which to run, we use a statement as follows:

./ 54001640PAN_heightasz.bpf "0:0.5 0.5:1 1:2 2:5 5:10 10:15 15:20 20:35 35:50 50:75 75:100 100:200"

I forgot all my data is in feet, but my height classes in meters, so we’ll just multiply our input values by a factor of 3 to get in the same ballpark (future refinement likely):

./ 54001640PAN_heightasz.bpf "0:1.5 1.5:3 3:6 6:15 15:30 30:45 45:60 60:105 105:150 150:200"

We could alternatively stick to our original categories (short shrub, tall shrub, sub-canopy, canopy) as our break points:

./ 54001640PAN_heightasz.bpf "0:1.5 1.5:6 6:30 30:45 45:200"

Finally, we can convert to laz, and load all our slices of point clouds in an animate between the slices

for OUTPUT in $(ls *.bpf); do docker run -v /home/gisuser/test/test:/data pdal/master pdal translate //data/$OUTPUT //data/$OUTPUT.laz; done

If you look closely, you should be able to see where a tornado in the 80s knocked down much of the forest here. It’s signature is in the tall shrub / sub-canopy layer:

Posted in 3D, LiDAR, Other, PDAL, pointcloud | Tagged: , , | 2 Comments »

Point Clouds in Sketchup?

Posted by smathermather on February 4, 2016

Wow. Thanks to my colleague Brandon Garman for this tip — Point Clouds in Sketchup. Mind you, I don’t think this is an inexpensive proposition, but a pretty cool looking tool:

Posted in 3D, LiDAR, Other | Tagged: , | Leave a Comment »

Taking Slices from LiDAR data: Part II

Posted by smathermather on February 3, 2016

Ok, with a little help from Bradley Chambers on the PDAL mailing list, we are back in business. If we want to filter our newly calculated heights into a new PDAL output, we can do that easily, say all points 100-500 above ground level:

pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_gt100.bpf -f range --filters.range.limits=&quot;Height[100:500]&quot;

A little sanity check to see if we are getting appropriate values:

pdal info --stats 54001640PAN_heightasz_gt100.bpf --dimensions &quot;Height&quot;
  &quot;filename&quot;: &quot;54001640PAN_heightasz_gt100.bpf&quot;,
  &quot;pdal_version&quot;: &quot;1.1.0 (git-version: 64c722)&quot;,
        &quot;average&quot;: 105.8738232,
        &quot;count&quot;: 179909,
        &quot;maximum&quot;: 194.5800018,
        &quot;minimum&quot;: 100,
        &quot;name&quot;: &quot;Height&quot;,
        &quot;position&quot;: 0

Ok, now I want to view this. I could convert to a *.laz file and view it with (as long as I use Chrome as my browser).
I have to switch back to docker, ’cause that’s where I have PDAL built with laszip:

docker run -v /home/gisuser/test:/data pdal/master pdal translate //data/54001640PAN_heightasz_gt100.bpf //data/54001640PAN_heightasz_gt100.laz

And now I can view in

Screen Shot 2016-02-03 at 9.57.23 PM

Posted in 3D, Other, PDAL | Tagged: , , | Leave a Comment »

Taking slices from LiDAR data

Posted by smathermather on January 28, 2016


Welp, my BASH skills could use some honing, but I’m just working on quick-and-dirty stuff for this series. PDAL as a utility is pretty interesting, so we’ll focus our learnings on PDAL. (Prior posts start here). I may have a build tutorial written for getting through the critical (as contrasted with minimal) parts of a PDAL build. I’m getting tired of running everything through docker. Besides, I think running through docker is causing me stability problems when I try to run on 40 cores or so… .

If you recall, PDAL is Point Data Abstraction Library. For those who know GDAL, this is the equivalent of GDAL. Much like GDAL simplifies and unifies access to raster data, PDAL is meant to allow people to write software that uses point clouds without having to bother too much with understanding the underlying structure of the data.

Once we have heights calculated for everything, now it’s time to slice up the data according to height above ground. I really wanted to do this in PDAL, but haven’t yet figured out how. I tried using a pcl filter. The syntax is simple enough:

pdal pcl -i /path/to/input/las -p /path/to/pcl/block/json -o /path/to/output/las

And so I create a filter that will filter based on elevation:

        "name": "PassThroughExample",
                "name": "PassThrough",
                "setFilterFieldName": "z",
                    "min": 2000,
                    "max": 2400

And then I use my filter with the pdal pcl command:


name=`basename $1 .bpf`

docker run -v $pathname:/data pdal/master pdal pcl -i //data/"$name".bpf -p //data/elev.filter.json -o //data/"$name"_5.0.bpf;

Sadly however, if I try to adapt for height, I get no points out:

        "name": "PassThroughExample",
                "name": "PassThrough",
                "setFilterFieldName": "Height",
                    "min": 0.5,
                    "max": 5

More investigating to do… .

Posted in 3D, LiDAR, Other, PDAL | Tagged: , | Comments Off on Taking slices from LiDAR data

parallel processing in PDAL

Posted by smathermather on January 28, 2016

Frankfurt Airport tunnel.JPG
Frankfurt Airport tunnel” by Peter IsotaloOwn work. Licensed under CC BY-SA 3.0 via Commons.

In my ongoing quest to process all the LiDAR data for Pennsylvania and Ohio into one gigantic usable dataset, I finally had to break down and learn how to do parallel processing in BASH. Yes, I still need to jump on the Python band wagon (the wagon is even long in the tooth, if we choose to mix metaphors), but BASH makes me soooo happy.

So, in a previous post, I wanted to process relative height in a point cloud. By relative height, I mean height relative to ground. PDAL has a nice utility for this, and it’s pretty easy to use, if you get PDAL installed successfully.

pdal translate 55001640PAN.las 55001640PAN_height.bpf height --writers.bpf.output_dims=&amp;quot;X,Y,Z,Height&amp;quot;;

Installing PDAL is not too easy, so I used the dockerized version of PDAL and it worked great. Problem is, the dockerized version complicates my commands on the command line, especially if I want to run it on a bunch of files.

Naturally, the next step is to run it on a whole bunch of LiDAR files. For that I need a little control script which I called, and then I need to run that in a “for” loop.

# Get the pathname from the input value

# Get the short name of the file, sans path and sans extension
name=`basename $1 .las`

docker run -v $pathname:/data pdal/master pdal translate //data/"$name".las //data/"$name"_height.bpf height --writers.bpf.output_dims="X,Y,Z,Intensity,ReturnNumber,NumberOfReturns,ScanDirectionFlag,EdgeOfFlightLine,Classification,ScanAngleRank,UserData,PointSourceId";

Now we need a basic for loop will take care of sending the las files into our, thus looping through all available las files:

for OUTPUT in $(ls *.las); do ~/./ $OUTPUT; done;

This is great, but I calculated it would take 13 days to complete on my 58366 LiDAR files. We’re talking approximately 41,000 square miles of non-water areas for Ohio, and approximately 44,000 square miles of non-water areas for Pennsylvania. I’m on no particular timeline, but I’m not really that patient. Quick duckduckgo search later, and I remember the GNU Parallel project. It’s wicked easy to use for this use case.

ls *.las | parallel -j24 ~/./

How simple! First, we list our las files, then we “pipe” them as a list to parallel, we tell parallel we want it to spawn 24 independent processes using that list as the input for our pdal_height script.

Now we can run it on 24 cores simultaneously. Sadly, I have slow disks :( so really I ran it like this:

ls *.las | parallel -j6 ~/./

Time to upgrade my disks! Finally, I want to process my entire LiDAR dataset irrespective of location. For this, we use the find command, name the starting directory location, and tell it we want to search based on name.

find /home/myspecialuser/LiDAR/ -name "*.las" | parallel -j6 ~/./

Estimated completion time: 2 days. I can live with that until I get better disks. Only one problem, I should make sure this doesn’t stop if my network drops for any reason. Let’s wrap this in nohup which will prevent network-based hangups:

nohup sh -c 'find /home/myspecialuser/LiDAR -name "*.las" | parallel -j6 ~/./ {}'

Posted in 3D, Analysis, parallel, PDAL, pointcloud | Tagged: , , , | Leave a Comment »

wget for downloading boatloads of data

Posted by smathermather on January 27, 2016

My current project to create a complete dataset of airborne LiDAR data for Ohio and Pennsylvania has been teaching me some simple, albeit really powerful tricks. We’ll just discuss one today — recursive use of wget. This allows us to download entire trees of web sites to mirror, or in our case download all the data. Additionally, wget works on ftp trees as well, with no real change in syntax.

wget -r


I recognize that curl can be a more powerful utility than wget, but for dead simple problems like this, a dead simple tool is more than adequate.

Posted in 3D, Docker, PDAL, pointcloud | Tagged: , , , | Leave a Comment »

PDAL and point cloud height

Posted by smathermather on January 20, 2016


PDAL now has the capacity to calculate heights from your point cloud data. With pre-classified LiDAR data, this means you can do this pretty easily:

pdal translate 55001640PAN.las 55001640PAN_height.bpf height --writers.bpf.output_dims="X,Y,Z,Height;"

A problem you might have is you may not have all the wonderful PDAL goodness built and installed. So you might get something like this:

PDAL: Couldn't create filter stage of type 'filters.height'

An easy way around this is to let docker do all the work. Once the docker version of pdal is in place, we need to change our command a little to suit our new configuration:

docker run -v /home/gisuser/test:/data pdal/master pdal translate //data/54001640PAN.las //data/54001640PAN_height.bpf height --writers.bpf.output_dims="X,Y,Z,Height"

Posted in 3D, Docker, Other, PDAL | Tagged: , , , | 1 Comment »

Point Clouds – the (re)start of a journey

Posted by smathermather on January 19, 2016

If you have followed this blog for a while, you will notice a continual returning to and refinement of ideas and topics. That’s how the blog started, and this role it has served, as a touch stone in my exploration of topics is critical to how I use it. I hope it is useful to you too, as a reader.

So, let’s talk about point clouds again. Point clouds are like ordinary geographic points — they have X, Y, and Z (typically), but there are a whole lot more of them than ordinary. How much more? Well, instead of thinking in the 100s of thousands, we think more along the lines of 100 of millions or billions of points.

The other thing we need to recognize about point clouds is that they may be n-dimensional — a LiDAR point cloud may have X, Y, Z, intensity, return number, class, etc..

Good Free and Open Source Software (FOSS) tools for dealing with point clouds include a spectrum of projects, including PDAL (Point Data Abstraction Library), and the Pointcloud extension for PostgreSQL (with PostGIS integration, of course).

Previously, I was attempting to process large amounts of LiDAR point clouds in order to model bird habitat. Now that I have the requisite storage, and PDAL has a height calculator, I am ready to dive back into this.

In the meantime, there are some cool things to check in this space. Let’s start with this presentation by Michael Smith on PDAL from FOSS4G in Seoul last year.

Screen Shot 2016-01-19 at 8.32.03 PM

Want to see some of the incredible things you can do in the browser, go no further than here:

Screen Shot 2016-01-19 at 8.34.52 PM.png

This goes deeper than simple vector tile approaches and employs level of detail optimizations that are necessary for this large of a use case. More on this later

postscript 1/20/2016:

Looks like I outed the above project… . Apologies to for this. For the record and from their site:

Hobu is a team of three software engineers located in Iowa City, Iowa and Cambridge, Iowa. We have been on the forefront of open source LiDAR software for over five years, and we have been building open source GIS software for even longer. We provide open source lidar software systems development, open source GIS software development in addition to contract research and development, systems design, evaluation and implementation.

Please contact us about any of the open source software for which we provide support.

We are the primary support organization for the following libraries:

  • PDAL — Point Data Abstraction Library
  • — WebGL point cloud rendering
  • Greyhound — Point cloud data streaming
  • libLAS — ASPRS LAS read/write

Welp, there we go, an unintentional public launch. Awesome work though… .

Posted in 3D, Database, FLANN, pointcloud, PostGIS, PostgreSQL, SQL | Tagged: , , | 1 Comment »

Moar generative PostGIS Art

Posted by smathermather on January 15, 2016

Screen Shot 2016-01-15 at 12.07.06 AM
Screen Shot 2016-01-15 at 12.09.33 AMScreen Shot 2016-01-15 at 12.09.18 AM
Screen Shot 2016-01-15 at 12.18.08 AM
Screen Shot 2016-01-15 at 12.11.17 AM
DROP TABLE IF EXISTS cuyahogaaaahh;
CREATE TABLE cuyahogaaaahh AS
SELECT row_number() OVER() As id, geom
SELECT ST_Buffer(ST_SubDivide(ST_Buffer(
     geom, 50, 100), generate_series), generate_series) AS geom
	FROM generate_series(8,50), cuyahoga
) subdivide

Posted in Other, PostGIS, PostgreSQL | Tagged: | Leave a Comment »


Get every new post delivered to your Inbox.

Join 1,639 other followers