Skip to content

Commit

Permalink
Merge pull request #158 from aquacropos/gdd_switch_type
Browse files Browse the repository at this point in the history
Gdd switch type
  • Loading branch information
tfoster88 authored Jun 2, 2024
2 parents 3066a2e + 695f117 commit e861bea
Show file tree
Hide file tree
Showing 20 changed files with 261 additions and 116 deletions.
8 changes: 8 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -573,3 +573,11 @@ v7_final_tidying.pptx
outputs/

test_test.py

*.png

no-stress-testing.py

temp-delete.py

gdd_compute_crop_cal_testing_200524.py
4 changes: 4 additions & 0 deletions aquacrop/entities/crop.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ def __init__(self, c_name, planting_date, harvest_date=None, **kwargs):
self.bface = (
0.001165 # WP co2 adjustment parameter given by FACE experiments
)
self.SwitchGDDType = 'mean' # calculate GDD phenology based on mean of CD phenology across entire simulation period (mean/median)

if c_name == "custom":

Expand Down Expand Up @@ -110,6 +111,7 @@ def __init__(self, c_name, planting_date, harvest_date=None, **kwargs):
"PlantMethod",
"CalendarType",
"SwitchGDD",
"SwitchGDDType",
"planting_date",
"harvest_date",
"Emergence",
Expand Down Expand Up @@ -258,6 +260,7 @@ def calculate_additional_params(
("PlantMethod", int64),
("CalendarType", int64),
("SwitchGDD", int64),
("SwitchGDDType", str),
("EmergenceCD", int64),
("Canopy10PctCD", int64),
("MaxRootingCD", int64),
Expand Down Expand Up @@ -377,6 +380,7 @@ def __init__(
2 # Calendar Type (1 = Calendar days, 2 = Growing degree days)
)
self.SwitchGDD = 0 # Convert calendar to gdd mode if inputs are given in calendar days (0 = No; 1 = Yes)
self.SwitchGDDType = 'mean' # calculate GDD phenology based on mean of CD phenology across entire simulation period (mean/median)

self.EmergenceCD = 0
self.Canopy10PctCD = 0
Expand Down
70 changes: 70 additions & 0 deletions aquacrop/entities/crops/crop_params.py
Original file line number Diff line number Diff line change
Expand Up @@ -1471,6 +1471,76 @@
'p_up2': 0.65,
'p_up3': 0.75,
'p_up4': 0.8},
'SugarBeetGDD_UK': {'Aer': 5.0,
'CCx': 0.98,
'CDC': 0.003857,
'CDC_CD': 0.07128,
'CGC': 0.010541,
'CGC_CD': 0.13227,
'CalendarType': 2,
'CropType': 2,
'Determinant': 0.0,
'ETadj': 1.0,
'Emergence': 23.0,
'EmergenceCD': 5.0,
'Flowering': 0.0,
'FloweringCD': 0.0,
'GDD_lo': 0,
'GDD_up': 9.0,
'GDDmethod': 3,
'HI0': 0.75,
'HIstart': 865.0,
'HIstartCD': 71.0,
'Kcb': 1.15,
'Maturity': 2203.0,
'MaturityCD': 142.0,
'MaxRooting': 408.0,
'MaxRootingCD': 43.0,
'Name': 'SugarBeetGDD_UK',
'PlantMethod': 1.0,
'PlantPop': 100000.0,
'PolColdStress': 1,
'PolHeatStress': 1,
'SeedSize': 1.0,
'Senescence': 1704.0,
'SenescenceCD': 116.0,
'SwitchGDD': 0,
'SxBotQ': 0.012,
'SxTopQ': 0.048,
'Tbase': 3.0,
'Tmax_lo': 45.0,
'Tmax_up': 40.0,
'Tmin_lo': 3.0,
'Tmin_up': 8.0,
'TrColdStress': 1,
'Tupp': 25.0,
'WP': 18.0,
'WPy': 100.0,
'YldForm': 1301.0,
'YldFormCD': 70.0,
'YldWC': 20,
'Zmax': 1.0,
'Zmin': 0.3,
'a_HI': 6.0,
'b_HI': -9.0,
'dHI0': 20.0,
'dHI_pre': 0.0,
'exc': -9.0,
'fage': 0.15,
'fshape_r': 1.5,
'fshape_w1': 3.0,
'fshape_w2': 3.0,
'fshape_w3': 3.0,
'fshape_w4': 1,
'fsink': 0.5,
'p_lo1': 0.6,
'p_lo2': 1,
'p_lo3': 1,
'p_lo4': 1,
'p_up1': 0.2,
'p_up2': 0.65,
'p_up3': 0.75,
'p_up4': 0.8},
'SugarCane': {'Aer': 5.0,
'CCx': 0.95,
'CDC': -9.0,
Expand Down
41 changes: 9 additions & 32 deletions aquacrop/initialize/compute_crop_calendar.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
import pandas as pd

from ..entities.modelConstants import ModelConstants
from ..utils.prepare_gdd import prepare_gdd
from typing import TYPE_CHECKING

if TYPE_CHECKING:
Expand All @@ -14,6 +15,7 @@ def compute_crop_calendar(
crop: "Crop",
clock_struct_planting_dates: "DatetimeIndex",
clock_struct_simulation_start_date: str,
clock_struct_simulation_end_date: str,
clock_struct_time_span: "DatetimeIndex",
weather_df: "DataFrame",
) -> "Crop":
Expand Down Expand Up @@ -148,36 +150,11 @@ def compute_crop_calendar(
Tmean = (temp_max + temp_min) / 2
Tmean = Tmean.clip(lower=crop.Tbase)
gdd = Tmean - crop.Tbase

gdd_cum = np.cumsum(gdd)
# Find gdd equivalent for each crop calendar variable
# 1. gdd's from sowing to emergence
crop.Emergence = gdd_cum.iloc[int(crop.EmergenceCD)]
# 2. gdd's from sowing to 10# canopy cover
crop.Canopy10Pct = gdd_cum.iloc[int(crop.Canopy10PctCD)]
# 3. gdd's from sowing to maximum rooting
crop.MaxRooting = gdd_cum.iloc[int(crop.MaxRootingCD)]
# 4. gdd's from sowing to maximum canopy cover
crop.MaxCanopy = gdd_cum.iloc[int(crop.MaxCanopyCD)]
# 5. gdd's from sowing to end of vegetative growth
crop.CanopyDevEnd = gdd_cum.iloc[int(crop.CanopyDevEndCD)]
# 6. gdd's from sowing to senescence
crop.Senescence = gdd_cum.iloc[int(crop.SenescenceCD)]
# 7. gdd's from sowing to maturity
crop.Maturity = gdd_cum.iloc[int(crop.MaturityCD)]
# 8. gdd's from sowing to start of yield_ formation
crop.HIstart = gdd_cum.iloc[int(crop.HIstartCD)]
# 9. gdd's from sowing to start of yield_ formation
crop.HIend = gdd_cum.iloc[int(crop.HIendCD)]
# 10. Duration of yield_ formation (gdd's)
crop.YldForm = crop.HIend - crop.HIstart

# 11. Duration of flowering (gdd's) - (fruit/grain crops only)
if crop.CropType == 3:
# gdd's from sowing to end of flowering
crop.FloweringEnd = gdd_cum.iloc[int(crop.FloweringEndCD)]
# Duration of flowering (gdd's)
crop.Flowering = crop.FloweringEnd - crop.HIstart

crop = prepare_gdd(weather_df,
clock_struct_simulation_start_date,
clock_struct_simulation_end_date,
gdd, crop, crop.SwitchGDDType)

# Convert CGC to gdd mode
# crop.CGC_CD = crop.CGC
Expand All @@ -194,15 +171,15 @@ def compute_crop_calendar(
if tCD <= 0:
tCD = 1

CCi = crop.CCx * (1 - 0.05 * (np.exp((crop.CDC_CD / crop.CCx) * tCD) - 1))
CCi = crop.CCx * (1 - 0.05 * (np.exp(((3.33 * crop.CDC_CD) / (crop.CCx + 2.29)) * tCD) - 1))
if CCi < 0:
CCi = 0

tGDD = crop.Maturity - crop.Senescence
if tGDD <= 0:
tGDD = 5

crop.CDC = (crop.CCx / tGDD) * np.log(1 + ((1 - CCi / crop.CCx) / 0.05))
crop.CDC = ((crop.CCx + 2.29) * np.log((((CCi/crop.CCx) - 1) / -0.05) + 1)) / (3.33 * tGDD)
# Set calendar type to gdd mode
crop.CalendarType = 2

Expand Down
1 change: 1 addition & 0 deletions aquacrop/initialize/compute_variables.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ def compute_variables(
crop,
clock_struct.planting_dates,
clock_struct.simulation_start_date,
clock_struct.simulation_end_date,
clock_struct.time_span,
weather_df,
)
Expand Down
1 change: 1 addition & 0 deletions aquacrop/initialize/read_model_parameters.py
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,7 @@ def read_model_parameters(
crop,
clock_struct.planting_dates,
clock_struct.simulation_start_date,
clock_struct.simulation_end_date,
clock_struct.time_span,
weather_df,
)
Expand Down
125 changes: 125 additions & 0 deletions aquacrop/utils/prepare_gdd.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
import numpy as np
import pandas as pd

def prepare_gdd(weather_df, sim_start, sim_end, gdd, crop, sum_fun):
"""
function to read in GDD and crop data to return
mean or median GDD crop growth stages based on all seasons
in the simulation period
Arguments:
weather_df(DataFrame): weather data across entire simulation timespan
sim_start (str): simulation start date
sim_end (str): simulation end date
gdd (numeric vector): daily gdd values across entire simulation time
crop (Crop object): Crop object containing crop parameters
sum_fun (str): 'mean' or 'median', select method of summarising
multiple years of GDD growth stages calculations
Returns:
crop (Crop object): updated crop object with new GDD growth stages
calculated based on entire simulation period
"""
# convert str to datetime
sim_start_date=pd.to_datetime(sim_start)
sim_end_date=pd.to_datetime(sim_end)

# add gdd as column
assert len(gdd) == len(weather_df), "The length of 'gdd' does not match the number of rows in 'weather_df', check planting date is on or after simulation start date in first year."
weather_df['gdd']=gdd

# Convert mm/dd formatted dates to datetime objects
def parse_mmdd_to_datetime(mmdd_date, year):
mm, dd = map(int, mmdd_date.split('/'))
return pd.to_datetime(f'{year}-{mm:02d}-{dd:02d}')

# Initialize the season counter
season_counter = 1

current_year = sim_start_date.year
planting_date=crop.planting_date

# Determine the first planting date
first_planting_date = parse_mmdd_to_datetime(planting_date, current_year)
if first_planting_date < sim_start_date:
current_year += 1
first_planting_date = parse_mmdd_to_datetime(planting_date, current_year)

# Loop through each year to assign seasons
while first_planting_date <= sim_end_date:
# Determine the end date of the season (a day before the next planting date)
next_planting_date = parse_mmdd_to_datetime(planting_date, current_year + 1)
season_end_date = next_planting_date - pd.Timedelta(days=1)

# If the season end date is beyond the simulation end date, truncate it
if season_end_date > sim_end_date:
season_end_date = sim_end_date

# Update the 'season' column for the current season
weather_df.loc[(weather_df['Date'] >= first_planting_date) & (weather_df['Date'] <= season_end_date), 'season'] = season_counter

# Increment the season counter
season_counter += 1

# Move to the next year
current_year += 1
first_planting_date = next_planting_date

# List of growth stages
growth_stages = [
'Emergence', 'Canopy10Pct', 'MaxRooting', 'MaxCanopy', 'CanopyDevEnd',
'Senescence', 'Maturity', 'HIstart', 'HIend', 'YieldFormation'
]
if crop.CropType == 3:
growth_stages.extend(['FloweringEnd', 'FloweringDuration'])

# Create dictionary of lists to store yearly GDD values for each growth stage
# i.e. create as many lists as there are elements in 'growth_stages',
# named with the values in 'growth_stages' in a dictionary
gdd_lists = {f'{stage}': [] for stage in growth_stages}

# get list of season numbers to iterate across
seasons=weather_df['season'].unique()

# iterate across seasons, calculating GDD to each CD growth stage
# per season and storing in the gdd_lists lists
for season in seasons:
# filter to current season
season_data = weather_df[weather_df['season']==season]

# get cumulative GDD for current season
gdd_cum=np.cumsum(season_data['gdd'])

# Find GDD equivalent for each crop calendar day growth stage
gdd_lists['Emergence'].append(gdd_cum.iloc[int(crop.EmergenceCD)])
gdd_lists['Canopy10Pct'].append(gdd_cum.iloc[int(crop.Canopy10PctCD)])
gdd_lists['MaxRooting'].append(gdd_cum.iloc[int(crop.MaxRootingCD)])
gdd_lists['MaxCanopy'].append(gdd_cum.iloc[int(crop.MaxCanopyCD)])
gdd_lists['CanopyDevEnd'].append(gdd_cum.iloc[int(crop.CanopyDevEndCD)])
gdd_lists['Senescence'].append(gdd_cum.iloc[int(crop.SenescenceCD)])
gdd_lists['Maturity'].append(gdd_cum.iloc[int(crop.MaturityCD)])
gdd_lists['HIstart'].append(gdd_cum.iloc[int(crop.HIstartCD)])
gdd_lists['HIend'].append(gdd_cum.iloc[int(crop.HIendCD)])
gdd_lists['YieldFormation'].append(crop.HIend - crop.HIstart)

# Duration of flowering (gdd's) - (fruit/grain crops only)
if crop.CropType == 3:
flowering_end=gdd_cum.iloc[int(crop.FloweringEndCD)]
# gdd's from sowing to end of flowering
gdd_lists['FloweringEnd'].append(flowering_end)
# Duration of flowering (gdd's)
gdd_lists['FloweringDuration'].append(flowering_end - crop.HIstart)

# calculate mean/median of GDD growth stages using dictionary logic,
# set the attribute to update the crop object
if sum_fun == 'mean':
for stage, values in gdd_lists.items():
setattr(crop, stage, np.mean(values))
elif sum_fun == 'median':
for stage, values in gdd_lists.items():
setattr(crop, stage, np.median(values))

return crop
2 changes: 1 addition & 1 deletion setup.cfg
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[metadata]
name = aquacrop
version = 3.0.7
version = 3.0.8
author = Tim Foster
author_email = timothy.foster@manchester.ac.uk
description = Soil-Crop-Water model based on AquaCrop-OS.
Expand Down
5 changes: 3 additions & 2 deletions tests/test_bug.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
"""
This test is for bug testing
"""
import os
os.environ['DEVELOPMENT'] = 'True'
import unittest

from aquacrop import AquaCropModel, Soil, Crop, InitialWaterContent, GroundWater
from aquacrop.utils import prepare_weather, get_filepath

import os
os.environ['DEVELOPMENT'] = 'True'



class TestModelExceptions(unittest.TestCase):
Expand Down
7 changes: 5 additions & 2 deletions tests/test_co2_timeseries.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,16 @@
"""
Test exeptions in the model.
"""
import os
os.environ['DEVELOPMENT'] = 'True'
import unittest
import pandas as pd


from aquacrop import AquaCropModel, Soil, Crop, InitialWaterContent, GroundWater, CO2
from aquacrop.utils import prepare_weather, get_filepath

import os
os.environ['DEVELOPMENT'] = 'True'


class TestModelExceptions(unittest.TestCase):
"""
Expand Down
Loading

0 comments on commit e861bea

Please sign in to comment.