<jupyter_start><jupyter_text>Forest structure using PDAL + PythonDr <NAME>, November 2019.This work is a set of python modules to replace MATLAB code for generating TERN forest metrics from airborne LIDAR. Fundamental ideasExisting code uses a series of nested loops, meaning we can't take advantage of array operations or easily reformat or paralellise functionalityThe approach used here defines a transportable function for each TERN product. These are applied to the data using a single loop (which could be chunked and parallelised).A simple process step-through looks like:1. Read LAS tile using PDAL. This removes an uncompression step. It also removes low outliers and computes normalised height for each point on the fly2. Read numpy labelled arrays from PDAL output into a GeoPandas dataframe, and apply a 2D spatial index3. From LAS file metadata, produce a fishnet grid with cells of size 'output resolution X output resolution'4. Iterate over grid cells, select valid points and generate TERN products for each grid cell5. Assemble an output array for each TERN product and write to GeoTIFFThis set of functions operates per-las-tile. An additional layer may be added to merge mutliple raster outputs into larger datasets to do:- snake_casify variable names<jupyter_code>NODATA_VALUE = -9999
#imports
import pdal
import numpy as np
import json

from shapely.geometry import Point
from shapely.geometry import MultiPolygon
from shapely.geometry import box
#from shapely.strtree import STRtree

import geopandas as gpd
import pandas as pd
import osmnx as ox

import os

# not using this, using geopandas instead
from rtree import index

# this is needed to create a raster from the output array
from osgeo import gdal
import osgeo.osr as osr
def writegeotiff(griddedpoints, outfile, parameters):
    """
    writes out a geotiff from a numpy array of forest metric
    results.
    
    inputs:
    - a numpy array of metrics [griddedpoints]
    - an outfile name [outfile]
    - a dictionary of parameters for the raster
    
    outputs:
    - a gdal dataset object
    - [outfile] written to disk
    """
    
    width = parameters["width"]
    height = parameters["height"]
    
    drv = gdal.GetDriverByName("GTiff")
    ds = drv.Create(outfile, width, height, 6, gdal.GDT_Float32)
    ds.SetGeoTransform(parameters["upperleft_y"],
                       parameters["resolution"],
                       0,
                       parameters["upperleft_y"],
                       0,
                       parameters["resolution"])
    ds.setProjection = parameters["projection"]
    ds.GetRasterBand(1).WriteArray(arr)

    return(ds)

def pdal2df(points):
    """
    Feed me a PDAL pipeline return array, get back a 
    GeoPandas dataframe 
    """

    arr = points[0]
    description = arr.dtype.descr
    cols = [col for col, __ in description]
    gdf = gpd.GeoDataFrame({col: arr[col] for col in cols})
    gdf.name = 'nodes'
    gdf['geometry'] = gdf.apply(lambda row: Point((row['X'], row['Y'])), axis=1)
    
    return(gdf_nodes)

def spatialindex(dataframe):
    sindex = dataframe.sindex
    return(sindex)

#get a pointview from PDAL
def readlasfile(lasfile):
    """
    Run a PDAL pipeline. Input is a JSON declaration to 
    deliver to PDAL. Output is a labelled numpy array.
    
    Data are filtered to:
    - label local minima as noise
    - compute height above ground using nearest ground point
      neighbours (TIN method arriving soon)
    - sort using a morton order (space filling curve) to 
      speed indexing later.
    
    """
    pipeline = {
        "pipeline": [
            {
                "type": "readers.las",
                "filename": lasfile
            },
            {
                "type": "filters.hag"
            }
        ]
    }
    
    pipeline = pdal.Pipeline(json.dumps(pipeline))
    count = pipeline.execute()
    
    #read points into labelled arrays
    arrays = pipeline.arrays

    #return a numpy array to operate on
    return(arrays)

def readlasmetadata(lasfile):
    pipeline = {
        "pipeline": [
            {
                "type": "readers.las",
                "filename": lasfile,
                "count": 1
            },
            {
                "type": "filters.info"
            }
        ]
    }
    
    pipeline = pdal.Pipeline(json.dumps(pipeline))
    pipeline.validate()
    pipeline.loglevel = 2  # stay quiet
    count = pipeline.execute()
    
    #extract metadata into a JSON blob
    metadata = json.loads(pipeline.metadata)
    
    return(metadata)


