WorldQual_Lite_TP.py — Hauptskript¶

Quelldatei

src/WorldQual_Lite_TP.py

Das zentrale Modellskript. Berechnet die monatlichen und jährlichen Gesamtphosphor-Frachten (TP) aus acht Eintragspfaden, aggregiert sie auf Einzugsgebietsebene und erzeugt die Ergebnis-CSV sowie den Stacked-Bar-Plot.

GLOSSAR FÜR DEUTSCHE ABKÜRZUNGEN (German Abbreviation Glossary)

ef              = Emissionsfaktor [g TP/Person/Tag] (Emission Factor)
con_*           = Anschlussrate (connection rate) [%]
  - con_prim    = Primär (Primary treatment connection rate)
  - con_sec     = Sekundär (Secondary treatment connection rate)
  - con_tert    = Tertiär (Tertiary treatment connection rate)
  - con_quat    = Quartär (Quaternary treatment connection rate)
  - con_untr    = Unbekannt (Unknown treatment connection rate)
rem_*           = Entfernungsrate (removal rate) [%]
  - rem_prim    = Primäre Eliminationsrate
  - rem_sec     = Sekundäre Eliminationsrate
  - rem_tert    = Tertiäre Eliminationsrate
  - rem_quat    = Quartäre Eliminationsrate
  - rem_untr    = Unbekannte Eliminationsrate
rtf_*           = Return-Flow (Rückfluss) [m³/a]
  - rtf_man     = Industrie-Rückfluss
Ld_*            = Last (Fracht) [t/a oder t/Monat]
  - Ld_ds       = Kanalgebundene Last
  - Ld_mf       = Industrie-Last
  - Ld_sc       = Streusiedlungs-Last
  - Ld_atm      = Atmosphärische Depositions-Last
  - Ld_cw       = Verwitterungs-Last
  - Ld_usr      = Urbane Oberflächenabfluss-Last
ds              = Kanalgebundenes häusliches Abwasser
sc              = Streusiedlungen (nicht kanalisiert)
mf              = Industrieabwässer
atm             = Atmosphärische Deposition
cw              = Chemische Verwitterung
usr             = Urbaner Oberflächenabfluss (Urban Surface Runoff)
Cell_GCRC       = Grid Cell Row Column ID (Rasterzellen-Kennung)
IDReg           = Regions-ID (1=eu, 2=af, 3=as, 4=au, 5=na, 6=sa)
IDScen          = Szenario-ID
HL              = Hydraulic Loading (Hydraulische Belastung) [m/a]
SI              = Boden-Erosions-Index
stp_failure     = Kläranlage-Ausfallrate [%] (STP failure rate)
SPO_treated     = Anteil behandelte Streusiedlungs-Abwässer [%] (Scattered settlement treatment percentage)
UNF             = Unformatted Binary File (WaterGAP-Binärformat)

import time as tm
start_run_time = tm.time()
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

from InputDataFetchFunctions import CountryEmmisionFactor
from InputDataFetchFunctions import CountryConcInReturnFlows
from InputDataFetchFunctions import CountryPopulation
from InputDataFetchFunctions import CountryReturnFlows
from InputDataFetchFunctions import CountryConnectionToTreatment
from InputDataFetchFunctions import RemovalRate
from InputDataFetchFunctions import Cell_ID_To_GCRC
from InputDataFetchFunctions import IDFaoReg_from_Country_Id
from BinaryFileHandler import ReadBin
from BinaryFileHandler import Path_Concatenate
from InputDataFetchFunctions import LivestockExcretionRate

import Paths_and_params as PP

Häusliches Abwasser (kanalgebunden)¶

def DomesticSewered(
    country_tot_pop: float,
    country_urb_pop: float,
    country_rur_pop: float,
    cell_urb_pop: float,
    cell_rur_pop: float,
    con_prim: float,
    con_sec: float,
    con_tert: float,
    con_untr: float,
    rem_prim: float,
    rem_sec: float,
    rem_tert: float,
    rem_untr: float,
    stp_failure: float | None,
    ef: float,
    Cell_monthly_correction_factor: float,
    con_urb: float,
    con_rur: float,
    con_quat: float,
    rem_quat: float
) -> tuple[float, float, float]:

Berechnet die Phosphorbelastung aus kanalgebundenem häuslichen Abwasser.

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus häuslichen Abwässern auf Landes- und Rasterzellen-Ebene unter Berücksichtigung von Anschlussraten, Behandlungsstufen und deren Eliminationsraten.

Parameter

Parameter	Typ	Beschreibung
`country_tot_pop`	`float`	Gesamtbevölkerung des Landes [Personen]
`country_urb_pop`	`float`	Urbane Bevölkerung des Landes [Personen]
`country_rur_pop`	`float`	Ländliche Bevölkerung des Landes [Personen]
`cell_urb_pop`	`float`	Urbane Bevölkerung der Rasterzelle [Personen]
`cell_rur_pop`	`float`	Ländliche Bevölkerung der Rasterzelle [Personen]
`con_prim`	`float`	Anschlussrate primäre Behandlung [%]
`con_sec`	`float`	Anschlussrate sekundäre Behandlung [%]
`con_tert`	`float`	Anschlussrate tertiäre Behandlung [%]
`con_untr`	`float`	Anschlussrate unbekannte Behandlung [%]
`rem_prim`	`float`	Eliminationsrate primäre Behandlung [%]
`rem_sec`	`float`	Eliminationsrate sekundäre Behandlung [%]
`rem_tert`	`float`	Eliminationsrate tertiäre Behandlung [%]
`rem_untr`	`float`	Eliminationsrate unbekannte Behandlung [%]
`stp_failure`	`float \\| None`	Kläranlage-Ausfallrate [%], oder None wenn keine Ausfälle angenommen werden
`ef`	`float`	Emissionsfaktor [g TP/Person/Tag]
`Cell_monthly_correction_factor`	`float`	Korrekturfaktor für monatliche Schwankungen [-]
`con_urb`	`float`	Anschlussrate urbane Bevölkerung [%]
`con_rur`	`float`	Anschlussrate ländliche Bevölkerung [%]
`con_quat`	`float`	Anschlussrate quartäre Behandlung [%]
`rem_quat`	`float`	Eliminationsrate quartäre Behandlung [%]

Rückgabe

tuple[float, float, float]

Ld_ds_cell_monthly : Monatliche Phosphorfracht pro Rasterzelle [t/Monat]
Ld_ds_cell_yearly : Jährliche Phosphorfracht pro Rasterzelle [t/a]
Ld_ds_treated_country_yearly : Jährliche behandelte Fracht auf Landesebene [t/a]

    total_connection = con_prim + con_sec + con_tert + con_untr + con_quat

    frac_prim = con_prim / total_connection * 100
    frac_sec = con_sec / total_connection * 100
    frac_tert = con_tert / total_connection * 100
    frac_quat = con_quat / total_connection * 100
    frac_treat_unknown = con_untr / total_connection * 100

    Ld_ds_before_treatment_country = ef * ((country_urb_pop * con_urb/100) + (country_rur_pop * con_rur/100)) / 1000  # In Tonnen pro Jahr

    if stp_failure is not None and not math.isnan(stp_failure):
        Leftover_percentage_after_Treatment = frac_prim / 100 * (100 - rem_prim*stp_failure/100) / 100 + \
                        frac_sec / 100 * (100 - rem_sec*stp_failure/100) / 100 + \
                        frac_tert / 100 * (100 - rem_tert*stp_failure/100) / 100 + \
                        frac_quat / 100 * (100 - rem_quat*stp_failure/100) / 100 + \
                        frac_treat_unknown / 100 * (100 - rem_untr) / 100
    else:
        Leftover_percentage_after_Treatment = frac_prim / 100 * (100 - rem_prim) / 100 + \
                        frac_sec / 100 * (100 - rem_sec) / 100 + \
                        frac_tert / 100 * (100 - rem_tert) / 100 + \
                        frac_quat / 100 * (100 - rem_quat) / 100 + \
                        frac_treat_unknown / 100 * (100 - rem_untr) / 100

    Ld_ds_treated_country_yearly = Ld_ds_before_treatment_country * Leftover_percentage_after_Treatment

    Ld_ds_cell_yearly = ef * ((cell_urb_pop * con_urb/100) + (cell_rur_pop * con_rur/100)) / 1000 * Leftover_percentage_after_Treatment * Cell_monthly_correction_factor
    Ld_ds_cell_monthly = Ld_ds_cell_yearly / 12

    return Ld_ds_cell_monthly, Ld_ds_cell_yearly, Ld_ds_treated_country_yearly

Industrieabwässer¶

def Manufacturing(
    country_rtf_man: float | None,
    cell_rtf_man: float,
    Conc_mf: float,
    con_sec: float | None,
    con_tert: float | None,
    rem_sec: float,
    rem_tert: float,
    stp_failure: float | None,
    Cell_monthly_correction_factor: float
) -> tuple[float, float]:

Berechnet die Phosphorbelastung aus Industrieabwässern.

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus Industrieeinleitungen basierend auf Rückflussraten und Konzentrationen mit Berücksichtigung von Behandlungsstufen.

Parameter

Parameter	Typ	Beschreibung
`country_rtf_man`	`float \\| None`	Return-Flow-Rate Industrie auf Landesebene [m³/a]
`cell_rtf_man`	`float`	Return-Flow-Rate Industrie pro Rasterzelle [m³/a]
`Conc_mf`	`float`	Phosphorkonzentration in Industrieabwässern [mg/l]
`con_sec`	`float \\| None`	Anschlussrate sekundäre Behandlung [%], oder None
`con_tert`	`float \\| None`	Anschlussrate tertiäre Behandlung [%], oder None
`rem_sec`	`float`	Eliminationsrate sekundäre Behandlung [%]
`rem_tert`	`float`	Eliminationsrate tertiäre Behandlung [%]
`stp_failure`	`float \\| None`	Kläranlage-Ausfallrate [%], oder None
`Cell_monthly_correction_factor`	`float`	Korrekturfaktor für monatliche Schwankungen [-]

