Emission inventory summary

VITObelgium · Jan 28, 2022 · cf04fb8 · cf04fb8
1 parent 98744d2
commit cf04fb8
Show file tree

Hide file tree

Showing 18 changed files with 333 additions and 169 deletions.
diff --git a/emap.natvis b/emap.natvis
@@ -4,7 +4,7 @@
         <DisplayString>{_name}</DisplayString>
     </Type>
 	<Type Name="emap::Country">
-		<DisplayString>{_isoCode}</DisplayString>
+		<DisplayString>{_isoCode.value_}</DisplayString>
 	</Type>
 	<Type Name="emap::Pollutant">
 		<DisplayString>{_code}</DisplayString>

diff --git a/logic/CMakeLists.txt b/logic/CMakeLists.txt
@@ -23,7 +23,7 @@ add_library(emaplogic
     include/emap/outputwriters.h outputwriters.cpp
     configurationutil.h
     enuminfo.h
-    emissioninventory.h
+    emissioninventory.h emissioninventory.cpp
     unitconversion.h
     geometry.h geometry.cpp
     runsummary.h runsummary.cpp

diff --git a/logic/emissioninventory.cpp b/logic/emissioninventory.cpp
@@ -0,0 +1,147 @@
+#pragma once
+
+#include "emap/emissions.h"
+#include "emap/scalingfactors.h"
+#include "infra/algo.h"
+#include "infra/log.h"
+#include "runsummary.h"
+
+#include <cassert>
+#include <numeric>
+
+namespace emap {
+
+using namespace inf;
+
+EmissionSector convert_sector_to_gnfr_level(const EmissionSector& sec)
+{
+    return EmissionSector(sec.gnfr_sector());
+}
+
+EmissionIdentifier convert_emission_id_to_gnfr_level(const EmissionIdentifier& id)
+{
+    EmissionIdentifier result = id;
+    result.sector             = convert_sector_to_gnfr_level(id.sector);
+    return result;
+}
+
+static std::unordered_map<EmissionIdentifier, double> create_gnfr_sums(const SingleEmissions& totalEmissionsNfr)
+{
+    std::unordered_map<EmissionIdentifier, double> result;
+
+    for (const auto& em : totalEmissionsNfr) {
+        if (em.country().is_belgium() || !em.value().amount().has_value()) {
+            continue;
+        }
+
+        assert(em.sector().type() == EmissionSector::Type::Gnfr);
+        assert(result.count(em.id()) == 0);
+
+        result.emplace(em.id(), *em.value().amount());
+    }
+
+    return result;
+}
+
+static std::unordered_map<EmissionIdentifier, double> create_nfr_sums(const SingleEmissions& totalEmissionsNfr)
+{
+    std::unordered_map<EmissionIdentifier, double> result;
+
+    for (const auto& em : totalEmissionsNfr) {
+        if (em.country().is_belgium() || !em.value().amount().has_value()) {
+            continue;
+        }
+
+        assert(em.sector().type() == EmissionSector::Type::Nfr);
+        EmissionIdentifier gnfrId = convert_emission_id_to_gnfr_level(em.id());
+
+        result[gnfrId] += *em.value().amount();
+    }
+
+    return result;
+}
+
+static std::unordered_map<EmissionIdentifier, double> create_nfr_correction_ratios(
+    const std::unordered_map<EmissionIdentifier, double>& nfrBasedTotals,
+    const std::unordered_map<EmissionIdentifier, double>& gnfrBasedTotals,
+    RunSummary& summary)
+{
+    std::unordered_map<EmissionIdentifier, double> result;
+
+    for (auto& [id, nfrBasedTotal] : nfrBasedTotals) {
+        const auto gnfrBasedTotal = inf::find_in_map_optional(gnfrBasedTotals, id);
+        double correction         = 1.0;
+
+        if (gnfrBasedTotal.has_value()) {
+            correction = *gnfrBasedTotal / nfrBasedTotal;
+            summary.add_gnfr_correction(id, *gnfrBasedTotal, nfrBasedTotal, correction);
+        } else {
+            summary.add_gnfr_correction(id, 0.0, nfrBasedTotal, 1.0);
+        }
+
+        result[id] = correction;
+    }
+
+    return result;
+}
+
+EmissionInventory create_emission_inventory(const SingleEmissions& totalEmissionsNfr,
+                                            const SingleEmissions& totalEmissionsGnfr,
+                                            const SingleEmissions& pointSourceEmissions,
+                                            const ScalingFactors& diffuseScalings,
+                                            const ScalingFactors& pointScalings,
+                                            RunSummary& runSummary)
+{
+    EmissionInventory result;
+
+    // Create gnfr sum of the actual reported nfr values
+    // Create gnfr sum of the validated gnfr values (put in map)
+    // Then calculate the ratio between the two
+    auto nfrCorrectionRatios = create_nfr_correction_ratios(create_nfr_sums(totalEmissionsNfr), create_gnfr_sums(totalEmissionsGnfr), runSummary);
+
+    for (const auto& em : totalEmissionsNfr) {
+        assert(em.sector().type() == EmissionSector::Type::Nfr);
+
+        double diffuseEmission  = em.value().amount().value_or(0.0);
+        double pointEmissionSum = 0.0;
+        std::vector<EmissionEntry> pointSourceEntries;
+
+        if (em.country().is_belgium()) {
+            // For belgian regions we calculate the diffuse emissions by subtracting the point source emissions
+            // from the total emissions
+
+            pointSourceEntries = pointSourceEmissions.emissions_with_id(em.id());
+            pointEmissionSum   = std::accumulate(pointSourceEntries.cbegin(), pointSourceEntries.cend(), 0.0, [](double total, const auto& current) {
+                return total + current.value().amount().value_or(0.0);
+            });
+
+            if (diffuseEmission > 0 && pointEmissionSum > diffuseEmission) {
+                // Check if the difference is caused by floating point rounding
+                if (std::abs(pointEmissionSum - diffuseEmission) < 1e-6) {
+                    // Minor difference caused by rounding, make them the same
+                    pointEmissionSum = diffuseEmission;
+                } else {
+                    throw RuntimeError("The sum of the point emissions ({}) for {} is bigger than the diffuse emissions ({}) for sector {}", pointEmissionSum, em.country(), diffuseEmission, em.sector());
+                }
+            } else {
+            }
+        } else {
+            // Rest of Europe
+            if (diffuseEmission < 0.0) {
+                inf::Log::warn("Negative emissions reported for {}", em.id());
+                diffuseEmission = 0.0;
+            }
+
+            diffuseEmission *= nfrCorrectionRatios[em.id()];
+        }
+
+        EmissionInventoryEntry entry(em.id(), diffuseEmission - pointEmissionSum, std::move(pointSourceEntries));
+        entry.set_diffuse_scaling(diffuseScalings.scaling_for_id(em.id()).value_or(1.0));
+        entry.set_point_scaling(pointScalings.scaling_for_id(em.id()).value_or(1.0));
+        result.add_emission(std::move(entry));
+    }
+
+    return result;
+}
+
+}
diff --git a/logic/emissioninventory.h b/logic/emissioninventory.h
@@ -12,95 +12,13 @@ namespace emap {
 
 using namespace inf;
 
-static std::unordered_map<CountryId, std::unordered_map<GnfrId, double>> create_gnfr_sums(const SingleEmissions& totalEmissionsNfr)
-{
-    std::unordered_map<CountryId, std::unordered_map<GnfrId, double>> result;
+class RunSummary;
 
-    for (const auto& em : totalEmissionsNfr) {
-        if (em.country().is_belgium() || !em.value().amount().has_value()) {
-            continue;
-        }
+EmissionInventory create_emission_inventory(const SingleEmissions& totalEmissionsNfr,
+                                            const SingleEmissions& totalEmissionsGnfr,
+                                            const SingleEmissions& pointSourceEmissions,
+                                            const ScalingFactors& diffuseScalings,
+                                            const ScalingFactors& pointScalings,
+                                            RunSummary& runSummary);
 
-        assert(em.sector().type() == EmissionSector::Type::Nfr);
-        result[em.country().id()][em.sector().gnfr_sector().id()] += *em.value().amount();
-    }
-
-    return result;
-}
-
-static std::unordered_map<CountryId, std::unordered_map<GnfrId, double>> create_nfr_correction_ratios(
-    const std::unordered_map<CountryId, std::unordered_map<GnfrId, double>>& nfrBasedTotals,
-    const std::unordered_map<CountryId, std::unordered_map<GnfrId, double>>& gnfrBasedTotals)
-{
-    std::unordered_map<CountryId, std::unordered_map<GnfrId, double>> result;
-
-    for (auto& [countryId, map] : nfrBasedTotals) {
-        for (auto& [gnfrSector, nfrBasedTotal] : map) {
-            const auto& sectorMap = inf::find_in_map_required(gnfrBasedTotals, countryId);
-            double gnfrBasedTotal = inf::find_in_map_required(sectorMap, gnfrSector);
-
-            result[countryId][gnfrSector] = gnfrBasedTotal / nfrBasedTotal;
-        }
-    }
-
-    return result;
-}
-
-inline EmissionInventory create_emission_inventory(const SingleEmissions& totalEmissionsNfr,
-                                                   const SingleEmissions& totalEmissionsGnfr,
-                                                   const SingleEmissions& pointSourceEmissions,
-                                                   const ScalingFactors& diffuseScalings,
-                                                   const ScalingFactors& pointScalings)
-{
-    EmissionInventory result;
-
-    // Create gnfr sum of the actual reported nfr values
-    // Create gnfr sum of the validated gnfr values (put in map)
-    // Then calculate the ratio between the two
-    auto nfrCorrectionRatios = create_nfr_correction_ratios(create_gnfr_sums(totalEmissionsNfr), create_gnfr_sums(totalEmissionsGnfr));
-
-    for (const auto& em : totalEmissionsNfr) {
-        assert(em.sector().type() == EmissionSector::Type::Nfr);
-
-        double diffuseEmission  = em.value().amount().value_or(0.0);
-        double pointEmissionSum = 0.0;
-        std::vector<EmissionEntry> pointSourceEntries;
-
-        if (em.country().is_belgium()) {
-            // For belgian regions we calculate the diffuse emissions by subtracting the point source emissions
-            // from the total emissions
-
-            pointSourceEntries = pointSourceEmissions.emissions_with_id(em.id());
-            pointEmissionSum   = std::accumulate(pointSourceEntries.cbegin(), pointSourceEntries.cend(), 0.0, [](double total, const auto& current) {
-                return total + current.value().amount().value_or(0.0);
-              });
-
-            if (diffuseEmission > 0 && pointEmissionSum > diffuseEmission) {
-                // Check if the difference is caused by floating point rounding
-                if (std::abs(pointEmissionSum - diffuseEmission) < 1e-6) {
-                    // Minor difference caused by rounding, make them the same
-                    pointEmissionSum = diffuseEmission;
-                } else {
-                    throw RuntimeError("The sum of the point emissions ({}) for {} is bigger than the diffuse emissions ({}) for sector {}", pointEmissionSum, em.country(), diffuseEmission, em.sector());
-                }
-            } else {
-            }
-        } else {
-            // Rest of Europe
-            if (diffuseEmission < 0.0) {
-                inf::Log::warn("Negative emissions reported for {}", em.id());
-                diffuseEmission = 0.0;
-            }
-
-            diffuseEmission *= nfrCorrectionRatios[em.country().id()][em.sector().gnfr_sector().id()];
-        }
-
-        EmissionInventoryEntry entry(em.id(), diffuseEmission - pointEmissionSum, std::move(pointSourceEntries));
-        entry.set_diffuse_scaling(diffuseScalings.scaling_for_id(em.id()).value_or(1.0));
-        entry.set_point_scaling(pointScalings.scaling_for_id(em.id()).value_or(1.0));
-        result.add_emission(std::move(entry));
-    }
-
-    return result;
-}
 }
diff --git a/logic/include/emap/emissions.h b/logic/include/emap/emissions.h
@@ -4,6 +4,7 @@
 #include "emap/pollutant.h"
 #include "emap/sector.h"
 #include "infra/algo.h"
+#include "infra/hash.h"
 #include "infra/point.h"
 #include "infra/span.h"
 
@@ -467,3 +468,16 @@ struct formatter<emap::EmissionSector::Type>
     }
 };
 }
