From ffdcc94897fa867a1e62dc59e4694aa7c68cb2b6 Mon Sep 17 00:00:00 2001 From: mati Date: Mon, 27 Jan 2025 13:16:06 -0300 Subject: [PATCH 1/3] Read CBD and CBH from raster layers in Kitral --- Cell2Fire/FuelModelKitral.cpp | 42 ++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/Cell2Fire/FuelModelKitral.cpp b/Cell2Fire/FuelModelKitral.cpp index de6a6e5d..16080e59 100644 --- a/Cell2Fire/FuelModelKitral.cpp +++ b/Cell2Fire/FuelModelKitral.cpp @@ -698,11 +698,26 @@ float crownfractionburn(inputs* data, main_outs* at, int FMC) { // generar output de cfb float a, cbd, ros, ros0, H, wa, i0, cbh, cfb; - cbh = cbhs[data->nftype][0]; + if (std::isnan(data->cbh)) + { + cbh = cbhs[data->nftype][0]; + } + else + { + cbh = data->cbh; + } i0 = pow((0.01 * cbh * (460 + 25.9 * FMC)), 1.5); H = hs[data->nftype][0]; wa = fls_david[data->nftype][0]; - cbd = cbds[data->nftype][0]; + if (std::isnan(data->cbd)) + { + cbd = cbds[data->nftype][0]; + } + else + { + cbd = data->cbd; + } + ros0 = 60 * i0 / (H * wa); ros = at->rss; if (cbd != 0) @@ -767,15 +782,26 @@ checkActive(inputs* data, main_outs* at, int FMC) // En KITRAL SE USA PL04 { float ros_critical, cbd, H, wa, i0, cbh; bool active; - cbh = cbhs[data->nftype][0]; + if (std::isnan(data->cbh)) + { + cbh = cbhs[data->nftype][0]; + } + else + { + cbh = data->cbh; + } i0 = pow((0.01 * cbh * (460 + 25.9 * FMC)), 1.5); H = hs[data->nftype][0]; wa = fls_david[data->nftype][0]; - cbd = cbds[data->nftype][0]; ros_critical = 60 * i0 / (H * wa); - - cbd = cbds[data->nftype][0]; - + if (std::isnan(data->cbd)) + { + cbd = cbds[data->nftype][0]; + } + else + { + cbd = data->cbd; + } active = cbd * ros_critical >= 3; return active; } @@ -870,7 +896,7 @@ calculate_k(inputs* data, // Step 10: Criterion for Crown Fire Initiation (no init if user does not // want to include it) - if (args->AllowCROS && cbhs[data->nftype][0] != 0) + if (args->AllowCROS && (data->cbh != 0 || cbhs[data->nftype][0] != 0)) { if (activeCrown) { // si el fuego esta activo en copas chequeamos condiciones From 882a18e6026a8ce05e4be2a682e425a1afefcbb8 Mon Sep 17 00:00:00 2001 From: mati Date: Tue, 28 Jan 2025 16:28:07 -0300 Subject: [PATCH 2/3] Calculate crown intensity and flamen length in kitral --- Cell2Fire/Cell2Fire.cpp | 12 ++--- Cell2Fire/Cells.cpp | 2 +- Cell2Fire/FuelModelKitral.cpp | 84 +++++++++++++++++++++++++++++++++-- Cell2Fire/FuelModelKitral.h | 6 +++ Cell2Fire/FuelModelSpain.h | 5 +++ 5 files changed, 99 insertions(+), 10 deletions(-) diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index 907907c9..12e7f9ce 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -682,7 +682,7 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) this->surfaceIntensityFolder = this->args.OutFolder + "SurfaceIntensity" + separator(); } // Crown Byram Intensity Folder - if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator == "S")) + if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator != "C")) { CSVWriter CSVFolder("", ""); this->crownIntensityFolder = this->args.OutFolder + "CrownIntensity"; @@ -698,7 +698,7 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) this->surfaceFlameLengthFolder = this->args.OutFolder + "SurfaceFlameLength" + separator(); } // Crown Flame Length Folder - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S")) + if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) { CSVWriter CSVFolder("", ""); this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength"; @@ -706,7 +706,7 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator(); } // max Flame Length Folder - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S")) + if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) { CSVWriter CSVFolder("", ""); this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength"; @@ -1963,7 +1963,7 @@ Cell2Fire::Results() } // Crown Intensity - if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator == "S")) + if ((this->args.OutIntensity) && (this->args.AllowCROS) && (this->args.Simulator != "C")) { this->crownIntensityFolder = this->args.OutFolder + "CrownIntensity" + separator(); std::string intensityName; @@ -1999,7 +1999,7 @@ Cell2Fire::Results() } // Crown Flame length - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S")) + if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) { this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator(); std::string fileName; @@ -2017,7 +2017,7 @@ Cell2Fire::Results() } // Max Flame length - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator == "S")) + if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) { this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength" + separator(); std::string fileName; diff --git a/Cell2Fire/Cells.cpp b/Cell2Fire/Cells.cpp index 696784a1..83d375f3 100644 --- a/Cell2Fire/Cells.cpp +++ b/Cell2Fire/Cells.cpp @@ -736,7 +736,7 @@ Cells::manageFire(int period, surfFraction[nb] = metrics.sfc; SurfaceFlameLengths[this->realId - 1] = mainstruct.fl; SurfaceFlameLengths[nb - 1] = metrics.fl; - if ((args->AllowCROS) && (args->Simulator == "S")) + if ((args->AllowCROS) && (args->Simulator != "C")) { float comp_zero = 0; MaxFlameLengths[this->realId - 1] diff --git a/Cell2Fire/FuelModelKitral.cpp b/Cell2Fire/FuelModelKitral.cpp index 16080e59..df624577 100644 --- a/Cell2Fire/FuelModelKitral.cpp +++ b/Cell2Fire/FuelModelKitral.cpp @@ -680,13 +680,78 @@ byram_intensity(inputs* data, main_outs* at) return ib; // unidad de medida } +/** + * Calculates byram fire intensity when there is active crown fire. + * In order for this to be calculated, the input folder must contain + * files with CBD, CBH and tree height data for each cell. + * @param at Structure containing the cell's output data. + * @param data Structure containing the cell's input data. + * @return Fire intensity. + */ + +float +crown_byram_intensity_k(main_outs* at, inputs* data) +{ + float cbd, cbh; + if (std::isnan(data->cbh)) + { + cbh = cbhs[data->nftype][0]; + } + else + { + cbh = data->cbh; + } + float canopy_height = cbh * 2; // TODO: use real tree height + if (canopy_height < 0) + { + throw std::runtime_error("Tree height is lower than canopy base height, " + "please provide valid files."); + } + if (std::isnan(data->cbd)) + { + cbd = cbds[data->nftype][0]; + } + else + { + cbd = data->cbd; + } + return std::ceil((hs[data->nftype][0] / 60) * cbd * canopy_height * at->ros_active * 100.0) / 100.0; +} + +/** + * @brief Calculates the flame length of a cell when there is crown fire. + * @param intensity Byram intensity for crown fires + * @return the flame length + */ + +float +crown_flame_length_k(float intensity) +{ + float fl = 0.1 * pow(intensity, 0.5); + if (fl < 0.01) + { + return 0; + } + else + { + return std::ceil(fl * 100.0) / 100.0; + } +} + bool fire_type(inputs* data, main_outs* at, int FMC) { float intensity, critical_intensity, cbh; bool crownFire = false; intensity = at->sfi; - cbh = cbhs[data->nftype][0]; + if (std::isnan(data->cbh)) + { + cbh = cbhs[data->nftype][0]; + } + else + { + cbh = data->cbh; + } // cbh = data->cbh; critical_intensity = pow((0.01 * cbh * (460 + 25.9 * FMC)), 1.5); if ((intensity > critical_intensity) && cbh != 0) @@ -854,7 +919,14 @@ calculate_k(inputs* data, FMC = args->FMC; ptr->nftype = data->nftype; ptr->fmc = fmcs[data->nftype][0]; - ptr->cbh = cbhs[data->nftype][0]; + if (isnan(data->cbh)) + { + ptr->cbh = cbhs[data->nftype][0]; + } + else + { + ptr->cbh = data->cbh; + } // cout << " cbh " << ptr->cbh << "\n"; ptr->fl = fls_david[data->nftype][0]; @@ -896,7 +968,7 @@ calculate_k(inputs* data, // Step 10: Criterion for Crown Fire Initiation (no init if user does not // want to include it) - if (args->AllowCROS && (data->cbh != 0 || cbhs[data->nftype][0] != 0)) + if (args->AllowCROS && (data->cbh > 0 || cbhs[data->nftype][0] != 0)) { if (activeCrown) { // si el fuego esta activo en copas chequeamos condiciones @@ -930,6 +1002,8 @@ calculate_k(inputs* data, at->rss = hptr->ros; bptr->ros = backfire_ros10_k(hptr, sec); fptr->ros = flankfire_ros_k(hptr->ros, bptr->ros, sec->lb); + at->crown_intensity = crown_byram_intensity_k(at, data); + at->crown_flame_length = crown_flame_length_k(at->crown_intensity); if (args->verbose) { @@ -952,6 +1026,8 @@ calculate_k(inputs* data, at->rss = hptr->ros; bptr->ros = backfire_ros10_k(hptr, sec); fptr->ros = flankfire_ros_k(hptr->ros, bptr->ros, sec->lb); + at->crown_intensity = crown_byram_intensity_k(at, data); + at->crown_flame_length = crown_flame_length_k(at->crown_intensity); if (args->verbose) { @@ -1038,6 +1114,8 @@ determine_destiny_metrics_k(inputs* data, fuel_coefs* ptr, arguments* args, main if (crownFire) { metrics->cfb = crownfractionburn(data, metrics, FMC); + metrics->crown_intensity = crown_byram_intensity_k(metrics, data); + metrics->crown_flame_length = crown_flame_length_k(metrics->crown_intensity); } if (args->verbose) { diff --git a/Cell2Fire/FuelModelKitral.h b/Cell2Fire/FuelModelKitral.h index c28cac5d..5e30c345 100644 --- a/Cell2Fire/FuelModelKitral.h +++ b/Cell2Fire/FuelModelKitral.h @@ -32,6 +32,9 @@ float backfire_ros_k(main_outs* at, snd_outs* sec); // Flame length [m]) float flame_length(inputs* data, main_outs* at); +// Crown Flame length [m]) +float crown_flame_length_k(float intensity); + // Angle of the flame w.r.t. horizontal surface (Putnam's) float angleFL(inputs* data, main_outs* at); @@ -41,6 +44,9 @@ float flame_height(inputs* data, main_outs* at); // byram intensity float byram_intensity(inputs* data, main_outs* at); +// crown byram intensity +float crown_byram_intensity_k(main_outs* at, inputs* data); + // Type of fire (Crown = 1, SROS = 0) bool fire_type(inputs* data, main_outs* atr, int FMC); // Active crown diff --git a/Cell2Fire/FuelModelSpain.h b/Cell2Fire/FuelModelSpain.h index 2aa0f8f6..98c9daef 100644 --- a/Cell2Fire/FuelModelSpain.h +++ b/Cell2Fire/FuelModelSpain.h @@ -35,6 +35,8 @@ float backfire_ros_s(main_outs* at, snd_outs* sec); // Flame length [m]) float flame_length(inputs* data, fuel_coefs* ptr); +// Crown Flame length [m]) +float crown_flame_length(float intensity); // Angle of the flame w.r.t. horizontal surface (Putnam's) float angleFL(inputs* data, fuel_coefs* ptr); @@ -45,6 +47,9 @@ float flame_height(inputs* data, fuel_coefs* ptr); // byram intensity float byram_intensity(inputs* data, fuel_coefs* ptr); +// crown byram intensity +float crown_byram_intensity(main_outs* at, inputs* data); + // Type of fire (Crown = 1, SROS = 0) bool fire_type(inputs* data, fuel_coefs* ptr); From a5251fcc897fd0f624f6ea9fa672cc8bf5b8b6af Mon Sep 17 00:00:00 2001 From: mati Date: Thu, 30 Jan 2025 16:17:45 -0300 Subject: [PATCH 3/3] [bug] Use real tree height --- Cell2Fire/FuelModelKitral.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Cell2Fire/FuelModelKitral.cpp b/Cell2Fire/FuelModelKitral.cpp index df624577..89f0cee2 100644 --- a/Cell2Fire/FuelModelKitral.cpp +++ b/Cell2Fire/FuelModelKitral.cpp @@ -701,7 +701,8 @@ crown_byram_intensity_k(main_outs* at, inputs* data) { cbh = data->cbh; } - float canopy_height = cbh * 2; // TODO: use real tree height + // CBH is 0,1012 * Height^(1,4822) + float canopy_height = std::pow(cbh / 0.1012, 1 / 1.4822) - cbh; if (canopy_height < 0) { throw std::runtime_error("Tree height is lower than canopy base height, "