Back to Project
Methodology: Mapping Dominant Tree Species in the Petawawa Research Forest with LiDAR

Detailed code and steps used in the analysis

1 Data Loading and Preparation

Loaded plot locations from a GeoPackage using geopandas, dropping unnecessary columns. Species codes were loaded from a CSV and plot-level dominant species data from an Excel file. The 67-band ALS metrics raster was opened using rasterio and metric values extracted at each plot location.

import geopandas as gpd
import pandas as pd
import numpy as np
import rasterio
from rasterio.transform import rowcol
import matplotlib.pyplot as plt
from sklearn.feature_selection import mutual_info_classif

# Load plot locations
plots = gpd.read_file("plot_locations.gpkg")
plots = plots.drop(columns=["Source", "Date"])

# Load species codes
sp_codes = pd.read_csv("sp_codes.csv")
sp_codes = sp_codes.rename(columns={"Name": "Species", "MNR Code": "Code"})

# Load plot-level dominant species data
plot_data = pd.read_excel("plot_data.xlsx")

# Open ALS raster and extract values at plot locations
with rasterio.open("als_metrics.tif") as src:
    transform = src.transform
    crs = src.crs
    bands = src.count
    band_names = [src.descriptions[i] for i in range(bands)]

    lidar_values = []
    for _, row in plots.iterrows():
        x, y = row.geometry.x, row.geometry.y
        r, c = rowcol(transform, x, y)
        try:
            vals = src.read(window=((r, r+1), (c, c+1))).flatten()
        except Exception:
            vals = np.full(bands, np.nan)
        lidar_values.append(vals)

lidar_df = pd.DataFrame(lidar_values, columns=band_names)

2 Data Cleaning

Rows with missing LiDAR values were identified using isnull() and removed from both the plot data and LiDAR dataframes. Species with fewer than 15 observations were filtered out, retaining only White Pine, Red Pine, and Red Oak.

# Identify and remove missing LiDAR values
valid_mask = ~lidar_df.isnull().any(axis=1)
lidar_df = lidar_df[valid_mask].reset_index(drop=True)
plot_data = plot_data[valid_mask].reset_index(drop=True)

# Merge plot data with species codes
merged = plot_data.merge(sp_codes, on="Code", how="left")

# Filter species with fewer than 15 observations
species_counts = merged["Species"].value_counts()
valid_species = species_counts[species_counts >= 15].index
filtered = merged[merged["Species"].isin(valid_species)].reset_index(drop=True)
lidar_filtered = lidar_df[merged["Species"].isin(valid_species)].reset_index(drop=True)

print(f"Retained species: {list(valid_species)}")
print(f"Total plots after filtering: {len(filtered)}")

3 Feature Selection

Species classes were encoded numerically for use with scikit-learn. mutual_info_classif was applied to rank all 67 ALS metrics by their correlation with species class labels. The top 5 most informative features were printed.

from sklearn.preprocessing import LabelEncoder

# Encode species labels
le = LabelEncoder()
y = le.fit_transform(filtered["Species"])
X = lidar_filtered.values

# Compute mutual information scores
mi_scores = mutual_info_classif(X, y, random_state=42)
mi_series = pd.Series(mi_scores, index=lidar_filtered.columns)
mi_series = mi_series.sort_values(ascending=False)

print("Top 5 most informative ALS features:")
print(mi_series.head(5))

4 Data Exploration and Visualization

A bar chart was generated showing the number of ground plots per dominant species. Box plot panels were created for the top 3 mutual information features plus Vertical Complexity Index (vci_0.5bin) to visually assess class separability.

# Bar chart: number of sites per dominant species
all_species_counts = merged["Species"].value_counts()
plt.figure(figsize=(12, 5))
all_species_counts.plot(kind="bar", color="#1f77b4")
plt.title("Number of Sites per Dominant Species")
plt.xlabel("")
plt.ylabel("Count")
plt.xticks(rotation=45, ha="right")
plt.tight_layout()
plt.savefig("bar_graph_petawawa.png", dpi=150)
plt.show()

# Box plots for top features + vci_0.5bin
top_features = list(mi_series.head(3).index) + ["vci_0.5bin"]
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
fig.suptitle("Box plot of LiDAR metrics for each species")

for ax, feat in zip(axes.flatten(), top_features):
    groups = [lidar_filtered[feat][filtered["Species"] == sp].values
              for sp in le.classes_]
    ax.boxplot(groups, labels=le.classes_, patch_artist=True,
               medianprops=dict(color="orange"))
    ax.set_title(feat)

plt.tight_layout()
plt.savefig("box_plot_petawawa.png", dpi=150)
plt.show()

5 Conclusions and Recommendations