initialize

aquacrop.initialize

aquacrop.initialize.calculate_HI_linear

calculate_HI_linear(crop_YldFormCD, crop_HIini, crop_HI0, crop_HIGC)

Function to calculate time to switch to linear harvest index build-up, and associated linear rate of build-up. Only for fruit/grain crops.

Reference Manual (pg. 112)

Arguments:

crop_YldFormCD (int):  length of yield formaiton period (calendar days)

crop_HIini (float):  initial harvest index

crop_HI0 (float):  reference harvest index

crop_HIGC (float):  harvest index growth coefficent

Returns:

crop_tLinSwitch (float): time to switch to linear harvest index build-up

crop_dHILinear (float): linear rate of HI build-up

Source code in aquacrop/initialize/calculate_HI_linear.py

def calculate_HI_linear(
    crop_YldFormCD: int,
    crop_HIini: float,
    crop_HI0: float,
    crop_HIGC: float,
) -> Tuple[float, float]:
    """
    Function to calculate time to switch to linear harvest index build-up,
    and associated linear rate of build-up. Only for fruit/grain crops.

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=121" target="_blank">Reference Manual</a> (pg. 112)


    Arguments:

        crop_YldFormCD (int):  length of yield formaiton period (calendar days)

        crop_HIini (float):  initial harvest index

        crop_HI0 (float):  reference harvest index

        crop_HIGC (float):  harvest index growth coefficent


    Returns:

        crop_tLinSwitch (float): time to switch to linear harvest index build-up

        crop_dHILinear (float): linear rate of HI build-up


    """
    # Determine linear switch point
    # Initialise variables
    ti = 0
    tmax = crop_YldFormCD
    HIest = 0
    HIprev = crop_HIini
    # Iterate to find linear switch point
    while (HIest <= crop_HI0) and (ti < tmax):
        ti = ti + 1
        HInew = (crop_HIini * crop_HI0) / (
            crop_HIini + (crop_HI0 - crop_HIini) * np.exp(-crop_HIGC * ti)
        )
        HIest = HInew + (tmax - ti) * (HInew - HIprev)
        HIprev = HInew

    tSwitch = ti - 1

    # Determine linear build-up rate
    if tSwitch > 0:
        HIest = (crop_HIini * crop_HI0) / (
            crop_HIini + (crop_HI0 - crop_HIini) * np.exp(-crop_HIGC * tSwitch)
        )
    else:
        HIest = 0

    dHILin = (crop_HI0 - HIest) / (tmax - tSwitch)

    crop_tLinSwitch = tSwitch
    crop_dHILinear = dHILin

    return crop_tLinSwitch, crop_dHILinear

aquacrop.initialize.calculate_HIGC

calculate_HIGC(crop_YldFormCD, crop_HI0, crop_HIini)

Function to calculate harvest index growth coefficient

Reference Manual (pg. 110)

Arguments:

crop_YldFormCD (int):  length of yield formation period (calendar days)

crop_HI0 (float):  reference harvest index

crop_HIini (float):  initial harvest index

Returns:

crop_HIGC (float): harvest index growth coefficient

Source code in aquacrop/initialize/calculate_HIGC.py

def calculate_HIGC(
    crop_YldFormCD: int,
    crop_HI0: float,
    crop_HIini: float,
) -> float:
    """
    Function to calculate harvest index growth coefficient

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=119" target="_blank">Reference Manual</a> (pg. 110)


    Arguments:

        crop_YldFormCD (int):  length of yield formation period (calendar days)

        crop_HI0 (float):  reference harvest index

        crop_HIini (float):  initial harvest index


    Returns:

        crop_HIGC (float): harvest index growth coefficient


    """
    # Determine HIGC
    # Total yield_ formation days
    tHI = crop_YldFormCD
    # Iteratively estimate HIGC
    HIGC = 0.001
    HIest = 0
    while HIest <= (0.98 * crop_HI0):
        HIGC = HIGC + 0.001
        HIest = (crop_HIini * crop_HI0) / (
            crop_HIini + (crop_HI0 - crop_HIini) * np.exp(-HIGC * tHI)
        )

    if HIest >= crop_HI0:
        HIGC = HIGC - 0.001

    crop_HIGC = HIGC

    return crop_HIGC

aquacrop.initialize.compute_crop_calendar

compute_crop_calendar(crop, clock_struct_planting_dates, clock_struct_simulation_start_date, clock_struct_time_span, weather_df)

Function to compute additional parameters needed to define crop phenological calendar

Reference Manual (pg. 19-20)

Arguments:

crop (Crop):  Crop object containing crop paramaters

clock_struct_planting_dates (DatetimeIndex):  list of planting dates

clock_struct_simulation_start_date (str):  sim start date

clock_struct_time_span (DatetimeIndex):  all dates between sim start and end dates

weather_df (DataFrame):  weather data for simulation period

Returns:

crop (Crop): updated Crop object

Source code in aquacrop/initialize/compute_crop_calendar.py

