From 2531574be113f6a2169b55779b155202bfcaa3a7 Mon Sep 17 00:00:00 2001
From: Dirk Vanden Boer <dirk.vandenboer@vito.be>
Date: Thu, 29 Aug 2024 14:31:00 +0200
Subject: [PATCH] Use a faster csv parser to parse the point sources

---
 bootstrap.py                                  |   6 +-
 deps/geodynamix                               |   2 +-
 .../FindFastCppCsvParser.cmake                |  20 ---
 .../fast-cpp-csv-parser/portfile.cmake        |  15 --
 .../fast-cpp-csv-parser/vcpkg.json            |   7 -
 deps/vcpkg                                    |   2 +-
 logic/configurationparser.cpp                 |   1 +
 logic/debugtools.cpp                          |   1 +
 logic/emissioninventory.cpp                   |   1 +
 logic/emissionscollector.cpp                  |   4 +-
 logic/gridprocessing.cpp                      |   1 +
 logic/include/emap/emissioninventory.h        |  38 +++++
 logic/include/emap/inputconversion.h          |   1 +
 logic/inputparsers.cpp                        | 161 +++++++++++++-----
 logic/modelrun.cpp                            |   3 +-
 logic/outputbuilderfactory.cpp                |   3 +-
 logic/spatialpatterninventory.cpp             |   1 +
 .../05_model_parameters/code_conversions.xlsx | Bin 18253 -> 18851 bytes
 logic/test/data/point_sources.csv             |  19 +++
 logic/test/data/point_sources_empty_coord.csv |   3 +
 logic/test/gridprocessingtest.cpp             |   1 +
 logic/test/inputparsertest.cpp                |  74 ++++++--
 22 files changed, 253 insertions(+), 111 deletions(-)
 delete mode 100644 deps/overlay-ports/fast-cpp-csv-parser/FindFastCppCsvParser.cmake
 delete mode 100644 deps/overlay-ports/fast-cpp-csv-parser/portfile.cmake
 delete mode 100644 deps/overlay-ports/fast-cpp-csv-parser/vcpkg.json
 create mode 100644 logic/test/data/point_sources.csv
 create mode 100644 logic/test/data/point_sources_empty_coord.csv

diff --git a/bootstrap.py b/bootstrap.py
index 860474c..f692b3a 100755
--- a/bootstrap.py
+++ b/bootstrap.py
@@ -36,14 +36,10 @@
                     "c:/DEV/bld"  # avoid long path issues by using a short build path
                 )
 
-        overlay_ports = os.path.abspath(
-            os.path.join(os.path.dirname(__file__), "deps", "overlay-ports")
-        )
-
         if args.clean:
             vcpkg.clean(triplet=triplet)
         else:
-            vcpkg.bootstrap(ports_dir=os.path.join(".", "deps"), triplet=triplet, build_root=build_root, install_root=install_root, overlay_ports=overlay_ports, clean_after_build=args.clean_after_build)
+            vcpkg.bootstrap(ports_dir=os.path.join(".", "deps"), triplet=triplet, build_root=build_root, install_root=install_root, clean_after_build=args.clean_after_build)
     except KeyboardInterrupt:
         print("\nInterrupted")
         sys.exit(-1)
