Initial work on reading and processing carbon flux and wood yield data.

Further forest work on GAMS code.

GAMS/IntExtOpt.gms

 
- 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);

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";

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) {

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;
 	}

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());

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);

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) {

src/ac/ed/lurg/landuse/CarbonFluxReader.java

+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"));
+		
+                                                                                
+	}
+}

src/ac/ed/lurg/landuse/WoodYieldRasterSet.java

+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;
+	}
+}

src/ac/ed/lurg/landuse/WoodYieldReader.java

+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"));
+		                                                                          
+	}
+}

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),