From e919dcfb6d62cf6d769efe61252aa8c3048d0522 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Thu, 8 Jul 2021 21:04:04 +0100
Subject: [PATCH] Further work on forest growth mechanism.
---
GAMS/IntExtOpt.gms | 32 ++--
debug_config.properties | 6 +-
src/ac/ed/lurg/ModelConfig.java | 9 +-
src/ac/ed/lurg/ModelMain.java | 31 ++-
src/ac/ed/lurg/country/CountryAgent.java | 58 +++++-
.../ed/lurg/country/CountryAgentManager.java | 7 +-
.../lurg/country/gams/GamsCountryInput.java | 9 +-
.../lurg/country/gams/GamsLocationInput.java | 8 +-
.../country/gams/GamsLocationOptimiser.java | 23 ++-
.../lurg/country/gams/GamsLocationOutput.java | 9 +-
.../country/gams/GamsRasterOptimiser.java | 98 ++++------
src/ac/ed/lurg/forestry/ForestGrowthItem.java | 24 ++-
src/ac/ed/lurg/forestry/ForestManager.java | 178 +++++++++++++++---
src/ac/ed/lurg/forestry/ForestStandItem.java | 28 ++-
src/ac/ed/lurg/landuse/LandUseItem.java | 52 +++--
src/ac/ed/lurg/utils/DoubleMap.java | 4 +-
16 files changed, 439 insertions(+), 137 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 44d288cc..8c61305e 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -70,7 +70,7 @@
SCALAR woodDemand;
SCALAR carbonPrice price of carbon - $1000 per tonne;
SCALAR woodPrice price of wood and timber - $1000 per tC-eq;
- SCALAR forestRotationPeriod;
+ PARAMETER forestRotationPeriod(location);
SCALAR woodMaxNetImport;
SCALAR woodMinNetImport;
SCALAR managedForestAreaMaxChangeRate;
@@ -161,7 +161,7 @@ $gdxin
totalConversionCost(land_cover, location) land cover conversion cost - $1000 per ha
newForestPotentialHarvest(location)
naturalDeforestedArea(location)
- woodHarvestRota(location) wood harvested from timber forest rotation - Mt C-eq
+ woodHarvestNat(location)
woodHarvestLuc(location)
woodExported total wood sold - Mt C-eq
woodImported
@@ -174,7 +174,7 @@ $gdxin
total_cost total cost of domestic supply including net imports;
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
- totalFeedDM, landCoverArea, landCoverChange, totalConversionCost, woodHarvestRota, woodHarvestLuc, woodExported, woodImported,
+ totalFeedDM, landCoverArea, landCoverChange, totalConversionCost, woodHarvestNat, woodHarvestLuc, woodExported, woodImported,
supply, woodSold, woodSupply, exportAmount, newForestPotentialHarvest, naturalDeforestedArea;
* POSITIVE VARIABLE A;
@@ -286,9 +286,9 @@ $gdxin
NEW_FOREST_POTENTIAL_HARVEST(location) .. newForestPotentialHarvest(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location));
- WOOD_HARVEST_NATURAL_CALC(location) .. woodHarvestNat(location) =E= naturalDeforestedArea(location) * woodYield('natural', location);
+ WOOD_HARVEST_NATURAL_CALC(location) .. woodHarvestNat(location) =E= naturalDeforestedArea(location) * woodYieldLuc('natural', 'pasture', location);
- WOOD_HARVEST_NATURAL_CONSTRAINT(location) .. naturalDeforestedArea(location) =L= landCoverChange('natural', 'natural' location);
+ WOOD_HARVEST_NATURAL_CONSTRAINT(location) .. naturalDeforestedArea(location) =L= landCoverChange('natural', 'natural', location);
WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((forested, land_cover_after), landCoverChange(forested, land_cover_after, location) * woodYieldLuc(forested, land_cover_after, location));
@@ -315,15 +315,17 @@ $gdxin
] * exp(-discountRate) -
sum((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) -
sum(import_crop, importPrices(import_crop) * importAmount(import_crop))
- } / (1 - exp(-discountRate) +
+ } / (1 - exp(-discountRate)) +
- {
- woodPrice * newForestPotentialHarvest(location) * exp(-discountRate * forestRotationPeriod) -
- sum(location, landCoverArea('timberForest', location) * forestEstablishmentCost) -
- woodPrice * woodImported
- } / (1 - exp(-discountRate * forestRotationPeriod)) -
-
- sum((land_cover, location), totalConversionCost(land_cover, location));
+ sum(location,
+ [
+ woodPrice * newForestPotentialHarvest(location) * exp(-discountRate * forestRotationPeriod(location)) -
+ landCoverArea('timberForest', location) * forestEstablishmentCost
+ ] / (1 - exp(-discountRate * forestRotationPeriod(location)))
+ ) -
+ sum((land_cover, location), totalConversionCost(land_cover, location)) +
+ woodPrice * woodSold -
+ woodPrice * woodImported;
@@ -385,8 +387,8 @@ $gdxin
netCarbonFlux = SUM(location, carbonFlux.L(location));
netWoodImport = woodImported.L - woodExported.L;
- netWoodImport = 0;
- totalWoodHarvest = sum(location, woodHarvestRota.L(location) + woodHarvestLuc.L(location));
+* netWoodImport = 0;
+ totalWoodHarvest = sum(location, woodHarvestNat.L(location) + woodHarvestLuc.L(location)) + managedWoodHarvest;
Scalar totalCostsLessLU;
diff --git a/debug_config.properties b/debug_config.properties
index 9b72a2af..abec5b0e 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -17,10 +17,10 @@ GENERATE_NEW_YIELD_CLUSTERS=false
YIELD_FILENAME=yield.out
DEBUG_LIMIT_COUNTRIES=true
-DEBUG_COUNTRY_NAME=Brazil
-GAMS_COUNTRY_TO_SAVE=Brazil
+DEBUG_COUNTRY_NAME=Central America
+GAMS_COUNTRY_TO_SAVE=Central America
-INIT_WOOD_PRICE=0.0
+INIT_WOOD_PRICE=0.05
INIT_WOOD_STOCK=1000
INIT_CARBON_PRICE=0.0
FOREST_ROTATION_PERIOD=30
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index f62a7b20..b372bdce 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -250,9 +250,9 @@ public class ModelConfig {
public static final String FOREST_DIR = SPATIAL_DATA_DIR + File.separator + "forestry";
//public static final String CARBON_FLUX_FILE = FOREST_DIR + File.separator + "carbon_flux_";
public static final String WOOD_YIELD_FILENAME = "wood_yield.out";
- public static final String TIMBER_FOREST_GROWTH_FILENAME = "x1.out";
- public static final String CARBON_FOREST_GROWTH_FILENAME = "x2.out";
- public static final String NATURAL_FOREST_GROWTH_FILENAME = "x3.out";
+ public static final String TIMBER_FOREST_GROWTH_FILENAME = "growth_timber.out";
+ public static final String CARBON_FOREST_GROWTH_FILENAME = "growth_carbon.out";
+ public static final String NATURAL_FOREST_GROWTH_FILENAME = "growth_natural.out";
// Output
public static final String LAND_COVER_OUTPUT_FILE = OUTPUT_DIR + File.separator + "lc.txt";
@@ -448,5 +448,6 @@ public class ModelConfig {
public static final String WOOD_NET_IMPORTS_FILENAME = getProperty("WOOD_NET_IMPORTS_FILENAME", "wood_net_imports.csv");
public static final String WOOD_NET_IMPORTS_FILE = getProperty("WOOD_NET_IMPORTS_FILE", DATA_DIR + File.separator + WOOD_NET_IMPORTS_FILENAME);
public static final double WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 2.688e-7); // m3 to MtC-eq p.16 [https://doi.org/10.5194/gmd-13-5425-2020]
- public static final int FOREST_INIT_AGE = getIntProperty("FOREST_INIT_AGE", 60);
+ public static final int FOREST_INIT_AGE = getIntProperty("FOREST_INIT_AGE", 30);
+ public static final double DISCOUNT_RATE = getDoubleProperty("DISCOUNT_RATE", 0.08);
}
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 487cfa38..88814413 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -28,6 +28,7 @@ import ac.ed.lurg.demand.DemandManagerFromFile;
import ac.ed.lurg.demand.DemandManagerSSP;
import ac.ed.lurg.demand.ElasticDemandManager;
import ac.ed.lurg.forestry.CarbonForestGrowthDataReader;
+import ac.ed.lurg.forestry.ForestGrowthItem;
import ac.ed.lurg.forestry.ForestGrowthRasterSet;
import ac.ed.lurg.forestry.NaturalForestGrowthDataReader;
import ac.ed.lurg.forestry.TimberForestGrowthDataReader;
@@ -85,6 +86,8 @@ public class ModelMain {
private IrrigationRasterSet currentIrrigationData;
private RasterSet<LandUseItem> globalLandUseRaster;
private RasterSet<IntegerRasterItem> clusterIdRaster;
+
+ private ForestGrowthRasterSet forestGrowthRaster;
public static void main(String[] args) {
ModelMain theModel = new ModelMain();
@@ -118,6 +121,10 @@ public class ModelMain {
globalLandUseRaster = new RasterSet<LandUseItem>(desiredProjection);
internationalMarket = new InternationalMarket();
+
+ forestGrowthRaster = new ForestGrowthRasterSet(desiredProjection);
+ getForestGrowthData(new Timestep(ModelConfig.START_TIMESTEP));
+
createCountryAgents(compositeCountryManager.getAll());
}
@@ -157,10 +164,12 @@ public class ModelMain {
CarbonFluxRasterSet currentCarbonFluxData = getCarbonFluxData(timestep);
WoodYieldRasterSet currentWoodYieldData = getWoodYieldData(timestep);
+ getForestGrowthData(timestep);
+
DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = new ConversionCostReader().read();
countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData,
- currentWoodYieldData, conversionCosts, carbonDemandIncrease);
+ currentWoodYieldData, conversionCosts, carbonDemandIncrease, forestGrowthRaster);
internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep);
// loop for iterations. Could check within a tolerance using internationalMarket.findMaxPriceDiff, not doing so as doesn't find a solution due to inelastic consumption
@@ -679,11 +688,18 @@ public class ModelMain {
new LandCoverReader(initialLC).getRasterDataFromFile(ModelConfig.INITAL_LAND_COVER_FILE); // Land cover fractions
RasterSet<LandUseItem> landUseRaster = new RasterSet<LandUseItem>(initialLC.getHeaderDetails());
+
+ int totalMissing = 0;
for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
- landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue()));
+ ForestGrowthItem growthFunction = forestGrowthRaster.get(entry.getKey());
+ landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue(), growthFunction));
+ if (growthFunction == null)
+ totalMissing +=1;
}
+ LogWriter.printlnWarning("Missing forest growth functions = " + totalMissing);
+
return landUseRaster;
}
@@ -720,13 +736,12 @@ public class ModelMain {
return woodYieldData;
}
- private ForestGrowthRasterSet getForestGrowthData(Timestep timestep) {
- ForestGrowthRasterSet forestGrowthData = new ForestGrowthRasterSet(desiredProjection);
+ private void getForestGrowthData(Timestep timestep) {
LogWriter.println("Reading forest growth data");
- new TimberForestGrowthDataReader(forestGrowthData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + ModelConfig.TIMBER_FOREST_GROWTH_FILENAME);
- new CarbonForestGrowthDataReader(forestGrowthData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + ModelConfig.CARBON_FOREST_GROWTH_FILENAME);
- new NaturalForestGrowthDataReader(forestGrowthData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + ModelConfig.NATURAL_FOREST_GROWTH_FILENAME);
- return forestGrowthData;
+ String rootDir = timestep.getYearSubDir(ModelConfig.FOREST_DIR);
+ new TimberForestGrowthDataReader(forestGrowthRaster).getRasterDataFromFile(rootDir + File.separator + ModelConfig.TIMBER_FOREST_GROWTH_FILENAME);
+ new CarbonForestGrowthDataReader(forestGrowthRaster).getRasterDataFromFile(rootDir + File.separator + ModelConfig.CARBON_FOREST_GROWTH_FILENAME);
+ new NaturalForestGrowthDataReader(forestGrowthRaster).getRasterDataFromFile(rootDir + File.separator + ModelConfig.NATURAL_FOREST_GROWTH_FILENAME);
}
/** Ugly in situ update of currentIrrigationData, better if IrrigationRasterSets were handled more immutably */
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 96a982ad..e032e2ac 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -16,6 +16,7 @@ import ac.ed.lurg.country.gams.GamsRasterInput;
import ac.ed.lurg.country.gams.GamsRasterOptimiser;
import ac.ed.lurg.country.gams.GamsRasterOutput;
import ac.ed.lurg.demand.AbstractDemandManager;
+import ac.ed.lurg.forestry.ForestGrowthItem;
import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.IrrigationItem;
@@ -108,10 +109,14 @@ public class CountryAgent extends AbstractCountryAgent {
public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData,
Map<CropType, GlobalPrice> worldPrices, double globalGen2EcIncrease, RasterSet<CarbonFluxItem> carbonFluxData,
RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
- double carbonDemandIncrease, GlobalPrice carbonPrice, GlobalPrice timberPrice) {
+ double carbonDemandIncrease, GlobalPrice carbonPrice, GlobalPrice woodPrice, RasterSet<ForestGrowthItem> forestGrowthData) {
// project demand
- calculateCountryPricesAndDemand(worldPrices, carbonPrice, timberPrice, false);
+ calculateCountryPricesAndDemand(worldPrices, carbonPrice, woodPrice, false);
+
+ updateForestGrowthData(forestGrowthData);
+ growForests();
+ calculateForestRotation(woodPrice.getExportPrice(), ModelConfig.DISCOUNT_RATE);
if (saveGamsGdxFiles && ModelConfig.PRICE_ELASTIC_DEMAND)
saveGDXFile("demand");
@@ -242,7 +247,7 @@ public class CountryAgent extends AbstractCountryAgent {
GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints,
previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates, currentCarbonPrice, carbonTradeConstraint, currentWoodDemand,
- currentTimberPrice, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData());
+ currentTimberPrice, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData(), getWoodHarvestFromRotation());
GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData,
carbonFluxData, woodYieldData, conversionCosts, countryLevelInputs);
@@ -339,4 +344,51 @@ public class CountryAgent extends AbstractCountryAgent {
public Map<WoodType, WoodUsageData> getWoodUsageData() {
return previousGamsRasterOutput.getWoodUsageData();
}
+
+ private void updateForestGrowthData(RasterSet<ForestGrowthItem> forestGrowthData) {
+ RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
+ for (Map.Entry<RasterKey, LandUseItem> entry : previousLandUses.entrySet()) {
+ RasterKey key = entry.getKey();
+ LandUseItem luItem = entry.getValue();
+ if (luItem == null)
+ continue;
+
+ ForestGrowthItem growthFunction = forestGrowthData.get(key);
+ if (growthFunction == null)
+ continue;
+
+ luItem.updateForestGrowthFunction(growthFunction);
+ }
+ }
+
+ private void growForests() {
+ RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
+ for (LandUseItem luItem : previousLandUses.values()) {
+ if (luItem == null)
+ continue;
+ luItem.growForests();
+ }
+ }
+
+ private void calculateForestRotation(double woodPrice, double discountRate) {
+ RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
+ for (LandUseItem luItem : previousLandUses.values()) {
+ if (luItem == null)
+ continue;
+
+ luItem.calculateOptimalRotation(woodPrice, discountRate);;
+ }
+ }
+
+ private double getWoodHarvestFromRotation() {
+ RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
+ double totalHarvest = 0.0;
+ for (LandUseItem luItem : previousLandUses.values()) {
+ if (luItem == null)
+ continue;
+
+ totalHarvest += luItem.getWoodHarvestFromRotation();
+ }
+ return totalHarvest;
+ }
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index 1683b9ad..5285d49c 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -15,6 +15,8 @@ import ac.ed.lurg.Timestep;
import ac.ed.lurg.country.crafty.CraftyCountryAgent;
import ac.ed.lurg.country.crafty.CraftyProdManager;
import ac.ed.lurg.demand.AbstractDemandManager;
+import ac.ed.lurg.forestry.ForestGrowthItem;
+import ac.ed.lurg.forestry.ForestGrowthRasterSet;
import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.CarbonFluxRasterSet;
import ac.ed.lurg.landuse.CropUsageData;
@@ -108,7 +110,7 @@ public class CountryAgentManager {
public void determineProductionForAll(Timestep timestep, YieldRaster yieldSurfaces, IrrigationRasterSet currentIrrigationData, double gen2Increase,
CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
- double carbonDemandIncrease) {
+ double carbonDemandIncrease, ForestGrowthRasterSet forestGrowthRaster) {
for (AbstractCountryAgent aca : countryAgents) {
aca.setCurrentTimestep(timestep);
@@ -121,10 +123,11 @@ public class CountryAgentManager {
RasterSet<IrrigationItem> irrigData = currentIrrigationData.createSubsetForKeys(countryKeys);
RasterSet<CarbonFluxItem> carbonFluxData = currentCarbonFluxData.createSubsetForKeys(countryKeys);
RasterSet<WoodYieldItem> woodYieldData = currentWoodYieldData.createSubsetForKeys(countryKeys);
+ RasterSet<ForestGrowthItem> forestGrowthData = forestGrowthRaster.createSubsetForKeys(countryKeys);
try {
ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase, carbonFluxData, woodYieldData,
- conversionCosts, carbonDemandIncrease, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice());
+ conversionCosts, carbonDemandIncrease, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice(), forestGrowthData);
// update global rasters
globalLandUseRaster.putAll(ca.getLandUses());
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 7b9f503c..df876f69 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -28,12 +28,14 @@ public class GamsCountryInput {
private CountryPrice woodPrice;
private Map<WoodType, TradeConstraint> woodTradeConstraints;
private Map<WoodType, WoodUsageData> previousWoodUsageData;
+ private double woodHarvestFromRotation;
public GamsCountryInput(CompositeCountry country, Map<CommodityType, Double> projectedDemand, Map<CropType, CountryPrice> countryPrices,
Map<CropType, TradeConstraint> importConstraints, Map<CropType, CropUsageData> previousCropUsageData,
Map<CommodityType, Map<CropType, Double>> minDemandFracts, Map<CropType, Double> subsidyRates,
CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint, Map<WoodType, Double> woodDemand,
- CountryPrice woodPrice, Map<WoodType, TradeConstraint> woodTradeConstraints, Map<WoodType, WoodUsageData> previousWoodUsageData) {
+ CountryPrice woodPrice, Map<WoodType, TradeConstraint> woodTradeConstraints, Map<WoodType, WoodUsageData> previousWoodUsageData,
+ double woodHarvestFromRotation) {
super();
this.country = country;
this.projectedDemand = projectedDemand;
@@ -48,6 +50,7 @@ public class GamsCountryInput {
this.woodPrice = woodPrice;
this.woodTradeConstraints = woodTradeConstraints;
this.previousWoodUsageData = previousWoodUsageData;
+ this.woodHarvestFromRotation = woodHarvestFromRotation;
}
public CompositeCountry getCountry() {
@@ -124,6 +127,10 @@ public class GamsCountryInput {
return previousWoodUsageData;
}
+ public double getWoodHarvestFromRotation() {
+ return woodHarvestFromRotation;
+ }
+
public void setTradeConstraints(Map<CropType, TradeConstraint> tradeConstraints) {
this.tradeConstraints = tradeConstraints;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index a5ac0d8e..b89f0200 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -21,11 +21,12 @@ public class GamsLocationInput {
private Map<Integer, ? extends WoodYieldItem> woodYields;
private DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts;
private GamsCountryInput countryInput;
+ private Map<Integer, Double> rotationPeriods;
public GamsLocationInput(Timestep timestep, Map<Integer, ? extends YieldResponsesItem> yields, Map<Integer, ? extends LandUseItem> previousLandUse,
Map<Integer, ? extends IrrigationItem> irrigationCosts, Map<Integer, ? extends CarbonFluxItem> carbonFluxes,
Map<Integer, ? extends WoodYieldItem> woodYields,
- DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts, GamsCountryInput countryInput) {
+ DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts, GamsCountryInput countryInput, Map<Integer, Double> rotationPeriods) {
super();
this.timestep = timestep;
this.yields = yields;
@@ -35,6 +36,7 @@ public class GamsLocationInput {
this.woodYields = woodYields;
this.conversionCosts = conversionCosts;
this.countryInput = countryInput;
+ this.rotationPeriods = rotationPeriods;
}
public Map<Integer, ? extends YieldResponsesItem> getYields() {
@@ -65,6 +67,10 @@ public class GamsLocationInput {
return countryInput;
}
+ public Map<Integer, Double> getRotationPeriods() {
+ return rotationPeriods;
+ }
+
public Timestep getTimestep() {
return timestep;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 6ef2102c..4b62389c 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -454,12 +454,20 @@ public class GamsLocationOptimiser {
setGamsParamValue(conversionCostP.addRecord(v), cost, 5);
}
}
-
- addScalar(inDB, "forestRotationPeriod", ModelConfig.FOREST_ROTATION_PERIOD, 3);
+
+ GAMSParameter forestRotationP = inDB.addParameter("forestRotationPeriod", 1);
+ for (Map.Entry<Integer, Double> entry : inputData.getRotationPeriods().entrySet()) {
+ int locationId = entry.getKey();
+ double rotPeriod = entry.getValue();
+ String locString = Integer.toString(locationId);
+ setGamsParamValue(forestRotationP.addRecord(locString), rotPeriod, 0);
+ }
+
double forestChangeRate = (ModelConfig.IS_CALIBRATION_RUN) ? 0.0 : ModelConfig.FOREST_MAX_CHANGE;
addScalar(inDB, "managedForestAreaMaxChangeRate", forestChangeRate, 3);
addScalar(inDB, "forestEstablishmentCost", ModelConfig.FOREST_ESTABLISHMENT_COST, 3);
addScalar(inDB, "naturalAreaValue", ModelConfig.NATURAL_AREA_VALUE, 3);
+ addScalar(inDB, "managedWoodHarvest", countryInput.getWoodHarvestFromRotation(), 3);
}
private void addScalar(GAMSDatabase gamsDb, String recordName, double val, int places) {
@@ -643,8 +651,17 @@ public class GamsLocationOptimiser {
}
WoodUsageData woodUsageData = new WoodUsageData(totalWoodHarvest, netWoodImport);
+
+ Map<Integer, Double> naturalForestCleared = new HashMap<Integer, Double>();
+ GAMSVariable varNatForestCleared = outDB.getVariable("naturalDeforestedArea");
+ for (GAMSVariableRecord rec : varNatForestCleared) {
+ String locationName = rec.getKeys()[0];
+ int locId = Integer.parseInt(locationName);
+ double areaCleared = rec.getLevel();
+ naturalForestCleared.put(locId, areaCleared);
+ }
- GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, gamsLandCoverChanges, netCarbonFlux, newWoodUsageMap);
+ GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, gamsLandCoverChanges, netCarbonFlux, newWoodUsageMap, naturalForestCleared);
return results;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 07637c1b..1d1b177f 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -23,13 +23,15 @@ public class GamsLocationOutput {
ArrayList<LandCoverChangeItem> gamsLandCoverChanges;
private double netCarbonFlux;
private Map<WoodType, WoodUsageData> woodUsageData;
+ private Map<Integer, Double> naturalForestCleared;
public GamsLocationOutput(ModelStat status,
Map<Integer, LandUseItem> landUses,
Map<CropType, CropUsageData> cropUsageData,
TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChange,
ArrayList<LandCoverChangeItem> gamsLandCoverChanges,
- double netCarbonFlux, Map<WoodType, WoodUsageData> woodUsageData) {
+ double netCarbonFlux, Map<WoodType, WoodUsageData> woodUsageData,
+ Map<Integer, Double> naturalForestCleared) {
super();
this.status = status;
this.landUses = landUses;
@@ -38,6 +40,7 @@ public class GamsLocationOutput {
this.gamsLandCoverChanges = gamsLandCoverChanges;
this.netCarbonFlux = netCarbonFlux;
this.woodUsageData = woodUsageData;
+ this.naturalForestCleared = naturalForestCleared;
}
public ModelStat getStatus() {
@@ -66,4 +69,8 @@ public class GamsLocationOutput {
public Map<WoodType, WoodUsageData> getWoodUsageData() {
return woodUsageData;
}
+
+ public Map<Integer, Double> getNaturalForestCleared() {
+ return naturalForestCleared;
+ }
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index a724510f..856bc27c 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -131,6 +131,16 @@ public class GamsRasterOptimiser {
}
}
}
+
+ // Deforested natural areas. No change in natural area but forest stand area decreased and new stand created
+ if (gamsOutput.getNaturalForestCleared().containsKey(locId)) {
+ double natForestCleared = gamsOutput.getNaturalForestCleared().get(locId);
+ if (natForestCleared > 0) {
+ int foo =1;
+ }
+ allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.NATURAL, LandCoverType.NATURAL, natForestCleared, locId);
+ }
+
for (RasterKey key : keys) {
LandUseItem newLandUseItem = newLandUseRaster.get(key);
@@ -241,6 +251,8 @@ public class GamsRasterOptimiser {
LazyTreeMap<Integer, WoodYieldItem> aggregatedWoodYields = new LazyTreeMap<Integer, WoodYieldItem>() {
protected WoodYieldItem createValue() { return new WoodYieldItem(); }
};
+
+ Map<Integer, Double> aggregatedForestRotation = new HashMap<Integer, Double>();
int countFound = 0, countMissing = 0;
@@ -294,21 +306,6 @@ public class GamsRasterOptimiser {
//LogWriter.println(id + ", " + key + ", " + avgIrigCost.getIrrigCost() + ", from " + areaSoFar + ", " + suitableAreaThisTime + ", " + irrigCost.getIrrigCost());
}
- /*
- if (carbonFluxItem == null) {
- LogWriter.println("" + rasterInputData.getCarbonFluxes().getXCoordin(key));
- LogWriter.println("" + rasterInputData.getCarbonFluxes().getYCoordin(key));
- }
- */
-
- /*
- if (clusterId == 16) {
- int foo =1;
- } else {
- continue;
- }
- */
-
// Aggregate carbon fluxes and wood yields
for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
@@ -329,52 +326,39 @@ public class GamsRasterOptimiser {
}
}
- // Wood yield
- for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
- for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
- double woodYield;
- if (woodYieldItem != null) {
- if (woodYieldItem.checkForWoodYieldLucKeys(prevLC, newLC)) {
- woodYield = woodYieldItem.getWoodYieldLuc(prevLC, newLC);
- } else {
- continue;
- }
- } else {
- woodYield = 0.0; // if missing data assume 0
- }
+ // LUC wood yield
+ for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
+ for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
+ double wYieldThisTime = 0;
+ if (fromLc.isNatural() && !fromLc.equals(toLc))
+ wYieldThisTime = landUseItem.getWoodYield(fromLc);
+
+ double wYieldSoFar = (aggWYield.checkForWoodYieldLucKeys(fromLc, toLc)) ? aggWYield.getWoodYieldLuc(fromLc, toLc) : wYieldThisTime;
+
+ // Want to average yield over forest area rather than suitable area
+ double lcAreaSoFar = aggLandUse.getLandCoverArea(fromLc);
+ double lcAreaThisTime = landUseItem.getLandCoverArea(fromLc);
- if (!aggWYield.checkForWoodYieldLucKeys(prevLC, newLC)) { // if not seen yet
- aggWYield.setWoodYieldLuc(prevLC, newLC, woodYield);
- } else {
- aggWYield.setWoodYieldLuc(prevLC, newLC, aggregateMean(aggWYield.getWoodYieldLuc(prevLC, newLC), suitableAreaSoFar,
- woodYield, suitableAreaThisTime));
- }
+ double wYieldAgg = aggregateMean(wYieldSoFar, lcAreaSoFar, wYieldThisTime, lcAreaThisTime);
+ aggWYield.setWoodYieldLuc(fromLc, toLc, wYieldAgg);
}
}
- // Timber yield
- for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
- for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
- double timberYield;
- if (woodYieldItem != null) {
- if (woodYieldItem.checkForWoodYieldRotaKeys(prevLC, newLC)) {
- timberYield = woodYieldItem.getWoodYieldRota(prevLC, newLC);
- } else {
- continue;
- }
- } else {
- timberYield = 0.0; // if missing data assume 0
- }
-
- if (!aggWYield.checkForWoodYieldRotaKeys(prevLC, newLC)) { // if not seen yet
- aggWYield.setWoodYieldRota(prevLC, newLC, timberYield);
- } else {
- aggWYield.setWoodYieldRota(prevLC, newLC, aggregateMean(aggWYield.getWoodYieldRota(prevLC, newLC), suitableAreaSoFar,
- timberYield, suitableAreaThisTime));
- }
- }
+ // Rotation wood yield
+ for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
+ LandCoverType toLc = LandCoverType.TIMBER_FOREST;
+ double wYieldThisTime = landUseItem.getWoodYieldAtRotation();
+ double wYieldSoFar = (aggWYield.checkForWoodYieldRotaKeys(fromLc, toLc)) ? aggWYield.getWoodYieldRota(fromLc, toLc) : wYieldThisTime;
+ double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
+
+ aggWYield.setWoodYieldRota(fromLc, toLc, wYieldAgg);
}
-
+
+ double forestRotaThisTime = landUseItem.getForestRotation();
+ double forestRotaSoFar = (aggregatedForestRotation.containsKey(clusterId)) ? aggregatedForestRotation.get(clusterId) : forestRotaThisTime;
+ double forestRotaAgg = aggregateMean(forestRotaSoFar, suitableAreaSoFar, forestRotaThisTime, suitableAreaThisTime);
+ aggregatedForestRotation.put(clusterId, forestRotaAgg);
+
// Crops yields and area fractions
for (CropType crop : CropType.getNonMeatTypes()) {
@@ -452,7 +436,7 @@ public class GamsRasterOptimiser {
checkedTotalAreas(aggregatedAreas, "after");
return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts,
- aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput());
+ aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput(), aggregatedForestRotation);
}
private void logWarningWithCoord(String message, RasterKey key, YieldRaster yieldRaster, CropType crop) {
diff --git a/src/ac/ed/lurg/forestry/ForestGrowthItem.java b/src/ac/ed/lurg/forestry/ForestGrowthItem.java
index 97cdfbf5..0fc926a6 100644
--- a/src/ac/ed/lurg/forestry/ForestGrowthItem.java
+++ b/src/ac/ed/lurg/forestry/ForestGrowthItem.java
@@ -1,17 +1,31 @@
package ac.ed.lurg.forestry;
+import java.io.Serializable;
+import java.util.Collections;
+import java.util.Map;
+
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.utils.Interpolator;
import ac.sac.raster.RasterItem;
-public class ForestGrowthItem implements RasterItem {
+public class ForestGrowthItem implements RasterItem, Serializable {
+ private static final long serialVersionUID = -4822328490070047362L;
+
private DoubleMap<LandCoverType, Integer, Double> forestGrowthFunction = new DoubleMap<LandCoverType, Integer, Double>();
public void setForestGrowthFunction(LandCoverType forestType, int time, double deltaC) {
forestGrowthFunction.put(forestType, time, deltaC);
}
+ public void setDefaultForMissingData() {
+ for (LandCoverType forestType : LandCoverType.getNaturalTypes()) {
+ for (int t = 0; t <= 100; t += 5) {
+ setForestGrowthFunction(forestType, t, 0);
+ }
+ }
+ }
+
public double getGrowth(LandCoverType forestType, int time) {
int lowerTime = (time/5) * 5;
int upperTime = lowerTime + 5;
@@ -29,4 +43,12 @@ public class ForestGrowthItem implements RasterItem {
}
return cumulGrowth;
}
+
+ public Map<Integer, Double> getGrowthFunction(LandCoverType forestType) {
+ return forestGrowthFunction.getInnerMap(forestType);
+ }
+
+ public int getMaxAge(LandCoverType forestType) {
+ return Collections.max(forestGrowthFunction.getInnerMap(forestType).keySet());
+ }
}
diff --git a/src/ac/ed/lurg/forestry/ForestManager.java b/src/ac/ed/lurg/forestry/ForestManager.java
index a8ac616c..34ac9869 100644
--- a/src/ac/ed/lurg/forestry/ForestManager.java
+++ b/src/ac/ed/lurg/forestry/ForestManager.java
@@ -1,57 +1,93 @@
package ac.ed.lurg.forestry;
+import java.io.Serializable;
import java.util.ArrayList;
+import java.util.Collections;
import java.util.HashMap;
-import java.util.List;
import java.util.Map;
-import java.util.Map.Entry;
+
import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.Timestep;
-import ac.ed.lurg.landuse.LandCoverChangeItem;
-import ac.ed.lurg.landuse.LandUseItem;
-import ac.ed.lurg.landuse.WoodYieldItem;
-import ac.ed.lurg.landuse.WoodYieldRasterSet;
import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.DoubleMap;
-import ac.ed.lurg.utils.LazyTreeMap;
import ac.ed.lurg.utils.LogWriter;
-import ac.sac.raster.IntegerRasterItem;
-import ac.sac.raster.RasterKey;
-import ac.sac.raster.RasterSet;
-public class ForestManager {
+public class ForestManager implements Serializable {
+ private static final long serialVersionUID = 3093028776147903435L;
+
private Map<LandCoverType, ArrayList<ForestStandItem>> forest;
+ private Map<LandCoverType, Double> forestArea;
+ private Map<LandCoverType, Double> unprotectedForestArea;
private ForestGrowthItem growthFunction;
+ private int optimalRotation;
- ForestManager(Map<LandCoverType, Double> landCoverAreas, double initialForestedNaturalArea, ForestGrowthItem growthFunction) {
- forest = new HashMap<LandCoverType, ArrayList<ForestStandItem>>();
- initialiseForests(landCoverAreas, initialForestedNaturalArea);
+ public ForestManager(Map<LandCoverType, Double> landCoverAreas, Map<LandCoverType, Double> unprotectedAreas, double initialForestedNaturalArea, ForestGrowthItem growthFunction) {
+ this.forest = new HashMap<LandCoverType, ArrayList<ForestStandItem>>();
+ this.forestArea = new HashMap<LandCoverType, Double>();
+ this.unprotectedForestArea = new HashMap<LandCoverType, Double>();
+ for (LandCoverType forestType : LandCoverType.getNaturalTypes()) {
+ forestArea.put(forestType, landCoverAreas.get(forestType));
+ unprotectedForestArea.put(forestType, unprotectedAreas.get(forestType));
+ }
+ if (growthFunction == null) {
+ growthFunction = new ForestGrowthItem();
+ growthFunction.setDefaultForMissingData();
+// LogWriter.printlnWarning("Null growthFunction");
+ }
+ this.growthFunction = growthFunction;
+ initialiseForests(landCoverAreas, unprotectedAreas, initialForestedNaturalArea);
}
- public void initialiseForests(Map<LandCoverType, Double> landCoverAreas, double initialForestedNaturalArea) {
+ public void initialiseForests(Map<LandCoverType, Double> landCoverAreas, Map<LandCoverType, Double> unprotectedAreas, double initialForestedNaturalArea) {
for (LandCoverType forestType : LandCoverType.getManagedForestTypes()) {
ArrayList<ForestStandItem> forestStands = new ArrayList<ForestStandItem>();
+
double yield = growthFunction.getCumulGrowth(forestType, 0, ModelConfig.FOREST_INIT_AGE);
- forestStands.add(new ForestStandItem(forestType, landCoverAreas.get(forestType), ModelConfig.FOREST_INIT_AGE, yield));
- forest.put(forestType, forestStands);
+ double unprotArea = unprotectedAreas.get(forestType);
+ double protArea = landCoverAreas.get(forestType) - unprotArea;
+
+ if (unprotArea > 0)
+ forestStands.add(new ForestStandItem(forestType, unprotArea, ModelConfig.FOREST_INIT_AGE, yield, false));
+
+ if (protArea > 0)
+ forestStands.add(new ForestStandItem(forestType, protArea, ModelConfig.FOREST_INIT_AGE, yield, true));
+
+ if (!forestStands.isEmpty())
+ forest.put(forestType, forestStands);
}
ArrayList<ForestStandItem> naturalStands= new ArrayList<ForestStandItem>();
+
+ double protectedFract = (landCoverAreas.get(LandCoverType.NATURAL) - unprotectedAreas.get(LandCoverType.NATURAL)) / landCoverAreas.get(LandCoverType.NATURAL);
double natForestedArea = initialForestedNaturalArea;
+ double natForestedProtectedArea = natForestedArea * protectedFract;
double natOtherArea = landCoverAreas.get(LandCoverType.NATURAL) - initialForestedNaturalArea;
-
+ double natOtherProtectedArea = natOtherArea * protectedFract;
double yield = growthFunction.getCumulGrowth(LandCoverType.NATURAL, 0, ModelConfig.FOREST_INIT_AGE);
- naturalStands.add(new ForestStandItem(LandCoverType.NATURAL, natForestedArea, ModelConfig.FOREST_INIT_AGE, yield));
- naturalStands.add(new ForestStandItem(LandCoverType.NATURAL, natOtherArea, 0, 0.0));
+ if (natForestedArea - natForestedProtectedArea > 0)
+ naturalStands.add(new ForestStandItem(LandCoverType.NATURAL, natForestedArea - natForestedProtectedArea, ModelConfig.FOREST_INIT_AGE, yield, false));
+
+ if (natForestedProtectedArea > 0)
+ naturalStands.add(new ForestStandItem(LandCoverType.NATURAL, natForestedProtectedArea, ModelConfig.FOREST_INIT_AGE, yield, true));
+
+ if (natOtherArea - natOtherProtectedArea > 0)
+ naturalStands.add(new ForestStandItem(LandCoverType.NATURAL, natOtherArea - natOtherProtectedArea, 0, 0.0, false));
- forest.put(LandCoverType.NATURAL, naturalStands);
+ if (natOtherProtectedArea > 0)
+ naturalStands.add(new ForestStandItem(LandCoverType.NATURAL, natOtherProtectedArea, 0, 0.0, true));
+
+ if (!naturalStands.isEmpty())
+ forest.put(LandCoverType.NATURAL, naturalStands);
+
+ // Initial optimal rotation
+ calculateOptimalRotation(ModelConfig.INIT_WOOD_PRICE, ModelConfig.DISCOUNT_RATE);
}
public void growForests() {
for (LandCoverType forestType : LandCoverType.getNaturalTypes()) {
ArrayList<ForestStandItem> forestStands = forest.get(forestType);
+ if (forestStands == null)
+ return;
for (ForestStandItem stand: forestStands) {
stand.updateAgeAndYield(growthFunction);
}
@@ -61,5 +97,99 @@ public class ForestManager {
public void setGrowthFunction(ForestGrowthItem growthFunction) {
this.growthFunction = growthFunction;
}
-
+
+ public void calculateOptimalRotation(double woodPrice, double discountRate) {
+ Map<Double, Integer> levMap = new HashMap<Double, Integer>(); // Land Expected Value
+ int maxAge = growthFunction.getMaxAge(LandCoverType.TIMBER_FOREST);
+
+ for (int age = 0; age <= maxAge; age++) {
+ double yield = growthFunction.getCumulGrowth(LandCoverType.TIMBER_FOREST, 0, age);
+ double lev = (woodPrice * yield * Math.exp(-discountRate * age) - ModelConfig.FOREST_ESTABLISHMENT_COST) / (1 - Math.exp(-discountRate * age));
+ levMap.put(lev, age);
+ }
+ // TODO potential edge case with tied LEV
+ double maxLev = Collections.max(levMap.keySet());
+ optimalRotation = levMap.get(maxLev);
+ }
+
+ public int getOptimalRotation() {
+ return optimalRotation;
+ }
+
+ public double harvestWood() {
+ ArrayList<ForestStandItem> timberForest = forest.get(LandCoverType.TIMBER_FOREST);
+ if (timberForest == null)
+ return 0;
+
+ double woodHarvested = 0;
+ for (ForestStandItem stand : timberForest) {
+ if (stand.getAge() >= optimalRotation) {
+ woodHarvested += stand.harvestWood();
+ }
+ }
+ return woodHarvested;
+ }
+
+ public void moveForestArea(LandCoverType toType, LandCoverType fromType, double area) {
+
+ ArrayList<ForestStandItem> toForest = forest.get(toType);
+ ArrayList<ForestStandItem> fromForest = forest.get(fromType);
+
+ // Edge case where both types are not forest
+ if (toForest == null && fromForest == null) {
+ return;
+ } else
+
+ // Moving from forest to non-forest i.e. removing forest
+ if (fromForest != null) {
+ double totalCurrentArea = unprotectedForestArea.get(fromType);
+ for (ForestStandItem s: fromForest) {
+ if (s.isProtected())
+ continue;
+ double areaToRemove = s.getArea() / totalCurrentArea * area;
+ s.removeArea(areaToRemove);
+ }
+ unprotectedForestArea.put(fromType, unprotectedForestArea.get(fromType) - area);
+ forestArea.put(fromType, forestArea.get(fromType) - area);
+ }
+
+ // Moving from non-forest to forest i.e. adding new forest
+ if (toForest != null) {
+ ForestStandItem newStand = new ForestStandItem(toType, area, 0, 0, false);
+ forest.get(toType).add(newStand);
+ unprotectedForestArea.put(toType, unprotectedForestArea.get(toType) + area);
+ forestArea.put(toType, forestArea.get(toType) + area);
+ }
+ }
+
+ public double getAverageYield(LandCoverType forestType) {
+ ArrayList<ForestStandItem> stands = forest.get(forestType);
+ if (stands == null)
+ return 0;
+ Double totalWoodStock = stands.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getWoodStock()).sum();
+ if (totalWoodStock == null || totalWoodStock == 0) {
+ return 0.0;
+ } else {
+ return totalWoodStock / unprotectedForestArea.get(forestType);
+ }
+ }
+
+ public double getYieldAtRotation() {
+ return growthFunction.getCumulGrowth(LandCoverType.TIMBER_FOREST, 0, optimalRotation);
+ }
+
+ public void runAreaCheck() {
+ for (Map.Entry<LandCoverType, ArrayList<ForestStandItem>> entry : forest.entrySet()) {
+ LandCoverType forestType = entry.getKey();
+ ArrayList<ForestStandItem> stands = entry.getValue();
+ double totalArea = stands.stream().mapToDouble(o -> o.getArea()).sum();
+ double totalUnprotectedArea = stands.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
+ if (totalArea != forestArea.get(forestType)) {
+ LogWriter.printlnError("ForestManager: total forest area and stand area different: " + forestArea.get(forestType) + ", " + totalArea);
+ }
+ if (totalUnprotectedArea != unprotectedForestArea.get(forestType)) {
+ LogWriter.printlnError("ForestManager: total forest area and stand area different: " + unprotectedForestArea.get(forestType) + ", " + totalUnprotectedArea);
+ }
+ }
+ }
}
diff --git a/src/ac/ed/lurg/forestry/ForestStandItem.java b/src/ac/ed/lurg/forestry/ForestStandItem.java
index f27b9e76..fee56e4d 100644
--- a/src/ac/ed/lurg/forestry/ForestStandItem.java
+++ b/src/ac/ed/lurg/forestry/ForestStandItem.java
@@ -1,20 +1,26 @@
package ac.ed.lurg.forestry;
+import java.io.Serializable;
+
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.LogWriter;
-public class ForestStandItem {
+public class ForestStandItem implements Serializable {
+ private static final long serialVersionUID = -5437305221519749016L;
+
private LandCoverType forestType;
private double area;
private int age;
private double yield;
+ private boolean isProtected;
- public ForestStandItem(LandCoverType forestType, double area, int age, double yield) {
+ public ForestStandItem(LandCoverType forestType, double area, int age, double yield, boolean isProtected) {
this.forestType = forestType;
this.area = area;
this.age = age;
this.yield = yield;
+ this.isProtected = isProtected;
}
public LandCoverType getForestType() {
@@ -29,10 +35,18 @@ public class ForestStandItem {
return age;
}
+ public void setAge(int age) {
+ this.age = age;
+ }
+
public double getYield() {
return yield;
}
+ public boolean isProtected() {
+ return isProtected;
+ }
+
public void removeArea(double a) {
if (a > area) {
area -= Math.min(area, a);
@@ -47,5 +61,15 @@ public class ForestStandItem {
yield += forestGrowthFunction.getGrowth(forestType, age);
}
+
+ public double harvestWood() {
+ double harvest = getArea() * getYield();
+ setAge(0);
+ return harvest;
+ }
+
+ public double getWoodStock() {
+ return area * yield;
+ }
}
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index 59372649..8898d1df 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.ForestGrowthItem;
import ac.ed.lurg.forestry.ForestManager;
import ac.ed.lurg.types.CropToDouble;
import ac.ed.lurg.types.CropType;
@@ -23,12 +24,11 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
private double protectedArea; //protected area in Mha
private double unavailableArea; //area unavailable due to altitude etc
private double suitableArea; //Mha
- private double initialForestedNaturalArea;
private ForestManager forestManager;
public LandUseItem() {}
- public LandUseItem(LandCoverItem landCover) {
+ public LandUseItem(LandCoverItem landCover, ForestGrowthItem growthFunction) {
this();
if (landCover != null) {
for (LandCoverType lcType : LandCoverType.values())
@@ -38,7 +38,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
setCropFraction(CropType.MAIZE, 0.5);
setUnavailableArea(landCover.getUnavailableFract());
updateSuitableArea(ModelConfig.BASE_YEAR);
- initialForestedNaturalArea = landCover.getInitialForestedNaturalArea();
+ forestManager = new ForestManager(landCoverAreas, unprotectedAreas, landCover.getInitialForestedNaturalArea(), growthFunction);
}
}
@@ -50,6 +50,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
unprotectedAreas.putAll(luItemToCopy.unprotectedAreas);
protectedArea = (luItemToCopy.protectedArea);
suitableArea = (luItemToCopy.suitableArea);
+ forestManager = (luItemToCopy.forestManager);
}
public void setProtectedFract(double protFrac) {
@@ -221,9 +222,11 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
double areaMoved;
areaMoved = Math.min(area, prevFrom);
-
- unprotectedAreas.put(toType, prevTo + areaMoved);
- unprotectedAreas.put(fromType, prevFrom - areaMoved);
+
+ if (!fromType.equals(toType)) {
+ unprotectedAreas.put(toType, prevTo + areaMoved);
+ unprotectedAreas.put(fromType, prevFrom - areaMoved);
+ }
/*
if (area - areaMoved > 0) {
@@ -231,6 +234,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
}
*/
+ if (toType.isNatural() || fromType.isNatural()) {
+ forestManager.moveForestArea(toType, fromType, areaMoved);
+ }
+
return area - areaMoved;
}
@@ -401,6 +408,34 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return changes;
}
+
+ public void updateForestGrowthFunction(ForestGrowthItem newGrowthFunction) {
+ forestManager.setGrowthFunction(newGrowthFunction);
+ }
+
+ public void growForests() {
+ forestManager.growForests();
+ }
+
+ public void calculateOptimalRotation(double woodPrice, double discountRate) {
+ forestManager.calculateOptimalRotation(woodPrice, discountRate);
+ }
+
+ public double getWoodYield(LandCoverType forestType) {
+ return forestManager.getAverageYield(forestType);
+ }
+
+ public double getWoodYieldAtRotation() {
+ return forestManager.getYieldAtRotation();
+ }
+
+ public double getWoodHarvestFromRotation() {
+ return forestManager.harvestWood();
+ }
+
+ public int getForestRotation() {
+ return forestManager.getOptimalRotation();
+ }
private boolean isZeroOrNull(Double d) {
return d == null || d == 0;
@@ -481,11 +516,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return total;
}
-
-
- public double getInitialForestedNaturalArea() {
- return initialForestedNaturalArea;
- }
@Override
public String toString() {
diff --git a/src/ac/ed/lurg/utils/DoubleMap.java b/src/ac/ed/lurg/utils/DoubleMap.java
index 8a50e7c0..ed3a2d10 100644
--- a/src/ac/ed/lurg/utils/DoubleMap.java
+++ b/src/ac/ed/lurg/utils/DoubleMap.java
@@ -1,9 +1,11 @@
package ac.ed.lurg.utils;
+import java.io.Serializable;
import java.util.HashMap;
import java.util.Map;
-public class DoubleMap<K, L, V> {
+public class DoubleMap<K, L, V> implements Serializable {
+ private static final long serialVersionUID = 309214428660247740L;
// Stores a double mapped by two consecutive keys
private Map<K, Map<L, Double>> theMap;
--
GitLab