diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 5b1afba8e6e4934fdf665400f145627b629e0dce..ec4b52772d7201aa8ae46d395d4a1d04af5c8526 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -36,7 +36,7 @@
PARAMETER importPrices(import_crop) prices for imports;
PARAMETER maxNetImport(import_crop) maximum net import for each crop based on world market;
PARAMETER minNetImport(import_crop) minimum net import for each crop based on world market;
- PARAMETER irrigCost(location) irrigation cost in cost per 1000 Mlitre or Mha for each litre per m2;
+ PARAMETER irrigCost(location) irrigation cost in cost per 1000 Mlitre or Mha for each litre per m2;
PARAMETER irrigMaxRate(crop, location) max water application rate irrigation in litre per m2;
PARAMETER irrigConstraint(location) max water available for irrigation in litre per m2;
PARAMETER minDemandPerCereal(cereal_crop) min demand for each cereal crop as factor of all cereals;
@@ -44,11 +44,11 @@
PARAMETER seedAndWasteRate(all_types) rate of use for seed and waste combined;
PARAMETER subsidyRate(crop) rates of subsidy compared to costs;
- PARAMETER agriExpansionCost(location) price for increase agricultural area varies based on amount of managed or unmanaged forest;
- SCALAR pastureIncCost price for increasing pasture area;
- SCALAR cropIncCost price for increasing crop;
- SCALAR cropDecCost price for decreasing crop;
- SCALAR pastureDecCost price for decreasing pasture;
+ PARAMETER agriExpansionCost(location) price for increase agricultural area varies based on amount of managed or unmanaged forest - 1000$ per ha;
+ SCALAR pastureIncCost price for increasing pasture area - 1000$ per ha;
+ SCALAR cropIncCost price for increasing crop - 1000$ per ha;
+ SCALAR cropDecCost price for decreasing crop - 1000$ per ha;
+ SCALAR pastureDecCost price for decreasing pasture - 1000$ per ha;
SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
SCALAR fertiliserUnitCost fert cost at max fert rate;
@@ -61,10 +61,10 @@
PARAMETER previousLandCoverArea(land_cover, location) land cover area in Mha;
* PARAMETER minimumLandCoverArea(land_cover, location) minimum land cover area to constrain conversion;
- PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location) carbon flux;
- PARAMETER woodYield(land_cover_before, land_cover_after, location) wood yield;
- SCALAR carbonPrice price of carbon;
- SCALAR woodPrice price of wood and timber;
+ 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;
+ SCALAR carbonPrice price of carbon - $ per tonne or million$ per Mt;
+ SCALAR woodPrice price of wood and timber - $ per tonne or million$ per Mt;
*$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
$gdxin %gdxincname%
@@ -77,9 +77,22 @@ $load previousLandCoverArea, carbonFluxRate, woodYield, carbonPrice, woodPrice
*$load minimumLandCoverArea
$gdxin
-* need to convert from $/kg to $/Mt
- carbonPrice = carbonPrice * 0.01;
- woodPrice = woodPrice * 0.01;
+* convert units from $/t to billion$/Mt
+ carbonPrice = carbonPrice * 0.001;
+ woodPrice = woodPrice * 0.001;
+
+* convert units from 1000$/ha to million$/Mha
+* agriExpansionCost(location) = agriExpansionCost(location) * 1000;
+* pastureIncCost = pastureIncCost * 1000;
+* cropIncCost = cropIncCost * 1000;
+* cropDecCost = cropDecCost * 1000;
+* pastureDecCost = pastureDecCost * 1000;
+
+* covert units from 1000$/t to million$/Mt
+* exportPrices(import_crop) = exportPrices(import_crop) * 1000;
+* importPrices(import_crop) = importPrices(import_crop) * 1000;
+
+
SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /;
@@ -149,8 +162,8 @@ $gdxin
landCoverArea, landCoverChange, woodHarvest;
EQUATIONS
- UNIT_COST_EQ(crop, location) cost per area
- YIELD_EQ(crop, location) yield given chosen intensity
+ UNIT_COST_EQ(crop, location) cost per area - million$ per Mha
+ YIELD_EQ(crop, location) yield given chosen intensity - tonnes per hectare
CROP_DEMAND_CONSTRAINT(crop) satisfy demand for individual crops
TOTAL_CEREAL_DEMAND_CONSTRAINT satisfy demand for combined cereals
TOTAL_OIL_PULSE_DEMAND_CONSTRAINT satisfy demand for combined oilcrops and pulses
@@ -178,7 +191,6 @@ $gdxin
CROPLAND_LAND_COVER_CALC(location) cropland area equals sum of all crop areas
PASTURE_LAND_COVER_CALC(location) pasture area (land cover) equals pasture area (land use)
* MINIMUM_LAND_COVER_CONSTRAINT(location) constraint on land cover conversion
-* TOTAL_LAND_COVER_CONSTRAINT(location) total land cover equals suitable area
LAND_COVER_CHANGE_CALC(land_cover, location)
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location)
LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location)
@@ -190,7 +202,7 @@ $gdxin
fertiliserUnitCost * fertI(crop, location) +
irrigCost(location) * irrigMaxRate(crop, location) * irrigI(crop, location) +
otherIntCost(crop) * otherIntensity(crop, location)
- ) ;
+ );
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yieldNone(crop, location) +
@@ -251,8 +263,6 @@ $gdxin
LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =E= 0;
-* TOTAL_LAND_COVER_CONSTRAINT(location) .. sum(land_cover, landCoverArea(land_cover, location)) =E= suitableLandArea(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));
CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
diff --git a/debug_config.properties b/debug_config.properties
index a59915593399115739bbf169389cc39e26bfb4eb..c161ff3a4daf04abe5e0566aa3d4a3d142b82c6b 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -1,24 +1,22 @@
-BASE_DIR=..
+BASE_DIR=C:/Users/Bart/git/plumv2
-YIELD_DIR=/Users/palexand/Documents/LURG/LPJ/LPJGPLUM_remap6p7_20190225/rcp26
+YIELD_DIR=C:/Users/Bart/Documents/PhD/LURG/rcp60
-END_TIMESTEP=0
+OUTPUT_DIR = C:/Users/Bart/git/plumv2/output
+
+END_TIMESTEP=10
END_FIRST_STAGE_CALIBRATION=10
TIMESTEP_SIZE=1
IS_CALIBRATION_RUN=true
-DEBUG_JUST_DEMAND_OUTPUT=true
-GENERATE_NEW_YIELD_CLUSTERS=false
+GENERATE_NEW_YIELD_CLUSTERS=true
-DEMAND_PRICE_IMPORT_AND_PROD_COST=true
+DEMAND_PRICE_IMPORT_AND_PROD_COST=false
-EXTRAPOLATE_YIELD_FERT_RESPONSE=false
+YIELD_FILENAME=yield.out
DEBUG_LIMIT_COUNTRIES=false
-DEBUG_COUNTRY_NAME=Italy
-
-ADJUST_DIET_PREFS=true
-
-TARGET_DIET_FILE=/Users/palexand/Downloads/Tau_diet.csv
\ No newline at end of file
+DEBUG_COUNTRY_NAME=Peru & Ecuador
+GAMS_COUNTRY_TO_SAVE=Peru & Ecuador
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index d8f42418536952d2ca5aa3093940e213cb134257..c33525081df0405ab433b74521c38ecda6201833 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -423,7 +423,7 @@ public class ModelConfig {
public static final String TARGET_DIET_FILE = getProperty("TARGET_DIET_FILE", DATA_DIR + File.separator + "TargetDiet.txt");
// Forestry parameters
- public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 20.0); // $/tC
- public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.175); // $/tC
+ 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
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 9f1d7b93f060c0ddef00a0288a924e1a1d5f27f6..26eaa91fa082a1d09e1c95647206adbd35cbd0fa 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -38,6 +38,7 @@ import ac.ed.lurg.types.Parameter;
import ac.ed.lurg.types.YieldType;
import ac.ed.lurg.utils.LazyHashMap;
import ac.ed.lurg.utils.LogWriter;
+import ac.ed.lurg.utils.TripleMap;
import ac.ed.lurg.yield.YieldResponsesItem;
public class GamsLocationOptimiser {
@@ -74,6 +75,9 @@ public class GamsLocationOptimiser {
long startTime = System.currentTimeMillis();
gamsJob.run(opt, inDB);
+ if (inputData.getCountryInput().getCountry().getName().equals("Peru & Ecuador"))
+ LogWriter.println("" + inputData.getPreviousLandUse().get(2).getUnprotectedLandCoverArea(LandCoverType.CROPLAND));
+
if (ModelConfig.CLEANUP_GAMS_DIR)
cleanup(ws.workingDirectory());
@@ -560,7 +564,7 @@ public class GamsLocationOptimiser {
// Land cover change
GAMSVariable varLandCoverChange = outDB.getVariable("landCoverChange");
- LandCoverChangeItem landCoverChange = new LandCoverChangeItem();
+ TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges = new TripleMap<Integer, LandCoverType, LandCoverType, Double>();
for (GAMSVariableRecord rec : varLandCoverChange) {
@@ -577,19 +581,20 @@ public class GamsLocationOptimiser {
if (fromLC.equals("natural")) {
double prevUnmanagedForestProp = getPrevUnmanagedForestProp(locId);
- landCoverChange.setLandCoverChange(LandCoverType.UNMANAGED_FOREST, LandCoverType.getForName(toLC), locId, change * prevUnmanagedForestProp);
- landCoverChange.setLandCoverChange(LandCoverType.OTHER_NATURAL, LandCoverType.getForName(toLC), locId, change * (1 - prevUnmanagedForestProp));
+ 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);
- landCoverChange.setLandCoverChange(LandCoverType.getForName(fromLC), LandCoverType.UNMANAGED_FOREST, locId, change * prevUnmanagedForestProp);
- landCoverChange.setLandCoverChange(LandCoverType.getForName(fromLC), LandCoverType.OTHER_NATURAL, locId, change * (1 - prevUnmanagedForestProp));
+ landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.UNMANAGED_FOREST, change * prevUnmanagedForestProp);
+ landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.OTHER_NATURAL, change * (1 - prevUnmanagedForestProp));
} else {
- landCoverChange.setLandCoverChange(LandCoverType.getForName(fromLC), LandCoverType.getForName(toLC), locId, change);
+ landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.getForName(toLC), change);
}
+ // Potential case where fromLC = natural and toLC = natural unaccounted for but should never happen due to GAMS constraint (LAND_COVER_SELF_CHANGE_CONTRAINT).
}
- GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChange);
+ GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges);
return results;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index b64a6d6ddb97826f6f64313cd44e042fe22a8f8b..258064949b97f96c58357e0115a45fdea64f659d 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -5,26 +5,29 @@ import java.util.Map;
import com.gams.api.GAMSGlobals.ModelStat;
import ac.ed.lurg.landuse.CropUsageData;
-import ac.ed.lurg.landuse.LandCoverChangeItem;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.TripleMap;
+import ac.ed.lurg.types.LandCoverType;
public class GamsLocationOutput {
ModelStat status;
Map<Integer, LandUseItem> landUses; // data mapped from id (not raster)
private Map<CropType, CropUsageData> cropUsageData;
- LandCoverChangeItem landCoverChange = new LandCoverChangeItem();
+ //Map<Integer, DoubleMap<LandCoverType, LandCoverType, Double>> landCoverChange;
+ //LandCoverChangeItem landCoverChange;
+ TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges;
public GamsLocationOutput(ModelStat status,
Map<Integer, LandUseItem> landUses,
Map<CropType, CropUsageData> cropUsageData,
- LandCoverChangeItem landCoverChange) {
+ TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChange) {
super();
this.status = status;
this.landUses = landUses;
this.cropUsageData = cropUsageData;
- this.landCoverChange = landCoverChange;
+ this.landCoverChanges = landCoverChange;
}
public ModelStat getStatus() {
@@ -38,7 +41,7 @@ public class GamsLocationOutput {
return cropUsageData;
}
- public LandCoverChangeItem getLandCoverChange() {
- return landCoverChange;
+ public TripleMap<Integer, LandCoverType, LandCoverType, Double> getLandCoverChanges() {
+ return landCoverChanges;
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index d4b6b9394d6222385dbc8c83945017bfb9454968..a520bf68552b8a25b8078899ab796e4d94964a49 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -5,6 +5,7 @@ import java.util.Map;
import java.util.Map.Entry;
import java.util.Set;
+import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.Intensity;
import ac.ed.lurg.landuse.IrrigationItem;
@@ -15,8 +16,10 @@ import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.GamsLandCoverType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.YieldType;
+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.ed.lurg.yield.YieldRaster;
import ac.ed.lurg.yield.YieldResponsesItem;
import ac.sac.raster.IntegerRasterItem;
@@ -101,7 +104,7 @@ public class GamsRasterOptimiser {
//checkedTotalAreas(rasterInputData.getPreviousAreas(), "old");
//checkedTotalAreas(newAreaRaster, "new");
- LandCoverChangeItem landCoverChangeSet = gamsOutput.getLandCoverChange();
+ //TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges = gamsOutput.getLandCoverChanges();
for (Map.Entry<Integer, LandUseItem> entry : gamsOutput.getLandUses().entrySet()) {
Integer locId = entry.getKey();
@@ -116,54 +119,39 @@ public class GamsRasterOptimiser {
keys.add(mapEntry.getKey());
}
- // Allocate land cover change to grid cells
+
+ RasterSet<LandUseItem> landUseItemsForLocation = newLandUseRaster.createSubsetForKeys(keys);
- Map<LandCoverType, Map<LandCoverType, Double>> landCoverChange = landCoverChangeSet.getLandCoverChangeForLocation(locId);
+ DoubleMap<LandCoverType, LandCoverType, Double> landCoverChange = gamsOutput.getLandCoverChanges().getInnerMap(locId);
/*
- for (LandCoverType fromLC : LandCoverType.getConvertibleTypes()) {
- for (LandCoverType toLC : LandCoverType.getConvertibleTypes()) {
- double change = landCoverChange.getLandCoverChange(fromLC, toLC, locId);
- allocAllLandCropsTop(newLandUseRaster, keys, toLC, fromLC, change, locId);
- }
+ if (landCoverChange == null) {
+ landCoverChange = new DoubleMap<LandCoverType, LandCoverType, Double>();
}
*/
- for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap: landCoverChange.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);
- }
- }
-
-
-/*
- double gamsPastureChange = newLandUseAggItem.getLandCoverArea(LandCoverType.PASTURE) - prevLandUseAggItem.getLandCoverArea(LandCoverType.PASTURE);
- double gamsCroplandChange = newLandUseAggItem.getLandCoverArea(LandCoverType.CROPLAND) - prevLandUseAggItem.getLandCoverArea(LandCoverType.CROPLAND);
- RasterSet<LandUseItem> landUseItemsForLocation = newLandUseRaster.createSubsetForKeys(keys);
-
- if (DEBUG) {
- checkedTotalAreas(landUseItemsForLocation, locId + " before");
- LogWriter.println("pastureChange " + gamsPastureChange);
- LogWriter.println("croplandChange " + gamsCroplandChange);
- }
-
- double prevManagedForest = 0, prevUnmanagedForest = 0, prevNatural = 0, protectedAreasPastureChange = 0, protectedAreasCroplandChange = 0;
+ double prevTimberForest = 0, prevCarbonForest = 0, prevUnmanagedForest = 0, prevNatural = 0;
+ double protectedAreasPastureChange = 0, protectedAreasCroplandChange = 0;
for (LandUseItem luItem: landUseItemsForLocation.values()) {
- prevManagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.MANAGED_FOREST);
+ prevTimberForest += luItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST);
+ prevCarbonForest += luItem.getUnprotectedLandCoverArea(LandCoverType.CARBON_FOREST);
prevUnmanagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST);
prevNatural += luItem.getTotalUnprotectedNatural();
}
- double prevForestTotal = prevManagedForest + prevUnmanagedForest;
+ double prevForestTotal = prevTimberForest + + prevCarbonForest + prevUnmanagedForest;
double prevForestToNaturalFraction = (prevNatural > 0) ? prevForestTotal / prevNatural : 0;
- double prevForestManagedFraction = (prevForestTotal > 0) ? prevManagedForest / prevForestTotal : 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(prevForestManagedFraction > 1.0) prevForestManagedFraction = 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) {
@@ -183,16 +171,20 @@ public class GamsRasterOptimiser {
double pastureShortfall = excessAgriculture * prevPastureFract;
double croplandShortfall = excessAgriculture * (1 - prevPastureFract);
- protectedLandShortfall += luItem.moveAreas(LandCoverType.MANAGED_FOREST, LandCoverType.PASTURE,
- pastureShortfall * prevForestToNaturalFraction * prevForestManagedFraction);
+ 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 * (1 - prevForestManagedFraction));
+ pastureShortfall * prevForestToNaturalFraction * prevForestUnmanagedFraction);
protectedLandShortfall += luItem.moveAreas(LandCoverType.OTHER_NATURAL, LandCoverType.PASTURE,
pastureShortfall * (1 - prevForestToNaturalFraction));
- protectedLandShortfall += luItem.moveAreas(LandCoverType.MANAGED_FOREST, LandCoverType.CROPLAND,
- croplandShortfall * prevForestToNaturalFraction * prevForestManagedFraction);
+ 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 * (1 - prevForestManagedFraction));
+ croplandShortfall * prevForestToNaturalFraction * prevForestUnmanagedFraction);
protectedLandShortfall += luItem.moveAreas(LandCoverType.OTHER_NATURAL, LandCoverType.CROPLAND,
croplandShortfall * (1 - prevForestToNaturalFraction));
@@ -209,58 +201,28 @@ public class GamsRasterOptimiser {
if (DEBUG) {
LogWriter.println("protectedAreasPastureChange " + protectedAreasPastureChange + " protectedAreasCroplandChange " + protectedAreasCroplandChange);
}
-
-
- }
-
- double totalPastureChange = protectedAreasPastureChange + gamsPastureChange;
- double totalCroplandChange = protectedAreasCroplandChange + gamsCroplandChange;
-
- double pastureFromCrop = 0;
- double pastureFromNatural = 0;
- double cropFromNatural = 0;
-
- // Gained pasture and lost cropland
- if (totalPastureChange > 0 && totalCroplandChange < 0) {
- pastureFromCrop = Math.min(totalPastureChange, -totalCroplandChange);
- // Gain in pasture greater than loss in cropland
- if (totalPastureChange > -totalCroplandChange)
- pastureFromNatural = totalPastureChange+totalCroplandChange; // take pasture gain not met by cropland loss from natural
- // Gain in pasture less than loss in cropland
- else
- cropFromNatural = totalPastureChange+totalCroplandChange; // (-ve) add cropland loss not met by pasture gain to natural
- }
- // Gained cropland and lost pasture
- else if (totalPastureChange < 0 && totalCroplandChange > 0) {
- pastureFromCrop = -Math.min(-totalPastureChange, totalCroplandChange);
-
- if (totalCroplandChange > -totalPastureChange)
- cropFromNatural = totalPastureChange+totalCroplandChange;
- else
- pastureFromNatural = totalPastureChange+totalCroplandChange;
- }
- else {
- pastureFromNatural = totalPastureChange;
- cropFromNatural = totalCroplandChange;
- }
+ double totalNaturalArea = newLandUseAggItem.getTotalNatural();
+
+
+ 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);
+ }
- if (DEBUG) {
- LogWriter.println("pastureFromCrop " + pastureFromCrop);
- LogWriter.println("pastureFromNatural " + pastureFromNatural);
- LogWriter.println("cropFromNatural " + cropFromNatural);
+ }
+
+ // 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);
+ }
}
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.PASTURE, LandCoverType.MANAGED_FOREST, pastureFromNatural * prevForestToNaturalFraction * prevForestManagedFraction, locId);
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.PASTURE, LandCoverType.UNMANAGED_FOREST, pastureFromNatural * prevForestToNaturalFraction * (1-prevForestManagedFraction), locId);
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.PASTURE, LandCoverType.OTHER_NATURAL, pastureFromNatural * (1-prevForestToNaturalFraction), locId);
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.PASTURE, LandCoverType.CROPLAND, pastureFromCrop, locId);
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.CROPLAND, LandCoverType.MANAGED_FOREST, cropFromNatural * prevForestToNaturalFraction * prevForestManagedFraction, locId);
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.CROPLAND, LandCoverType.UNMANAGED_FOREST, cropFromNatural * prevForestToNaturalFraction * (1-prevForestManagedFraction), locId);
- allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.CROPLAND, LandCoverType.OTHER_NATURAL, cropFromNatural * (1-prevForestToNaturalFraction), locId);
-
- if (DEBUG) checkedTotalAreas(landUseItemsForLocation, locId + " after");
-*/
for (RasterKey key : keys) {
LandUseItem newLandUseItem = newLandUseRaster.get(key);
if (newLandUseItem != null) {
@@ -272,10 +234,9 @@ public class GamsRasterOptimiser {
}
}
-
+ if (DEBUG) checkedTotalAreas(landUseItemsForLocation, locId + " after");
}
-
return newLandUseRaster;
}
@@ -459,8 +420,6 @@ public class GamsRasterOptimiser {
// Carbon flux
if (!aggCFlux.checkForKeys(prevLC, newLC)) { // if not seen yet
- //LogWriter.println("prevLC: " + prevLC + ", newLC: " + newLC);
- //LogWriter.println("Data: " + carbonFluxItem.getCarbonFluxes());
aggCFlux.setCarbonFlux(prevLC, newLC, carbonFluxItem.getCarbonFlux(prevLC, newLC));
} else {
aggCFlux.setCarbonFlux(prevLC, newLC, aggregateMean(aggCFlux.getCarbonFlux(prevLC, newLC), suitableAreaSoFar,
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index f6774523f6f5c2082b5512e23dae52df8e5fab92..dff462d47d8a7cabf0baa1baf370673572bb9940 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -1,9 +1,11 @@
package ac.ed.lurg.landuse;
import java.io.Serializable;
+import java.util.ArrayList;
import java.util.Collection;
import java.util.HashMap;
import java.util.Map;
+import java.util.Map.Entry;
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.types.CropToDouble;
@@ -19,7 +21,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>();
private Map<LandCoverType, Double> unprotectedAreas = new HashMap<LandCoverType, Double>();
- //private Map<GamsLandCoverType, Double> gamsLandCoverAreas = new HashMap<GamsLandCoverType, Double>();
+ private Map<LandCoverType, Map<Integer, Double>> minimumLandCover = new HashMap<LandCoverType, Map<Integer, Double>>(); // Minimum land cover by year of expiry
private double protectedArea; //protected area in Mha
private double unavailableArea; //area unavailable due to altitude etc
private double suitableArea;
@@ -433,37 +435,17 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return total;
}
- /*
- public boolean checkForUnprotected(LandCoverType landCover) {
- return unprotectedAreas.containsKey(landCover);
- }
- */
- /*
- public double getGamsLandCoverArea(GamsLandCoverType landCoverType) {
- return gamsLandCoverAreas.get(landCoverType);
- }
-
- public void setGamsLandCoverArea(GamsLandCoverType c, double d) {
- if (Double.isNaN(d) || Double.isInfinite(d))
- throw new RuntimeException("AreasItem for " + c + " is " + d);
-
- double landCover = (d < 0.0) ? 0.0 : d;
-
- gamsLandCoverAreas.put(c, landCover);
+ public double getMinimumLandCover(LandCoverType lcType, Integer year) {
+ double area = 0;
+ Map<Integer, Double> yearMap = minimumLandCover.get(lcType);
+ for (Entry<Integer, Double> entry : yearMap.entrySet()) {
+ if (entry.getKey() > year) {
+ area += entry.getValue();
+ }
+ }
+ return area;
}
-
- // move areas from one land cover to another, return any residual not possible
- public void moveGamsAreas(GamsLandCoverType toType, GamsLandCoverType fromType, double changeReq) {
-
- double prevTo = getGamsLandCoverArea(toType);
- double prevFrom = getGamsLandCoverArea(fromType);
-
- setGamsLandCoverArea(toType, prevTo + changeReq);
- setGamsLandCoverArea(fromType, prevFrom - changeReq);
-
- }
- */
@Override
public String toString() {
return "LandUseItem: [landCoverAreas=" + landCoverAreas + ", protectedArea=" + protectedArea + ", unavailableArea=" + unavailableArea + "]";
diff --git a/src/ac/ed/lurg/output/LpjgOutputer.java b/src/ac/ed/lurg/output/LpjgOutputer.java
index 5ba5f785f1d232c49ebfad835ff38d3565f6bfeb..2562d1e2826644ebe319c99500e8283b8702a2cd 100644
--- a/src/ac/ed/lurg/output/LpjgOutputer.java
+++ b/src/ac/ed/lurg/output/LpjgOutputer.java
@@ -172,7 +172,7 @@ public class LpjgOutputer extends AbstractLandUseOutputer {
try {
String landCoverFileName = outputDir.getPath() + File.separator + "GamsLandCoverFract.txt";
landCoverWriter = new BufferedWriter(new FileWriter(landCoverFileName, false));
- landCoverWriter.write("Lon Lat Year CROPLAND PASTURE TIMBER CARBON NATURAL BARREN URBAN");
+ landCoverWriter.write("Lon Lat Year CROPLAND PASTURE TIMBER CARBON UNMANAGED NATURAL BARREN URBAN");
landCoverWriter.newLine();
for (Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) {
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index 47c1563360ac3c5321c0bccbc8f9f70169da7a81..84277884348507591164cb63df228a39c471e8b3 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -65,5 +65,17 @@ public enum LandCoverType {
return convertibleTypes;
}
+
+ public static Collection<LandCoverType> getProtectibleTypes() {
+
+ Collection<LandCoverType> protectibleTypes = new HashSet<LandCoverType>();
+
+ for (LandCoverType c : values())
+ if (c.isProtectable)
+ protectibleTypes.add(c);
+
+ return protectibleTypes;
+
+ }
}
diff --git a/src/ac/ed/lurg/utils/DoubleMap.java b/src/ac/ed/lurg/utils/DoubleMap.java
new file mode 100644
index 0000000000000000000000000000000000000000..9c2762b3a7ea918cd876a173527655af034d5e2c
--- /dev/null
+++ b/src/ac/ed/lurg/utils/DoubleMap.java
@@ -0,0 +1,51 @@
+package ac.ed.lurg.utils;
+
+import java.util.HashMap;
+import java.util.Map;
+
+public class DoubleMap<K, L, V extends Double> {
+
+ Map<K, Map<L, Double>> theMap = new HashMap<K, Map<L, Double>>();
+
+ public DoubleMap() {
+
+ }
+
+ public void put(K key1, L key2, double value) {
+ if (theMap.containsKey(key1)) { // check if inner map already exists
+ theMap.get(key1).put(key2, value); // adds new entry or overwrites existing if exists
+ } else { // add new complete entry
+ Map<L, Double> innerMap = new HashMap<L, Double>();
+ innerMap.put(key2, value);
+ theMap.put(key1, innerMap);
+ }
+ }
+
+ public Double get(K key1, L key2) {
+ if (theMap.containsKey(key1)) { // check if first key exits
+ return theMap.get(key1).get(key2); // will return null if second key doesn't exist
+ } else {
+ return null; // will return null if first key doesn't exist
+ }
+ }
+
+ public Map<L, Double> getInnerMap(K key) {
+ return theMap.get(key);
+
+ }
+
+ public Map<K, Map<L, Double>> getMap() {
+ return theMap;
+ }
+
+ public void addTo(K key1, L key2, double value) {
+ if (theMap.containsKey(key1)) {
+ double current = get(key1, key2);
+ put(key1, key2, current + value);
+ } else {
+ put(key1, key2, value);
+ }
+
+
+ }
+}
diff --git a/src/ac/ed/lurg/utils/NestedMap.java b/src/ac/ed/lurg/utils/NestedMap.java
new file mode 100644
index 0000000000000000000000000000000000000000..54f72c9b9a326241846770277ffe9878286debe8
--- /dev/null
+++ b/src/ac/ed/lurg/utils/NestedMap.java
@@ -0,0 +1,51 @@
+package ac.ed.lurg.utils;
+
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.Map;
+import java.util.List;
+
+public class NestedMap {
+
+ Map<Object, Object> theMap;
+ Integer depth;
+
+ NestedMap() {
+ theMap = new HashMap<Object, Object>();
+ }
+
+ public void put(List<Object> keys, Object value) {
+
+
+
+ Collections.reverse(keys);
+ theMap.put(keys.get(keys.size() - 1), createNestedMap(keys.subList(0, keys.size() - 1), value));
+ }
+
+ private Map<Object, Object> createNestedMap(List<Object> keys, Object value) {
+ Map<Object, Object> valueMap = new HashMap<Object, Object>();
+ valueMap.put(keys.get(0), value);
+ Map<Object, Object> higherMap;
+ if (keys.size() > 1) {
+ higherMap = createNestedMap(keys.subList(1, keys.size()), valueMap);
+ return higherMap;
+ } else {
+ return valueMap;
+ }
+ }
+ /*
+ private Map<Object, Object> getLowestLevel(List<Object> keys) {
+ if (theMap.containsKey(keys.get(0))) {
+ getLowestLevel(keys.subList(1, keys.size()));
+ } else {
+ return theMap.get(key);
+ }
+ }
+ */
+ /*
+ private Map<Object, Object> getMapForKeys(List<Object> keys) {
+ Map<Object, Object> map = new HashMap<Object, Object>();
+
+ }
+ */
+}
diff --git a/src/ac/ed/lurg/utils/TripleMap.java b/src/ac/ed/lurg/utils/TripleMap.java
new file mode 100644
index 0000000000000000000000000000000000000000..78a0dbb9136a1ec786771786055c6f024e4d85b7
--- /dev/null
+++ b/src/ac/ed/lurg/utils/TripleMap.java
@@ -0,0 +1,37 @@
+package ac.ed.lurg.utils;
+
+import java.util.HashMap;
+import java.util.Map;
+
+public class TripleMap<K, L, M, V extends Double> {
+
+ Map<K, DoubleMap<L, M, V >> theMap = new HashMap<K, DoubleMap<L, M, V>>();
+
+ public TripleMap() {
+
+ }
+
+ public void put(K key1, L key2, M key3, double value) {
+ if (theMap.containsKey(key1)) {
+ theMap.get(key1).put(key2, key3, value);
+ } else {
+ DoubleMap<L, M, V> innerMap = new DoubleMap<L, M, V>();
+ innerMap.put(key2, key3, value);
+ theMap.put(key1, innerMap);
+ }
+ }
+
+
+ public Double get(K key1, L key2, M key3) {
+ if (theMap.containsKey(key1)) {
+ return theMap.get(key1).get(key2, key3);
+ } else {
+ return null;
+ }
+ }
+
+ public DoubleMap<L, M, V> getInnerMap(K key) {
+ return theMap.get(key);
+
+ }
+}