Lidar

Dependencies

%conda install gdown

Setup

# Import modules
import sys
import subprocess
from pathlib import Path
from contextlib import suppress
import gdown

# Find GRASS Python packages
sys.path.append(
  subprocess.check_output(
    ["grass", "--config", "python_path"],
    text=True
    ).strip()
  )

# Import GRASS packages
import grass.script as gs
import grass.jupyter as gj
from grass.tools import Tools

# Create and initialize GRASS project
try:
    gj.init("hilltop")
except ValueError:
    with suppress(ValueError):
        gs.create_project("hilltop", epsg="2801")
        gj.init("hilltop")
tools = Tools(overwrite=True)

# Print projection
projection = tools.g_proj(format="shell", flags="p").keyval
print(projection.get("name"))

NAD83(HARN) / Louisiana South

Download parcel

# Set url
url = "https://drive.google.com/file/d/1Okkjea42jx4omRXGFFujWadvOjQnDtQO/view?usp=drive_link"

# Set filename
output = "parcel.gpkg"

# Download data
gdown.download(url, output, quiet=False)

# Import parcel
tools.v_import(input="parcel.gpkg", output="parcel")

Downloading...
From: https://drive.google.com/uc?id=1Okkjea42jx4omRXGFFujWadvOjQnDtQO
To: /Users/baharmon/academia/computational-ecology/parcel.gpkg
100%|███████████████████████████████████████████████████████████████████████████████████| 98.3k/98.3k [00:00<00:00, 41.7MB/s]

Download lidar

# Set url
url = "https://drive.google.com/file/d/1FQPGd3MNjqB3itYwUBCZqqQjoN2mIz19/view?usp=drive_link"

# Set filename
output = "hilltop.laz"

# Download data
gdown.download(url, output, quiet=False)

Downloading...
From (original): https://drive.google.com/uc?id=1FQPGd3MNjqB3itYwUBCZqqQjoN2mIz19
From (redirected): https://drive.google.com/uc?id=1FQPGd3MNjqB3itYwUBCZqqQjoN2mIz19&confirm=t&uuid=09f3a12f-b061-4e95-8aca-e6d944399d84
To: /Users/baharmon/academia/computational-ecology/hilltop.laz
100%|█████████████████████████████████████████████████████████████████████████████████████| 203M/203M [01:03<00:00, 3.20MB/s]

'hilltop.laz'

Set region

# Set region
tools.g_region(vector="parcel", res=0.3)

Bin lidar

%%time

# Bin lidar
tools.r_in_pdal(
    input="hilltop.laz",
    output="surface",
    method="max",
    resolution=0.3,
    flags="oe"
    )

# Set color gradient
tools.r_colors(map="surface", color="viridis")

# Visualize
m = gj.Map(width=1600)
m.d_rast(map="surface")
m.d_legend(raster="surface", color="white", at=(75, 95, 2, 4))
m.show()

CPU times: user 13 ms, sys: 385 ms, total: 398 ms
Wall time: 5.68 s

Set region

# Set region
tools.g_region(vector="parcel", res=0.3)

Process surface

%%time

# Set mask
with gs.MaskManager():

    # Set mask
    tools.r_mask(vector="parcel")

    # Clip
    tools.r_mapcalc(expression="clip = surface")
    
    # Fill nulls
    tools.r_fillnulls(input="clip", output="surface", method="rst", tension=40)

    # Remove temporary map
    tools.g_remove(type="raster", name="clip", flags="f")

    # Set color gradient
    tools.r_colors(map="surface", color="viridis")

# Visualize
m = gj.Map(width=1600)
m.d_rast(map="surface")
m.d_legend(raster="surface", color="white", at=(75, 95, 2, 4))
m.show()

CPU times: user 16.1 ms, sys: 593 ms, total: 610 ms
Wall time: 58.2 s

Import lidar

%%time

# Import lidar
tools.v_in_pdal(
    input="hilltop.laz",
    output="lidar",
    flags="ro"
    )

# Print point count
info = tools.v_info(map="lidar", format="json", flags="gt").json
print(info["points"])

15130452
CPU times: user 5.17 ms, sys: 141 ms, total: 146 ms
Wall time: 58.8 s

Decimate lidar

%%time

# Set region
tools.g_region(vector="lidar", res=0.3)

# Decimate lidar
tools.v_decimate(
    input="lidar",
    output="points",
    layer=1,
    zdiff=1.0,
    flags="gz"
    )

# Print point count
info = tools.v_info(map="points", format="json", flags="gt").json
print(info["points"])

1813159
CPU times: user 7.98 ms, sys: 359 ms, total: 367 ms
Wall time: 1min 30s

Classify ground

# Install addon
tools.g_extension(extension="v.lidar.mcc")

%%time

# Classify lidar
tools.v_lidar_mcc(
    input="points",
    ground="ground",
    nonground="nonground",
    nl=3,
    t=0.2,
    j=0.1,
    f=2,
    s=15
    )

CPU times: user 4.75 ms, sys: 72.6 ms, total: 77.3 ms
Wall time: 1min 27s

Interpolate ground

%%time

# Interpolate ground
tools.v_surf_rst(
    input="ground",
    elevation="interpolation",
    nprocs=12
    )

# Clip ground
with gs.MaskManager():

    # Set mask
    tools.r_mask(vector="parcel")   

    # Clip to region
    tools.r_mapcalc(expression="elevation = interpolation")
    

# Set color gradient
tools.r_colors(map="elevation", color="viridis")

# Visualize
m = gj.Map(width=1600)
m.d_rast(map="elevation")
m.d_legend(raster="elevation", color="white", at=(75, 95, 2, 4))
m.show()

CPU times: user 20.2 ms, sys: 665 ms, total: 685 ms
Wall time: 2min 3s

Calculate canopy height

# Calculate difference
tools.r_mapcalc(expression="canopy = surface - elevation")

# Set color gradient
tools.r_colors(map="canopy", color="viridis")

# Visualize
m = gj.Map(width=1600)
m.d_rast(map="canopy")
m.d_legend(raster="canopy", color="white", at=(75, 95, 2, 4))
m.show()