diff --git a/Cell2Fire/Cell2Fire.cpp b/Cell2Fire/Cell2Fire.cpp index adab7611..ce1f64db 100644 --- a/Cell2Fire/Cell2Fire.cpp +++ b/Cell2Fire/Cell2Fire.cpp @@ -49,6 +49,11 @@ std::unordered_map> HarvestedCells; std::vector NFTypesCells; std::unordered_map IgnitionHistory; std::unordered_map WeatherHistory; +std::unordered_map> StatisticsPerCell; +std::unordered_map meanSurfaceFlameLength; +std::unordered_map meanCrownFlameLength; +std::unordered_map meanMaxFlameLength; +std::unordered_map> StatisticsPerSim; /****************************************************************************** Utils @@ -176,9 +181,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 +343,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; @@ -673,7 +674,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,15 +684,19 @@ 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"); } + if (this->args.Stats) + { + this->statsFolder = Cell2Fire::createOutputFolder("Statistics"); + } // Crown Folder if (this->args.OutCrown && this->args.AllowCROS) { @@ -1749,7 +1754,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 +1790,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 +1808,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; @@ -1820,6 +1825,58 @@ Cell2Fire::Results() this->rows, this->cols, this->xllcorner, this->yllcorner, this->cellSide, this->maxFlameLengths); } + // 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 cellFlameLengthMeans; + + float surfaceFlameLength = this->surfaceFlameLengths[cell] / args.TotalSims; + meanSurfaceFlameLength[cell] += surfaceFlameLength; + cellFlameLengthMeans.push_back(meanSurfaceFlameLength[cell]); + totalSurface += this->surfaceFlameLengths[cell]; + maxPerSim = max(maxPerSim, this->surfaceFlameLengths[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; + cellFlameLengthMeans.push_back(meanCrownFlameLength[cell]); + cellFlameLengthMeans.push_back(meanMaxFlameLength[cell]); + totalCrown += this->maxFlameLengths[cell]; + maxPerSim = max(maxPerSim, this->crownFlameLengths[cell]); + } + 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 + "statisticsPerCell" + oss.str() + ".csv"; + CSVWriter statsFile(Statsname); + statsFile.printStats(StatisticsPerCell, "cell", (this->args.AllowCROS) && (this->args.Simulator != "C")); + } + } + // Intensity if ((this->args.OutCrownConsumption) && (this->args.AllowCROS)) { 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/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 582608a0..09832e4e 100644 --- a/Cell2Fire/FuelModelKitral.cpp +++ b/Cell2Fire/FuelModelKitral.cpp @@ -705,6 +705,51 @@ 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; + 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) + { + throw std::runtime_error("Tree height is lower than canopy base height, " + "please provide valid files."); + } + 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) { @@ -935,6 +980,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) { @@ -957,6 +1004,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) { @@ -1045,6 +1094,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); diff --git a/Cell2Fire/ReadArgs.cpp b/Cell2Fire/ReadArgs.cpp index e2fc172c..90268667 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("StatisticsPerCell: %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 @@ -628,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; @@ -660,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 46b980e2..eaa70be7 100644 --- a/Cell2Fire/WriteCSV.cpp +++ b/Cell2Fire/WriteCSV.cpp @@ -287,6 +287,48 @@ CSVWriter::printIgnitions(std::unordered_map ignitionsHistory, ofs.close(); } +void +CSVWriter::printStats(std::unordered_map> statistics, std::string type, bool hasCrown) +{ + std::ofstream ofs(fileName, std::ofstream::out); + + // Print column titles + 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) + { + 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..74e77af3 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, 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);