Rückgabe

tuple[float, float]

Ld_mf_cell_monthly : Monatliche Phosphorfracht pro Rasterzelle [t/Monat]
Ld_mf_cell_yearly : Jährliche Phosphorfracht pro Rasterzelle [t/a]

    if con_sec == None or con_tert == None:
        frac_sec = 50
        frac_tert = 50
    else:
        frac_sec = con_sec / (con_sec + con_tert)  * 100
        frac_tert = con_tert / (con_sec + con_tert)  * 100

    Ld_mf_untreated_country = Conc_mf * country_rtf_man / 1000000  # Umrechnung in Tonnen pro Jahr

    if stp_failure is not None and not math.isnan(stp_failure):
        Leftover_percentage_after_Treatment =  \
        frac_sec/100*(100 - rem_sec*stp_failure/100)/100 + frac_tert/100*(100 - rem_tert*stp_failure/100)/100
    else:
        Leftover_percentage_after_Treatment =  \
        frac_sec/100*(100 - rem_sec)/100 + frac_tert/100*(100 - rem_tert)/100

Nur Sekundär- und Tertiärbehandlung berücksichtigt. Für europäische Industrien wird 100% Anschluss angenommen (Williams-Paper)

    Ld_mf_treated_country = Ld_mf_untreated_country * Leftover_percentage_after_Treatment
    if country_rtf_man is not None and country_rtf_man != 0:
        Ld_mf_cell_yearly = Ld_mf_treated_country * cell_rtf_man / country_rtf_man * Cell_monthly_correction_factor
    else:
        country_rtf_man = 1000000000
        Ld_mf_cell_yearly = Ld_mf_treated_country * cell_rtf_man / country_rtf_man * Cell_monthly_correction_factor
        print("Error: Invalid country_rtf_man value encountered. A value of 10e9 is assumed.")
    Ld_mf_cell_monthly =  Ld_mf_cell_yearly / 12
    return Ld_mf_cell_monthly , Ld_mf_cell_yearly

Erosionsfaktor nach Fink et al.¶

def Cell_Yearly_ErodedPortion(
    cell: int,
    IDReg: int,
    SoilErosionData: list,
    MeanRunoff_Data: list,
    Yearly_Runoff_Data: list,
    Lmax_value: float,
    a: float,
    b: float,
    c: float,
    Cell_GCRC: int
) -> float:

Berechnet den jährlichen Erosionsfaktor für eine Rasterzelle nach Fink et al.

Diese Funktion ermittelt den Anteil der erodierten Fläche basierend auf hydrologischen Parametern, Bodenerosionsindex und kalibrierbaren Fink-Parametern (Lmax, a, b, c).

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`IDReg`	`int`	Regions-ID (1=EU, 2=AF, 3=AS, 4=AU, 5=NA, 6=SA)
`SoilErosionData`	`list`	Bodenerosionsindex für alle Rasterzellen (SI) [-]
`MeanRunoff_Data`	`list`	Mittlerer Abfluss für alle Rasterzellen [mm/a]
`Yearly_Runoff_Data`	`list`	Jährlicher Abfluss für alle Rasterzellen [mm/a]
`Lmax_value`	`float`	Fink-Parameter: Maximale Abschwemmungslänge [m]
`a`	`float`	Fink-Parameter a [-]
`b`	`float`	Fink-Parameter b [-]
`c`	`float`	Fink-Parameter c [-]
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID der Rasterzelle

Rückgabe

float

Jährlicher Erosionsfaktor (eroded portion) [-]

Cell_GCRC = Cell_ID_To_GCRC(cell, IDReg)

    Cell_yearRunoff = Yearly_Runoff_Data[Cell_GCRC-1]
    Cell_meanRunoff = MeanRunoff_Data[Cell_GCRC-1]                              # Aus waterGAP-Hydrologiemodell
    SI = SoilErosionData[Cell_GCRC-1]

    Cell_eroded_portion_yearly = ((Lmax_value / (1 + (Cell_yearRunoff/a)**(b) )) +
                                  (SI * c * Cell_yearRunoff / Cell_meanRunoff)) # * Cell_net_Area

    return Cell_eroded_portion_yearly

def Cell_Yearly_to_monthly_Load_Converter(
    cell: int,
    country_id: int,
    IDReg: int,
    month: int,
    Actual_Surface_Runoff: list,
    Yearly_Actual_Surface_Runoff: list,
    Cell_GCRC: int
) -> float:

Konvertiert Jahresfrachten zu Monatsfrachten mittels Abflussratio.

Diese Funktion berechnet das Verhältnis zwischen monatlichem und jährlichem Oberflächenabfluss, um die Umrechnung von Jahres- zu Monatsfrachten zu ermöglichen.

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`country_id`	`int`	Länder-ID
`IDReg`	`int`	Regions-ID
`month`	`int`	Monatsnummer (1-12)
`Actual_Surface_Runoff`	`list`	Monatliche Oberflächenabflüsse (12 Monate pro Rasterzelle)
`Yearly_Actual_Surface_Runoff`	`list`	Jährlicher Oberflächenabfluss pro Rasterzelle
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID

Rückgabe

float

Verhältnis Monats- zu Jahresabfluss (monatlich_zu_jährlich_Ratio) [-]

Monatlicher Oberflächenabfluss aus 12er-Array (Monatsdaten)

    Cell_monthly_Actual_Surface_Runoff = Actual_Surface_Runoff[12*Cell_GCRC + month - 13]

Berechne Verhältnis zwischen monatlichem und jährlichem Abfluss

    Cell_monthly_to_yearly_ratio = Cell_monthly_Actual_Surface_Runoff / Yearly_Actual_Surface_Runoff[Cell_GCRC-1]
    return Cell_monthly_to_yearly_ratio

Anorganischer Dünger¶

def Inorganic_Fertilizer_new_method(
    cell: int,
    IDReg: int,
    P_rate_ton_Data: list,
    CropLand_Corrected_Data: list,
    BuiltUPRatioData: list,
    Eroded_portion: float,
    yearly_to_monthly_coeficient: float,
    Cell_GCRC: int
) -> tuple[float, float, float]:

Berechnet die Phosphorbelastung aus anorganischem Dünger.

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus anorganischem Dünger unter Berücksichtigung von Erosionsfaktoren, Anbauflächen und Flächennutzungsmustern.

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`IDReg`	`int`	Regions-ID
`P_rate_ton_Data`	`list`	Phosphordüngeapplikation [t/km²/a] pro Rasterzelle
`CropLand_Corrected_Data`	`list`	Korrigierte Anbaufläche [km²] pro Rasterzelle
`BuiltUPRatioData`	`list`	Anteil bebauter Fläche [0-1] pro Rasterzelle
`Eroded_portion`	`float`	Jährlicher Erosionsfaktor [-]
`yearly_to_monthly_coeficient`	`float`	Jahres-zu-Monats-Konversionsfaktor [-]
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID

Rückgabe

tuple[float, float, float]

Cell_monthly_inorg_load : Monatliche Last [t/Monat]
Cell_Ld_inorg_yearly : Jährliche Last [t/a]
P_rate_ton_km2 : Düngerapplikationsrate [t/km²/a]

Lade Düngeapplikationsrate für Rasterzelle

    P_rate_ton_km2 = P_rate_ton_Data[Cell_GCRC-1]

Lade Anbaufläche

    Cell_CropLand_area = CropLand_Corrected_Data[Cell_GCRC-1]

Lade bebauten Flächenanteil

    BuiltUP_Ratio_Cell = BuiltUPRatioData[Cell_GCRC-1]

Berechne jährliche Last: Düngerrate × Erosion × Anbaufläche × (1 - bebaute Fläche)

    Cell_Ld_inorg_yearly = P_rate_ton_km2 * Eroded_portion * Cell_CropLand_area * (1 - BuiltUP_Ratio_Cell)

Konvertiere zu monatlicher Last

    Cell_monthly_inorg_load = Cell_Ld_inorg_yearly * yearly_to_monthly_coeficient
    return Cell_monthly_inorg_load, Cell_Ld_inorg_yearly, P_rate_ton_km2

Viehwirtschaftlicher TP-Eintrag¶

def AgricultureLivestock(
    cell: int,
    country_id: int,
    IDReg: int,
    Livestock_densityData: list,
    ls_exr_rate: list,
    Eroded_portion: float,
    yearly_to_monthly_coeficient: float,
    Cell_GCRC: int
) -> tuple[float, float]:

Berechnet die Phosphorbelastung aus Viehwirtschaft und Gülleeintrag.

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus Viehbestandsdichten unter Berücksichtigung von monatlichen Ausscheidungsraten und Erosionsfaktoren.

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`country_id`	`int`	Länder-ID
`IDReg`	`int`	Regions-ID
`Livestock_densityData`	`list`	Monatliche Viehbestandsdichten [Tiere/km²]
`ls_exr_rate`	`list`	Monatliche Ausscheidungsraten [t TP/Tier/Monat]
`Eroded_portion`	`float`	Jährlicher Erosionsfaktor [-]
`yearly_to_monthly_coeficient`	`float`	Jahres-zu-Monats-Konversionsfaktor [-]
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID

Rückgabe

tuple[float, float]

Cell_monthly_org_load : Monatliche Last [t/Monat]
Cell_Ld_org_yearly : Jährliche Last [t/a]

Berechne monatliche Viehausscheidungen für alle 12 Monate

    ex = []
    j = 0
    for i in range(12*(Cell_GCRC-1), 12*Cell_GCRC):

