diff --git a/experiments/rdmsol_analyse.py b/experiments/rdmsol_analyse.py new file mode 100644 index 0000000..51527f3 --- /dev/null +++ b/experiments/rdmsol_analyse.py @@ -0,0 +1,52 @@ +""" +Analyse a PCSE ensemble generated by rdmsol_generate.py. + +Example: + %run experiments/rdmsol_analyse.py outputs/rdmsol -v -c barley +""" +from pathlib import Path + +from matplotlib import pyplot as plt +from tqdm import tqdm + +import fpcup + +# Parse command line arguments +import argparse +parser = argparse.ArgumentParser(description="Run a PCSE ensemble for multiple location with one sowing date.") +parser.add_argument("output_dir", help="folder to load PCSE outputs from", type=Path) +parser.add_argument("-c", "--crop", help="crop to run simulations on", default="barley", choices=("barley", "maize", "sorghum", "soy", "wheat"), type=str.lower) +parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") +args = parser.parse_args() + +# Load the ensemble summary +summary = fpcup.io.load_ensemble_summary_from_folder(args.output_dir) + +# Get the associated RDMSOLs +rdmsol_from_run_id = lambda run_id: int(run_id.split("RDMSOL-")[1].split("_")[0]) +summary["RDMSOL"] = summary.index.to_series().apply(rdmsol_from_run_id) +summary.sort_values("RDMSOL", inplace=True) + +# Plot summary outputs as a function of RDMSOL for each soiltype and site +summary_outputs = ("RD", "DOM", "LAIMAX", "TWSO") + +fig, axs = plt.subplots(nrows=len(summary_outputs), ncols=3, sharex=True, sharey="row", figsize=(10, 10), squeeze=False) + +# Loop over the columns (soil types) first +for ax_col, (soiltype, summary_by_soiltype) in zip(axs.T, summary.groupby("soiltype")): + ax_col[0].set_title(f"Soil type: {soiltype}") + + # Loop over the rows (summary outputs) next + for ax, key in zip(ax_col, summary_outputs): + # Plot a line for each site + summary_by_soiltype.groupby("geometry").plot.line("RDMSOL", key, ax=ax, alpha=0.5, legend=False) + +axs[0, 0].set_xlim(0, 150) +for ax, key in zip(axs[:, 0], summary_outputs): + ax.set_ylabel(key) + +fig.align_xlabels() +fig.align_ylabels() + +plt.show() +plt.close() diff --git a/experiments/rdmsol_generate.py b/experiments/rdmsol_generate.py new file mode 100644 index 0000000..d664ca1 --- /dev/null +++ b/experiments/rdmsol_generate.py @@ -0,0 +1,86 @@ +""" +Run a PCSE ensemble for different values of RDMSOL (maximum rooting depth) to determine its effect. +All available soil types are tested. +A single site, roughly central to the Netherlands, is used. + +Example: + %run experiments/rdmsol_generate.py -v -c barley -n 9 +""" +from pathlib import Path + +from tqdm import tqdm + +import fpcup + +# Parse command line arguments +import argparse +parser = argparse.ArgumentParser(description="Run a PCSE ensemble for different values of RDMSOL (maximum rooting depth) to determine its effect.") +parser.add_argument("-d", "--data_dir", help="folder to load PCSE data from", type=Path, default=fpcup.settings.DEFAULT_DATA) +parser.add_argument("-o", "--output_dir", help="folder to save PCSE outputs to", type=Path, default=fpcup.settings.DEFAULT_OUTPUT / "rdmsol") +parser.add_argument("-c", "--crop", help="crop to run simulations on", default="barley", choices=("barley", "maize", "sorghum", "soy", "wheat"), type=str.lower) +parser.add_argument("-n", "--number", help="number of locations; result may be lower due to rounding", type=int, default=25) +parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") +args = parser.parse_args() + +# Constants +YEAR = 2022 +RDMSOLs = list(range(10, 150, 1)) + +# Generate the coordinates +crs = fpcup.constants.WGS84 +coords = fpcup.site.generate_sites_space(latitude=(51, 53), longitude=(4, 6), n=args.number) +if args.verbose: + print(f"Generated {len(coords)} coordinates") + +### Generate constants +# Crop data +cropdata = fpcup.crop.default +if args.verbose: + print("Loaded crop data") + +# Agromanagement calendars +sowdate = fpcup.agro.sowdate_range(args.crop, YEAR)[0] +agromanagement = fpcup.agro.load_agrotemplate(args.crop, sowdate=sowdate) +if args.verbose: + print("Loaded agro management data:") + print(agromanagement) + +### Loop over sites (coordinates) +failed_runs = [] +for c in tqdm(coords, desc="Sites", unit="site", leave=args.verbose): + ### Generate site-specific data + # Get site data + sitedata = fpcup.site.example(c) + + # Get weather data + weatherdata = fpcup.weather.load_weather_data_NASAPower(c) + + ### Loop over soil data and depths + for soildata in tqdm(fpcup.soil.soil_types.values(), desc="Soil types", unit="type", leave=False, disable=fpcup.RUNNING_IN_IPYTHON): + for depth in tqdm(RDMSOLs, desc="Running models", unit="RDMSOL", leave=False, disable=fpcup.RUNNING_IN_IPYTHON): + soildata["RDMSOL"] = depth + run_id = f"{args.crop}_{soildata.name}_RDMSOL-{depth:d}_dos{sowdate:%Y%j}_lat{c[0]:.7f}-lon{c[1]:.7f}" + + # Combine input data + run = fpcup.model.RunData(sitedata=sitedata, soildata=soildata, cropdata=cropdata, weatherdata=weatherdata, agromanagement=agromanagement, geometry=c, crs=crs, run_id=run_id) + + # Run model + output = fpcup.model.run_pcse_single(run) + + # Save the results to file + try: + output.to_file(args.output_dir) + # If the run failed, saving to file will also fail, so we instead note that this run failed + except AttributeError: + failed_runs.append(run_id) + +# Feedback on failed runs: if any failed, let the user know. If none failed, only let the user know in verbose mode. +print() +if len(failed_runs) > 0: + print(f"Number of failed runs: {len(failed_runs)}/{len(coords)}") +else: + if args.verbose: + print("No runs failed.") + +# Save an ensemble summary +fpcup.io.save_ensemble_summary(args.output_dir, verbose=args.verbose) diff --git a/fpcup/__init__.py b/fpcup/__init__.py index 4f8aa17..bce7a84 100644 --- a/fpcup/__init__.py +++ b/fpcup/__init__.py @@ -1,3 +1,4 @@ from . import aggregate, constants, geo, io, plotting, settings from . import agro, crop, site, soil, weather from .model import run_pcse_single, run_pcse_ensemble +from .tools import RUNNING_IN_IPYTHON diff --git a/fpcup/crop.py b/fpcup/crop.py index c9d3eb3..e9a8b8b 100644 --- a/fpcup/crop.py +++ b/fpcup/crop.py @@ -5,6 +5,19 @@ from ._brp_dictionary import brp_crops_NL2EN +# Parameter names are from "A gentle introduction to WOFOST" (De Wit & Boogaard 2021) and from the YAML crop parameter files +parameter_names = {"TSUMEM": "Temperature sum from sowing to emergence [°C day]", + "TSUM1": "Temperature sum from emergence to anthesis [°C day]", + "TSUM2": "Temperature sum from anthesis to maturity [°C day]", + "TBASEM": "Lower threshold temperature for emergence [°C]", + "TEFFMX": "Maximum effective temperature for emergence [°C]", + "SLATB": "Specific leaf area as a function of DVS [; ha/kg]", + "AMAXTB": "Maximum leaf CO2 assimilation rate as function of DVS [;kg/ha/hr]", + "RDI": "Initial rooting depth [cm]", + "RRI": "Maximum daily increase in rooting depth [cm/day]", + "RDMCR": "Maximum rooting depth [cm]", + } + default = YAMLCropDataProvider() def main_croptype(crop: str) -> str: diff --git a/fpcup/plotting.py b/fpcup/plotting.py index dba8228..f7a638d 100644 --- a/fpcup/plotting.py +++ b/fpcup/plotting.py @@ -329,11 +329,11 @@ def wofost_summary_histogram(summary: Summary, keys: Iterable[str]=KEYS_AGGREGAT plt.close() -def _remove_area_from_keys(keys: Iterable[str]) -> tuple[str]: +def _remove_key_from_keys(keys: Iterable[str], key_to_remove: str) -> tuple[str]: """ - Remove the "area" key from a list of keys. + Remove the key_to_remove from a list of keys. """ - return tuple(k for k in keys if k != "area") + return tuple(k for k in keys if k != key_to_remove) def wofost_summary_geo(data_geo: gpd.GeoDataFrame, keys: Iterable[str]=KEYS_AGGREGATE_PLOT, *, @@ -348,9 +348,11 @@ def wofost_summary_geo(data_geo: gpd.GeoDataFrame, keys: Iterable[str]=KEYS_AGGR """ NEW_FIGURE = (axs is None) - # Don't plot area if this is not available - if NEW_FIGURE and "area" in keys and "area" not in data_geo.columns: - keys = _remove_area_from_keys(keys) + # Don't plot n, area if these are not available + if NEW_FIGURE: + for key_to_remove in ("n", "area"): + if key_to_remove in keys and key_to_remove not in data_geo.columns: + keys = _remove_key_from_keys(keys, key_to_remove) # Configure figure and axes if none were provided if NEW_FIGURE: @@ -403,7 +405,7 @@ def plot_wofost_summary(summary: Summary, keys: Iterable[str]=KEYS_AGGREGATE_PLO f"Available columns: {summary.columns}") else: if "area" in keys: - keys = _remove_area_from_keys(keys) + keys = _remove_key_from_keys(keys, "area") weights = None # Create figure and panels diff --git a/fpcup/soil.py b/fpcup/soil.py index 7e1c0f4..24e2b75 100644 --- a/fpcup/soil.py +++ b/fpcup/soil.py @@ -13,6 +13,25 @@ DEFAULT_SOIL_DATA = DEFAULT_DATA / "soil" +# Parameter names are from "A gentle introduction to WOFOST" (De Wit & Boogaard 2021) and from the CABO file descriptions +parameter_names = {"SOLNAM": "Soil name", + "CRAIRC": "Critical soil air content for aeration [cm^3 / cm^3]", + "SM0": "Soil moisture content of saturated soil [cm^3 / cm^3]", + "SMTAB": "Vol. soil moisture content as function of pF [log(cm) ; cm^3 / cm^3]", + "SMFCF": "Soil moisture content at field capacity [cm^3 / cm^3]", + "SMW": "Soil moisture content at wilting point [cm^3 / cm^3]", + "RDMSOL": "Maximum rootable depth of soil [cm]", + "K0": "Hydraulic conductivity of saturated soil [cm / day]", + "KSUB": "Maximum percolation rate of water to subsoil [cm / day]", + "SOPE": "Maximum percolation rate of water through the root zone [cm / day]", + "SPADS": "1st topsoil seepage parameter deep seedbed", + "SPODS": "2nd topsoil seepage parameter deep seedbed", + "SPASS": "1st topsoil seepage parameter shallow seedbed", + "SPOSS": "2nd topsoil seepage parameter shallow seedbed", + "DEFLIM": "required moisture deficit deep seedbed", + "CONTAB": "10-log hydraulic conductivity as function of pF [log(cm) ; log(cm/day)]", + } + # Dummy with default settings dummy = DummySoilDataProvider() dummy.name = "dummy" diff --git a/fpcup/tools.py b/fpcup/tools.py index 1ccbf09..10fc22a 100644 --- a/fpcup/tools.py +++ b/fpcup/tools.py @@ -5,6 +5,13 @@ from ._typing import Iterable +try: + get_ipython() +except NameError: + RUNNING_IN_IPYTHON = False +else: + RUNNING_IN_IPYTHON = True + def make_iterable(x: object, exclude: Iterable[type]=[str]) -> Iterable: """ Check if x is an iterable.