compute_carbon_credits.py

import numpy as np
import rasterio

# Function to compute pixel area from CHM raster
def compute_pixel_area(chm_path):
    """Compute pixel area from a raster file."""
    with rasterio.open(chm_path) as src:
        pixel_width, pixel_height = src.res
        pixel_area = pixel_width * pixel_height
    return pixel_area

# Function to compute carbon credits
def compute_carbon_credits(chm_path, BEF=0.5, CF=0.47, carbon_price_per_ton=50, a=0.05, b=2.5):
    """
    Compute carbon credits based on CHM raster file.
    
    Parameters:
    - chm_path: Path to the Canopy Height Model (CHM) raster file.
    - BEF: Biomass expansion factor (default is 0.5).
    - CF: Carbon fraction of biomass (default is 0.47).
    - carbon_price_per_ton: Carbon price in USD per tonne (default is 50).
    - a: Coefficient for AGB calculation (default is 0.05).
    - b: Exponent for AGB calculation (default is 2.5).
    
    Returns:
    - Estimated carbon stock in tonnes.
    - Potential carbon credits in USD.
    """
    # Compute pixel area for the raster
    pixel_area = compute_pixel_area(chm_path)
    print(f"Pixel area: {pixel_area} square meters")

    with rasterio.open(chm_path) as src:
        chm_data = src.read(1)  # Read the first band

        # Mask out no-data values
        chm_data = np.ma.masked_equal(chm_data, src.nodata)

    # Function to calculate AGB
    def calculate_agb(chm_data, pixel_area):
        """Estimate above-ground biomass (AGB) from the Canopy Height Model (CHM)."""
        agb = a * (chm_data ** b) * BEF  # Corrected to multiply by BEF
        return agb * pixel_area  # Convert to kg per pixel area

    # Function to calculate carbon stock
    def calculate_carbon_stock(agb, CF):
        """Convert above-ground biomass (AGB) to carbon stock."""
        return agb * CF / 1000  # Convert kg to tonnes

    # Function to calculate carbon credits
    def calculate_carbon_credits(carbon_stock):
        """Calculate the potential carbon credits in USD."""
        return carbon_stock * carbon_price_per_ton

    # Calculate AGB
    agb_data = calculate_agb(chm_data, pixel_area)
    
    # Calculate carbon stock
    carbon_stock = calculate_carbon_stock(agb_data.sum(), CF)
    
    # Calculate carbon credits
    carbon_credits = calculate_carbon_credits(carbon_stock)

    print(f"Estimated Carbon Stock: {carbon_stock:.2f} tonnes")
    print(f"Potential Carbon Credits: ${carbon_credits:.2f}")

    return carbon_stock, carbon_credits