Multipliziere Viehbestandsdichte mit Ausscheidungsrate für jeden Monat

        ex.append(Livestock_densityData[i] * ls_exr_rate[j])
        j = j + 1

Summiere alle monatlichen Ausscheidungen (Jahressumme)

    F_org = sum(ex)

Berechne jährliche Last mit Erosionsfaktor

    Cell_Ld_org_yearly = F_org * Eroded_portion

Konvertiere zu monatlicher Last

    Cell_monthly_org_load = Cell_Ld_org_yearly * yearly_to_monthly_coeficient

    return Cell_monthly_org_load, Cell_Ld_org_yearly

Streusiedlungen (nicht kanalisiert)¶

def DomesticNonsewered(
    ef: float,
    cell_rur_pop: float,
    cell_urb_pop: float,
    country_urb_pop: float,
    country_rur_pop: float,
    country_tot_pop: float,
    to_treat_and_unknown: float | None,
    to_hanging_t: float | None,
    to_open_def: float | None,
    rem_untr: float,
    con_prim: float,
    rem_prim: float,
    con_sec: float,
    rem_sec: float,
    con_tert: float,
    rem_tert: float,
    con_untr: float,
    stp_failure: float | None,
    SPO_treated: float | None,
    yearly_to_monthly_coeficient: float,
    treat_failure: float | None,
    con_urb: float,
    con_rur: float,
    con_quat: float,
    rem_quat: float
) -> tuple[float, float, float, float, float]:

Berechnet die Phosphorbelastung aus nicht kanalgebundenem häuslichen Abwasser (Streusiedlungen).

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus nicht an Kläranlagen angeschlossener Bevölkerung mit verschiedenen Entsorgungsvarianten (Behandlung, Hängetsysteme, offene Defäkation).

Parameter

Parameter	Typ	Beschreibung
`ef`	`float`	Emissionsfaktor [g TP/Person/Tag]
`cell_rur_pop`	`float`	Ländliche Bevölkerung der Rasterzelle [Personen]
`cell_urb_pop`	`float`	Urbane Bevölkerung der Rasterzelle [Personen]
`country_urb_pop`	`float`	Urbane Bevölkerung des Landes [Personen]
`country_rur_pop`	`float`	Ländliche Bevölkerung des Landes [Personen]
`country_tot_pop`	`float`	Gesamtbevölkerung des Landes [Personen]
`to_treat_and_unknown`	`float \\| None`	Anteil Streusiedlungen mit Behandlung/unbekannt [%]
`to_hanging_t`	`float \\| None`	Anteil Streusiedlungen mit Hängetsystemen [%]
`to_open_def`	`float \\| None`	Anteil Streusiedlungen mit offener Defäkation [%]
`rem_untr`	`float`	Eliminationsrate unbekannte Behandlung [%]
`con_prim, rem_prim`	`float`	Anschluss- und Eliminationsrate primäre Behandlung [%]
`con_sec, rem_sec`	`float`	Anschluss- und Eliminationsrate sekundäre Behandlung [%]
`con_tert, rem_tert`	`float`	Anschluss- und Eliminationsrate tertiäre Behandlung [%]
`con_untr`	`float`	Anschlussrate unbekannte Behandlung [%]
`stp_failure`	`float \\| None`	Kläranlage-Ausfallrate [%]
`SPO_treated`	`float \\| None`	Anteil behandelte Streusiedlungs-Abwässer [%]
`yearly_to_monthly_coeficient`	`float`	Jahres-zu-Monats-Konversionsfaktor [-]
`treat_failure`	`float \\| None`	Behandlungs-Ausfallrate [%]
`con_urb, con_rur`	`float`	Anschlussraten urbane/ländliche Bevölkerung [%]
`con_quat, rem_quat`	`float`	Anschluss- und Eliminationsrate quartäre Behandlung [%]

Rückgabe

tuple[float, float, float, float, float]

ld_treat_sc_cell_monthly : Monatliche behandelte Last [t/Monat]
ld_treat_sc_cell : Jährliche behandelte Last [t/a]
ld_untr_sc_cell : Jährliche unbehandelte Last [t/a]
ld_hanging_l_cell : Jährliche Last Hängetsysteme [t/a]
ld_diff_untr_sc_cell : Jährliche Last offene Defäkation [t/a]

Validiere und ersetze None/NaN-Werte mit Nullen

    if to_treat_and_unknown == None or math.isnan(to_treat_and_unknown):
        to_treat_and_unknown = 0
    if to_hanging_t == None or math.isnan(to_hanging_t):
        to_hanging_t = 0
    if to_open_def == None or math.isnan(to_open_def):
        to_open_def = 0
    if SPO_treated == None or math.isnan(SPO_treated):
        SPO_treated = 0
    if treat_failure == None or math.isnan(treat_failure):
        treat_failure = 0

Setze STP-Ausfallrate auf Null für Streusiedlungen

    if stp_failure is not None and not math.isnan(stp_failure):
        stp_failure = 0

Berechne Gesamtanschlussrate über alle Behandlungsstufen

    total_connection = con_prim + con_sec + con_tert + con_untr + con_quat

Berechne Anteil nicht angeschlossener urbaner Bevölkerung (behebe Floating-Point-Fehler)

    if ( (100-con_urb)>-0.0001 and (100-con_urb)<0.0001 ):
        x = 0.0
    else:
        x = 100-con_urb

Berechne Anteil nicht angeschlossener ländlicher Bevölkerung

    if ( (100-con_rur)>-0.002 and (100-con_rur)<0.002 ):
        y = 0.0
    else:
        y = 100-con_rur

Berechne Gesamtzahl nicht angeschlossener Personen

    NotConnectedPeople = ((country_urb_pop * x/100) + (country_rur_pop * y/100))

Berechne Basislast für nicht angeschlossene Bevölkerung

    ld_untr_sc = ef * NotConnectedPeople / 1000

Wenn detaillierte Informationen über Entsorgungswege vorhanden sind

    if to_hanging_t != None and to_open_def != None and to_treat_and_unknown != None:

Summiere alle Entsorgungswege (müssen 100% ergeben)

        sum_ld_sc = to_treat_and_unknown + to_hanging_t + to_open_def
        if (100 - (total_connection + sum_ld_sc)) > -0.0001 and (100 - (total_connection + sum_ld_sc)) < 0.0001:
            miss_con_rate = 0
        else:
            miss_con_rate = 100 - total_connection - sum_ld_sc
        sum_ld_sc += miss_con_rate
        if sum_ld_sc != 0:
            ld_diff_untr_sc = to_open_def / (sum_ld_sc / 10000)
            ld_treat_sc = (SPO_treated /(sum_ld_sc/100)) * (100 - rem_sec*treat_failure/100) +\
                    ((to_treat_and_unknown + miss_con_rate -SPO_treated)/(sum_ld_sc/100))*(100.0-treat_failure)
        else:
            ld_treat_sc = 0

    else:

Weniger detaillierte Informationen: Annahme nur Primär- und Sekundärbehandlung

        if total_connection == 0:
            frac_sec = 100
            frac_prim = 0
            frac_tert = 0
            frac_quat = 0
            frac_treat_unknown = 0
        else:
            frac_prim = con_prim / total_connection * 100
            frac_sec = (con_sec  + con_tert + frac_quat)/ total_connection * 100
            frac_tert = 0
            frac_quat = 0
            frac_treat_unknown = con_untr / total_connection * 100

        ld_diff_untr_sc = frac_treat_unknown * (100 - rem_untr) / 100
        if treat_failure != 0:
            ld_treat_sc = frac_prim * (100 - rem_prim*treat_failure/100) / 100 + \
                        frac_sec / 100 * (100 - rem_sec*treat_failure/100) / 100 + \
                        frac_tert / 100 * (100 - rem_tert*treat_failure/100) / 100 + \
                        frac_quat / 100 * (100 - rem_quat*treat_failure/100) / 100
        else:
            ld_treat_sc = frac_prim * (100 - rem_prim) / 100 + \
                        frac_sec / 100 * (100 - rem_sec) / 100 + \
                        frac_tert / 100 * (100 - rem_tert) / 100 + \
                        frac_quat / 100 * (100 - rem_quat) / 100

    ld_treat_sc *= ld_untr_sc / 10000
    ld_diff_untr_sc *= ld_untr_sc / 10000
    ld_hanging_l = (ld_untr_sc * to_hanging_t / (sum_ld_sc / 100)) / 100

Skaliere auf Zellenebene herunter

    ld_treat_sc_cell = ld_treat_sc * ((cell_urb_pop *( 1-con_urb / 100)) + (cell_rur_pop * (1-con_rur / 100))) / (country_tot_pop*(1 - total_connection/100))

    ld_untr_sc_cell = ld_untr_sc * ((cell_urb_pop * x / 100) + (cell_rur_pop * y / 100)) / (country_tot_pop*(1 - total_connection/100))

    ld_diff_untr_sc_cell = ld_diff_untr_sc * ((cell_urb_pop * x / 100) + (cell_rur_pop * y / 100)) / (country_tot_pop*(1 - total_connection/100))

    ld_hanging_l_cell = ld_hanging_l * ((cell_urb_pop * x / 100) + (cell_rur_pop * y / 100)) / (country_tot_pop*(1 - total_connection/100))
    ld_treat_sc_cell_monthly = yearly_to_monthly_coeficient * ld_treat_sc_cell

    return ld_treat_sc_cell_monthly, ld_treat_sc_cell, ld_untr_sc_cell, ld_hanging_l_cell, ld_diff_untr_sc_cell

Atmosphärische Deposition¶

