From 22095f91ea830dcfdcf14da08110e377d839185a Mon Sep 17 00:00:00 2001
From: Bart Arendarczyk <s1924442@ed.ac.uk>
Date: Sat, 4 Mar 2023 15:26:57 +0000
Subject: [PATCH] Forestry and carbon fixes & some parameter adjustments.
---
src/ac/ed/lurg/ModelConfig.java | 34 ++--
src/ac/ed/lurg/ModelMain.java | 2 +-
src/ac/ed/lurg/carbon/CarbonFluxItem.java | 17 +-
src/ac/ed/lurg/country/CountryAgent.java | 3 +-
.../country/gams/GamsRasterOptimiser.java | 44 ++++-
src/ac/ed/lurg/demand/WoodDemandManager.java | 9 +-
src/ac/ed/lurg/forestry/ForestManager.java | 180 +++++++++---------
src/ac/ed/lurg/forestry/WoodHarvestItem.java | 34 ++++
src/ac/ed/lurg/landuse/LandCoverTile.java | 3 -
src/ac/ed/lurg/landuse/LandUseItem.java | 53 +++++-
src/ac/ed/lurg/landuse/WoodUsageReader.java | 4 +-
src/ac/ed/lurg/types/CropType.java | 20 +-
src/ac/ed/lurg/types/WoodType.java | 2 +-
13 files changed, 258 insertions(+), 147 deletions(-)
create mode 100644 src/ac/ed/lurg/forestry/WoodHarvestItem.java
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 6fe441a2..b53c18e5 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -200,14 +200,14 @@ public class ModelConfig {
public static final String YIELD_FILENAME = getProperty("YIELD_FILENAME", "yield.out");
public static final boolean PASTURE_FERT_RESPONSE_FROM_LPJ = getBooleanProperty("PASTURE_FERT_RESPONSE_FROM_LPJ", false);;
- public static final double CALIB_FACTOR_CEREAL_C3 = getDoubleProperty("CALIB_FACTOR_CEREAL_C3", 0.904);
- public static final double CALIB_FACTOR_CEREAL_C4 = getDoubleProperty("CALIB_FACTOR_CEREAL_C4", 0.504);
- public static final double CALIB_FACTOR_MISCANTHUS = getDoubleProperty("CALIB_FACTOR_MISCANTHUS", 2.047);
- public static final double CALIB_FACTOR_RICE = getDoubleProperty("CALIB_FACTOR_RICE", 0.978);
- public static final double CALIB_FACTOR_OILCROPS = getDoubleProperty("CALIB_FACTOR_OILCROPS", 1.0);
- public static final double CALIB_FACTOR_PULSES = getDoubleProperty("CALIB_FACTOR_PULSES", 0.347);
- public static final double CALIB_FACTOR_STARCHY_ROOTS = getDoubleProperty("CALIB_FACTOR_STARCHY_ROOTS",4.679);
- public static final double CALIB_FACTOR_FRUITVEG = getDoubleProperty("CALIB_FACTOR_FRUITVEG",3.377);
+ public static final double CALIB_FACTOR_CEREAL_C3 = getDoubleProperty("CALIB_FACTOR_CEREAL_C3", 1.07);
+ public static final double CALIB_FACTOR_CEREAL_C4 = getDoubleProperty("CALIB_FACTOR_CEREAL_C4", 0.707);
+ public static final double CALIB_FACTOR_MISCANTHUS = getDoubleProperty("CALIB_FACTOR_MISCANTHUS", 2.481);
+ public static final double CALIB_FACTOR_RICE = getDoubleProperty("CALIB_FACTOR_RICE", 1.464);
+ public static final double CALIB_FACTOR_OILCROPS = getDoubleProperty("CALIB_FACTOR_OILCROPS", 1.47);
+ public static final double CALIB_FACTOR_PULSES = getDoubleProperty("CALIB_FACTOR_PULSES", 0.446);
+ public static final double CALIB_FACTOR_STARCHY_ROOTS = getDoubleProperty("CALIB_FACTOR_STARCHY_ROOTS",5.656);
+ public static final double CALIB_FACTOR_FRUITVEG = getDoubleProperty("CALIB_FACTOR_FRUITVEG",4.812);
public static final double CALIB_FACTOR_SUGAR = getDoubleProperty("CALIB_FACTOR_SUGAR", 1.0);
public static final double CALIB_FACTOR_PASTURE = getDoubleProperty("CALIB_FACTOR_PASTURE", 1.0);
@@ -363,7 +363,7 @@ public class ModelConfig {
public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2019);
// Import export limits
- public static final double ANNUAL_MAX_IMPORT_CHANGE = getDoubleProperty("ANNUAL_MAX_IMPORT_CHANGE", 0.02);
+ public static final double ANNUAL_MAX_IMPORT_CHANGE = getDoubleProperty("ANNUAL_MAX_IMPORT_CHANGE", 0.05);
public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", ANNUAL_MAX_IMPORT_CHANGE*TIMESTEP_SIZE);
public static final double IMPORT_CHANGE_ELASTICITY = getDoubleProperty("IMPORT_CHANGE_ELASTICITY", 3.0); // higher value = steeper trade imbalance adjustment
@@ -394,7 +394,7 @@ public class ModelConfig {
public static final double PASTURE_HARVEST_FRACTION = getDoubleProperty("PASTURE_HARVEST_FRACTION", 0.5);
public static final double MEAT_EFFICIENCY = getDoubleProperty("MEAT_EFFICIENCY", 1.0); // 'meat' is includes feed conversion ratio already, this is tech. change or similar
- public static final double IRRIGIATION_EFFICIENCY = getDoubleProperty("IRRIGIATION_EFFICIENCY", 0.6);
+ public static final double IRRIGIATION_EFFICIENCY = getDoubleProperty("IRRIGIATION_EFFICIENCY", 0.33);
public static final int ELLIOTT_BASEYEAR = 2010;
public static final double ENVIRONMENTAL_WATER_CONSTRAINT = getDoubleProperty("ENVIRONMENTAL_WATER_CONSTRAINT", 0.5); // change with care, as due to normalisation it might not have the impact you first imagine
public static final double OTHER_WATER_USE_FACTOR = getDoubleProperty("OTHER_WATER_USE_FACTOR", 1.0);
@@ -425,8 +425,8 @@ public class ModelConfig {
public static final boolean RESET_ENERGYCROP_PRICE = getBooleanProperty("RESET_ENERGYCROP_PRICE", true); // Resets price after calibration to avoid problems due to low initial demand
// public static final double BIOENERGY_HEATING_VALUE_GJ_PER_T = getDoubleProperty("BIOENERGY_HEATING_VALUE_GJ_PER_T", 17.5); // GJ per t DM
- public static final double MARKET_LAMBDA = getDoubleProperty("MARKET_LAMBDA", 0.2); // controls international market price adjustment rate
- public static final double MARKET_DAMPING = getDoubleProperty("MARKET_DAMPING", 2.0); // controls international market price adjustment rate
+ public static final double MARKET_LAMBDA = getDoubleProperty("MARKET_LAMBDA", 0.04); // controls international market price adjustment rate
+ public static final double MARKET_DAMPING = getDoubleProperty("MARKET_DAMPING", 0.4); // controls international market price adjustment rate
public static final String PRICE_UPDATE_METHOD = getProperty("PRICE_UPDATE_METHOD", "StockUseRatio"); // Options: StockUseRatio, MarketImbalance, StockChange
public static final double MAX_PRICE_INCREASE = getDoubleProperty("MAX_PRICE_INCREASE", 1.5);
@@ -445,10 +445,10 @@ public class ModelConfig {
public static final double UNHANDLED_CROP_RATE = getDoubleProperty("UNHANDLED_CROP_RATE", 0.0498); // mostly forage crops
public static final double SETASIDE_RATE = getDoubleProperty("SETASIDE_RATE", 0.103); // includes aside, fallow and failed cropland areas
- public static final double OTHER_INTENSITY_COST = getDoubleProperty("OTHER_INTENSITY_COST", 0.6);
+ public static final double OTHER_INTENSITY_COST = getDoubleProperty("OTHER_INTENSITY_COST", 0.8);
public static final double OTHER_INTENSITY_PARAM = getDoubleProperty("OTHER_INTENSITY_PARAM", 3.22);
- public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 0.0002);
+ public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 0.00035);
public static final double IRRIG_COST_MULTIPLIER = getDoubleProperty("IRRIG_COST_MULTIPLIER", 1.0);
public static final double FERTILISER_COST_PER_T = getDoubleProperty("FERTILISER_COST_PER_T", 1.6); // $500/t, 18% N/t
public static final double FERTILISER_MAX_COST = FERTILISER_COST_PER_T * MAX_FERT_AMOUNT/1000;
@@ -522,11 +522,11 @@ public class ModelConfig {
public static final double FUELWOOD_DEMAND_ELASTICITY = getDoubleProperty("FUELDWOOD_DEMAND_ELASTICITY", -0.3598551);
public static final double WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 0.3); // m3 to tC-eq p.16 [https://doi.org/10.5194/gmd-13-5425-2020]
public static final double VEGETATION_CLEARING_COST = getDoubleProperty("VEGETATION_CLEARING_COST", 0.001); //$1000/tC
- public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_ESTABLISHMENT_COST", 3.0) : 0.0; // establishment, management etc. $1000/ha
+ public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_ESTABLISHMENT_COST", 5.0) : 0.0; // establishment, management etc. $1000/ha
public static final double FOREST_BASE_COST = getDoubleProperty("FOREST_BASE_COST", 0.05);
public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.2); //$1000/tC
- public static final double INIT_ROUNDWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_ROUNDWOOD_PRICE", 0.36) : 0.0; // $1000/tC-eq
- public static final double INIT_FUELWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_FUELWOOD_PRICE", 0.24) : 0.0; // $1000/tC-eq
+ public static final double INIT_ROUNDWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_ROUNDWOOD_PRICE", 0.1) : 0.0; // $1000/tC-eq
+ public static final double INIT_FUELWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_FUELWOOD_PRICE", 0.05) : 0.0; // $1000/tC-eq
public static final int MAX_FOREST_ROTATION_PERIOD = getIntProperty("CARBON_WOOD_MAX_TIME", 160); // upper data limit, years
public static final int LAND_AGE_GROUP_SIZE = getIntProperty("LAND_AGE_GROUP_SIZE", 5); // grouping size for land cover age
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 4b725b71..b17122b7 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -370,7 +370,7 @@ public class ModelMain {
continue;
double prod = woodUsage.getHarvest();
- double netImports = woodUsage.getNetImport() - woodUsage.getLucHarvest(); // does not include transport losses
+ double netImports = woodUsage.getNetImport();
CountryPrice px = country.getCurrentCountryWoodPrices().get(woodType);
double importPrice = px.getImportPrice();
diff --git a/src/ac/ed/lurg/carbon/CarbonFluxItem.java b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
index d020dde9..cc164a0e 100644
--- a/src/ac/ed/lurg/carbon/CarbonFluxItem.java
+++ b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
@@ -53,18 +53,17 @@ public class CarbonFluxItem implements RasterItem, Serializable {
}
double meanFlux = totalFlux / totalArea;
carbonFluxRate.put(key, meanFlux);
-
- // Carbon credit rate for unchanged land cover
-
-
}
-
+
+ // Potential carbon flux from this existing land cover
double fromLcCumulFlux = 0;
for (int age : landTile.getAgeKeys()) {
for (int i = age; i <= age + ModelConfig.CARBON_HORIZON; i++) {
fromLcCumulFlux += cFluxes[i];
}
}
+
+ // Potential carbon flux if converting to this land cover
double toLcCumulFlux = 0;
for (int age = 0; age <= ModelConfig.CARBON_HORIZON; age++) {
toLcCumulFlux += cFluxes[age];
@@ -77,8 +76,12 @@ public class CarbonFluxItem implements RasterItem, Serializable {
for (LandCoverType lc : LandCoverType.getConvertibleTypes()) {
LccKey fromKey = new LccKey(lcType, lc);
LccKey toKey = new LccKey(lc, lcType);
- carbonCreditRate.merge(fromKey, fromLcCumulFlux, (a, b) -> a - b);
- carbonCreditRate.merge(toKey, toLcCumulFlux, (a, b) -> a + b);
+ if (lcType.equals(lc)) { // no net difference if land cover stays same
+ carbonCreditRate.put(fromKey, 0.0);
+ } else {
+ carbonCreditRate.merge(fromKey, -fromLcCumulFlux, Double::sum);
+ carbonCreditRate.merge(toKey, toLcCumulFlux, Double::sum);
+ }
}
}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index ee700c56..ded725a0 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -261,8 +261,7 @@ public class CountryAgent extends AbstractCountryAgent {
if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
changeUp = changeDown = 0.0;
} else {
- changeUp = maxOfProdOrSupply * globalWoodPrices.get(woodType).getMaxImportChange();
- changeDown = maxOfProdOrSupply * globalWoodPrices.get(woodType).getMaxExportChange();
+ changeUp = changeDown = maxOfProdOrSupply * ModelConfig.MAX_IMPORT_CHANGE;
}
woodTradeConstraints.put(woodType, new TradeConstraint(baseTrade - changeDown, baseTrade + changeUp));
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 235ab9c0..fa726076 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -7,6 +7,7 @@ import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.carbon.CarbonFluxItem;
import ac.ed.lurg.country.LandCoverChangeItem;
import ac.ed.lurg.forestry.ForestManager;
+import ac.ed.lurg.forestry.WoodHarvestItem;
import ac.ed.lurg.forestry.WoodYieldData;
import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.*;
@@ -97,7 +98,7 @@ public class GamsRasterOptimiser {
private RasterSet<LandUseItem> allocAreas(Map<Integer, ? extends LandUseItem> prevAreasAgg, GamsLocationOutput gamsOutput, int year) {
RasterSet<LandUseItem> newLandUseRaster = createWithSameLandCovers(rasterInputData.getPreviousLandUses());
- // Order is important here
+ // *** Order is important here ***
// Disaggregate cluster land cover changes to raster grid
Map<RasterKey, List<LandCoverChangeItem>> lcChanges = disaggregateLandCoverChanges(gamsOutput.getLandCoverChanges(),
@@ -110,11 +111,46 @@ public class GamsRasterOptimiser {
// Allocate land cover changes
allocLandCoverChanges(newLandUseRaster, lcChanges);
- // Calculate land cover changes for timber forest rotation and clear-cut
- Map<RasterKey, List<LandCoverChangeItem>> lcChangesRota = ForestManager.handleRotaWood(
- newLandUseRaster, rasterInputData.getCountryInput().getWoodDemand(),
+ Map<LandCoverType, Double> foo = new HashMap<>();
+ for (LandUseItem i : newLandUseRaster.values()) {
+ for (LandCoverType c : LandCoverType.getConvertibleTypes()) {
+ foo.merge(c, i.getLandCoverArea(c, LandProtectionType.CONVERTIBLE), Double::sum);
+ }
+ }
+
+ // Set wood type forest fractions
+ Map<Integer, Map<WoodType, Double>> woodAreas = gamsOutput.getForestAreas();
+ for (int locId : woodAreas.keySet()) {
+ Set<RasterKey> keys = new HashSet<RasterKey>();
+ for (Entry<RasterKey, IntegerRasterItem> mapEntry : mapping.entrySet()) {
+ IntegerRasterItem iri = mapEntry.getValue();
+ if (locId == iri.getInt())
+ keys.add(mapEntry.getKey());
+ }
+ Map<WoodType, Double> areaMap = woodAreas.get(locId);
+ double totalArea = areaMap.values().stream().reduce(0.0, Double::sum);
+ for (WoodType wType : WoodType.values()) {
+ double fraction = totalArea > 0 ? areaMap.get(wType) / totalArea : 0;
+ for (RasterKey key : keys) {
+ LandUseItem luItem = newLandUseRaster.get(key);
+ luItem.setForestFraction(wType, fraction);
+ }
+ }
+ }
+
+ // Calculate land cover changes for timber forest rotation and clear-cut. Resets harvested stands to age 0
+ Map<RasterKey, WoodHarvestItem> rotaHarvestMap = ForestManager.handleRotaWood(newLandUseRaster,
rasterInputData.getWoodYields(), gamsOutput, mapping);
+ // Set forest rotations
+ for (RasterKey key : rotaHarvestMap.keySet()) {
+ WoodHarvestItem whItem = rotaHarvestMap.get(key);
+ LandUseItem luItem = newLandUseRaster.get(key);
+ luItem.setWoodHarvests(whItem);
+ }
+
+
+
// Add wood produced from LCC to output
Map<WoodType, WoodUsageData> woodUsageMap = gamsOutput.getWoodUsageData();
for (WoodType wType : lccWoodAmount.keySet()) {
diff --git a/src/ac/ed/lurg/demand/WoodDemandManager.java b/src/ac/ed/lurg/demand/WoodDemandManager.java
index d9fcac3f..2ef60933 100644
--- a/src/ac/ed/lurg/demand/WoodDemandManager.java
+++ b/src/ac/ed/lurg/demand/WoodDemandManager.java
@@ -1,5 +1,7 @@
package ac.ed.lurg.demand;
+import java.io.Serial;
+import java.io.Serializable;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
@@ -11,7 +13,9 @@ import ac.ed.lurg.types.WoodType;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.utils.StringTabularReader;
-public class WoodDemandManager {
+public class WoodDemandManager implements Serializable {
+ @Serial
+ private static final long serialVersionUID = 171124405614624397L;
private Map<SingleCountry, Map<WoodType, Double>> baseWoodDemand = new HashMap<SingleCountry, Map<WoodType, Double>>();
public WoodDemandManager() {
@@ -26,8 +30,7 @@ public class WoodDemandManager {
for (Map<String, String> row : rowList) {
String countryCode = row.get("Iso3");
String woodTypeStr = row.get("Type");
- double demand = Double.parseDouble(row.get("Demand"));
- demand *= ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; // converting from m3 to MtC-eq
+ double demand = Double.parseDouble(row.get("Demand")) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR;
SingleCountry sc = CountryManager.getInstance().getForCode(countryCode);
WoodType woodType = WoodType.getForName(woodTypeStr);
Map<WoodType, Double> wMap = baseWoodDemand.computeIfAbsent(sc, k -> new HashMap<WoodType, Double>());
diff --git a/src/ac/ed/lurg/forestry/ForestManager.java b/src/ac/ed/lurg/forestry/ForestManager.java
index 44a4bbc5..4fbb77c4 100644
--- a/src/ac/ed/lurg/forestry/ForestManager.java
+++ b/src/ac/ed/lurg/forestry/ForestManager.java
@@ -23,11 +23,11 @@ public class ForestManager {
RasterSet<LandUseItem> inputLandUseRaster) {
Map<WoodType, Double> woodAmount = new HashMap<>();
+ for (WoodType w : WoodType.values()) {
+ woodAmount.put(w, 0.0);
+ }
if (!ModelConfig.IS_FORESTRY_ON) {
- for (WoodType w : WoodType.values()) {
- woodAmount.put(w, 0.0);
- }
return woodAmount;
}
@@ -41,9 +41,10 @@ public class ForestManager {
LandCoverType toLc = lccItem.getToLandCover();
double changedArea = lccItem.getArea();
double totalArea = tiles.getTotalArea(fromLc, LandProtectionType.CONVERTIBLE);
- if (changedArea > totalArea) {
+ if (changedArea > totalArea && Math.abs(changedArea - totalArea) > 1e-10) {
LogWriter.printlnError("Error: ForestManager changed area exceeds total area.");
}
+ changedArea = Math.min(changedArea, totalArea); // deal with floating point errors
double fraction = changedArea / totalArea;
// Assume stands are converted proportionally by area.
// Type of wood harvested depends on stand age.
@@ -65,114 +66,113 @@ public class ForestManager {
}
// Calculates area of timber forest cleared to meet demand
- public static Map<RasterKey, List<LandCoverChangeItem>> handleRotaWood(RasterSet<LandUseItem> outputLandUseRaster,
- Map<WoodType, Double> woodDemand,
+ public static Map<RasterKey, WoodHarvestItem> handleRotaWood(RasterSet<LandUseItem> outputLandUseRaster,
RasterSet<WoodYieldItem> woodYieldRaster,
- GamsLocationOutput locationOutput,
- RasterSet<IntegerRasterItem> mapping) {
+ GamsLocationOutput locationOutput,
+ RasterSet<IntegerRasterItem> mapping) {
// Since some trees could be used for either roundwood or fuelwood, we harvest the more valuable roundwood first
final List<WoodType> harvestOrder = new ArrayList<>(Arrays.asList(WoodType.IND_ROUNDWOOD, WoodType.FUELWOOD));
- Map<RasterKey, List<LandCoverChangeItem>> allChanges = new HashMap<>();
+ Map<RasterKey, WoodHarvestItem> harvestMap = new HashMap<>();
if (!ModelConfig.IS_FORESTRY_ON) {
- return allChanges;
+ return harvestMap;
}
for (WoodType wType : harvestOrder) {
- Map<RasterKey, List<LandCoverChangeItem>> changes = handleRotaWoodForType(outputLandUseRaster, woodDemand,
- woodYieldRaster, locationOutput, mapping, wType);
- for (RasterKey key : changes.keySet()) {
- List<LandCoverChangeItem> changesList = changes.get(key);
- List<LandCoverChangeItem> changesListAll = allChanges.computeIfAbsent(key, k -> new ArrayList<>());
- changesListAll.addAll(changesList);
- }
+ handleRotaWoodForType(outputLandUseRaster, woodYieldRaster, locationOutput, mapping, wType, harvestMap);
}
- return allChanges;
+ return harvestMap;
}
- private static Map<RasterKey, List<LandCoverChangeItem>> handleRotaWoodForType(RasterSet<LandUseItem> outputLandUseRaster,
- Map<WoodType, Double> woodDemand,
- RasterSet<WoodYieldItem> woodYieldRaster,
- GamsLocationOutput locationOutput,
- RasterSet<IntegerRasterItem> mapping,
- WoodType wType) {
- Map<Integer, Map<WoodType, Double>> rotationPeriods = locationOutput.getRotationPeriods();
- Map<RasterKey, List<LandCoverChangeItem>> lcChanges = new HashMap<>();
-
- double demand = woodDemand.get(wType);
- if (demand == 0)
- return lcChanges;
-
- // First calculate how much wood is available
- Map<RasterKey, Double> stockMap = new HashMap<>();
- for (RasterKey key : outputLandUseRaster.keySet()) {
- LandUseItem luItem = outputLandUseRaster.get(key);
- LandCoverTile tiles = luItem.getLandCoverTiles();
- WoodYieldItem wyItem = woodYieldRaster.get(key);
- Integer locId = mapping.get(key).getInt();
- double stock = 0;
- double rotationAge = rotationPeriods.get(locId).get(wType);
- for (int age : tiles.getAgeKeys()) {
- if (age >= rotationAge) {
- double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age);
- double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age);
- stock += yield * area;
- }
- }
- stockMap.put(key, stock);
- }
+ private static void handleRotaWoodForType(RasterSet<LandUseItem> outputLandUseRaster,
+ RasterSet<WoodYieldItem> woodYieldRaster,
+ GamsLocationOutput locationOutput,
+ RasterSet<IntegerRasterItem> mapping,
+ WoodType wType,
+ Map<RasterKey, WoodHarvestItem> harvestMap) {
+
+ double demand = locationOutput.getWoodUsageData().get(wType).getHarvest()
+ - locationOutput.getWoodUsageData().get(wType).getLucHarvest();
+ if (demand <= 0) {
+ return;
+ }
- double totalStock = stockMap.values().stream().reduce(0.0, Double::sum);
- if (demand > totalStock) {
- LogWriter.printlnWarning("Insufficient " + wType.getName() + " to meet demand. Shortfall: " + (demand - totalStock));
- }
- double harvestFraction = Math.min(demand / totalStock, 1);
- // Now harvest wood to meet demand.
- // Side effect! Resets land cover age to 0 for harvested stands
- for (RasterKey key : outputLandUseRaster.keySet()) {
- LandUseItem luItem = outputLandUseRaster.get(key);
- LandCoverTile tiles = luItem.getLandCoverTiles();
- WoodYieldItem wyItem = woodYieldRaster.get(key);
- Integer locId = mapping.get(key).getInt();
- double rotationAge = rotationPeriods.get(locId).get(wType);
- double totalAmountToRemove = stockMap.get(key) * harvestFraction;
- List<Integer> ageList = new ArrayList<>(tiles.getAgeKeys());
- List<Integer> rotationAgeList = new ArrayList<>();
- for (int i : ageList) {
- if (i >= rotationAge) {
- rotationAgeList.add(i);
- }
- }
- // Harvest in descending age order
- rotationAgeList.sort(Collections.reverseOrder());
- for (int age : rotationAgeList) {
- if (totalAmountToRemove == 0)
- break;
+ // First calculate how much wood is available
+ Map<RasterKey, Double> stockMap = new HashMap<>();
+ for (RasterKey key : outputLandUseRaster.keySet()) {
+ LandUseItem luItem = outputLandUseRaster.get(key);
+ LandCoverTile tiles = luItem.getLandCoverTiles();
+ WoodYieldItem wyItem = woodYieldRaster.get(key);
+ Integer locId = mapping.get(key).getInt();
+ double stock = 0;
+ double rotationAge = wType.getMinRotation();
+ for (int age : tiles.getAgeKeys()) {
+ if (age >= rotationAge) {
double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age);
- double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age);
- if (area == 0 || yield == 0)
+ if (area == 0)
continue;
- double standStock = area * yield;
- double amountToRemove = Math.min(totalAmountToRemove, standStock);
- double convertedArea = Math.min(area, amountToRemove / yield);
- LandCoverChangeItem lccItem = new LandCoverChangeItem(LandCoverType.TIMBER_FOREST,
- LandCoverType.TIMBER_FOREST, age, convertedArea);
- List<LandCoverChangeItem> mapForKey = lcChanges.computeIfAbsent(key, k -> new ArrayList<>());
- mapForKey.add(lccItem);
- if (!ModelConfig.IS_CALIBRATION_RUN)
- luItem.doForestRotation(age, convertedArea); // Side effect!
- totalAmountToRemove -= amountToRemove;
+ double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age);
+ stock += yield * area;
}
+ }
+ stockMap.put(key, stock);
+ }
- if (totalAmountToRemove > 0)
- LogWriter.printlnWarning("ForestManager: Could not meet wood demand.");
+ double totalStock = stockMap.values().stream().reduce(0.0, Double::sum);
+ LogWriter.println("Total stock for " + wType + " = " + totalStock);
+ if (demand > totalStock) {
+ LogWriter.printlnWarning("Insufficient " + wType.getName() + " to meet demand. Shortfall: " + (demand - totalStock));
+ }
+ if (totalStock < 1e-10) {
+ return;
+ }
+ double harvestFraction = Math.min(demand / totalStock, 1);
+ // Now harvest wood to meet demand.
+ for (RasterKey key : outputLandUseRaster.keySet()) {
+ LandUseItem luItem = outputLandUseRaster.get(key);
+ LandCoverTile tiles = luItem.getLandCoverTiles();
+ WoodYieldItem wyItem = woodYieldRaster.get(key);
+ Integer locId = mapping.get(key).getInt();
+ double rotationAge = wType.getMinRotation();
+ double totalAmountToRemove = stockMap.get(key) * harvestFraction;
+ if (totalAmountToRemove < 1e-10) {
+ continue;
+ }
+ // Assuming that oldest trees get harvested first
+ List<Integer> ageList = new ArrayList<>(tiles.getAgeKeys());
+ List<Integer> rotationAgeList = new ArrayList<>();
+ for (int i : ageList) {
+ if (i >= rotationAge) {
+ rotationAgeList.add(i);
+ }
+ }
+ // Harvest in descending age order
+ rotationAgeList.sort(Collections.reverseOrder());
+
+ for (int age : rotationAgeList) {
+ double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age);
+ double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age);
+ if (area == 0 || yield == 0)
+ continue;
+ double standStock = area * yield;
+ double amountToRemove = Math.min(totalAmountToRemove, standStock);
+ double convertedArea = Math.min(area, amountToRemove / yield);
+ LandCoverChangeItem lccItem = new LandCoverChangeItem(LandCoverType.TIMBER_FOREST,
+ LandCoverType.TIMBER_FOREST, age, convertedArea);
+ WoodHarvestItem whItem = harvestMap.computeIfAbsent(key, k -> new WoodHarvestItem());
+ whItem.addAreasHarvested(lccItem);
+ whItem.addWoodAmount(wType, amountToRemove);
+ luItem.doForestRotation(age, convertedArea); // WARNING! Side effect
+ totalAmountToRemove -= amountToRemove;
}
- return lcChanges;
+ if (totalAmountToRemove > 0)
+ LogWriter.printlnWarning("ForestManager: Could not meet wood demand.");
+ }
}
}
diff --git a/src/ac/ed/lurg/forestry/WoodHarvestItem.java b/src/ac/ed/lurg/forestry/WoodHarvestItem.java
new file mode 100644
index 00000000..228c118a
--- /dev/null
+++ b/src/ac/ed/lurg/forestry/WoodHarvestItem.java
@@ -0,0 +1,34 @@
+package ac.ed.lurg.forestry;
+
+import ac.ed.lurg.country.LandCoverChangeItem;
+import ac.ed.lurg.types.WoodType;
+
+import java.io.Serial;
+import java.io.Serializable;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+public class WoodHarvestItem implements Serializable {
+ @Serial
+ private static final long serialVersionUID = -79494635617251642L;
+ private List<LandCoverChangeItem> areasHarvested = new ArrayList<>();
+ private Map<WoodType, Double> woodAmount = new HashMap<>();
+
+ public List<LandCoverChangeItem> getAreasHarvested() {
+ return areasHarvested;
+ }
+
+ public void addAreasHarvested(LandCoverChangeItem lccItem) {
+ areasHarvested.add(lccItem);
+ }
+
+ public Map<WoodType, Double> getWoodAmount() {
+ return woodAmount;
+ }
+
+ public void addWoodAmount(WoodType woodType, double amount) {
+ woodAmount.merge(woodType, amount, Double::sum);
+ }
+}
diff --git a/src/ac/ed/lurg/landuse/LandCoverTile.java b/src/ac/ed/lurg/landuse/LandCoverTile.java
index 246adc27..3fb40566 100644
--- a/src/ac/ed/lurg/landuse/LandCoverTile.java
+++ b/src/ac/ed/lurg/landuse/LandCoverTile.java
@@ -39,9 +39,6 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
public void addArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) {
LandKey key = new LandKey(lcType, lpType, age);
- if (Double.isNaN(a)) {
- int foo = 1;
- }
addArea(key, a);
}
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index cdb21b97..d343e173 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -5,6 +5,7 @@ import java.util.Collection;
import java.util.HashMap;
import java.util.Map;
import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.forestry.WoodHarvestItem;
import ac.ed.lurg.types.*;
import ac.ed.lurg.utils.Interpolator;
import ac.sac.raster.InterpolatingRasterItem;
@@ -16,7 +17,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
private LandCoverTile landCoverAreas = new LandCoverTile();
private Map<WoodType, Double> forestFractions = new HashMap<>();
- private Map<LccKey, Double> landCoverChange = new HashMap<LccKey, Double>(); // previous timestep LCC
+ private Map<LccKey, Double> landCoverChange = new HashMap<LccKey, Double>(); // land cover change for previous timestep
+ private WoodHarvestItem woodHarvests; // forest rotations
private double protectedFraction; // protected fraction of total area
public LandUseItem() {}
@@ -112,21 +114,40 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
if (area == 0)
return 0;
- double prevFrom = getLandCoverArea(fromType, LandProtectionType.CONVERTIBLE);
+ double prevFromArea = getLandCoverArea(fromType, LandProtectionType.CONVERTIBLE);
double areaMoved;
- areaMoved = Math.min(area, prevFrom);
+ areaMoved = Math.min(area, prevFromArea);
if (Math.abs(area - areaMoved) / area > 1e-6) {
throw new RuntimeException("Move unprotected area shortfall. Should not happen due to GAMS constraint.");
}
-
- // Note if toType == fromType, area will be removed from fromType and added back in with age 0
- landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMoved);
if (ModelConfig.IS_CALIBRATION_RUN) {
- landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, areaMoved); // preserve age distribution
+ // Preserve age distribution of source land cover.
+ // This is mainly so that forest <> non-forest transitions don't lose trees.
+ if (fromType.equals(LandCoverType.NATURAL) && toType.equals(LandCoverType.TIMBER_FOREST)) {
+ double moveFraction = areaMoved / prevFromArea;
+ for (int age : landCoverAreas.getAgeKeys()) {
+ double fromArea = landCoverAreas.getArea(fromType, LandProtectionType.CONVERTIBLE, age);
+ if (fromArea == 0)
+ continue;
+ double partialArea = fromArea * moveFraction;
+ landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, age, partialArea);
+ landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, age, partialArea);
+ }
+ } else {
+
+ double areaMovedCalib = areaMoved;
+ if (fromType.isNatural()) {
+ areaMovedCalib *= 0.95; // keep some area during calibration to preserve age distribution
+ }
+ landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMovedCalib);
+ landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, areaMovedCalib);
+ }
+
} else {
+ landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMoved);
landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, 0, areaMoved); // new area starts at age=0
}
@@ -137,9 +158,27 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
}
public void doForestRotation(int age, double area) {
+ if (ModelConfig.IS_CALIBRATION_RUN)
+ return;
landCoverAreas.removeArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age, area);
landCoverAreas.addArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, 0, area);
}
+
+ public void setForestFraction(WoodType woodType, double fraction) {
+ forestFractions.put(woodType, fraction);
+ }
+
+ public double getForestFraction(WoodType woodType) {
+ return forestFractions.get(woodType);
+ }
+
+ public void setWoodHarvests(WoodHarvestItem woodHarvests) {
+ this.woodHarvests = woodHarvests;
+ }
+
+ public WoodHarvestItem getWoodHarvests() {
+ return woodHarvests;
+ }
public void resetLccAreas() {
landCoverChange.clear();
diff --git a/src/ac/ed/lurg/landuse/WoodUsageReader.java b/src/ac/ed/lurg/landuse/WoodUsageReader.java
index 5c95651d..a9d0862d 100644
--- a/src/ac/ed/lurg/landuse/WoodUsageReader.java
+++ b/src/ac/ed/lurg/landuse/WoodUsageReader.java
@@ -50,8 +50,8 @@ public class WoodUsageReader {
for (Map<String, String> row : rowList) {
String countryCode = row.get("Iso3");
String woodTypeName = row.get("Type");
- double production = Double.parseDouble(row.get("Production"));
- double netImport = Double.parseDouble(row.get("NetImports"));
+ double production = Double.parseDouble(row.get("Production")) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR;
+ double netImport = Double.parseDouble(row.get("NetImports")) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR;
SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
WoodType woodType = WoodType.getForName(woodTypeName);
CompositeCountry cc = CountryManager.getInstance().getForSingleCountry(country);
diff --git a/src/ac/ed/lurg/types/CropType.java b/src/ac/ed/lurg/types/CropType.java
index 4441e9f7..06457535 100644
--- a/src/ac/ed/lurg/types/CropType.java
+++ b/src/ac/ed/lurg/types/CropType.java
@@ -11,18 +11,18 @@ import ac.ed.lurg.utils.LogWriter;
public enum CropType {
- WHEAT("WheatBarleyOats", "wheat", true, false, 9.5, ModelConfig.INITAL_PRICE_WHEAT, 0.35),
- MAIZE("MaizeMilletSorghum", "maize", true, false, 5.1, ModelConfig.INITAL_PRICE_MAIZE, 0.35),
- RICE("Rice (Paddy Equivalent)", "rice", true, false, 8.3, ModelConfig.INITAL_PRICE_RICE, 0.35),
- OILCROPS("Oilcrops", "oilcrops", true, false, 4.4, ModelConfig.INITAL_PRICE_OILCROPS, 0.25),
- PULSES("Pulses", "pulses", true, false, 10.8, ModelConfig.INITAL_PRICE_PULSES, 0.35),
- STARCHY_ROOTS("Starchy Roots", "starchyRoots", true, false, 14.3, ModelConfig.INITAL_PRICE_STARCHYROOTS, 0.15),
+ WHEAT("WheatBarleyOats", "wheat", true, false, 9.5, ModelConfig.INITAL_PRICE_WHEAT, 0.4),
+ MAIZE("MaizeMilletSorghum", "maize", true, false, 5.1, ModelConfig.INITAL_PRICE_MAIZE, 0.4),
+ RICE("Rice (Paddy Equivalent)", "rice", true, false, 8.3, ModelConfig.INITAL_PRICE_RICE, 0.4),
+ OILCROPS("Oilcrops", "oilcrops", true, false, 4.4, ModelConfig.INITAL_PRICE_OILCROPS, 0.3),
+ PULSES("Pulses", "pulses", true, false, 10.8, ModelConfig.INITAL_PRICE_PULSES, 0.5),
+ STARCHY_ROOTS("Starchy Roots", "starchyRoots", true, false, 14.3, ModelConfig.INITAL_PRICE_STARCHYROOTS, 0.2),
ENERGY_CROPS("Energy crops", "energycrops", true, false, 5, ModelConfig.INITAL_PRICE_ENERGYCROPS, 0.2),
SETASIDE("setaside", "setaside", false, false, 0, 0, 0),
- MONOGASTRICS("Monogastrics", "monogastrics", true, true, 3.1, ModelConfig.INITAL_PRICE_MONOGASTRICS, 0.1),
- RUMINANTS("Ruminants", "ruminants", true, true, 2.2, ModelConfig.INITAL_PRICE_RUMINANTS, 0.1),
- FRUITVEG("FruitVeg", "fruitveg", true, false, 8.9, ModelConfig.INITAL_PRICE_FRUITVEG, 0.1),
- SUGAR("Sugar", "sugar", true, false, 0.1, ModelConfig.INITAL_PRICE_SUGAR, 0.4),
+ MONOGASTRICS("Monogastrics", "monogastrics", true, true, 3.1, ModelConfig.INITAL_PRICE_MONOGASTRICS, 0.2),
+ RUMINANTS("Ruminants", "ruminants", true, true, 2.2, ModelConfig.INITAL_PRICE_RUMINANTS, 0.2),
+ FRUITVEG("FruitVeg", "fruitveg", true, false, 8.9, ModelConfig.INITAL_PRICE_FRUITVEG, 0.2),
+ SUGAR("Sugar", "sugar", true, false, 0.1, ModelConfig.INITAL_PRICE_SUGAR, 0.5),
PASTURE("pasture", "pasture", false, false, 0, 0, 0); // confusing having a land cover and a crop type. Needed here for yields (but not used for cropland area fractions).
private String faoName;
diff --git a/src/ac/ed/lurg/types/WoodType.java b/src/ac/ed/lurg/types/WoodType.java
index 4e4e7c68..73c6f370 100644
--- a/src/ac/ed/lurg/types/WoodType.java
+++ b/src/ac/ed/lurg/types/WoodType.java
@@ -5,7 +5,7 @@ import java.util.Map;
import ac.ed.lurg.ModelConfig;
public enum WoodType {
- IND_ROUNDWOOD("roundwood", ModelConfig.IND_ROUNDWOOD_DEMAND_ELASTICITY, 40, ModelConfig.INIT_ROUNDWOOD_PRICE),
+ IND_ROUNDWOOD("roundwood", ModelConfig.IND_ROUNDWOOD_DEMAND_ELASTICITY, 30, ModelConfig.INIT_ROUNDWOOD_PRICE),
FUELWOOD("fuelwood", ModelConfig.FUELWOOD_DEMAND_ELASTICITY, 10, ModelConfig.INIT_FUELWOOD_PRICE);
private String name;
--
GitLab