Don't calculcate PMCoarse when pm10 and pm2.5 are not included in the…

… run
VITObelgium · Mar 16, 2023 · a06a47f · a06a47f
1 parent e48acf0
commit a06a47f
Show file tree

Hide file tree

Showing 4 changed files with 62 additions and 36 deletions.
diff --git a/logic/emissioninventory.cpp b/logic/emissioninventory.cpp
@@ -403,47 +403,49 @@ SingleEmissions read_country_point_sources(const RunConfiguration& cfg, const Co
             }
         }
 
-        try {
-            const auto pm10     = cfg.pollutants().try_pollutant_from_string(constants::pollutant::PM10);
-            const auto pm2_5    = cfg.pollutants().try_pollutant_from_string(constants::pollutant::PM2_5);
-            const auto pmCoarse = cfg.pollutants().try_pollutant_from_string(constants::pollutant::PMCoarse);
-
-            if (pm10.has_value() && pm2_5.has_value() && pmCoarse.has_value()) {
-                // Calculate PMcoarse data from pm10 and pm2.5 data
-                const auto pm10Emissions = read_country_pollutant_point_sources(pointEmissionsDir, *pm10, cfg, runSummary);
-                auto pm2_5Emissions      = read_country_pollutant_point_sources(pointEmissionsDir, *pm2_5, cfg, runSummary);
-
-                std::sort(pm2_5Emissions.begin(), pm2_5Emissions.end(), [](const EmissionEntry& lhs, const EmissionEntry& rhs) {
-                    return lhs.source_id() < rhs.source_id();
-                });
-
-                for (auto& pm10Entry : pm10Emissions) {
-                    auto pm10Val = pm10Entry.value().amount();
-                    if (pm10Val.has_value()) {
-                        auto iter = std::lower_bound(pm2_5Emissions.begin(), pm2_5Emissions.end(), pm10Entry.source_id(), [](const EmissionEntry& entry, std::string_view srcId) {
-                            return entry.source_id() < srcId;
-                        });
-
-                        double pm2_5Val = 0.0;
-                        if (iter != pm2_5Emissions.end() && iter->source_id() == pm10Entry.source_id()) {
-                            pm2_5Val = iter->value().amount().value_or(0.0);
-                            assert(iter->coordinate() == pm10Entry.coordinate());
-                        }
+        if (cfg.pmcoarse_calculation_needed()) {
+            try {
+                const auto pm10     = cfg.pollutants().try_pollutant_from_string(constants::pollutant::PM10);
+                const auto pm2_5    = cfg.pollutants().try_pollutant_from_string(constants::pollutant::PM2_5);
+                const auto pmCoarse = cfg.pollutants().try_pollutant_from_string(constants::pollutant::PMCoarse);
+
+                if (pm10.has_value() && pm2_5.has_value() && pmCoarse.has_value()) {
+                    // Calculate PMcoarse data from pm10 and pm2.5 data
+                    const auto pm10Emissions = read_country_pollutant_point_sources(pointEmissionsDir, *pm10, cfg, runSummary);
+                    auto pm2_5Emissions      = read_country_pollutant_point_sources(pointEmissionsDir, *pm2_5, cfg, runSummary);
+
+                    std::sort(pm2_5Emissions.begin(), pm2_5Emissions.end(), [](const EmissionEntry& lhs, const EmissionEntry& rhs) {
+                        return lhs.source_id() < rhs.source_id();
+                    });
+
+                    for (auto& pm10Entry : pm10Emissions) {
+                        auto pm10Val = pm10Entry.value().amount();
+                        if (pm10Val.has_value()) {
+                            auto iter = std::lower_bound(pm2_5Emissions.begin(), pm2_5Emissions.end(), pm10Entry.source_id(), [](const EmissionEntry& entry, std::string_view srcId) {
+                                return entry.source_id() < srcId;
+                            });
+
+                            double pm2_5Val = 0.0;
+                            if (iter != pm2_5Emissions.end() && iter->source_id() == pm10Entry.source_id()) {
+                                pm2_5Val = iter->value().amount().value_or(0.0);
+                                assert(iter->coordinate() == pm10Entry.coordinate());
+                            }
 
-                        try {
-                            if (auto pmCoarseVal = EmissionValue(pmcoarse_from_pm25_pm10(pm2_5Val, pm10Val)); pmCoarseVal.amount().has_value()) {
-                                EmissionEntry entry(EmissionIdentifier(country, pm10Entry.id().sector, *pmCoarse), pmCoarseVal);
-                                entry.set_coordinate(pm10Entry.coordinate().value());
-                                result.add_emission(std::move(entry));
+                            try {
+                                if (auto pmCoarseVal = EmissionValue(pmcoarse_from_pm25_pm10(pm2_5Val, pm10Val)); pmCoarseVal.amount().has_value()) {
+                                    EmissionEntry entry(EmissionIdentifier(country, pm10Entry.id().sector, *pmCoarse), pmCoarseVal);
+                                    entry.set_coordinate(pm10Entry.coordinate().value());
+                                    result.add_emission(std::move(entry));
+                                }
+                            } catch (const std::exception& e) {
+                                throw RuntimeError("Sector {} with EIL nr {} ({})", pm10Entry.id().sector, pm10Entry.source_id(), e.what());
                             }
-                        } catch (const std::exception& e) {
-                            throw RuntimeError("Sector {} with EIL nr {} ({})", pm10Entry.id().sector, pm10Entry.source_id(), e.what());
                         }
                     }
                 }
+            } catch (const std::exception& e) {
+                Log::debug("Failed to create PMcoarse point sources: {}", e.what());
             }
-        } catch (const std::exception& e) {
-            Log::debug("Failed to create PMcoarse point sources: {}", e.what());
         }
 
         const auto flandersMeta   = grid_data(GridDefinition::Flanders1km).meta;

diff --git a/logic/include/emap/runconfiguration.h b/logic/include/emap/runconfiguration.h
@@ -65,6 +65,8 @@ class RunConfiguration
     fs::path emission_output_raster_path(date::year year, const EmissionIdentifier& emissionId) const;
     fs::path emission_brn_output_path(date::year year, const Pollutant& pol, const EmissionSector& sector) const;
 
+    bool pmcoarse_calculation_needed() const noexcept;
+
     const fs::path& data_root() const noexcept;
     void set_data_root(const fs::path& root);
 
@@ -92,6 +94,7 @@ class RunConfiguration
     std::optional<int32_t> max_concurrency() const noexcept;
 
     std::vector<Pollutant> included_pollutants() const;
+    bool pollutant_is_included(std::string_view pollutant) const noexcept;
 
     const SectorInventory& sectors() const noexcept;
     const PollutantInventory& pollutants() const noexcept;

diff --git a/logic/inputparsers.cpp b/logic/inputparsers.cpp
@@ -387,7 +387,9 @@ SingleEmissions parse_emissions(EmissionSector::Type sectorType, const fs::path&
         }
 
         SingleEmissions emissions(requestYear, std::move(entries));
-        merge_unique_emissions(emissions, calculate_pmcoarse_emissions(pollutantInv, emissions));
+        if (cfg.pmcoarse_calculation_needed()) {
+            merge_unique_emissions(emissions, calculate_pmcoarse_emissions(pollutantInv, emissions));
+        }
         return emissions;
     } catch (const std::exception& e) {
         throw RuntimeError("Error parsing {} ({})", emissionsCsv, e.what());

diff --git a/logic/runconfiguration.cpp b/logic/runconfiguration.cpp
@@ -1,6 +1,7 @@
 #include "emap/runconfiguration.h"
 
 #include "configurationutil.h"
+#include "emap/constants.h"
 #include "infra/exception.h"
 #include "infra/string.h"
 
@@ -80,6 +81,12 @@ fs::path RunConfiguration::emission_brn_output_path(date::year year, const Pollu
     return _paths.emission_brn_output_path(year, pol, sector);
 }
 
+bool RunConfiguration::pmcoarse_calculation_needed() const noexcept
+{
+    return pollutant_is_included(constants::pollutant::PM10) &&
+           pollutant_is_included(constants::pollutant::PM2_5);
+}
+
 const fs::path& RunConfiguration::emission_scalings_path() const noexcept
 {
     return _emissionScalingsPath;
@@ -179,6 +186,18 @@ std::vector<Pollutant> RunConfiguration::included_pollutants() const
     return _includedPollutants;
 }
 
+bool RunConfiguration::pollutant_is_included(std::string_view pollutantName) const noexcept
+{
+    auto pol = _pollutantInventory.try_pollutant_from_string(pollutantName);
+
+    if (_includedPollutants.empty()) {
+        // All pollutants are included, if we know the pollutant return true
+        return pol.has_value();
+    }
+
+    return container_contains(_includedPollutants, *pol);
+}
+
 const SectorInventory& RunConfiguration::sectors() const noexcept
 {
     return _sectorInventory;