def BackgroundAtm(
    cell: int,
    country_id: int,
    IDReg: int,
    P_atm_Dep_Data: list,
    Eroded_portion: float,
    yearly_to_monthly_coeficient: float,
    GR_Data: list,
    Area_Data: list,
    Land_Area_PercentageData: list,
    BuiltUPRatioData: list,
    Cell_GCRC: int
) -> tuple[float, float]:

Berechnet die Phosphorbelastung aus atmosphärischer Deposition.

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus atmosphärischer Deposition unter Berücksichtigung von Landflächen und Erosionsfaktoren.

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`country_id`	`int`	Länder-ID
`IDReg`	`int`	Regions-ID
`P_atm_Dep_Data`	`list`	Atmosphärische Phosphordeposition [t/km²/a]
`Eroded_portion`	`float`	Jährlicher Erosionsfaktor [-]
`yearly_to_monthly_coeficient`	`float`	Jahres-zu-Monats-Konversionsfaktor [-]
`GR_Data`	`list`	Reihen-Indizes für Rasterzellen
`Area_Data`	`list`	Rasterzellenflächen [km²]
`Land_Area_PercentageData`	`list`	Anteil Landfläche [%]
`BuiltUPRatioData`	`list`	Anteil bebauter Fläche [0-1]
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID

Rückgabe

tuple[float, float]

Cell_monthly_atm_load : Monatliche Last [t/Monat]
Cell_Ld_atm_yearly : Jährliche Last [t/a]

Lade Depositionsrate für Rasterzelle

    Fatm = P_atm_Dep_Data[Cell_GCRC-1]

Lade Rasterzellengröße

    Cell_GR = GR_Data[Cell_GCRC-1]
    Cell_Area = Area_Data[Cell_GR-1]

Berechne Landfläche (ohne Wasserflächen)

    Cell_Land_Area = Cell_Area * int(Land_Area_PercentageData[Cell_GCRC-1]) / 100

Lade bebauten Flächenanteil

    BuiltUP_Ratio_Cell = BuiltUPRatioData[Cell_GCRC-1]

Berechne natürliche Landfläche (ohne bebaute Fläche)

    Cell_net_Area = Cell_Land_Area * (1 - BuiltUP_Ratio_Cell)

Berechne jährliche Last: Depositionsrate × Erosion × natürliche Fläche

    Cell_Ld_atm_yearly = Fatm * Eroded_portion * Cell_net_Area / 1000

Konvertiere zu monatlicher Last

    Cell_monthly_atm_load = Cell_Ld_atm_yearly * yearly_to_monthly_coeficient

    return Cell_monthly_atm_load, Cell_Ld_atm_yearly

Chemische Verwitterung¶

def BackgroundCW(
    cell: int,
    country_id: int,
    IDReg: int,
    PWeathering_Data: list,
    Area_Data: list,
    Land_Area_PercentageData: list,
    BuiltUPRatioData: list,
    Yearly_Runoff_Data: list,
    MeanRunoff_Data: list,
    yearly_to_monthly_coeficient: float,
    Cell_GCRC: int,
    GR_Data: list
) -> tuple[float, float]:

Berechnet die Phosphorbelastung aus chemischer Verwitterung.

Diese Funktion ermittelt die monatliche und jährliche Phosphorfracht aus Gesteinsverwitterung unter Berücksichtigung von Abflussvariabilität und Erosionsfaktoren.

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`country_id`	`int`	Länder-ID
`IDReg`	`int`	Regions-ID
`PWeathering_Data`	`list`	Phosphor-Verwitterungsrate [t/km²/a]
`Area_Data`	`list`	Rasterzellenflächen [km²]
`Land_Area_PercentageData`	`list`	Anteil Landfläche [%]
`BuiltUPRatioData`	`list`	Anteil bebauter Fläche [0-1]
`Yearly_Runoff_Data`	`list`	Jährlicher Abfluss [mm/a]
`MeanRunoff_Data`	`list`	Mittlerer Abfluss [mm/a]
`yearly_to_monthly_coeficient`	`float`	Jahres-zu-Monats-Konversionsfaktor [-]
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID
`GR_Data`	`list`	Reihen-Indizes für Rasterzellen

Rückgabe

tuple[float, float]

Ld_cw_monthly : Monatliche Last [t/Monat]
Ld_cw_yearly : Jährliche Last [t/a]

Lade Verwitterungsrate für Rasterzelle

    Ld_cw = PWeathering_Data[Cell_GCRC-1]

Lade Rasterzellengröße

    Cell_GR = GR_Data[Cell_GCRC-1]
    Cell_Area = Area_Data[Cell_GR-1]

Berechne Landfläche

    Cell_Land_Area = Cell_Area * int(Land_Area_PercentageData[Cell_GCRC-1]) / 100

Lade bebauten Flächenanteil

    BuiltUP_Ratio_Cell = BuiltUPRatioData[Cell_GCRC-1]

Berechne natürliche Landfläche

    Cell_net_Area = Cell_Land_Area * (1 - BuiltUP_Ratio_Cell)

Lade jährliche und mittlere Abflussmengen

    Cell_yearRunoff = Yearly_Runoff_Data[Cell_GCRC-1]
    Cell_meanRunoff = MeanRunoff_Data[Cell_GCRC-1]

Berechne jährliche Last mit Abflussratio

    Ld_cw_yearly = Ld_cw * Cell_yearRunoff / Cell_meanRunoff * Cell_net_Area / 1000

Konvertiere zu monatlicher Last

    Ld_cw_monthly = Ld_cw_yearly * yearly_to_monthly_coeficient
    return Ld_cw_monthly, Ld_cw_yearly

Urbaner Oberflächenabfluss¶

def UrbanSurfaceRunoff(
    cell: int,
    country_id: int,
    IDReg: int,
    month: int,
    UrbanRunoffData: list,
    EventMeanConc: float,
    Area_Data: list,
    Land_Area_PercentageData: list,
    BuiltUPRatioData: list,
    con_prim: float,
    con_sec: float,
    con_tert: float,
    con_untr: float,
    rem_prim: float,
    rem_sec: float,
    rem_tert: float,
    rem_untr: float,
    stp_failure: float | None,
    Cell_GCRC: int,
    GR_Data: list,
    con_quat: float,
    rem_quat: float
) -> float:

Berechnet die Phosphorbelastung aus urbanen Oberflächenabflüssen (Niederschlagsabfluss).

Diese Funktion ermittelt die monatliche Phosphorfracht aus urbanen Oberflächenabflüssen unter Berücksichtigung von Behandlungsstufen, Ereignismittelkonzentrationen und bebauten Flächenanteilen.

Parameter

Parameter	Typ	Beschreibung
`cell`	`int`	Rasterzellen-ID
`country_id`	`int`	Länder-ID
`IDReg`	`int`	Regions-ID
`month`	`int`	Monatsnummer (1-12)
`UrbanRunoffData`	`list`	Monatliche urbane Abflussmengen [mm/a]
`EventMeanConc`	`float`	Ereignismittelkonzentration Phosphor [mg/l]
`Area_Data`	`list`	Rasterzellenflächen [km²]
`Land_Area_PercentageData`	`list`	Anteil Landfläche [%]
`BuiltUPRatioData`	`list`	Anteil bebauter Fläche [0-1]
`con_prim, con_sec, con_tert, con_untr, con_quat`	`float`	Anschlussraten verschiedener Behandlungsstufen [%]
`rem_prim, rem_sec, rem_tert, rem_untr, rem_quat`	`float`	Eliminationsraten verschiedener Behandlungsstufen [%]
`stp_failure`	`float \\| None`	Kläranlage-Ausfallrate [%], oder None
`Cell_GCRC`	`int`	Gitter-Reihe-Spalte-ID
`GR_Data`	`list`	Reihen-Indizes für Rasterzellen

Rückgabe

float

Monatliche Phosphorfracht aus urbanen Oberflächenabflüssen [t/Monat]

Berechne Gesamtanschlussrate und Fraktionen nach Behandlungsstufe

    total_connection = con_prim + con_sec + con_tert + con_untr + con_quat
    frac_prim = con_prim / total_connection * 100
    frac_sec = con_sec / total_connection * 100
    frac_tert = con_tert / total_connection * 100
    frac_quat = con_quat / total_connection * 100
    frac_treat_unknown = con_untr / total_connection * 100

Lade monatliche urbane Abflussmenge für Rasterzelle und Monat

    Cell_monthly_Urban_Surface_Runoff = UrbanRunoffData[12*Cell_GCRC + month - 13]

Lade Rasterzellengröße

    Cell_GR = GR_Data[Cell_GCRC-1]
    Cell_Area = Area_Data[Cell_GR-1]

Berechne Landfläche

    Cell_Land_Area = Cell_Area * Land_Area_PercentageData[Cell_GCRC-1] / 100

Lade bebauten Flächenanteil

    BuiltUP_Ratio_Cell = BuiltUPRatioData[Cell_GCRC-1]

Berechne unbehandelte urbane Oberflächenlast

    Ld_usr = Cell_monthly_Urban_Surface_Runoff * EventMeanConc * Cell_Land_Area * BuiltUP_Ratio_Cell / 1000

Berechne behandelte Last abhängig von STP-Ausfallrate

    if stp_failure is not None and not math.isnan(stp_failure):

Mit Ausfallrate: Eliminationsrate wird mit Ausfallrate multipliziert

        Ld_usr_cell = Ld_usr * frac_prim / 100 * (100 - rem_prim*stp_failure/100) / 100 + \
                     Ld_usr * frac_sec / 100 * (100 - rem_sec*stp_failure/100) / 100 + \
                     Ld_usr * frac_tert / 100 * (100 - rem_tert*stp_failure/100) / 100 + \
                     Ld_usr * frac_quat / 100 * (100 - rem_quat*stp_failure/100) / 100 + \
                     Ld_usr * frac_treat_unknown / 100 * (100 - rem_untr) / 100
    else:

