From 2752f62159b63cf150c56576e47da522a6a66648 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Tue, 9 Feb 2021 15:41:56 +0000
Subject: [PATCH] Initial work on reading and processing carbon flux and wood
yield data. Further forest work on GAMS code.
---
GAMS/IntExtOpt.gms | 59 +++++++++++--------
src/ac/ed/lurg/ModelConfig.java | 4 +-
src/ac/ed/lurg/ModelMain.java | 25 +++++++-
src/ac/ed/lurg/country/CountryAgent.java | 12 ++--
.../ed/lurg/country/CountryAgentManager.java | 11 +++-
.../country/gams/GamsLocationOptimiser.java | 4 +-
.../country/gams/GamsRasterOptimiser.java | 41 +++++++++++--
src/ac/ed/lurg/landuse/CarbonFluxReader.java | 47 +++++++++++++++
.../ed/lurg/landuse/WoodYieldRasterSet.java | 28 +++++++++
src/ac/ed/lurg/landuse/WoodYieldReader.java | 45 ++++++++++++++
src/ac/ed/lurg/types/LandCoverType.java | 3 +-
11 files changed, 238 insertions(+), 41 deletions(-)
create mode 100644 src/ac/ed/lurg/landuse/CarbonFluxReader.java
create mode 100644 src/ac/ed/lurg/landuse/WoodYieldRasterSet.java
create mode 100644 src/ac/ed/lurg/landuse/WoodYieldReader.java
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 33495b6e..97955b75 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -1,5 +1,5 @@
- SET all_types / monogastrics, ruminants, cereals, oilcropspulses, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside, carbonForest, timberForest, natural /;
+ SET all_types / monogastrics, ruminants, cereals, oilcropspulses, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside/;
SET crop(all_types) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside /;
SET crop_less_pasture(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, setaside/;
@@ -8,15 +8,19 @@
SET feed_crop(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, pasture/;
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 wood_producing(all_types) / timberForest, carbonForest, natural/;
-* SET forest(all_types) / timberForest, carbonForest /;
- ALIAS (all_types, lu_type_previous);
+ SET land_cover / cropland, pasture, timberForest, carbonForest, natural /;
+ SET wood_producing / timberForest, carbonForest, natural /;
+ ALIAS (land_cover, land_cover_before);
+ ALIAS (land_cover, land_cover_start);
+ ALIAS (land_cover, land_cover_after);
+
SET location;
PARAMETER suitableLandArea(location) areas of land in Mha;
PARAMETER previousCropArea(crop, location) areas for previous timestep in Mha;
- PARAMETER previousNonCropArea(wood_producing, location);
+ PARAMETER previousLandCoverArea(land_cover_before, land_cover_start, location) land cover area at timestep t-1 split by land cover at t-2;
+ PARAMETER minimumLandCoverArea(land_cover, location) minimum land cover area to constrain conversion;
PARAMETER previousFertIntensity(crop, location);
PARAMETER previousIrrigIntensity(crop, location);
PARAMETER previousOtherIntensity(crop, location);
@@ -59,8 +63,8 @@
SCALAR domesticPriceMarkup factor price increased from cost of production;
SCALAR maxLandExpansionRate max rate of country land expansion;
- PARAMETER carbonFluxRate(all_types, location) carbon flux;
- PARAMETER woodYield(wood_producing, location) wood yield;
+ PARAMETER carbonFluxRate(land_cover_before, land_cover_start, location) carbon flux;
+ PARAMETER woodYield(land_cover_before, land_cover_start, location) wood yield;
SCALAR carbonPrice price of carbon;
SCALAR woodPrice price of wood and timber;
@@ -115,7 +119,7 @@ $gdxin
VARIABLES
cropArea(crop, location) total area for crops - Mha
- nonCropArea(wood_producing, location) forest and natural area
+ landCoverArea(land_cover_start, land_cover_after, location) forest and natural area
fertI(crop, location) fertilizer intensity for each crop - factor between 0 and 1
irrigI(crop, location) irrigation intensity for each crop - factor between 0 and 1
otherIntensity(crop, location)
@@ -159,13 +163,7 @@ $gdxin
MIN_NET_IMPORT_CONSTRAINT(import_crop) constraint on min net imports
PASTURE_IMPORT_CONSTRAINT constraint to not import pasture
PASTURE_EXPORT_CONSTRAINT
- TIMBER_IMPORT_CONSTRAINT constraint to not import timber (for now)
- TIMBER_EXPORT_CONSTRAINT
- CARBON_IMPORT_CONSTRAINT constraint to not import carbon
- CARBON_EXPORT_CONSTRAINT
- NATURAL_IMPORT_CONSTRAINT constraint to not import natural
- NATURAL_EXPORT_CONSTRAINT
- IRRIGATION_CONSTRAINT(location) constraint no water usage
+ IRRIGATION_CONSTRAINT(location) constraint on water usage
NET_SUPPLY_EQ(crop) calc net supply for crops
AGRI_LAND_EXPANSION_CALC(location) calc agriLandExpansion
CROP_INCREASE_CALC(location)
@@ -173,6 +171,11 @@ $gdxin
PASTURE_INCREASE_CONV_CALC(location)
PASTURE_DECREASE_CONV_CALC(location)
PASTURE_TOTAL_CHANGE_CONSTRAINT(location)
+ PREVIOUS_LAND_COVER_CONSTRAINT(land_cover_start, location) constraint so that previous land cover at this timestep equals land cover at previous timestep
+ CROPLAND_CONSTRAINT(location) cropland area equals sum of all crop areas
+ PASTURE_CONSTRAINT(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
WOOD_HARVEST_CALC(location) calc wood harvested
CARBON_FLUX_CALC(location) calc carbon flux
COST_EQ total cost objective function;
@@ -212,8 +215,9 @@ $gdxin
SETASIDE_AREA_CALC(location) .. cropArea('setaside', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) * setAsideRate;
- TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate) + cropArea('pasture', location) +
- sum(wood_producing, nonCropArea(wood_producing, location));
+ TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate) + cropArea('pasture', location);
+
+
MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
@@ -223,13 +227,6 @@ $gdxin
MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
- TIMBER_IMPORT_CONSTRAINT .. importAmount('timberForest') =E= 0;
- TIMBER_EXPORT_CONSTRAINT .. exportAmount('timberForest') =E= 0;
- CARBON_IMPORT_CONSTRAINT .. importAmount('carbonForest') =E= 0;
- CARBON_EXPORT_CONSTRAINT .. exportAmount('carbonForest') =E= 0;
- NATURAL_IMPORT_CONSTRAINT .. importAmount('natural') =E= 0;
- NATURAL_EXPORT_CONSTRAINT .. exportAmount('natural') =E= 0;
-
IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) * (1.0 - unhandledCropRate) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * cropArea(crop, location));
AGRI_LAND_EXPANSION_CALC(location) .. agriLandExpansion(location) =G= sum(crop, cropArea(crop, location) - previousCropArea(crop, location));
@@ -240,9 +237,19 @@ $gdxin
PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(cropArea('pasture', location) - previousCropArea('pasture', location));
PASTURE_TOTAL_CHANGE_CONSTRAINT(location) .. pastureIncrease(location) -pastureDecrease(location) =G= cropArea('pasture', location) - previousCropArea('pasture', location);
- WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum(wood_producing, nonCropArea(wood_producing, location) * woodYield(wood_producing, location));
+ PREVIOUS_LAND_COVER_CONSTRAINT(land_cover_start, location) .. sum(land_cover_after, landCoverArea(land_cover_start, land_cover_after, location)) =E= sum(land_cover_before, previousLandCoverArea(land_cover_before, land_cover_start, location));
+
+ CROPLAND_CONSTRAINT(location) .. sum(land_cover_start, landCoverArea(land_cover_start, 'cropland', location)) =E= sum(crop, cropArea(crop, location));
+
+ PASTURE_CONSTRAINT(location) .. sum(land_cover_start, landCoverArea(land_cover_start, 'pasture', location)) =E= cropArea('pasture', location);
+
+ MINIMUM_LAND_COVER_CONSTRAINT(location, land_cover_after) .. sum(land_cover_start, landCoverArea(land_cover_start, land_cover_after, location)) =G= minimumLandCover(land_cover_after, location);
+
+ TOTAL_LAND_COVER_CONSTRAINT(location) .. sum((land_cover_start, land_cover_after), landCoverArea(land_cover_start, land_cover_after, location)) =E= suitableArea(location);
+
+ WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((land_cover_start, land_cover_after), landCoverArea(land_cover_start, land_cover_after, location) * woodYield(land_cover_start, land_cover_after, location));
- CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum(wood_producing, nonCropArea(wood_producing, location) * carbonFluxRate(wood_producing, location));
+ CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_start, land_cover_after), landCoverArea(land_cover_start, land_cover_after, location) * carbonFluxRate(land_cover_start, land_cover_after, location));
* CROPLAND_INCREASE_CONSTRAINT(location) .. cropIncrease(location)/(1.0 - unhandledCropRate) =L= maxLandChangeRate * ( suitableLandArea(location) -sum(crop_less_pasture, previousCropArea(crop_less_pasture, location))/*(1.0 - unhandledCropRate) - previousCropArea('pasture', location) );
* PASTURE_DECREASE_CONSTRAINT(location) .. pastureDecrease(location) =L= maxLandChangeRate * previousCropArea('pasture', location);
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index e289c2a0..99f61ca6 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -243,7 +243,9 @@ public class ModelConfig {
public static final String YIELDSHOCK_MAP_DIR = SPATIAL_DATA_DIR + File.separator + "yieldshockmaps";
public static final String YIELDSHOCKS_PARAMETER_FILE = getProperty("YIELDSHOCKS_PARAMETER_FILE", OUTPUT_DIR + File.separator+ "yieldshocks.csv");
public static final String PRICESHOCKS_PARAMETER_FILE = getProperty("PRICESHOCKS_PARAMETER_FILE", OUTPUT_DIR + File.separator+ "priceshocks.csv");
- public static final String EXPORT_RESTRICTIONS_FILE = getProperty("EXPORT_RESTRICTIONS_FILE", OUTPUT_DIR + File.separator+ "exportrestictions.csv");;
+ public static final String EXPORT_RESTRICTIONS_FILE = getProperty("EXPORT_RESTRICTIONS_FILE", OUTPUT_DIR + File.separator+ "exportrestictions.csv");
+ public static final String CARBON_FLUX_FILE = SPATIAL_DATA_DIR + File.separator + "carbon_flux_";
+ public static final String WOOD_YIELD_FILE = SPATIAL_DATA_DIR + File.separator + "wood_yield_";
// Output
public static final String LAND_COVER_OUTPUT_FILE = OUTPUT_DIR + File.separator + "lc.txt";
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 070cfc7d..5a678e9f 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -26,6 +26,8 @@ import ac.ed.lurg.demand.CalorieManager;
import ac.ed.lurg.demand.DemandManagerFromFile;
import ac.ed.lurg.demand.DemandManagerSSP;
import ac.ed.lurg.demand.ElasticDemandManager;
+import ac.ed.lurg.landuse.CarbonFluxRasterSet;
+import ac.ed.lurg.landuse.CarbonFluxReader;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.CropUsageReader;
import ac.ed.lurg.landuse.FPUManager;
@@ -40,6 +42,8 @@ import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.MaxCropAreaReader;
import ac.ed.lurg.landuse.ProtectedAreasReader;
import ac.ed.lurg.landuse.RunOffReader;
+import ac.ed.lurg.landuse.WoodYieldRasterSet;
+import ac.ed.lurg.landuse.WoodYieldReader;
import ac.ed.lurg.output.LandUseOutputer;
import ac.ed.lurg.output.LpjgOutputer;
import ac.ed.lurg.types.CommodityType;
@@ -69,6 +73,8 @@ public class ModelMain {
private InternationalMarket internationalMarket;
private IrrigationRasterSet currentIrrigationData;
+ private CarbonFluxRasterSet currentCarbonFluxData;
+ private WoodYieldRasterSet currentWoodYieldData;
private RasterSet<LandUseItem> globalLandUseRaster;
private RasterSet<IntegerRasterItem> clusterIdRaster;
@@ -136,7 +142,10 @@ public class ModelMain {
double gen2EcDDemand = demandManager.getSecondGenBioenergyDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
double gen2Increase = (gen2EcDDemand>previousGen2EcDDemand) ? gen2EcDDemand - previousGen2EcDDemand : 0.0;
- countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase);
+ CarbonFluxRasterSet currentCarbonFluxData = getCarbonFluxData(timestep);
+ WoodYieldRasterSet currentWoodYieldData = getWoodYieldData(timestep);
+
+ countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, currentWoodYieldData);
internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, timestep);
// loop for iterations. Could check within a tolerance using internationalMarket.findMaxPriceDiff, not doing so as doesn't find a solution due to inelastic consumption
@@ -619,6 +628,20 @@ public class ModelMain {
fixedIrrigData.calcIrrigConstraintOffsets(); // should have everything we need to calc offset between Elliott and LPJ data
return fixedIrrigData;
}
+
+ /** Get carbon flux data */
+ private CarbonFluxRasterSet getCarbonFluxData(Timestep timestep) {
+ CarbonFluxRasterSet carbonFluxData = new CarbonFluxRasterSet(desiredProjection);
+ new CarbonFluxReader(carbonFluxData).getRasterDataFromFile(ModelConfig.CARBON_FLUX_FILE + timestep.getYieldYear() + ".out");
+ return carbonFluxData;
+ }
+
+ /** Get carbon wood yield data */
+ private WoodYieldRasterSet getWoodYieldData(Timestep timestep) {
+ WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
+ new WoodYieldReader(woodYieldData).getRasterDataFromFile(ModelConfig.WOOD_YIELD_FILE + timestep.getYieldYear() + ".out");
+ return woodYieldData;
+ }
/** Ugly in situ update of currentIrrigationData, better if IrrigationRasterSets were handled more immutably */
private void getUpdateIrrigationData(Timestep timestep, YieldRaster yieldSurfaces) {
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index c3accb6c..a723454c 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -17,9 +17,11 @@ import ac.ed.lurg.country.gams.GamsRasterInput;
import ac.ed.lurg.country.gams.GamsRasterOptimiser;
import ac.ed.lurg.country.gams.GamsRasterOutput;
import ac.ed.lurg.demand.AbstractDemandManager;
+import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.WoodYieldItem;
import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
@@ -102,7 +104,7 @@ public class CountryAgent extends AbstractCountryAgent {
}
public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData,
- Map<CropType, GlobalPrice> worldPrices, double globalGen2EcIncrease) {
+ Map<CropType, GlobalPrice> worldPrices, double globalGen2EcIncrease, RasterSet<CarbonFluxItem> carbonFluxData, RasterSet<WoodYieldItem> woodYieldData) {
// project demand
calculateCountryPricesAndDemand(worldPrices, false);
@@ -123,7 +125,7 @@ public class CountryAgent extends AbstractCountryAgent {
yieldClusters = calcYieldClusters(irrigData, countryYieldSurfaces); // this should only be on the first timestep
// optimize areas and intensity
- GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease);
+ GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, carbonFluxData, woodYieldData);
GamsRasterOptimiser opti = new GamsRasterOptimiser(input, yieldClusters);
LogWriter.println("Running " + country.getName() + ", currentTimestep " + currentTimestep);
@@ -155,7 +157,8 @@ public class CountryAgent extends AbstractCountryAgent {
}
}
- private GamsRasterInput getGamsRasterInput(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces, double gen2EcIncrease) {
+ private GamsRasterInput getGamsRasterInput(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces, double gen2EcIncrease,
+ RasterSet<CarbonFluxItem> carbonFluxData, RasterSet<WoodYieldItem> woodYieldData) {
double allowedImportChange;
if (currentTimestep.isInitialTimestep() || (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION)) { // initialisation time-step
@@ -200,7 +203,8 @@ public class CountryAgent extends AbstractCountryAgent {
GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints,
previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates);
- GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, countryLevelInputs);
+ GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData,
+ carbonFluxData, woodYieldData, countryLevelInputs);
return input;
}
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index 8340d27a..a7395b36 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -15,10 +15,14 @@ import ac.ed.lurg.Timestep;
import ac.ed.lurg.country.crafty.CraftyCountryAgent;
import ac.ed.lurg.country.crafty.CraftyProdManager;
import ac.ed.lurg.demand.AbstractDemandManager;
+import ac.ed.lurg.landuse.CarbonFluxItem;
+import ac.ed.lurg.landuse.CarbonFluxRasterSet;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.IrrigationRasterSet;
import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.WoodYieldItem;
+import ac.ed.lurg.landuse.WoodYieldRasterSet;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.YieldRaster;
@@ -96,7 +100,8 @@ public class CountryAgentManager {
return countryAgents;
}
- public void determineProductionForAll(Timestep timestep, YieldRaster yieldSurfaces, IrrigationRasterSet currentIrrigationData, double gen2Increase) {
+ public void determineProductionForAll(Timestep timestep, YieldRaster yieldSurfaces, IrrigationRasterSet currentIrrigationData, double gen2Increase,
+ CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData) {
for (AbstractCountryAgent aca : countryAgents) {
aca.setCurrentTimestep(timestep);
@@ -107,9 +112,11 @@ public class CountryAgentManager {
Collection<RasterKey> countryKeys = countryBoundaryRaster.getKeysFor(ca.getCountry());
YieldRaster countryYieldSurfaces = yieldSurfaces.createSubsetForKeys(countryKeys);
RasterSet<IrrigationItem> irrigData = currentIrrigationData.createSubsetForKeys(countryKeys);
+ RasterSet<CarbonFluxItem> carbonFluxData = currentCarbonFluxData.createSubsetForKeys(countryKeys);
+ RasterSet<WoodYieldItem> woodYieldData = currentWoodYieldData.createSubsetForKeys(countryKeys);
try {
- ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase);
+ ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase, carbonFluxData, woodYieldData);
// update global rasters
globalLandUseRaster.putAll(ca.getLandUses());
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 55c27ada..256a88a5 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -87,7 +87,9 @@ public class GamsLocationOptimiser {
}
if (DEBUG) LogWriter.println("\nPrevious crop and land areas");
- GAMSParameter prevCropP = inDB.addParameter("previousArea", 2);
+ GAMSParameter prevCropP = inDB.addParameter("previousCropArea", 2);
+ GAMSParameter prevLandCoverP = inDB.addParameter("previousLandCoverArea", 3);
+ GAMSParameter minLandCoverP = inDB.addParameter("minimumLandCoverArea", 3);
GAMSParameter prevFertIP = inDB.addParameter("previousFertIntensity", 2);
GAMSParameter prevIrrigIP = inDB.addParameter("previousIrrigIntensity", 2);
GAMSParameter prevOtherIP = inDB.addParameter("previousOtherIntensity", 2);
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index d75c384e..292a38af 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -6,9 +6,11 @@ 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;
import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.WoodYieldItem;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.YieldType;
@@ -331,7 +333,11 @@ public class GamsRasterOptimiser {
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));
@@ -356,6 +362,12 @@ public class GamsRasterOptimiser {
LazyTreeMap<Integer, IrrigationItem> aggregatedIrrigCosts = new LazyTreeMap<Integer, IrrigationItem>() {
protected IrrigationItem createValue() { return new IrrigationItem(); }
};
+ LazyTreeMap<Integer, CarbonFluxItem> aggregatedCarbonFluxes = new LazyTreeMap<Integer, CarbonFluxItem>() {
+ protected CarbonFluxItem createValue() { return new CarbonFluxItem(); }
+ };
+ LazyTreeMap<Integer, WoodYieldItem> aggregatedWoodYields = new LazyTreeMap<Integer, WoodYieldItem>() {
+ protected WoodYieldItem createValue() { return new WoodYieldItem(); }
+ };
int countFound = 0, countMissing = 0;
@@ -379,24 +391,42 @@ public class GamsRasterOptimiser {
}
IrrigationItem irrigItem = irrigRaster.get(key);
+
+ CarbonFluxItem carbonFluxItem = carbonFluxRaster.get(key);
+ WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
int clusterId = mapping.get(key).getInt();
YieldResponsesItem aggYResp = aggregatedYields.lazyGet(clusterId);
LandUseItem aggLandUse = aggregatedAreas.lazyGet(clusterId);
IrrigationItem aggIrig = aggregatedIrrigCosts.lazyGet(clusterId);
-
- // Irrigation cost
+
+ CarbonFluxItem aggCFlux = aggregatedCarbonFluxes.lazyGet(clusterId);
+ WoodYieldItem aggWYield = aggregatedWoodYields.lazyGet(clusterId);
+
landUseItem.updateSuitableArea(rasterInputData.getTimestep().getYear());
double suitableAreaThisTime = landUseItem.getSuitableArea();
double suitableAreaSoFar = aggLandUse.getSuitableArea();
-
+
+ // Irrigation cost
if (irrigItem!= null) {
aggIrig.setIrrigCost( aggregateMean(aggIrig.getIrrigCost(), suitableAreaSoFar, irrigItem.getIrrigCost(), suitableAreaThisTime));
aggIrig.setIrrigConstraint(aggregateMean(aggIrig.getIrrigConstraint(), suitableAreaSoFar, irrigItem.getIrrigConstraint(), suitableAreaThisTime));
//LogWriter.println(id + ", " + key + ", " + avgIrigCost.getIrrigCost() + ", from " + areaSoFar + ", " + suitableAreaThisTime + ", " + irrigCost.getIrrigCost());
}
+
+ // Aggregate carbon fluxes and wood yields
+ for (LandCoverType lc_previous : LandCoverType.values()) {
+ for (LandCoverType lc_new : LandCoverType.values()) {
+ aggCFlux.setCarbonFlux(lc_previous, lc_new, aggregateMean(aggCFlux.getCarbonFlux(lc_previous, lc_new), suitableAreaSoFar,
+ carbonFluxItem.getCarbonFlux(lc_previous, lc_new), suitableAreaThisTime));
+
+ aggWYield.setWoodYield(lc_previous, lc_new, aggregateMean(aggWYield.getWoodYield(lc_previous, lc_new), suitableAreaSoFar,
+ woodYieldItem.getWoodYield(lc_previous, lc_new), suitableAreaThisTime));
+ }
+ }
+
// Crops yields and area fractions
for (CropType crop : CropType.getNonMeatTypes()) {
if (!crop.equals(CropType.SETASIDE)) {
@@ -480,7 +510,8 @@ public class GamsRasterOptimiser {
checkedTotalAreas(landUseRaster, "before");
checkedTotalAreas(aggregatedAreas, "after");
- return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, rasterInputData.getCountryInput());
+ return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts,
+ aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getCountryInput());
}
private void logWarningWithCoord(String message, RasterKey key, YieldRaster yieldRaster, CropType crop) {
diff --git a/src/ac/ed/lurg/landuse/CarbonFluxReader.java b/src/ac/ed/lurg/landuse/CarbonFluxReader.java
new file mode 100644
index 00000000..d23a771c
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/CarbonFluxReader.java
@@ -0,0 +1,47 @@
+package ac.ed.lurg.landuse;
+
+import java.util.Map;
+
+import ac.ed.lurg.types.LandCoverType;
+import ac.sac.raster.AbstractTabularRasterReader;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+
+public class CarbonFluxReader extends AbstractTabularRasterReader<CarbonFluxItem> {
+
+ private static final int MIN_COLS = 18;
+
+ public CarbonFluxReader(RasterSet<CarbonFluxItem> carbonFlux) {
+ super("[ |\t]+", MIN_COLS, carbonFlux);
+ }
+
+ @Override
+ protected void setData(RasterKey key, CarbonFluxItem item, Map<String, Double> rowValues) {
+ item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "c2n"));
+ item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "c2f"));
+ item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "c2x"));
+ item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.PASTURE, getValueForCol(rowValues, "c2p"));
+
+ item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "p2n"));
+ item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "p2f"));
+ item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "p2x"));
+ item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.CROPLAND, getValueForCol(rowValues, "p2c"));
+
+ item.setCarbonFlux(LandCoverType.OTHER_NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "n2p"));
+ item.setCarbonFlux(LandCoverType.OTHER_NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "n2f"));
+ item.setCarbonFlux(LandCoverType.OTHER_NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "n2x"));
+ item.setCarbonFlux(LandCoverType.OTHER_NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "n2p"));
+
+ item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "f2p"));
+ item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "f2n"));
+ item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "f2x"));
+ item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "f2c"));
+
+ item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "x2p"));
+ item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "x2n"));
+ item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "x2f"));
+ item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "x2c"));
+
+
+ }
+}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldRasterSet.java b/src/ac/ed/lurg/landuse/WoodYieldRasterSet.java
new file mode 100644
index 00000000..e330aefa
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/WoodYieldRasterSet.java
@@ -0,0 +1,28 @@
+package ac.ed.lurg.landuse;
+
+import java.util.Collection;
+
+import ac.sac.raster.RasterHeaderDetails;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+
+public class WoodYieldRasterSet extends RasterSet<WoodYieldItem> {
+
+ private static final long serialVersionUID = 3310840203660415050L;
+
+ public WoodYieldRasterSet(RasterHeaderDetails header) {
+ super(header);
+ }
+
+ protected WoodYieldItem createRasterData() {
+ return new WoodYieldItem();
+ }
+
+ // not very efficient, we could keep the mapping of country to area somewhere.
+ @Override
+ public WoodYieldRasterSet createSubsetForKeys(Collection<RasterKey> keys) {
+ WoodYieldRasterSet subsetWoodYieldRaster = new WoodYieldRasterSet(getHeaderDetails());
+ popSubsetForKeys(subsetWoodYieldRaster, keys);
+ return subsetWoodYieldRaster;
+ }
+}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldReader.java b/src/ac/ed/lurg/landuse/WoodYieldReader.java
new file mode 100644
index 00000000..8a2522fa
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/WoodYieldReader.java
@@ -0,0 +1,45 @@
+package ac.ed.lurg.landuse;
+
+import java.util.Map;
+
+import ac.ed.lurg.types.LandCoverType;
+import ac.sac.raster.AbstractTabularRasterReader;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+
+public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem> {
+
+ private static final int MIN_COLS = 18;
+
+ public WoodYieldReader(RasterSet<WoodYieldItem> woodYield) {
+ super("[ |\t]+", MIN_COLS, woodYield);
+ }
+
+ protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
+ item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "c2n"));
+ item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "c2f"));
+ item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "c2x"));
+ item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.PASTURE, getValueForCol(rowValues, "c2p"));
+
+ item.setWoodYield(LandCoverType.PASTURE, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "p2n"));
+ item.setWoodYield(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "p2f"));
+ item.setWoodYield(LandCoverType.PASTURE, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "p2x"));
+ item.setWoodYield(LandCoverType.PASTURE, LandCoverType.CROPLAND, getValueForCol(rowValues, "p2c"));
+
+ item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "n2p"));
+ item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "n2f"));
+ item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "n2x"));
+ item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "n2p"));
+
+ item.setWoodYield(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "f2p"));
+ item.setWoodYield(LandCoverType.TIMBER_FOREST, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "f2n"));
+ item.setWoodYield(LandCoverType.TIMBER_FOREST, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "f2x"));
+ item.setWoodYield(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "f2c"));
+
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "x2p"));
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "x2n"));
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "x2f"));
+ item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "x2c"));
+
+ }
+}
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index 6c0c8428..1bbc1c6f 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -2,7 +2,8 @@ package ac.ed.lurg.types;
public enum LandCoverType {
- MANAGED_FOREST ("managedForest", true),
+ TIMBER_FOREST ("timberForest", true),
+ CARBON_FOREST ("carbonForest", true),
UNMANAGED_FOREST ("unmanagedForest", true),
OTHER_NATURAL ("otherNatural", true),
CROPLAND ("cropland", false),
--
GitLab