def readlaschunk(lasfile, poly):
    """
    Run a PDAL pipeline to collect a chunk of points from a LASfile.
    
    Inputs
    - a PDAL readable file name
    - a geometry bounding the chunk of data to be clipped
    
    Output
    - a labelled numpy array.

    
    """
    xmin = str(np.int(np.floor(poly.bounds[0])))
    ymin = str(np.int(np.floor(poly.bounds[1])))
    xmax = str(np.int(np.ceil(poly.bounds[2])))
    ymax = str(np.int(np.ceil(poly.bounds[3])))
    
    cropbox = "([" + xmin + "," + ymin + "],[" + xmax + "," + ymax + "])"
    print(cropbox)
    
    pipeline = {
        "pipeline": [
            {
                "type": "readers.las",
                "filename": lasfile
            },
            {
                "type": "filters.crop",
                "bounds": cropbox
            }
        ]
    }
    
    pipeline = pdal.Pipeline(json.dumps(pipeline))
    count = pipeline.execute()

    #read points into labelled arrays
    array = pipeline.arrays

    #return a numpy array to operate on
    return(array)


def extract_vars(df):
    """
    extract relevant variables
    do we need to do this now? or wait till we've grabbed the indexed chunk?
    lets write it anyway, then the index chunkifier can call it...
    
    inputs:
    - a numpy labelled array resulting from a PDAL LAS/LAZ file read
    
    outputs:
    - 1D arrays containing relevant variables
    
    """
    classification = df["Classification"].values
    intensity = df["Intensity"].values
    returnnumber = df["ReturnNumber"].values
    numberofreturns = df["NumberOfReturns"].values
    elevation = df["Z"].values
    hag = df["HeightAboveGround"].values

    return(intensity, returnnumber, numberofreturns, elevation, hag)


def gen_raster_cells(metadata, resolution):
    """
    Generate cells of 'resolution x resolution' for point querying
    
    input:
    - PDAL metadata
    
    output:
    - shapely geometry containing polygons defining 'resolution x resolution'
      boxes covering the LAS tile extent
      
    """
    bbox = box(metadata["metadata"]["readers.las"][0]["minx"],
               metadata["metadata"]["readers.las"][0]["miny"],
               metadata["metadata"]["readers.las"][0]["maxx"],
               metadata["metadata"]["readers.las"][0]["maxy"])
    
    tiledBBox = ox.quadrat_cut_geometry(bbox, quadrat_width=resolution)
    
    return(tiledBBox)

def get_cell_points(poly, df, sindex):
    
    poly = poly.buffer(1e-14).buffer(0)
    possible_matches_index = list(sindex.intersection(poly.bounds))
    possible_matches = df.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(poly)]
    
    return(precise_matches)
# Vegetation cover fraction: (Nfirst - Nsingle) / Nfirst
def comp_vcf(points):
    """
    Computes vegetation cover fraction according to the TERN product manual.
    
    inputs:
    - a labelled array of points from an input LAS tile
    
    outputs:
    - a numpy array of grid cells containing the result of:
    
    (Nfirst - Nsingle) / Nfirst
    
    ...where:
    Nfirst = count of first returns
    Nsingle = count of single returns
    
    ...per grid cell.
    """
    # collect all the first and single return indices
    nSingle = np.size(np.where(points["NumberOfReturns"].values == 1))
    nFirst = np.size(np.where(points["ReturnNumber"].values == 1))
    if (nFirst > 0):
        vcf = (nFirst - nSingle) / nFirst
    else:
        print('no first returns, set vcf to {}'.format(NODATA_VALUE))
        vcf = -9999
        
    return(vcf)
# Canopy layering index:

# R = total returns
# 
# vegetation layer cover fraction: LCF