Ohne Ausfallrate: volle Eliminationsraten

        Ld_usr_cell = Ld_usr * frac_prim / 100 * (100 - rem_prim) / 100 + \
                     Ld_usr * frac_sec / 100 * (100 - rem_sec) / 100 + \
                     Ld_usr * frac_tert / 100 * (100 - rem_tert) / 100 + \
                     Ld_usr * frac_quat / 100 * (100 - rem_quat) / 100 + \
                     Ld_usr * frac_treat_unknown / 100 * (100 - rem_untr) / 100

    return Ld_usr_cell


def Load_After_Retention_factor(
    HL: float,
    a_ret: float = PP.a_ret,
    b_ret: float = PP.b_ret
) -> float:

Berechnet den Retentionsfaktor nach Vollenweider für Seen und Flussstrecken.

Diese Funktion ermittelt den Anteil der Phosphorfracht, die die Seeretention passiert, basierend auf der hydraulischen Belastung (HL) und kalibrierten Retentionsparametern (a_ret, b_ret) nach Vollenweider.

Parameter

Parameter	Typ	Beschreibung
`HL`	`float`	Hydraulische Belastung (Durchfluss/Fläche) [m/a]
`a_ret`	`float, optional`	Kalibrier-Parameter a für Retentionsberechnung (aus Paths_and_params.PP.a_ret)
`b_ret`	`float, optional`	Kalibrier-Parameter b für Retentionsberechnung (aus Paths_and_params.PP.b_ret)

Rückgabe

float

Retentionsfaktor (0-1): Anteil der Fracht, der nicht zurückgehalten wird [-] 1 = keine Retentionierung, 0 = vollständige Retentionierung

Berechne Retentionsfaktor nach Vollenweider-Formel

    return 1 / (1 + a_ret*HL**b_ret) if HL != 0 else 1

print("Functions definition completed.")

Filterparameter aus Konfiguration¶

initial_year = PP.initial_year
final_year_included = PP.final_year_included
years = range(initial_year,final_year_included+1)
parameter_id = PP.parameter_id                                  # Phosphor = 60 (wq_load_general)
months = PP.months                                              # Februar = 2 (Monate 1-12)
country_id = PP.country_id                                      # Länder-ID (wq_general)
IDScen = PP.IDScen                                              # Szenario-ID = 9
dbname_cell = PP.dbname_cell                                    # Zellen-Datenbankname (Europa)
dbname1 = PP.dbname1                                            # Database name where country inputs and parameters are stored
continent_index = PP.continent_index                            # Index-Nummer des Kontinents aus der 'name'-Liste unten
IDReg = 1                                                       # eu=1, af=2, as=3, au=4,na=5,sa=6
Point_load_corr = PP.Point_load_corr


name = PP.name
ng =  PP.ng                                                     # Total cells in the continent (Default)
nrow =  PP.nrow                                                 # Number of rows forming the matrix including water outside the continent (Default)
ncol =  PP.ncol                                                 # Number of columns forming the matrix including water outside the continent (Default)


Basin_cells_with_percentage = pd.read_csv(PP.Basin_cells_list_csv_path)                                                                   # Cell_id (Grid cell)
Basin_grid_cells = list(Basin_cells_with_percentage['Cell_ID'])
years = range(initial_year,final_year_included+PP.time_step, PP.time_step)
IDFAOReg = IDFaoReg_from_Country_Id(country_id)                                 # OECD from Database= "wq_general" and Table= "_fao_reg"


GLCC_path = PP.Other_UNF_files_folder + "/GLCC2000.UNF2"
GC_path = PP.Other_UNF_files_folder +  "/GC.UNF2"
GR_path = PP.Other_UNF_files_folder +  "/GR.UNF2"
CellAreaPath = PP.Other_UNF_files_folder + "/GAREA.UNF0"
LandAreaPercentagePath = PP.Other_UNF_files_folder + "/G_LAND_AREA.UNF1"
BuiltUPRatioPath = PP.Other_UNF_files_folder + "/GBUILTUP.UNF0"
SoilErosionPath = PP.Other_UNF_files_folder + "/G_SOILEROS.UNF0"
MeanRunoff_path = PP.Surface_Runoff_folder + "/G_SURFACE_RUNOFF_MEAN.UNF0"
Gfreqw_path = PP.Other_UNF_files_folder + "/GFREQW.UNF1"                             # Wasserflächen-Anteil in den Zellen [%]
P_atm_Dep_Path = PP.Other_UNF_files_folder + "/G_PATMDEPOS.UNF0"
PWeathering_Path = PP.Other_UNF_files_folder + "/G_PWEATHERING.UNF0"

GLCC_data = ReadBin(GLCC_path, ng[continent_index])

GR_Data = ReadBin(GR_path, ng[continent_index])                                 # Liste der Reihen jeder Zelle, sortiert nach GCRC
MeanRunoff_Data = ReadBin(MeanRunoff_path, ng[continent_index])
Area_Data = ReadBin(CellAreaPath, nrow[continent_index])                        # Sortiert nach Zellenreihe
Land_Area_PercentageData = ReadBin(LandAreaPercentagePath, ng[continent_index])
BuiltUPRatioData = ReadBin(BuiltUPRatioPath, ng[continent_index])
SoilErosionData = ReadBin(SoilErosionPath, ng[continent_index])
Gfreqw_Data = ReadBin(Gfreqw_path, ng[continent_index])                  # Wasserflächen-Anteil in jeder Zelle [%]
P_atm_Dep_Data = ReadBin(P_atm_Dep_Path, ng[continent_index])
PWeathering_Data = ReadBin(PWeathering_Path, ng[continent_index])

Modellberechnung¶

Die Funktion Model() ist die zentrale Berechnungsroutine. Sie nimmt die sechs Erosions-/Kalibrierungsparameter (Lmax, a, b, c, sc_corr, bg_corr) entgegen und iteriert ueber alle Jahre, Monate und Rasterzellen im Einzugsgebiet. Pro Zelle werden die acht TP-Quellen berechnet und auf Einzugsgebietsebene aggregiert.

def Model(X: list) -> float:

Hauptberechnungsroutine für das TP-Belastungsmodell.

Diese Funktion ist das Kernstück des WorldQual Lite-Modells. Sie iteriert über alle Jahre, Monate und Rasterzellen im Einzugsgebiet und berechnet die acht Phosphor-Quellen: (1) kanalgebundene Abwässer, (2) Industrie, (3) anorganische Dünger, (4) Viehwirtschaft, (5) Streusiedlungen, (6) atmosphärische Deposition, (7) chemische Verwitterung und (8) urbane Oberflächenabflüsse. Die Ergebnisse werden auf Einzugs- gebietsebene aggregiert und mit Retentionsfaktoren korrigiert.

Parameter

Parameter	Typ	Beschreibung
`X`	`list`	Kalibrier-Parameter als Liste mit 6 Elementen: - X[0] : Lmax (Fink-Parameter, maximale Abschwemmungslänge) [m] - X[1] : a (Fink-Parameter) [-] - X[2] : b (Fink-Parameter) [-] - X[3] : c (Fink-Parameter) [-] - X[4] : sc_corr (Streusiedlungs-Korrekturfaktor) [-] - X[5] : bg_corr (Hintergrund-Korrekturfaktor) [-]

Rückgabe

float

Mean Squared Error (MSE) zwischen simulierten und beobachteten Werten (oder ähnliche Gütemaßzahl für Kalibrierung)

    Mean_Squared_Error_list = []

    column_names_monthly_sum = ['Year', 'DomesticSeweredLoad', 'Scattered Load', 'ManufacturingLoad',
                                'InorganicFertilizerLoad', 'AgricultureLivestockLoad', 'AtmBackgroundLoad',
                                'CWBackgroundLoad', 'UrbanSurfaceRunoffLoad', 'SummedLoadings']

    monthly_sum_by_year = pd.DataFrame(columns = column_names_monthly_sum)

    for time in years:

Schritt 1: Zellinput-Dateien laden¶

Lade CSV-Dateien mit Bevölkerung, Rückflüssen und anderen Zelleninputs

        if PP.run_type == 'Future':
            df_continent_cell_input = pd.read_csv(Path_Concatenate(PP.cell_input_folder + "/" + PP.Scenario + "/europe_cell_input_", time, ".csv"))
        elif PP.run_type == 'Historical':
            df_continent_cell_input = pd.read_csv(Path_Concatenate(PP.cell_input_folder + "/europe_cell_input_", time, ".csv"))

Schritt 2: UNF-Dateipfade zusammensetzen¶

