"""Hägglund site index models for Swedish Norway spruce and Scots pine.
Three separate publications, one per model. The ``Hagglund_1970`` class name is a
legacy identifier: there is no Hägglund (1970) site-index publication.
Sources:
Hägglund, B. (1972). *Om övre höjdens utveckling för gran i norra Sverige.*
Skogshögskolan, institutionen för skogsproduktion, Rapporter och uppsatser nr 21.
Hägglund, B. (1973). *Om övre höjdens utveckling för gran i södra Sverige.*
Skogshögskolan, institutionen för skogsproduktion, Rapporter och uppsatser nr 24.
Hägglund, B. (1974). *Övre höjdens utveckling i tallbestånd.* Skogshögskolan,
institutionen för skogsproduktion, Rapporter och uppsatser nr 31.
"""
import math
import warnings
from enum import Enum
from functools import wraps
from typing import Callable, Union, cast
from numpy import exp, log
from pyforestry.base.contracts import FormulaDescriptor, SourceReference
from pyforestry.base.helpers import Age, AgeMeasurement, SiteIndexValue, TreeSpecies
[docs]
class HagglundSpruceModel:
"""Height trajectory models for Norway spruce."""
[docs]
@staticmethod
def northern_sweden(
dominant_height: float,
age: Union[float, AgeMeasurement],
age2: Union[float, AgeMeasurement],
latitude: float,
culture: bool = True,
) -> tuple[SiteIndexValue, float]:
"""
Calculates the height trajectory of Norway Spruce in northern Sweden
based on Hägglund (1972).
In this version, we allow the input parameter 'age' to be provided either
as ``Age.DBH`` (the effective DBH age for the stand) or as ``Age.TOTAL``.
In the latter case we use Newton–Raphson to solve for the effective DBH
age ``x`` that satisfies ``x + T13(x) = age_total``.
For the prediction age 'age2':
- If provided as Age.DBH, we add the computed T13 to obtain the total age.
- Otherwise (if Age.TOTAL) it is used directly.
"""
# Extract numeric values from age and its code.
if isinstance(age, AgeMeasurement):
age_value = float(age)
age_type = age.code
else:
age_value = float(age)
age_type = Age.DBH.value # default if plain number
# For age2, extract numeric value and preserve type.
if isinstance(age2, AgeMeasurement):
age2_value = float(age2)
age2_type = age2.code
else:
age2_value = float(age2)
age2_type = Age.TOTAL.value # default to total if no code provided
P = 0.9175 if culture else 1.0
# Convert dominant height (meters) to decimeters and subtract 13 (the base value)
top_height_dm = dominant_height * 10 - 13
if age_value > (407 - 1.167 * top_height_dm):
warnings.warn("Too old stand, outside of the material.", stacklevel=2)
if latitude >= 67 or latitude <= 60:
warnings.warn(
"Outside of latitudinal range, 60°<= L <= 67° N, using function 8.4 ",
stacklevel=2,
)
# Parameters for function 8.4
B = 3.4501
C = 0.77518
D = -0.42579
E = 1.33935
else:
# Parameters for function 8.7
B = 3.3816
C = 0.77896
D = -1.24207 + 0.0014629 * latitude * 10
E = 1.25998
# Define the bonitering subroutine that uses the effective DBH age (eff_age)
def subroutineBonitering(eff_age: float):
"""Return productivity parameters for an effective DBH age."""
AI1 = 10.0
AI2 = 600.0
AI3 = 0.0
A2 = 0.0
RK = 0.0
RM2 = 0.0
while abs(AI1 - AI2) > 1:
AI3 = (AI1 + AI2) / 2.0
RK = 0.001936 + 0.00004100 * AI3**1.0105
A2 = B * AI3**C
RM2 = D + E / (0.56721 + 0.000008 * AI3**1.8008)
DIF = top_height_dm - A2 * (1 - exp(-eff_age * RK)) ** RM2
if DIF <= 0:
AI2 = AI3
else:
AI1 = AI3
T26 = (-1 / RK) * log(1 - (13 / A2) ** (1 / RM2))
T13_local = P * (7.0287 + 0.66118 * T26)
return A2, RK, RM2, T26, T13_local
# If age is provided as DBH, use it directly; otherwise (if TOTAL)
# solve for the effective DBH age.
if age_type == Age.DBH.value:
eff_age = age_value
else:
# Use Newton–Raphson to solve: f(x) = x + T13(x) - age_value = 0
def f(x):
"""Root function for DBH-age estimation."""
# Compute T13 from the bonitering subroutine for effective age x.
_, _, _, _, T13_local = subroutineBonitering(x)
return x + T13_local - age_value
def fprime(x, h=0.001):
"""Numerical derivative for the Newton solver."""
return (f(x + h) - f(x - h)) / (2 * h)
# Choose an initial guess. (Empirically, effective DBH age is lower than total age.)
x = age_value * 0.35
for _ in range(30):
fx = f(x)
fpx = fprime(x)
if abs(fpx) < 1e-8:
break
x_new = x - fx / fpx
if abs(x_new - x) < 1e-6:
x = x_new
break
x = x_new
eff_age = x
# Now compute productivity parameters at the effective DBH age.
A2, RK, RM2, T26, T13 = subroutineBonitering(eff_age)
if A2 > 336:
warnings.warn("Too high productivity, outside of material.", stacklevel=2)
if A2 < 189:
warnings.warn("Too low productivity, outside of material.", stacklevel=2)
# Determine the effective DBH age for the height prediction:
if age2_type == Age.DBH.value:
# If the target age 'age2' is already DBH, use its value directly for the formula.
effective_age2_dbh = age2_value
# The corresponding total age for output reference (if needed) would be:
output_total_age = age2_value + T13
else: # age2_type is TOTAL
# If the target age 'age2' is TOTAL, calculate the corresponding effective DBH age.
# We use the T13 calculated based on the input conditions (age, dominant_height).
effective_age2_dbh = age2_value - T13
# Check for potentially non-physical results (e.g., total age less than T13)
if effective_age2_dbh <= 0:
warnings.warn(
"Calculated effective DBH age for prediction "
f"({effective_age2_dbh:.2f}) is non-positive. "
f"This might happen if age2 (Total Age={age2_value}) is less "
f"than T13 ({T13:.2f}). Using a small positive value (e.g., 1.0) instead.",
stacklevel=2,
)
effective_age2_dbh = 1.0 # Use a small positive fallback or handle as an error
# The corresponding total age for output reference is the input total age itself:
output_total_age = age2_value
# Calculate height using the derived effective DBH age
height = (13 + A2 * (1 - exp(-effective_age2_dbh * RK)) ** RM2) / 10
return (
SiteIndexValue(
height,
reference_age=Age.TOTAL(output_total_age),
species={TreeSpecies.Sweden.picea_abies},
fn=Hagglund_1970.height_trajectory.picea_abies.northern_sweden,
),
T13,
)
[docs]
@staticmethod
def southern_sweden(
dominant_height: float,
age: Union[float, AgeMeasurement],
age2: Union[float, AgeMeasurement],
) -> tuple[SiteIndexValue, float]:
"""
Calculates the height trajectory of Norway Spruce in southern Sweden
based on Hägglund (1973).
As in the northern function, the input 'age' can be provided either as DBH or TOTAL.
If provided as TOTAL, we solve using Newton–Raphson for the effective DBH age.
For the prediction age 'age2', if provided as DBH then T13 is added.
"""
if isinstance(age, AgeMeasurement):
age_value = float(age)
age_type = age.code
else:
age_value = float(age)
age_type = Age.DBH.value
if isinstance(age2, AgeMeasurement):
age2_value = float(age2)
age2_type = age2.code
else:
age2_value = float(age2)
age2_type = Age.TOTAL.value
top_height_dm = dominant_height * 10 - 13
def subroutineBonitering(eff_age: float):
"""Return productivity parameters for an effective DBH age."""
AI1 = 10.0
AI2 = 600.0
A2 = 0.0
RK = 0.0
RM2 = 0.0
while abs(AI1 - AI2) > 1:
AI3 = (AI1 + AI2) / 2.0
RK = 0.042624 - 7.1145 / AI3**1.0068
A2 = 1.0017 * AI3**0.99808
RM = 0.15933 + 3.7e6 / AI3**3.156
if RM > 0.95:
RM = 0.95
RM2 = 0.98822 / (1 - RM)
if RK < 0.0001:
RK = 0.0001
DIF = top_height_dm - A2 * (1 - exp(-eff_age * RK)) ** RM2
if DIF <= 0:
AI2 = AI3
else:
AI1 = AI3
T26 = (-1 / RK) * log(1 - (13 / A2) ** (1 / RM2))
T13_local = 4.9546 + 0.63934 * T26 + 0.031992 * T26 * T26
return A2, RK, RM2, T26, T13_local
if age_type == Age.DBH.value:
eff_age = age_value
else:
def f(x):
"""Root function for DBH-age estimation."""
_, _, _, _, T13_local = subroutineBonitering(x)
return x + T13_local - age_value
def fprime(x, h=0.001):
"""Numerical derivative for the Newton solver."""
return (f(x + h) - f(x - h)) / (2 * h)
x = age_value * 0.35
for _ in range(30):
fx = f(x)
fpx = fprime(x)
if abs(fpx) < 1e-8:
break
x_new = x - fx / fpx
if abs(x_new - x) < 1e-6:
x = x_new
break
x = x_new
eff_age = x
A2, RK, RM2, T26, T13 = subroutineBonitering(eff_age)
if A2 > 400:
warnings.warn("Too high productivity, outside of the material.", stacklevel=2)
if A2 < 250:
warnings.warn("Too low productivity, outside of the material.", stacklevel=2)
if A2 > 375 and top_height_dm > 267:
warnings.warn("Too old stand, outside of the material.", stacklevel=2)
if eff_age > 90:
warnings.warn("Too old stand, outside of the material.", stacklevel=2)
# Determine the effective DBH age for the height prediction:
if age2_type == Age.DBH.value:
# If the target age 'age2' is already DBH, use its value directly for the formula.
effective_age2_dbh = age2_value
# The corresponding total age for output reference (if needed) would be:
output_total_age = age2_value + T13
else: # age2_type is TOTAL
# If the target age 'age2' is TOTAL, calculate the corresponding effective DBH age.
# We use the T13 calculated based on the input conditions (age, dominant_height).
effective_age2_dbh = age2_value - T13
# Check for potentially non-physical results (e.g., total age less than T13)
if effective_age2_dbh <= 0:
warnings.warn(
"Calculated effective DBH age for prediction "
f"({effective_age2_dbh:.2f}) is non-positive. "
f"This might happen if age2 (Total Age={age2_value}) is less "
f"than T13 ({T13:.2f}). Using a small positive value (e.g., 1.0) instead.",
stacklevel=2,
)
effective_age2_dbh = 1.0 # Use a small positive fallback or handle as an error
# The corresponding total age for output reference is the input total age itself:
output_total_age = age2_value
# Calculate height using the derived effective DBH age
height = (13 + A2 * (1 - exp(-effective_age2_dbh * RK)) ** RM2) / 10
return (
SiteIndexValue(
height,
reference_age=Age.TOTAL(output_total_age),
species={TreeSpecies.Sweden.picea_abies},
fn=Hagglund_1970.height_trajectory.picea_abies.southern_sweden,
),
T13,
)
[docs]
class HagglundPineRegeneration(Enum):
"""Regeneration method categories for the pine model."""
CULTURE = "culture"
NATURAL = "natural"
UNKNOWN = "unknown"
def __str__(self):
"""Return the enum value as a string."""
return self.value
[docs]
class HagglundPineModel:
"""Height trajectory model for Scots pine."""
[docs]
@staticmethod
def sweden(
dominant_height_m: float,
age: Union[float, AgeMeasurement],
age2: Union[float, AgeMeasurement],
regeneration: HagglundPineRegeneration,
) -> tuple[SiteIndexValue, float]:
"""
Hägglund 1974: Height growth of Scots Pine in Sweden.
The forward-only solver is used.
- For ``age``: if provided as ``Age.TOTAL``, we solve using
Newton–Raphson for the effective DBH age.
- For ``age2``: if ``Age.DBH`` is provided, ``T13`` is added.
"""
if not isinstance(regeneration, HagglundPineRegeneration):
raise TypeError("regeneration argument must be of type HagglundPineRegeneration")
if isinstance(age, AgeMeasurement):
age_value = float(age)
age_type = age.code
else:
age_value = float(age)
age_type = Age.DBH.value
if isinstance(age2, AgeMeasurement):
age2_value = float(age2)
age2_type = age2.code
else:
age2_value = float(age2)
age2_type = Age.TOTAL.value
top_height_dm = dominant_height_m * 10 - 13
if age_value > 120:
warnings.warn("Too old stand, outside of the material.", stacklevel=2)
def subroutineBonitering(eff_age: float):
"""Return productivity parameters for an effective DBH age."""
AI1 = 10.0
AI2 = 600.0
A2 = 0.0
RK = 0.0
RM2 = 0.0
while abs(AI1 - AI2) > 1:
AI3 = (AI1 + AI2) / 2.0
RM = 0.066074 + 4.4189e5 / AI3**2.9134
RM = min(RM, 0.95)
RM2 = 1.0 / (1 - RM)
RK = 1.0002e-4 + 9.5953 * AI3**1.3755 / 1e6
RK = max(RK, 0.0001)
A2 = 1.0075 * AI3
DIF = top_height_dm - A2 * (1 - math.exp(-eff_age * RK)) ** RM2
if DIF <= 0:
AI2 = AI3
else:
AI1 = AI3
T26 = (-1 / RK) * math.log(1 - (13 / A2) ** (1 / RM2))
T262 = T26**2
if regeneration == HagglundPineRegeneration.NATURAL:
T13_local = 7.4624 + 0.11672 * T262
elif regeneration == HagglundPineRegeneration.UNKNOWN:
T13_local = 6.8889 + 0.12405 * T262
elif regeneration == HagglundPineRegeneration.CULTURE:
T13_local = 7.4624 + 0.11672 * T262 - 0.39276 * T26
else: # pragma: no cover - defensive
raise ValueError("Unexpected regeneration class.")
T13_local = min(T13_local, 50)
return A2, RK, RM2, T13_local
if age_type == Age.DBH.value:
eff_age = age_value
else:
def f(x):
"""Root function for DBH-age estimation."""
_, _, _, T13_local = subroutineBonitering(x)
return x + T13_local - age_value
def fprime(x, h=0.001):
"""Numerical derivative for the Newton solver."""
return (f(x + h) - f(x - h)) / (2 * h)
x = age_value * 0.35
for _ in range(30):
fx = f(x)
fpx = fprime(x)
if abs(fpx) < 1e-8:
break
x_new = x - fx / fpx
if abs(x_new - x) < 1e-6:
x = x_new
break
x = x_new
eff_age = x
A2, RK, RM2, T13 = subroutineBonitering(eff_age)
if A2 > 311:
warnings.warn("Too high productivity, outside of the material.", stacklevel=2)
if A2 < 180:
warnings.warn("Too low productivity, outside of the material.", stacklevel=2)
if A2 > 250 and eff_age > 100:
warnings.warn("Too old stand, outside of material.", stacklevel=2)
# Determine the effective DBH age for the height prediction:
if age2_type == Age.DBH.value:
# If the target age 'age2' is already DBH, use its value directly for the formula.
effective_age2_dbh = age2_value
# The corresponding total age for output reference (if needed) would be:
output_total_age = age2_value + T13
else: # age2_type is TOTAL
# If the target age 'age2' is TOTAL, calculate the corresponding effective DBH age.
# We use the T13 calculated based on the input conditions (age, dominant_height).
effective_age2_dbh = age2_value - T13
# Check for potentially non-physical results (e.g., total age less than T13)
if effective_age2_dbh <= 0:
warnings.warn(
"Calculated effective DBH age for prediction "
f"({effective_age2_dbh:.2f}) is non-positive. "
f"This might happen if age2 (Total Age={age2_value}) is less "
f"than T13 ({T13:.2f}). Using a small positive value (e.g., 1.0) instead.",
stacklevel=2,
)
effective_age2_dbh = 1.0 # Use a small positive fallback or handle as an error
# The corresponding total age for output reference is the input total age itself:
output_total_age = age2_value
# Calculate height using the derived effective DBH age
height = (13 + A2 * (1 - exp(-effective_age2_dbh * RK)) ** RM2) / 10
# Return the SiteIndexValue, referencing the corresponding TOTAL age.
# This standardizes the output reference age type.
return (
SiteIndexValue(
height,
reference_age=Age.TOTAL(output_total_age), # Always reference total age
species={TreeSpecies.Sweden.pinus_sylvestris},
fn=Hagglund_1970.height_trajectory.pinus_sylvestris.sweden,
),
T13,
)
# =============================================================================
# Container wrappers to select return value
# =============================================================================
[docs]
class HeightTrajectoryWrapper:
"""
Wrapper that calls the underlying model function and returns only the SiteIndexValue.
"""
def __init__(self, model):
"""Wrap a model class to return SiteIndexValue results."""
self._model = model
def __getattr__(self, name):
"""Proxy attributes, converting callables to SI-only wrappers."""
attr = getattr(self._model, name)
if callable(attr):
model_attr = cast(Callable[..., tuple[SiteIndexValue, float]], attr)
@wraps(model_attr)
def wrapper(*args, **kwargs):
"""Return only the SiteIndexValue from the model call."""
si_value, _ = model_attr(*args, **kwargs)
return si_value
return wrapper
return attr
[docs]
class TimeToBreastHeightWrapper:
"""
Wrapper that calls the underlying model function and returns only the T13 value.
"""
def __init__(self, model):
"""Wrap a model class to return T13 results."""
self._model = model
def __getattr__(self, name):
"""Proxy attributes, converting callables to T13-only wrappers."""
attr = getattr(self._model, name)
if callable(attr):
model_attr = cast(Callable[..., tuple[SiteIndexValue, float]], attr)
@wraps(model_attr)
def wrapper(*args, **kwargs):
"""Return only the T13 value from the model call."""
_, T13 = model_attr(*args, **kwargs)
return T13
return wrapper
return attr
# =============================================================================
# Container class to expose species-specific models
# =============================================================================
[docs]
class HeightTrajectoryContainer:
"""Expose height trajectory functions by species."""
picea_abies: HeightTrajectoryWrapper
pinus_sylvestris: HeightTrajectoryWrapper
def __init__(self) -> None:
"""Initialize species-specific height trajectory wrappers."""
self.picea_abies = HeightTrajectoryWrapper(HagglundSpruceModel)
self.pinus_sylvestris = HeightTrajectoryWrapper(HagglundPineModel)
[docs]
class TimeToBreastHeightContainer:
"""Expose T13 functions by species."""
picea_abies: TimeToBreastHeightWrapper
pinus_sylvestris: TimeToBreastHeightWrapper
def __init__(self) -> None:
"""Initialize species-specific T13 wrappers."""
self.picea_abies = TimeToBreastHeightWrapper(HagglundSpruceModel)
self.pinus_sylvestris = TimeToBreastHeightWrapper(HagglundPineModel)
[docs]
class Hagglund_1970:
"""
Provides a class-based interface to Hägglund's site index functions.
Callers can choose to obtain either the height trajectory (site index)
or the time to breast height (T13).
"""
height_trajectory: HeightTrajectoryContainer = HeightTrajectoryContainer()
time_to_breast_height: TimeToBreastHeightContainer = TimeToBreastHeightContainer()
regeneration = HagglundPineRegeneration
# ---------------------------------------------------------------------------
# Introspection
# ---------------------------------------------------------------------------
DESCRIPTOR = FormulaDescriptor(
component_id="hagglund_1970_siteindex",
source=SourceReference(
author="Hägglund, B.",
year=1972,
title="Om övre höjdens utveckling för gran i norra Sverige",
note=(
"Skogshögskolan, institutionen för skogsproduktion, Rapporter och uppsatser "
"nr 21. This module spans three of the author's publications and a single "
"reference cannot carry them all; each model is named for its own year. "
"Norway spruce, northern Sweden: Hägglund (1972), above. Norway spruce, "
"southern Sweden: Hägglund (1973) 'Om övre höjdens utveckling för gran i "
"södra Sverige', Rapporter och uppsatser nr 24. Scots pine: Hägglund (1974) "
"'Övre höjdens utveckling i tallbestånd', Rapporter och uppsatser nr 31. "
"There is no Hägglund (1970) site-index publication; the class name "
"Hagglund_1970 is a legacy identifier, not a citation."
),
),
species_groups={"pine": frozenset(), "spruce": frozenset()},
units={
"dominant_height_m": "m",
"age_years": "years",
"return": "SiteIndexValue (m at age 100)",
},
kernel_names=("Hagglund_1970", "HagglundSpruceModel", "HagglundPineModel"),
)