def comp_lcf(points, heights, vcf):
    """
    Compute LCF as per the TERN product manual:
    
    LCF = VCF * (((veg returns below H2) - (veg returns below H1)) / (veg returns below H2))
    
    Inputs:
    - a set of points to compute LCF over
    - a height threshold pair, containing H1 and H2 as an array [h1, h2]
    - a precomputed VCF
    
    Outputs:
    - a floating point number denoting LCF
    
    Conditions:
    
    The LCF *must* be computed over the same set of points as the VCF used as input.
    
    """
    
    h1 = heights[0]
    h2 = heights[1]
    
    #find veg returns - ASPRS classes 3,4,5
    veg_returns = np.where(np.logical_or(points["Classification"].values == 3,
                             points["Classification"].values == 4,
                             points["Classification"].values == 5))
    # how many veg returns have height below the first threshold?
    vegbelowh1 = np.size(np.where(points["HeightAboveGround"][vegreturns] < h1))
    
    # how many veg returns have height below the second threshold?
    vegbelowh2 = np.size(np.where(points["HeightAboveGround"][vegreturns] < h2))
    
    # compute the LCF
    lcf = vcf * ( (vegbelowh2 - vegbelowh1) / vegbelowh2)
    
    return(lcf)

#CTH
def comp_cth(points):
    # compute the highest vegetation point in each grid cell
    
    veg_returns = np.where(np.logical_or(points["Classification"].values == 3,
                             points["Classification"].values == 4,
                             points["Classification"].values == 5)) 

    vegpoints = points["HeightAboveGround"].values[veg_returns]

    cth = np.max(vegpoints)

    return(cth)
def comp_dem(points):
    # interpolate ground returns in a grid and output a raster
    
    
    return()
def comp_fbf(points):
    # if building classes exist, compute a fractional conver per grid cell...
    
    return()
def read_data(lasfile):
    """
    wrapper to read in LAS data and produce a dataframe + spatial index
    """
    metadata, points = readlasfile(lasfile)
    
    dataframe, spatial_index = pdal2df(points)
    
    return(metadata, dataframe, spatial_index)
    
def compute_tern_products(metadata, points, sindex, resolution):
    """
    Wrapper to iterate over the input data and generate rasters for each product.
    
    *note this part could be paralellised - maybe per-product, or per-cell
    
    Each grid square processed in this loop corresponds to one pixel in an output raster.
    
    """
    
    #set up an 'output resolution' sized grid - like a fishnet grid.
    # each polygon in the resulting set covers an area of 'resolution X resolution'
    pixel_grid = gen_raster_cells(metadata, resolution)
    
    #set up output rasters
    
    # get tile width and height
    tile_width = metadata["metadata"]["readers.las"][0]["maxx"] - metadata["metadata"]["readers.las"][0]["minx"]
    tile_height = metadata["metadata"]["readers.las"][0]["maxy"] - metadata["metadata"]["readers.las"][0]["miny"]

    raster_xsize = int(np.ceil(tile_width) / resolution)
    raster_ysize = int(np.ceil(tile_height) / resolution)
    
    print(tile_width)
    print(raster_xsize)
    
    vcf_raster = np.zeros((raster_xsize, raster_ysize))
    
    print(np.shape(vcf_raster))

    lcf_raster = np.zeros((raster_xsize, raster_ysize))
    cth_raster = np.zeros((raster_xsize, raster_ysize))
    
    
    for pixel in pixel_grid:
        
        #compute output array index for this cell:
        
        poly_x, poly_y = pixel.centroid.xy
        
        poly_base_x = poly_x[0] - metadata["metadata"]["readers.las"][0]["minx"]
        poly_base_y = poly_y[0] - metadata["metadata"]["readers.las"][0]["miny"]
        
        print(poly_base_x)
        print(poly_base_y)
        
        array_x = int(np.floor((poly_base_x / (resolution)) ))
        array_y = int(np.floor((poly_base_y / (resolution)) ))
        
        #print('array X: {}; array Y: {}'.format(array_x, array_y))
        
        #get points for this cell
        matches = get_cell_points(pixel, points, sindex)
        
        #compute in order
        #VCF
        vcf_raster[array_x, array_y] = comp_vcf(matches)
        
        #LCF - need stuff about levels here...
        #lcf_raster[array_x, array_y] = comp_lcf(points)
        
        #CTH
        try:
            cth_raster[array_x, array_y] = comp_cth(matches)
        except ValueError:
            print('no vegetation returns were present, CTH set to {} for array index {} {}'.format(NODATA_VALUE, array_x, array_y))
            cth_raster[array_x, array_y] = NODATA_VALUE
            
    #end of computing stuff
    
    #extract EPSG code from LAS:

    
    return(vcf_raster, cth_raster)<jupyter_output><empty_output><jupyter_text>Testing functionality using a local fileThe following section generates metrics from a local LAZ file. Plugging in download mechanics from ELVIS will be added later<jupyter_code>#lidar test file - Mt Ainslie, chosen for varied vegetation cover and topography