Konstruiere vollständige Pfade zu Binärdateien basierend auf Jahr und Lauftypus

        if PP.run_type == 'Historical':
            Runoff_path = Path_Concatenate(PP.Surface_Runoff_folder + "/G_SURFACE_RUNOFF_", time, ".12.UNF0")
            UrbanRunoffPath = Path_Concatenate(PP.Urban_Runoff_folder + "/G_URBAN_RUNOFF_", time, ".12.UNF0")
            Livestock_Density_Path = Path_Concatenate(PP.Livestock_Density_folder + "/G_LIVESTOCK_NR_", time, ".12.UNF0")
            Correction_Factor_Path = Path_Concatenate(PP.Correction_Factor_folder + "/G_CORR_FACT_RTF_", time, ".12.UNF0")
            P_Rate_ton_path = Path_Concatenate(PP.P_Rate_ton_folder + "/P_RATE_TON_KM2_", time, ".UNF0")
            CropLand_Corrected_Path = Path_Concatenate(PP.CropLand_Corrected_folder + "/CROPLAND_CORR_KM2_", time, ".UNF0")
        elif PP.run_type == 'Future':
            Runoff_path = Path_Concatenate(PP.Future_UNF_files + '/Hydrology/' + PP.name[continent_index] + '/' + PP.rcp + '_' + PP.GCM + "/G_SURFACE_RUNOFF_", time, ".12.UNF0")
            UrbanRunoffPath = Path_Concatenate(PP.Future_UNF_files + '/Hydrology/' + PP.name[continent_index] + '/' + PP.rcp + '_' + PP.GCM + "/G_URBAN_RUNOFF_", time, ".12.UNF0")
            Livestock_Density_Path = Path_Concatenate(PP.Future_UNF_files + '/LIVESTOCK_NR/' + PP.Scenario + '/' + PP.name[continent_index] + '/' + "/G_LIVESTOCK_NR_", time, ".12.UNF0")
            Correction_Factor_Path = Path_Concatenate(PP.Future_UNF_files + '/correction_factors/' + PP.rcp + '_' + PP.Scenario + '/' + PP.name[continent_index] + '/' + PP.GCM + "/G_CORR_FACT_RTF_", time, ".12.UNF0")
            P_Rate_ton_path = Path_Concatenate(PP.Future_UNF_files + '/P_RATE_TON_KM2/' + PP.Scenario + '/' + PP.name[continent_index] + "/P_RATE_TON_KM2_"+PP.Scenario+"_",time,"_"+PP.name[continent_index]+".UNF0")
            CropLand_Corrected_Path = Path_Concatenate(PP.Future_UNF_files + '/CROPLAND_AREA_KM2/' + PP.Scenario + '/' + PP.name[continent_index] + "/"+PP.Scenario+"_",time, "_5arcmin_cropland_"+PP.name[continent_index]+".UNF0")

Schritt 3: Erosionsparameter entpacken¶

Extrahiere Fink-Parameter aus Kalibrier-Vektor X

        Lmax = X[0]
        a = X[1]
        b = X[2]
        c = X[3]
        sc_corr = X[4]
        bg_corr = X[5]

Schritt 4: UNF-Rasterdaten einlesen¶

Lade binäre Rasterkarten: Abfluss, Vieh, Dünger, Korrekturfaktoren etc.

        Runoff_Data = ReadBin(Runoff_path, ng[continent_index])                         # Monatliche Daten
        UrbanRunoffData = ReadBin(UrbanRunoffPath, ng[continent_index])                 # Monatliche Daten
        Livestock_densityData = ReadBin(Livestock_Density_Path, ng[continent_index])
        Correction_Factor_Data = ReadBin(Correction_Factor_Path, ng[continent_index])   # Monatliche Daten
        P_rate_ton_Data = ReadBin(P_Rate_ton_path, ng[continent_index])
        CropLand_Corrected_Data = ReadBin(CropLand_Corrected_Path, ng[continent_index])

        print("Input data from UNF files have been fetched without error.")
        print("Pre-processing of input data in progress . . .")

Schritt 5: Abflussdaten vorverarbeiten¶

Berechne aktuellen Oberflächenabfluss und aggregiere zu Jahressummen

        Actual_Surface_Runoff = []
        zip_object = zip(Runoff_Data, UrbanRunoffData)
        for list1_i, list2_i in zip_object:
            difference = list1_i-list2_i
            if difference < 0:
                difference = 0
            Actual_Surface_Runoff.append(difference)

        Yearly_Runoff_Data = []
        Yearly_Actual_Surface_Runoff = []
        ActSurfaceRunoffSummer = []
        for i in range(0 , len(Runoff_Data) , 12):
            Yearly_Runoff_Data.append(sum(Runoff_Data[i:(i+12)]))
            Yearly_Actual_Surface_Runoff.append(sum(Actual_Surface_Runoff[i:(i+12)]))
            ActSurfaceRunoffSummer.append(sum(Actual_Surface_Runoff[(i+1):(i+9)]))     # Zweiter bis neunter Monat

Schritt 6: Länderdaten laden (Datenbank oder CSV)¶

Lade nationale Parameter: Bevölkerung, Anschlussraten, Eliminationsraten

        if PP.data_source  == 'DB':
            country_urb_pop = CountryPopulation(dbname1, IDScen, time, country_id)[1]
            country_rur_pop = CountryPopulation(dbname1, IDScen, time, country_id)[2]
            country_tot_pop = CountryPopulation(dbname1, IDScen, time, country_id)[0]
            country_urb_pop_percentage = country_urb_pop / country_tot_pop * 100

country_rur_pop_percentage = country_rur_pop / country_tot_pop * 100

            ef = CountryEmmisionFactor(dbname1, parameter_id, time, country_id)[0]            # g TP/Person/Tag

ef = 0.4

            con_prim = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[0]                 # In Prozent
            con_sec = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[1]                  # In Prozent
            con_tert = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[2]                 # In Prozent
            con_untr = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[3]                 # In Prozent
            stp_failure = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[4]              # In Prozent
            con_quat = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[11]
            con_quat = (con_quat if con_quat is not None else 0)
            con_tot = con_prim + con_sec + con_tert + con_untr + con_quat                                           # Gesamtanschlussrate des Landes [%]


            rem_prim = RemovalRate(dbname1, IDScen, time, parameter_id)[0]
            rem_sec = RemovalRate(dbname1, IDScen, time, parameter_id)[1]
            rem_tert = RemovalRate(dbname1, IDScen, time, parameter_id)[2]
            rem_quat = RemovalRate(dbname1, IDScen, time, parameter_id)[8]
            rem_quat = (rem_quat if rem_quat is not None else 0)
            rem_untr = RemovalRate(dbname1, IDScen, time, parameter_id)[3]
            treat_failure = RemovalRate(dbname1, IDScen, time, parameter_id)[4]
            to_treat_and_unknown = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[5]
            to_hanging_t = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[6]
            to_open_def = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[7]
            UrbSewerConn = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[8]
            RurSewerConn = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[9]
            SPO_treated = CountryConnectionToTreatment(dbname1, IDScen, time, country_id)[10]

            country_rtf_man = CountryReturnFlows(dbname1, IDScen, time, country_id)[0]        # mg/l
            Conc_mf = CountryConcInReturnFlows(dbname1, parameter_id, time, country_id)[6]    # conc_man_nd (mg/l)
            EventMeanConc = CountryConcInReturnFlows(dbname1, parameter_id, time, country_id)[8]     # Event mean concentration (conc_urb)

        elif PP.data_source  == 'Excel':
            df_country_input = pd.read_csv(str(PP.data_path) + "/country_input/{}_country_input.csv".format(PP.Scenario))
            df_country_parameter_input = pd.read_csv(str(PP.data_path) + "/country_parameter_input/{}_country_parameter_input.csv".format(PP.Scenario))
            df_parameter_input = pd.read_csv(str(PP.data_path) + "/parameter_input/{}_parameter_input.csv".format(PP.Scenario))

Laenderinput: Bevoelkerung und Klaeranlagenanschluesse

            criteria = (df_country_input['IDScen'] == PP.Future_IDScen) & (df_country_input['country_id'] == country_id) & (df_country_input['time'] == time)

Access the value using the .loc[] method

            country_urb_pop = df_country_input.loc[criteria, 'pop_urb'].iloc[0]
            country_rur_pop = df_country_input.loc[criteria, 'pop_rur'].iloc[0]
            country_tot_pop = df_country_input.loc[criteria, 'pop_tot'].iloc[0]
            country_urb_pop_percentage = country_urb_pop / country_tot_pop * 100


            con_prim = df_country_input.loc[criteria, 'con_prim'].iloc[0]
            con_sec = df_country_input.loc[criteria, 'con_sec'].iloc[0]
            con_tert = df_country_input.loc[criteria, 'con_tert'].iloc[0]
            con_quat = df_country_input.loc[criteria, 'con_quat'].iloc[0]
            con_untr = df_country_input.loc[criteria, 'con_untr'].iloc[0]
            stp_failure = df_country_input.loc[criteria, 'stp_failure'].iloc[0]
            country_rtf_man = df_country_input.loc[criteria, 'rtf_man'].iloc[0]
            con_tot = con_prim + con_sec + con_tert + con_untr

            to_treat_and_unknown = df_country_input.loc[criteria, 'to_treat_and_unknown'].iloc[0]
            to_hanging_t = df_country_input.loc[criteria, 'to_hanging_t'].iloc[0]
            to_open_def = df_country_input.loc[criteria, 'to_open_def'].iloc[0]
            UrbSewerConn = df_country_input.loc[criteria, 'UrbSewerConn'].iloc[0]
            RurSewerConn = df_country_input.loc[criteria, 'RurSewerConn'].iloc[0]
            SPO_treated = df_country_input.loc[criteria, 'SPO_treat'].iloc[0]

Laenderparameter: Emissionsfaktor, Konzentrationen

            criteria2 = (df_country_parameter_input['IDScen'] == PP.Future_IDScen) & (df_country_parameter_input['country_id'] == country_id)\
                        & (df_country_parameter_input['parameter_id'] == parameter_id) & (df_country_parameter_input['time'] == time)

            ef = df_country_parameter_input.loc[criteria2, 'ef'].iloc[0]
            Conc_mf = df_country_parameter_input.loc[criteria2, 'conc_man_nd'].iloc[0]
            EventMeanConc = df_country_parameter_input.loc[criteria2, 'conc_urb'].iloc[0]

Eliminationsraten nach Behandlungsstufe

            criteria3 = (df_parameter_input['IDScen'] == PP.Future_IDScen) & (df_parameter_input['parameter_id'] == parameter_id) & (df_parameter_input['time'] == time)

            rem_prim = df_parameter_input.loc[criteria3, 'rem_prim'].iloc[0]
            rem_sec = df_parameter_input.loc[criteria3, 'rem_sec'].iloc[0]
            rem_tert = df_parameter_input.loc[criteria3, 'rem_tert'].iloc[0]
            rem_untr = df_parameter_input.loc[criteria3, 'rem_untr'].iloc[0]
            treat_failure = df_parameter_input.loc[criteria3, 'treat_failure'].iloc[0]
            rem_quat = df_parameter_input.loc[criteria3, 'rem_quat'].iloc[0]
        else:
            "Define the source of sanitation and pop data in the Paths_and_params.py"

