These are chat archives for PDAL/PDAL

7th
Nov 2017
matt wilkie
@maphew
Nov 07 2017 20:50
Hello World! I'm investigating use of PDAL to extract Bare Earth points from a LAS source and then feed to ArcGIS TopoToRaster (or ANUDEM) in order to end up with hydrologically correct DEM. .... I've figured out how to get Bare Earth using pdal ground, following the tutorial, but when I then convert to shapefile there aren't any Z values. What am I missing?
pdal translate -w ogr -i thing_01_BE.laz -o thing_01_BE.shp
pdal ground --verbose 4 -i thing_01.laz -o thing_01_BE.laz (posted in wrong order)
Howard Butler
@hobu
Nov 07 2017 22:00
@maphew what do you want written "to shapefile"?
rasterize your bare earth with writers.gdal and then go do GDAL things to it :)
note that PDAL doesn't do hydro enforcement of anything, so you'll need to clip out your waterbodies and flatten them using GDAL/whatever
matt wilkie
@maphew
Nov 07 2017 22:02
I want the rasterization part handled by TopoToRaster/Anudem. It's the hydro enforcement part that it's good at. I don't want to throw made up values (interpolated values) at it.
Howard Butler
@hobu
Nov 07 2017 22:03
so what do you need PDAL to be writing for output? LAS?
matt wilkie
@maphew
Nov 07 2017 22:03
basically I'm trying get a point cloud of X,Y,Z
where Z is ground
Howard Butler
@hobu
Nov 07 2017 22:05
So if you write a shapefile right now with PDAL (is this OSGeo4W?), you don't get Z's in the multipoint shapefile?
(also, good to see you again!)
matt wilkie
@maphew
Nov 07 2017 22:05
if i use TopoToRaster the input has to be shapefile points with an elevation attribute. If Anudem it's just ascii (but feeding the parameters is a pain. It's a fortran program)
(nice to chat with you too!)
Howard Butler
@hobu
Nov 07 2017 22:06
so every point as its own shapefile entry, with an attribute for Z?
matt wilkie
@maphew
Nov 07 2017 22:06
affirmative
Howard Butler
@hobu
Nov 07 2017 22:06
not ESRI-style multipointZ
I must say that's a very stupid point cloud format :P
matt wilkie
@maphew
Nov 07 2017 22:06
right no multipoint (which is stupid since that's Esri's own format but that's the way it is)
it all boils down to fortran
Howard Butler
@hobu
Nov 07 2017 22:07
so ASCII output is ok too? or do you have a workflow for that?
matt wilkie
@maphew
Nov 07 2017 22:08
I can work with ascii. I was just trying a slightly more efficient format first ;-)
Howard Butler
@hobu
Nov 07 2017 22:08
I doubt one-point-per-shape-record + attribute is that much cheaper :)
You might be able to get what you want with some OGR VRT magic, but it isnt' going to be worth a half of day effort
especially if you have ASCII available as an option.
Even for all of NWT
Andrew Bell
@abellgithub
Nov 07 2017 22:11
Why would LAS not work but text would? Just curious...
Howard Butler
@hobu
Nov 07 2017 22:12
The next tools in the chain only support ASCII or the silly shapefile contraption
matt wilkie
@maphew
Nov 07 2017 22:12
It boils down to workflow. When the input is XYZ.txt I have to create an Event layer before feeding T2R. Which is dumb, because then T2R converts it back to txt before feeding the backend engine. But. I hate writing the straight-front-load-to-engine parameter files by hand, because it takes so many tries and retries to get right.
Howard Butler
@hobu
Nov 07 2017 22:13
I think you could get the shapefiles you want pretty quickly with a little python https://pypi.python.org/pypi/PDAL
PDAL Python extension doesn't exist on OSGeo4W though at the moment. I could package it up for you though. It works there
ooh, I could use the py wrapper I think
(direct from pypi i sfine)
Howard Butler
@hobu
Nov 07 2017 22:14
are you linux, docker, or windows?
matt wilkie
@maphew
Nov 07 2017 22:14
win
Howard Butler
@hobu
Nov 07 2017 22:14
ok, I'll go make a wheel
matt wilkie
@maphew
Nov 07 2017 22:15
Oh wow, thanks!!
Howard Butler
@hobu
Nov 07 2017 22:16
would be python3 only though.
You should be able to OGR+python3 in OSGeo4W64
matt wilkie
@maphew
Nov 07 2017 22:17
You know, if you ever get bored (hah!) implementing Hutchinson's drainage enforcement algorithms in a modern language would be huge
I can do py3 using miniconda; haven't tried to load osgeo stuff into that yet though. ... I saw py3 get selected as a requirement when I installed osgeo pdal. That should work
Howard Butler
@hobu
Nov 07 2017 22:19
yep. PDAL/OSGeo4W is python3 only for now. It embeds Python3 for being able to do inline python processing on data
matt wilkie
@maphew
Nov 07 2017 22:20
What's practical diff between SMRF ground filtering and PMF? I'm not sure how to choose
Howard Butler
@hobu
Nov 07 2017 22:21
SMRF is probably quite a bit faster and is going to work better in flatter terrain with default settings
PMF might do better in urban situations
if you just want to strip the trees off, use SMRF. Make sure to denoise first though
matt wilkie
@maphew
Nov 07 2017 22:22
hmm, I'm working with mixed urban/natural. I'll start with the fast one!
denoise?
Howard Butler
@hobu
Nov 07 2017 22:23
Does your lidar have a lot of atmospheric noise in it? Any of the ground algos are going to be quite sensitive to it if the noise has not already been classified
https://www.pdal.io/workshop/exercises/analysis/denoising/denoising.html might help a bit, but there are many options in PDAL, and it's a quite deep topic
If your vendor provided data with classifications for noise already, that's probably a more convenient place to start
matt wilkie
@maphew
Nov 07 2017 22:24
I think so. Noise is in the default surface created by Esri EzLAS tools is what I'm trying to fix.
image.png
Howard Butler
@hobu
Nov 07 2017 22:25
all those little bumps are :christmas_tree: :)
matt wilkie
@maphew
Nov 07 2017 22:25
big pimple in centre has gotta go, plus removing all the trees and buildings
yup!
Howard Butler
@hobu
Nov 07 2017 22:26
big pimple could be removed with filters.statisticaloutlier or even a simple filters.range that chops based on Dimension=Z
matt wilkie
@maphew
Nov 07 2017 22:27
I'm a LAS newbie. There are classifications in the points, but I'm not sure what to do with them (or have any idea how good/reliable they are)
Howard Butler
@hobu
Nov 07 2017 22:30
I would start with their classifications before attempting to make your own with PDAL
matt wilkie
@maphew
Nov 07 2017 22:30
is public data (just not generally available)
Howard Butler
@hobu
Nov 07 2017 22:32
The classifications in all of that data appear to be 0.
pdal info Carmacks_01.laz
matt wilkie
@maphew
Nov 07 2017 22:34
ooooh.
Howard Butler
@hobu
Nov 07 2017 22:34
if the Classification section there had any range to it, you could count on there being some values that might be helpful.
matt wilkie
@maphew
Nov 07 2017 22:39
I see. I grabbed from the folder labelled "raw". When I run pdal info on classified\bare_earth\carmacks_001.laz I get a Classified min/max of 1-7. ... Silly me assumed raw meant "has everything"
Howard Butler
@hobu
Nov 07 2017 22:39
raw == "unclassified"
so, 7 == noise
you want filters.range with Classification[2:2] as the limits
that will likely get rid of your pimple
matt wilkie
@maphew
Nov 07 2017 22:40
...anyway, thank you so much for your time and advice. I've been spinning in circles for 2 days.
Howard Butler
@hobu
Nov 07 2017 22:40
how'd you get to PDAL?
matt wilkie
@maphew
Nov 07 2017 22:41
a comment to an answer somewhere on gis.stackexchange.com. I'm surprised I haven't heard more about it.
(I did do a write up on plas.io quite a while ago, but hadn't noticed the tools behind it)
Howard Butler
@hobu
Nov 07 2017 22:42
yeah I try to answer there a little bit. Feels like 80% marketing an GIS.se
Howard Butler
@hobu
Nov 07 2017 22:43
We've gotten better at that :)
matt wilkie
@maphew
Nov 07 2017 22:43
a whopping readership of 3 (counting @hobu)
Howard Butler
@hobu
Nov 07 2017 22:43
Go to http://potree.entwine.io/ and look at all of Netherland's lidar in a browser (30,000+ laz files)
(or minnesota, or Iowa, or ...)
matt wilkie
@maphew
Nov 07 2017 22:44
Damn
ww
wow
Howard Butler
@hobu
Nov 07 2017 22:45
shit, I forgot to DLL export some crap and need to refresh the whole OSGeo4W64 PDAL stack to account for it (this is to get the Python extension working on windows).
I'll set that all in motion (it's almost all automatic), but it takes a few hours
matt wilkie
@maphew
Nov 07 2017 22:46
no rush on my part. I have a lot of work to do in prep anyway
is there a page that describes how to translate syntax between json usage examples and command line parameters?
Howard Butler
@hobu
Nov 07 2017 22:49
do you mean for setting options and such?
matt wilkie
@maphew
Nov 07 2017 22:50
yeah
Howard Butler
@hobu
Nov 07 2017 22:50
you can set any option via command line by setting its stage name + option
--filters.range.limits=Classification[2:2]
(er, maybe need quotes around the Classification bit -- can't remember)
so what we typically do is write a pipeline json that sketches out the basic things we want to do -- filter noise, run smrf, write las -- and then pump that through an xargs or something and substitute
--readers.las.filename={} and --writers.las.filename=output/{}
Mike covered this in more detail at FOSS4G this year http://2017.foss4g.org/post_conference/PDAL.pdf
matt wilkie
@maphew
Nov 07 2017 22:53
I got an error when attempting to use multiple writer as per https://www.pdal.io/apps/translate.html, but when using gdal
Howard Butler
@hobu
Nov 07 2017 22:54
Multiple writers is iff'y
matt wilkie
@maphew
Nov 07 2017 22:54
pdal translate -w gdal --writer gdal.resolution=1 -i ... -o ...
PDAL: kernels.translate: Attempted to set value twice for argument 'writer'
Howard Butler
@hobu
Nov 07 2017 22:55
pdal translate in.las out.tif --writers.gdal.resolution=1
it'll use writers.gdal based on the tif extension
Mike's slide deck has some good material for you too
matt wilkie
@maphew
Nov 07 2017 22:56
ahhh, plural plus dot for the options. ... very nice about using extensions. That's a great feature. ... Thanks for the slide deck
Howard Butler
@hobu
Nov 07 2017 22:56
yeah, I wish GDAL had something like it
PDAL's Pipeline == GDAL VRT version 2
switch from pdal translate to pdal pipeline once you start making really messy command line invocations
having the pipeline around as a record of what you did is quite valuable
you can even burn that pipeline into the LAS file if you use the writers.las.pdal_metadata=true option
matt wilkie
@maphew
Nov 07 2017 22:59
(!) make that the default
Howard Butler
@hobu
Nov 07 2017 22:59
Lots of good content in here too https://www.pdal.io/workshop/exercises/index.html
don't want to default to that because it might adulterate files unexpectedly. We try to make users be explicit about what they want
matt wilkie
@maphew
Nov 07 2017 23:01
I'm tickled pink. Thanks for taking the time to educate. I've been second guessing myself all day about "wasting" time investigating yet another 3rd party toolkit when we already have all this other stuff (that cost so much to invest in). No more 2nding needed!
Howard Butler
@hobu
Nov 07 2017 23:01
There's definitely some gaps here in PDAL land
but we're first-rate workflow, translation, and processing for point cloud data
After I get your wheel, you can write some Fiona or ogr.py that reads the LAS data after it has been classified as ground and then writes it into your craptacular shapefile format
The same pipeline you write to define your range filtering can be used in python land too. The stuff is all designed to link and be used together
Another option, if you are pure LAS+Python (not LAZ) is laspy.
speed-wise it is likely to be similar to PDAL for this particular task, but you're not going to have access to all of PDAL's goodies, of course.
matt wilkie
@maphew
Nov 07 2017 23:07
+1
Michael's slide deck: very good intro for where I am right now.
Howard Butler
@hobu
Nov 07 2017 23:21
Have you ever worked with ANUDEM, @adamsteer ?
matt wilkie
@maphew
Nov 07 2017 23:22
I'm willing to answer any questions you have on it -- from a user's perspective. I've spent a lot of time on/in it over the years (though not for awhile now).
Howard Butler
@hobu
Nov 07 2017 23:25
such AML!
matt wilkie
@maphew
Nov 07 2017 23:26
&you &know &it!
matt wilkie
@maphew
Nov 07 2017 23:50
here's the better tasting snake-based coffee: https://github.com/envygeo/elev/tree/master/anudem (a mess)