# this is pretty big, may need splitting up!

#lasfile = "/Volumes/Antares/ACT-lidar/8ppm/callingelvis-testdata/ACT2015_8ppm-C3-AHD_6966094_55.las"

lasfile = "/Volumes/Antares/fire-test/NSW Government - Spatial Services-2/Point Clouds/AHD/StAlbans201709-LID2-C3-AHD_2866308_56_0002_0002/StAlbans201709-LID2-C3-AHD_2866308_56_0002_0002.las"
lasfile = "./Berridale201802-LID2-C3-AHD_6585974_55_0002_0002.las"
%%time

metadata = readlasmetadata(lasfile)
metadata
resolution = 25
pixel_grid = gen_raster_cells(metadata, resolution)
pixel_grid[0]
boxbounds = pixel_grid[0].bounds[0]
boxbounds
#get a pointview from PDAL
def readlaschunk(lasfile, poly):
    """
    Run a PDAL pipeline to collect a chunk of points from a LASfile.
    
    Inputs
    - a PDAL readable file name
    - a geometry bounding the chunk of data to be clipped
    
    Output
    - a labelled numpy array.

    
    """
    xmin = str(np.int(np.floor(poly.bounds[0])))
    ymin = str(np.int(np.floor(poly.bounds[1])))
    xmax = str(np.int(np.ceil(poly.bounds[2])))
    ymax = str(np.int(np.ceil(poly.bounds[3])))
    
    cropbox = "([" + xmin + "," + ymin + "],[" + xmax + "," + ymax + "])"
    print(cropbox)
    
    pipeline = {
        "pipeline": [
            {
                "type": "readers.las",
                "filename": lasfile
            },
            {
                "type": "filters.crop",
                "bounds": cropbox
            }
        ]
    }
    
    pipeline = pdal.Pipeline(json.dumps(pipeline))
    count = pipeline.execute()

    #read points into labelled arrays
    array = pipeline.arrays

    #return a numpy array to operate on
    return(array)
%%time

laschunk = readlaschunk(lasfile, pixel_grid[0])
laschunk
%%time

laspoints = readlasfile(lasfile)
%%time

#here we read points into a GeoDataFrame and dump the labelled array.
# this is a pretty expensive step RAM wise, we're duplicating all the points...

df = pdal2df(laspoints)
# set points to None, we don't use them anymore
laspoints = None
%%time

sindex = spatialindex(df)
# dump everything from memory
points = None
df = None
vcf_raster = None
cth_raster = None
%%time
# this part of the process is simply reading from the source file. No analysis yet.

metadata, points = readlasfile(lasfile)
%%time

#here we read points into a GeoDataFrame and dump the labelled array.
# this is a pretty expensive step RAM wise, we're duplicating all the points...

df = pdal2df(points)

# set points to None, we don't use them anymore
points = None
%%time

# here we generate an RTree index on the dataframe using GeoPandas.
# also pretty expensive... 

sindex = spatialindex(df)
%%time
## rtree index building straight from the point dataset...

idx = index.Index()
for pid, point in enumerate(points[0]):
    idx.insert(pid, (point[0], point[1],point[0], point[1]), point)
# set an output resolution

resolution = 100
%%time

vcf_raster, cth_raster = compute_tern_products(metadata, df, sindex, resolution)
from matplotlib import pyplot as plt

%matplotlib inline
plt.imshow(vcf_rasters)
plt.colorbar()
plt.imshow(cth_rasters)
plt.colorbar()
wktcrs = metadata["metadata"]["readers.las"][0]["comp_spatialreference"]
type(wktcrs)
srs = osr.SpatialReference()
srs.SetFromUserInput("EPSG:28356")
srs.ImportFromWkt(wktcrs)
srs.GetAuthorityCode(None)
srs.GetAuthorityName(None)
tiledBBox = gen_raster_cells(metadata,resolution)
tiledBBox[50]
%%time
idxpoints, tree = create_spatial_index(points)
%%time
thepoints = tree.query(geometry_cut[50])
query<jupyter_output><empty_output><jupyter_text>Run a sample workflow on one square<jupyter_code>pixel_grid[0]
matches = get_cell_points(pixel_grid[10],df, sindex)
def comp_cth1(points):
    # compute the highest vegetation point in each grid cell
    
    veg_returns = np.where(np.logical_or(points["Classification"].values == 3,
                             points["Classification"].values == 4,
                             points["Classification"].values == 5)) 

    vegpoints = points["HeightAboveGround"].values[veg_returns]

    cth = np.max(vegpoints)
    
    return(cth)