Schritt 7: Anschlussraten Stadt/Land berechnen¶

Berechne getrennte Anschlussraten für urbane und ländliche Bevölkerung

if country_urb_pop_percentage < con_tot:                                                # To give priority for urban population
    con_urb = 100                                                                       # Percentage of Urban population connected to sewer
    con_rur = country_rur_pop/country_tot_pop * (con_tot - country_urb_pop_percentage)  # Percentage of Rural population connected to sewer
else:
    con_rur = 0                                                                         # Percentage of Rural population connected to sewer
    con_urb = country_tot_pop/country_urb_pop * con_tot                                 # Percentage of Urban population connected to sewer

        if country_tot_pop ==0:
            con_rur = 0
            con_urb = 0
        elif country_rur_pop == 0:
            con_rur = 0
            con_urb = UrbSewerConn
        elif country_urb_pop ==0:
            con_rur = RurSewerConn
            con_urb = 0
        else:
            con_rur = RurSewerConn
            con_urb = UrbSewerConn

        ls_exr_rate = LivestockExcretionRate(dbname1, parameter_id, IDFAOReg)

        print("Preprocessing of the input data from UNF files has been done without error.")

Ergebnis-DataFrame für monatliche und jährliche TP-Frachten

        column_names = ['Cell_ID', 'Month', 'DomesticSeweredLoad', 'Scattered Load','ManufacturingLoad', 'InorganicFertilizerLoad',
                          'AgricultureLivestockLoad', 'AtmBackgroundLoad', 'CWBackgroundLoad', 'UrbanSurfaceRunoffLoad','SummedLoadings']

        df_summary = pd.DataFrame(columns = column_names)

Schritt 8: Hydraulische Belastung und Retention berechnen¶

Berechne Wasserfläche und HL für das gesamte Einzugsgebiet

Berechne Retentionsfaktor nach Vollenweider für Seeretention
Hydraulische Belastung (HL) = Jahresabfluss / Wasserfläche [m/a]

        Basin_Water_Surface_Area = 0
        Basin_annual_runoff = 0
        for cell in Basin_grid_cells:
            Cell_GCRC = Cell_ID_To_GCRC(cell, IDReg)
            Cell_GR = GR_Data[Cell_GCRC-1]
            percentage_part_of_cell_in_the_basin = float(Basin_cells_with_percentage[Basin_cells_with_percentage['Cell_ID'] == cell]['Portion of Cell in Basin (%)'])
            Cell_Water_Area = Gfreqw_Data[Cell_GCRC-1] / 100 * Area_Data[Cell_GR-1] * percentage_part_of_cell_in_the_basin/ 100.823 # km²
            Basin_Water_Surface_Area = Basin_Water_Surface_Area + Cell_Water_Area
            cell_annual_runoff_Km3 = Yearly_Runoff_Data[Cell_GCRC-1] /(10**6) * Area_Data[Cell_GR-1] * percentage_part_of_cell_in_the_basin /100.823 # km³
            Basin_annual_runoff = Basin_annual_runoff + cell_annual_runoff_Km3
        HL = Basin_annual_runoff/ Basin_Water_Surface_Area * 1000  if Basin_Water_Surface_Area != 0 else 0             # km³ / km² = km = 1000m
        retention_factor  = Load_After_Retention_factor(HL)

Schritt 9: Zelleniteration (8 TP-Quellen berechnen)¶

Iteriere über alle Rasterzellen im Einzugsgebiet

Pro Zelle: Berechne 8 Phosphor-Quellfrachten für jeden Monat

        for cell in Basin_grid_cells:
            is_cell = df_continent_cell_input['cell']== cell
            cell_urb_pop = df_continent_cell_input[is_cell]["pop_urb"].item()
            cell_rur_pop = df_continent_cell_input[is_cell]["pop_rur"].item()

cell_tot_pop = df_continent_cell_input[is_cell]["pop_tot"].item()

            cell_rtf_man = df_continent_cell_input[is_cell]["rtf_man"].item()                # (m³/Jahr)

print("Zelle Industrie-Rückfluss ist {}".format(cell_rtf_man))

            Cell_GCRC = Cell_ID_To_GCRC(cell, IDReg)
            urbanloadyearly = 0  # initialization


            for month in months:

                Cell_monthly_correction_factor = Correction_Factor_Data[12*Cell_GCRC + month - 13]

DomesticSeweredLoad = DomesticSewered(country_tot_pop, country_urb_pop, country_rur_pop, cell_urb_pop, cell_rur_pop, con_prim, con_sec, con_tert, con_untr, rem_prim, rem_sec, rem_tert, rem_untr, stp_failure, ef, Cell_monthly_correction_factor, con_urb, con_rur)

                DomesticSeweredLoad = DomesticSewered(country_tot_pop, country_urb_pop, country_rur_pop, cell_urb_pop, cell_rur_pop, con_prim, con_sec, con_tert, con_untr,
                      rem_prim, rem_sec, rem_tert, rem_untr, stp_failure, ef, Cell_monthly_correction_factor, con_urb, con_rur, con_quat, rem_quat)

print(Ld_ds_Monatslast[0])

                ManufacturingLoad = Manufacturing(country_rtf_man, cell_rtf_man, Conc_mf, con_sec,
                                          con_tert, rem_sec, rem_tert, stp_failure, Cell_monthly_correction_factor)

                Eroded_portion = Cell_Yearly_ErodedPortion(cell, IDReg, SoilErosionData, MeanRunoff_Data,
                                Yearly_Runoff_Data, Lmax, a, b, c, Cell_GCRC)              # Fink-Parameter überprüfen: Lmax, a, b, c
                yearly_to_monthly_coeficient = Cell_Yearly_to_monthly_Load_Converter(cell, country_id, IDReg, month,
                                                                             Actual_Surface_Runoff, Yearly_Actual_Surface_Runoff, Cell_GCRC)
                InorganicFertilizerLoad = Inorganic_Fertilizer_new_method(cell, IDReg, P_rate_ton_Data, CropLand_Corrected_Data, BuiltUPRatioData,
                                 Eroded_portion, yearly_to_monthly_coeficient, Cell_GCRC)

                AgricultureLivestockLoad = AgricultureLivestock(cell, country_id, IDReg, Livestock_densityData, ls_exr_rate,
                                                        Eroded_portion, yearly_to_monthly_coeficient, Cell_GCRC)

                Scattered__treated_Load = DomesticNonsewered(ef, cell_rur_pop, cell_urb_pop, country_urb_pop, country_rur_pop, country_tot_pop,
                               to_treat_and_unknown, to_hanging_t, to_open_def, rem_untr,
                               con_prim, rem_prim, con_sec, rem_sec, con_tert, rem_tert,
                               con_untr, stp_failure, SPO_treated, yearly_to_monthly_coeficient, treat_failure, con_urb, con_rur, con_quat, rem_quat)[0:2]

                BackgroundAtmLoad = BackgroundAtm(cell, country_id, IDReg, P_atm_Dep_Data,
                          Eroded_portion, yearly_to_monthly_coeficient,
                          GR_Data, Area_Data, Land_Area_PercentageData, BuiltUPRatioData, Cell_GCRC)

                BackgroundCWLoad = BackgroundCW(cell, country_id, IDReg, PWeathering_Data, Area_Data, Land_Area_PercentageData, BuiltUPRatioData,
                         Yearly_Runoff_Data, MeanRunoff_Data, yearly_to_monthly_coeficient, Cell_GCRC, GR_Data)

                UrbanSurfaceRunoffLoad =  UrbanSurfaceRunoff(cell, country_id, IDReg, month, UrbanRunoffData, EventMeanConc, Area_Data, Land_Area_PercentageData, BuiltUPRatioData,
                               con_prim, con_sec, con_tert, con_untr, rem_prim, rem_sec, rem_tert, rem_untr, stp_failure, Cell_GCRC, GR_Data, con_quat, rem_quat)
                SummedLoadings = DomesticSeweredLoad[0]*Point_load_corr + Scattered__treated_Load[0] * sc_corr + ManufacturingLoad[0]*Point_load_corr +\
                        InorganicFertilizerLoad[0] + AgricultureLivestockLoad[0] + BackgroundAtmLoad[0] +\
                            BackgroundCWLoad[0] * bg_corr + UrbanSurfaceRunoffLoad
                urbanloadyearly += UrbanSurfaceRunoffLoad

Monatliche Ergebnisse an df_summary anfuegen

                resulting_data = {'Cell_ID': [cell], 'Month':[month],'DomesticSeweredLoad': [DomesticSeweredLoad[0]*Point_load_corr], 'Scattered Load' : [Scattered__treated_Load[0] * sc_corr],
                          'ManufacturingLoad': [ManufacturingLoad[0]*Point_load_corr], 'InorganicFertilizerLoad': [InorganicFertilizerLoad[0]],
                          'AgricultureLivestockLoad': [AgricultureLivestockLoad[0]], 'AtmBackgroundLoad': [BackgroundAtmLoad[0]], 'CWBackgroundLoad': [BackgroundCWLoad[0] *bg_corr],
                          'UrbanSurfaceRunoffLoad': [UrbanSurfaceRunoffLoad],
                          'SummedLoadings': [SummedLoadings] }
                df = pd.DataFrame(resulting_data)

df_summary = df_summary.append(df)

                df_summary = pd.concat([df_summary, df])


                print("Loadings Calculation for cell " + str(cell) + " for year " + str(time) +" for month " + str(month) + " has been completed.")

