coordinate_extraction.py

import laspy
import numpy as np
import os
import csv

# Set directory containing LAZ files
folder_path = r"S:\qfl_australia_tls\Curtain Fig\output"

# Output CSV file
output_csv = os.path.join(r"S:\qfl_australia_tls\Curtain Fig\\", "tree_locations.csv")

# Load already processed filenames if CSV exists
processed_files = set()
if os.path.exists(output_csv):
    with open(output_csv, mode="r", newline="") as file:
        reader = csv.DictReader(file)
        for row in reader:
            processed_files.add(row["filename"])

# Initialize CSV writing (append mode if already exists, else write mode)
with open(output_csv, mode="a" if os.path.exists(output_csv) else "w", newline="") as file:
    writer = csv.writer(file)

    # If new file, write the header
    if os.stat(output_csv).st_size == 0:
        writer.writerow(["filename", "center_x", "center_y"])

    # Loop through all LAS files
    for filename in os.listdir(folder_path):
        if filename.endswith(".las") and filename not in processed_files:
            file_path = os.path.join(folder_path, filename)

            try:
                # Load LAS file
                las = laspy.read(file_path, laz_backend=laspy.LazBackend.Lazrs)

                # Extract coordinates
                x = np.array(las.x)
                y = np.array(las.y)
                z = np.array(las.z)

                # Define cross-section range (e.g., bottom 0.5 meters)
                z_min = np.min(z)
                z_threshold = z_min + 0.5

                # Filter points in the cross-section
                mask = (z >= z_min) & (z <= z_threshold)
                x_section = x[mask]
                y_section = y[mask]

                if len(x_section) > 0:
                    # Compute centroid
                    center_x = np.mean(x_section)
                    center_y = np.mean(y_section)

                    # Write to CSV
                    writer.writerow([filename, center_x, center_y])
                    print(f"Processed {filename} -> Location: ({center_x}, {center_y})")
                else:
                    print(f"Skipped {filename} (no valid points in cross-section).")

            except Exception as e:
                print(f"Error processing {filename}: {e}")
        elif filename in processed_files:
            print(f"Skipping {filename} (already processed)")

print(f"Results saved to {output_csv}")