comp_cth1(matches)<jupyter_output><empty_output><jupyter_text>create a dataframe for pretty querying purposes..<jupyter_code>points
#import pandas as pd




#points = None
# find the points that intersect with each subpolygon and add them to points_within_geometry
points_within_geometry = pd.DataFrame()
def get_points(poly, df, sindex):
    poly = poly.buffer(1e-14).buffer(0)
    possible_matches_index = list(sindex.intersection(poly.bounds))
    possible_matches = df.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(poly)]
    
    return(precise_matches)
polyX, polyY = poly.centroid.xy
polyX
polyY
polyBaseXCoord = polyX[0] - metadata["metadata"]["readers.las"][0]["minx"]
arrayXindex = (polyBaseXCoord / (resolution/2 )) -1
arrayXindex
polyBaseYCoord = polyY[0] - metadata["metadata"]["readers.las"][0]["miny"]
arrayYindex = (polyBaseYCoord / (resolution/2 )) - 1 
arrayYindex
arrayX = PolyBaseCoord - numberofcells
polyY[0] - metadata["metadata"]["readers.las"][0]["miny"]
tilewidth = metadata["metadata"]["readers.las"][0]["maxx"] - metadata["metadata"]["readers.las"][0]["minx"]
tilewidth
tileheight = metadata["metadata"]["readers.las"][0]["maxy"] - metadata["metadata"]["readers.las"][0]["miny"]
tileheight
vcfRaster = np.zeros((tilewidth, tileheight))
vcfRaster
%%time 

matches = get_points(poly, df)
matches
intensity, returnnumber, numberofreturns, elevation, hag  = extract_vars(matches)
np.size(np.where(matches["NumberOfReturns"].values == 1))
## OK now we can make magic - extracting each grid cell, we can rasterify it...
vcf = comp_vcf(matches)
vcf
vcfRaster[0,138] = vcf
vcfRaster
this = matches["Classification"].values
veg = np.where(np.logical_or(matches["Classification"].values == 3,matches["Classification"].values == 4, matches["Classification"].values == 5)) 
this = matches["Classification"].values
np.where(matches["HeightAboveGround"].values[veg] < 1)
lcf0105 = comp_lcf(matches, [1, 5], vcf)<jupyter_output><empty_output><jupyter_text>code purgatorystuff here might be useful, or not<jupyter_code># this will likely evolve to take a 'what to grid' input
# actually likely not needed...
def grid_setup(pointmetadata, resolution):
    """
    Sets up array indexes for an incoming las file, using an input resolution
    
    input:
    - las file metadata
    - a scalar resolution
    
    output:
    - a numpy array of indexing values to divide the input points into
      'resolution' x 'resolution' bins.
      
    implicit assumptions:
    - 'resolution' is always set in native LAS file units
    """
    
    return()
# not used yet - geopandas is doing this part!
def create_spatial_index(arrays):
    """
    task here is to map geospatial space to numpy array space.
    
    - using the point metadata set up a grid of [resolution] x [resolution]
    - scan the point coordinates array to see which points live in which cell
    - create an index which maps point array indexes to grid indexes
    
    idly wondering if a fast point-in-polygon does this job. Shapely / OGR to the rescue?
    ...this would be parallelisable...
    """
    points = []
    
    for thepoint in arrays[0]:
        #print(thepoint[0])
        points.append(Point(thepoint[0], thepoint[1]))
    
    tree = STRtree(points)
    
    return(points, tree)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1033518205931999, 'fraction_numerical': 0.011960453339763685, 'mean_word_length': 2.9414559969587533, 'pattern_counts': {'":': 19, '<': 25, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '56247', 'license': 'Apache-2.0', 'n_tokens_mistral': 6644, 'n_tokens_neox': 6411, 'n_words': 1861, 'path': 'notebooks/forest-structure-pdal-crop-approach.ipynb', 'repo_name': 'adamsteer/forestmetrics', 'size': '20739'}