Jaehrliche Gesamtfracht pro Zelle berechnen und anfuegen

            SummedLoadingsYear = DomesticSeweredLoad[1]*Point_load_corr + Scattered__treated_Load[1] * sc_corr + ManufacturingLoad[1]*Point_load_corr +\
                        InorganicFertilizerLoad[1] + AgricultureLivestockLoad[1] + BackgroundAtmLoad[1] +\
                            BackgroundCWLoad[1] * bg_corr + urbanloadyearly
            resulting_data_yearly = {'Cell_ID': [cell], 'Month':["Yearly load"],'DomesticSeweredLoad': [DomesticSeweredLoad[1]*Point_load_corr], 'Scattered Load' : [Scattered__treated_Load[1] * sc_corr],
                          'ManufacturingLoad': [ManufacturingLoad[1]*Point_load_corr], 'InorganicFertilizerLoad': [InorganicFertilizerLoad[1]],
                          'AgricultureLivestockLoad': [AgricultureLivestockLoad[1]], 'AtmBackgroundLoad': [BackgroundAtmLoad[1]], 'CWBackgroundLoad':[BackgroundCWLoad[1] * bg_corr],
                          'UrbanSurfaceRunoffLoad': [urbanloadyearly],'SummedLoadings': [SummedLoadingsYear] }
            df2 = pd.DataFrame(resulting_data_yearly)

df_summary = df_summary.append(df2)

            df_summary = pd.concat([df_summary, df2])

Schritt 10: Flächenanteil-Korrektur¶

Rasterzellen, die nur teilweise im Einzugsgebiet liegen, werden

anteilig gewichtet (ggf. Kalibrierung für Rastergitter-Fehler)

        lookup_percentage = np.dot(
                        (df_summary['Cell_ID'].values[:,None] == Basin_cells_with_percentage['Cell_ID'].values), # Wenn dies wahr ist,
                        Basin_cells_with_percentage['Portion of Cell in Basin (%)']    # Gib dies als Ausgabe an
                        )
        df_summary['Percentage of cell in the basin'] = lookup_percentage

df_summary = df_summary[['DomesticSeweredLoad', 'Scattered Load','ManufacturingLoad', 'InorganicFertilizerLoad', 'AgricultureLivestockLoad', 'BackgroundLoad', 'UrbanSurfaceRunoffLoad','SummedLoadings']].multiply(df_summary['Percentage of cell in the basin'], axis="index")

        df_summary.loc[:,['DomesticSeweredLoad', 'Scattered Load','ManufacturingLoad', 'InorganicFertilizerLoad',
                     'AgricultureLivestockLoad', 'AtmBackgroundLoad', 'CWBackgroundLoad','UrbanSurfaceRunoffLoad','SummedLoadings']] = df_summary.loc[:,['DomesticSeweredLoad', 'Scattered Load','ManufacturingLoad', 'InorganicFertilizerLoad',
                     'AgricultureLivestockLoad', 'AtmBackgroundLoad','CWBackgroundLoad', 'UrbanSurfaceRunoffLoad','SummedLoadings']].multiply(df_summary.loc[:, 'Percentage of cell in the basin'] / 100.823, axis="index") # Rasterzellen sind zu 100,823% vertreten (Fehler in den Daten)

Schritt 11: Retentionsfaktor anwenden¶

Berücksichtige Phosphor-Rückhaltung (Sedimentation) in stehenden Gewässern nach Vollenweider basierend auf hydraulischer Belastung

        df_summary.loc[:,['DomesticSeweredLoad', 'Scattered Load','ManufacturingLoad', 'InorganicFertilizerLoad',
             'AgricultureLivestockLoad', 'AtmBackgroundLoad', 'CWBackgroundLoad', 'UrbanSurfaceRunoffLoad','SummedLoadings']] = df_summary.loc[:,['DomesticSeweredLoad', 'Scattered Load','ManufacturingLoad', 'InorganicFertilizerLoad',
             'AgricultureLivestockLoad', 'AtmBackgroundLoad', 'CWBackgroundLoad', 'UrbanSurfaceRunoffLoad','SummedLoadings']].multiply(retention_factor, axis="index")

Schritt 12: Aggregation auf Einzugsgebietsebene¶

Aggregiere Zellergebnisse zu monatlichen Einzugsgebiets-Frachten

(Summation aller Zellen) und speichere in Jahres-Dataframe
df_summary_copy = df_summary.copy()

        df_summary = df_summary.drop(['Cell_ID','Percentage of cell in the basin'], axis = 1)

        monthly_sum = df_summary.groupby(['Month']).sum(numeric_only= False)

        monthly_sum['Year'] = time

shift column 'Year' to first position

        first_column = monthly_sum.pop('Year')

insert column using insert(position,column_name,first_column) function

        monthly_sum.insert(0, 'Year', first_column)
        monthly_sum = monthly_sum.drop('Yearly load')    # Entferne letzte Reihe (Jahresfracht)

Jahresergebnisse an Gesamt-DataFrame anfuegen

        monthly_sum_by_year = pd.concat([monthly_sum_by_year, monthly_sum])

        monthly_sum_by_year['Month'] = monthly_sum_by_year.index
        second_column = monthly_sum_by_year.pop('Month')
        monthly_sum_by_year.insert(1, 'Month', second_column)

    return monthly_sum_by_year

Kalibrierungsparameter aus config.yaml laden

Lmax_calib = PP.Lmax_calib
a_calib = PP.a_calib
b_calib = PP.b_calib
c_calib = PP.c_calib
sc_corr_calib = PP.sc_corr_calib
bg_corr_calib = PP.bg_corr_calib

calib_params = [Lmax_calib, a_calib, b_calib, c_calib, sc_corr_calib, bg_corr_calib]
Predicted_df2 = Model(calib_params)

end_run_time = tm.time()
print("The time it took to run the model is ", round(end_run_time - start_run_time, 2), " Seconds")

Ergebnisse als CSV exportieren¶

Predicted_df2[['DomesticSeweredLoad', 'Scattered Load',
        'ManufacturingLoad', 'InorganicFertilizerLoad',
        'AgricultureLivestockLoad', 'AtmBackgroundLoad', 'CWBackgroundLoad',
        'UrbanSurfaceRunoffLoad', 'SummedLoadings']
              ] = Predicted_df2[
        ['DomesticSeweredLoad', 'Scattered Load',
        'ManufacturingLoad', 'InorganicFertilizerLoad',
        'AgricultureLivestockLoad', 'AtmBackgroundLoad', 'CWBackgroundLoad',
        'UrbanSurfaceRunoffLoad', 'SummedLoadings']].astype(float)

Predicted_df2_yearly = Predicted_df2.groupby('Year').sum()
Predicted_df2_yearly = Predicted_df2_yearly.drop(['Month', 'SummedLoadings'], axis=1)

Spalten in logische Reihenfolge bringen

Predicted_df2_yearly = Predicted_df2_yearly.reindex(columns=['DomesticSeweredLoad', 'ManufacturingLoad', 'Scattered Load',
       'InorganicFertilizerLoad', 'AgricultureLivestockLoad',
       'AtmBackgroundLoad', 'CWBackgroundLoad', 'UrbanSurfaceRunoffLoad'])

import os
os.makedirs(PP.output_folder, exist_ok=True)
Predicted_df2_yearly.to_csv(PP.output_folder + '/TP_yearly_results.csv')

Stacked-Bar-Plot der jaehrlichen TP-Frachten¶

df_to_plot = Predicted_df2_yearly
ax = df_to_plot.plot(kind='bar', stacked=True)

Diagrammtitel und Achsenbeschriftungen

ax.set_title('TP Load')
ax.set_xlabel('Year')
ax.set_ylabel('TP load (ton)')

Diagramm anzeigen

plt.show()

Modellvalidierung (RMSE und R²) — nur wenn Messdaten vorhanden¶

_measured_data_path = os.path.join(PP.output_folder, "Measured_Loading.xlsx")
if os.path.exists(_measured_data_path):
    Mesaured_Data = pd.read_excel(_measured_data_path)
    merged_df = pd.merge(Predicted_df2, Mesaured_Data, on=['Year', 'Month'])

    from sklearn.metrics import mean_squared_error
    RMSE = np.sqrt(mean_squared_error(merged_df['TP Load (tons)'], merged_df['SummedLoadings']))

    import scipy.stats as stats
    r, p = stats.pearsonr(merged_df['TP Load (tons)'], merged_df['SummedLoadings'])
    Rsquared = r**2
    print(f"Validierung: RMSE = {RMSE:.4f}, R² = {Rsquared:.4f}")
else:
    print(f"Keine Messdaten gefunden unter {_measured_data_path} — Validierung übersprungen.")

print('RMSE =', np.round(RMSE,2), 'R2 =', np.round(Rsquared, 2))

Vergleich: Messwerte vs. Simulation¶

Mess- und Simulationsdaten zusammenfuehren

merged_df['YearMonth'] = merged_df['Year'].astype(str) + '-' + merged_df['Month'].astype(str)

Linienplot: Simuliert vs. Gemessen

plt.plot(merged_df['YearMonth'], merged_df['SummedLoadings'], label='Simulated')
plt.plot(merged_df['YearMonth'], merged_df['TP Load (tons)'], label='Measured')

Achsenbeschriftung: nur jede 6. Markierung anzeigen

plt.xticks(merged_df['YearMonth'][::6], rotation=45)

Achsenbeschriftung und Legende

plt.xlabel('Year')
plt.ylabel('TP Load (ton/month)')
plt.legend()

plt.text(2, 6.5, f'R2 = {Rsquared:.2f}\nRMSE = {RMSE:.2f}', fontsize=10, color='green')

Plot anzeigen

plt.show()