def compute_crop_calendar(
    crop: "Crop",
    clock_struct_planting_dates: "DatetimeIndex",
    clock_struct_simulation_start_date: str,
    clock_struct_time_span: "DatetimeIndex",
    weather_df: "DataFrame",
) -> "Crop":
    """
    Function to compute additional parameters needed to define crop phenological calendar

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=28" target="_blank">Reference Manual</a> (pg. 19-20)


    Arguments:

        crop (Crop):  Crop object containing crop paramaters

        clock_struct_planting_dates (DatetimeIndex):  list of planting dates

        clock_struct_simulation_start_date (str):  sim start date

        clock_struct_time_span (DatetimeIndex):  all dates between sim start and end dates

        weather_df (DataFrame):  weather data for simulation period


    Returns:

        crop (Crop): updated Crop object



    """

    if len(clock_struct_planting_dates) == 0:
        plant_year = pd.DatetimeIndex([clock_struct_simulation_start_date]).year[0]
        if (
            pd.to_datetime(str(plant_year) + "/" + crop.planting_date)
            < clock_struct_simulation_start_date
        ):
            pl_date = str(plant_year + 1) + "/" + crop.planting_date
        else:
            pl_date = str(plant_year) + "/" + crop.planting_date
    else:
        pl_date = clock_struct_planting_dates[0]

    # Define crop calendar mode
    Mode = crop.CalendarType

    # Calculate variables %%
    if Mode == 1:  # Growth in calendar days

        # Time from sowing to end of vegatative growth period
        if crop.Determinant == 1:
            crop.CanopyDevEndCD = round(crop.HIstartCD + (crop.FloweringCD / 2))
        else:
            crop.CanopyDevEndCD = crop.SenescenceCD

        # Time from sowing to 10% canopy cover (non-stressed conditions)
        crop.Canopy10PctCD = round(
            crop.EmergenceCD + (np.log(0.1 / crop.CC0) / crop.CGC_CD)
        )

        # Time from sowing to maximum canopy cover (non-stressed conditions)
        crop.MaxCanopyCD = round(
            crop.EmergenceCD
            + (
                np.log(
                    (0.25 * crop.CCx * crop.CCx / crop.CC0)
                    / (crop.CCx - (0.98 * crop.CCx))
                )
                / crop.CGC_CD
            )
        )

        # Time from sowing to end of yield_ formation
        crop.HIendCD = crop.HIstartCD + crop.YldFormCD

        # Duplicate calendar values (needed to minimise if
        # statements when switching between gdd and CD runs)
        crop.Emergence = crop.EmergenceCD
        crop.Canopy10Pct = crop.Canopy10PctCD
        crop.MaxRooting = crop.MaxRootingCD
        crop.Senescence = crop.SenescenceCD
        crop.Maturity = crop.MaturityCD
        crop.MaxCanopy = crop.MaxCanopyCD
        crop.CanopyDevEnd = crop.CanopyDevEndCD
        crop.HIstart = crop.HIstartCD
        crop.HIend = crop.HIendCD
        crop.YldForm = crop.YldFormCD
        if crop.CropType == 3:
            crop.FloweringEndCD = crop.HIstartCD + crop.FloweringCD
            # crop.FloweringEndCD = crop.FloweringEnd
            # crop.FloweringCD = crop.Flowering
        else:
            crop.FloweringEnd = ModelConstants.NO_VALUE
            crop.FloweringEndCD = ModelConstants.NO_VALUE
            crop.FloweringCD = ModelConstants.NO_VALUE

        # Check if converting crop calendar to gdd mode
        if crop.SwitchGDD == 1:
            #             # Extract weather data for first growing season that crop is planted
            #             for i,n in enumerate(ParamStruct.CropChoices):
            #                 if n == crop.Name:
            #                     idx = i
            #                     break
            #                 else:
            #                     idx = -1
            #             assert idx > -1

            date_range = pd.date_range(pl_date, clock_struct_time_span[-1])
            weather_df = weather_df.copy()
            weather_df.index = weather_df.Date
            weather_df = weather_df.loc[date_range]
            temp_min = weather_df.MinTemp
            temp_max = weather_df.MaxTemp

            # Calculate gdd's
            if crop.GDDmethod == 1:

                Tmean = (temp_max + temp_min) / 2
                Tmean = Tmean.clip(lower=crop.Tbase, upper=crop.Tupp)
                gdd = Tmean - crop.Tbase

            elif crop.GDDmethod == 2:

                temp_max = temp_max.clip(lower=crop.Tbase, upper=crop.Tupp)
                temp_min = temp_min.clip(lower=crop.Tbase, upper=crop.Tupp)
                Tmean = (temp_max + temp_min) / 2
                gdd = Tmean - crop.Tbase

            elif crop.GDDmethod == 3:

                temp_max = temp_max.clip(lower=crop.Tbase, upper=crop.Tupp)
                temp_min = temp_min.clip(upper=crop.Tupp)
                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

            # Convert CGC to gdd mode
            # crop.CGC_CD = crop.CGC
            crop.CGC = (
                np.log(
                    (((0.98 * crop.CCx) - crop.CCx) * crop.CC0)
                    / (-0.25 * (crop.CCx**2))
                )
            ) / (-(crop.MaxCanopy - crop.Emergence))

            # Convert CDC to gdd mode
            # crop.CDC_CD = crop.CDC
            tCD = crop.MaturityCD - crop.SenescenceCD
            if tCD <= 0:
                tCD = 1

            CCi = crop.CCx * (1 - 0.05 * (np.exp((crop.CDC_CD / crop.CCx) * 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))
            # Set calendar type to gdd mode
            crop.CalendarType = 2

        else:
            crop.CDC = crop.CDC_CD
            crop.CGC = crop.CGC_CD

        # print(crop.__dict__)
    elif Mode == 2:
        # Growth in growing degree days
        # Time from sowing to end of vegatative growth period
        if crop.Determinant == 1:
            crop.CanopyDevEnd = round(crop.HIstart + (crop.Flowering / 2))
        else:
            crop.CanopyDevEnd = crop.Senescence

        # Time from sowing to 10# canopy cover (non-stressed conditions)
        crop.Canopy10Pct = round(crop.Emergence + (np.log(0.1 / crop.CC0) / crop.CGC))

        # Time from sowing to maximum canopy cover (non-stressed conditions)
        crop.MaxCanopy = round(
            crop.Emergence
            + (
                np.log(
                    (0.25 * crop.CCx * crop.CCx / crop.CC0)
                    / (crop.CCx - (0.98 * crop.CCx))
                )
                / crop.CGC
            )
        )

        # Time from sowing to end of yield_ formation
        crop.HIend = crop.HIstart + crop.YldForm

        # Time from sowing to end of flowering (if fruit/grain crop)
        if crop.CropType == 3:
            crop.FloweringEnd = crop.HIstart + crop.Flowering

        # Extract weather data for first growing season that crop is planted
        #         for i,n in enumerate(ParamStruct.CropChoices):
        #             if n == crop.Name:
        #                 idx = i
        #                 break
        #             else:
        #                 idx = -1
        #         assert idx> -1
        date_range = pd.date_range(pl_date, clock_struct_time_span[-1])
        weather_df = weather_df.copy()
        weather_df.index = weather_df.Date

        weather_df = weather_df.loc[date_range]
        temp_min = weather_df.MinTemp
        temp_max = weather_df.MaxTemp

        # Calculate gdd's
        if crop.GDDmethod == 1:

            Tmean = (temp_max + temp_min) / 2
            Tmean = Tmean.clip(lower=crop.Tbase, upper=crop.Tupp)
            gdd = Tmean - crop.Tbase

        elif crop.GDDmethod == 2:

            temp_max = temp_max.clip(lower=crop.Tbase, upper=crop.Tupp)
            temp_min = temp_min.clip(lower=crop.Tbase, upper=crop.Tupp)
            Tmean = (temp_max + temp_min) / 2
            gdd = Tmean - crop.Tbase

        elif crop.GDDmethod == 3:

            temp_max = temp_max.clip(lower=crop.Tbase, upper=crop.Tupp)
            temp_min = temp_min.clip(upper=crop.Tupp)
            Tmean = (temp_max + temp_min) / 2
            Tmean = Tmean.clip(lower=crop.Tbase)
            gdd = Tmean - crop.Tbase

        gdd_cum = np.cumsum(gdd).reset_index(drop=True)

        assert (
            gdd_cum.values[-1] > crop.Maturity
        ), f"not enough growing degree days in simulation ({gdd_cum.values[-1]}) to reach maturity ({crop.Maturity})"

        crop.MaturityCD = (gdd_cum > crop.Maturity).idxmax() + 1

        assert crop.MaturityCD < 365, "crop will take longer than 1 year to mature"

        # 1. gdd's from sowing to maximum canopy cover
        crop.MaxCanopyCD = (gdd_cum > crop.MaxCanopy).idxmax() + 1
        # 2. gdd's from sowing to end of vegetative growth
        crop.CanopyDevEndCD = (gdd_cum > crop.CanopyDevEnd).idxmax() + 1
        # 3. Calendar days from sowing to start of yield_ formation
        crop.HIstartCD = (gdd_cum > crop.HIstart).idxmax() + 1
        # 4. Calendar days from sowing to end of yield_ formation
        crop.HIendCD = (gdd_cum > crop.HIend).idxmax() + 1
        # 5. Duration of yield_ formation in calendar days
        crop.YldFormCD = crop.HIendCD - crop.HIstartCD
        if crop.CropType == 3:
            # 1. Calendar days from sowing to end of flowering
            FloweringEnd = (gdd_cum > crop.FloweringEnd).idxmax() + 1
            # 2. Duration of flowering in calendar days
            crop.FloweringCD = FloweringEnd - crop.HIstartCD
        else:
            crop.FloweringCD = ModelConstants.NO_VALUE

    return crop

aquacrop.initialize.compute_variables

compute_variables(param_struct, weather_df, clock_struct, acfp=dirname(dirname(abspath(__file__))))

Function to compute additional variables needed to run the model eg. CO2 Creates cropstruct jit class objects

Arguments:

param_struct (ParamStruct):  Contains model paramaters

weather_df (DataFrame):  weather data

clock_struct (ClockStruct):  time params

acfp (Path):  path to aquacrop directory containing co2 data

Returns:

param_struct (ParamStruct):  updated model params

Source code in aquacrop/initialize/compute_variables.py

def compute_variables(
    param_struct: "ParamStruct",
    weather_df: "DataFrame",
    clock_struct: "ClockStruct",
    acfp: str = dirname(dirname(abspath(__file__))),
) -> "ParamStruct":
    """
    Function to compute additional variables needed to run the model eg. CO2
    Creates cropstruct jit class objects

    Arguments:

        param_struct (ParamStruct):  Contains model paramaters

        weather_df (DataFrame):  weather data

        clock_struct (ClockStruct):  time params

        acfp (Path):  path to aquacrop directory containing co2 data

    Returns:

        param_struct (ParamStruct):  updated model params


    """

    if param_struct.water_table == 1:

        param_struct.Soil.add_capillary_rise_params()

    # Calculate readily evaporable water in surface layer
    if param_struct.Soil.adj_rew == 0:
        param_struct.Soil.rew = round(
            (
                1000
                * (
                    param_struct.Soil.profile.th_fc.iloc[0]
                    - param_struct.Soil.profile.th_dry.iloc[0]
                )
                * param_struct.Soil.evap_z_surf
            ),
            2,
        )

    if param_struct.Soil.calc_cn == 1:
        # adjust curve number
        ksat = param_struct.Soil.profile.Ksat.iloc[0]
        if ksat > 864:
            param_struct.Soil.cn = 46
        elif ksat > 347:
            param_struct.Soil.cn = 61
        elif ksat > 36:
            param_struct.Soil.cn = 72
        elif ksat > 0:
            param_struct.Soil.cn = 77

        assert ksat > 0

    for i in range(param_struct.NCrops):

        crop = param_struct.CropList[i]
        # crop.calculate_additional_params()

        # Crop calander
        crop = compute_crop_calendar(
            crop,
            clock_struct.planting_dates,
            clock_struct.simulation_start_date,
            clock_struct.time_span,
            weather_df,
        )

        # Harvest index param_struct.Seasonal_Crop_List[clock_struct.season_counter].Paramsgrowth coefficient
        crop.HIGC = calculate_HIGC(
            crop.YldFormCD,
            crop.HI0,
            crop.HIini,
        )

        # Days to linear harvest_index switch point
        if crop.CropType == 3:
            # Determine linear switch point and HIGC rate for fruit/grain crops
            crop.tLinSwitch, crop.dHILinear = calculate_HI_linear(
                crop.YldFormCD, crop.HIini, crop.HI0, crop.HIGC
            )
        else:
            # No linear switch for leafy vegetable or root/tiber crops
            crop.tLinSwitch = 0
            crop.dHILinear = 0.0

        param_struct.CropList[i] = crop

    # Calculate WP adjustment factor for elevation in CO2 concentration
    # Load CO2 data
    co2Data = param_struct.CO2.co2_data

    # Years
    start_year, end_year = pd.DatetimeIndex(
        [clock_struct.simulation_start_date, clock_struct.simulation_end_date]
    ).year
    sim_years = np.arange(start_year, end_year + 1)

    # Interpolate data
    CO2conc_interp = np.interp(sim_years, co2Data.year, co2Data.ppm)

    # Store data
    param_struct.CO2.co2_data_processed = pd.Series(CO2conc_interp, index=sim_years)  # maybe get rid of this

    # Get CO2 concentration for first year
    CO2conc = param_struct.CO2.co2_data_processed.iloc[0]

    # param_struct.CO2 = param_struct.co2_concentration_adj

    # if user specified constant concentration
    if  param_struct.CO2.constant_conc is True:
        if param_struct.CO2.current_concentration > 0.:
            CO2conc = param_struct.CO2.current_concentration
        else:
            CO2conc = param_struct.CO2.co2_data_processed.iloc[0]

    param_struct.CO2.current_concentration = CO2conc

    CO2ref = param_struct.CO2.ref_concentration

    # Get CO2 weighting factor for first year
    if CO2conc <= CO2ref:
        fw = 0
    else:
        if CO2conc >= 550:
            fw = 1
        else:
            fw = 1 - ((550 - CO2conc) / (550 - CO2ref))

    # Determine adjustment for each crop in first year of simulation
    for i in range(param_struct.NCrops):
        crop = param_struct.CropList[i]
        # Determine initial adjustment
        fCO2old = (CO2conc / CO2ref) / (
            1
            + (CO2conc - CO2ref)
            * (
                (1 - fw) * crop.bsted
                + fw * ((crop.bsted * crop.fsink) + (crop.bface * (1 - crop.fsink)))
            )
        )
        # New adjusted correction coefficient for CO2 (version 7 of AquaCrop)
    if (CO2conc > CO2ref):
        # Calculate shape factor
        fshape = -4.61824 - 3.43831*crop.fsink - 5.32587*crop.fsink*crop.fsink
        # Determine adjustment for CO2
        if (CO2conc >= 2000):
            fCO2new = 1.58  # Maximum CO2 adjustment 
        else:
            CO2rel = (CO2conc-CO2ref)/(2000-CO2ref)
            fCO2new = 1 + 0.58 * ((np.exp(CO2rel*fshape) - 1)/(np.exp(fshape) - 1))


    # Select adjusted coefficient for CO2
    if (CO2conc <= CO2ref):
        fCO2 = fCO2old
    elif ((CO2conc <= 550) and (fCO2old < fCO2new)):
        fCO2 = fCO2old
    else:
        fCO2 = fCO2new

        # Consider crop type
    if crop.WP >= 40:
        # No correction for C4 crops
        ftype = 0
    elif crop.WP <= 20:
        # Full correction for C3 crops
        ftype = 1
    else:
        ftype = (40 - crop.WP) / (40 - 20)

        # Total adjustment
    crop.fCO2 = 1 + ftype * (fCO2 - 1)

    param_struct.CropList[i] = crop


    # change this later
    if param_struct.NCrops == 1:
        crop_list = [
            deepcopy(param_struct.CropList[0])
            for i in range(len(param_struct.CropChoices))
        ]
        # param_struct.Seasonal_Crop_List = [deepcopy(param_struct.CropList[0]) for i in range(len(param_struct.CropChoices))]

    else:
        crop_list = param_struct.CropList

    # add crop for out of growing season
    # param_struct.Fallow_Crop = deepcopy(param_struct.Seasonal_Crop_List[0])
    Fallow_Crop = deepcopy(crop_list[0])

    param_struct.Seasonal_Crop_List = []
    for crop in crop_list:
        crop_struct = CropStruct()
        for a, v in crop.__dict__.items():
            if hasattr(crop_struct, a):
                crop_struct.__setattr__(a, v)

        param_struct.Seasonal_Crop_List.append(crop_struct)

    fallow_struct = CropStruct()
    for a, v in Fallow_Crop.__dict__.items():
        if hasattr(fallow_struct, a):
            fallow_struct.__setattr__(a, v)

    param_struct.Fallow_Crop = fallow_struct

    return param_struct

aquacrop.initialize.create_soil_profile

create_soil_profile(param_struct)

funciton to create soil profile namedTuple to store soil info. Its much faster to access the info when its in a namedTuple compared to a dataframe

Arguments:

param_struct (ParamStruct):  Contains model crop and soil paramaters

Returns:

param_struct (ParamStruct):  updated with soil profile

Source code in aquacrop/initialize/create_soil_profile.py

def create_soil_profile(param_struct: "ParamStruct") -> "ParamStruct":
    """
    funciton to create soil profile namedTuple to store soil info.
    Its much faster to access the info when its in a namedTuple
    compared to a dataframe

    Arguments:

        param_struct (ParamStruct):  Contains model crop and soil paramaters

    Returns:

        param_struct (ParamStruct):  updated with soil profile


    """

    profile = SoilProfile(int(param_struct.Soil.profile.shape[0]))

    pdf = param_struct.Soil.profile.astype("float64")

    profile.dz = pdf.dz.values
    profile.dzsum = pdf.dzsum.values
    profile.zBot = pdf.zBot.values
    profile.z_top = pdf.z_top.values
    profile.zMid = pdf.zMid.values

    profile.Comp = np.int64(pdf.Comp.values)
    profile.Layer = np.int64(pdf.Layer.values)
    # profile.Layer_dz = pdf.Layer_dz.values
    profile.th_wp = pdf.th_wp.values
    profile.th_fc = pdf.th_fc.values
    profile.th_s = pdf.th_s.values

    profile.Ksat = pdf.Ksat.values
    profile.Penetrability = pdf.penetrability.values
    profile.th_dry = pdf.th_dry.values
    profile.tau = pdf.tau.values
    profile.th_fc_Adj = pdf.th_fc_Adj.values

    if param_struct.water_table == 1:
        profile.aCR = pdf.aCR.values
        profile.bCR = pdf.bCR.values
    else:
        profile.aCR = pdf.dz.values * 0.0
        profile.bCR = pdf.dz.values * 0.0

    # param_struct.Soil.profile = profile

    param_struct.Soil.Profile = SoilProfileNT(
        dz=profile.dz,
        dzsum=profile.dzsum,
        zBot=profile.zBot,
        z_top=profile.z_top,
        zMid=profile.zMid,
        Comp=profile.Comp,
        Layer=profile.Layer,
        th_wp=profile.th_wp,
        th_fc=profile.th_fc,
        th_s=profile.th_s,
        Ksat=profile.Ksat,
        Penetrability=profile.Penetrability,
        th_dry=profile.th_dry,
        tau=profile.tau,
        th_fc_Adj=profile.th_fc_Adj,
        aCR=profile.aCR,
        bCR=profile.bCR,
    )

    return param_struct

aquacrop.initialize.read_clocks_parameters

Inititalize clocks parameters

check_max_simulation_days(sim_start_time, sim_end_time)

Check that the date range of the simulation is less than 580 years. In pandas this cannot happen due to the size of the variable

Arguments:

sim_start_time (str): simulation start date YYYY/MM/DD

sim_end_time (str): simulation start date YYYY/MM/DD

Source code in aquacrop/initialize/read_clocks_parameters.py

def check_max_simulation_days(
    sim_start_time: str,
    sim_end_time: str):
    """
    Check that the date range of the simulation is less than 580 years.
    In pandas this cannot happen due to the size of the variable

    Arguments:

        sim_start_time (str): simulation start date YYYY/MM/DD

        sim_end_time (str): simulation start date YYYY/MM/DD

    """
    start_year = int(sim_start_time.split("/")[0])
    end_year = int(sim_end_time.split("/")[0])
    if (end_year - start_year) > 580:
        raise ValueError("Simulation period must be less than 580 years.")

read_clock_parameters(sim_start_time, sim_end_time, off_season=False)

Function to read in start and end simulation time and return a ClockStruct object

Arguments:

sim_start_time (str): simulation start date

sim_end_time (str): simulation start date

off_season (bool): True, simulate off season
                  False, skip ahead to next season post-harvest

Returns:

clock_struct (ClockStruct): simulation time paramaters

Source code in aquacrop/initialize/read_clocks_parameters.py

def read_clock_parameters(
    sim_start_time: str,
    sim_end_time: str,
    off_season: bool=False) -> ClockStruct:
    """
    Function to read in start and end simulation time and return a ClockStruct object

    Arguments:

        sim_start_time (str): simulation start date

        sim_end_time (str): simulation start date

        off_season (bool): True, simulate off season
                          False, skip ahead to next season post-harvest

    Returns:

        clock_struct (ClockStruct): simulation time paramaters


    """
    check_max_simulation_days(sim_start_time, sim_end_time)

    # Extract data and put into pandas datetime format
    pandas_sim_start_time = pd.to_datetime(sim_start_time)
    pandas_sim_end_time = pd.to_datetime(sim_end_time)

    # create ClockStruct object
    clock_struct = ClockStruct()

    # Add variables
    clock_struct.simulation_start_date = pandas_sim_start_time
    clock_struct.simulation_end_date = pandas_sim_end_time

    clock_struct.n_steps = (pandas_sim_end_time - pandas_sim_start_time).days + 1
    clock_struct.time_span = pd.date_range(
        freq="D", start=pandas_sim_start_time, end=pandas_sim_end_time
    )

    clock_struct.step_start_time = clock_struct.time_span[0]
    clock_struct.step_end_time = clock_struct.time_span[1]

    clock_struct.sim_off_season = off_season

    return clock_struct

aquacrop.initialize.read_field_managment

read_field_management(ParamStruct, FieldMngt, FallowFieldMngt)

store field management variables as FieldMngtStruct object

Arguments:

ParamStruct (ParamStruct):  Contains model crop and soil paramaters

FieldMngt (FieldMngt):  field mngt params

FallowFieldMngt (FieldMngt): fallow field mngt params

Returns:

ParamStruct (ParamStruct):  updated ParamStruct with field management info

Source code in aquacrop/initialize/read_field_managment.py

def read_field_management(
    ParamStruct: "ParamStruct",
    FieldMngt: "FieldMngt",
    FallowFieldMngt: "FieldMngt") -> "ParamStruct":

    """
    store field management variables as FieldMngtStruct object

    Arguments:

        ParamStruct (ParamStruct):  Contains model crop and soil paramaters

        FieldMngt (FieldMngt):  field mngt params

        FallowFieldMngt (FieldMngt): fallow field mngt params

    Returns:

        ParamStruct (ParamStruct):  updated ParamStruct with field management info


    """

    field_mngt_struct = FieldMngtStruct()
    for a, v in FieldMngt.__dict__.items():
        if hasattr(field_mngt_struct, a):
            field_mngt_struct.__setattr__(a, v)

    fallow_field_mngt_struct = FieldMngtStruct()
    for a, v in FallowFieldMngt.__dict__.items():
        if hasattr(fallow_field_mngt_struct, a):
            fallow_field_mngt_struct.__setattr__(a, v)

    ParamStruct.FieldMngt = field_mngt_struct
    ParamStruct.FallowFieldMngt = fallow_field_mngt_struct

    return ParamStruct

aquacrop.initialize.read_groundwater_table

read_groundwater_table(ParamStruct, GwStruct, ClockStruct)

Function to initialise groundwater parameters

Arguments:

ParamStruct (ParamStruct): Contains model paramaters

GwStruct (GroundWater): groundwater params

ClockStruct (ClockStruct): time params

Returns:

ParamStruct (ParamStruct): updated with GW info

Source code in aquacrop/initialize/read_groundwater_table.py

def read_groundwater_table(
    ParamStruct: "ParamStruct",
    GwStruct: "GroundWater",
    ClockStruct: "ClockStruct") -> "ParamStruct":
    """
    Function to initialise groundwater parameters

    Arguments:

        ParamStruct (ParamStruct): Contains model paramaters

        GwStruct (GroundWater): groundwater params

        ClockStruct (ClockStruct): time params

    Returns:

        ParamStruct (ParamStruct): updated with GW info

    """

    # assign water table value and method
    WT = GwStruct.water_table
    WTMethod = GwStruct.method

    # check if water table present
    if WT == "N":
        ParamStruct.water_table = 0
        ParamStruct.z_gw = 999 * np.ones(len(ClockStruct.time_span))
        ParamStruct.zGW_dates = ClockStruct.time_span
        ParamStruct.WTMethod = "None"
    elif WT == "Y":
        ParamStruct.water_table = 1

        df = pd.DataFrame([GwStruct.dates, GwStruct.values]).T
        df.columns = ["Date", "Depth(mm)"]

        # get date in correct format
        df.Date = pd.DatetimeIndex(df.Date)
        # print(f'DF length: {len(df)}')
        # print(f'Index length: {len(df.index)}')

        if len(df) == 1:

            # if only 1 watertable depth then set that value to be constant
            # accross whole simulation            
            z_gw = pd.DataFrame(
                data=df["Depth(mm)"].iloc[0]*np.ones(len(ClockStruct.time_span)),
                index=pd.to_datetime(ClockStruct.time_span),
                columns=['Depth(mm)']
            )['Depth(mm)']

        elif len(df) > 1:
            # check water table method
            if WTMethod == "Constant":

                # No interpolation between dates

                # create daily depths for each simulation day
                z_gw = pd.Series(
                    np.nan * np.ones(len(ClockStruct.time_span)), index=ClockStruct.time_span
                )

                # assign constant depth for all dates in between
                for row in range(len(df)):
                    date = df.Date.iloc[row]
                    depth = df["Depth(mm)"].iloc[row]
                    z_gw.loc[z_gw.index >= date] = depth
                    if row == 0:
                        z_gw.loc[z_gw.index <= date] = depth

            elif WTMethod == "Variable":

                # Linear interpolation between dates

                # create daily depths for each simulation day
                # fill unspecified days with NaN
                z_gw = pd.Series(
                    np.nan * np.ones(len(ClockStruct.time_span)), index=ClockStruct.time_span
                )

                for row in range(len(df)):
                    date = df.Date.iloc[row]
                    depth = df["Depth(mm)"].iloc[row]
                    z_gw.loc[date] = depth

                # Interpolate daily groundwater depths
                z_gw = z_gw.interpolate()

        # assign values to Paramstruct object
        ParamStruct.z_gw = z_gw.values
        ParamStruct.zGW_dates = z_gw.index.values
        ParamStruct.WTMethod = WTMethod

    return ParamStruct

aquacrop.initialize.read_irrigation_management

read_irrigation_management(ParamStruct, IrrMngt, ClockStruct)

initilize irrigation management and store as IrrMngtStruct object

Arguments:

ParamStruct (ParamStruct):  Contains model crop and soil paramaters

IrrMngt (IrrigationManagement):  irr mngt params object

ClockStruct (ClockStruct):  time paramaters

Returns:

ParamStruct (ParamStruct):  updated model paramaters

Source code in aquacrop/initialize/read_irrigation_management.py

def read_irrigation_management(
    ParamStruct: "ParamStruct",
    IrrMngt: "IrrigationManagement",
    ClockStruct: "ClockStruct") -> "ParamStruct":
    """
    initilize irrigation management and store as IrrMngtStruct object

    Arguments:

        ParamStruct (ParamStruct):  Contains model crop and soil paramaters

        IrrMngt (IrrigationManagement):  irr mngt params object

        ClockStruct (ClockStruct):  time paramaters


    Returns:

        ParamStruct (ParamStruct):  updated model paramaters



    """
    # If specified, read input irrigation time-series
    if IrrMngt.irrigation_method == 3:

        df = IrrMngt.Schedule.copy()
        # change the index to the date
        df.index = pd.DatetimeIndex(df.Date)

        try:
            # create a dateframe containing the daily irrigation to
            # be applied for every day in the simulation
            df = df.reindex(ClockStruct.time_span, fill_value=0).drop("Date", axis=1)

            IrrMngt.Schedule = np.array(df.values, dtype=float).flatten()

        except TypeError:
            # older version of pandas with not reindex

            # create new dataframe for whole simulation
            # populate new dataframe with old values
            new_df = pd.DataFrame(data=np.zeros(len(ClockStruct.time_span)),
                index=pd.to_datetime(ClockStruct.time_span),
                columns=['Depth']
                )

            # fill in the new dataframe with irrigation schedule
            new_df.loc[df.index]=df.Depth.values

            IrrMngt.Schedule = np.array(new_df.values, dtype=float).flatten()

    else:

        IrrMngt.Schedule = np.zeros(len(ClockStruct.time_span))

    IrrMngt.SMT = np.array(IrrMngt.SMT, dtype=float)

    irr_mngt_struct = IrrMngtStruct(len(ClockStruct.time_span))
    for a, v in IrrMngt.__dict__.items():
        if hasattr(irr_mngt_struct, a):
            irr_mngt_struct.__setattr__(a, v)

    ParamStruct.IrrMngt = irr_mngt_struct
    ParamStruct.FallowIrrMngt = IrrMngtStruct(len(ClockStruct.time_span))

    return ParamStruct

aquacrop.initialize.read_model_initial_conditions

read_model_initial_conditions(ParamStruct, ClockStruct, InitWC, crop)

Function to set up initial model conditions

Arguments:

ParamStruct (ParamStruct):  Contains model paramaters

ClockStruct (ClockStruct):  time paramaters

InitWC (InitialWaterContent):  initial water content

crop (Crop): crop parameters

Returns:

ParamStruct (ParamStruct):  updated ParamStruct object

InitCond (InitialCondition):  containing initial model conditions/counters

Source code in aquacrop/initialize/read_model_initial_conditions.py

def read_model_initial_conditions(
    ParamStruct: "ParamStruct",
    ClockStruct: "ClockStruct",
    InitWC: "InitialWaterContent",
    crop: "Crop") -> Tuple["ParamStruct", "InitialCondition"]:
    """
    Function to set up initial model conditions

    Arguments:

        ParamStruct (ParamStruct):  Contains model paramaters

        ClockStruct (ClockStruct):  time paramaters

        InitWC (InitialWaterContent):  initial water content

        crop (Crop): crop parameters


    Returns:

        ParamStruct (ParamStruct):  updated ParamStruct object

        InitCond (InitialCondition):  containing initial model conditions/counters

    """

    ###################
    # creat initial condition class
    ###################

    InitCond = InitialCondition(len(ParamStruct.Soil.profile))

    # class_args = {key:value for key, value in InitCond_class.__dict__.items() if not key.startswith('__') and not callable(key)}
    # InitCond = InitCondStruct(**class_args)

    if ClockStruct.season_counter == -1:
        InitCond.z_root = 0.
        InitCond.cc0_adj = 0.

    elif ClockStruct.season_counter == 0:
        InitCond.z_root = ParamStruct.Seasonal_Crop_List[0].Zmin
        InitCond.cc0_adj = ParamStruct.Seasonal_Crop_List[0].CC0

    # Set HIfinal to crop's reference harvest index
    InitCond.HIfinal = crop.HI0

    ##################
    # save field management
    ##################

    # Initial surface storage between any soil bunds
    if ClockStruct.season_counter == -1:
        # First day of simulation is in fallow period
        if (ParamStruct.FallowFieldMngt.bunds) and (
            float(ParamStruct.FallowFieldMngt.z_bund) > 0.001
        ):
            # Get initial storage between surface bunds
            InitCond.surface_storage = float(ParamStruct.FallowFieldMngt.bund_water)
            if InitCond.surface_storage > float(ParamStruct.FallowFieldMngt.z_bund):
                InitCond.surface_storage = float(ParamStruct.FallowFieldMngt.z_bund)
        else:
            # No surface bunds
            InitCond.surface_storage = 0

    elif ClockStruct.season_counter == 0:
        # First day of simulation is in first growing season
        # Get relevant field management structure parameters
        FieldMngtTmp = ParamStruct.FieldMngt
        if (FieldMngtTmp.bunds) and (float(FieldMngtTmp.z_bund) > 0.001):
            # Get initial storage between surface bunds
            InitCond.surface_storage = float(FieldMngtTmp.bund_water)
            if InitCond.surface_storage > float(FieldMngtTmp.z_bund):
                InitCond.surface_storage = float(FieldMngtTmp.z_bund)
        else:
            # No surface bunds
            InitCond.surface_storage = 0

    ############
    # watertable
    ############

    profile = ParamStruct.Soil.profile

    # Check for presence of groundwater table
    if ParamStruct.water_table == 0:  # No water table present
        # Set initial groundwater level to dummy value
        InitCond.z_gw = ModelConstants.NO_VALUE
        InitCond.wt_in_soil = False
        # Set adjusted field capacity to default field capacity
        InitCond.th_fc_Adj = profile.th_fc.values
    elif ParamStruct.water_table == 1:  # Water table is present
        # Set initial groundwater level
        InitCond.z_gw = float(ParamStruct.z_gw[ClockStruct.time_step_counter])
        # Find compartment mid-points
        zMid = profile.zMid
        # Check if water table is within modelled soil profile
        if InitCond.z_gw >= 0:
            idx = zMid[zMid >= InitCond.z_gw].index
            if idx.shape[0] == 0:
                InitCond.wt_in_soil = False
            else:
                InitCond.wt_in_soil = True
        else:
            InitCond.wt_in_soil = False

        # Adjust compartment field capacity
        compi = int(len(profile)) - 1
        thfcAdj = np.zeros(compi + 1)
        while compi >= 0:
            # get soil layer of compartment
            compdf = profile.loc[compi]
            if compdf.th_fc <= 0.1:
                Xmax = 1
            else:
                if compdf.th_fc >= 0.3:
                    Xmax = 2
                else:
                    pF = 2 + 0.3 * (compdf.th_fc - 0.1) / 0.2
                    Xmax = (np.exp(pF * np.log(10))) / 100

            if (InitCond.z_gw < 0) or ((InitCond.z_gw - zMid.iloc[compi]) >= Xmax):
                for ii in range(compi+1):
                    compdfii = profile.loc[ii]
                    thfcAdj[ii] = compdfii.th_fc

                compi = -1
            else:
                if compdf.th_fc >= compdf.th_s:
                    thfcAdj[compi] = compdf.th_fc
                else:
                    if zMid.iloc[compi] >= InitCond.z_gw:
                        thfcAdj[compi] = compdf.th_s
                    else:
                        dV = compdf.th_s - compdf.th_fc
                        dFC = (dV / (Xmax ** 2)) * ((zMid.iloc[compi] - (InitCond.z_gw - Xmax)) ** 2)
                        thfcAdj[compi] = compdf.th_fc + dFC

                compi = compi - 1

        # Store adjusted field capacity values
        InitCond.th_fc_Adj = np.round(thfcAdj, 3)

    profile["th_fc_Adj"] = np.round(InitCond.th_fc_Adj, 3)

    # create hydrology df to group by layer instead of compartment
    ParamStruct.Soil.Hydrology = profile.groupby("Layer").mean().drop(["dz", "dzsum"], axis=1)
    ParamStruct.Soil.Hydrology["dz"] = profile.groupby("Layer").sum().dz

    ###################
    # initial water contents
    ###################

    typestr = InitWC.wc_type
    methodstr = InitWC.method

    depth_layer = InitWC.depth_layer
    datapoints = InitWC.value

    values = np.zeros(len(datapoints))

    hydf = ParamStruct.Soil.Hydrology

    # Assign data
    if typestr == "Num":
        # Values are defined as numbers (m3/m3) so no calculation required
        depth_layer = np.array(depth_layer, dtype=float)
        values = np.array(datapoints, dtype=float)

    elif typestr == "Pct":
        # Values are defined as percentage of taw. Extract and assign value for
        # each soil layer based on calculated/input soil hydraulic properties
        depth_layer = np.array(depth_layer, dtype=float)
        datapoints = np.array(datapoints, dtype=float)

        for ii in range(len(values)):
            if methodstr == "Depth":
                depth = depth_layer[ii]
                value = datapoints[ii]

                # Find layer at specified depth
                if depth < profile.dzsum.iloc[-1]:
                    layer = profile.query(f"{depth}<dzsum").Layer.iloc[0]
                else:
                    layer = profile.Layer.iloc[-1]

                compdf = hydf.loc[layer]

                # Calculate moisture content at specified depth
                values[ii] = compdf.th_wp + ((value / 100) * (compdf.th_fc - compdf.th_wp))
            elif methodstr == "Layer":
                # Calculate moisture content at specified layer
                layer = depth_layer[ii]
                value = datapoints[ii]

                compdf = hydf.loc[layer]

                values[ii] = compdf.th_wp + ((value / 100) * (compdf.th_fc - compdf.th_wp))

    elif typestr == "Prop":
        # Values are specified as soil hydraulic properties (SAT, FC, or WP).
        # Extract and assign value for each soil layer
        depth_layer = np.array(depth_layer, dtype=float)
        datapoints = np.array(datapoints, dtype=str)

        for ii in range(len(values)):
            if methodstr == "Depth":
                # Find layer at specified depth
                depth = depth_layer[ii]
                value = datapoints[ii]

                # Find layer at specified depth
                if depth < profile.dzsum.iloc[-1]:
                    layer = profile.query(f"{depth}<dzsum").Layer.iloc[0]
                else:
                    layer = profile.Layer.iloc[-1]

                compdf = hydf.loc[layer]

                # Calculate moisture content at specified depth
                if value == "SAT":
                    values[ii] = compdf.th_s
                if value == "FC":
                    values[ii] = compdf.th_fc
                if value == "WP":
                    values[ii] = compdf.th_wp

            elif methodstr == "Layer":
                # Calculate moisture content at specified layer
                layer = depth_layer[ii]
                value = datapoints[ii]

                compdf = hydf.loc[layer]

                if value == "SAT":
                    values[ii] = compdf.th_s
                if value == "FC":
                    values[ii] = compdf.th_fc
                if value == "WP":
                    values[ii] = compdf.th_wp

    # Interpolate values to all soil compartments

    thini = np.zeros(int(profile.shape[0]))
    if methodstr == "Layer":
        for ii in range(len(values)):
            layer = depth_layer[ii]
            value = values[ii]

            idx = profile.query(f"Layer=={int(layer)}").index

            thini[idx] = value

        InitCond.th = thini

    elif methodstr == "Depth":
        depths = depth_layer

        # Add zero point
        if depths[0] > 0:
            depths = np.append([0], depths)
            values = np.append([values[0]], values)

        # Add end point (bottom of soil profile)
        if depths[-1] < ParamStruct.Soil.zSoil:
            depths = np.append(depths, [ParamStruct.Soil.zSoil])
            values = np.append(values, [values[-1]])

        # Find centroids of compartments
        SoilDepths = profile.dzsum.values
        comp_top = np.append([0], SoilDepths[:-1])
        comp_bot = SoilDepths
        comp_mid = (comp_top + comp_bot) / 2
        # Interpolate initial water contents to each compartment
        thini = np.interp(comp_mid, depths, values)
        InitCond.th = thini

    # If groundwater table is present and calculating water contents based on
    # field capacity, then reset value to account for possible changes in field
    # capacity caused by capillary rise effects
    if ParamStruct.water_table == 1:
        if (typestr == "Prop") and (datapoints[-1] == "FC"):
            InitCond.th = InitCond.th_fc_Adj

    # If groundwater table is present in soil profile then set all water
    # contents below the water table to saturation
    if InitCond.wt_in_soil is True:
        # Find compartment mid-points
        SoilDepths = profile.dzsum.values
        comp_top = np.append([0], SoilDepths[:-1])
        comp_bot = SoilDepths
        comp_mid = (comp_top + comp_bot) / 2
        idx = np.where(comp_mid >= InitCond.z_gw)[0][0]
        for ii in range(idx, len(profile)):
            layeri = profile.loc[ii].Layer
            InitCond.th[ii] = hydf.th_s.loc[layeri]

    InitCond.thini = InitCond.th

    ParamStruct.Soil.profile = profile
    ParamStruct.Soil.Hydrology = hydf

    return ParamStruct, InitCond

aquacrop.initialize.read_model_parameters

read_model_parameters(clock_struct, soil, crop, weather_df)

Finalise soil and crop paramaters including planting and harvest dates save to new object param_struct

Arguments:

clock_struct (ClockStruct):  time params

soil (Soil):  soil object

crop (Crop):  crop object

weather_df (DataFrame): list of datetimes

Returns:

clock_struct (ClockStruct): updated time paramaters

param_struct (ParamStruct):  Contains model crop and soil paramaters

Source code in aquacrop/initialize/read_model_parameters.py

def read_model_parameters(
    clock_struct: "ClockStruct",
    soil: "Soil",
    crop: "Crop",
    weather_df: "DataFrame"):

    """
    Finalise soil and crop paramaters including planting and harvest dates
    save to new object param_struct


    Arguments:

        clock_struct (ClockStruct):  time params

        soil (Soil):  soil object

        crop (Crop):  crop object

        weather_df (DataFrame): list of datetimes

    Returns:

        clock_struct (ClockStruct): updated time paramaters

        param_struct (ParamStruct):  Contains model crop and soil paramaters

    """
    # create param_struct object
    param_struct = ParamStruct()

    soil.fill_nan()

    # Assign soil object to param_struct
    param_struct.Soil = soil

    while soil.zSoil < crop.Zmax + 0.1:
        for i in soil.profile.index[::-1]:
            if soil.profile.loc[i, "dz"] < 0.25:
                soil.profile.loc[i, "dz"] += 0.1
                soil.fill_nan()
                break

    # TODO: Why all these commented lines? The model does not allow rotations now?
    ###########
    # crop
    ###########

    #     if isinstance(crop, Iterable):
    #         cropList=list(crop)
    #     else:
    #         cropList = [crop]

    #     # assign variables to paramstruct
    #     paramStruct.nCrops = len(cropList)
    #     if paramStruct.nCrops > 1:
    #         paramStruct.SpecifiedPlantcalendar = 'yield_'
    #     else:
    #         paramStruct.SpecifiedPlantcalendar = 'N'

    #     # add crop list to paramStruct
    #     paramStruct.cropList = cropList

    ############################
    # plant and harvest times
    ############################

    #     # find planting and harvest dates
    #     # check if there is more than 1 crop or multiple plant dates in sim year
    #     if paramStruct.SpecifiedPlantcalendar == "yield_":
    #         # if here than crop rotation occours during same period

    #         # create variables from dataframe
    #         plantingDates = pd.to_datetime(planting_dates)
    #         harvestDates = pd.to_datetime(harvest_dates)

    #         if (paramStruct.nCrops > 1):

    #             cropChoices = [crop.name for crop in paramStruct.cropList]

    #         assert len(cropChoices) == len(plantingDates) == len(harvestDates)

    # elif paramStruct.nCrops == 1:
    # Only one crop type considered during simulation - i.e. no rotations
    # either within or between years
    crop_list = [crop]
    param_struct.CropList = crop_list
    param_struct.NCrops = 1

    # Get start and end years for full simulation
    sim_start_date = clock_struct.simulation_start_date
    sim_end_date = clock_struct.simulation_end_date

    if crop.harvest_date is None:
        crop = compute_crop_calendar(
            crop,
            clock_struct.planting_dates,
            clock_struct.simulation_start_date,
            clock_struct.time_span,
            weather_df,
        )
        mature = int(crop.MaturityCD + 30)
        plant = pd.to_datetime("1990/" + crop.planting_date)
        harv = plant + np.timedelta64(mature, "D")
        new_harvest_date = str(harv.month) + "/" + str(harv.day)
        crop.harvest_date = new_harvest_date

    # extract years from simulation start and end date
    start_end_years = [sim_start_date.year, sim_end_date.year]

    # check if crop growing season runs over calander year
    # Planting and harvest dates are in days/months format so just add arbitrary year
    single_year = pd.to_datetime("1990/" + crop.planting_date) < pd.to_datetime(
        "1990/" + crop.harvest_date
    )

    if single_year:
        # if normal year

        # Check if the simulation in the following year does not exceed planting date.
        mock_simulation_end_date = pd.to_datetime("1990/" + f'{sim_end_date.month}' + "/" + f'{sim_end_date.day}')
        mock_simulation_start_date = pd.to_datetime("1990/" + crop.planting_date)
        last_simulation_year_does_not_start = mock_simulation_end_date <= mock_simulation_start_date

        if last_simulation_year_does_not_start:
            start_end_years[1] = start_end_years[1] - 1

        # specify the planting and harvest years as normal
        plant_years = list(range(start_end_years[0], start_end_years[1] + 1))
        harvest_years = plant_years
    else:
        # if it takes over a year then the plant year finishes 1 year before end of sim
        # and harvest year starts 1 year after sim start

        if (
            pd.to_datetime(str(start_end_years[1] + 2) + "/" + crop.harvest_date)
            < sim_end_date
        ):

            # specify shifted planting and harvest years
            plant_years = list(range(start_end_years[0], start_end_years[1] + 1))
            harvest_years = list(range(start_end_years[0] + 1, start_end_years[1] + 2))
        else:

            plant_years = list(range(start_end_years[0], start_end_years[1]))
            harvest_years = list(range(start_end_years[0] + 1, start_end_years[1] + 1))

    # Correct for partial first growing season (may occur when simulating
    # off-season soil water balance)
    if (
        pd.to_datetime(str(plant_years[0]) + "/" + crop.planting_date)
        < clock_struct.simulation_start_date
    ):
        # shift everything by 1 year
        plant_years = plant_years[1:]
        harvest_years = harvest_years[1:]

    # ensure number of planting and harvest years are the same
    assert len(plant_years) == len(harvest_years)

    # create lists to hold variables
    planting_dates = []
    harvest_dates = []
    crop_choices = []

    # save full harvest/planting dates and crop choices to lists
    for i, _ in enumerate(plant_years):
        planting_dates.append(
            str(plant_years[i]) + "/" + param_struct.CropList[0].planting_date
        )
        harvest_dates.append(
            str(harvest_years[i]) + "/" + param_struct.CropList[0].harvest_date
        )
        crop_choices.append(param_struct.CropList[0].Name)

    # save crop choices
    param_struct.CropChoices = list(crop_choices)

    # save clock paramaters
    clock_struct.planting_dates = pd.to_datetime(planting_dates)
    clock_struct.harvest_dates = pd.to_datetime(harvest_dates)
    clock_struct.n_seasons = len(planting_dates)

    # Initialise growing season counter
    if pd.to_datetime(clock_struct.step_start_time) == clock_struct.planting_dates[0]:
        clock_struct.season_counter = 0
    else:
        clock_struct.season_counter = -1

    # return the FileLocations object as i have added some elements
    return clock_struct, param_struct

aquacrop.initialize.read_weather_inputs

Initialize weather data

read_weather_inputs(clock_sctruct, weather_df)

Clip weather to start and end simulation dates

Arguments:

clock_sctruct (ClockStruct): ClockStruct object

weather_df (DataFrame): weather dataframe

Returns:

weather_df (DataFrame): clipped weather dataframe

Source code in aquacrop/initialize/read_weather_inputs.py

def read_weather_inputs(
    clock_sctruct: "ClockStruct",
    weather_df: "DataFrame") -> "DataFrame":
    """
    Clip weather to start and end simulation dates

    Arguments:

        clock_sctruct (ClockStruct): ClockStruct object

        weather_df (DataFrame): weather dataframe

    Returns:

        weather_df (DataFrame): clipped weather dataframe

    """

    # get the start and end dates of simulation
    start_date = clock_sctruct.simulation_start_date
    end_date = clock_sctruct.simulation_end_date

    if weather_df.Date.iloc[0] > start_date:
        raise ValueError(
            "The first date of the climate data cannot be longer than the start date of the model."
        )

    if weather_df.Date.iloc[-1] < end_date:
        raise ValueError(
            "The model end date cannot be longer than the last date of climate data."
        )

    # remove weather data outside of simulation dates
    weather_df = weather_df[weather_df.Date >= start_date]
    weather_df = weather_df[weather_df.Date <= end_date]

    return weather_df