Skip to content

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 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
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 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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_simulation_end_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 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
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":
    """
    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

            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
            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(((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 + 2.29) * np.log((((CCi/crop.CCx) - 1) / -0.05) + 1)) / (3.33 * tGDD) 
            # 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 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
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.simulation_end_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 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
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 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
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 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
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 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
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 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
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 
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
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.simulation_end_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 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
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