diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 3fd99c4d56922feb671289b1d880827d58116fd7..4fa6d2381d6b4f5624865a90dd6845dad7d7fefd 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -142,7 +142,9 @@ $gdxin
pastureDecrease(location)
landCoverArea(land_cover, location) land cover area in Mha
landCoverChange(land_cover_before, land_cover_after, location) land cover change
- prematureDeforestationCost(location)
+ deforestedArea(forest, location) area of forest cut down
+ reservedAreaDeforested(forest, location) area of reserved (minimumLandCoverArea) forest cut down
+ prematureDeforestationCost
woodHarvest(location) total wood harvested
carbonFlux(location) total carbon flux
totalFeedDM
@@ -187,10 +189,12 @@ $gdxin
PASTURE_LAND_COVER_CALC(location) pasture area (land cover) equals pasture area (land use)
* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) constraint on land cover conversion
LAND_COVER_CHANGE_CALC(land_cover, location)
- LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) can't convert more land than was previously available
- LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) calculate no change area (e.g. cropland to cropland)
+ LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
+* LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) calculate no change area (e.g. cropland to cropland)
FOREST_EXPANSION_CONSTRAINT
- PREMATURE_DEFORESTATION_COST_CALC(location)
+ DEFORESTED_AREA_CALC(forest, location)
+ RESERVED_AREA_DEFORESTED_CALC(forest, location)
+ PREMATURE_DEFORESTATION_COST_CALC
WOOD_HARVEST_CALC(location) calc wood harvested
CARBON_FLUX_CALC(location) calc carbon flux
COST_EQ total cost objective function;
@@ -254,20 +258,29 @@ $gdxin
PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
+* landCoverChange(A, A, location) defined as unchanged land cover
LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= previousLandCoverArea(land_cover, location) +
- sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) - sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
+ sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) -
+ sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
- LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =L= previousLandCoverArea(land_cover, location);
+ LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= previousLandCoverArea(land_cover, location);
- LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =E= previousLandCoverArea(land_cover, location) -
- sum(land_cover_after$[not sameAs(land_cover, land_cover_after)], landCoverChange(land_cover, land_cover_after, location));
- FOREST_EXPANSION_CONSTRAINT .. sum((non_forest, forest, location), landCoverChange(non_forest, forest, location)) =L= sum(location, suitableLandArea(location)) * forestExpansionMaxRate;
+ FOREST_EXPANSION_CONSTRAINT .. sum((non_forest, forest, location), landCoverChange(non_forest, forest, location)) =L= sum((non_forest, location), previousLandCoverArea(non_forest, location)) * forestExpansionMaxRate;
* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =G= minimumLandCoverArea(land_cover, location);
- PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum(forest, (sum((land_cover_after)$[not sameAs(forest, land_cover_after)], landCoverChange(forest, land_cover_after, location)) -
- previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location))) * 1000;
+* PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum(forest, (sum((land_cover_after)$[not sameAs(forest, land_cover_after)], landCoverChange(forest, land_cover_after, location)) -
+* previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)));
+
+ DEFORESTED_AREA_CALC(forest, location) .. deforestedArea(forest, location) =E= sum((land_cover_after)$[not sameAs(forest, land_cover_after)], landCoverChange(forest, land_cover_after, location));
+
+* Equivalent to MAX(0, deforestedArea(forest, location) - minimumLandCoverArea(forest, location)) https://www.gams.com/latest/docs/UG_NLP_GoodFormulations.html#UG_NLP_GoodFormulations_ReformulatingDNLPModels
+ RESERVED_AREA_DEFORESTED_CALC(forest, location) .. reservedAreaDeforested(forest, location) =E= [(deforestedArea(forest, location) - [previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)] + delta) +
+ sqrt(sqr(deforestedArea(forest, location) - [previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)] + delta)) +
+ sqr(delta)] / 2;
+
+ PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, location), reservedAreaDeforested(forest, location) * 100);
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));
@@ -292,7 +305,8 @@ $gdxin
SUM(location, carbonFlux(location) * carbonPrice) - SUM(location, woodHarvest(location) * woodPrice) +
SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
- SUM(location, prematureDeforestationCost(location))
+
+ prematureDeforestationCost
);
MODEL LAND_USE /ALL/ ;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 0bf741dd8d49f38122b66723f647f51c8f81be49..ef6786e06d55f47001e901026750d7fb349d4a55 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -425,6 +425,6 @@ public class ModelConfig {
// Forestry parameters
public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 10.0); // $/tC-eq
public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 50.0); // $/tC-eq
- public static final int FOREST_LOCKIN_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 20); // cannot convert forest after planting for this many years
+ public static final int FOREST_LOCKIN_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 30); // cannot convert forest after planting for this many years
}
diff --git a/src/ac/ed/lurg/Timestep.java b/src/ac/ed/lurg/Timestep.java
index fe6b49ecb9fa126d1d59ba2e9e393c796752eec4..48b7c7b2e40b63c9bf24e53fcf8b438f8b505bcb 100644
--- a/src/ac/ed/lurg/Timestep.java
+++ b/src/ac/ed/lurg/Timestep.java
@@ -21,6 +21,14 @@ public class Timestep {
return ModelConfig.TIMESTEP_SIZE * (timestep - ModelConfig.START_TIMESTEP) + ModelConfig.BASE_YEAR;
}
+ public int getPreviousYear() {
+ if (timestep == ModelConfig.START_TIMESTEP)
+ throw new RuntimeException("Can't getPreviousYear for " + timestep);
+
+ return ModelConfig.TIMESTEP_SIZE * (getPreviousTimestep().getTimestep() - ModelConfig.START_TIMESTEP) + ModelConfig.BASE_YEAR;
+ }
+
+
public int getYieldYear() {
if (isInitialTimestep()) {
return ModelConfig.BASE_YEAR;
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 6c55218a66b3c0a6d2609957a404f2578bfc9666..421a98a332e418afbb07e18fe0910315f618f4bf 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -5,11 +5,14 @@ import java.io.IOException;
import java.nio.file.FileSystems;
import java.nio.file.Files;
import java.nio.file.StandardCopyOption;
+import java.util.ArrayList;
+import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
+import java.util.Set;
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.country.gams.GamsCountryInput;
@@ -225,11 +228,28 @@ public class CountryAgent extends AbstractCountryAgent {
}
}
}
+
+ // Logged timberForest area - to be converted to otherNatural
+ DoubleMap<Integer, LandCoverType, Double> loggedForest = new DoubleMap<Integer, LandCoverType, Double>(); //<Location, LandCoverType, Area>
+
+ if (currentTimestep.getTimestep() != ModelConfig.START_TIMESTEP) {
+ for (Entry<Integer, DoubleMap<Integer, LandCoverType, Double>> yearMap : minimumLandCover.getMap().entrySet()) {
+ if (yearMap.getKey() <= currentTimestep.getYear() && yearMap.getKey() > currentTimestep.getPreviousYear()) {
+ for (Entry<Integer, Map<LandCoverType, Double>> locMap : yearMap.getValue().getMap().entrySet()) {
+ for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
+ loggedForest.put(locMap.getKey(), coverMap.getKey(), coverMap.getValue());
+ }
+ }
+ }
+ }
+ }
+
+
GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints,
previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates);
GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData,
- carbonFluxData, woodYieldData, minimumLandCoverForTimestep, countryLevelInputs);
+ carbonFluxData, woodYieldData, minimumLandCoverForTimestep, loggedForest, countryLevelInputs);
return input;
}
@@ -349,5 +369,30 @@ public class CountryAgent extends AbstractCountryAgent {
minimumLandCover.put(expiryYear, locId, lcChange.getKey(), lcChange.getValue());
}
}
+ // Remove deforested area from minimum area requirement
+ DoubleMap<Integer, LandCoverType, Double> reservedDeforested = previousGamsRasterOutput.getReservedDeforested();
+ if (!reservedDeforested.getMap().isEmpty()) {
+ for (Entry<Integer, Map<LandCoverType, Double>> locMap : reservedDeforested.getMap().entrySet()) {
+ int locId = locMap.getKey();
+ for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
+ LandCoverType forestType = coverMap.getKey();
+ double deforestedArea = coverMap.getValue();
+ List<Integer> expiryYearsFound = new ArrayList<Integer>(minimumLandCover.getMap().keySet()) ;
+ Collections.sort(expiryYearsFound); // ascending order, remove area from most mature forest
+ while (deforestedArea > 0) {
+ for (int year : expiryYearsFound) {
+ if (year > currentTimestep.getYear()) {
+ double currentMinimumArea = minimumLandCover.get(year, locId, forestType);
+ if (currentMinimumArea <= 0)
+ LogWriter.printlnError("currentMinimumArea is <=0");
+ double areaToSubtract = Math.min(currentMinimumArea, deforestedArea);
+ minimumLandCover.put(year, locId, forestType, currentMinimumArea - areaToSubtract);
+ deforestedArea -= areaToSubtract;
+ }
+ }
+ }
+ }
+ }
+ }
}
}
\ 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 a7395b36a0018b1cc50abb63610d71cce5ccbf28..66eb021b7abb48b15ddbcb29c3803517ad7492ee 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -126,7 +126,7 @@ public class CountryAgentManager {
clusterIdRaster.putAll(ca.getYieldClusters());
} catch (Exception e) {
- LogWriter.printlnError("Exception processing " + ca.getCountry() + " will continue with other countries");
+ LogWriter.printlnError("Exception processing " + ca.getCountry() + " will continue with other countries");
LogWriter.print(e);
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 6844c309a2e982dfde3329a5413c80c2f8183e94..815eb0979f33dbf9b3111cb8e3c83bc558fb30bc 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -535,42 +535,7 @@ public class GamsLocationOptimiser {
LogWriter.println(String.format("\n%s %s: Total area= %.1f (crop=%.1f, pasture %.1f)",
inputData.getCountryInput().getCountry(), inputData.getTimestep().getYear(),
totalCropArea+totalPastureArea, totalCropArea, totalPastureArea));
-
-
-/*
- // Land cover areas (exc. cropland and pasture)
- GAMSVariable varLandCoverArea = outDB.getVariable("landCoverArea");
- double landCoverArea;
-
- for (GAMSVariableRecord rec : varLandCoverArea) {
- String lc = rec.getKeys()[0];
- String locationName = rec.getKeys()[1];
- int locId = Integer.parseInt(locationName);
-
- landCoverArea = rec.getLevel();
- LandUseItem landUseItem = landUses.lazyGet(locId);
-
- switch(lc) {
- case "forestT":
- landUseItem.setLandCoverArea(LandCoverType.TIMBER_FOREST, landCoverArea);
- break;
- case "forestC":
- landUseItem.setLandCoverArea(LandCoverType.CARBON_FOREST, landCoverArea);
- break;
- case "natural":
-
- double prevOtherNaturalArea = inputData.getPreviousLandUse().get(locId).getLandCoverArea(LandCoverType.OTHER_NATURAL);
- double prevUnmanagedForestArea = inputData.getPreviousLandUse().get(locId).getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
- double prevUnmanagedForestProp = (prevOtherNaturalArea + prevUnmanagedForestArea) != 0 ? prevUnmanagedForestArea / (prevUnmanagedForestArea + prevOtherNaturalArea) : 0;
-
- landUseItem.setLandCoverArea(LandCoverType.UNMANAGED_FOREST, landCoverArea * prevUnmanagedForestProp); // TODO Split into unmanaged forest and other natural
- landUseItem.setLandCoverArea(LandCoverType.OTHER_NATURAL, landCoverArea * (1 - prevUnmanagedForestProp)); // TODO Split into unmanaged forest and other natural
- break;
- }
-
- }
-*/
// Land cover change
GAMSVariable varLandCoverChange = outDB.getVariable("landCoverChange");
@@ -587,27 +552,6 @@ public class GamsLocationOptimiser {
if (!fromLC.equals(toLC)) { //Important as don't want to include unchanged land cover
landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.getForName(toLC), change);
}
- /*
- if (change == 0.0) {
- continue;
- }
- */
-
- /*
- if (fromLC.equals("natural") && toLC.equals("natural")) {
- continue;
- } else if (fromLC.equals("natural")) {
- double prevUnmanagedForestProp = getPrevUnmanagedForestProp(locId);
- landCoverChanges.put(locId, LandCoverType.UNMANAGED_FOREST, LandCoverType.getForName(toLC), change * prevUnmanagedForestProp);
- landCoverChanges.put(locId, LandCoverType.OTHER_NATURAL, LandCoverType.getForName(toLC), change * (1 - prevUnmanagedForestProp));
- } else if (toLC.equals("natural")) {
- double prevUnmanagedForestProp = getPrevUnmanagedForestProp(locId);
- landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.UNMANAGED_FOREST, change * prevUnmanagedForestProp);
- landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.OTHER_NATURAL, change * (1 - prevUnmanagedForestProp));
- } else {
- landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.getForName(toLC), change);
- }
- */
}
@@ -628,8 +572,24 @@ public class GamsLocationOptimiser {
*/
}
}
+
+ // reserved Forest which was deforested
+ GAMSVariable varReservedDeforested = outDB.getVariable("reservedAreaDeforested");
+ DoubleMap<Integer, LandCoverType, Double> reservedDeforested = new DoubleMap<Integer, LandCoverType, Double>();
+
+ for (GAMSVariableRecord rec : varReservedDeforested) {
+ String landCover = rec.getKeys()[0];
+ String locationName = rec.getKeys()[1];
+ int locId = Integer.parseInt(locationName);
+
+ double deforestedArea = rec.getLevel();
+
+ if (deforestedArea > 0.00001) {
+ reservedDeforested.put(locId, LandCoverType.getForName(landCover), deforestedArea);
+ }
+ }
- GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, minimumLandCoverAdditions);
+ GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, minimumLandCoverAdditions, reservedDeforested);
return results;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index a30c9d2deb5988673fabe9e7191318d7fec4f7e7..05a7a54f845b4a81d56c8ceec43a03bc290f80fb 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -18,18 +18,21 @@ public class GamsLocationOutput {
private Map<CropType, CropUsageData> cropUsageData;
TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges;
DoubleMap<Integer, LandCoverType, Double> minimumLandCover;
+ DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
public GamsLocationOutput(ModelStat status,
Map<Integer, LandUseItem> landUses,
Map<CropType, CropUsageData> cropUsageData,
TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChange,
- DoubleMap<Integer, LandCoverType, Double> minimumLandCover) {
+ DoubleMap<Integer, LandCoverType, Double> minimumLandCover,
+ DoubleMap<Integer, LandCoverType, Double> reservedDeforested) {
super();
this.status = status;
this.landUses = landUses;
this.cropUsageData = cropUsageData;
this.landCoverChanges = landCoverChange;
this.minimumLandCover = minimumLandCover;
+ this.reservedDeforested = reservedDeforested;
}
public ModelStat getStatus() {
@@ -50,4 +53,8 @@ public class GamsLocationOutput {
public DoubleMap<Integer, LandCoverType, Double> getMinimumLandCover() {
return minimumLandCover;
}
+
+ public DoubleMap<Integer, LandCoverType, Double> getReservedDeforested() {
+ return reservedDeforested;
+ }
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index 759107b3a3ddaf1836a4fedf684af1ac540310a8..e4e1171b2666853e2d71cba009cf5b7f1ce1661d 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -19,11 +19,12 @@ public class GamsRasterInput {
private RasterSet<CarbonFluxItem> carbonFluxes;
private RasterSet<WoodYieldItem> woodYields;
private DoubleMap<Integer, LandCoverType, Double> minimumLandCover;
+ private DoubleMap<Integer, LandCoverType, Double> loggedForest;
private GamsCountryInput countryInput;
public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses, RasterSet<IrrigationItem> irrigationCost,
RasterSet<CarbonFluxItem> carbonFluxes, RasterSet<WoodYieldItem> woodYields, DoubleMap<Integer, LandCoverType, Double> minimumLandCover,
- GamsCountryInput countryInput) {
+ DoubleMap<Integer, LandCoverType, Double> loggedForest, GamsCountryInput countryInput) {
super();
this.timestep = timestep;
this.yields = yields;
@@ -32,6 +33,7 @@ public class GamsRasterInput {
this.carbonFluxes = carbonFluxes;
this.woodYields = woodYields;
this.minimumLandCover = minimumLandCover;
+ this.loggedForest = loggedForest;
this.countryInput = countryInput;
}
@@ -59,6 +61,10 @@ public class GamsRasterInput {
return minimumLandCover;
}
+ public DoubleMap<Integer, LandCoverType, Double> getLoggedForest() {
+ return loggedForest;
+ }
+
public GamsCountryInput getCountryInput() {
return countryInput;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 74fe24d22d6fa318166433c231d91087c1cce6d5..b7dcecdafb4266639cf34ecfa994d8ccf40e1e0d 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -53,7 +53,7 @@ public class GamsRasterOptimiser {
private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) {
RasterSet<LandUseItem> newIntensityRaster = allocAreas(gamsInput.getPreviousLandUse(), gamsOutput, gamsInput.getTimestep().getYear());
- return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getMinimumLandCover());
+ return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getMinimumLandCover(), gamsOutput.getReservedDeforested());
}
private RasterSet<LandUseItem> createWithSameLandCovers(RasterSet<LandUseItem> toCopy) {
@@ -98,19 +98,9 @@ public class GamsRasterOptimiser {
private RasterSet<LandUseItem> allocAreas(Map<Integer, ? extends LandUseItem> prevAreasAgg, GamsLocationOutput gamsOutput, int year) {
RasterSet<LandUseItem> newLandUseRaster = createWithSameLandCovers(rasterInputData.getPreviousLandUses());
- //checkedTotalAreas(rasterInputData.getPreviousAreas(), "old");
- //checkedTotalAreas(newAreaRaster, "new");
-
for (Map.Entry<Integer, LandUseItem> entry : gamsOutput.getLandUses().entrySet()) {
Integer locId = entry.getKey();
LandUseItem newLandUseAggItem = entry.getValue();
-// LandUseItem prevLandUseAggItem = prevAreasAgg.get(locId);
-
- /*
- for (Entry<RasterKey, IntegerRasterItem> foo : mapping.entrySet()) {
- LogWriter.println(foo.getKey() + "|" + foo.getValue().getInt());
- }
- */
Set<RasterKey> keys = new HashSet<RasterKey>();
for (Entry<RasterKey, IntegerRasterItem> mapEntry : mapping.entrySet()) {
@@ -120,104 +110,10 @@ public class GamsRasterOptimiser {
keys.add(mapEntry.getKey());
}
-
RasterSet<LandUseItem> landUseItemsForLocation = newLandUseRaster.createSubsetForKeys(keys);
- if (rasterInputData.getCountryInput().getCountry().getName().equals("Peru & Ecuador")) {
- LogWriter.println("foo");
- }
-
DoubleMap<LandCoverType, LandCoverType, Double> landCoverChange = gamsOutput.getLandCoverChanges().getInnerMap(locId);
- /*
- if (landCoverChange == null) {
- landCoverChange = new DoubleMap<LandCoverType, LandCoverType, Double>();
- }
- */
-
- double prevTimberForest = 0, prevCarbonForest = 0, prevUnmanagedForest = 0, prevNatural = 0;
- double protectedAreasPastureChange = 0, protectedAreasCroplandChange = 0;
-
- for (LandUseItem luItem: landUseItemsForLocation.values()) {
- prevTimberForest += luItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST);
- prevCarbonForest += luItem.getUnprotectedLandCoverArea(LandCoverType.CARBON_FOREST);
- prevUnmanagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST);
- prevNatural += luItem.getTotalUnprotectedNatural();
- }
-
- double prevForestTotal = prevTimberForest + + prevCarbonForest + prevUnmanagedForest;
- double prevForestToNaturalFraction = (prevNatural > 0) ? prevForestTotal / prevNatural : 0;
- double prevForestTimberFraction = (prevForestTotal > 0) ? prevTimberForest / prevForestTotal : 0;
- double prevForestCarbonFraction = (prevForestTotal > 0) ? prevCarbonForest / prevForestTotal : 0;
- double prevForestUnmanagedFraction = (prevForestTotal > 0) ? prevUnmanagedForest / prevForestTotal : 0;
-
- // Floating point errors
- if(prevForestToNaturalFraction > 1.0) prevForestToNaturalFraction = 1.0;
- if(prevForestTimberFraction > 1.0) prevForestTimberFraction = 1.0;
- if(prevForestCarbonFraction > 1.0) prevForestCarbonFraction = 1.0;
- if(prevForestUnmanagedFraction > 1.0) prevForestUnmanagedFraction = 1.0;
-
-
- //this tries to honour spatial distribution of protected areas, force removal of pasture and cropland from protected areas
- if (ModelConfig.FORCE_PROTECTED_AREAS && year >= ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR) {
- double protectedLandShortfall = 0;
-
- for (LandUseItem luItem: landUseItemsForLocation.values()) {
-
- double suitableLand = luItem.getSuitableArea();
- double prevCropland = luItem.getLandCoverArea(LandCoverType.CROPLAND);
- double prevPasture = luItem.getLandCoverArea(LandCoverType.PASTURE);
-
- //if we have more agricultural land in this cell than we should, remove it and add it to land needing reallocated
- if (suitableLand < (prevCropland + prevPasture)) {
-
- double excessAgriculture = prevCropland + prevPasture - suitableLand;
- double prevPastureFract = prevPasture / (prevCropland + prevPasture);
- double pastureShortfall = excessAgriculture * prevPastureFract;
- double croplandShortfall = excessAgriculture * (1 - prevPastureFract);
-
- protectedLandShortfall += luItem.moveAreas(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE,
- pastureShortfall * prevForestToNaturalFraction * prevForestTimberFraction);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE,
- pastureShortfall * prevForestToNaturalFraction * prevForestCarbonFraction);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.UNMANAGED_FOREST, LandCoverType.PASTURE,
- pastureShortfall * prevForestToNaturalFraction * prevForestUnmanagedFraction);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.OTHER_NATURAL, LandCoverType.PASTURE,
- pastureShortfall * (1 - prevForestToNaturalFraction));
- protectedLandShortfall += luItem.moveAreas(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND,
- croplandShortfall * prevForestToNaturalFraction * prevForestTimberFraction);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND,
- croplandShortfall * prevForestToNaturalFraction * prevForestCarbonFraction);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.UNMANAGED_FOREST, LandCoverType.CROPLAND,
- croplandShortfall * prevForestToNaturalFraction * prevForestUnmanagedFraction);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.OTHER_NATURAL, LandCoverType.CROPLAND,
- croplandShortfall * (1 - prevForestToNaturalFraction));
-
- protectedAreasPastureChange+= pastureShortfall;
- protectedAreasCroplandChange+= croplandShortfall;
- }
-
- }
-
- if(protectedLandShortfall > 0.00001) {
- LogWriter.printlnWarning("locID " + locId + " Not able to incorporate all protected area changes " + protectedLandShortfall + " remains unprotected.");
- }
-
- if (DEBUG) {
- LogWriter.println("protectedAreasPastureChange " + protectedAreasPastureChange + " protectedAreasCroplandChange " + protectedAreasCroplandChange);
- }
-
- double totalNaturalArea = newLandUseAggItem.getTotalNatural();
-
- // Add protected area changes to GAMS changes
- for (LandCoverType lc : LandCoverType.getProtectibleTypes()) {
- double fraction = newLandUseAggItem.getLandCoverArea(lc) / totalNaturalArea;
- landCoverChange.addTo(lc, LandCoverType.CROPLAND, protectedAreasCroplandChange * fraction);
- landCoverChange.addTo(lc, LandCoverType.PASTURE, protectedAreasPastureChange * fraction);
- }
-
- }
-
// Allocate land cover change to grid cells
for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap: landCoverChange.getMap().entrySet()) {
LandCoverType fromLC = fromMap.getKey();
@@ -235,7 +131,6 @@ public class GamsRasterOptimiser {
newLandUseItem.setCropFraction(crop, newLandUseAggItem.getCropFraction(crop));
newLandUseItem.setIntensity(crop, newLandUseAggItem.getIntensity(crop)); // intensities constant over single aggregated land category
}
-
}
}
@@ -278,12 +173,13 @@ public class GamsRasterOptimiser {
}
for (RasterKey key : keys) {
- //if (DEBUG) LogWriter.println(" Processing raster key " + key);
+ //if (DEBUG) LogWriter.println("Processing raster key " + key);
LandUseItem newLandUseItem = newLandUseRaster.get(key);
if (newLandUseItem!=null) {
- double cellChange = (fromType.isProtectable() && totalUnprotectedLC > 0) ? change * newLandUseItem.getUnprotectedLandCoverArea(fromType)/totalUnprotectedLC : change / keys.size();
+ double cellChange = (totalUnprotectedLC > 0) ? change * newLandUseItem.getUnprotectedLandCoverArea(fromType)/totalUnprotectedLC : change / keys.size();
double shortfall = newLandUseItem.moveAreas(toType, fromType, cellChange);
+ shortfall += newLandUseItem.moveUnprotectedArea(toType, fromType, cellChange);
if (shortfall == 0)
keysWithSpace.add(key);
else
@@ -293,12 +189,6 @@ public class GamsRasterOptimiser {
}
}
- // Now we check for unprotected areas, I think this is redundant. However, it does cause a problem so leaving just in case it is needed.
- if (totalShortfall > 0 & keysWithSpace.size() > 0) { // more to allocate and some free areas to allocate into
- if (DEBUG) LogWriter.println("A total land cover of " + totalShortfall + " is unallocated, so trying again with areas with space");
- totalShortfall = allocAllLandCrops(newLandUseRaster, keysWithSpace, toType, fromType, totalShortfall);
- }
-
if (DEBUG) LogWriter.println("allocAllLandCrops: from " + fromType + " to " + toType + " change " + change + " totalShortfall " + totalShortfall);
return totalShortfall;
}
@@ -309,43 +199,11 @@ public class GamsRasterOptimiser {
YieldRaster yieldRaster = rasterInputData.getYields();
RasterSet<LandUseItem> landUseRaster = rasterInputData.getPreviousLandUses();
- // LandUseItem test = landUseRaster.getFromCoordinates(71.25, 35.25);
-
- /* new RasterOutputer<AreasItem>(cropAreaRaster, "BeforeLandCover") {
- @Override
- public Double getValue(RasterKey location) {
- AreasItem area = results.get(location);
- if (area == null)
- return null;
-
- int count = 0;
- double maxLC = 0;
- int maxId = 0;
- for (LandCoverType l : LandCoverType.values()) {
- count++;
- if (maxLC < area.getLandCoverArea(l)) {
- maxLC = area.getLandCoverArea(l);
- maxId = count;
- }
- }
- return (double)maxId;
- }
-
- @Override
- public int convertToPixelValue(double value) {
- return (int) (1 + value * 2);
- }
- }.writeOutput(true); */
-
-
RasterSet<IrrigationItem> irrigRaster = rasterInputData.getIrrigationData();
RasterSet<CarbonFluxItem> carbonFluxRaster = rasterInputData.getCarbonFluxes();
RasterSet<WoodYieldItem> woodYieldRaster = rasterInputData.getWoodYields();
-
- // YieldResponsesItem yresp = yieldRaster.get(new RasterKey(40,38));//getFromCoordinates(-118.0, -33.0);
- // LogWriter.printlnError("Test key2: " + yresp.getYieldMax(CropType.WHEAT) + ", " + yresp.getYieldFertOnly(CropType.WHEAT) + ", " + yresp.getYieldIrrigOnly(CropType.WHEAT));
// look for inconsistencies
for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
@@ -359,6 +217,7 @@ public class GamsRasterOptimiser {
}
}
+ // Aggregation
LazyTreeMap<Integer, YieldResponsesItem> aggregatedYields = new LazyTreeMap<Integer, YieldResponsesItem>() {
protected YieldResponsesItem createValue() { return new YieldResponsesItem(); }
};
@@ -375,7 +234,6 @@ public class GamsRasterOptimiser {
protected WoodYieldItem createValue() { return new WoodYieldItem(); }
};
-
int countFound = 0, countMissing = 0;
for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
@@ -410,7 +268,14 @@ public class GamsRasterOptimiser {
CarbonFluxItem aggCFlux = aggregatedCarbonFluxes.lazyGet(clusterId);
WoodYieldItem aggWYield = aggregatedWoodYields.lazyGet(clusterId);
- landUseItem.updateSuitableArea(rasterInputData.getTimestep().getYear());
+ /*
+ if (clusterId == 4 && rasterInputData.getTimestep().getYear() == 2011) {
+ int foo = 1;
+ }
+ */
+
+ // TODO disabled protected area update
+ //landUseItem.updateSuitableArea(rasterInputData.getTimestep().getYear());
double suitableAreaThisTime = landUseItem.getSuitableArea();
double suitableAreaSoFar = aggLandUse.getSuitableArea();
@@ -514,47 +379,33 @@ public class GamsRasterOptimiser {
}
}
-/*
- for (Entry<RasterKey, LandUseItem> entry : rasterInputData.getPreviousLandUses().entrySet()) {
- LandUseItem landUseItem = entry.getValue();
- RasterKey key = entry.getKey();
- int clusterId = mapping.get(key).getInt();
- LandUseItem aggLandUse = aggregatedAreas.lazyGet(clusterId);
-
- // Land covers areas
- for (LandCoverType landType : LandCoverType.values()) {
- double areaThisTime = landUseItem.getUnprotectedLandCoverArea(landType);
- double areaSoFar = aggLandUse.getUnprotectedLandCoverArea(landType);
-
- aggLandUse.setLandCoverArea(landType, areaSoFar + areaThisTime);
- }
- }
-*/
-/*
- if (rasterInputData.getCountryInput().getCountry().getName().equals("Central America")) {
- LogWriter.println("foo");
- }
-*/
+
LogWriter.println("YieldResponsesItem: " + rasterInputData.getCountryInput().getCountry() + ", countFound=" + countFound + ", countMissing=" + countMissing);
-
-
- // for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) {
- // LogWriter.println(e.getKey() + " zone has " + e.getValue().size() + " raster areas");
-
- /* CropType[] cs = {CropType.WHEAT, CropType.MAIZE};
- for (CropType c : cs) {
- for (RasterKey key : e.getValue()) {
- YieldResponsesItem yri = yieldRaster.get(key);
- LogWriter.println(String.format("\t%s: x=%.1f y=%.1f: %s yields (none,fmid,fmax,irrig,both) %.1f, %.1f, %.1f, %.1f, %.1f",
- key, yieldRaster.getXCoordin(key), yieldRaster.getYCoordin(key), c.getGamsName(),
- yri.getYieldNone(c), yri.getYieldMidFertOnly(c), yri.getYieldFertOnly(c), yri.getYieldIrrigOnly(c), yri.getYieldMax(c)));
- }
- LogWriter.println("");
- } */
- // }
checkedTotalAreas(landUseRaster, "before");
checkedTotalAreas(aggregatedAreas, "after");
+
+ // Auto convert timberForest to otherNatural once matured and cut down
+ if (!rasterInputData.getLoggedForest().getMap().isEmpty()) {
+ for (Entry<Integer, Map<LandCoverType, Double>> locMap : rasterInputData.getLoggedForest().getMap().entrySet()) {
+ int locId = locMap.getKey();
+
+ Set<RasterKey> keys = new HashSet<RasterKey>();
+ for (Entry<RasterKey, IntegerRasterItem> mapEntry : mapping.entrySet()) {
+ IntegerRasterItem iri = mapEntry.getValue();
+ if (iri != null && locId == iri.getInt())
+ keys.add(mapEntry.getKey());
+ }
+
+ for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
+ // move land for GAMS input
+ aggregatedAreas.get(locId).moveUnprotectedArea(LandCoverType.OTHER_NATURAL, coverMap.getKey(), coverMap.getValue());
+ // move land for raster
+ allocAllLandCropsTop(rasterInputData.getPreviousLandUses(), keys, LandCoverType.OTHER_NATURAL, coverMap.getKey(), coverMap.getValue(), locId);
+ }
+ }
+ }
+
return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts,
aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getMinimumLandCover(), rasterInputData.getCountryInput());
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
index 78bbbfefd903a86d6dee4ce44e3c4f69d496166e..b141c061d57eb6cbb8c882b5a7e44ff8114843e1 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
@@ -17,6 +17,7 @@ public class GamsRasterOutput {
private RasterSet<LandUseItem> landUses;
private Map<CropType, CropUsageData> cropUsageData;
private DoubleMap<Integer, LandCoverType, Double> minimumLandCoverAdditions;
+ DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
public GamsRasterOutput(RasterSet<LandUseItem> landUses, Map<CropType, CropUsageData> cropUsageData) {
super();
@@ -25,10 +26,11 @@ public class GamsRasterOutput {
}
public GamsRasterOutput(ModelStat status, RasterSet<LandUseItem> intensityRaster, Map<CropType, CropUsageData> cropUsageData,
- DoubleMap<Integer, LandCoverType, Double> minimumLandCoverAdditions) {
+ DoubleMap<Integer, LandCoverType, Double> minimumLandCoverAdditions, DoubleMap<Integer, LandCoverType, Double> reservedDeforested) {
this(intensityRaster, cropUsageData);
this.status = status;
this.minimumLandCoverAdditions = minimumLandCoverAdditions;
+ this.reservedDeforested = reservedDeforested;
}
public ModelStat getStatus() {
@@ -46,4 +48,8 @@ public class GamsRasterOutput {
public DoubleMap<Integer, LandCoverType, Double> getMinimumLandCoverAdditions() {
return minimumLandCoverAdditions;
}
+
+ public DoubleMap<Integer, LandCoverType, Double> getReservedDeforested() {
+ return reservedDeforested;
+ }
}
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index 07541d8f566be62968d43de8518629314fb25c8b..1e6c5cb692c667277810c0a3c40bc8dbeb0d3d0f 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -9,6 +9,7 @@ import ac.ed.lurg.types.CropToDouble;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.Interpolator;
+import ac.ed.lurg.utils.LogWriter;
import ac.sac.raster.InterpolatingRasterItem;
public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serializable {
@@ -80,7 +81,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
public void setUnavailableArea(double unavailF) {
- unavailableArea = unavailF*getTotalLandCoverArea();
+ // subtract BARREN from unavailableArea, set to 0 if negative
+ unavailableArea = Math.max(unavailF*getTotalLandCoverArea() - getLandCoverArea(LandCoverType.BARREN), 0.0);
}
public double getUnavailableArea() {
@@ -194,15 +196,39 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
double prevTo = getLandCoverArea(toType);
double prevFrom = getLandCoverArea(fromType);
- double availLC = fromType.isProtectable() ? Math.min(prevFrom, unprotectedAreas.get(fromType)) : prevFrom; // if protected can not go past protected area
- double actuallyChanged = Math.min(availLC, changeReq);
+ double actuallyChanged = Math.min(prevFrom, changeReq);
+
+ /*
+ if (changeReq - actuallyChanged > 0) {
+ throw new RuntimeException("Move area shortfall. Should not happen due to GAMS constraint.");
+ }*/
setLandCoverArea(toType, prevTo + actuallyChanged);
setLandCoverArea(fromType, prevFrom - actuallyChanged);
-
+
return changeReq - actuallyChanged;
}
+
+ public double moveUnprotectedArea(LandCoverType toType, LandCoverType fromType, double area) {
+
+ double prevTo = getUnprotectedLandCoverArea(toType);
+ double prevFrom = getUnprotectedLandCoverArea(fromType);
+ double areaMoved;
+
+ areaMoved = Math.min(area, prevFrom);
+
+ unprotectedAreas.put(toType, prevTo + areaMoved);
+ unprotectedAreas.put(fromType, prevFrom - areaMoved);
+
+ /*
+ if (area - areaMoved > 0) {
+ throw new RuntimeException("Move unprotected area shortfall. Should not happen due to GAMS constraint.");
+ }
+ */
+
+ return area - areaMoved;
+ }
public double getCropArea(CropType c) {
Double d = getLandCoverArea(LandCoverType.CROPLAND) * getCropFraction(c);
@@ -279,7 +305,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
}
public double getTotalNatural() {
- double totalNatural = getLandCoverArea(LandCoverType.OTHER_NATURAL) + getLandCoverArea(LandCoverType.TIMBER_FOREST) + getLandCoverArea(LandCoverType.CARBON_FOREST) + getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
+ double totalNatural = getLandCoverArea(LandCoverType.OTHER_NATURAL) + getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
return totalNatural;
}
@@ -304,18 +330,25 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
}
private double getProtectedandUnavailable() {
- double totalSuitableLandCover = getTotalLandCoverArea()-getLandCoverArea(LandCoverType.BARREN)-getLandCoverArea(LandCoverType.URBAN);
- if (totalSuitableLandCover <=0)
- return 0;
- double minAndProtectedA = getUnavailableArea() + getProtectedArea() + totalSuitableLandCover * ModelConfig.MIN_NATURAL_RATE;
- return Math.min(totalSuitableLandCover, minAndProtectedA);
+ double totalSuitableLandCover = getTotalLandCoverArea() - getLandCoverArea(LandCoverType.BARREN) - getLandCoverArea(LandCoverType.URBAN);
+ if (totalSuitableLandCover < 0)
+ throw new RuntimeException("totalSuitableLandCover < 0");
+ //return 0;
+ // Protected area should be at least MIN_NATURAL_RATE
+ double desiredProtectedArea = Math.max(getProtectedArea(), totalSuitableLandCover * ModelConfig.MIN_NATURAL_RATE);
+ // if desired protected is more than suitable then protect all suitable
+ if (desiredProtectedArea + getUnavailableArea() > totalSuitableLandCover) {
+ LogWriter.printlnWarning("Desired protected area exceeds total suitable area");
+ }
+ return Math.min(totalSuitableLandCover, desiredProtectedArea + getUnavailableArea());
}
/** averages protected areas across land cover types */
private void updateUnprotectedAreas() {
double desiredProtected = getProtectedandUnavailable();
double totalNatural = getTotalNatural();
- double unprotectedNat = Math.max(0, totalNatural - desiredProtected); // if we are already using more than is protected then the unprotectedArea is 0
+ double unprotectedNat = Math.max(0, totalNatural - desiredProtected); // if desiredProtected more than total Natural then protect all
+
double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0;
for (LandCoverType landType : LandCoverType.values()) {