Getting source bounds when filtering LAS using PDAL?

Question

I'm using PDAL to filter out non ground points to create a DTM along with a DSM also made through PDAL to create a Canopy Height Model using this method.

My workflow looks like this.

#!/bin/sh
# name=output name, i=filetype (LAS or LAZ)

name="$1"
i="$2"

mkdir DTM
mkdir DSM

for f in *.${i}
do
JSONFILE="{
  \"pipeline\":[
    \"${f}\",
     {
      \"type\":\"filters.range\",
      \"limits\":\"Classification[2:2], Classification[9:9]\"
    },
    {
      \"type\":\"writers.gdal\",
      \"resolution\": 0.5,
      \"radius\": 1.0,
      \"filename\":\"DTM/${f}.tif\",
      \"output_type\":\"min\"
    }
  ]
}"


echo "$JSONFILE" > ${f}_dtm.json

done


ls *.json | parallel pdal pipeline {}

rm *.json

for f in *.${i}
do
JSONFILE="{
  \"pipeline\":[
    \"${f}\",
    {
      \"resolution\": 0.5,
      \"radius\": 1.0,
      \"filename\":\"DSM/${f}.tif\",
      \"output_type\":\"max\"
    }
  ]
}"

echo "$JSONFILE" > ${f}_dsm.json

done

ls *.json | parallel pdal pipeline {}

rm *.json

ls DTM/*.tif | sed -e "s'DTM/''g" | parallel gdal_calc.py -A DSM/{} -B DTM/{} --outfile=DTM/{}_CHM.tif --calc="A-B" --NoDataValue=0

Note, I'm using parallel to be able speed up things using more cpu cores.

This workflow works well with most files, but sometimes when there aren't any ground points (Classification 2) all the way out to the file's extents, the output TIF gets cropped and will not have the same dimensions as the original LAS file (and DSM). When the DSM and DTM extents doesn't match I get an error in gdal_calc.py.

Now I wonder how I can solve this.

Either if PDAL could write a TIF with the same bounds as the original LAS file regardless of the extents of the filtered points. I've seen the bounds flag in writers.gdal but I have no idea how I extract these from the LAS file in a simple way.

Another option would perhaps be changing the extent of the TIFFs afterwards so it fits with the DSM file's extent. Not sure how to achieve this either. Perhaps by extracting extents from the DSM and then "applying" to the DTM tif.

Example images covering the same area. As you can see the bottom DTM is cropped because there isn't any non ground points in the southern part.

Answer 1

I finally came up with a solution thanks to a colleague who told me to use GREP.

For each LAS file (in a loop I used these lines of code):

pdal info ${f} > ${f}_lasinfo.json

MAXX=`grep "maxx" ${f}_lasinfo.json | sed -e 's/[^0-9.]//g'`
MAXY=`grep "maxy" ${f}_lasinfo.json | sed -e 's/[^0-9.]//g'`
MINY=`grep "miny" ${f}_lasinfo.json | sed -e 's/[^0-9.]//g'`
MINX=`grep "minx" ${f}_lasinfo.json | sed -e 's/[^0-9.]//g'`

pdal info creates a json-file including all kinds of information about the file. Some lines covers the bounds of the file which I then extracted with grep and sed.

After that I could add a new line to my pdal pipe:

"bounds":"([$MINX,$MAXX],[$MINY,$MAXY])"