diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 883e614220e3f6facbc4050c0cdf0c7fb2620976..f7b925658b8655bfc952286dfe6b43bb4baae961 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -9,10 +9,12 @@
SET feed_crop_less_pasture(feed_crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg /;
SET import_crop(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops/;
- SET land_cover / cropland, pasture, timberForest, carbonForest, natural /;
+ SET land_cover / cropland, pasture, timberForest, carbonForest, natural, naturalNew /;
SET forest(land_cover) / timberForest, carbonForest, natural /;
+ SET natural_forest(land_cover) / carbonForest, natural /;
SET non_forest(land_cover) / cropland, pasture /;
- SET non_timber_forest_before(land_cover) / cropland, pasture, carbonForest, natural /;
+ SET non_timber_forest_before(land_cover) / cropland, pasture, carbonForest, natural, naturalNew /;
+ SET exc_natural(land_cover) / cropland, pasture, timberForest, carbonForest, naturalNew /;
ALIAS(non_timber_forest_before, non_timber_forest_after);
ALIAS (land_cover, land_cover_before);
ALIAS (land_cover, land_cover_after);
@@ -61,6 +63,7 @@
PARAMETER unrestrictedArea(land_cover, location)
PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location) carbon flux - MtC-eq per Mha;
PARAMETER woodYield(land_cover_before, land_cover_after, location) wood yield - MtC-eq per Mha;
+ PARAMETER timberYield(land_cover_before, land_cover_after, location)
SCALAR carbonPrice price of carbon - $1000 per tonne;
SCALAR woodPrice price of wood and timber - $1000 per tC-eq;
PARAMETER timberForestYield(location);
@@ -75,7 +78,7 @@ $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
$load previousLandCoverArea, carbonFluxRate, woodYield, carbonPrice, woodPrice
-$load minimumLandCoverArea, conversionCost, timberForestYield
+$load minimumLandCoverArea, conversionCost, timberForestYield, timberYield
$gdxin
SCALAR delta "use to smooth power function see 7.5 www.gams.com/dd/docs/solversconopt.pdf" / 0.00000000001 /;
@@ -117,6 +120,10 @@ $gdxin
SCALAR forestExpansionMaxRate / 1 /;
SCALAR restrictedDeforestationCost / 10000000 /;
unrestrictedArea(land_cover, location) = previousLandCoverArea(land_cover, location) - minimumLandCoverArea(land_cover, location);
+ woodYield('carbonForest', 'naturalNew', location) = woodYield('carbonForest', 'timberForest', location);
+ woodYield('natural', 'naturalNew', location) = woodYield('natural', 'pasture', location);
+ conversionCost(land_cover, 'naturalNew') = conversionCost(land_cover, 'natural');
+ conversionCost('natural', 'naturalNew') = 0.05;
VARIABLES
cropArea(crop, location) total area for crops - Mha
@@ -134,7 +141,7 @@ $gdxin
landCoverChange(land_cover_before, land_cover_after, location) land cover change
restrictedLandCoverArea(land_cover, location)
restrictedLandCoverChange(land_cover, land_cover_after, location)
- prematureDeforestationCost
+ prematureDeforestationCost(location)
totalLandCoverArea(land_cover, location)
totalLandCoverChange(land_cover, land_cover_after, location)
woodHarvest(location) total wood harvested
@@ -143,7 +150,7 @@ $gdxin
timberForwardContract(location)
totalConversionCost(location)
totalFeedDM
- A "artificial variable for debugging https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
+* A "artificial variable for debugging https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
total_cost total cost of domestic supply including net imports;
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
@@ -152,7 +159,7 @@ $gdxin
landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost,
restrictedLandCoverArea, restrictedLandCoverChange, totalLandCoverArea, totalLandCoverChange;
- POSITIVE VARIABLE A;
+* POSITIVE VARIABLE A;
EQUATIONS
UNIT_COST_EQ(crop, location) cost per area - $1000 per ha or $billion per Mha
@@ -182,7 +189,7 @@ $gdxin
RESTRICTED_LAND_COVER_CHANGE_CALC(land_cover, location)
RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location)
FOREST_EXPANSION_CONSTRAINT
- PREMATURE_DEFORESTATION_COST_CALC
+ PREMATURE_DEFORESTATION_COST_CALC(location)
TOTAL_LAND_COVER_CALC(land_cover, location)
TOTAL_LAND_COVER_CHANGE_CALC(land_cover, land_cover_after, location)
* TOTAL_LAND_COVER_CONSTRAINT(location)
@@ -191,6 +198,7 @@ $gdxin
TIMBER_HARVEST_CALC(location)
TIMBER_FORWARD_CONTRACT_CALC(location)
CONVERSION_COST(location)
+ NATURAL_CONVERSION_CONSTRAINT(exc_natural, location)
COST_EQ total cost objective function;
UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E= ( baseCost(crop) +
@@ -263,24 +271,25 @@ $gdxin
RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, restrictedLandCoverChange(land_cover, land_cover_after, location)) =E= minimumLandCoverArea(land_cover, location);
+ NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) .. totalLandCoverChange(exc_natural, 'natural', location) =E= 0;
+
FOREST_EXPANSION_CONSTRAINT .. sum((non_forest, forest, location), landCoverChange(non_forest, forest, location)) =L= sum((non_forest, location), previousLandCoverArea(non_forest, location)) * forestExpansionMaxRate;
- PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, land_cover, location)$[not sameAs(forest, land_cover)], restrictedLandCoverChange(forest, land_cover, location)) * restrictedDeforestationCost;
+ PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum((forest, land_cover)$[not sameAs(forest, land_cover)], restrictedLandCoverChange(forest, land_cover, location)) * restrictedDeforestationCost;
* Timber from land cover change
- WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((non_timber_forest_before, non_timber_forest_after), landCoverChange(non_timber_forest_before, non_timber_forest_after, location) * woodYield(non_timber_forest_before, non_timber_forest_after, location));
+ WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYield(land_cover_before, land_cover_after, location));
* Timber from timberForest rotation
TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E= (sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location))) * timberForestYield(location);
* Future timber from newly planted timberForest
- TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(non_timber_forest_before, landCoverChange(non_timber_forest_before, 'timberForest', location) * woodYield(non_timber_forest_before, 'timberForest', location)) +
- (landCoverChange('timberForest', 'timberForest', location)) * woodYield('timberForest', 'timberForest', location);
+ TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * timberYield(land_cover, 'timberForest', location));
CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
* Subtracting cost of not fully grown timberForest reserved under minimumLandCoverArea
- CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after));
+ CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * (conversionCost(land_cover_before, land_cover_after) + ((previousLandCoverArea('natural', location) / suitableLandArea(location)) * 0.5 + 0.1)));
COST_EQ .. total_cost =E=
(
@@ -290,7 +299,7 @@ $gdxin
SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
- prematureDeforestationCost + SUM(location, totalConversionCost(location)) + A * 1000000
+ SUM(location, prematureDeforestationCost(location)) + SUM(location, totalConversionCost(location))
);
MODEL LAND_USE /ALL/ ;
@@ -298,8 +307,9 @@ $gdxin
irrigI.L(crop, location) = previousIrrigIntensity(crop, location);
otherIntensity.L(crop, location) = previousOtherIntensity(crop, location);
cropArea.L(crop, location) = previousCropArea(crop, location);
- landCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, location) - minimumLandCoverArea(land_cover, location);
+ landCoverArea.L(land_cover, location) = unrestrictedArea(land_cover, location);
restrictedLandCoverArea.L(land_cover, location) = minimumLandCoverArea(land_cover, location);
+ totalLandCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, location);
ruminantFeed.L(feed_crop) = previousRuminantFeed(feed_crop);
monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop);
importAmount.L(all_types) = previousImportAmount(all_types);
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index d1015b72be9666a8bbc5f6ea351c4c7e246b010c..73000f35db33750ea862a34e547c12f296279040 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -153,7 +153,7 @@ public class CountryAgent extends AbstractCountryAgent {
saveGDXFile("landuse");
previousGamsRasterOutput = result;
- forestManager.updateForestHistory(result.getNewForest(), result.getReservedDeforested(), currentTimestep);;
+ forestManager.updateForestHistory(result.getNewForest(), result.getReservedDeforested(), result.getNaturalChanges(), currentTimestep);;
return result;
}
diff --git a/src/ac/ed/lurg/country/ForestManager.java b/src/ac/ed/lurg/country/ForestManager.java
index 1f3efd83b8f7860e08e40efcb5f33129a9f4f9e5..6e0b6e45d612ae4988cacf26f70603915931cba1 100644
--- a/src/ac/ed/lurg/country/ForestManager.java
+++ b/src/ac/ed/lurg/country/ForestManager.java
@@ -16,6 +16,7 @@ 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.ed.lurg.utils.TripleMap;
import ac.sac.raster.IntegerRasterItem;
import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterSet;
@@ -23,9 +24,13 @@ import ac.ed.lurg.landuse.WoodYieldItem;
public class ForestManager {
private List<ForestRecord> forestHistory;
+ private Map<Integer, Double> naturalMeanAge;
+ private Map<Integer, Double> naturalArea;
ForestManager() {
forestHistory = new ArrayList<ForestRecord>();
+ naturalMeanAge = new HashMap<Integer, Double>();
+ naturalArea = new HashMap<Integer, Double>();
}
public void addForestRecord(ForestRecord record) {
@@ -67,13 +72,13 @@ public class ForestManager {
//LogWriter.println(key.toString());
WoodYieldItem wYItem = initWoodYieldData.get(key);
double timberYieldThisTime;
- // TODO Assuming forest was pasture before
- timberYieldThisTime = (wYItem != null) ? initWoodYieldData.get(key).getWoodYield(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST) : 0.0;
-
- unmanagedForestArea.put(clusterId, unmanagedForestAreaSoFar + luItem.getInitialForestedNaturalArea());
- otherNaturalArea.put(clusterId, otherNaturalAreaSoFar + luItem.getLandCoverArea(LandCoverType.NATURAL) - luItem.getInitialForestedNaturalArea()); // Natural is other natural + unmanaged forest so need to subtract the latter
- timberForestArea.put(clusterId, timberForestAreaSoFar + luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST));
- timberYieldTimesArea.put(clusterId, timberYieldSoFar + luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST) * timberYieldThisTime);
+ // TODO Assuming forest was forest before
+ timberYieldThisTime = (wYItem != null) ? initWoodYieldData.get(key).getTimberYield(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST) : 0.0;
+ double naturalForestFraction = (luItem.getLandCoverArea(LandCoverType.NATURAL) > 0) ? luItem.getInitialForestedNaturalArea() / luItem.getLandCoverArea(LandCoverType.NATURAL) : 0;
+ unmanagedForestArea.put(clusterId, unmanagedForestAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.NATURAL) * naturalForestFraction);
+ otherNaturalArea.put(clusterId, otherNaturalAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.NATURAL) * (1 - naturalForestFraction)); // Natural is other natural + unmanaged forest so need to subtract the latter
+ timberForestArea.put(clusterId, timberForestAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST));
+ timberYieldTimesArea.put(clusterId, timberYieldSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST) * timberYieldThisTime);
/*
if (clusterId == 3) {
@@ -87,8 +92,11 @@ public class ForestManager {
// Assume unmanaged forest is fully mature (max wood yield) and other natural is new (no wood yield).
int otherNaturalPseudoYear = ModelConfig.BASE_YEAR;
- int forestPseudoYear = ModelConfig.BASE_YEAR - ModelConfig.FOREST_LOCKIN_PERIOD;
+ int forestYearPlanted = ModelConfig.BASE_YEAR - ModelConfig.FOREST_LOCKIN_PERIOD;
+ int forestAge = ModelConfig.FOREST_LOCKIN_PERIOD;
+ int otherNaturalAge = 0;
+ /*
for (Entry<Integer, Double> entry : unmanagedForestArea.entrySet()) {
if (entry.getValue() > 0)
addForestRecord(new ForestRecord(forestPseudoYear, entry.getKey(), LandCoverType.NATURAL, entry.getValue(), 0));
@@ -98,13 +106,23 @@ public class ForestManager {
if (entry.getValue() > 0)
addForestRecord(new ForestRecord(otherNaturalPseudoYear, entry.getKey(), LandCoverType.NATURAL, entry.getValue(), 0));
}
+ */
+
+ for (Entry<Integer, Double> entry : unmanagedForestArea.entrySet()) {
+ int clusterId = entry.getKey();
+ double forestArea = entry.getValue();
+ double ntrlArea = otherNaturalArea.get(clusterId);
+ double meanAge = (forestArea + ntrlArea > 0) ? (forestAge * forestArea + otherNaturalAge * ntrlArea) / (forestArea + ntrlArea) : 0.0;
+ naturalMeanAge.put(clusterId, meanAge);
+ naturalArea.put(clusterId, forestArea + ntrlArea);
+ }
for (Entry<Integer, Double> entry : timberForestArea.entrySet()) {
int clusterId = entry.getKey();
double area = entry.getValue();
if (area > 0) {
double timberYield = timberYieldTimesArea.get(clusterId) / area;
- addForestRecord(new ForestRecord(forestPseudoYear, clusterId, LandCoverType.TIMBER_FOREST, area, timberYield));
+ addForestRecord(new ForestRecord(forestYearPlanted, clusterId, LandCoverType.TIMBER_FOREST, area, timberYield));
}
}
@@ -154,11 +172,12 @@ public class ForestManager {
return yield;
}
- // Calculates wood yield multiplier for natural areas based on age
- double naturalWoodYieldFunction(double age) {
+ // Calculates wood yield factor for natural areas based on age
+ private double naturalWoodYieldFunction(double age) {
return Math.min(age, ModelConfig.FOREST_LOCKIN_PERIOD) / ModelConfig.FOREST_LOCKIN_PERIOD;
}
+ /*
public Map<Integer, Double> getNaturalYieldFactorForTimestep(Timestep timestep) {
DoubleMap<Integer, LandCoverType, Double> ageTimesArea = new DoubleMap<Integer, LandCoverType, Double>();
DoubleMap<Integer, LandCoverType, Double> area = new DoubleMap<Integer, LandCoverType, Double>();
@@ -187,6 +206,18 @@ public class ForestManager {
return yieldFactor;
}
+ */
+
+ public Map<Integer, Double> getNaturalYieldFactorForTimestep(Timestep timestep) {
+ Map<Integer, Double> yieldFactor = new HashMap<Integer, Double>();
+ for (Entry<Integer, Double> entry : naturalMeanAge.entrySet()) {
+ int clusterId = entry.getKey();
+ double age = entry.getValue();
+ yieldFactor.put(clusterId, naturalWoodYieldFunction(age));
+ }
+ return yieldFactor;
+ }
+
public double getForestPlantedArea(int year, int location, LandCoverType forestType) {
List<ForestRecord> recordsFound = forestHistory
@@ -233,12 +264,96 @@ public class ForestManager {
return recordsFound.get(0).getPotentialYield();
}
- public void updateForestHistory(List<ForestRecord> newForest, DoubleMap<Integer, LandCoverType, Double> reservedDeforested, Timestep timestep) {
+ public void updateForestHistory(List<ForestRecord> newForest, DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
+ TripleMap<Integer, String, String, Double> naturalChanges, Timestep timestep) {
+
int yearPlanted = timestep.getYear();
for (ForestRecord item : newForest) {
- forestHistory.add(new ForestRecord(yearPlanted, item));
+ if (item.getForestType().isManagedForest())
+ forestHistory.add(new ForestRecord(yearPlanted, item));
+ }
+
+ for (ForestRecord record : forestHistory) {
+ if (record.getArea() < 0) {
+ LogWriter.printlnError("area < 0");
+ }
}
+ Map<Integer, Double> newNaturalArea = new HashMap<Integer, Double>();
+ Map<Integer, Double> clearedNaturalArea = new HashMap<Integer, Double>();
+ Map<Integer, Double> removedNaturalArea = new HashMap<Integer, Double>();
+
+ for (Entry<Integer, DoubleMap<String, String, Double>> locMap : naturalChanges.getMap().entrySet()) {
+ int locId = locMap.getKey();
+ for (Entry<String, Map<String, Double>> fromMap : locMap.getValue().getMap().entrySet()) {
+ String fromLC = fromMap.getKey();
+ for (Entry<String, Double> toMap : fromMap.getValue().entrySet()) {
+ String toLC = toMap.getKey();
+ double area = toMap.getValue();
+
+ if (!fromLC.equals("natural") && toLC.equals("naturalNew")) {
+ double areaSoFar = (newNaturalArea.containsKey(locId)) ? newNaturalArea.get(locId) : 0;
+ newNaturalArea.put(locId, areaSoFar + area);
+ }
+
+ if (fromLC.equals("natural") && toLC.equals("naturalNew")) {
+ double areaSoFar = (clearedNaturalArea.containsKey(locId)) ? clearedNaturalArea.get(locId) : 0;
+ clearedNaturalArea.put(locId, areaSoFar + area);
+ }
+
+ if (fromLC.equals("natural") && !toLC.equals("natural") && !toLC.equals("naturalNew")) {
+ double areaSoFar = (removedNaturalArea.containsKey(locId)) ? removedNaturalArea.get(locId) : 0;
+ removedNaturalArea.put(locId, areaSoFar + area);
+ }
+ }
+ }
+ }
+
+ // Update removed natural areas
+ for (Entry<Integer, Double> entry : removedNaturalArea.entrySet()) {
+ int locId = entry.getKey();
+ double removedArea = entry.getValue();
+ if (!naturalArea.containsKey(locId)) {
+ LogWriter.printlnError("ForestManager.updateForestHistory(): cannot decrement natural area as area is 0");
+ continue;
+ }
+ double currentArea = naturalArea.get(locId);
+ if (removedArea > currentArea) {
+ LogWriter.printlnError("ForestManager.updateForestHistory(): natural area to remove > current area: diff = " + (removedArea - currentArea));
+ removedArea = currentArea;
+ }
+ naturalArea.put(locId, currentArea - removedArea);
+
+ }
+
+ // Updated deforested natural areas
+ for (Entry<Integer, Double> entry : clearedNaturalArea.entrySet()) {
+ int locId = entry.getKey();
+ double clearedArea = entry.getValue();
+ if (!naturalArea.containsKey(locId)) {
+ LogWriter.printlnError("ForestManager.updateForestHistory(): cannot clear natural area as area is 0");
+ continue;
+ }
+ double currentArea = naturalArea.get(locId);
+ if (clearedArea > currentArea) {
+ LogWriter.printlnError("ForestManager.updateForestHistory(): natural area to clear > current area : diff = " + (clearedArea - currentArea));
+ }
+ double currentAge = naturalMeanAge.get(locId);
+ double newAge = ((currentArea - clearedArea) * currentAge) / currentArea;
+ naturalMeanAge.put(locId, newAge);
+ }
+
+ // Add new natural areas
+ for (Entry<Integer, Double> entry : newNaturalArea.entrySet()) {
+ int locId = entry.getKey();
+ double newArea = entry.getValue();
+ double currentArea = naturalArea.get(locId);
+ double currentAge = naturalMeanAge.get(locId);
+ double newAge = (currentArea * currentAge) / (currentArea + newArea);
+ naturalArea.put(locId, currentArea + newArea);
+ naturalMeanAge.put(locId, newAge);
+ }
+
/*
for (ForestRecord record: forestHistory) {
LogWriter.println(""+record.getLocation()+","+record.getYearPlanted()+","+record.getForestType()+","+record.getArea()+","+record.getPotentialYield());
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 078e720bded4265be2dd96c58c4dad585c356c25..e12392c6bcf6557850002d7bd084e92df37d1ad2 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -374,12 +374,15 @@ public class GamsLocationOptimiser {
double factor = entry.getValue();
for (LandCoverType toLC : LandCoverType.getConvertibleTypes()) {
WoodYieldItem item = woodYieldData.get(locId);
- double maxYield = item.getWoodYield(LandCoverType.NATURAL, toLC);
- item.setWoodYield(LandCoverType.NATURAL, toLC, maxYield * factor);
+ if (item.checkForWoodKeys(LandCoverType.NATURAL, toLC)) {
+ double maxYield = item.getWoodYield(LandCoverType.NATURAL, toLC);
+ item.setWoodYield(LandCoverType.NATURAL, toLC, maxYield * factor);
+ }
}
}
GAMSParameter woodYieldP = inDB.addParameter("woodYield", 3);
+ GAMSParameter timberYieldP = inDB.addParameter("timberYield", 3);
for (Entry<Integer, ? extends WoodYieldItem> entry : woodYieldData.entrySet()) {
Integer locationId = entry.getKey();
@@ -392,11 +395,13 @@ public class GamsLocationOptimiser {
v.add(prevLC.getName());
v.add(newLC.getName());
v.add(locString);
- if (wYield.checkForKeys(prevLC, newLC))
+ if (wYield.checkForWoodKeys(prevLC, newLC))
setGamsParamValue(woodYieldP.addRecord(v), wYield.getWoodYield(prevLC, newLC), 3);
+ if (wYield.checkForTimberKeys(prevLC, newLC))
+ setGamsParamValue(timberYieldP.addRecord(v), wYield.getTimberYield(prevLC, newLC), 3);
}
}
- }
+ }
// Timber yields
GAMSParameter timberForestYieldP = inDB.addParameter("timberForestYield", 1);
@@ -432,6 +437,9 @@ public class GamsLocationOptimiser {
v.add(landCover.getName());
v.add(locString);
double minCover = coverMap.getValue();
+ if (minCover < 0) {
+ LogWriter.printlnError("minCover < 0");
+ }
double previousCover = inputData.getPreviousLandUse().get(locationId).getUnprotectedLandCoverArea(landCover);
// Correction for floating point errors
double minCoverCorrected = Math.min(minCover, previousCover);
@@ -583,29 +591,61 @@ public class GamsLocationOptimiser {
totalCropArea+totalPastureArea, totalCropArea, totalPastureArea));
// Land cover change
- GAMSVariable varLandCoverChange = outDB.getVariable("landCoverChange");
+ GAMSVariable varLandCoverChange = outDB.getVariable("totalLandCoverChange");
TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges = new TripleMap<Integer, LandCoverType, LandCoverType, Double>();
+ DoubleMap<Integer, LandCoverType, Double> totalArea = new DoubleMap<Integer, LandCoverType, Double>();
+ DoubleMap<Integer, LandCoverType, Double> areaTimesYield = new DoubleMap<Integer, LandCoverType, Double>();
+ TripleMap<Integer, String, String, Double> naturalChanges = new TripleMap<Integer, String, String, Double>();
for (GAMSVariableRecord rec : varLandCoverChange) {
- String fromLC = rec.getKeys()[0];
- String toLC = rec.getKeys()[1];
+ String fromLcStr = rec.getKeys()[0];
+ String toLcStr = rec.getKeys()[1];
String locationName = rec.getKeys()[2];
int locId = Integer.parseInt(locationName);
-
double change = rec.getLevel();
+ if (change == 0)
+ continue;
+
+ // Changes in natural areas
+ if (fromLcStr.equals("natural") || toLcStr.equals("naturalNew")) {
+ naturalChanges.put(locId, fromLcStr, toLcStr, change);
+ }
+
+ // All land cover changes
+ fromLcStr = (fromLcStr.equals("naturalNew")) ? "natural" : fromLcStr;
+ toLcStr = (toLcStr.equals("naturalNew")) ? "natural" : toLcStr;
+ LandCoverType fromLC = LandCoverType.getForName(fromLcStr);
+ LandCoverType toLC = LandCoverType.getForName(toLcStr);
+
+
if (!fromLC.equals(toLC)) { //Important as don't want to include unchanged land cover
- landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.getForName(toLC), change);
- }
+ landCoverChanges.put(locId, fromLC, toLC, change);
+ }
+
+ // New forest area
+ if (toLC.isManagedForest()) {
+ Double minArea = inputData.getMinimumLandCover().get(locId, fromLC);
+ minArea = (minArea == null) ? 0 : minArea;
+ double newForestArea = change - minArea; // need to subtract unchanged area between forest rotations
+ if (newForestArea < 0) {
+ LogWriter.printlnError("GamsLocationOptimiser.handleResults() newForestArea < 0: " + newForestArea);
+ newForestArea = 0;
+ }
+ if (newForestArea != 0) {
+ double woodYield = (toLC.equals(LandCoverType.TIMBER_FOREST)) ? inputData.getWoodYields().get(locId).getTimberYield(fromLC, toLC) : 0.0;
+ totalArea.addTo(locId, toLC, newForestArea);
+ areaTimesYield.addTo(locId, toLC, woodYield * newForestArea);
+ }
+ }
}
// Minimum land cover additions. Need to keep track of new forest areas to restrict conversion
// Also keep track of natural areas
List<ForestRecord> newForest = new ArrayList<ForestRecord>();
- DoubleMap<Integer, LandCoverType, Double> totalArea = new DoubleMap<Integer, LandCoverType, Double>();
- DoubleMap<Integer, LandCoverType, Double> areaTimesYield = new DoubleMap<Integer, LandCoverType, Double>();
+/*
for (Entry<Integer, DoubleMap<LandCoverType, LandCoverType, Double>> locMap : landCoverChanges.getMap().entrySet()) {
int locId = locMap.getKey();
DoubleMap<LandCoverType, LandCoverType, Double> changeMap = locMap.getValue();
@@ -614,18 +654,18 @@ public class GamsLocationOptimiser {
for (Entry<LandCoverType, Double> toMap : fromMap.getValue().entrySet()) {
LandCoverType toLC = toMap.getKey();
double newArea = toMap.getValue();
- if (!toLC.isNatural() || newArea == 0.0)
+ if (!toLC.isManagedForest() || newArea == 0.0)
continue;
Double minArea = inputData.getMinimumLandCover().get(locId, fromLC);
- minArea = (minArea == null) ? 0.0 : minArea;
- double newForestArea = newArea - minArea; // subtract unchanged forest area
+ minArea = (minArea == null || !fromLC.equals(toLC)) ? 0.0 : minArea;
+ double newForestArea = newArea - minArea; // subtract unchanged forest area if going from land cover A to A
double woodYield = (toLC.equals(LandCoverType.TIMBER_FOREST)) ? inputData.getWoodYields().get(locId).getWoodYield(fromLC, toLC) : 0.0;
totalArea.addTo(locId, toLC, newForestArea);
areaTimesYield.addTo(locId, toLC, woodYield * newForestArea);
}
}
}
-
+*/
// Calculate weighted average yield
for (Entry<Integer, Map<LandCoverType, Double>> locMap : totalArea.getMap().entrySet()) {
Integer locId = locMap.getKey();
@@ -637,8 +677,7 @@ public class GamsLocationOptimiser {
newForest.add(new ForestRecord(locId, lcType, a, y));
}
}
-
-
+
// reserved Forest which was deforested
GAMSParameter varReservedDeforested = outDB.getParameter("reservedAreaDeforested");
DoubleMap<Integer, LandCoverType, Double> reservedDeforested = new DoubleMap<Integer, LandCoverType, Double>();
@@ -649,10 +688,7 @@ public class GamsLocationOptimiser {
int locId = Integer.parseInt(locationName);
double deforestedArea = rec.getValue();
-
- if (deforestedArea > 0.00001) {
- reservedDeforested.put(locId, LandCoverType.getForName(landCover), deforestedArea);
- }
+ reservedDeforested.put(locId, LandCoverType.getForName(landCover), deforestedArea);
}
// Carbon flux
@@ -661,7 +697,7 @@ public class GamsLocationOptimiser {
// Timber harvest
double totalTimberProduction = outDB.getParameter("totalTimberProduction").getFirstRecord().getValue();
- GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, newForest, reservedDeforested, netCarbonFlux, totalTimberProduction);
+ GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, newForest, reservedDeforested, naturalChanges, netCarbonFlux, totalTimberProduction);
return results;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index f0b45e884670bf8654c84ab786a4b366eca36109..8da054d5e7abd0fab47556d7bb87e125e96e505e 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -21,6 +21,7 @@ public class GamsLocationOutput {
private TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges;
private List<ForestRecord> newForest;
private DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
+ private TripleMap<Integer, String, String, Double> naturalChanges;
private double netCarbonFlux;
private double timberProduction;
@@ -30,6 +31,7 @@ public class GamsLocationOutput {
TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChange,
List<ForestRecord> newForest,
DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
+ TripleMap<Integer, String, String, Double> naturalChanges,
double netCarbonFlux, double timberProduction) {
super();
this.status = status;
@@ -38,6 +40,7 @@ public class GamsLocationOutput {
this.landCoverChanges = landCoverChange;
this.newForest = newForest;
this.reservedDeforested = reservedDeforested;
+ this.naturalChanges = naturalChanges;
this.netCarbonFlux = netCarbonFlux;
this.timberProduction = timberProduction;
}
@@ -65,6 +68,10 @@ public class GamsLocationOutput {
return reservedDeforested;
}
+ public TripleMap<Integer, String, String, Double> getNaturalChanges() {
+ return naturalChanges;
+ }
+
public double getNetCarbonFlux() {
return netCarbonFlux;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 2eb23fe87e4cc89845da9affcc20be061219a70b..cf8da9493e33ae0f73a1f70b90b8c065a7f73a51 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -54,7 +54,7 @@ public class GamsRasterOptimiser {
RasterSet<LandUseItem> newIntensityRaster = allocAreas(gamsInput.getPreviousLandUse(), gamsOutput, gamsInput.getTimestep().getYear());
return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getNewForest(),
- gamsOutput.getReservedDeforested(), gamsOutput.getNetCarbonFlux(), gamsOutput.getTimberProduction());
+ gamsOutput.getReservedDeforested(), gamsOutput.getNaturalChanges(), gamsOutput.getNetCarbonFlux(), gamsOutput.getTimberProduction());
}
private RasterSet<LandUseItem> createWithSameLandCovers(RasterSet<LandUseItem> toCopy) {
@@ -115,13 +115,15 @@ public class GamsRasterOptimiser {
DoubleMap<LandCoverType, LandCoverType, Double> landCoverChange = gamsOutput.getLandCoverChanges().getInnerMap(locId);
- // Allocate land cover change to grid cells
- for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap: landCoverChange.getMap().entrySet()) {
- LandCoverType fromLC = fromMap.getKey();
- for (Entry<LandCoverType, Double> toMap: fromMap.getValue().entrySet()) {
- LandCoverType toLC = toMap.getKey();
- double change = toMap.getValue();
- allocAllLandCropsTop(newLandUseRaster, keys, toLC, fromLC, change, locId);
+ if (landCoverChange != null) {
+ // Allocate land cover change to grid cells
+ for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap: landCoverChange.getMap().entrySet()) {
+ LandCoverType fromLC = fromMap.getKey();
+ for (Entry<LandCoverType, Double> toMap: fromMap.getValue().entrySet()) {
+ LandCoverType toLC = toMap.getKey();
+ double change = toMap.getValue();
+ allocAllLandCropsTop(newLandUseRaster, keys, toLC, fromLC, change, locId);
+ }
}
}
@@ -319,15 +321,24 @@ public class GamsRasterOptimiser {
aggCFlux.setCarbonFlux(prevLC, newLC, aggregateMean(aggCFlux.getCarbonFlux(prevLC, newLC), suitableAreaSoFar,
carbonFlux, suitableAreaThisTime));
}
- // Wood yield
+ }
+ }
+
+ // Wood yield
+ for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
+ for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
double woodYield;
if (woodYieldItem != null) {
- woodYield = woodYieldItem.getWoodYield(prevLC, newLC);
+ if (woodYieldItem.checkForWoodKeys(prevLC, newLC)) {
+ woodYield = woodYieldItem.getWoodYield(prevLC, newLC);
+ } else {
+ continue;
+ }
} else {
woodYield = 0.0; // if missing data assume 0
}
- if (!aggWYield.checkForKeys(prevLC, newLC)) { // if not seen yet
+ if (!aggWYield.checkForWoodKeys(prevLC, newLC)) { // if not seen yet
aggWYield.setWoodYield(prevLC, newLC, woodYield);
} else {
aggWYield.setWoodYield(prevLC, newLC, aggregateMean(aggWYield.getWoodYield(prevLC, newLC), suitableAreaSoFar,
@@ -335,6 +346,29 @@ public class GamsRasterOptimiser {
}
}
}
+
+ // Timber yield
+ for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
+ for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
+ double timberYield;
+ if (woodYieldItem != null) {
+ if (woodYieldItem.checkForTimberKeys(prevLC, newLC)) {
+ timberYield = woodYieldItem.getTimberYield(prevLC, newLC);
+ } else {
+ continue;
+ }
+ } else {
+ timberYield = 0.0; // if missing data assume 0
+ }
+
+ if (!aggWYield.checkForTimberKeys(prevLC, newLC)) { // if not seen yet
+ aggWYield.setTimberYield(prevLC, newLC, timberYield);
+ } else {
+ aggWYield.setTimberYield(prevLC, newLC, aggregateMean(aggWYield.getTimberYield(prevLC, newLC), suitableAreaSoFar,
+ timberYield, suitableAreaThisTime));
+ }
+ }
+ }
// Crops yields and area fractions
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
index c6875491caa2ee6e1f779ccac611223ffa934eba..01dfd917b748226b3c8c22ea838072aa272ecd79 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
@@ -11,6 +11,7 @@ import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.DoubleMap;
+import ac.ed.lurg.utils.TripleMap;
import ac.sac.raster.RasterSet;
public class GamsRasterOutput {
@@ -20,6 +21,7 @@ public class GamsRasterOutput {
private Map<CropType, CropUsageData> cropUsageData;
private List<ForestRecord> newForest;
DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
+ private TripleMap<Integer, String, String, Double> naturalChanges;
double netCarbonFlux;
double timberProduction;
@@ -31,11 +33,12 @@ public class GamsRasterOutput {
public GamsRasterOutput(ModelStat status, RasterSet<LandUseItem> intensityRaster, Map<CropType, CropUsageData> cropUsageData,
List<ForestRecord> newForest, DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
- double netCarbonFlux, double timberProduction) {
+ TripleMap<Integer, String, String, Double> naturalChanges, double netCarbonFlux, double timberProduction) {
this(intensityRaster, cropUsageData);
this.status = status;
this.newForest = newForest;
this.reservedDeforested = reservedDeforested;
+ this.naturalChanges = naturalChanges;
this.netCarbonFlux = netCarbonFlux;
this.timberProduction = timberProduction;
}
@@ -60,6 +63,10 @@ public class GamsRasterOutput {
return reservedDeforested;
}
+ public TripleMap<Integer, String, String, Double> getNaturalChanges() {
+ return naturalChanges;
+ }
+
public double getNetCarbonFlux() {
return netCarbonFlux;
}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldItem.java b/src/ac/ed/lurg/landuse/WoodYieldItem.java
index 872143690a7ce535ab22b87cc4699c35c5a0dcc4..f259c1707ab582ee84e0f00eff1cd275eda84288 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldItem.java
+++ b/src/ac/ed/lurg/landuse/WoodYieldItem.java
@@ -38,13 +38,17 @@ public class WoodYieldItem implements RasterItem {
return timberYields.get(previousLandCover).get(newLandCover);
}
- public boolean checkForKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ public boolean checkForWoodKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
if (woodYields.containsKey(previousLandCover)) {
- if (woodYields.get(previousLandCover).containsKey(newLandCover)) {
- return true;
- } else {
- return false;
- }
+ return woodYields.get(previousLandCover).containsKey(newLandCover);
+ } else {
+ return false;
+ }
+ }
+
+ public boolean checkForTimberKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ if (timberYields.containsKey(previousLandCover)) {
+ return timberYields.get(previousLandCover).containsKey(newLandCover);
} else {
return false;
}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldReader.java b/src/ac/ed/lurg/landuse/WoodYieldReader.java
index a274adf236d9def7ace773e2409216961028af70..cfbaa2aa79faad52e2d4c256d90cfe64b45ce327 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldReader.java
+++ b/src/ac/ed/lurg/landuse/WoodYieldReader.java
@@ -9,7 +9,7 @@ import ac.sac.raster.RasterSet;
public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem> {
- private static final int MIN_COLS = 27;
+ private static final int MIN_COLS = 14;
private static final double conversionFactor = 10.0; // convert kgC/m2 to tC/ha
public WoodYieldReader(RasterSet<WoodYieldItem> woodYield) {
@@ -17,6 +17,21 @@ public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem>
}
protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "forC_to_crop") * conversionFactor);
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT") * conversionFactor);
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "forC_to_past") * conversionFactor);
+ item.setWoodYield(LandCoverType.NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "ntrl_to_crop") * conversionFactor);
+ item.setWoodYield(LandCoverType.NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "ntrl_to_forC") * conversionFactor);
+ item.setWoodYield(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT") * conversionFactor);
+ item.setWoodYield(LandCoverType.NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "ntrl_to_past") * conversionFactor);
+
+ item.setTimberYield(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "crop_to_forT_rot") * conversionFactor);
+ item.setTimberYield(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "past_to_forT_rot") * conversionFactor);
+ item.setTimberYield(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forT_to_forT_rot") * conversionFactor);
+ item.setTimberYield(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT_rot") * conversionFactor);
+ item.setTimberYield(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT_rot") * conversionFactor);
+
+ /*
item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.NATURAL, getValueForCol(rowValues, "crop_to_ntrl") * conversionFactor);
item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "crop_to_forT") * conversionFactor);
item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "crop_to_forC") * conversionFactor);
@@ -45,6 +60,7 @@ public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem>
item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.NATURAL, getValueForCol(rowValues, "forC_to_ntrl") * conversionFactor);
item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT") * conversionFactor);
item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "forC_to_crop") * conversionFactor);
- item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "forC_to_forC") * conversionFactor);
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "forC_to_forC") * conversionFactor);
+ */
}
}