diff --git a/deps/geodynamix b/deps/geodynamix
index 0ab6194..7366eff 160000
--- a/deps/geodynamix
+++ b/deps/geodynamix
@@ -1 +1 @@
-Subproject commit 0ab61942c9f0787541d24f7b165f6cf37eb74b85
+Subproject commit 7366eff9d2b032225d36e5380d6881ef7704080e
diff --git a/deps/overlay-ports/fast-cpp-csv-parser/FindFastCppCsvParser.cmake b/deps/overlay-ports/fast-cpp-csv-parser/FindFastCppCsvParser.cmake
deleted file mode 100644
index 7ba8809..0000000
--- a/deps/overlay-ports/fast-cpp-csv-parser/FindFastCppCsvParser.cmake
+++ /dev/null
@@ -1,20 +0,0 @@
-include(FindPackageHandleStandardArgs)
-
-find_path(FastCppCsvParser_INCLUDE_DIR
-    NAMES csv.h
-    HINTS ${FastCppCsvParser_ROOT_DIR}/include ${FastCppCsvParser_INCLUDEDIR}
-)
-
-find_package_handle_standard_args(FastCppCsvParser
-    FOUND_VAR FastCppCsvParser_FOUND
-    REQUIRED_VARS FastCppCsvParser_INCLUDE_DIR
-)
-
-mark_as_advanced(FastCppCsvParser_ROOT_DIR FastCppCsvParser_INCLUDE_DIR)
-
-if(FastCppCsvParser_FOUND AND NOT TARGET FastCppCsvParser::csv)
-    add_library(FastCppCsvParser::csv INTERFACE IMPORTED)
-    set_target_properties(FastCppCsvParser::csv PROPERTIES
-        INTERFACE_INCLUDE_DIRECTORIES "${FastCppCsvParser_INCLUDE_DIR}"
-    )
-endif()
diff --git a/deps/overlay-ports/fast-cpp-csv-parser/portfile.cmake b/deps/overlay-ports/fast-cpp-csv-parser/portfile.cmake
deleted file mode 100644
index 25d840b..0000000
--- a/deps/overlay-ports/fast-cpp-csv-parser/portfile.cmake
+++ /dev/null
@@ -1,15 +0,0 @@
-# header-only library
-
-vcpkg_from_github(
-    OUT_SOURCE_PATH SOURCE_PATH
-    REPO ben-strasser/fast-cpp-csv-parser
-    REF 75600d0b77448e6c410893830df0aec1dbacf8e3
-    SHA512 aab418e98eb895dabd6369b186b7a55beddb84b89e358395a9f125829074916eff9086d80f9cd342d1bfd91acacc7103875c970a84164b75fff259cc93729285
-    HEAD_REF master
-)
-
-file(INSTALL ${SOURCE_PATH}/csv.h DESTINATION ${CURRENT_PACKAGES_DIR}/include)
-file(INSTALL ${CMAKE_CURRENT_LIST_DIR}/FindFastCppCsvParser.cmake DESTINATION ${CURRENT_PACKAGES_DIR}/share/cmake)
-
-# Handle copyright
-configure_file(${SOURCE_PATH}/LICENSE ${CURRENT_PACKAGES_DIR}/share/${PORT}/copyright COPYONLY)
diff --git a/deps/overlay-ports/fast-cpp-csv-parser/vcpkg.json b/deps/overlay-ports/fast-cpp-csv-parser/vcpkg.json
deleted file mode 100644
index 8e2640a..0000000
--- a/deps/overlay-ports/fast-cpp-csv-parser/vcpkg.json
+++ /dev/null
@@ -1,7 +0,0 @@
-{
-  "name": "fast-cpp-csv-parser",
-  "version-string": "2021-01-03",
-  "port-version": 1,
-  "description": "A small, easy-to-use and fast header-only library for reading comma separated value (CSV) files",
-  "homepage": "https://github.com/ben-strasser/fast-cpp-csv-parser"
-}
diff --git a/deps/vcpkg b/deps/vcpkg
index 53238eb..57adc34 160000
--- a/deps/vcpkg
+++ b/deps/vcpkg
@@ -1 +1 @@
-Subproject commit 53238eb7cdb557aa018585f4b12e288727fc9bdb
+Subproject commit 57adc34a5f058ef12ebcba10a78fec9026da9b11
diff --git a/logic/configurationparser.cpp b/logic/configurationparser.cpp
index d0ad5e8..4614f67 100644
--- a/logic/configurationparser.cpp
+++ b/logic/configurationparser.cpp
@@ -14,6 +14,7 @@ namespace emap {
 
 using namespace inf;
 using namespace std::string_view_literals;
+namespace gdal = inf::gdal;
 
 static EmissionDestination emission_destination_from_string(std::string_view str)
 {
diff --git a/logic/debugtools.cpp b/logic/debugtools.cpp
index 76ec229..dbb6a0d 100644
--- a/logic/debugtools.cpp
+++ b/logic/debugtools.cpp
@@ -18,6 +18,7 @@ namespace emap {
 
 using namespace inf;
 using namespace std::string_literals;
+namespace gdal = inf::gdal;
 
 class VectorBuilder
 {
diff --git a/logic/emissioninventory.cpp b/logic/emissioninventory.cpp
index 004bc33..0f9cdf2 100644
--- a/logic/emissioninventory.cpp
+++ b/logic/emissioninventory.cpp
@@ -15,6 +15,7 @@ namespace emap {
 
 using namespace inf;
 using namespace date::literals;
+namespace gdal = inf::gdal;
 
 static fs::path throw_if_not_exists(const fs::path& path)
 {
diff --git a/logic/emissionscollector.cpp b/logic/emissionscollector.cpp
index 7587547..ddf66e5 100644
--- a/logic/emissionscollector.cpp
+++ b/logic/emissionscollector.cpp
@@ -77,9 +77,7 @@ void EmissionsCollector::add_emissions(const CountryCellCoverage& countryInfo, c
     }
 
     for (auto& entry : pointEmissions) {
-        if (entry.value().amount() > 0.0) {
-            _outputBuilder->add_point_output_entry(entry);
-        }
+        _outputBuilder->add_point_output_entry(entry);
     }
 
     if (diffuseEmissions.empty() && !pointEmissions.empty()) {
diff --git a/logic/gridprocessing.cpp b/logic/gridprocessing.cpp
index 3bf8621..2ead359 100644
--- a/logic/gridprocessing.cpp
+++ b/logic/gridprocessing.cpp
@@ -32,6 +32,7 @@ namespace emap {
 
 using namespace inf;
 using namespace std::string_literals;
+namespace gdal = inf::gdal;
 
 gdal::VectorDataSet transform_vector(const fs::path& vectorPath, const GeoMetadata& destMeta)
 {
diff --git a/logic/include/emap/emissioninventory.h b/logic/include/emap/emissioninventory.h
index 690a13b..5e01f11 100644
--- a/logic/include/emap/emissioninventory.h
+++ b/logic/include/emap/emissioninventory.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include "emap/emissions.h"
+#include "infra/math.h"
 
 namespace emap {
 
@@ -223,6 +224,26 @@ class EmissionCollection
         throw inf::RuntimeError("No emission found with id: {}", id);
     }
 
+    const TEmission& emission_with_id_at_coordinate(const EmissionIdentifier& id, Coordinate coord) const
+    {
+        auto iter = std::find_if(_emissions.begin(), _emissions.end(), [&id, &coord](const TEmission& em) {
+            if (em.id() == id) {
+                if (auto coordOpt = em.coordinate(); coordOpt.has_value()) {
+                    return inf::math::approx_equal(coord.x, coordOpt->x, 1e-4) &&
+                           inf::math::approx_equal(coord.y, coordOpt->y, 1e-4);
+                }
+            }
+
+            return false;
+        });
+
+        if (iter == _emissions.end()) {
+            throw inf::RuntimeError("No emission found with id: {} at coordinate {}", id, coord);
+        }
+
+        return *iter;
+    }
+
     std::optional<TEmission> try_emission_with_id(const EmissionIdentifier& id) const noexcept
     {
         auto emissionIter = find_sorted(id);
@@ -243,6 +264,23 @@ class EmissionCollection
         return result;
     }
 
+    std::vector<TEmission> emissions_with_id_at_coordinate(const EmissionIdentifier& id, Coordinate coord) const
+    {
+        std::vector<TEmission> result;
+        std::copy_if(_emissions.begin(), _emissions.end(), std::back_inserter(result), [&id, &coord](const TEmission& em) {
+            if (em.id() == id) {
+                if (auto coordOpt = em.coordinate(); coordOpt.has_value()) {
+                    return inf::math::approx_equal(coord.x, coordOpt->x, 1e-4) &&
+                           inf::math::approx_equal(coord.y, coordOpt->y, 1e-4);
+                }
+            }
+
+            return false;
+        });
+
+        return result;
+    }
+
     size_t empty() const noexcept
     {
         return _emissions.empty();
diff --git a/logic/include/emap/inputconversion.h b/logic/include/emap/inputconversion.h
index f7b42db..3ee75d7 100644
--- a/logic/include/emap/inputconversion.h
+++ b/logic/include/emap/inputconversion.h
@@ -3,6 +3,7 @@
 #include "infra/exception.h"
 #include "infra/string.h"
 
+#include <algorithm>
 #include <string>
 #include <string_view>
 #include <utility>
diff --git a/logic/inputparsers.cpp b/logic/inputparsers.cpp
index 5aaccee..ae6eca5 100644
--- a/logic/inputparsers.cpp
+++ b/logic/inputparsers.cpp
@@ -18,6 +18,7 @@
 
 #include <cassert>
 #include <cpl_port.h>
+#include <csv.h>
 #include <exception>
 #include <limits>
 #include <string>
@@ -25,18 +26,10 @@
 #include <utility>
 #include <vector>
 
-#if defined(_MSC_VER)
-#pragma warning(push)
-#pragma warning(disable : 4267 4244)
-#endif
-#include <csv.h>
-#if defined(_MSC_VER)
-#pragma warning(pop)
-#endif
-
 namespace emap {
 
 using namespace inf;
+namespace gdal = inf::gdal;
 
 static int32_t required_csv_column(const inf::CsvReader& csv, const std::string& columnName)
 {
@@ -209,64 +202,131 @@ SingleEmissions parse_point_sources(const fs::path& emissionsCsv, const RunConfi
 
     size_t lineNr = 2;
 
-    std::unordered_map<PointSourceIdentifier, double> pointSourceEmissions;
-
     try {
         Log::debug("Parse emissions: {}", emissionsCsv);
 
-        /*
+        std::vector<EmissionEntry> pointSources;
+        std::unordered_map<PointSourceIdentifier, double> pointSourceEmissions;
+
+#if 1
         using namespace io;
-        CSVReader<21, trim_chars<' ', '\t'>, no_quote_escape<';'>, throw_on_overflow> in(str::from_u8(emissionsCsv.u8string()));
+        CSVReader<25, trim_chars<' ', '\t'>, no_quote_escape<';'>, throw_on_overflow> in(file::u8string(emissionsCsv));
+
+        in.read_header(ignore_missing_column | ignore_extra_column,
+                       "type",
+                       "scenario",
+                       "year",
+                       "reporting_country",
+                       "nfr_sector",
+                       "gnfr_sector",
+                       "pollutant",
+                       "emission",
+                       "unit",
+                       "x",
+                       "y",
+                       "hoogte_m",
+                       "diameter_m",
+                       "temperatuur_C",
+                       "warmteinhoud_MW",
+                       "debiet_Nm3/u",
+                       "Debiet_Nm3/u",
+                       "dv",
+                       "type_emissie",
+                       "EIL_nummer",
+                       "exploitatie_naam",
+                       "EIL_Emissiepunt_Jaar_Naam",
+                       "Activiteit_type",
+                       "subtype",
+                       "pointsource_index");
+
+        if (!in.has_column("nfr_sector") && !in.has_column("gnfr_sector")) {
+            throw RuntimeError("Missing nfr_sector or gnfr_sector column");
+        }
 
+        auto sectorType             = in.has_column("nfr_sector") ? EmissionSector::Type::Nfr : EmissionSector::Type::Gnfr;
+        bool hasSubTypeCol          = in.has_column("subtype");
+        bool hasDvCol               = in.has_column("dv");
+        bool hasPointSourceIndexCol = in.has_column("pointsource_index");
+        bool hasCoordinates         = in.has_column("x") && in.has_column("y");
+
+        int32_t year = 0, dv = 0;
+        double value, height, diameter, temp, warmthContents, flowRate;
+        char *type, *scenario, *countryStr, *nfrSectorName, *gnfrSectorName, *pollutantName, *unit, *emissionType, *eilNr, *explName, *eilYearName, *activityType;
+        char *subtype = nullptr, *psIndex = nullptr;
+        char *x = nullptr, *y = nullptr;
+        while (in.read_row(type, scenario, year, countryStr, nfrSectorName, gnfrSectorName, pollutantName, value, unit, x, y, height, diameter, temp, warmthContents, flowRate, flowRate, dv, emissionType, eilNr, explName, eilYearName, activityType, subtype, psIndex)) {
+            auto sectorName    = std::string_view(sectorType == EmissionSector::Type::Nfr ? nfrSectorName : gnfrSectorName);
+            const auto country = countryInv.try_country_from_string(countryStr);
 
-        int32_t year;
-        double value, height, diameter, temp, warmthContents, x, y, flowRate;
-        char *type, *scenario, *countryStr, *sectorName, *pollutantName, *unit, *excType, *eilNr, *explName, *naceCode, *eilYearName, *activityType, *subtype;
-        while (in.read_row(type, scenario, year, countryStr, sectorName, pollutant, value, x, y, unit, height, diameter, temp, warmthContents, flowRate, excType, eilNr, explName, naceCode, eilYearName, activityType, subtype)) {
-            if (sectorInv.is_ignored_sector(sectorType, sectorName) ||
-                pollutantInv.is_ignored_pollutant(pollutantName)) {
+            if (sectorName.empty() ||
+                sectorInv.is_ignored_sector(sectorType, sectorName, *country) ||
+                pollutantInv.is_ignored_pollutant(pollutantName, *country)) {
                 continue;
             }
 
-            double emissionValue = to_giga_gram(to_double(line.get_string(colEmission), lineNr), line.get_string(colUnit));
+            double emissionValue = to_giga_gram(value, unit);
+            if (emissionValue == 0.0) {
+                continue;
+            }
 
             auto sector    = sectorInv.try_sector_from_string(sectorType, sectorName);
-            auto country   = countryInv.try_country_from_string(line.get_string(colCountry));
             auto pollutant = pollutantInv.try_pollutant_from_string(pollutantName);
-            if (sector.has_value() && country.has_value() && pollutant.has_value()) {
-                EmissionEntry info(
-                    EmissionIdentifier(*country, *sector, *pollutant),
-                    EmissionValue(emissionValue));
-
-                info.set_height(line.get_double(colHeight).value_or(0.0));
-                info.set_diameter(line.get_double(colDiameter).value_or(0.0));
-                info.set_temperature(line.get_double(colTemperature).value_or(-9999.0));
-                info.set_warmth_contents(line.get_double(colWarmthContents).value_or(0.0));
-                info.set_flow_rate(line.get_double(colFlowRate).value_or(0.0));
-
-                std::string subType = "none";
-                if (colSubType.has_value()) {
-                    subType = line.get_string(*colSubType);
+
+            if (sector.has_value() && pollutant.has_value()) {
+                PointSourceIdentifier ps;
+                ps.sector         = *sector;
+                ps.country        = *country;
+                ps.pollutant      = *pollutant;
+                ps.height         = height;
+                ps.diameter       = diameter;
+                ps.temperature    = temp;
+                ps.warmthContents = warmthContents;
+                ps.flowRate       = flowRate;
+
+                if (hasSubTypeCol) {
+                    ps.subType = subtype ? subtype : "none";
+                } else if (hasPointSourceIndexCol) {
+                    ps.subType = psIndex ? psIndex : "none";
+                } else {
+                    ps.subType = "none";
                 }
 
-                info.set_source_id(fmt::format("{}_{}_{}_{}_{}_{}_{}_{}", info.height(), info.diameter(), info.temperature(), info.warmth_contents(), info.flow_rate(), line.get_string(colEilPoint), line.get_string(colEil), subType));
+                ps.eilNumber = eilNr;
+                ps.eilPoint  = eilYearName;
 
-                if (colX.has_value() && colY.has_value()) {
-                    auto x = line.get_double(*colX);
-                    auto y = line.get_double(*colY);
-                    if (x.has_value() && y.has_value()) {
-                        info.set_coordinate(Coordinate(*x, *y));
+                if (hasCoordinates) {
+                    auto xVal = str::to_double(x);
+                    auto yVal = str::to_double(y);
+                    if (xVal.has_value() && yVal.has_value()) {
+                        ps.coordinate = Coordinate(*xVal, *yVal);
                     } else {
-                        throw RuntimeError("Invalid coordinate in point sources: {}", line.get_string(*colX), line.get_string(*colY));
+                        throw RuntimeError("Invalid coordinate in point sources: x='{}' y='{}'", x, y);
                     }
                 }
 
-                result.add_emission(std::move(info));
+                if (hasDvCol) {
+                    ps.dv = dv;
+                }
+
+                if (combineIdentical) {
+                    pointSourceEmissions[ps] += emissionValue;
+                } else {
+                    pointSources.push_back(ps.to_emission_entry(emissionValue));
+                }
+            } else {
+                if (!pollutant.has_value()) {
+                    Log::warn("Unknown pollutant name: {}", pollutantName);
+                }
+
+                if (!sector.has_value()) {
+                    Log::warn("Unknown sector name: {}", sectorName);
+                }
             }
+
             ++lineNr;
-        }*/
+        }
+#else
 
-        SingleEmissions result(cfg.year());
         inf::CsvReader csv(emissionsCsv);
 
         auto colCountry   = required_csv_column(csv, "reporting_country");
@@ -304,6 +364,9 @@ SingleEmissions parse_point_sources(const fs::path& emissionsCsv, const RunConfi
             }
 
             double emissionValue = to_giga_gram(to_double(line.get_string(colEmission), lineNr), line.get_string(colUnit));
+            if (emissionValue == 0.0) {
+                continue;
+            }
 
             auto sector    = sectorInv.try_sector_from_string(sectorType, sectorName);
             auto pollutant = pollutantInv.try_pollutant_from_string(pollutantName);
@@ -339,7 +402,7 @@ SingleEmissions parse_point_sources(const fs::path& emissionsCsv, const RunConfi
                 if (combineIdentical) {
                     pointSourceEmissions[ps] += emissionValue;
                 } else {
-                    result.add_emission(ps.to_emission_entry(emissionValue));
+                    pointSources.push_back(ps.to_emission_entry(emissionValue));
                 }
             } else {
                 if (!pollutant.has_value()) {
@@ -354,12 +417,16 @@ SingleEmissions parse_point_sources(const fs::path& emissionsCsv, const RunConfi
             ++lineNr;
         }
 
+#endif
+
         if (combineIdentical) {
             for (const auto& [ps, emission] : pointSourceEmissions) {
-                result.add_emission(ps.to_emission_entry(emission));
+                pointSources.push_back(ps.to_emission_entry(emission));
             }
         }
 
+        SingleEmissions result(cfg.year());
+        result.set_emissions(std::move(pointSources));
         return result;
     } catch (const std::exception& e) {
         throw RuntimeError("Error parsing {} line {} ({})", emissionsCsv, lineNr, e.what());
diff --git a/logic/modelrun.cpp b/logic/modelrun.cpp
index cbaa9d8..eca7760 100644
--- a/logic/modelrun.cpp
+++ b/logic/modelrun.cpp
@@ -32,6 +32,7 @@
 namespace emap {
 
 using namespace inf;
+namespace gdal = inf::gdal;
 
 struct SpatialPatternProcessInfo
 {
@@ -453,4 +454,4 @@ int run_model(const RunConfiguration& cfg, const ModelProgress::Callback& progre
         return EXIT_FAILURE;
     }
 }
-}
\ No newline at end of file
+}
diff --git a/logic/outputbuilderfactory.cpp b/logic/outputbuilderfactory.cpp
index 385e15c..ea53d8e 100644
--- a/logic/outputbuilderfactory.cpp
+++ b/logic/outputbuilderfactory.cpp
@@ -13,6 +13,7 @@
 namespace emap {
 
 using namespace inf;
+namespace gdal = inf::gdal;
 
 static std::unordered_map<CountryId, int32_t> parse_chimere_country_mapping(const fs::path& mappingPath, const CountryInventory& countryInv)
 {
@@ -90,4 +91,4 @@ std::unique_ptr<IOutputBuilder> make_output_builder(const RunConfiguration& cfg)
 
     throw RuntimeError("No known output builder for the specified grid definition");
 }
-}
\ No newline at end of file
+}
diff --git a/logic/spatialpatterninventory.cpp b/logic/spatialpatterninventory.cpp
index 011241c..8172e7c 100644
--- a/logic/spatialpatterninventory.cpp
+++ b/logic/spatialpatterninventory.cpp
@@ -20,6 +20,7 @@
 namespace emap {
 
 using namespace inf;
+namespace gdal = inf::gdal;
 
 static std::set<date::year> scan_available_years(const fs::path& spatialPatternPath)
 {
diff --git a/logic/test/data/_input/05_model_parameters/code_conversions.xlsx b/logic/test/data/_input/05_model_parameters/code_conversions.xlsx
index 1f8f1c680a2636a16fa3ceb8a190a9889d9172dd..2a928b1017144660449c030efe3d5c069c71747f 100644
GIT binary patch
delta 9002
zcmb_?1ymf{wk_`N?(XivU4kTda0~95gbME3I0Schg1b8er-LOyf(A%%UXye0J?Gr<
z$Nle(_j**<SY2z*wdUScwQH~1JsAsK0D`W4KmcTvY=RSEp`hxJprA0JprCvlIK7?S
z?9H8>?KylL9V)fgoYsV}-kF#_N^Q2&Xz8kiYYo)LpwHk*WvH2;5k<@D%hGA8Xx0Vg
zK6c66;Iwm~IJbqWWqjE?otQ|U{aWJr!jMz-=c`6$8bTud%1OE2HLJ)6|Fjia%ge}S
z_$okdT<c=Y?_&R8T7pAYT3iD8%mG?*t)x_|T!pM;nsssJLE57-1f(?*%Az+`H1C}_
zn_olpLha}yQP^`C-zWj32r)#NB>L$U*)KBv22IcjmBlD~zwl3P$PGX_uuB@(Vjf-`
z85mUAK9dU-&~iE;M!|6%#haHwYxX#L_Q=2nmOoto=CM;?0jPo`TV*y%%LuBmf%@IH
zL@x%l0<~Vw(54vOM#5ppcd}#s(bJ&HzGpgD81kpi1=#T^v&t~Nw^RnTyWjFF#uJwg
zU1^P?yugI3d#Vc~<+e4&@Vb?Xns6$0<_e^4dqc1_Tc6Nr<EVD+ICn*D@hb~uswe=-
z5mBGgXs&3mq+uXewf&%_3jVlkg=Cmvx4*7Cf@lRC(|ViDBibt-Cbk*ED-_yE+?!7u
zn%eq7l`oR$&u`c4cSgwUsxz0&N6%YfOm^Y5J-Rj>1xtl72x6G@I`>M0v7j&g{i%%o
z445CKy7-p3U8$Hzw=*9l5Rk6MEg6AXQF`_(fcF~B_D0!+%9bj|+;Cnn9)FvC^%aBV
znf%*StiZZ=M+kLC>aM249IlQ$kqQ(Pan`PB?+|S))4l{^l&>C0$&`$QeaRsM&s`n|
zZ(42+ipB_b+%Nk*m@R*Jw}kP1>et9hIiMSBZ1CN=gMENSPM~*2*d#8M)%^;HAuLx;
z2(!RsacKBOxHFDI1mcSHSW(jDRhcDq1;4+7TcaDJS;Or8rI}>3)@vzWtWTs0j<6vT
zModzTgSWG8Hh4k!q4Ya6No$~gwP9@pZ@p!HO@d%ZSztfnsuo3>bo*Lnn+Qx;@wkQI
zNzWqhk7evJ+Yd4!`=jfWxQ9c)YGBpEu5mi?z+}^ocvq~;;hN6@e6vJkuj_{U(L~P2
zd8&-(#XFKV#@+7{9}5sONmRUQ_$TZ0wVo+vNc4}NW+z*$2EA0vBHMhU=tUMiv<w^h
z5Iy0LGq6k%Cu3qWib}&XTHJg1q0+1S&J>YM_?xZ#0u)$E5TW7uh0xkD&#~b6<d6=u
zg7$}Xb&vHH^~bSORv(iF-WZ!W8xm_1CrYOn0sn&@h8{n@dKWL-)*C`l0ux@t0E#pW
zN*j`Rs}MXC6d5r{ikukWe!5|^Kx?s?Tv(wete50$8$khBiMezYWF($N#xI{4i&`V;
za1}SZYc-vG<`zCw-j{2{*@V7*d>|I1kRaDq`)OMFX)<`oOct3>ny`LFb#I!o!c;NB
zvxC{EiZpEoSMP@J{>tp_jI}*0ltTQ>GICXxz57sQsiZ6E0Rcdg#oVO1zE!cpwMMcj
z$Ed}xlFwujIUbsg-0B;gH?D=EXIJ-P<|FTb<RtkwXV#TfR<#HPm&nYcYV=?ocACZd
zl;clQ%$hX7XTH@_&lpe2og`E@J|B9YU)2KOU!sikA+tKE9I`y3yKhY%!efXC1?S$~
zW14Lz5ihHB2h;$iU*0hf-jlNWhFpF|9(Uq;Y*+094F*5N+y-6Ow~Pe(cPIvWP+$g{
zHD<uOUf7DK$?5Y@%5!>%aL^ss)<hC9F?tZ}p7Rwrh#Zsbr3MGR!^It5+Q>uLf>JYA
z&T{IO+$YdjGS1T67WJySTvW+%)eL562KnG}WUcr6UUvY3J=Kr%*Ed9JE0dRwZzoL`
zJg;?SP?)Z(y*4hkct_N}awYSr*7<FS`??wg;xD(nruZqQvqwUFMW=jo9_Bb42lj~%
zso)yx^Ou+m4HBm}xA~yc&528RtsmaLfb80V1Y(OW(479hXbM)odG#=CYJ~9!8L=P#
z*d+mJV5c=Mv=9^eQ&=|YvWGtO+I&tNOK7hJbk7+_X9-(!I*L^?uj5q;7!!I>4%_tj
z()iaKx3=uq^R*Vr_LdT*V(^dqi?&ZW6rQX1Np2z%{@0go`ctV(-{(E9EazRtTjr*l
znjfsny1OaI%|3Ti?%rQJU!mn)xR?l&>p%mK=oL-m1o~^lxzUmftFa`H*E#$W&NlB^
z;K~=hl0wPXr=UD#a_t-yz;H>^PAxeL>;dIX<erooYs8fYu}>5~q|zfajb9LaF;mXI
zStYS4CS2uWEJd^$($;Iz_oD&nUkt=t=wsayY2tXoT7m|glIY|pXfe-1<7c-3O5P7e
zBNk{hxe@IjbCe>*<Y~8!s9nF*<M~Fb44}$vuUj$yAj%7G39XLJfY2=P&kkmU{PiI~
zc9@57M5z>-xVk8xpwhtyRjp0tom3sMsn)$zHaUtN!b^nNJwjD_>{2afMq@5}s=`s|
zifB`Vv*j5kB3015B5`lyHqaXd8&?+X-^VyYB;RS~0e?}3PDvHV9qkvUm*u0VEe}f@
zuNmz(SQ!v+O+LmX-)YVkD^|R6Lo`^L8eS!jHD^f^t*Od-Om26EG@<3=Kg3jwk7?{M
z0})_BbdFJnkBPSPxfDwg94N`X70yM;aP9<=;Oh4*z-0tuX=oW3?17BW;;Qn!L@kZm
zZyndDt7EXEd_4zRgEaWDqqxm&2lqRXq?KeHqje2j{dxQNY$y7Bitu+k*g}SNiE1D%
zINcl^0$XB@gcrTnGw(Oem@qc$n2S8ug)|*jYeTQYK!;H?K#<hiyN}%JP^*Wdf;{At
zQ8v~HQgIJF6W}i&IAfG8?sI#Fek%Cqg{c}Vg%K!;CRa0e0R@9oz#Ypf$f2q#)XM}<
zlJH!4Hh5k%E%)%o7+3|6x{9NgGTJE&$|ua)R%w`)CdD#?Yrs_N1FG4){9jb&BHV(N
z2doui!}gt|1mKXw;|atP`yD!Eu&fRyNF`NvxxzF21cDiEG`I7g!F#2ZF!x1$y6WO<
z+nvkH9!J~3J#R^=Hm(Y5M>Ar}+Uo1$?984Um7v&Szu6<D5%mY<ap9vVDyHarqd-9w
zs6CCTAZHu|0A%Ex$%Pks!gwU1*M(?7r6fhI{hFl8`$bdLQY?2bN^^o^hV1;qEk>ad
zqufEhau0X&#_roIAAdlD-4B0JsQ^{AABi@vwXE<n!PRb@r%PDExvFj@zX9GA6IqwZ
z+k2i3(suV4EDC6;b{N6W-on(eg3{cC8P8wqmBnaJ0R1&!oF)C<a?L^m9EI5=CI)_S
zY_~V4{hyk?m|OYO1pj0V$7ess>bAGso#GpG?q#<2en;q1AR_5t$$IY^Y+mT|(6}&E
zuZHSN`Fc!9M9;_E2h@=1`Y;AcTqm)yql8J%GVgjWK}&3QC*3yn$bhC6fnVb$)Z_i~
zt~Gw|3(yz2VK&NWHjnDSlsfrFAOGr?UrJHZoK-EG?suvfTf8KzYoZx9f=P<8Fw3@l
z;MCk$Nob3j)l$xAlYWxYmS@nCn^JTy74OZjX(Dv=R&HxSHnD($*UQS)KdXmB@0>L0
z?G{@eRt~$7TOGwb0Fq*?J(oFO&sWI9d4o`R3;=9oCc#}VA+zu+{5>o$q9bFSt^iGf
zIVKfr+-RSZ4oyzybLRsP#(~H|PqwF$P<r8f`y*0bXaUxtewzPgeDv7^<Qe(TQ8Meu
zX-ytW?rD(ZBh$KL$%48p#5DhiaKp;JP35@J<?SewmEOeZ{zxyZvK5>X%VY-(x2OZ+
z;{K0rkM|GfM=jkC*QZDGG2pj8H@6SpaUU+ezdheFX<wJPyBP*dC~_z$y8RzLc<*??
zZ@b>TRss4A?#6K+9(EeoZ;ys4bDSl_JLcMDm#AOFdGNK*#hrc}iw;CvG`8Zfa5bA5
zPJNUV|MJUcM5xjk^?4xGB1ItZ0^4$Cjd+;Y-*AhD$@D1Cg&76yCWLv@_s9G+apg0-
zZ}~mB)2e3Jc(p$@^dez`y0}8LHF{QMwtSqPbypq}T_iI*1i~y*pco3HqTLMPc}YRr
zV}@j4BD}Lf5nw{^->w*V#p1a{abR?i`z6FyOW6W*cc<v#NAzn<&)0YdK(4AxMCstc
zWK({U2-CcVmMwRfyCZ!O#u#7HD+AT-S8tTcf_TvB(+L`cmaMf7An?)RLB6S`%qMvG
z{y<!|iN}E0RMUr{(jcCgwCNWyX(O#a=ZQQuim`G1%CUKcy;Du2+h!wU>+)?#9MAYQ
zRH~SYFqEf;h*%+RK(zH!WOJAZoxPk?x_sqc8RR=>6vcBuTebHaz}oiPKP8`sF*fxN
zfEKF=4@3#A&c#bw9gi}6a^k0KPvM@>vC&TrVG$X6&hd{%{k!c&Y02!Vk<ZPH?+8TI
zaE{(YT*qo5T9GKrl`M}o=_LKX?c4RK(jC3xyq6x4A70tZK_dJ|8{g<p%-`+rO-~)b
z$Du5@whh?JwfM0sKlO1@eyZMV)jzXy*5w;Mbrj@0L>ikjwiswClPB3S;dc&1I)TYG
z<jIHsGehogjDL27Uj#nb?n!z+H!HqlR1H1={+G&F{QrjiVg0Z8_LD>+PqJ}{#HqUc
zH~*CT|Ff0(2#WF@&zkZDK1Bm=hdrvY`+JF9F_d<S!d|8f`AHrT^i6riNfwCKG|DQe
z%UYy%s=pndJR<m;^1l7yP*2rUCAHJ|?ZESlkZa1jjm$P5o>i5#J;2GdGs7JjU#ce)
z<-*l0ZyP9I)o2*2gXP68s!%yBN9V7EII#7T>zfAw1^ud&<lm|EmtHQ`g7!Um5eH)n
z(!;eD`I8sy^#73Lsb<w=pc6};7Qtj?^y04+oeHb3%K^XPMe~1!JE7wJhdAL5%KB^D
z1(kpWJL@pUgltm16Y~d)PV0ZfVo&z0{Ab+%v0dc;SMFhRn%1FE21GbmOMm!xBf<~=
zNdAJI75F3haakqg$N#a-{NZmgTS9Um=SLRKTcuy%p~2B}x+1RG7`TbKE9*?QV3pVy
zmtgX3Y8DysSi2SuL)(~@V2T)gwn#Vcc(6;x^VZ}s*OK6uK1j8tDaU=}6nr0)Xs~Xt
z#sfwGrl|XWU12Iaxpwo9daT_&k@qE-MhBlQLC8O!$lVY!;V&Zy`Rfze=rJhKc>N88
zJOun^KVT~RyFNa#A0cc`7m`8|c9)FgdjF(IcUA5WsxOq_pX0vxIV7pg;_)?p3*HH&
zv6L`;&I>|WkJBz{WvL+J>F9B!P>xrEKMt#y^GMdlAU56<z^5N6g!x6QV>>j)@)R~R
z2FFq>;8v=F&(<ngCY11_nj(PAQ`U3;Q9GlhevG;wb^Ti*T>8sR|2-BpBmLi*=CuM<
zSzILilC5hh-9}DfELfBiAXyM=R}08>1D~Skk@VH`*;s*z7LkXc7&&YCBCTv!$J)#)
zaWT?CWjbX&Y2_#qp-yz;LJ$qLd?{A8KgQa6A>Jp6fac|SO)(!}1@=jC{e49DNwl|&
zK}vDs?}RP;EK&^kA4wPL`6~WUJ*oF38r+jo6v6&S%jbT#e5jdi((w<~=Q@F^Yz~jb
z_jjh#-HTr`w(5@i17}r2<EDr9m*Ug-1}{5q(iTs|MS*7$di)30HzvESpZOE3L+2-S
z9)E&eZbLZp5-h;!A+0y5RNeM;@WL_9mnU23-twOE3&tfKPuQsWEZs->&n`AEkBK(P
zIt#CjXqhK%ODh0D=a0NBE$F5fbGzO#h2kaG)er2QTT`WCk7gG*`;bG5C7LyLx3SAY
zVH3bw!h_G9dEq1Uzh)7t(S|5}kST2oPCe8b4NX$>qdyW@`k;_khLSlwe9apV;l6#>
z<s(=s`=n6(7+EK~YWT)OZe~wrg<*Kf7>j`$RHb3+R!J#BJ3*t9OCnNW`!(Scj@Xyh
z-5P{Pw%4RaaenUT-`Z}z5*z2Mz@bbCbbzp_@sVkTx?8_N-qm@4B+01&kV+h+&>ucY
zzJV0_y+lY@c=afw)om^12RnyExRs`=j=M$W8=bJekUF2*068sS_Rc!=`pr*|m+Y&k
z7DdM>s`pT|v#v&?PA@hoH%vOvXus>}*6DcX_b{5B&EK7utVm^6IR?8RA$?$X*?E<N
zYwk;n*U)t=hnQeMI!Fd|WT4k|Vpo!YUaIs_DTbTZ_uEK2QJC$urM#--+hCq;<sAr5
z5#HMnsBh8D*pA!8<6Sd4ajo)`th8bU+y0)Gw(U2XeJA)PK9}a|!RHd<qOioX^Q&E1
z*{z_e`U9q=+wryk9dkO)kaOYfHD%vS8c|(pKg@>4*04c72q0yB8ze^=OfR^6`c&v^
z2RzDt4U}K;*p~w(`f;n}bn@tny{TAmYeLEwdfsr#nLW*dOiEEKE&mDLc&4L}nc;ke
zDK|%p>`;^v-S3rsCpk=w1|w83e${3W4;E|#p~`GCJO^(4#R)bF2h?cH)#vUx6M1-m
zzZ8dmt?4)u5Pys9HZLT^QAgQt8N~PFpb;k}h#{|$5?!y_(6k4aSvZS59#1s(Cl%|`
zd%|m1)g3BiMxwJ1o(A`AvK-gDs(OQI$UhT_i~UpRnJ4EStrP1mB>t?Ll9Fegg1FG+
z0Qj%42nefm+L4Xr5`jc>*==RHPPezh>_GWZbnbQ4QRg;a!aJn-TcN34|KoO-`zw}e
z1=!*=@!Rp>9E>LBuDqkL{<Ue}tk5$aO0~<mG!C@7n@FQOtG<&Kj8~*%hJ<5gO6ok|
zXY)8p>SQzzG5s9<N4Qyv?3Ftt@EGlAHW1sfGu9qo=#VE%#f7G6QpBmG$IKfR7@j97
z_gps)fDo>~c+fQQDVo2mgb4$qlc|4)yaVeIhR}o0J&b<w3rBpiCg^2#S+)`)xkP+k
zn6fCFg~2kJZIKRbyG+@0#TLo#7)e(?LGxwH_A_CAI6He3o!u9hKFWu<i9%x@URaXu
z0iJg4jHs{Ca}w#~a53iQ>Da)Fps{5uJQR}g0Up|7oaDN3Jv_er2rGHjXDk+}*vzkA
z;f@=&kmPs*HT7B|67-$nz+Z0Q|9wqr_1`K=|1WA%P?-!?g8+303dX}8N^3G46x0DO
z6cip5NJ^Ffc$;SEye)(`yrGbBVtOW9)chXy$PPZ*Ej{HqXKvn`*)G{oqEj3k2G`f3
zIsWRVHtjziEh)}L2`O{j50@$Y{lG=rxlwapR!SGz51WIhQhd$pFH0}i&L;;s*)<K4
z=~*UupV@Wp*I(M+o%PqURZ*=~k~Ik|9s@@nV_)S@fS11MC!X?&+Qx@TMbqYB=3gN?
zO0Bc2$8)iVUNs`3R`n@*-?Xh2s{B)js(8+;=C&df_nU_L=c2>B$exFhi&A(U)ZqGN
z_`Ih8azZ)YPE~1@`OvIn@k$E4gYVq0cwM0-ruLPXQ^j}dux`l{-8$}vTThSsd-v$z
z_^YeDSiov(KucvU(tIm;`aE6K4r(QX{nZC9*HOzOoNl-Z7)AJWUXFD8)VpM<;y$Lf
zw2)bJ#=UR@X=G*92fO-uS{?1YdYp0sCdDCXn~n2k5!~ArX&Wr1+V|M%1W~+o$`x6m
zb)-{HYo;yttze^+JvMJf-~4>AXzK3gK_1&kFW{ScF2kj=qTES7X8QRFe2)^5X(+5q
z+w^>6{$%e^qlM($=c4pgTdw2lXJ0X-(fs0<aW3gBScwfsSEO}+y+z)_HRRr@#zZhk
zeQ7`8qMqIyfN*#{KnjOCH6D1WmZU3gy8byZtP<8Mk^-c}nWtZSuChGCCz$o^<!z)U
z{Rp7RVf9YH>Ym=F+p-hKDhW|_1HJYHrnm+6k8O0KGgMD1pUiUGY*C8Jj9jsMw>kEL
zp62RSnG#n94Tc)Qk?x+Ah`L7zCYlzqyGj^)?+Q#tFNe&H<t_~}2NM~F*M$%nBHdP^
zaj#?bqt|dZo@PWGl7L8!Es|+5fqOZMn<$sb@sXXo-eV=d+BY~s^RRX3pJ;xVejPQO
z$uK2n&_fhs9usxBgdcepkC^H+LBoh(dr|BD?C=%mo{VO=RDHxBTNQPb!m^|<+8S!>
z`#|4>Km)9PzyqCI0xr*436}tiTTjY(;L4)IjDl9xB3M88JQESj3!H+dh&F%n_b?Cz
zUxBDXfd<;(A6ZOsd3H@an8b*+;Rx7C-dNlQ8g;!kg8aJil_i-qhi7%)4tr<gY!+zv
z$MqlbQI`O$dFECY$FpzBeLl)I7v@gMMeWrrpZne-cH#~SF0AUgzyEecHhq4P%ZrM*
zl6hrcow$2bUN?>U9f*BHIxT43U4q)gnm+K*OIF}Y?|tUc#4=LBdWG|tr;OSQ9)4)!
z!o>5kJ^#Ac<R=svH}5SMS5hfCqo&>J$bp&!!?GDT0zbo1d$%db6|gQf^$O~ZsKZfQ
zbp+wZ;Hgv8d1Gyx%1y>Z>)E_fqlVmBV)89v=Z3y7A$6;I+yD-z2-waqP^r;pG-@ES
z$rvwr#lWn7VDIpXqj5crZ^47|8A}SC!FQQthVAxiQQRZ7?3&Sq=ij9dAKm@EI0(X2
zG!<qxZ^B3g-&{!I^r!E=)9%~8-6;RI=nk%W4236!=_HUX5rh>Cr}Tc4?VR_8w}WGz
z7R|Y8W;Ki=_L#63X!|yP7k~88ddDgILhrJAHLs=ZOrjx=O92+eP6`AZz(7HnGJ_n1
znIH=k=C)3U24Fi}DZB@3(A=VA5!017?rCu@zN5wLI2L{sH@aG5E7PZBqp?UOx;%kS
ztT$UPF0k=(dJ#|YJJ^>?!gnW*9dH{m<;`6zrh|#SY^D#huBWD_gwq`zcIX#12NA}W
z3oCGpHU{m~CH5Y~QGDJr`b6W44e%_+mjL3DMS1?4^K*mFwRjk=Nu!g(8fk-$MW&jW
zG%QplYomO8nhyJH(8+Ff){39kul(;KZ$3WCUC&Mp=`fQSa^Tqy^6I@zQ*b(Faz}ko
z-<U}=C%9@u`xu#X(^;(EODz`ax*}MXFWtjZDOk>ayJ^g(Ag)h&rp(5z6!v*Vp#vaM
z$s(bPB{cPm|EO0w)N#vL*6-{f)<XK`=Ni)upJ(U*?S8EOrUZk)JZpWhWm4g30+!IW
zBg;M-3V<vvS$W`$64Xtb|2gnIl}i_P?Fjm*{M5#?1Qk7@es)=hs<4<%8)l8o{0JY=
zXe=X9o+$ZSs5y^O?S{E~q^#3sJwaf*tXNG3KRcgvv@t<!{ISvQ2z7w$)op-7AHmeY
zN<nvqg$bzNDW}leSFG4Oz~J5{FqwE9xytzve|I(zku-t(v!e5fe3X;wlwNL1Y$LNp
zuapS;2OE)?RiZZ3f@z+;u%wU-?S`?kiY#ViG4$zzljuDmPjk~8f&L5{e^{W1Fc=<M
zN)n3;={0AEOy<<v@u(>ogrtyf#^dgWD~aut5ylG@?T-E#V28z{At<%4xesl^UuZvm
z<@X{lcT9+(I;EuH0=2qb`(3oj^L6TJ3|tJD!v#<yxp;U>a}TLSN{hAmss~K(M#wix
z-IboI_70@%8p%dJ@Ox{Ny#DAYbsZl0x{`{^3=St*XGVb*WB8|tV@JWOIq{5_OM8~1
z<RAAm_`hQE2;>8SJ$gvB?}*<Os2>bO2?Lf>1g>3)-_O``n%jcm+2&4EUjHPYoE2;O
z@fBM;2Lq1GoL4bLJ$+<(>dtM_xDYqAhx6%eGTfeYj7`H1pCmM}1?Quq-#4jyc~kqM
zrWy>bS2JR%J^rz1^>*i2&X5?9Sh){nn*p8LsR7SgJFBnMDG<pjlZW|SaI;1oSCds7
z*|y`k+8hCu8*wbiEs`=Otijd8D?7JbPBXKocyt%f^=J=2P;+_FH*HZfeLLzr>N}%g
z-bGE8f|ti6Mk1G<Ixm^z%U$*XSF(OP&@R{Wafb2li|;lprnLol*Socxro<@@1Lofd
zm%8&(>;{M;NX%!e^&xXre`KE@$`+fp5+4#z$edY6saR12OL%m7`g!1YBI!TE{a%g4
z5cGjZ-f_%Jeg&Cgxi!Rz(Lw|EL0ou)8;G|UiQX}F;JrYDJ>?u7qKSJmXX;i*8MP$3
z4@WH~pBD86ex_x)1H!wB9{w+_PT_~flx@Ogxlts~l@&WDXYA`xY?x|eZF|R*WkR2u
z&Wd$+D*3AUT5Ybp<R@b6MaV8>|LPg)HMM|eG47q47Skc~#`Oj9V3q`s!Cp<g(ZIe*
z7kg3%GYsdB*3wf;5xy6V_f~Qd@r4&uFLzR|q-d0q>}lLmfacHXfak=*$sQ%E-Xg=~
z4e%#}z3&>`xCGR2`9U0^)9R<gH~U2uB+B69NZ7F#z086=udIFu7JG-yNfTG5GvAkM
ziD|8=rA)LnncT<+>Sg`VF5<iGQAeKX&4&SPZ%{l+E|mcAvYWVbUDI;DIbs}9_+ua;
zYHFpvV7xD}&ux-aCDxF1&mFfPq(21RsY7P#Y^oHIKJq)rYTMXnFeB~jOWkyj2!Hgp
zge)xbzy7xYGL(J}T>~nYW+(e|Js1j#2vSP?`8h)a0W$a?BpH0NKNmRu9b|&@8{~`w
z55kwlBK=pjiT3Y6G?2fH95fSXKt>X}1OzWjMf&%N>7RL_pe{g?vQjYZT%ga~q@WEh
z6p*qU7U-ia3CsdF<PiPupU83~FmpV=5t4EwWPd(2pr9E48Vl-(7vv&L2g;DcCwscH
z?jcuO6{H@<`x}T%5Hv0)38NwmLYJ3>$q|K6m_a`BEHEizph9^DvVUEtkOvd?-<bX4
Vpl|ZeV05HF^a?Zx6tYjw{{l|R?WO<#

delta 8186
zcmZ`;WmHt%8XizUx+J6#1f;uDq`MpG9J&QLNH@sPjVRq6(lUSoN+TnJGz=2beFwhl
zx2`|#th3JS^S;mXym8K%*=z4Y2+~MAQq46gV9tB~EaM&s)b#)a!UcgqzF-a?7k5Vs
z7Z*o%Ung*dj+0B4@WX5I&1<xeJ|x4A!I`R(k-GL7Wvp=Cr<yB4=(<QAarT-z7rP=1
z94d46pME(f&(^;3z0CH=^=P?nsZ{*+COUQh<-27Ttr@$V-?3BhH`h+$Quj<8pOo@w
z0>3lk4;kmSs9mnlTswX)nfITCpn|+a#xb8}b#u-&<P+If`${_s4z-p<hQ8M)ccNn&
znhVx6NC~&j%+<ixg^Pqr8VC0rWWV<uo~U};`2F<*Ljv?eSTNr&q)J7O*X5DE8)V79
z%a9_8STX{(q`WE>LttJbC<ZBu$HeV(x<L0#2sOWX4URz1SXMc6D#nlB#->OOwRBfd
z08B=j7iM<BkD5j+6rP`;yt+52=(}sRN?H4XgcE$tM}}&DWQb1VyZhWD;XCWqny~GY
zjs*ohI+cSkZDy_aRVnXM6rl7YvTe9Zodve!Zyn#pr?@(NoBi#zyvDRHQee%?XaYQ4
zr&+ltsvW+L&8BA%%-O4{sX?e<X{z+~=a1^HDdPl65{idaBO6_mgby(=4^cXM<;R06
z=0qN2RuZnt-{7-$PARk7PPlSm?-GB0`o+#*{N)>(M*jKLigH+g(v%A0RGlX3=kN!S
z4}HI~tMV%UARY$PxfFElup&p|Z+30~xZrE@$AU~psu<YU3-J=tO`S<{gOp74u(V?~
zdfS=+?>#fSu-cesFvUd6aP;V<f{3@T!`$rT3M3u?!6dzynySJE+ei7Fa&61kYj%oF
z5rkk=%Fm`d0RIW~yEVPOpWDU1wY@t1Wd(l4xKC=r0&TsZ0(6w-m>>LvUpSuweY>m>
z6nPgrVY52e#q)1%An5GLB!j{e|92jB!IZ1P?|tJBZNUcY-ZfD*)E9&MHDRcM6pwvG
zdla@?qxhQq1EJ_Woy*8$fzw-^ie#p#lnq|!^Co;&M9wV=_qWL0vdmXkHt%_{t11_q
zQL0qV6u~`=E<R50=(y7)Cvk@XFSWFQ8@#M(b%U>!m)}1HqL%Wk3cmXx5P`!@2$R8y
zM8>9oJ*(wuBA8KwDy%S5miTDw0iG+JOTm3gs<o_<i5^npXuf=-c2sYOAHoKPK%K=_
zbUt_){OX<zzQGityc!Sv;y)$*I`^eu+*t{$aF+mVB`8R*CXxqnUWPdc5nc(aFI#o_
ze;@WC_eaSl^5G;eZ`sYehWymq>wHH-rnbYu8{4J4K4D7U<z+@<#&=xDiDTb-h5{*I
zdRX6q@eCQqR8j5_4+R7oBn5#;L68I)B7p1m-OV;|a#<3@yAB+`ieOrxttHHt=`Kqy
zuNu#3Xq)f2*hamlOPKAU{KDJo<_cP%B8sZ#`+huB6yhg(JiKg`{_17vr;Vs`*}*GO
zhCbf-)6PW4+J(_&mlwOJdTx_;O)u&-cg{ouS5+?U_u(A%wxv{Al<UICxED55z-~iT
zvBLO65z|^5CdW=RD>|Ki-T)ogyHX-s_;FF0vR!68Bs%I<z#*qj+G%1|9Mh>HaYXd)
z(OA^dlL^ji8^*1*c_SU;r7DZ9eRW(6gQ>+OspMx}O{&W!%ta@#4(;uN(3ThSqtP2O
zJ=9@p$Ep*i66R&k3#{k_UOvc$0L$rWeO;<mgu0%&s>@x2*h7T<d7%c~KZvecY2Rq9
z#~hE=<$3b&WM^z@PK4-xS*Et)HK?kC?VCM!-lTbvG_=Uv(f%m$uz7WbfpTJ1=lsQo
z^h%V9nkZ(?0G<SwLC=hYB`YpR!VQ<6JaV|`c-?h~%OM+wlBH-o9i#wP6A&kx@HbQK
zvF}qdO*f7rrN})U3rI1MPQ8KsPWPZu5GKC3%-OFR+4{n4cz`XEe9m!pDMA9v2`3<v
zs$eiEg3l1SU<_*JQ;mCXgt7LZTqZD^f5NVd?PwUVt(_JTB*&G5@Y0!45k<&;-pUi`
z#wu&aCob;5k9QFm-%niO1pHhV?ztT1a7qiU(+R~FV7v758Vx50#M7K~x0Pm;T5<Em
z0Q{(=A+;SymNa%|vU$XQ!b-VjsnG(iQIg}Y(eb>%-L699q&{HwZ+=*AQpYBvC{Rj0
ziDUk)RmAlwt+Cm6Hw{~YFP4#*)91fix?7gXgpB$5#k)+^)zUnG*KPw|tFAc%>c$W}
zcFvuaD$1$@(j{&ftlt0^YRXPtXiv>&!NcvrMNWE>u~GCG`Qp12_&)zZP!cLUm&3<#
zF2tgj6T8+cxJ5`ol!*Z=3E#pIWs~qq|GW4@LR5QnsbDU*;z<OEMvAt{0k8Qm00Gh@
zs)rpZEKld*0{x*<;}-nKraWriRJhUTA%|hJkDtwR9~Pu@%3D3qB=9lXyB*p_RzBhN
zsGkgZ{EX}!FSPB2w4VK)5}AV|x>pKn;YYe6-tDOV)J!ee(3jrEWb}N{wmdC82hOc;
z*>;JOr10FYl>xCy;-!1csP=JwO;<r<vh}8<0H=B4w7QWBck>~>qCy@tol72E!p2`5
z``~ERK5<sv2+ZAl%(WH6d243;rttLfq!y_`;w&kg-I3cz>9D{KyPl?&Yg(+pi1e;t
z6`OJe?`M$>(3iu#vF5w-^5AT2x_T+KpSDt(mK7bl><_lRaOFte@jM!TBp^%<9BMx&
ze#Z@+qn`J}lOv($`X<Mr{D7Vce<b)JRR7|ICFj=XRGKiZj%?fqVoaH1qVndkjd^$l
zlP}@ju~GyOx+pi9hlROGgCX<?cpf{aV!)l5NC+JtJ@M#m2qRy`qYx@PRZ@J#mnr?R
zQqS)FQg+mxzs5x_HHsOy23Bfc>(#T}KuRlT1vJpb)F)g*ne_(3t^9H3$BxVFaTT$h
z`cGU$HU>R5Hr`98kH;mSUtu<o$1QLeQ*V{n@aety=wla3c%V*`64Yq}^UND))HdWQ
zShho(=~4oH(qUV3OMPbD?Sg6e@R*2?@0FQX7ei5TsS0q#qtzCLLF&gjqKRzMzvD1z
zFww*Y%aR_by}-8Emw?dH65;s9xjQyufIz~k|NXrr1G2M?)Ry&$Wf$Z!Z}3(rmn;Va
z?l(lMq!dxfI*c3%-M^>tgiG<632)}jj|GL2)9hE#9NVeqn*DV-M9USVkYAq{I;&F)
ztLJUnH)ixKeRAXE=3o9;;FwO(TLy0669zm~PCoTVSKx2e?dAIk&LY>p&Z~je%W0QR
za9*3U0MJW34J+TDa6Sp($^fsc%{w%S?dG-f4N7+zjCh~I*WBt}{qU^6xKb=%rT<dr
zP?yOzz1&C4&cFLj_GYgp-KY8L(t4`WCn65qwhBhO`8BcRcj4h!e&)~Loaxm0K{RjV
zVv2{~ty!PlTjH3UGC}ts=1T%Rkf+^f$q(?E0KNqrz6t@{FO4n((ywM9a)xYcS%EN#
z?a|uRp&p@E=fB?N46hBaz5jA!lI1orRCLj{x6o~_b(W{)zUpRlR#UsA=3F7{JiJ>q
z(*Jqs+TpTG%}c~<^p)U+tH;GAIJ13wdu?rPIitO_-q877iSzu;z)ac+W<bPMJ(*Km
z1W>+@=CeOWH0A65s#PNRx_+b6DMTW{yS^Rp{_NZd-(zwX{f)L^vH<L3213)gySC^J
zPA0BFSFd>5cByPQKF>>>y8v+F1T_)xyU?ba&KKWrT7EToP`J(x`<b4zja)pTqPn51
zB)(1wxo2~HwbvbUF;93sncHr!YrZ*l30Swo^=@8Wox|rMj%u<_;6HB`O#R$Urwx9p
zRz#Y#y}-Y@X<bSdU0hjg2hN-8=eN(@Cjj!Fi6N)Qerqqfca}!hV;Hy6u7bj9F`G_z
z0y?~P{Dsb<O@x8>0{E5n>5B3m>+zn*Zwj%!`!1Yn_Qq^7SstN2kq`~>O-cC44EX8Z
z%;8;JT!riVa<V4VhG;VN*>R|_;pqFirZc0VT7RN&t@r>)nEmW2(8j@19HL;3xkpSz
z!Mxg?uRUDQB#3aB3w9rKRFR{;x52p<JD;yMXB8FXb&%j>!uGJc6Z-ndT-|VO(e>CO
z@=>a|VKTFxK%kL}CepmX9xxN<+Efy5%^JZ;6_3$X*(0riWKC6+-Bts<-+B5T`ID&f
zn5Q|Zl8K5UO0r*i?`E%HYFl<P2Gap&|3qFmmRwL3;~fgg%bZ@xMEw>Oihr~o-6gsy
zoEsWi!0IJKS%{G7P0T%hU}><GBrJOi;S&3Wk6N#+3Vb3ZlJn1u5EI+{nI=wEoN^mf
z@mP-P_^?4ve^~Kk(#5f6QVhkK%3Mv&KUAzT=|D^NNwL^SUoJw=-*7Yc+zI>QHt7qx
z?M?fH8ia!VXh{$^VpJ^Fo2%Pad_7QQ4lpt25C09o)lxgCNeAu%BC13vS|Gr8l^yZo
z#$;S3DwPHiB=k;dbsg@@PfDmRs3SD`CMO=V+En^0mIiSn6k=-98~@uTF}7Zd=^?eF
zSd3{z19OSZxg6)U?wa<}FtoqIOYJF<`kMLHQ+6EGY7&IUZS2j<^p}7kqNY1vlY}xt
zB6s>v4h@Lv;0Lz~zRD_MM@rQGa~pS;e(5iTu@M)!n*ZcuKRY3NlB{`q+j3R=EcpKj
zc#<7AC{!hLGwY}Qe+qER$^DBnCiY%iG9#k2c%r5_?(TMoDr%`+jSu8D{gH0K@}JxK
z6Q6tAzWh%=KGzgi{@af}i@yjAU9T*(YuFDskN$@aOiuiN{(l4Rmnh%Wzr{_d9gmWp
zksfU*(>3(PQ&Wi8%*MjOfrx*~TDXh?laq%_rxBxki#tc>4%TKEX5y8B!!m$8tgLKZ
z+3Ve=T%2NwWW!L@q-$<ul6v_v&jRR~8~&(0!(360DK6uF>cxHsepT7q*`|z^dI#>+
zH8(dYBK%wnB81FM>USM})h{MTT=5WUS;9+dd)5m~C9s~=nmHkI0aNl#Z$jN!yN5qc
z&I&tJOvJ6~2W3Q+)_P6;MN$>`M;i)EYdp_!X1$9VS6K*cTM5yAh&heO5jw45`M)0F
zd3{7x>AXeI^Y;<ZwYJ*3g87=F9%sT2Ze1DaWU%LPBVLyOv;|Da$3OEYyvpRP@t=_V
z8kQb+TDJt?FNO^7|1w0+2ATVV8*kJ9!)IWqEUgXPET|3ple1GHl=$CV;>e6P{#4KJ
z0515`CRoY;BLM%8!Swd85JRpYL@H6}>}9C;kQpWNn_mQ!Z1Ym$M>1Ckh2K9-mbAgB
zJcvQ#AIb-<*ZizU&Uo6?LQMXRei=EtFI)}Fs%UZPCD;mrh3$$%E-sLm@ixY>?u`n!
zFwiXoK0QRRIjTO(@(+QK5Nw|N8BaS~_z`T&`?qWYt0K_y`7N99maUeU@i^A80m0V1
zC2LFBU~pAUAjoF7WVZP5M%8A)%mMe!7UFf%WMQsaYezTUl+{442BJnokJ&xq3-j`2
z8koqj%H)zK;K`CX@1EHg-aMONYyf^sTzt=X?0(t=1^<Wsz`aivgpGof%h78s*=q>t
z<7GYKtdf4~@g=*=&-OUvKP)6dt!twUv*SHwcuZ99v^7{`Zgp!|GPnPc#kIPVT{amD
z_(zuX1NXP*Y^o6-F4w|60O=9)J|i&sq5BwXnv!gXsw-B^Z(<m<!<ct$HR*LVAq799
zjhz!4bVYYBM;_Wqaa&?jTwtqR+vg{JFBk>EzA;x{Q-AUijrg`&-5++l%c|`zE3DvO
zF5O%oig+WiU23Tp>t+#Ap*_aDxNHAdmJt?!O*OpNn3w3krnNheY@abN`F|~&b+e(B
zDhs_QMz96FRo}fuUAc1`rs@!jzJ!r6Z1HR-O>o!g{@D1*{?&AYraQ~)nT%S26tL~!
zn6du3%}#CCx-#BpW8eek#AO6;SE-ef2k+z5&O-lerp?#-C&GZjaZ%DQW1-9^N(<9s
zs>ze9_qX3mzgy6*lMyDPP+EW?Y=`L;{y7}}T0^>9cHykALQVuVG)ahS(r{V>wT4pR
zcG;|IfsNa{+I-ol^c4JEqenllwIj%R#lBFrcYWbch&l3Z*rwV$A02(OI`VT(`$A&(
zeEmwV_TJ7p<&HGh(SeIv%a~z=2*35ZU@+Br;K>0;iwBW&a^rSjo|89h>;tfW((t>w
zZdp$*3lqqz8d>*m?y{V&5!^ZQ`~Zg?QQu8%HCZrvTM?^yp;-Se1@PQ01?ZVn&k5qa
zdSQPPlHwzueLsoj8A_V2KecXI6Ktox+OoS%A-fU#JIvwcBA!s*LY!F>eQ+T5%l@n_
zwC}Rd=uP20NOTOh*WRoL9!u$<Z>vg9H2{$H<7j3kKN^x5(p=AZy56^Rvalvc`(i4d
zD-KkkU0OfUme~qay%o{L<0B?w>8Zo7qamYOoz_XQqAy1yYK8m%m}Gmvk%MZ^K;46v
zg<TA0+vqfV?<`6GIV$|DoU(fIQDuHYH_`LXopO9485BGY(e?8VR4;WY+b=Lp`QJ&>
z-(3XhA9g65p6+3f@*Z8V6a{NNO`hH0qG=z=uI^L&{_#lzAk!9qLrY50mPFa|tsCdl
zP2<y_3NfobW<4CiDn^|3#_mn)23wa8bd@MR4oxtnSL`u;92*#OXuA|Ntj?QO;&{-<
zbbRAMd=T4*lU0&vl{IOw^bJEUzCC>g`5gL9<7ns2*m4YBHFOV4luD8*mB%6Z^#fy^
z?Wk41t%@R=O#pCIZ$Xi4mf2p&8<UvAxMtOW@>}J>yFTvYLYqafH@RQ8#B$nvyU2Ws
z5q*EDR)VWUdj*8lRkFy?D5<eh*-7PB<l{;bD-96zlV0M_v#>=2s;LvFHuT%i#R#g3
z-vsi;GeTU+<OJl}u)Qd!GVSWpmE4}waHNNnMJniu0z4P~9$TEY%kB4VBV~KIsfYz}
z)_0I}BXTQPV8Zs5?COmH-TQ*nf+wd%Bmt>>z0GHmHVFr>JrX$D7V9hC09hOqy=bFi
zou|QH>STMI#AkmEEZ<wt+dHZhKXyLIc<!gqvG{#iiG8WQLZg^LiDL{;-1Qp$#ln{V
zJeXD7185SXpceEAP-wRBYN}v$^kk$?q76VW;13t}45mvy6SIh86`uu#qU};y7N3QW
z1$nKND3NBl>&Ml8k{=(ser}&p$ZiY%BxQ|iltSNH=IC9}N=G4DGETbtqt|;$P{v;G
ztw?Jm8WH*OeBp<Rq9&nEE!B~Ieo1lN9tCzxpbQDV3eSm-MoXbso{#xmS2VYFs;UuI
zWFeod`%R&Wcsghb&B?_y?Xkmo+?rL@VCt2^N@^O$D)?UE2QIV0qMnl7#*(e(njns6
z(rK=Nb5a5$izUa<{QNk)^GnctM9)?4VdpcV()7#TE@{G;S4>*Z+&8qg!iBt;X&dQ{
z#YF7TcDpb+SK!gbs*6%>eK*J>_BV0a8C0s`wzCaPbBEBps7VZsT7ePVaPNXwF@`Q(
z&h^4R8}xHGxG|7CCgOU2+F=azm~D3Pr28Py3<Q&v5`eiQ7WTq-iyeMQ^%IpUA`jAj
zn@27RNDo!UQI|ee_g^N^58{VffHUY14q@RrIVM*l#^>-5!+AD?@x?wk1Br-<E|j=W
zclQH~piO^v_+hy@e5M#)pJ({iy=+k18B+}0u~~Lxet12fUzSY+drQ_dJ?aE(dY1JY
zx@da=&dqH!L*uFBs~h#kUS&?DI=YpO-dP{$+UGxy%4WQMLfpjuy@u#gxo<M;OcQ=8
z`7XI~ET6mZHF-ZoY9J0ZXALe%HGq^KvA|G^n^&R5Op-vK<%arzu{vMCkH|?UY9LDV
zr?n7Aw)kONXV0Nj?#Uxx4A<~Ccm%r(e6_>rfIgYMkp^q_o<0AOTpW$jL~Avnt@xwr
z00Lbz9pqA44=5wgFNg%@pmsLB?Sx?K+I`9~c_yih+UlYbVGDz_vXXp<H^ElJUi(u`
zQ&F>>DQ!XWZ-rugYtr)~txb){pp`=H>vY~<B%0dUxJ~qdBsAu^4lO^jY_>YfG1WdU
zVDg;-dh14xiR4<Z2$Cy>T>8%k+0IV6GYvv5$=+D<b3XP3{A98I=4xjlcp+5uO_E|m
zvx-n%_olB?9qP|cKc)`1l;^yjelms1BdT8Hl#k?3yH6J^v~E3m&hUw4mvhH#RC*do
zq(5lLDshZ;_RJb=mr1VbQAQH2@nleg9k3<=e_|)jz6w54OxD=pc{MZ88I%4a{-1-M
z=LFRPCCWW9D#8HmO{f=5U+#i1cXH(7E{n7>7T1GN(uOx!9Oj`9mLoai5vN5tW;*B<
z5sw9lB1@quNim)!m<vjbIBnN+LtQPg1Uh@aF0nAwu_!wSgz2avAEp8KH#mCuI3Pg>
zXzh%fJv$0ms<<Dg)nn*`QW#w%9b>&nEh+cjF}d_iG1B-*?z|%4jL9Z_$JITGUE>k5
zfgJq|D@o;_1B_zuy^k(=pr#`clSx03sRV`VQceKxKpO8nmDd~zKZNKZeXInOe2uGz
zmmyt7unXm@*wt5HH&)#F#Vc@dV5|hX#i;;Hn(52y-=t;JOb3Z>pqP?8Nxuv9LKGr6
zaH!>hKJ@a#^)MAR`6v}dpQqs$kNWgE5|i`m!erGUuhmZiF9?bWg$vVoiYey3bk2xj
z*oUe~rZX(d!#J7CmcsRCsw6l#bu*i`<#oktWwnj(4+c5_1zX4Z`SGy1wAcw(+vBm*
zB$#KNaq7wxNq{1tmlCqLGw!OB>5tLw?!CWlG{4(_v3c<}LRyz*!&0YJ$ds4~MY|g<
z3OC=x?5xMd468LH`wVpKJv&Do0HOJfp7V3PsG~wnzXz9oK&IVtCmu6ixR6eC-N}&5
zh^#n{P<kJDdtj@g#+gCP(KMW`zv*X?Wy`4i=8?(vLKv|lWIc1I(7!+WIZH?RfpP#J
zK=ihw$y=t1-!RSMqWS=$#kG(V!w^hY)z!LRFNaJ2G&`!ZT2VjX#$f9EOXs$uyc#p<
z{F&gDJ)1?e{a?;s4s4mu9#(%oe9&ja(tuMC0f7w>P<Nj7?nGnHfN#3|yf@694+U<g
zAN50=-1i`rBHWNYITXl>Nca7Cn%IbXb^%Ewz}#noS*hi$%#_-DYpPKM1SiU<aCp?L
z?GJqKuwz-Kkyg-``u>jgqf=^cMYI;_lnPz}Sc}2>8ujwX{;@+pBW>JFslpucbM`(r
zP1XIveOLbY(NRL}MkWWsS>Chk$5oUa&p0t)!nhN{%vu3%fo@kT49V9i2ZQld-=qlE
z^z6g$MQ(iv6O6aiclu2Mwf|=2H{Z7Fzuo%DCj3J0)5y1h5HZ)69B5w7&}C8AfkoT7
zdIdowC{eTZ+$?0hb$DQw(VZ~U$efS7^de%kSv8^REeS0=@2+eEAzu{q4Bh-Z4I~iA
zqkc}D4hc}}qj}iBpmK?Wz03AhWnLck#Nl{_fDoQ{D!b_H2!4uW)9k|G)NEMaFDJNu
zEw~<TRa^aOkvx{oA>}a{{=*=GM%GCNnd&p{?1FCLDtZn_l&8W$GBMld%=np6bKUi&
ze586o*rEys44%<VqxRvoL`EeoaUZ1p+)T2fMF4=8yrT7ViM*SOwz&m);E4}KbGJI~
zbI<3BTr@|c$&`{9SWmR`v-^49p3cC#wU3u1i(|8Oe&Wpe2tz6B-*5EybTs#yfioE#
z7)-G8@|F^&CY1E}Shd?(keHO-ne~Ev8i#(g%<U**TGuQU9hNKZ$vX)cxWjb=H{r-D
zS3ecTyIY64z7lQaAw4>ctl;p}{K1Hua7JLw6=g5K%vmJVVLtJwaN-MpU)~wNJD@=c
zxm1gle%5!&&gdcgpzr{H|NbCBz_Ci*-MqHmO6vJN^nbUiA+^s{k(!cZ*eL$q`a*19
z|NA+`y4?tan92~5|GnP(F9I<Ths4XUQ2ae<{$Kx7$~%!Q84e^;2$`%TQX0famWJZr
z|6)QA=;8n3K`Sf}JXT^zr7SaY2`6NqiyU$&%ZuF0g&018sL3%ScX8iYQ{~7YQSxL+
zbCAz+!W6e}$ra+!sYJ}&|JEDhh0w}NB7+1VmWVHHA;_>CJ!C-y1u`nnjI1OKIhSXk
g`1ip>27&PZ%Eng|BB~&U43dEOD?COWl)82P2hvui8vp<R

diff --git a/logic/test/data/point_sources.csv b/logic/test/data/point_sources.csv
new file mode 100644
index 0000000..c9d2467
--- /dev/null
+++ b/logic/test/data/point_sources.csv
@@ -0,0 +1,19 @@
+type;scenario;year;reporting_country;nfr_sector;pollutant;emission;unit;x;y;hoogte_m;diameter_m;temperatuur_C;warmteinhoud_MW;debiet_Nm3/u;dv;type_emissie;EIL_nummer;exploitatie_naam;EIL_Emissiepunt_Jaar_Naam;Activiteit_type
+historic;;2021;BEF;3B1b;NO2;0.360260113;kg/yr;116332.9766;166636.8709;5;0;0;0;0;4;Runderen jonger dan 1 jaar;57978963_57978963_24628802823;57978963_57978963_24628802823;57978963_57978963_24628802823;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;116332.9766;166636.8709;5;0;0;0;0;4;Runderen jonger dan 1 jaar;57978963_57978963_24628802823;57978963_57978963_24628802823;57978963_57978963_24628802823;Mengmest
+historic;;2021;BEF;3B1b;NO2;1.489021862;kg/yr;116332.9766;166636.8709;5;0;0;0;0;4;Zoogkoeien;57978963_57978963_24628802823;57978963_57978963_24628802823;57978963_57978963_24628802823;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;116332.9766;166636.8709;5;0;0;0;0;4;Zoogkoeien;57978963_57978963_24628802823;57978963_57978963_24628802823;57978963_57978963_24628802823;Mengmest
+historic;;2021;BEF;3B1b;NO2;1.049682095;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Mestkalveren;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Mestkalveren;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
+historic;;2021;BEF;3B1b;NO2;5.038105693;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Runderen (van 1 tot 2 jaar) ;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Runderen (van 1 tot 2 jaar) ;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
+historic;;2021;BEF;3B1b;NO2;34.85470149;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Andere runderen;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Andere runderen;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
+historic;;2021;BEF;3B1b;NO2;42.96467575;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Vervangingsvee (1-2 j.);86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0.015456732;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Vervangingsvee (1-2 j.);86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
+historic;;2021;BEF;3B1b;NO2;0.842354917;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Runderen jonger dan 1 jaar;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Runderen jonger dan 1 jaar;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
+historic;;2021;BEF;3B1b;NO2;1.239607175;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Zoogkoeien;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Zoogkoeien;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
+historic;;2021;BEF;3B1b;NO2;35.77269872;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Vervangingsvee < 1 j.;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Vaste mest
+historic;;2021;BEF;3B1b;NO2;0.042744747;kg/yr;95419.65;196533.8;5;0;0;0;0;4;Vervangingsvee < 1 j.;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;86422801341_86422801341_26422801437;Mengmest
diff --git a/logic/test/data/point_sources_empty_coord.csv b/logic/test/data/point_sources_empty_coord.csv
new file mode 100644
index 0000000..6b316f3
--- /dev/null
+++ b/logic/test/data/point_sources_empty_coord.csv
@@ -0,0 +1,3 @@
+type;scenario;year;reporting_country;nfr_sector;pollutant;emission;unit;x;y;hoogte_m;diameter_m;temperatuur_C;warmteinhoud_MW;debiet_Nm3/u;dv;type_emissie;EIL_nummer;exploitatie_naam;EIL_Emissiepunt_Jaar_Naam;Activiteit_type
+historic;;2021;BEF;3B1b;NO2;0.360260113;kg/yr;116332.9766;166636.8709;5;0;0;0;0;4;Runderen jonger dan 1 jaar;57978963_57978963_24628802823;57978963_57978963_24628802823;57978963_57978963_24628802823;Vaste mest
+historic;;2021;BEF;3B1b;NO2;1.489021862;kg/yr;;;5;0;0;0;0;4;Zoogkoeien;57978963_57978963_24628802823;57978963_57978963_24628802823;57978963_57978963_24628802823;Vaste mest
diff --git a/logic/test/gridprocessingtest.cpp b/logic/test/gridprocessingtest.cpp
index d27a112..95eb95c 100644
--- a/logic/test/gridprocessingtest.cpp
+++ b/logic/test/gridprocessingtest.cpp
@@ -25,6 +25,7 @@ namespace emap::test {
 using namespace inf;
 using namespace doctest;
 using namespace std::string_view_literals;
+namespace gdal = inf::gdal;
 
 TEST_CASE("create_geometry_extent")
 {
diff --git a/logic/test/inputparsertest.cpp b/logic/test/inputparsertest.cpp
index 82ebf67..d2f8e24 100644
--- a/logic/test/inputparsertest.cpp
+++ b/logic/test/inputparsertest.cpp
@@ -5,7 +5,9 @@
 #include "gdx/algo/sum.h"
 #include "gdx/rasteriterator.h"
 #include "infra/algo.h"
+#include "infra/chrono.h"
 #include "infra/test/printsupport.h"
+#include "unitconversion.h"
 
 #include "testconfig.h"
 #include "testconstants.h"
@@ -75,6 +77,7 @@ TEST_CASE("Input parsers")
             cfg.set_year(1990_y);
 
             auto emissions = parse_emissions(EmissionSector::Type::Gnfr, file::u8path(TEST_DATA_DIR) / "_input" / "01_data_emissions" / "inventory" / "reporting_2021" / "totals" / "gnfr_allyears_2021.txt", cfg.year(), cfg, RespectIgnoreList::Yes);
+
             REQUIRE(emissions.size() == 4);
 
             for (auto& em : emissions) {
@@ -176,10 +179,10 @@ TEST_CASE("Input parsers")
 
     SUBCASE("Load point source emissions")
     {
-        cfg.set_combine_identical_point_sources(false);
-
         SUBCASE("nfr sectors")
         {
+            cfg.set_combine_identical_point_sources(false);
+
             const auto emissions = parse_point_sources(file::u8path(TEST_DATA_DIR) / "_input" / "01_data_emissions" / "inventory" / "reporting_2021" / "pointsources" / "pointsource_emissions_2021.csv", cfg);
             REQUIRE(emissions.size() == 4);
 
@@ -190,18 +193,69 @@ TEST_CASE("Input parsers")
             }
 
             {
-                auto noxEmissions = emissions.emissions_with_id(EmissionIdentifier(country::BEF, EmissionSector(sectors::nfr::Nfr1A1a), pollutants::NOx));
-                REQUIRE(noxEmissions.size() == 2);
-                CHECK(noxEmissions[0].coordinate().value() == Coordinate(148450, 197211));
-                CHECK(noxEmissions[1].coordinate().value() == Coordinate(95820, 173080));
+                EmissionIdentifier emId(country::BEF, EmissionSector(sectors::nfr::Nfr1A1a), pollutants::NOx);
+                REQUIRE(emissions.emissions_with_id(emId).size() == 2);
+                REQUIRE(emissions.emissions_with_id_at_coordinate(emId, Coordinate(148450, 197211)).size() == 1);
+                REQUIRE(emissions.emissions_with_id_at_coordinate(emId, Coordinate(95820, 173080)).size() == 1);
             }
             {
-                auto nmvocEmissions = emissions.emissions_with_id(EmissionIdentifier(country::BEF, EmissionSector(sectors::nfr::Nfr1A2c), pollutants::NMVOC));
-                REQUIRE(nmvocEmissions.size() == 2);
-                CHECK(nmvocEmissions[0].coordinate().value() == Coordinate(130643, 159190));
-                CHECK(nmvocEmissions[1].coordinate().value() == Coordinate(205000, 209000));
+                EmissionIdentifier emId(country::BEF, EmissionSector(sectors::nfr::Nfr1A2c), pollutants::NMVOC);
+                REQUIRE(emissions.emissions_with_id(emId).size() == 2);
+                REQUIRE(emissions.emissions_with_id_at_coordinate(emId, Coordinate(130643, 159190)).size() == 1);
+                REQUIRE(emissions.emissions_with_id_at_coordinate(emId, Coordinate(205000, 209000)).size() == 1);
+            }
+        }
+
+        SUBCASE("combine point sources")
+        {
+            const EmissionIdentifier emissionId(countries::BEF, EmissionSector(sectors::nfr::Nfr3B1b), pollutants::NOx);
+
+            {
+                cfg.set_combine_identical_point_sources(true);
+                const auto emissions = parse_point_sources(file::u8path(TEST_DATA_DIR) / "point_sources.csv", cfg);
+                REQUIRE(emissions.size() == 2);
+
+                CHECK(emissions.emission_with_id_at_coordinate(emissionId, Coordinate(116332.9766, 166636.8709)).value().amount() == Approx(to_giga_gram(1.849281975, "kg/yr")));
+                CHECK(emissions.emission_with_id_at_coordinate(emissionId, Coordinate(95419.65, 196533.8)).value().amount() == Approx(to_giga_gram(121.820027319, "kg/yr")));
+            }
+
+            {
+                cfg.set_combine_identical_point_sources(false);
+                const auto emissions = parse_point_sources(file::u8path(TEST_DATA_DIR) / "point_sources.csv", cfg);
+                CHECK(emissions.size() == 11);
+
+                CHECK(emissions.emissions_with_id_at_coordinate(emissionId, Coordinate(116332.9766, 166636.8709)).size() == 2);
+                CHECK(emissions.emissions_with_id_at_coordinate(emissionId, Coordinate(95419.65, 196533.8)).size() == 9);
             }
         }
+
+        SUBCASE("point sources empty coordinate")
+        {
+            cfg.set_combine_identical_point_sources(true);
+            chrono::DurationRecorder duration;
+            CHECK_THROWS_AS(parse_point_sources(file::u8path(TEST_DATA_DIR) / "point_sources_empty_coord.csv", cfg), inf::RuntimeError);
+        }
+
+        /*SUBCASE("point sources perf")
+        {
+            {
+                cfg.set_combine_identical_point_sources(true);
+                chrono::DurationRecorder duration;
+                const auto emissions = parse_point_sources(file::u8path(TEST_DATA_DIR) / "emap_NOx_2021_2024_opslag.csv", cfg);
+                Log::info("CSV Point source parsing with combining took: {}", duration.elapsed_time_string());
+
+                REQUIRE(emissions.size() == 28124);
+            }
+
+            {
+                cfg.set_combine_identical_point_sources(false);
+                chrono::DurationRecorder duration;
+                const auto emissions = parse_point_sources(file::u8path(TEST_DATA_DIR) / "emap_NOx_2021_2024_opslag.csv", cfg);
+                Log::info("CSV Point source parsing without combining  took: {}", duration.elapsed_time_string());
+
+                REQUIRE(emissions.size() == 89358);
+            }
+        }*/
     }
 
     SUBCASE("Load scaling factors")