[WIP] Expanded planted area

scenarios/README.md

@@ -90,7 +90,6 @@ simulations:
   - **scenario**:
     - `no_resilient_foods` - Uses current food production only
     - `all_resilient_foods` - Enables scaled up seaweed, cellulosic sugar, methane SCP
-    - `all_resilient_foods_and_more_area` - Enables resilient foods and expands cropland area
     - `seaweed` - Enables scaled up seaweed production
     - `methane_scp` - Enables scaled up methane SCP production
     - `cellulosic_sugar` - Enables scaled up cellulosic sugar production
@@ -151,6 +150,10 @@ simulations:
   - **meat_strategy**:
     - `reduce_breeding` - Reduce breeding in accordance to "reduced" column in data/no_food_trade/animal_feed_data/species_options.csv. Essentiall, nearly all breeding is stopped immediately after the nuclear winter.
     - `baseline_breeding` - Normal breeding levels to maintain livestock levels, "baseline" column in data/no_food_trade/animal_feed_data/species_options.csv
+  - **expanded_area**:
+    - `none` - no expanded planted area shall be considered
+    - `no_trade` - expanded planted area with no equipment trade: each country relies solely on the equipment within its borders to clear its respective productive area, assuming that machines are fully mobile within the country they belong to and that productive area within a country is one combined region. It is assumed a country’s machinery is focused on land clearing, but since the equipment is assumed to only be operating 56 hours/week, other tasks could be performed with multiple shifts.
+    - `export_pool` - expanded planted area with export pool: a machinery inventory constituted by machines loaned by several countries without productive area that other countries with less machinery can tap into to increase their fleet, making the land clearing process faster. We stipulated countries without productive area loan 30% and 40% of their construction and forestry machines to the export pool, respectively. It is also assumed the machinery distribution is maintained constant after this loan, according to the fraction occupied by each equipment type in the global machinery stock.
 
 ## Detailed Notes
 

scenarios/USA_expanded.yaml

@@ -0,0 +1,25 @@
+settings:
+  countries: USA
+  NMONTHS: 120
+simulations:  
+  # Baseline Model By Country (converted from the provided code)
+  baseline_model_by_country:
+    title: Baseline Model by Country
+    scale: country
+    seasonality: country
+    grasses: baseline
+    crop_disruption: country_nuclear_winter
+    scenario: no_resilient_foods
+    fish: baseline
+    waste: baseline_in_country
+    nutrition: baseline
+    intake_constraints: enabled
+    buffer: baseline
+    shutoff: continued
+    cull: do_eat_culled
+    fat: not_required
+    protein: not_required
+    meat_strategy: baseline_breeding
+    stored_food: baseline
+    end_simulation_stocks_ratio: baseline
+    expanded_area: no_trade 

scenarios/argentina.yaml

@@ -109,9 +109,9 @@ simulations:
 
   # Nuclear winter with resilient foods & more area
   argentina_net_nuclear_resilient_more_area:
-    title: Argentina Net Food Production, Resilient Foods
+    title: Argentina Net Food Production, Resilient Foods, Expanded Area
     scale: country
-    scenario: all_resilient_foods_and_more_area
+    scenario: all_resilient_foods
     seasonality: country
     grasses: country_nuclear_winter
     crop_disruption: country_nuclear_winter 
@@ -126,3 +126,4 @@ simulations:
     shutoff: long_delayed_shutoff
     cull: do_eat_culled
     meat_strategy: reduce_breeding
+    expanded_area: export_pool 

scenarios/eu_countries.yaml

@@ -145,9 +145,9 @@ simulations:
 
   # Nuclear winter with resilient foods & more area
   net_nuclear_resilient_more_area:
-    title: Net Food Production, Resilient Foods
+    title: Net Food Production, Resilient Foods, Expanded Area
     scale: country
-    scenario: all_resilient_foods_and_more_area
+    scenario: all_resilient_foods
     seasonality: country
     grasses: country_nuclear_winter
     crop_disruption: country_nuclear_winter 
