Skip to content
Snippets Groups Projects
Commit 28ff5bc7 authored by Bianca Wentzel's avatar Bianca Wentzel
Browse files

removed unused expression script

parent f8b9f765
No related branches found
No related tags found
1 merge request!2Plugin template
Pipeline #92652 failed
# SPDX-FileCopyrightText: 2024 Brandenburgische Technische Universität Cottbus-Senftenberg
#
# SPDX-License-Identifier: EUPL-1.2
"""
PET plugin API wrapper.
"""
######################################################
# import commonly used libraries. You can add/remove
# libraries as needed
######################################################
import os
import json
from pathlib import Path
# from evaluation_system.api.exceptions import PluginError
from typing import List, Optional
from evaluation_system.api import plugin
from evaluation_system.api.parameters import (
Directory,
Integer,
ParameterDictionary,
File,
SolrField,
String,
Bool,
Date
)
from evaluation_system.model.user import User
from evaluation_system.misc import logger
import freva
class PET(plugin.PluginAbstract):
######################################################
# Define the details for the plugin
######################################################
__version__ = (2025, 1, 9)
__short_description__: str = (
"Calculation of daily Potential Evatransporation"
)
__category__: Optional[str] = "postproc"
__tags__: Optional[List[str]] = ["Potential Evatranspiration"]
tool_developer = {
"name": "Christian Beier",
"email": "christian.beier@b-tu.de",
}
__long_description__: Optional[
str
] = "This plugin calculates daily Potential Evapotranspiration (PET) based on the methodology outlined in FAO Irrigation and Drainage Paper No. 56 (Crop Evapotranspiration) by Allen, Pereira et al., 1998.\n\nFor more details, refer to the official publication: https://www.fao.org/4/x0490e/x0490e00.htm"
######################################################
# Define the parameters for the plugin * Mandatory
######################################################
__parameters__ = ParameterDictionary(
SolrField(
name="project",
facet="project",
mandatory=True,
help="Project facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
},
),
SolrField(
name="product",
facet="product",
mandatory=True,
help="Product facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
}),
SolrField(
name="institute",
facet="institute",
mandatory=True,
help="Institute facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
},
),
SolrField(
name="model",
facet="model",
mandatory=True,
help="Model facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
},
),
SolrField(
name="experiment",
facet="experiment",
mandatory=True,
help="Experiment facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
},
),
SolrField(
name="realm",
facet="realm",
mandatory=True,
help="Realm facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
},
),
SolrField(
name="ensemble",
facet="ensemble",
mandatory=True,
help="ensemble facet for input data search.",
predefined_facets={
"time_frequency": ["1day"],
"variable": ["tasmax", "tasmin", "sfcWind", "pvap", "rsds", "orog"],
},
),
Date(
name="start_date",
mandatory=True,
help="Insert start date you want to select (YYYY-MM-DD)"
),
Date(
name="end_date",
mandatory=True,
help="Insert end date you want to select (YYYY-MM-DD)"
),
String(
name="region",
default=None,
help="Region box to be processed. The format have to be W,E,S,N (e.g. -180,180,0,90 for NH).\nLeave empty to select the region from input grid."
),
Directory(
name="output_dir",
default="$USER_OUTPUT_DIR",
help="Directory where output will be stored.",
),
Bool(
name="debug",
default=False,
help="Turn on to get better detailed debug information in case you face any issue."
)
)
def find_files(self, config_dict, variables: List[str], time_frequency: str, time: Optional[str] = None):
'''
Selects data files based on the given parameters using the freva databrowser
Parameters:
config_dict (dict): Dictionary containing values used for searching the database
variables (list): List of variables used to search data files within the database
time_frequency (str): Abbreviation for time frequency the data should cover
time (str): Time period the data should cover
Returns:
file_dict (dict): Dictionary containing list of files for each variable
'''
search_keys = (
"project",
"product",
"institute",
"model",
"experiment",
"realm",
"ensemble"
)
search_args = {key: config_dict[key] for key in search_keys}
search_args["time_frequency"] = time_frequency
search_args["time"] = time
file_dict = {}
for variable in variables:
search_args["variable"] = variable
if variable == "orog":
del search_args["ensemble"] ## orog files are part of another ensemble than the other variables
files = list(freva.databrowser(**search_args))
file_dict[variable] = files
return file_dict
def run_tool(self, config_dict=None):
config_dict["files"] = self.find_files(config_dict, ["tasmax", "tasmin", "sfcWind", "pvap", "rsds"], "day", config_dict["start_date"] + " to " + config_dict["end_date"])
config_dict["alti_file"] = self.find_files(config_dict, ["orog"], "fx")["orog"]
if "alti_file" not in config_dict or not config_dict["alti_file"]:
logger.error("Can't find altitude file")
out_dir = self._special_variables.substitute({"d": "$USER_OUTPUT_DIR"})
out_dir = Path(out_dir["d"]) / str(self.rowid)
out_dir.mkdir(exist_ok=True, parents=True)
config_dict["output_dir"] = str(out_dir)
user = User()
scratch = os.getenv("SCRATCH", f"/scratch/{user.getName()[0]}/{user.getName()}")
cache_dir = Path(f"{scratch}/pet/{self.rowid}")
cache_dir.mkdir(exist_ok=True, parents=True)
config_dict["cache_dir"] = str(cache_dir)
json_file = Path(config_dict["cache_dir"]) / "out.json"
error_file = Path(config_dict["cache_dir"]) / "error.log"
with json_file.open("w") as json_f:
json.dump(config_dict, json_f, indent=4)
tool_path = Path(__file__).parent / "src" / "pet" / "calculate_pet.sh"
with error_file.open("w") as error_f:
self.call(
f"{tool_path} {json_file}",
stderr=error_f,
)
return self.prepare_output(config_dict["output_dir"])
######################################################
#
# Calculation of daily Potential Evapotranspiration
# according to: FAO Irrigation and Drainage Paper No. 56 (Crop Evapotranspiration), Allan, Pereira 1998
#
# settings for cdo expressions by Christian Beier, February 2024
# missing/fill value must be set to NaN: cdo setmissval,nan
# no checking for units here
#
######################################################
# R script: 2 min/year
# cdo expr: 34s/year
### DECLARATION OF VARIABLES from combined input file ###
# _... local variable, not in output file
_var_tmax=tasmax;
_var_tmin=tasmin;
e_a=pvap;
_var_wind=sfcWind;
R_s=rsds;
_var_altitude=HSURF;
_var_wind_height=Wind_height;
### NATURAL CONSTANTS ###
c_p=1.013*10^(-3); # specific heat at constant pressure [MJ kg-1 C-1]
epsilon=0.622; # ratio molecular weight of water vapour/dry air
lambda=2.45; # latent heat of vaporization [MJ kg-1]
G_sc=0.0820; # solar constant [MJ m-2 d-1]
sigma=4.903*10^-9; # Stefan-Boltzmann constant [MJ K-4 m-2 d-1]
a=0.23; # albedo reflection coefficient for the hypothetical grass reference crop [-]
pi=3.14159265358979323846; # Pi
### CALCULATIONS ###
# MEAN TEMPERATURE: daily mean temperature [C]
T_mean=(_var_tmax+_var_tmin)/2; #(Eq.9)
# VAPOR PRESSURE DEFICIT [kPa]
pressure=101.3*(((293-0.0065*_var_altitude)/293)^5.26); # atmospheric pressure at altitude z [kPa] (Eq.7) 2dim
gamma=(c_p*pressure/(epsilon*lambda)); # psychrometric constant [kPa C-1] (Eq.8) 3dim? daily values equal
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
##where are time und time_bnds in output file from? input file
e_sat=(0.6108*exp(17.27*_var_tmax/(_var_tmax+237.3)) + 0.6108*exp(17.27*_var_tmin/(_var_tmin+237.3)))/2; # mean saturation vapour pressure [kPa] (Eq.11,12)
e_deficit=e_sat-e_a; # vapour pressure deficit
# SLOPE OF SATURATION VAPOUR PRESSURE CURVE [kPa C-1]
delta=(4098*(0.6108*exp(17.27*T_mean/(T_mean+237.3))))/((T_mean+237.3)^2); #(Eq.13)
# WIND SPEED [m s-1]
u2=_var_wind*4.87/ln(67.8*_var_wind_height-5.42); # wind speed 10m -> 2m above ground (Eq.47)
# COMPONENTS OF RADIATION [MJ m-2 d-1]
# net solar (=shortwave) radiation
R_ns=(1-a)*R_s; #(Eq.38)
# extraterrestrial radiation
lat=clat(_var_tmax); # get latitudes as variable
##cdi warning (cdfCheckVarname): Changed double entry of variable name 'lat' to 'lat_2'!
lat_rad=pi/180*(lat); # convert latitude from [] to [rad] (Eq.22)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
#day_of_year=ctimestep(); # timestep number = day of year (for yearly files, not necassary to set with cdo: days since 1.1.) not in file
year=cyear();
month=cmonth();
day=cday();
N1=floor(275 * month / 9);
N2=floor((month + 9) / 12);
N3=(1 + floor((year - 4 * floor(year / 4) + 2) / 3));
day_of_year=N1-(N2 * N3) + day - 30;
### fehler -1 Tag fr 1900,2100 (kein schaltjahr) ###
dist_earth_sun=1+0.033*cos((2*pi/365)*day_of_year); # inverse relative distance Earth-Sun (Eq.23) not in file
solar_dec=0.409*sin((2*pi/365)*day_of_year-1.39); # solar declination [rad], must be -23.45 <= x <= 23.45 (Eq.24) not in file
test_var=(-tan(lat_rad)*tan(solar_dec)); #(Eq.25)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
# calculate sunset_angle only for valid values of cosine function
#sunset_angle_zw=(test_var>=-1 && test_var<=1)*acos(test_var); # acos (Eq.25)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
#sunset_angle=sunset_angle_zw + (test_var<-1 || test_var>1)*missval(_var_tmax); # NA
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
#sunset_angle_test=acos(test_var); #(Eq.25)
test_var_na=(test_var<-1 || test_var>1); # 1==not valid areas
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
sunset_angle=acos(test_var); #acos(>1)=NA (Eq.25)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
R_a=(24*60)/pi*G_sc*dist_earth_sun*(sunset_angle*sin(lat_rad)*sin(solar_dec)+cos(lat_rad)*cos(solar_dec)*sin(sunset_angle)); # extraterrestrial radiation (Eq.21)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
# clear sky solar radiation [MJ m-2 d-1]
R_so=(0.75+2*10^(-5)*_var_altitude)*R_a; # _var_altitude (HSURF, z_station) (Eq.37)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
# net long wave radiation [MJ m-2 d-1]
T_max_K=_var_tmax+273.15;
T_min_K=_var_tmin+273.15;
test_var2=R_s/R_so; # is limited to be <=1, if greater: will be set to 1 (Eq.39)
#test_var3=(test_var2>1);
#test_var4=test_var3+(test_var2<=1)*test_var2;
test_var3=(test_var2>1)+(test_var2<=1)*test_var2;
R_nl=sigma*((T_max_K^4+T_min_K^4)/2)*(0.34-0.14*sqrt(e_a))*(1.35*test_var3-0.35); #(Eq.39)
# net radiation [MJ m-2 d-1]
R_n=R_ns-R_nl; #(Eq.40)
G=0; # estimated soil heat flux for daily periods [MJ m-2] (Eq.42)
# daily potential evapotranspiration ET0 [mm d-1]
ET0=(0.408*delta*(R_n-G)+gamma*(900/(T_mean+273))*u2*e_deficit)/(delta+gamma*(1+0.34*u2)); #(Eq.6)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment