Skip to content

Commit

Permalink
PMCoarse tests added
Browse files Browse the repository at this point in the history
  • Loading branch information
@dirkvdb
dirkvdb committed Apr 4, 2023
1 parent 974fbf0 commit 54308ee
Show file tree
Hide file tree
Showing 6 changed files with 59 additions and 18 deletions.
2 changes: 1 addition & 1 deletion deps/vcpkg
Submodule vcpkg updated 22 files
+1 −5 ports/gdal/CONTROL
+12 −0 ports/gdal/intel.patch
+2 −2 ports/gdal/portfile.cmake
+0 −34 ports/libiconv/0001-Add-export-definitions.patch
+14 −1 ports/libiconv/0002-Config-for-MSVC.patch
+12 −0 ports/libiconv/0003-Add-export.patch
+0 −25 ports/libiconv/0003-Fix-uwp.patch
+13 −0 ports/libiconv/0004-ModuleFileName.patch
+0 −86 ports/libiconv/CMakeLists.txt
+0 −4 ports/libiconv/CONTROL
+64 −0 ports/libiconv/clang-fortify.patch
+53 −27 ports/libiconv/portfile.cmake
+0 −8 ports/libiconv/unofficial-iconv-config.cmake
+4 −0 ports/libiconv/usage
+20 −11 ports/libiconv/vcpkg-cmake-wrapper.cmake
+8 −0 ports/libiconv/vcpkg.json
+14 −0 ports/proj/inteloneapi.patch
+4 −7 ports/proj/portfile.cmake
+0 −24 ports/proj/version-reject.patch
+10 −4 scripts/toolchains/windows.cmake
+11 −0 triplets/x64-windows-intel.cmake
+4 −1 triplets/x64-windows-static-intel.cmake
33 changes: 18 additions & 15 deletions logic/emissioninventory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,9 @@ template <typename Callback>
static void zip_point_emissions(std::span<const EmissionEntry> pol1Points, std::span<const EmissionEntry> pol2Points, Callback&& callback)
{
auto pol2PointsCopy = container_as_vector(pol2Points);
std::sort(pol2PointsCopy.begin(), pol2PointsCopy.end(), [](const EmissionEntry& entry1, const EmissionEntry& entry2) {
return entry1.source_id() < entry2.source_id();
});

for (auto& pol1Entry : pol1Points) {
if (pol1Entry.value().amount().has_value()) {
Expand All @@ -187,19 +190,19 @@ static void zip_point_emissions(std::span<const EmissionEntry> pol1Points, std::
}
}

static void validate_pm10_pm25(const EmissionInventoryEntry& pm10, const EmissionInventoryEntry& pm25)
static void validate_pm10_pm25(const EmissionInventoryEntry& pm10Emissions, const EmissionInventoryEntry& pm25Emissions)
{
const double diffThreshold = 1e-5;

{
// validate the point emissions, perform validation after the autoscaling, do not take user scalings into account
const auto pm10AutoScale = pm10.point_auto_scaling_factor();
const auto pm25AutoScale = pm25.point_auto_scaling_factor();
const auto pm10AutoScale = pm10Emissions.point_auto_scaling_factor();
const auto pm25AutoScale = pm25Emissions.point_auto_scaling_factor();

const auto pm10UserScale = pm10.point_user_scaling_factor();
const auto pm25UserScale = pm25.point_user_scaling_factor();
const auto pm10UserScale = pm10Emissions.point_user_scaling_factor();
const auto pm25UserScale = pm25Emissions.point_user_scaling_factor();

zip_point_emissions(pm10.point_emissions(), pm25.point_emissions(), [diffThreshold, pm10AutoScale, pm25AutoScale, pm10UserScale, pm25UserScale](const EmissionEntry& pm10, const EmissionEntry& pm25) {
zip_point_emissions(pm10Emissions.point_emissions(), pm25Emissions.point_emissions(), [diffThreshold, pm10AutoScale, pm25AutoScale, pm10UserScale, pm25UserScale](const EmissionEntry& pm10, const EmissionEntry& pm25) {
const auto pm10AutoScaled = pm10.value().amount().value() * pm10AutoScale;
const auto pm25AutoScaled = pm25.value().amount().value() * pm25AutoScale;
if (pm10AutoScaled < pm25AutoScaled) {
Expand All @@ -222,20 +225,20 @@ static void validate_pm10_pm25(const EmissionInventoryEntry& pm10, const Emissio

{
// validate the diffuse emissions, do not take user scalings into account
const auto pm10AutoScale = pm10.diffuse_auto_scaling_factor();
const auto pm25AutoScale = pm25.diffuse_auto_scaling_factor();
const auto pm10AutoScale = pm10Emissions.diffuse_auto_scaling_factor();
const auto pm25AutoScale = pm25Emissions.diffuse_auto_scaling_factor();

const auto pm10UserScale = pm10.diffuse_user_scaling_factor();
const auto pm25UserScale = pm25.diffuse_user_scaling_factor();
const auto pm10UserScale = pm10Emissions.diffuse_user_scaling_factor();
const auto pm25UserScale = pm25Emissions.diffuse_user_scaling_factor();

const auto pm10Diffuse = pm10.diffuse_emissions() * pm10AutoScale;
const auto pm25Diffuse = pm25.diffuse_emissions() * pm25AutoScale;
const auto pm10Diffuse = pm10Emissions.diffuse_emissions() * pm10AutoScale;
const auto pm25Diffuse = pm25Emissions.diffuse_emissions() * pm25AutoScale;

bool autoScaledValueValid = true;
if (pm10Diffuse < pm25Diffuse) {
if (pm25Diffuse - pm10Diffuse > diffThreshold) {
autoScaledValueValid = false;
Log::warn("Invalid PM diffuse data for {} (PM10: {}, PM2.5 {})", pm10.id(), pm10Diffuse, pm25Diffuse);
Log::warn("Invalid PM diffuse data for {} (PM10: {}, PM2.5 {})", pm10Emissions.id(), pm10Diffuse, pm25Diffuse);
}
}

Expand All @@ -246,7 +249,7 @@ static void validate_pm10_pm25(const EmissionInventoryEntry& pm10, const Emissio
if (pm25UserScaled - pm10UserScaled > diffThreshold) {
if (autoScaledValueValid) {
// Invalid value is caused by the user scaling
throw RuntimeError("Invalid PM diffuse data after user scaling for {} (PM10: {} (auto scale = {} user scale = {}), PM2.5 {} (auto scale = {} user scale = {}))", pm10.id(), pm10Diffuse, pm10AutoScale, pm10UserScale, pm25Diffuse, pm25AutoScale, pm25UserScale);
throw RuntimeError("Invalid PM diffuse data after user scaling for {} (PM10: {} (auto scale = {} user scale = {}), PM2.5 {} (auto scale = {} user scale = {}))", pm10Emissions.id(), pm10Diffuse, pm10AutoScale, pm10UserScale, pm25Diffuse, pm25AutoScale, pm25UserScale);
}
}
}
Expand Down Expand Up @@ -339,7 +342,7 @@ static EmissionInventory create_emission_inventory_impl(const SingleEmissions& t
return total + current.value().amount().value_or(0.0);
});

if (diffuseEmission > 0 && pointEmissionSum > diffuseEmission) {
if (diffuseEmission > 0 && pointEmissionSum > (diffuseEmission * diffuseEmissionAutoScale)) {
// Check if the difference is caused by floating point rounding
auto scalingFactor = diffuseEmission / pointEmissionSum;
if (pointEmissionSum - diffuseEmission < 1e-4) {
Expand Down
12 changes: 12 additions & 0 deletions logic/include/emap/emissions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,18 @@ class EmissionEntry
_coordinate = std::optional<Coordinate>(coordinate);
}

EmissionEntry& with_source_id(std::string_view srcId) & noexcept
{
set_source_id(srcId);
return *this;
}

EmissionEntry&& with_source_id(std::string_view srcId) && noexcept
{
set_source_id(srcId);
return std::move(*this);
}

void set_source_id(std::string_view srcId) noexcept
{
_sourceId = srcId;
Expand Down
2 changes: 1 addition & 1 deletion logic/runsummary.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ void RunSummary::sources_to_spreadsheet(lxw_workbook* wb, const std::string& tab
ColumnInfo{"Point scaling user", 15.0},
ColumnInfo{"Point scaling auto", 15.0},
ColumnInfo{"Path", 125.0},
ColumnInfo{"Total emissions", 17.0},
ColumnInfo{"Diffuse emissions", 17.0},
ColumnInfo{"Emissions within grid", 17.0},
ColumnInfo{"Point Emissions", 17.0},
};
Expand Down
26 changes: 26 additions & 0 deletions logic/test/emissioninventorytest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,6 +264,32 @@ TEST_CASE("Emission inventory")
CHECK(emission.point_user_scaling_factor() == 2.0);
}
}

SUBCASE("PMCoarse calculation")
{
// empty NFR emissions
SingleEmissions totalEmissions(date::year(2019)), pointEmissions(date::year(2019));
ScalingFactors scalings;

// Two point emissions of which the sum is larger then the total reported emission (threshold = 90% -> 150 / 170 = 88%
pointEmissions.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM10), EmissionValue(120.0), Coordinate(10, 10)).with_source_id("id1"));
pointEmissions.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM10), EmissionValue(50.0), Coordinate(10, 20)).with_source_id("id2"));

pointEmissions.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM2_5), EmissionValue(100.0), Coordinate(10, 10)).with_source_id("id1"));
pointEmissions.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM2_5), EmissionValue(40.0), Coordinate(10, 20)).with_source_id("id2"));

totalEmissions.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM10), EmissionValue(300.0)));
totalEmissions.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM2_5), EmissionValue(150.0)));

SingleEmissions gnfrTotals(date::year(2019));
gnfrTotals.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::gnfr::Shipping), pollutants::PM10), EmissionValue(300.0)));
gnfrTotals.add_emission(EmissionEntry(EmissionIdentifier(countries::BEF, EmissionSector(sectors::gnfr::Shipping), pollutants::PM2_5), EmissionValue(100.0)));

const auto inv = create_emission_inventory(totalEmissions, gnfrTotals, {}, pointEmissions, scalings, cfg, summary);
checkEmission(inv, EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM10), 300.0 - 170.0 /* 130 */, 170.0);
checkEmission(inv, EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PM2_5), 150.0 - 140.0 /* 10 */, 140.0);
checkEmission(inv, EmissionIdentifier(countries::BEF, EmissionSector(sectors::nfr::Nfr1A3bi), pollutants::PMcoarse), 120.0, 30.0);
}
}

}
2 changes: 1 addition & 1 deletion vcpkg.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
"name": "gdal",
"features": [
"expat",
"spatialite",
"sqlite",
"geos"
]
},
Expand Down

0 comments on commit 54308ee

Please sign in to comment.