@@ -162,3 +162,4 @@ simulations:
     shutoff: long_delayed_shutoff
     cull: do_eat_culled
     meat_strategy: reduce_breeding
+    expanded_area: export_pool 

src/food_system/outdoor_crops.py

@@ -7,8 +7,9 @@
 """
 
 import numpy as np
-from src.utilities.plotter import Plotter
+
 from src.food_system.food import Food
+from src.utilities.plotter import Plotter
 
 
 class OutdoorCrops:
@@ -324,8 +325,8 @@ def calculate_monthly_production(self, constants_for_params):
 
         self.assign_reduction_from_climate_impact(constants_for_params)
 
-        if constants_for_params["RATIO_INCREASED_CROP_AREA"] > 1:
-            self.assign_increase_from_increased_cultivated_area(constants_for_params)
+        if constants_for_params["EXPANDED_AREA"] != "none":
+            self.assign_increase_from_expanded_planted_area(constants_for_params)
 
         PLOT_WITH_SEASONALITY_FLAG = False
         if PLOT_WITH_SEASONALITY_FLAG:
@@ -336,31 +337,43 @@ def calculate_monthly_production(self, constants_for_params):
             )
             Plotter.plot_monthly_reductions_seasonally(ratios)
 
-    def assign_increase_from_increased_cultivated_area(self, constants_for_params):
-        # Constants
-
-        N = constants_for_params["INITIAL_HARVEST_DURATION_IN_MONTHS"]
-        max_value = constants_for_params["RATIO_INCREASED_CROP_AREA"]
-        total_months = (
-            constants_for_params["NUMBER_YEARS_TAKES_TO_REACH_INCREASED_AREA"] * 12
-        )
-
-        # Create the linspace
-        linspace = np.ones(self.NMONTHS)
-
-        # Calculate the increment per month after the delay
-        increment = (max_value - 1) / (total_months - N)
-
-        # Update the linspace values
-        for i in range(N, total_months):
-            linspace[i] = 1 + (i - N) * increment
+    def assign_increase_from_expanded_planted_area(
+        self, constants_for_params: dict
+    ) -> None:
+        """
+        Add expanded planted area crop yield to outdoor growing crop Calories count.
 
-        # Maintain the value after the specified number of years
-        linspace[total_months:] = max_value
+        Arguments:
+            constants_for_params (dict): a dictionary containing the constants for parameters.
 
