{ "cells": [ { "cell_type": "markdown", "id": "7fb27b941602401d91542211134fc71a", "metadata": {}, "source": [ "# Young stand demo (NYSKOG + Näslund 1986)\n" ] }, { "cell_type": "markdown", "id": "acae54e37e7d407bbb7b55eff062a284", "metadata": {}, "source": [ "## Notebook Objectives\n", "- Build a young stand workflow from initialization through growth with reproducible inputs.\n", "- Provide runnable, copy-safe snippets that work in the docs build environment.\n", "\n", "## Prerequisites\n", "- Python environment with `pyforestry` installed from this repository.\n", "- Execute cells in order; random components should use fixed seeds where shown.\n", "\n", "## Sources\n", "- NYSKOG reconstruction and Swedish growth utilities in `pyforestry.sweden.adapters`.\n", "- Näslund (1986) height/diameter relationships where used in the workflow.\n" ] }, { "cell_type": "markdown", "id": "9a63283cbaf04dbcab1f6479b197f3a8", "metadata": {}, "source": [ "This notebook shows a small, reproducible workflow for constructing a Swedish site,\n", "reconstructing a young stand using the NYSKOG equations (Appendix 2), and running\n", "a short forward simulation with the young-stand equations and Näslund damage model.\n", "\n", "**Note:** The original HUGIN/NYSKOG database is not part of pyforestry, so this\n", "demo synthesizes an initial tree list from the reconstruction equations only.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "8dd0d8092fe74a7c96281538738b07e2", "metadata": { "execution": { "iopub.execute_input": "2026-02-19T23:10:09.552430Z", "iopub.status.busy": "2026-02-19T23:10:09.552259Z", "iopub.status.idle": "2026-02-19T23:10:10.755350Z", "shell.execute_reply": "2026-02-19T23:10:10.753899Z" } }, "outputs": [], "source": [ "import copy\n", "import math\n", "import random\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from pyforestry.base.helpers.primitives.sitebase import SiteBase\n", "from pyforestry.base.helpers.tree import Tree\n", "from pyforestry.base.helpers.tree_species import TreeSpecies\n", "from pyforestry.sweden.adapters.elfving_1982 import (\n", " HuginMeanHeightModel,\n", " NfiRegion,\n", " NyskogReconstruction,\n", " RegenerationType,\n", " nyskog_indicators_from_site,\n", ")\n", "from pyforestry.sweden.mortality.naslund_1986 import (\n", " Naslund1986DamageModel,\n", " SaplingSpeciesGroup,\n", ")\n", "from pyforestry.sweden.systems.nystrom_soderberg_1987 import NystromSoderberg1987\n", "from pyforestry.sweden.site import Sweden, SwedishSite\n", "from pyforestry.sweden.siteindex.sis.generated_site_category_trees import (\n", " predict_site_categories_county_tree,\n", ")\n", "from pyforestry.sweden.siteindex.sis.hagglund_lundmark_1977 import Hagglund_Lundmark_1977_SIS\n", "\n", "\n", "class SwedishSiteDemo(SwedishSite):\n", " \"\"\"Concrete wrapper for notebooks (implements SiteBase abstract method).\"\"\"\n", "\n", " def compute_attributes(self) -> None:\n", " SwedishSite.__post_init__(self)\n", "\n", " def __post_init__(self) -> None:\n", " SiteBase.__post_init__(self)\n", "\n", "\n", "random.seed(42)\n", "np.random.seed(42)\n" ] }, { "cell_type": "markdown", "id": "72eea5119410473aa328ad9291626812", "metadata": {}, "source": [ "## 1. Define a Swedish site\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "8edb47106e1a46a883d545849b8ab81b", "metadata": { "execution": { "iopub.execute_input": "2026-02-19T23:10:10.758883Z", "iopub.status.busy": "2026-02-19T23:10:10.758385Z", "iopub.status.idle": "2026-02-19T23:10:11.014865Z", "shell.execute_reply": "2026-02-19T23:10:11.013186Z" } }, "outputs": [ { "data": { "text/plain": [ "{'sis_closeness': {'requested_sis': 20.0,\n", " 'achieved_sis': 20.698654683125625,\n", " 'sis_abs_error': 0.6986546831256248,\n", " 'sis_rel_error_pct': 3.493273415628124},\n", " 'predicted_site_categories': {'field_layer': ,\n", " 'bottom_layer': ,\n", " 'soil_texture': ,\n", " 'soil_moisture': ,\n", " 'soil_depth': ,\n", " 'soil_water': ,\n", " 'ditched': False}}" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "SwedishSiteDemo(latitude=60.5, longitude=15.0, altitude=150.0, field_layer=, bottom_layer=, soil_texture=, soil_moisture=, soil_depth=, soil_water=, aspect=None, incline_percent=None, ditched=False, temperature_sum_odin1983=1215.1999999999998, county=, humidity=np.float64(75.0), distance_to_coast=np.float64(254.68051713105228), climate_zone=, sis_spruce_100=None, sis_pine_100=None, sis_birch_50=SiteIndexValue(19.889499999999984, reference_age=AgeMeasurement(50.0, code=2 [DBH]), species={TreeName(genus=TreeGenus(name='Betula', code='BETULA'), species_name='pendula', code='BPEN'), TreeName(genus=TreeGenus(name='Betula', code='BETULA'), species_name='pubescens', code='BPUB')}, fn=eriksson_1997_height_trajectory_sweden_birch), n_of_limes_norrlandicus=True)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "site_index_pine_m = 20.0\n", "site_index_spruce_m = 22.0\n", "county = Sweden.County.KOPPARBERG_OVRIGA\n", "site_category_species = \"Pinus sylvestris\"\n", "requested_sis = site_index_pine_m\n", "\n", "predicted_site_categories = predict_site_categories_county_tree(\n", " sis_hagglund_1979=requested_sis,\n", " species=site_category_species,\n", " Direktlan=county,\n", ")\n", "\n", "site = SwedishSiteDemo(\n", " latitude=60.5,\n", " longitude=15.0,\n", " altitude=150.0,\n", " field_layer=predicted_site_categories[\"field_layer\"],\n", " bottom_layer=predicted_site_categories[\"bottom_layer\"],\n", " soil_texture=predicted_site_categories[\"soil_texture\"],\n", " soil_moisture=predicted_site_categories[\"soil_moisture\"],\n", " soil_depth=predicted_site_categories[\"soil_depth\"],\n", " soil_water=predicted_site_categories[\"soil_water\"],\n", " ditched=predicted_site_categories[\"ditched\"],\n", ")\n", "\n", "achieved_sis = Hagglund_Lundmark_1977_SIS(\n", " species=site_category_species,\n", " latitude=site.latitude,\n", " altitude=site.altitude or 0.0,\n", " soil_moisture=site.soil_moisture,\n", " ground_layer=site.bottom_layer or Sweden.BottomLayer.FRESH_MOSS,\n", " vegetation=site.field_layer,\n", " soil_texture=site.soil_texture or Sweden.SoilTextureTill.SANDY,\n", " climate_code=site.climate_zone or Sweden.ClimateZone.K1,\n", " lateral_water=site.soil_water or Sweden.SoilWater.SELDOM_NEVER,\n", " soil_depth=site.soil_depth or Sweden.SoilDepth.DEEP,\n", " incline_percent=site.incline_percent or 0.0,\n", " aspect=site.aspect or 0.0,\n", " nfi_adjustments=True,\n", " dlan=site.county or county,\n", " ditched=bool(site.ditched),\n", " peat=False,\n", " gotland=False,\n", " coast=(site.distance_to_coast or 9999.0) < 50.0,\n", " limes_norrlandicus=bool(site.n_of_limes_norrlandicus),\n", ")\n", "\n", "sis_closeness = {\n", " \"requested_sis\": requested_sis,\n", " \"achieved_sis\": float(achieved_sis),\n", " \"sis_abs_error\": abs(float(achieved_sis) - requested_sis),\n", " \"sis_rel_error_pct\": 100.0 * abs(float(achieved_sis) - requested_sis) / requested_sis,\n", "}\n", "\n", "display({\"sis_closeness\": sis_closeness, \"predicted_site_categories\": predicted_site_categories})\n", "site" ] }, { "cell_type": "markdown", "id": "10185d26023b46108eb7d9f57d49d2b3", "metadata": {}, "source": [ "## 2. Reconstruct a young stand using NYSKOG equations\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8763a12b2bbd4a93a75aff182afb95dc", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'pine': 197.21557785033062,\n", " 'spruce': 435.7513070153287,\n", " 'contorta': 0.0,\n", " 'larch': 0.0,\n", " 'birch': 1159.603329004732,\n", " 'other_broadleaf': 48.31680537519716}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "regen_type = RegenerationType.NATURAL\n", "species_to_plant = TreeSpecies.Sweden.pinus_sylvestris\n", "nfi_region = NfiRegion.REG3\n", "\n", "age_years = 12.0\n", "\n", "mean_height_main_m = HuginMeanHeightModel.mean_height(\n", " age_years=age_years,\n", " species=species_to_plant,\n", " site_index_pine_m=site_index_pine_m,\n", " site_index_spruce_m=site_index_spruce_m,\n", ")\n", "\n", "indicators = nyskog_indicators_from_site(\n", " field_layer=site.field_layer,\n", " soil_moisture=site.soil_moisture,\n", ")\n", "\n", "# NOTE: ASINW is typically derived from the published Elfving regeneration functions (Appendix 2).\n", "asinw = 110.0\n", "q = NyskogReconstruction.production_potential_q(asinw)\n", "ln_q = math.log(q)\n", "ln_si = math.log(site_index_pine_m)\n", "\n", "under_dimension_prob = NyskogReconstruction.udim_probability(q)\n", "height_indicator_dm = max(15.0, 10.0 * mean_height_main_m)\n", "\n", "stem_total = NyskogReconstruction.total_stems(\n", " regeneration_type=regen_type,\n", " mean_height_main_m=mean_height_main_m,\n", " q=q,\n", " ln_q=ln_q,\n", " ln_si=ln_si,\n", " under_dimension_prob=under_dimension_prob,\n", " wet=indicators[\"wet\"],\n", " dry=indicators[\"dry\"],\n", " height_indicator_dm=height_indicator_dm,\n", " deterministic=True,\n", ")\n", "\n", "prop_conifer = NyskogReconstruction.proportion_conifer(\n", " regeneration_type=regen_type,\n", " q=q,\n", " ln_qind=ln_q,\n", " stem_total=stem_total,\n", " ln_si=ln_si,\n", " wet=indicators[\"wet\"],\n", " dry=indicators[\"dry\"],\n", " rich=indicators[\"rich\"],\n", " poor=indicators[\"poor\"],\n", " deterministic=True,\n", ")\n", "\n", "prop_dom_conifer = NyskogReconstruction.dominant_conifer_share(\n", " regeneration_type=regen_type,\n", " qind=q,\n", " ln_si=ln_si,\n", " wet=indicators[\"wet\"],\n", " dry=indicators[\"dry\"],\n", " rich=indicators[\"rich\"],\n", " poor=indicators[\"poor\"],\n", " hwod=indicators[\"hwod\"],\n", " hwd=indicators[\"hwd\"],\n", " shrubs=indicators[\"shrubs\"],\n", " lichen=indicators[\"lichen\"],\n", " deterministic=True,\n", ")\n", "\n", "stems = NyskogReconstruction.stems_per_species(\n", " regeneration_type=regen_type,\n", " species_to_plant=species_to_plant,\n", " stem_total=stem_total,\n", " prop_conifer=prop_conifer,\n", " prop_dom_conifer=prop_dom_conifer,\n", " site_index_m=site_index_pine_m,\n", " nfi_region=nfi_region,\n", ")\n", "\n", "stems\n" ] }, { "cell_type": "markdown", "id": "7623eae2785240b9bd12b16a66d81610", "metadata": {}, "source": [ "## 3. Sample a tree list from reconstructed stand metrics\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7cdc8c89c7104fffa095e18ddfef8986", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "199" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def mean_height_secondary(species: TreeSpecies) -> float:\n", " site_index_m = site_index_spruce_m if species in {\n", " TreeSpecies.Sweden.picea_abies,\n", " TreeSpecies.Sweden.picea_sitchensis,\n", " TreeSpecies.Sweden.picea_mariana,\n", " } else site_index_pine_m\n", "\n", " return NyskogReconstruction.secondary_mean_height(\n", " regeneration_type=regen_type,\n", " secondary_species=species,\n", " site_index_m=site_index_m,\n", " mean_height_main_m=mean_height_main_m,\n", " herb=indicators[\"herb\"],\n", " dry=indicators[\"dry\"],\n", " wet=indicators[\"wet\"],\n", " deterministic=True,\n", " )\n", "\n", "\n", "species_map = {\n", " \"pine\": TreeSpecies.Sweden.pinus_sylvestris,\n", " \"spruce\": TreeSpecies.Sweden.picea_abies,\n", " \"contorta\": TreeSpecies.Sweden.pinus_contorta,\n", " \"larch\": TreeSpecies.Sweden.larix_sibirica,\n", " \"birch\": TreeSpecies.Sweden.betula_pendula,\n", " \"other_broadleaf\": TreeSpecies.Sweden.populus_tremula,\n", "}\n", "\n", "mean_heights = {\n", " \"pine\": mean_height_main_m,\n", " \"spruce\": mean_height_secondary(TreeSpecies.Sweden.picea_abies),\n", " \"contorta\": mean_height_secondary(TreeSpecies.Sweden.pinus_contorta),\n", " \"larch\": mean_height_secondary(TreeSpecies.Sweden.larix_sibirica),\n", " \"birch\": mean_height_secondary(TreeSpecies.Sweden.betula_pendula),\n", " \"other_broadleaf\": mean_height_secondary(TreeSpecies.Sweden.populus_tremula),\n", "}\n", "\n", "sample_n = 200\n", "total_stems = sum(stems.values())\n", "trees: list[Tree] = []\n", "\n", "for label, stems_per_ha in stems.items():\n", " if stems_per_ha <= 0.0:\n", " continue\n", " share = stems_per_ha / total_stems\n", " n_trees = max(1, int(round(sample_n * share)))\n", "\n", " species = species_map[label]\n", " mean_height = mean_heights[label]\n", " cvh = NyskogReconstruction.height_variation(\n", " species=species,\n", " species_height_m=mean_height,\n", " q=q,\n", " ln_q=ln_q,\n", " self_rejuvenated=1,\n", " deterministic=True,\n", " )\n", " beta, shape = NyskogReconstruction.weibull_parameters(\n", " species=species,\n", " cvh=cvh,\n", " mean_height_m=mean_height,\n", " )\n", "\n", " heights = beta * np.random.weibull(shape, size=n_trees)\n", " weight = stems_per_ha / n_trees\n", "\n", " for height in heights:\n", " trees.append(\n", " Tree(\n", " species=species,\n", " height_m=float(height),\n", " weight_n=weight,\n", " )\n", " )\n", "\n", "base_trees = copy.deepcopy(trees)\n", "\n", "len(trees)\n" ] }, { "cell_type": "markdown", "id": "b118ea5561624da68c537baed56e602f", "metadata": {}, "source": [ "## 4. DBH and Näslund damage proportions\n" ] }, { "cell_type": "code", "execution_count": null, "id": "938c804e27f84196a10c8828c723f798", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{: 0.7583307366236001,\n", " : 0.2994340399604417,\n", " : 0.46655281692975714,\n", " : 0.3419758275552662,\n", " : 0.7987225971416196}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def stand_metrics(tree_list: list[Tree]) -> dict[str, float]:\n", " heights = np.array([t.height_m or 0.0 for t in tree_list])\n", " weights = np.array([t.weight_n or 0.0 for t in tree_list])\n", " total_weight = weights.sum()\n", " if total_weight <= 0.0:\n", " return {\n", " \"mean_height_m\": 0.0,\n", " \"total_height_sqr_m2_per_100m2\": 0.0,\n", " \"broadleaf_height_sqr_share\": 0.0,\n", " \"h_max_m\": 0.0,\n", " }\n", " mean_height = float((heights * weights).sum() / total_weight)\n", "\n", " height_sqr_sum = float(((heights**2) * weights).sum())\n", " total_height_sqr_m2_per_100m2 = height_sqr_sum * 0.01\n", "\n", " broadleaf_species = {\n", " TreeSpecies.Sweden.betula_pendula,\n", " TreeSpecies.Sweden.betula_pubescens,\n", " TreeSpecies.Sweden.populus_tremula,\n", " }\n", " broadleaf_mask = np.array([t.species in broadleaf_species for t in tree_list])\n", " broadleaf_sqr_sum = float(((heights**2) * weights * broadleaf_mask).sum())\n", " if height_sqr_sum > 0.0:\n", " broadleaf_share = broadleaf_sqr_sum / height_sqr_sum\n", " else:\n", " broadleaf_share = 0.0\n", "\n", " top_heights = sorted(heights, reverse=True)[:3]\n", " h_max_m = float(np.mean(top_heights)) if top_heights else 0.0\n", "\n", " return {\n", " \"mean_height_m\": mean_height,\n", " \"total_height_sqr_m2_per_100m2\": total_height_sqr_m2_per_100m2,\n", " \"broadleaf_height_sqr_share\": broadleaf_share,\n", " \"h_max_m\": h_max_m,\n", " }\n", "\n", "\n", "def sapling_group(species: TreeSpecies) -> SaplingSpeciesGroup:\n", " if species in {\n", " TreeSpecies.Sweden.pinus_sylvestris,\n", " TreeSpecies.Sweden.larix_sibirica,\n", " TreeSpecies.Sweden.larix_decidua,\n", " TreeSpecies.Sweden.larix_europaea_x_leptolepis,\n", " TreeSpecies.Sweden.larix_sukaczewii,\n", " }:\n", " return SaplingSpeciesGroup.PINE\n", " if species == TreeSpecies.Sweden.pinus_contorta:\n", " return SaplingSpeciesGroup.CONTORTA\n", " if species == TreeSpecies.Sweden.picea_abies:\n", " return SaplingSpeciesGroup.SPRUCE\n", " if species in {\n", " TreeSpecies.Sweden.betula_pendula,\n", " TreeSpecies.Sweden.betula_pubescens,\n", " }:\n", " return SaplingSpeciesGroup.BIRCH\n", " if species in {\n", " TreeSpecies.Sweden.populus_tremula,\n", " TreeSpecies.Sweden.populus_tremula_x_tremuloides,\n", " }:\n", " return SaplingSpeciesGroup.ASPEN\n", " return SaplingSpeciesGroup.OTHER_BROADLEAF\n", "\n", "\n", "def apply_dbh(tree_list: list[Tree]) -> None:\n", " metrics = stand_metrics(tree_list)\n", " if metrics[\"mean_height_m\"] <= 0.0:\n", " return\n", " for tree in tree_list:\n", " if tree.height_m is None or tree.species is None:\n", " continue\n", " tree.diameter_cm = NystromSoderberg1987.dbh_from_height(\n", " height_dm=tree.height_m * 10.0,\n", " species=tree.species,\n", " total_height_sqr_m2_per_100m2=metrics[\"total_height_sqr_m2_per_100m2\"],\n", " broadleaf_height_sqr_share=metrics[\"broadleaf_height_sqr_share\"],\n", " natural_regeneration=1,\n", " cleaning_indicator=0,\n", " years_since_cleaning=0,\n", " altitude_m=site.altitude or 0.0,\n", " latitude_deg=site.latitude,\n", " shrubs=indicators[\"shrubs\"],\n", " herb_grass=indicators[\"herb\"],\n", " near_coast=1 if (site.distance_to_coast or 99.0) < 5.0 else 0,\n", " site_index_pine_m=site_index_pine_m,\n", " h_max_m=metrics[\"h_max_m\"],\n", " veg=0,\n", " )\n", "\n", "\n", "def damage_proportions(tree_list: list[Tree]) -> dict[SaplingSpeciesGroup, float]:\n", " stems = {sg: 0.0 for sg in SaplingSpeciesGroup}\n", " height_sums = {sg: 0.0 for sg in SaplingSpeciesGroup}\n", "\n", " for tree in tree_list:\n", " if tree.species is None or tree.height_m is None:\n", " continue\n", " group = sapling_group(tree.species)\n", " stems[group] += tree.weight_n or 0.0\n", " height_sums[group] += (tree.weight_n or 0.0) * tree.height_m\n", "\n", " mean_heights = {\n", " sg: (height_sums[sg] / stems[sg] if stems[sg] > 0.0 else 0.0)\n", " for sg in SaplingSpeciesGroup\n", " }\n", "\n", " return Naslund1986DamageModel.damage_proportions(\n", " stems=stems,\n", " mean_heights=mean_heights,\n", " site_index_pine_m=site_index_pine_m,\n", " site_index_spruce_m=site_index_spruce_m,\n", " latitude_deg=site.latitude,\n", " altitude_m=site.altitude or 0.0,\n", " )\n", "\n", "\n", "def apply_damage_mortality(\n", " tree_list: list[Tree],\n", " damage_props: dict[SaplingSpeciesGroup, float],\n", ") -> None:\n", " for tree in tree_list:\n", " if tree.species is None:\n", " continue\n", " group = sapling_group(tree.species)\n", " damage_prop = damage_props.get(group, 0.0)\n", " damage_prop = min(1.0, max(0.0, damage_prop))\n", " if tree.weight_n is not None:\n", " tree.weight_n *= 1.0 - damage_prop\n", "\n", "\n", "apply_dbh(trees)\n", "damage_proportions(trees)\n" ] }, { "cell_type": "markdown", "id": "504fb2a444614c0babb325280ed9130a", "metadata": {}, "source": [ "## 5. Simple forward simulation\n" ] }, { "cell_type": "code", "execution_count": null, "id": "59bbdb311c014d738909a11f9e486628", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{'age_years': 12.0,\n", " 'mean_height_m': 1.6807176071804812,\n", " 'mean_dbh_cm': 1.2446607824795128,\n", " 'pine_damage_prop': 0.7583307366236001,\n", " 'spruce_damage_prop': 0.2994340399604417,\n", " 'birch_damage_prop': 0.3419758275552662},\n", " {'age_years': 17.0,\n", " 'mean_height_m': 2.8046456524991306,\n", " 'mean_dbh_cm': 2.919427409360174,\n", " 'pine_damage_prop': 0.7417234871043763,\n", " 'spruce_damage_prop': 0.2990151912828163,\n", " 'birch_damage_prop': 0.2680638355674637},\n", " {'age_years': 22.0,\n", " 'mean_height_m': 4.391840787559593,\n", " 'mean_dbh_cm': 5.2571613473828025,\n", " 'pine_damage_prop': 0.7256788035489241,\n", " 'spruce_damage_prop': 0.300457614773755,\n", " 'birch_damage_prop': 0.19750103092563362},\n", " {'age_years': 27.0,\n", " 'mean_height_m': 6.089057804888556,\n", " 'mean_dbh_cm': 7.574174717815942,\n", " 'pine_damage_prop': 0.7178690972371207,\n", " 'spruce_damage_prop': 0.30467528492400175,\n", " 'birch_damage_prop': 0.1546302859552517}]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def stand_summary(tree_list: list[Tree], age_years: float) -> dict[str, float]:\n", " heights = np.array([t.height_m or 0.0 for t in tree_list])\n", " diameters = np.array([t.diameter_cm or 0.0 for t in tree_list])\n", " weights = np.array([t.weight_n or 0.0 for t in tree_list])\n", " total_weight = weights.sum()\n", " if total_weight <= 0.0:\n", " return {\n", " \"age_years\": age_years,\n", " \"mean_height_m\": 0.0,\n", " \"mean_dbh_cm\": 0.0,\n", " \"basal_area_m2_ha\": 0.0,\n", " \"stems_per_ha\": 0.0,\n", " \"pine_damage_prop\": 0.0,\n", " \"spruce_damage_prop\": 0.0,\n", " \"birch_damage_prop\": 0.0,\n", " }\n", " mean_height = float((heights * weights).sum() / total_weight)\n", " mean_dbh = float((diameters * weights).sum() / total_weight)\n", " basal_area = float((math.pi * (diameters / 200.0) ** 2 * weights).sum())\n", "\n", " damages = damage_proportions(tree_list)\n", "\n", " return {\n", " \"age_years\": age_years,\n", " \"mean_height_m\": mean_height,\n", " \"mean_dbh_cm\": mean_dbh,\n", " \"basal_area_m2_ha\": basal_area,\n", " \"stems_per_ha\": float(total_weight),\n", " \"pine_damage_prop\": damages.get(SaplingSpeciesGroup.PINE, 0.0),\n", " \"spruce_damage_prop\": damages.get(SaplingSpeciesGroup.SPRUCE, 0.0),\n", " \"birch_damage_prop\": damages.get(SaplingSpeciesGroup.BIRCH, 0.0),\n", " }\n", "\n", "\n", "def advance_years(tree_list: list[Tree], age_years: float, dt: float) -> float:\n", " current_mean_height = stand_metrics(tree_list)[\"mean_height_m\"]\n", " new_mean_height = HuginMeanHeightModel.mean_height(\n", " age_years=age_years + dt,\n", " species=species_to_plant,\n", " site_index_pine_m=site_index_pine_m,\n", " site_index_spruce_m=site_index_spruce_m,\n", " )\n", " if current_mean_height > 0.0:\n", " scale = new_mean_height / current_mean_height\n", " else:\n", " scale = 1.0\n", " for tree in tree_list:\n", " if tree.height_m is not None:\n", " tree.height_m *= scale\n", " apply_dbh(tree_list)\n", " return age_years + dt\n", "\n", "\n", "def run_simulation(\n", " base_tree_list: list[Tree],\n", " age_start: float,\n", " steps: int,\n", " dt: float,\n", " *,\n", " apply_damage: bool,\n", ") -> list[dict[str, float]]:\n", " trees = copy.deepcopy(base_tree_list)\n", " apply_dbh(trees)\n", " age = age_start\n", " results: list[dict[str, float]] = []\n", " for step in range(steps + 1):\n", " results.append(stand_summary(trees, age))\n", " if step == steps:\n", " break\n", " age = advance_years(trees, age, dt=dt)\n", " if apply_damage:\n", " damage_props = damage_proportions(trees)\n", " apply_damage_mortality(trees, damage_props)\n", " return results\n", "\n", "\n", "baseline_results = run_simulation(\n", " base_tree_list=base_trees,\n", " age_start=age_years,\n", " steps=4,\n", " dt=5.0,\n", " apply_damage=False,\n", ")\n", "damage_results = run_simulation(\n", " base_tree_list=base_trees,\n", " age_start=age_years,\n", " steps=4,\n", " dt=5.0,\n", " apply_damage=True,\n", ")\n", "\n", "baseline_results, damage_results\n" ] }, { "cell_type": "markdown", "id": "b43b363d81ae4b689946ece5c682cd59", "metadata": {}, "source": [ "## 6. Damage vs no-damage comparison\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8a65eabff63a45729fe45fb5ade58bdc", "metadata": {}, "outputs": [], "source": [ "def plot_comparison(\n", " baseline: list[dict[str, float]],\n", " damaged: list[dict[str, float]],\n", ") -> None:\n", " ages = [row[\"age_years\"] for row in baseline]\n", " fig, axes = plt.subplots(2, 2, figsize=(10, 8), sharex=True)\n", " axes = axes.flatten()\n", "\n", " metrics = [\n", " (\"mean_height_m\", \"Mean height (m)\"),\n", " (\"mean_dbh_cm\", \"Mean DBH (cm)\"),\n", " (\"basal_area_m2_ha\", \"Basal area (m²/ha)\"),\n", " (\"stems_per_ha\", \"Stems (n/ha)\"),\n", " ]\n", "\n", " for ax, (key, label) in zip(axes, metrics):\n", " ax.plot(ages, [row[key] for row in baseline], label=\"No damage\")\n", " ax.plot(ages, [row[key] for row in damaged], label=\"Näslund damage\")\n", " ax.set_ylabel(label)\n", " ax.grid(True, alpha=0.3)\n", "\n", " axes[-2].set_xlabel(\"Age (years)\")\n", " axes[-1].set_xlabel(\"Age (years)\")\n", " axes[0].legend(loc=\"best\")\n", " fig.suptitle(\"Young stand development with and without damage\", fontsize=12)\n", " fig.tight_layout()\n", "\n", "\n", "plot_comparison(baseline_results, damage_results)\n" ] } ], "metadata": { "kernelspec": { "display_name": "pyforestry", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.2" } }, "nbformat": 4, "nbformat_minor": 5 }