From ffdcc94897fa867a1e62dc59e4694aa7c68cb2b6 Mon Sep 17 00:00:00 2001 From: mati Date: Mon, 27 Jan 2025 13:16:06 -0300 Subject: [PATCH 01/10] 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 02/10] 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 03/10] [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, " From 3af4895621d1568e2947afcc89073b1e397da46b Mon Sep 17 00:00:00 2001 From: "matias.vilches@usach.cl" Date: Wed, 5 Feb 2025 17:17:05 -0300 Subject: [PATCH 04/10] new-out-fl-format --- Cell2Fire/Cell2Fire.cpp | 101 +++++++++++++--------------------------- Cell2Fire/Cell2Fire.h | 1 + Cell2Fire/WriteCSV.cpp | 33 +++++++++++++ Cell2Fire/WriteCSV.h | 1 + 4 files changed, 67 insertions(+), 69 deletions(-) diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index 12e7f9ce..45dc0964 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -47,6 +47,10 @@ std::unordered_map> BBOFactors; std::unordered_map> HarvestedCells; std::vector NFTypesCells; std::unordered_map IgnitionHistory; +std::unordered_map> Statistics; +std::unordered_map meanSurfaceFlameLength; +std::unordered_map meanCrownFlameLength; +std::unordered_map meanMaxFlameLength; /****************************************************************************** Utils @@ -689,30 +693,7 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) CSVFolder.MakeDir(this->crownIntensityFolder); this->crownIntensityFolder = this->args.OutFolder + "CrownIntensity" + separator(); } - // Surface Flame Length Folder - if (this->args.OutFl) - { - CSVWriter CSVFolder("", ""); - this->surfaceFlameLengthFolder = this->args.OutFolder + "SurfaceFlameLength"; - CSVFolder.MakeDir(this->surfaceFlameLengthFolder); - this->surfaceFlameLengthFolder = this->args.OutFolder + "SurfaceFlameLength" + separator(); - } - // Crown Flame Length Folder - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) - { - CSVWriter CSVFolder("", ""); - this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength"; - CSVFolder.MakeDir(this->crownFlameLengthFolder); - this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator(); - } - // max Flame Length Folder - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) - { - CSVWriter CSVFolder("", ""); - this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength"; - CSVFolder.MakeDir(this->maxFlameLengthFolder); - this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength" + separator(); - } + // Crown Folder if (this->args.OutCrown && this->args.AllowCROS) { @@ -1980,58 +1961,40 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->crownIntensities); } - // Surface Flame Length + // Flame Length Statistics + // This will be used to store accumulated flame lengths across simulations if (this->args.OutFl) { - this->surfaceFlameLengthFolder = this->args.OutFolder + "SurfaceFlameLength" + separator(); - std::string flName; - std::ostringstream oss; - oss.str(""); - oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; - flName = this->surfaceFlameLengthFolder + "SurfaceFlameLength" + oss.str() + ".asc"; - if (this->args.verbose) - { - std::cout << "We are generating the Flame Length to a asc file " << flName << std::endl; - } - CSVWriter CSVPloter(flName, " "); - CSVPloter.printASCII( - this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->surfaceFlameLengths); - } - // Crown Flame length - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) - { - this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator(); - std::string fileName; - std::ostringstream oss; - oss.str(""); - oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; - fileName = this->crownFlameLengthFolder + "CrownFlameLength" + oss.str() + ".asc"; - if (this->args.verbose) + for (int cell : this->burntCells) { - std::cout << "We are generating the Crown Flame Length to a asc file " << fileName << std::endl; - } - CSVWriter CSVPloter(fileName, " "); - CSVPloter.printASCII( - this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->crownFlameLengths); - } + std::vector flameLengthMeans; - // Max Flame length - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) - { - this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength" + separator(); - std::string fileName; - std::ostringstream oss; - oss.str(""); - oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; - fileName = this->maxFlameLengthFolder + "MaxFlameLength" + oss.str() + ".asc"; - if (this->args.verbose) + float surfaceFlameLength = this->surfaceFlameLengths[cell] / args.TotalSims; + meanSurfaceFlameLength[cell] += surfaceFlameLength; + flameLengthMeans.push_back(meanSurfaceFlameLength[cell]); + + if ((this->args.AllowCROS) && (this->args.Simulator != "C")) + { + float crownFlameLength = this->crownFlameLengths[cell] / args.TotalSims; + float maxFlameLength = this->maxFlameLengths[cell] / args.TotalSims; + meanCrownFlameLength[cell] += crownFlameLength; + meanMaxFlameLength[cell] += maxFlameLength; + flameLengthMeans.push_back(meanCrownFlameLength[cell]); + flameLengthMeans.push_back(meanMaxFlameLength[cell]); + } + + Statistics[cell] = { flameLengthMeans }; + } + if (currentSim == args.TotalSims) { - std::cout << "We are generating the Max Flame Length to a asc file " << fileName << std::endl; + std::string filename = "statistics.csv"; + CSVWriter statisticsFolder("", ""); + this->statsFolder = this->args.OutFolder + "Statistics" + separator(); + statisticsFolder.MakeDir(this->statsFolder); + CSVWriter statsFile(this->statsFolder + filename); + statsFile.printStats(Statistics); } - CSVWriter CSVPloter(fileName, " "); - CSVPloter.printASCII( - this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->maxFlameLengths); } // Intensity diff --git a/Cell2Fire/Cell2Fire.h b/Cell2Fire/Cell2Fire.h index 660e009d..83e43442 100644 --- a/Cell2Fire/Cell2Fire.h +++ b/Cell2Fire/Cell2Fire.h @@ -85,6 +85,7 @@ class Cell2Fire string crownFlameLengthFolder; string historyFolder; string ignitionsFolder; + string statsFolder; // Vectors std::vector fire_period; diff --git a/Cell2Fire/WriteCSV.cpp b/Cell2Fire/WriteCSV.cpp index 4535c7d8..59ec6ac4 100644 --- a/Cell2Fire/WriteCSV.cpp +++ b/Cell2Fire/WriteCSV.cpp @@ -287,6 +287,39 @@ CSVWriter::printIgnitions(std::unordered_map ignitionsHistory) ofs.close(); } +void +CSVWriter::printStats(std::unordered_map> statistics) +{ + std::ofstream ofs(fileName, std::ofstream::out); + + // Print column titles + ofs << "cell,surfaceFlameLength,crownFlameLength,maxFlameLength\n"; + + // Iterate through the unordered_map using iterator + for (auto it = statistics.begin(); it != statistics.end(); ++it) + { + int sim = it->first; // The simulation ID (key) + const std::vector& flameLengths = it->second; // The associated flame lengths (value) + + ofs << sim; // Print the simulation number (key) + + // Print all the flame lengths associated with the current simulation + for (size_t j = 0; j < flameLengths.size(); j++) + { + ofs << "," << flameLengths[j]; + if (j < flameLengths.size() - 1) + { + ofs << " "; // Optional: space between multiple flame lengths + } + } + + ofs << "\n"; // Move to the next line for each entry + } + + // Close the file + ofs.close(); +} + void CSVWriter::printCSV_V2(int rows, int cols, std::vector statusCells) { diff --git a/Cell2Fire/WriteCSV.h b/Cell2Fire/WriteCSV.h index ea20f493..029a8162 100644 --- a/Cell2Fire/WriteCSV.h +++ b/Cell2Fire/WriteCSV.h @@ -44,6 +44,7 @@ class CSVWriter void asciiHeader(int rows, int cols, double xllcorner, double yllcorner, int cellside); void printWeather(std::vector weatherHistory); void printIgnitions(std::unordered_map ignitionsHistory); + void printStats(std::unordered_map> statistics); void printCSV_V2(int rows, int cols, std::vector statusCells); // Function to create a directory void MakeDir(std::string pathPlot); From f26d3d7f9cdd3e2f7a1adf4e2036ca740ffbbe99 Mon Sep 17 00:00:00 2001 From: mati Date: Mon, 27 Jan 2025 13:16:06 -0300 Subject: [PATCH 05/10] Read CBD and CBH from raster layers in Kitral --- Cell2Fire/FuelModelKitral.cpp | 93 ++++++++++++++++------------------- 1 file changed, 43 insertions(+), 50 deletions(-) diff --git a/Cell2Fire/FuelModelKitral.cpp b/Cell2Fire/FuelModelKitral.cpp index 582608a0..16080e59 100644 --- a/Cell2Fire/FuelModelKitral.cpp +++ b/Cell2Fire/FuelModelKitral.cpp @@ -559,23 +559,6 @@ setup_const() cbds.insert(std::make_pair(DX02, cbd_dx02)); } -/** - * @brief Sets default value for crown constants when no raster files are provided for them. - * @param data Cell data - */ -void -setup_crown_const(inputs* data) -{ - if (data->cbd == -9999) - { - data->cbd = cbds[data->nftype][0]; - } - if (data->cbh == -9999) - { - data->cbh = cbhs[data->nftype][0]; - } -} - float rate_of_spread_k(inputs* data, fuel_coefs* ptr, @@ -612,7 +595,6 @@ rate_of_spread_k(inputs* data, return at->rss * (at->rss >= 0); } -// TODO: citation needed float flankfire_ros_k(float ros, float bros, float lb) { @@ -620,7 +602,6 @@ flankfire_ros_k(float ros, float bros, float lb) } /* ----------------- Length-to-Breadth --------------------------*/ -// TODO: citation needed float l_to_b(float ws, fuel_coefs* ptr) { @@ -632,7 +613,6 @@ l_to_b(float ws, fuel_coefs* ptr) } /* ----------------- Back Rate of Spread --------------------------*/ -// TODO: citation needed float backfire_ros_k(main_outs* at, snd_outs* sec) { @@ -646,7 +626,6 @@ backfire_ros_k(main_outs* at, snd_outs* sec) return bros * (bros >= 0); } -// TODO: citation needed float slope_effect(float elev_i, float elev_j, int cellsize) { @@ -657,7 +636,6 @@ slope_effect(float elev_i, float elev_j, int cellsize) return se; } -// TODO: citation needed float flame_length(inputs* data, main_outs* at) // REVISAR ESTA ECUACI�N { @@ -669,7 +647,6 @@ flame_length(inputs* data, main_outs* at) // REVISAR ESTA ECUACI�N return fl; } -// TODO: citation needed float angleFL(inputs* data, main_outs* at) { @@ -682,7 +659,6 @@ angleFL(inputs* data, main_outs* at) return angle; } -// TODO: citation needed float flame_height(inputs* data, main_outs* at) { @@ -692,7 +668,6 @@ flame_height(inputs* data, main_outs* at) return fh; } -// TODO: citation needed float byram_intensity(inputs* data, main_outs* at) { @@ -711,28 +686,41 @@ fire_type(inputs* data, main_outs* at, int FMC) float intensity, critical_intensity, cbh; bool crownFire = false; intensity = at->sfi; - cbh = data->cbh; + cbh = cbhs[data->nftype][0]; + // cbh = data->cbh; critical_intensity = pow((0.01 * cbh * (460 + 25.9 * FMC)), 1.5); - if ((intensity > critical_intensity) && cbh > 0) + if ((intensity > critical_intensity) && cbh != 0) crownFire = true; return crownFire; } -// TODO: citation needed float crownfractionburn(inputs* data, main_outs* at, int FMC) -{ - // generar output de cfb +{ // generar output de cfb float a, cbd, ros, ros0, H, wa, i0, cbh, cfb; - cbh = data->cbh; + 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 = data->cbd; + 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) + if (cbd != 0) { a = -log(0.1) / (0.9 * (3.0 / cbd - ros0)); } @@ -745,7 +733,6 @@ crownfractionburn(inputs* data, main_outs* at, int FMC) return cfb; } -// TODO: citation needed float active_rate_of_spreadPL04(inputs* data, main_outs* at) // En KITRAL SE USA PL04 @@ -776,11 +763,10 @@ active_rate_of_spreadPL04(inputs* data, rospl04 = fmc * fch * (fv + ps); } ros_active = 3.34 * rospl04; // if rac*cbd>3.0, aplicar - // ros_final=3.34*rospl04 + // ros_final=3.34*rospl04 return ros_active; } -// TODO: citation needed float final_rate_of_spreadPL04(main_outs* at) // En KITRAL SE USA PL04 { @@ -796,19 +782,30 @@ checkActive(inputs* data, main_outs* at, int FMC) // En KITRAL SE USA PL04 { float ros_critical, cbd, H, wa, i0, cbh; bool active; - cbh = data->cbh; - + 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]; ros_critical = 60 * i0 / (H * wa); - cbd = data->cbd; - + if (std::isnan(data->cbd)) + { + cbd = cbds[data->nftype][0]; + } + else + { + cbd = data->cbd; + } active = cbd * ros_critical >= 3; return active; } -// TODO: citation needed float backfire_ros10_k(fire_struc* hptr, snd_outs* sec) { @@ -841,7 +838,7 @@ calculate_k(inputs* data, { // Hack: Initialize coefficients setup_const(); - setup_crown_const(data); + // Aux float ros, bros, lb, fros; int FMC; @@ -857,8 +854,7 @@ calculate_k(inputs* data, FMC = args->FMC; ptr->nftype = data->nftype; ptr->fmc = fmcs[data->nftype][0]; - ptr->cbh = data->cbh; - + ptr->cbh = cbhs[data->nftype][0]; // cout << " cbh " << ptr->cbh << "\n"; ptr->fl = fls_david[data->nftype][0]; @@ -900,11 +896,10 @@ 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)) + if (args->AllowCROS && (data->cbh != 0 || cbhs[data->nftype][0] != 0)) { if (activeCrown) - { - // si el fuego esta activo en copas chequeamos condiciones + { // si el fuego esta activo en copas chequeamos condiciones at->ros_active = active_rate_of_spreadPL04(data, at); if (!checkActive(data, at, FMC)) { @@ -1021,9 +1016,7 @@ determine_destiny_metrics_k(inputs* data, fuel_coefs* ptr, arguments* args, main { // Hack: Initialize coefficients setup_const(); - setup_crown_const(data); - ptr->cbh = data->cbh; // Aux float ros = 0, bros = 0, lb = 0, fros = 0; int FMC = args->FMC; @@ -1039,7 +1032,7 @@ determine_destiny_metrics_k(inputs* data, fuel_coefs* ptr, arguments* args, main metrics->fl = flame_length(data, metrics); // Step 10: Criterion for Crown Fire Initiation (no init if user does not // want to include it) - if (args->AllowCROS && data->cbh > 0) + if (args->AllowCROS) { crownFire = fire_type(data, metrics, FMC); if (crownFire) From 8044884c0d11d84d5be5783dfe5f2af080b08002 Mon Sep 17 00:00:00 2001 From: mati Date: Tue, 28 Jan 2025 16:28:07 -0300 Subject: [PATCH 06/10] 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 adab7611..19253712 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -673,7 +673,7 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) this->surfaceIntensityFolder = Cell2Fire::createOutputFolder("SurfaceIntensity"); } // 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")) { this->crownIntensityFolder = Cell2Fire::createOutputFolder("CrownIntensity"); } @@ -683,12 +683,12 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) this->surfaceFlameLengthFolder = Cell2Fire::createOutputFolder("SurfaceFlameLength"); } // 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")) { this->crownFlameLengthFolder = Cell2Fire::createOutputFolder("CrownFlameLength"); } // 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")) { this->maxFlameLengthFolder = Cell2Fire::createOutputFolder("MaxFlameLength"); } @@ -1749,7 +1749,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; @@ -1785,7 +1785,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; @@ -1803,7 +1803,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 d58d8ae7..6720f41b 100644 --- a/Cell2Fire/Cells.cpp +++ b/Cell2Fire/Cells.cpp @@ -677,7 +677,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 7170c390..7e742ed1 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 c00ba317e876e6ddf0429b7fc1de08e7d06cb8f6 Mon Sep 17 00:00:00 2001 From: mati Date: Thu, 30 Jan 2025 16:17:45 -0300 Subject: [PATCH 07/10] [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, " From f5c063a867ec46e87caa3bc8069b535603217b96 Mon Sep 17 00:00:00 2001 From: "matias.vilches@usach.cl" Date: Wed, 5 Feb 2025 17:17:05 -0300 Subject: [PATCH 08/10] new-out-fl-format --- Cell2Fire/Cell2Fire.cpp | 66 +++++++++++++++++++---------------------- Cell2Fire/Cell2Fire.h | 1 + Cell2Fire/WriteCSV.cpp | 33 +++++++++++++++++++++ Cell2Fire/WriteCSV.h | 1 + 4 files changed, 66 insertions(+), 35 deletions(-) diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index 19253712..0231fc0a 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -49,6 +49,10 @@ std::unordered_map> HarvestedCells; std::vector NFTypesCells; std::unordered_map IgnitionHistory; std::unordered_map WeatherHistory; +std::unordered_map> Statistics; +std::unordered_map meanSurfaceFlameLength; +std::unordered_map meanCrownFlameLength; +std::unordered_map meanMaxFlameLength; /****************************************************************************** Utils @@ -176,9 +180,6 @@ Cell2Fire::Cell2Fire(arguments _args) : CSVForest(_args.InFolder + "fuels", " ") this->args = _args; this->args_ptr = &this->args; - // this->wdf_ptr = new weatherDF; - // this->wdf[100000]; - /******************************************************************** * * Initialization @@ -341,7 +342,6 @@ Cell2Fire::Cell2Fire(arguments _args) : CSVForest(_args.InFolder + "fuels", " ") this->harvestCells.insert(i + 1); } - // POTENCIALMENTE AQUI ESTA MALO /* Weather DataFrame */ // this->WeatherData = this->CSVWeather.getData(); // std::cout << "\nWeather DataFrame from instance " << this->CSVWeather.fileName << std::endl; @@ -1766,7 +1766,8 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->crownIntensities); } - // Surface Flame Length + // Flame Length Statistics + // This will be used to store accumulated flame lengths across simulations if (this->args.OutFl) { this->surfaceFlameLengthFolder = this->args.OutFolder + "SurfaceFlameLength" + separator(); @@ -1784,40 +1785,35 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->surfaceFlameLengths); } - // Crown Flame length - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) - { - this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator(); - std::string fileName; - std::ostringstream oss; - oss.str(""); - oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; - fileName = this->crownFlameLengthFolder + "CrownFlameLength" + oss.str() + ".asc"; - if (this->args.verbose) + for (int cell : this->burntCells) { - std::cout << "We are generating the Crown Flame Length to a asc file " << fileName << std::endl; - } - CSVWriter CSVPloter(fileName, " "); - CSVPloter.printASCII( - this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->crownFlameLengths); - } + std::vector flameLengthMeans; - // Max Flame length - if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) - { - this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength" + separator(); - std::string fileName; - std::ostringstream oss; - oss.str(""); - oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; - fileName = this->maxFlameLengthFolder + "MaxFlameLength" + oss.str() + ".asc"; - if (this->args.verbose) + float surfaceFlameLength = this->surfaceFlameLengths[cell] / args.TotalSims; + meanSurfaceFlameLength[cell] += surfaceFlameLength; + flameLengthMeans.push_back(meanSurfaceFlameLength[cell]); + + if ((this->args.AllowCROS) && (this->args.Simulator != "C")) + { + float crownFlameLength = this->crownFlameLengths[cell] / args.TotalSims; + float maxFlameLength = this->maxFlameLengths[cell] / args.TotalSims; + meanCrownFlameLength[cell] += crownFlameLength; + meanMaxFlameLength[cell] += maxFlameLength; + flameLengthMeans.push_back(meanCrownFlameLength[cell]); + flameLengthMeans.push_back(meanMaxFlameLength[cell]); + } + + Statistics[cell] = { flameLengthMeans }; + } + if (currentSim == args.TotalSims) { - std::cout << "We are generating the Max Flame Length to a asc file " << fileName << std::endl; + std::string filename = "statistics.csv"; + CSVWriter statisticsFolder("", ""); + this->statsFolder = this->args.OutFolder + "Statistics" + separator(); + statisticsFolder.MakeDir(this->statsFolder); + CSVWriter statsFile(this->statsFolder + filename); + statsFile.printStats(Statistics); } - CSVWriter CSVPloter(fileName, " "); - CSVPloter.printASCII( - this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->maxFlameLengths); } // Intensity diff --git a/Cell2Fire/Cell2Fire.h b/Cell2Fire/Cell2Fire.h index 7e91d382..8d869bd2 100644 --- a/Cell2Fire/Cell2Fire.h +++ b/Cell2Fire/Cell2Fire.h @@ -89,6 +89,7 @@ class Cell2Fire string crownFlameLengthFolder; string historyFolder; string ignitionsFolder; + string statsFolder; // Vectors std::vector fire_period; diff --git a/Cell2Fire/WriteCSV.cpp b/Cell2Fire/WriteCSV.cpp index 46b980e2..5fd77df7 100644 --- a/Cell2Fire/WriteCSV.cpp +++ b/Cell2Fire/WriteCSV.cpp @@ -287,6 +287,39 @@ CSVWriter::printIgnitions(std::unordered_map ignitionsHistory, ofs.close(); } +void +CSVWriter::printStats(std::unordered_map> statistics) +{ + std::ofstream ofs(fileName, std::ofstream::out); + + // Print column titles + ofs << "cell,surfaceFlameLength,crownFlameLength,maxFlameLength\n"; + + // Iterate through the unordered_map using iterator + for (auto it = statistics.begin(); it != statistics.end(); ++it) + { + int sim = it->first; // The simulation ID (key) + const std::vector& flameLengths = it->second; // The associated flame lengths (value) + + ofs << sim; // Print the simulation number (key) + + // Print all the flame lengths associated with the current simulation + for (size_t j = 0; j < flameLengths.size(); j++) + { + ofs << "," << flameLengths[j]; + if (j < flameLengths.size() - 1) + { + ofs << " "; // Optional: space between multiple flame lengths + } + } + + ofs << "\n"; // Move to the next line for each entry + } + + // Close the file + ofs.close(); +} + void CSVWriter::printCSV_V2(int rows, int cols, std::vector statusCells) { diff --git a/Cell2Fire/WriteCSV.h b/Cell2Fire/WriteCSV.h index 356ffd81..dda4aacf 100644 --- a/Cell2Fire/WriteCSV.h +++ b/Cell2Fire/WriteCSV.h @@ -45,6 +45,7 @@ class CSVWriter void printWeather(std::vector weatherHistory); void printIgnitions(std::unordered_map ignitionsHistory, std::unordered_map weatherHistory); + void printStats(std::unordered_map> statistics); void printCSV_V2(int rows, int cols, std::vector statusCells); // Function to create a directory void MakeDir(std::string pathPlot); From 6b7ded73ef9856c254d51fc01fd3b8c7f425fd5e Mon Sep 17 00:00:00 2001 From: mati Date: Tue, 1 Jul 2025 21:25:08 -0400 Subject: [PATCH 09/10] Add stats flag --- Cell2Fire/Cell2Fire.cpp | 56 +++++++++++++-- Cell2Fire/FuelModelKitral.cpp | 125 ++++++++++++++-------------------- Cell2Fire/ReadArgs.cpp | 11 +-- 3 files changed, 107 insertions(+), 85 deletions(-) diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index 0231fc0a..c8fb7b97 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -692,6 +692,10 @@ Cell2Fire::reset(int rnumber, double rnumber2, int simExt = 1) { this->maxFlameLengthFolder = Cell2Fire::createOutputFolder("MaxFlameLength"); } + if (this->args.Stats) + { + this->statsFolder = Cell2Fire::createOutputFolder("Statistics"); + } // Crown Folder if (this->args.OutCrown && this->args.AllowCROS) { @@ -1766,8 +1770,7 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->crownIntensities); } - // Flame Length Statistics - // This will be used to store accumulated flame lengths across simulations + // Surface Flame Length if (this->args.OutFl) { this->surfaceFlameLengthFolder = this->args.OutFolder + "SurfaceFlameLength" + separator(); @@ -1785,6 +1788,47 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->surfaceFlameLengths); } + // Crown Flame length + if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) + { + this->crownFlameLengthFolder = this->args.OutFolder + "CrownFlameLength" + separator(); + std::string fileName; + std::ostringstream oss; + oss.str(""); + oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; + fileName = this->crownFlameLengthFolder + "CrownFlameLength" + oss.str() + ".asc"; + if (this->args.verbose) + { + std::cout << "We are generating the Crown Flame Length to a asc file " << fileName << std::endl; + } + CSVWriter CSVPloter(fileName, " "); + CSVPloter.printASCII( + this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->crownFlameLengths); + } + + // Max Flame length + if ((this->args.OutFl) && (this->args.AllowCROS) && (this->args.Simulator != "C")) + { + this->maxFlameLengthFolder = this->args.OutFolder + "MaxFlameLength" + separator(); + std::string fileName; + std::ostringstream oss; + oss.str(""); + oss << std::setfill('0') << std::setw(this->widthSims) << this->sim; + fileName = this->maxFlameLengthFolder + "MaxFlameLength" + oss.str() + ".asc"; + if (this->args.verbose) + { + std::cout << "We are generating the Max Flame Length to a asc file " << fileName << std::endl; + } + CSVWriter CSVPloter(fileName, " "); + CSVPloter.printASCII( + this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->maxFlameLengths); + } + + // Flame Length Statistics + // This will be used to store accumulated flame lengths across simulations + if (this->args.Stats) + { + for (int cell : this->burntCells) { std::vector flameLengthMeans; @@ -1807,11 +1851,9 @@ Cell2Fire::Results() } if (currentSim == args.TotalSims) { - std::string filename = "statistics.csv"; - CSVWriter statisticsFolder("", ""); - this->statsFolder = this->args.OutFolder + "Statistics" + separator(); - statisticsFolder.MakeDir(this->statsFolder); - CSVWriter statsFile(this->statsFolder + filename); + std::ostringstream oss; + std::string Statsname = this->statsFolder + "statistics" + oss.str() + ".csv"; + CSVWriter statsFile(Statsname); statsFile.printStats(Statistics); } } diff --git a/Cell2Fire/FuelModelKitral.cpp b/Cell2Fire/FuelModelKitral.cpp index 89f0cee2..09832e4e 100644 --- a/Cell2Fire/FuelModelKitral.cpp +++ b/Cell2Fire/FuelModelKitral.cpp @@ -559,6 +559,23 @@ setup_const() cbds.insert(std::make_pair(DX02, cbd_dx02)); } +/** + * @brief Sets default value for crown constants when no raster files are provided for them. + * @param data Cell data + */ +void +setup_crown_const(inputs* data) +{ + if (data->cbd == -9999) + { + data->cbd = cbds[data->nftype][0]; + } + if (data->cbh == -9999) + { + data->cbh = cbhs[data->nftype][0]; + } +} + float rate_of_spread_k(inputs* data, fuel_coefs* ptr, @@ -595,6 +612,7 @@ rate_of_spread_k(inputs* data, return at->rss * (at->rss >= 0); } +// TODO: citation needed float flankfire_ros_k(float ros, float bros, float lb) { @@ -602,6 +620,7 @@ flankfire_ros_k(float ros, float bros, float lb) } /* ----------------- Length-to-Breadth --------------------------*/ +// TODO: citation needed float l_to_b(float ws, fuel_coefs* ptr) { @@ -613,6 +632,7 @@ l_to_b(float ws, fuel_coefs* ptr) } /* ----------------- Back Rate of Spread --------------------------*/ +// TODO: citation needed float backfire_ros_k(main_outs* at, snd_outs* sec) { @@ -626,6 +646,7 @@ backfire_ros_k(main_outs* at, snd_outs* sec) return bros * (bros >= 0); } +// TODO: citation needed float slope_effect(float elev_i, float elev_j, int cellsize) { @@ -636,6 +657,7 @@ slope_effect(float elev_i, float elev_j, int cellsize) return se; } +// TODO: citation needed float flame_length(inputs* data, main_outs* at) // REVISAR ESTA ECUACI�N { @@ -647,6 +669,7 @@ flame_length(inputs* data, main_outs* at) // REVISAR ESTA ECUACI�N return fl; } +// TODO: citation needed float angleFL(inputs* data, main_outs* at) { @@ -659,6 +682,7 @@ angleFL(inputs* data, main_outs* at) return angle; } +// TODO: citation needed float flame_height(inputs* data, main_outs* at) { @@ -668,6 +692,7 @@ flame_height(inputs* data, main_outs* at) return fh; } +// TODO: citation needed float byram_intensity(inputs* data, main_outs* at) { @@ -693,14 +718,8 @@ 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; - } + cbd = data->cbd; + cbh = data->cbh; // CBH is 0,1012 * Height^(1,4822) float canopy_height = std::pow(cbh / 0.1012, 1 / 1.4822) - cbh; if (canopy_height < 0) @@ -708,14 +727,6 @@ crown_byram_intensity_k(main_outs* at, inputs* data) 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; } @@ -745,48 +756,28 @@ fire_type(inputs* data, main_outs* at, int FMC) float intensity, critical_intensity, cbh; bool crownFire = false; intensity = at->sfi; - if (std::isnan(data->cbh)) - { - cbh = cbhs[data->nftype][0]; - } - else - { - cbh = data->cbh; - } - // cbh = data->cbh; + cbh = data->cbh; critical_intensity = pow((0.01 * cbh * (460 + 25.9 * FMC)), 1.5); - if ((intensity > critical_intensity) && cbh != 0) + if ((intensity > critical_intensity) && cbh > 0) crownFire = true; return crownFire; } +// TODO: citation needed float crownfractionburn(inputs* data, main_outs* at, int FMC) -{ // generar output de cfb +{ + // generar output de cfb float a, cbd, ros, ros0, H, wa, i0, cbh, cfb; - if (std::isnan(data->cbh)) - { - cbh = cbhs[data->nftype][0]; - } - else - { - cbh = data->cbh; - } + 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]; - if (std::isnan(data->cbd)) - { - cbd = cbds[data->nftype][0]; - } - else - { - cbd = data->cbd; - } + cbd = data->cbd; ros0 = 60 * i0 / (H * wa); ros = at->rss; - if (cbd != 0) + if (cbd > 0) { a = -log(0.1) / (0.9 * (3.0 / cbd - ros0)); } @@ -799,6 +790,7 @@ crownfractionburn(inputs* data, main_outs* at, int FMC) return cfb; } +// TODO: citation needed float active_rate_of_spreadPL04(inputs* data, main_outs* at) // En KITRAL SE USA PL04 @@ -829,10 +821,11 @@ active_rate_of_spreadPL04(inputs* data, rospl04 = fmc * fch * (fv + ps); } ros_active = 3.34 * rospl04; // if rac*cbd>3.0, aplicar - // ros_final=3.34*rospl04 + // ros_final=3.34*rospl04 return ros_active; } +// TODO: citation needed float final_rate_of_spreadPL04(main_outs* at) // En KITRAL SE USA PL04 { @@ -848,30 +841,19 @@ checkActive(inputs* data, main_outs* at, int FMC) // En KITRAL SE USA PL04 { float ros_critical, cbd, H, wa, i0, cbh; bool active; - if (std::isnan(data->cbh)) - { - cbh = cbhs[data->nftype][0]; - } - else - { - cbh = data->cbh; - } + 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]; ros_critical = 60 * i0 / (H * wa); - if (std::isnan(data->cbd)) - { - cbd = cbds[data->nftype][0]; - } - else - { - cbd = data->cbd; - } + cbd = data->cbd; + active = cbd * ros_critical >= 3; return active; } +// TODO: citation needed float backfire_ros10_k(fire_struc* hptr, snd_outs* sec) { @@ -904,7 +886,7 @@ calculate_k(inputs* data, { // Hack: Initialize coefficients setup_const(); - + setup_crown_const(data); // Aux float ros, bros, lb, fros; int FMC; @@ -920,14 +902,8 @@ calculate_k(inputs* data, FMC = args->FMC; ptr->nftype = data->nftype; ptr->fmc = fmcs[data->nftype][0]; - if (isnan(data->cbh)) - { - ptr->cbh = cbhs[data->nftype][0]; - } - else - { - ptr->cbh = data->cbh; - } + ptr->cbh = data->cbh; + // cout << " cbh " << ptr->cbh << "\n"; ptr->fl = fls_david[data->nftype][0]; @@ -969,10 +945,11 @@ 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)) { if (activeCrown) - { // si el fuego esta activo en copas chequeamos condiciones + { + // si el fuego esta activo en copas chequeamos condiciones at->ros_active = active_rate_of_spreadPL04(data, at); if (!checkActive(data, at, FMC)) { @@ -1093,7 +1070,9 @@ determine_destiny_metrics_k(inputs* data, fuel_coefs* ptr, arguments* args, main { // Hack: Initialize coefficients setup_const(); + setup_crown_const(data); + ptr->cbh = data->cbh; // Aux float ros = 0, bros = 0, lb = 0, fros = 0; int FMC = args->FMC; @@ -1109,7 +1088,7 @@ determine_destiny_metrics_k(inputs* data, fuel_coefs* ptr, arguments* args, main metrics->fl = flame_length(data, metrics); // Step 10: Criterion for Crown Fire Initiation (no init if user does not // want to include it) - if (args->AllowCROS) + if (args->AllowCROS && data->cbh > 0) { crownFire = fire_type(data, metrics, FMC); if (crownFire) diff --git a/Cell2Fire/ReadArgs.cpp b/Cell2Fire/ReadArgs.cpp index e2fc172c..0a308dad 100644 --- a/Cell2Fire/ReadArgs.cpp +++ b/Cell2Fire/ReadArgs.cpp @@ -206,16 +206,17 @@ parseArgs(int argc, char* argv[], arguments* args_ptr) //} //--statistics - // if (cmdOptionExists(argv, argv + argc, "--statistics")) { - // out_stats = true; - // printf("Statistics: %d \n", out_stats); - //} + if (cmdOptionExists(argv, argv + argc, "--statistics")) + { + out_stats = true; + printf("Statistics: %d \n", out_stats); + } //--bbo if (cmdOptionExists(argv, argv + argc, "--bbo")) { bbo_tuning = true; - printf("BBOTuning: %s \n", btoa(out_stats)); + printf("BBOTuning: %s \n", btoa(bbo_tuning)); } //--cros From bee942970bf739004a8ef8e86a7f8f39f57b7c6c Mon Sep 17 00:00:00 2001 From: mati Date: Mon, 7 Jul 2025 16:26:35 -0400 Subject: [PATCH 10/10] Add flame length summaries --- Cell2Fire/Cell2Fire.cpp | 41 ++++++++++++++++++++++++++++++----------- Cell2Fire/ReadArgs.cpp | 6 +++--- Cell2Fire/WriteCSV.cpp | 13 +++++++++++-- Cell2Fire/WriteCSV.h | 2 +- 4 files changed, 45 insertions(+), 17 deletions(-) diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index c8fb7b97..ce1f64db 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -49,10 +49,11 @@ std::unordered_map> HarvestedCells; std::vector NFTypesCells; std::unordered_map IgnitionHistory; std::unordered_map WeatherHistory; -std::unordered_map> Statistics; +std::unordered_map> StatisticsPerCell; std::unordered_map meanSurfaceFlameLength; std::unordered_map meanCrownFlameLength; std::unordered_map meanMaxFlameLength; +std::unordered_map> StatisticsPerSim; /****************************************************************************** Utils @@ -1824,18 +1825,23 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->maxFlameLengths); } - // Flame Length Statistics + // Flame Length StatisticsPerCell // This will be used to store accumulated flame lengths across simulations if (this->args.Stats) { - + float totalCrown = 0; + float totalSurface = 0; + float maxPerSim = 0; + std::vector simStats; for (int cell : this->burntCells) { - std::vector flameLengthMeans; + std::vector cellFlameLengthMeans; float surfaceFlameLength = this->surfaceFlameLengths[cell] / args.TotalSims; meanSurfaceFlameLength[cell] += surfaceFlameLength; - flameLengthMeans.push_back(meanSurfaceFlameLength[cell]); + cellFlameLengthMeans.push_back(meanSurfaceFlameLength[cell]); + totalSurface += this->surfaceFlameLengths[cell]; + maxPerSim = max(maxPerSim, this->surfaceFlameLengths[cell]); if ((this->args.AllowCROS) && (this->args.Simulator != "C")) { @@ -1843,18 +1849,31 @@ Cell2Fire::Results() float maxFlameLength = this->maxFlameLengths[cell] / args.TotalSims; meanCrownFlameLength[cell] += crownFlameLength; meanMaxFlameLength[cell] += maxFlameLength; - flameLengthMeans.push_back(meanCrownFlameLength[cell]); - flameLengthMeans.push_back(meanMaxFlameLength[cell]); + cellFlameLengthMeans.push_back(meanCrownFlameLength[cell]); + cellFlameLengthMeans.push_back(meanMaxFlameLength[cell]); + totalCrown += this->maxFlameLengths[cell]; + maxPerSim = max(maxPerSim, this->crownFlameLengths[cell]); } - - Statistics[cell] = { flameLengthMeans }; + StatisticsPerCell[cell] = { cellFlameLengthMeans }; + } + simStats.push_back(totalSurface / BCells); + if ((this->args.AllowCROS) && (this->args.Simulator != "C")) + { + simStats.push_back(totalCrown / BCells); } + simStats.push_back(maxPerSim); + + StatisticsPerSim[this->sim] = { simStats }; + std::ostringstream oss; + std::string Stats = this->statsFolder + "statisticsPerSim" + oss.str() + ".csv"; + CSVWriter simStatsFile(Stats); + simStatsFile.printStats(StatisticsPerSim, "sim", (this->args.AllowCROS) && (this->args.Simulator != "C")); if (currentSim == args.TotalSims) { std::ostringstream oss; - std::string Statsname = this->statsFolder + "statistics" + oss.str() + ".csv"; + std::string Statsname = this->statsFolder + "statisticsPerCell" + oss.str() + ".csv"; CSVWriter statsFile(Statsname); - statsFile.printStats(Statistics); + statsFile.printStats(StatisticsPerCell, "cell", (this->args.AllowCROS) && (this->args.Simulator != "C")); } } diff --git a/Cell2Fire/ReadArgs.cpp b/Cell2Fire/ReadArgs.cpp index 0a308dad..90268667 100644 --- a/Cell2Fire/ReadArgs.cpp +++ b/Cell2Fire/ReadArgs.cpp @@ -209,7 +209,7 @@ parseArgs(int argc, char* argv[], arguments* args_ptr) if (cmdOptionExists(argv, argv + argc, "--statistics")) { out_stats = true; - printf("Statistics: %d \n", out_stats); + printf("StatisticsPerCell: %d \n", out_stats); } //--bbo @@ -629,7 +629,7 @@ printArgs(arguments args) std::cout << "FinalGrid: " << args.FinalGrid << std::endl; std::cout << "PromTuned: " << args.PromTuned << std::endl; std::cout << "BBOTuning: " << args.BBOTuning << std::endl; - std::cout << "Statistics: " << args.Stats << std::endl; + std::cout << "StatisticsPerCell: " << args.Stats << std::endl; std::cout << "noOutput: " << args.NoOutput << std::endl; std::cout << "verbose: " << args.verbose << std::endl; std::cout << "seed: " << args.seed << std::endl; @@ -661,7 +661,7 @@ printArgs(arguments args) std::cout << "FinalGrid: " << args.FinalGrid << std::endl; std::cout << "PromTuned: " << args.PromTuned << std::endl; std::cout << "BBOTuning: " << args.BBOTuning << std::endl; - std::cout << "Statistics: " << args.Stats << std::endl; + std::cout << "StatisticsPerCell: " << args.Stats << std::endl; std::cout << "noOutput: " << args.NoOutput << std::endl; std::cout << "verbose: " << args.verbose << std::endl; std::cout << "Ignition Points Log: " << args.IgnitionsLog << std::endl; diff --git a/Cell2Fire/WriteCSV.cpp b/Cell2Fire/WriteCSV.cpp index 5fd77df7..eaa70be7 100644 --- a/Cell2Fire/WriteCSV.cpp +++ b/Cell2Fire/WriteCSV.cpp @@ -288,12 +288,21 @@ CSVWriter::printIgnitions(std::unordered_map ignitionsHistory, } void -CSVWriter::printStats(std::unordered_map> statistics) +CSVWriter::printStats(std::unordered_map> statistics, std::string type, bool hasCrown) { std::ofstream ofs(fileName, std::ofstream::out); // Print column titles - ofs << "cell,surfaceFlameLength,crownFlameLength,maxFlameLength\n"; + ofs << type << ",surfaceFlameLengthMean"; + if (hasCrown) + { + ofs << type << ",crownFlameLengthMean,maxFlameLength"; + } + else if (type == "sim") + { + ofs << ",maxFlameLength"; + } + ofs << std::endl; // Iterate through the unordered_map using iterator for (auto it = statistics.begin(); it != statistics.end(); ++it) diff --git a/Cell2Fire/WriteCSV.h b/Cell2Fire/WriteCSV.h index dda4aacf..74e77af3 100644 --- a/Cell2Fire/WriteCSV.h +++ b/Cell2Fire/WriteCSV.h @@ -45,7 +45,7 @@ class CSVWriter void printWeather(std::vector weatherHistory); void printIgnitions(std::unordered_map ignitionsHistory, std::unordered_map weatherHistory); - void printStats(std::unordered_map> statistics); + void printStats(std::unordered_map> statistics, std::string type, bool hasCrown); void printCSV_V2(int rows, int cols, std::vector statusCells); // Function to create a directory void MakeDir(std::string pathPlot);