-        # increase the kcals grown using the increased ratio due to increased area
-        for i in range(self.NMONTHS):
-            self.KCALS_GROWN[i] = self.KCALS_GROWN[i] * linspace[i]
+        Returns:
+            None.
+        """
+        assert constants_for_params["EXPANDED_AREA"] != "none"
+        expanded_area_yearly_yield = [
+            constants_for_params[
+                f"expanded_area_{constants_for_params["EXPANDED_AREA"]}_kcals_year{year}"
+            ]
+            for year in range(1, 11)
+        ]
+        expanded_area_monthly_yield = np.outer(
+            constants_for_params["SEASONALITY"], expanded_area_yearly_yield
+        ).T.ravel()[: self.NMONTHS]
+
+        # the expanded area model assumes that:
+        # ```land cleared in a month is immediately cultivated,
+        # and wheat is harvested from that area 8 months after cultivation,
+        # which we considered a conservative average
+        # since harvest could happen 4 to 12 months after clearing.```
+        # therefore we zero out first 8 months
+        expanded_area_monthly_yield[:9] = 0
+
+        self.KCALS_GROWN = list(
+            np.array(self.KCALS_GROWN) + expanded_area_monthly_yield
+        )
+        self.NO_RELOCATION_KCALS_GROWN = list(
+            np.array(self.NO_RELOCATION_KCALS_GROWN) + expanded_area_monthly_yield
+        )
 
     def assign_reduction_from_climate_impact(self, constants_for_params):
         self.KCALS_GROWN = []

src/import_scripts_no_food_trade/create_expanded_area_csv.py

@@ -0,0 +1,202 @@
+from pathlib import Path
+
+import git
+import numpy as np
+import pandas as pd
+from git.types import PathLike
+
+
+def get_wheat_nutrients(repo_root: PathLike) -> pd.Series:
+    """
+    Retrieve nutritional information for wheat from
+    the "Supplemental_Data.xlsx" file.
+    NOTE: other crops are currently not available in the
+    expanded area model therefore we specificlly fetch
+    wheat only.
+    We can refactor this when once we have more data.
+
+    Arguments:
+        repo_root (git.PathLike): the working directory of the repository,
+            typically fetched via:
+            `git.Repo(".", search_parent_directories=True).working_dir`
+
+    Returns:
+        pandas.Series: Calories, Fat and Protein per kg of wheat.
+    """
+    nutrition = pd.read_excel(
+        Path(repo_root) / "data" / "Supplemental_Data.xlsx",
+        sheet_name="Nutrition",
+        usecols="A:D",
+    )
+    nutrition = nutrition.set_index(nutrition.columns[0])
+    nutrition.columns = ["Calories", "Protein", "Fat"]
+    wheat_nutrients = nutrition.loc["Wheat", :]
+    return wheat_nutrients
+
+
+def get_expanded_area_monthly(
+    repo_root: PathLike, equipment_trade="no trade"
+) -> pd.DataFrame:
+    """
+    Reads the expanded area model output CSV and
+    converts appropriate data to a dataframe.
+
+
+    Arguments:
+        repo_root (git.PathLike): the working directory of the repository,
+            typically fetched via:
+            `git.Repo(".", search_parent_directories=True).working_dir`
+        equipment_trade (str, optional): specifies expanded are model scenario.
+            Can either be `no_trade` or `export_pool`. Defaults to `no_trade`.
+
+    Returns:
+        pandas.DataFrame: additional wheat yields per month per country
+            from expanded area during a nuclear winter.
+    """
+    start_col = {
+        "no_trade": 17,
+        "export_pool": 128,
+    }
+    if equipment_trade not in start_col:
+        print(f"Unknown equipment trade scenario: {equipment_trade}.")
+        print("Using `no_trade`.")
+        equipment_trade = "no_trade"
+    expanded_area_yields = (
+        pd.read_csv(
+            Path(repo_root)
+            / "data"
+            / "no_food_trade"
+            / "raw_data"
+            / "luisa_expanded_planted_area.csv",
+            index_col=False,
+            dtype=str,
+        )
+        .iloc[16:-1, start_col[equipment_trade] : start_col[equipment_trade] + 87]
+        .reset_index(drop=True)
+    )
+    expanded_area_yields.columns = expanded_area_yields.iloc[0]
+    expanded_area_yields = expanded_area_yields.drop(expanded_area_yields.index[0])
+    # NOTE: using "iloc" and other shenanigans here
+    # because when using read_csv parameters `skiprows`
+    # and `usecols`, the results were sometimes incorrect;
+    # pandas would try to force str to be float for some reason
+    # similar to these cases:
+    # https://stackoverflow.com/questions/55503007/how-to-remove-the-0-when-reading-csv-with-pandas
+    # https://stackoverflow.com/questions/54132508/pandas-adding-decimal-points-when-using-read-csv
+    # however, my attempts to fix this failed.
+    expanded_area_yields = expanded_area_yields.drop(
+        columns=["Annual yield (kg per ha)"]
+    ).set_index("Country (ISO 3DIGIT)")
+    expanded_area_yields.index.name = "iso3"
+    # we want 0=May; the start of the nuclear winter
+    expanded_area_yields.columns = np.array(expanded_area_yields.columns, dtype=int) - 1
+    expanded_area_yields = expanded_area_yields.drop(columns=-1)
+    # we don't need countries with all NaNs
+    expanded_area_yields = expanded_area_yields.apply(pd.to_numeric, errors="coerce")
+    expanded_area_yields = expanded_area_yields.replace(0, np.nan)
+    expanded_area_yields = expanded_area_yields.dropna(axis=0, how="all")
+    expanded_area_yields = expanded_area_yields.replace(np.nan, 0)
+    return expanded_area_yields
+
+
+def calculate_yearly_yields(expanded_area_yields_monthly: pd.DataFrame) -> pd.DataFrame:
+    """
+    Convert monthly values to yearly yields.
+    As we want to apply seasonality later on in the code,
+    all we need is aggregate yearly values of the additional yields.
+
+    Arguments:
+        expanded_area_yields_monthly (pandas.DataFrame): Data frame containing
+            monthly yields for each country.
+
+    Returns:
+        pandas.DataFrame: a data frame with yearly yields for each country.
+    """
+    expanded_area_yields_yearly = expanded_area_yields_monthly.copy()
+    expanded_area_yields_yearly.columns = pd.to_datetime(
+        expanded_area_yields_yearly.columns, unit="M"
+    )
+    expanded_area_yields_yearly = expanded_area_yields_yearly.T.resample("A").sum().T
+    expanded_area_yields_yearly.columns = range(
+        len(expanded_area_yields_yearly.columns)
+    )
+    return expanded_area_yields_yearly
+
+
+def add_missing_years(expanded_area_yields_yearly: pd.DataFrame) -> pd.DataFrame:
+    """
+    Complete the data to conform with other data sets.
+    As of wrighting this, we have 7 years of expanded area predictions,
+    and assuming it plateaus after, we add columns with final 3 years
+    to match the 10-year time spand of the other datasets.
+    This can (and probably should) be modified to reflect actual
+    yield predictions in the final 3 years.
+
+    Arguments:
+        expanded_area_yields_monthly (pandas.DataFrame): Data frame containing
+            monthly yields for each country.
+
+    Returns:
+        pandas.DataFrame: a data frame with yearly yields for each country
+            with additional years appended.
+    """
+    expanded_area_yields_yearly = pd.concat(
+        [expanded_area_yields_yearly] + [expanded_area_yields_yearly.iloc[:, -1]] * 3,
+        axis=1,
+    )
+    expanded_area_yields_yearly.columns = range(
+        len(expanded_area_yields_yearly.columns)
+    )
+    return expanded_area_yields_yearly
+
+
+def main():
+    """
+    Process expanded planted area model output data and create
+    a CSV that can be later merged into the integrated model input files.
+    """
+    repo_root = git.Repo(".", search_parent_directories=True).working_dir
+    print("Processing expanded planted area data.")
+    wn = get_wheat_nutrients(repo_root)
+    expanded_area_yields = {}
+    for scenario in ["no_trade", "export_pool"]:
+        eam = get_expanded_area_monthly(repo_root, scenario)
+        eay = calculate_yearly_yields(eam)
+        eay = add_missing_years(eay)
+        # kg to kcal to dry caloric tonne
+        eay_kcal = eay.mul(wn["Calories"]).div(4000 * 1000)
+        eay_kcal.columns = [
+            f"expanded_area_{scenario}_kcals_year" + str(c + 1) for c in eay.columns
+        ]
+        eay_fat = eay.mul(wn["Fat"]).div(1e3)  # in tonnes
+        eay_fat.columns = [
+            f"expanded_area_{scenario}_fat_year" + str(c + 1) for c in eay.columns
+        ]
+        eay_protein = eay.mul(wn["Protein"]).div(1e3)  # in tonnes
+        eay_protein.columns = [
+            f"expanded_area_{scenario}_protein_year" + str(c + 1) for c in eay.columns
+        ]
+        eay = pd.concat([eay_kcal, eay_fat, eay_protein], axis=1)
+        eay = eay.sort_index()
+        expanded_area_yields[scenario] = eay
+    expanded_area_yields = pd.concat(expanded_area_yields.values(), axis=1)
+    expanded_area_yields.to_csv(
+        Path(repo_root)
+        / "data"
+        / "no_food_trade"
+        / "processed_data"
+        / "expanded_area.csv"
+    )
+    print("Saved expanded area macros to:")
+    print(
+        Path(repo_root)
+        / "data"
+        / "no_food_trade"
+        / "processed_data"
+        / "expanded_area.csv"
+    )
+    assert ~expanded_area_yields.isnull().to_numpy().any()
+
+
+if __name__ == "__main__":
+    main()