Test effect of different RDMSOL values

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()

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)

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

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:

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

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"

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.