Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for January, 2016

Taking slices from LiDAR data

Posted by smathermather on January 28, 2016

maxresdefault

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:

{
    "pipeline":
    {
        "name": "PassThroughExample",
        "filters":
        [
            {
                "name": "PassThrough",
                "setFilterFieldName": "z",
                "setFilterLimits":
                {
                    "min": 2000,
                    "max": 2400
                }
            }
        ]
    }
}

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

#!/bin/bash

pathname="${1%/*}"
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:

{
    "pipeline":
    {
        "name": "PassThroughExample",
        "filters":
        [
            {
                "name": "PassThrough",
                "setFilterFieldName": "Height",
                "setFilterLimits":
                {
                    "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="X,Y,Z,Height";

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 pdal_height.sh, and then I need to run that in a “for” loop.

#!/bin/bash
# Get the pathname from the input value
pathname="${1%/*}";

# 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 pdal_height.sh, thus looping through all available las files:


for OUTPUT in $(ls *.las); do ~/./pdal_height.sh $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 ~/./pdal_height.sh

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 ~/./pdal_height.sh

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 ~/./pdal_height.sh

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 ~/./pdal_height.sh {}'

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 ftp://pamap.pasda.psu.edu/pamap_lidar/cycle1/LAS/

Capture

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

point_cloud

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 https://hobu.co 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
  • plas.io — 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
FROM
(
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 »

Generative art in PostGIS

Posted by smathermather on January 14, 2016

Screen Shot 2016-01-14 at 11.46.55 PMScreen Shot 2016-01-14 at 11.45.55 PMScreen Shot 2016-01-14 at 11.45.41 PM

SELECT ST_Buffer(ST_SubDivide(ST_Buffer(
     ST_GeomFromText(
      'POINT(0 0)'
     ), 50000, 100), generate_series), generate_series * 100) AS geom
	FROM generate_series(8,100)

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

Yet another approach to ST_Buffer on geography

Posted by smathermather on January 13, 2016

Another approach to ST_Buffer would be to subdivide the geometries before buffering, and put them all together at the end. ST_SubDivide can do this for us. We can tell it how may vertices we want in each geometry (minimum of 8). Since _ST_BestSRID will try UTM first, we’ll add enough nodes to ensure we always have 8 nodes within the 1,000,000 meter width of a UTM zone by segmentizing at 62,500 meters.

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$

WITH subdivided AS (
 SELECT
 CASE WHEN ST_GeometryType(g1::geometry) = 'ST_Point' THEN g1
 ELSE (SELECT ST_SubDivide(ST_Segmentize(g1::geometry, 62500), 8) AS geom)
 END
),
transformed_local AS (
 SELECT ST_Transform(g1::geometry, _ST_BestSRID(geom)) AS geom FROM subdivided
),
buffered AS (
 SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
 FROM transformed_local
),
transformed_4326 AS (
 SELECT (ST_Dump(ST_Transform(geom, 4326))).geom AS geom FROM buffered
),
checksrid AS (
 SELECT geom FROM transformed_4326
)

SELECT ST_Union(geom)::geography FROM checksrid;
 
$BODY$
LANGUAGE sql VOLATILE;

Aaand, a somewhat unrelated image (just the buffered subdivision of a buffered point, not really an output from our function above. More on all this later.

Screen Shot 2016-01-13 at 6.38.15 PM

Posted in Database, PostGIS, PostgreSQL, SQL | Tagged: , , | 1 Comment »

ST_Buffer on Geography — Iteration

Posted by smathermather on January 12, 2016

In my previous post, Imagico said:

This is nice but still fails in one important aspect – it will create significant inaccuracies once your buffering distance reaches values where buffering in a single local coordinate system is getting problematic which, depending on your precision requirements, will likely happen at distances of a few degrees. A solution to this would require either an iterative split and merge technique of a completely different approach on a very basic level (i.e. attempting truly spherical buffering).

Ok, ok. I kinda forgot about this aspect. The proper way to do this is, as always is to do it at a lower level thank in an SQL function, preferably calculating the buffer on the spheroid, ellipsoid, or geoid. But, ouch. I think we can munge something good enough for now.

We’ll still do the buffer:

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
...

But now that we’ve performed the buffer, we need to determine that the buffered geometry is not too large for our local coordinate system. To do that, we’ll compare the best SRID for our buffered geometry to the original best SRID that we calculated based on the original input geometry.

SELECT
CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography

If the two SRIDs aren’t equal, then we need to use the new best SRID based on buffered geometry and buffer using that new best SRID.

ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, num_seg_quarter_circle), 4326)::geography

Now let’s put it all together:

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$
 
WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT
 CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography 
 ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, num_seg_quarter_circle), 4326)::geography
 END
 FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;


CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
buffer_style_parameters text)
RETURNS geography AS
$BODY$
 
WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, buffer_style_parameters) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT
 CASE WHEN _ST_BestSRID(geom) = _St_BestSRID(g1) THEN geom::geography 
 ELSE ST_Transform(ST_Buffer(ST_Transform(g1::geometry, _ST_BestSRID(geom)), radius_of_buffer, buffer_style_parameters), 4326)::geography
 END
 FROM transformed_4326
 
$BODY$
LANGUAGE sql VOLATILE;

So, what’s the affect? Pretty significant really.

SELECT ST_Buffer(
	 ST_GeomFromText(
	  'POINT(0 0)'
	 )::geography, 5000000, 50);
Blue line showing correct buffer, dashed red line showing distorted buffer from using wrong local coordinate system.

Blue line showing correct buffer, dashed red line showing distorted buffer from using wrong local coordinate system.

Posted in Database, PostGIS, PostgreSQL, SQL | Tagged: , , | 2 Comments »

Final?: ST_Buffer on Geography, or My love for @NullIsland

Posted by smathermather on January 11, 2016

Thanks to Paul Norman’s reminder in a previous post, we now have all the pieces we need to complete an ST_Buffer function that exposes all the wonderful goodness of buffer style parameters to geography users and also chooses the best local coordinate system automatically.

We use _ST_BestSRID to choose our local coordinate system.

Remember my previous three disadvantages:

  1. It didn’t do this at a low enough level to automatically use the best available local coordinate system. This leaves it to the user to choose the best available local coordinate system
  2. It used a different function name than you otherwise use for buffering.
  3. Finally, it wouldn’t let you use all the other parameters that ST_Buffer exposes through the use of buffer style parameters.

They are now all solved.


CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
num_seg_quarter_circle integer)
RETURNS geography AS
$BODY$

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, num_seg_quarter_circle) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT geom::geography FROM transformed_4326

$BODY$
LANGUAGE sql VOLATILE;

 

CREATE OR REPLACE FUNCTION ST_Buffer(g1 geography, radius_of_buffer float,
buffer_style_parameters text)
RETURNS geography AS
$BODY$

WITH transformed AS (
SELECT ST_Transform(g1::geometry, _ST_BestSRID(g1)) AS geom
),
buffered AS (
SELECT ST_Buffer(geom, radius_of_buffer, buffer_style_parameters) AS geom
FROM transformed
),
transformed_4326 AS (
SELECT ST_Transform(geom, 4326) AS geom FROM buffered
)
SELECT geom::geography FROM transformed_4326

$BODY$
LANGUAGE sql VOLATILE;

Let’s try this!

SELECT ST_Buffer(
ST_GeomFromText(
'LINESTRING(5 2.5,0 -2.50,-5 2.5)'
)::geography, 500000, 'join=mitre');
Heart shaped buffer over Null Island.

Heart shaped buffer over Null Island.

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

(poor editing earlier tonight. My apologies.

Posted in Database, PostGIS, PostgreSQL, SQL | Tagged: , , | 2 Comments »