+
+namespace std {
+template <>
+struct hash<emap::EmissionIdentifier>
+{
+    size_t operator()(const emap::EmissionIdentifier& id) const
+    {
+        size_t seed = 0;
+        inf::hash_combine(seed, id.country.id(), id.pollutant.code(), id.sector.id());
+        return seed;
+    }
+};
+}
diff --git a/logic/include/emap/runconfiguration.h b/logic/include/emap/runconfiguration.h
@@ -62,6 +62,8 @@ class RunConfiguration
     ValidationType validation_type() const noexcept;
 
     date::year year() const noexcept;
+    void set_year(date::year year) noexcept;
+
     date::year reporting_year() const noexcept;
 
     std::string_view scenario() const noexcept;

diff --git a/logic/include/emap/runsummary.h b/logic/include/emap/runsummary.h
diff --git a/logic/include/emap/sector.h b/logic/include/emap/sector.h
@@ -160,6 +160,8 @@ class EmissionSector
      */
     const GnfrSector& gnfr_sector() const noexcept;
 
+    int64_t id() const noexcept;
+
     bool is_land_sector() const noexcept;
     /* Returns the nfr sector this sector overrides if it is applicable */
     // std::optional<NfrSector> is_sector_override() const noexcept;

diff --git a/logic/inputparsers.cpp b/logic/inputparsers.cpp
@@ -142,9 +142,14 @@ SingleEmissions parse_emissions(EmissionSector::Type sectorType, const fs::path&
         using namespace io;
         CSVReader<6, trim_chars<' ', '\t'>, no_quote_escape<';'>, throw_on_overflow, single_line_comment<'#'>> in(emissionsCsv.u8string());
 
-        int year;
+        int32_t year;
         char *countryStr, *sector, *pollutant, *unit, *value;
         while (in.read_row(countryStr, year, sector, pollutant, unit, value)) {
+            if (year != static_cast<int32_t>(cfg.year())) {
+                // not the year we want
+                continue;
+            }
+
             auto emissionValue = str::to_double(value);
             if (emissionValue.has_value()) {
                 *emissionValue = to_giga_gram(*emissionValue, unit);

diff --git a/logic/modelrun.cpp b/logic/modelrun.cpp
@@ -352,13 +352,13 @@ void run_model(const RunConfiguration& cfg, const ModelProgress::Callback& progr
 
     Log::debug("Generate emission inventory");
     chrono::DurationRecorder dur;
-    const auto inventory = create_emission_inventory(nfrTotalEmissions, gnfrTotalEmissions, pointSourcesFlanders, scalingsDiffuse, scalingsPointSource);
+    const auto inventory = create_emission_inventory(nfrTotalEmissions, gnfrTotalEmissions, pointSourcesFlanders, scalingsDiffuse, scalingsPointSource, summary);
     Log::debug("Generate emission inventory took {}", dur.elapsed_time_string());
 
-    Log::debug("Spread emissions");
+    /*Log::debug("Spread emissions");
     dur.reset();
     spread_emissions(inventory, spatPatInv, cfg, progressCb);
-    Log::debug("Spread emissions took {}", dur.elapsed_time_string());
+    Log::debug("Spread emissions took {}", dur.elapsed_time_string());*/
 
     // Write the summary
     {
@@ -373,7 +373,7 @@ void run_model(const RunConfiguration& cfg, const ModelProgress::Callback& progr
             fmt::fprintf(fp, summary.spatial_pattern_usage_table());
         }
 
-        summary.write_spatial_pattern_spreadsheet(cfg.run_summary_spreadsheet_path());
+        summary.write_summary_spreadsheet(cfg.run_summary_spreadsheet_path());
     }
 }
 }
diff --git a/logic/runconfiguration.cpp b/logic/runconfiguration.cpp
@@ -151,6 +151,11 @@ date::year RunConfiguration::year() const noexcept
     return _year;
 }
 
+void RunConfiguration::set_year(date::year year) noexcept
+{
+    _year = year;
+}
+
 date::year RunConfiguration::reporting_year() const noexcept
 {
     return _reportYear;