From 67df0894a83cda358191c74c273ebb127f8b1da4 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Tue, 16 Dec 2014 19:02:21 +0000
Subject: [PATCH] no message

---
 GAMS/IntExtOpt.gms                            |   2 +-
 src/ac/ed/lurg/ModelConfig.java               |   8 +-
 src/ac/ed/lurg/ModelContext.java              |  16 ---
 src/ac/ed/lurg/ModelMain.java                 |  47 +++++---
 src/ac/ed/lurg/country/CountryAgent.java      | 102 +++++++-----------
 .../lurg/country/gams/GamsCountryInput.java   |  10 +-
 .../country/gams/GamsLocationOptimiser.java   |  10 +-
 .../lurg/country/gams/GamsLocationOutput.java |  17 ++-
 .../country/gams/GamsRasterOptimiser.java     |   3 +-
 .../lurg/country/gams/GamsRasterOutput.java   |  21 ++--
 src/ac/ed/lurg/landuse/CommodityData.java     |  78 ++++++++++++++
 src/ac/sac/raster/RasterKey.java              |  13 +++
 12 files changed, 193 insertions(+), 134 deletions(-)
 delete mode 100644 src/ac/ed/lurg/ModelContext.java
 create mode 100644 src/ac/ed/lurg/landuse/CommodityData.java

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index fb7c287b..013b3b9c 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -89,7 +89,7 @@ $gdxin
 *               (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-5*irrigI(crop, location)) +
 *               (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) * (1 - exp(-fertParam(crop, location)*fertI(crop, location))) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location)))
 *            ) * otherIntensity(crop, location);
- 
+  
  CROP_DEMAND_CONSTRAINT(crop_less_pasture) .. sum(location, area(crop_less_pasture, location) * yield(crop_less_pasture, location)) - feedAmount(crop_less_pasture) =G= 
             demand(crop_less_pasture) - netImportAmount(crop_less_pasture);
  
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 560e9e10..c1eb057b 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -61,7 +61,6 @@ public class ModelConfig {
 		return v==null ? null : Long.valueOf(v);
 	}
 	
-	@SuppressWarnings("unused")
 	private static Double getDoubleProperty(String prop, Double defaultDouble) {
 		Double propValue = getModelConfig().getDoubleProp(prop);
 		return propValue == null ? defaultDouble : propValue;
@@ -84,6 +83,7 @@ public class ModelConfig {
 	public static final boolean SUPPRESS_STD_OUTPUT = getBooleanProperty("SUPPRESS_STD_OUTPUT", Boolean.FALSE);
 
 	public static final String BASE_DIR = getProperty("OUTPUT_DIR", "/Users/peteralexander/Documents/R_Workspace");
+	
 	public static final String TEMP_DIR = BASE_DIR + File.separator + "temp";
 	public static final String OUTPUT_DIR = BASE_DIR + File.separator + "output";
 	public static final String DATA_DIR = BASE_DIR + File.separator + "data";
@@ -97,11 +97,15 @@ public class ModelConfig {
 
 	public static final String COUNTRY_BOUNDARY_FILE = DATA_DIR + File.separator + "country_boundaries.test.asc";
 	public static final String COUNTRY_CODES_FILE = DATA_DIR + File.separator + "country_codes3.csv";
-	
 	public static final String COUNTRY_DATA_FILE = DATA_DIR + File.separator + "country_data.csv";
+	public static final String COMMODITY_DATA_FILE = DATA_DIR + File.separator + "con_prod_c.csv";
 
 	public static final int START_TIMESTEP = getIntProperty("START_TIMESTEP", 0);
 	public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 1);
 	public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2010);
 	
+	public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", 0.05);
+	
+	public static final String SSP_SCENARIO = getProperty("SSP_SCENARIO", "SSP2_v9_130325");
+	
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelContext.java b/src/ac/ed/lurg/ModelContext.java
deleted file mode 100644
index 29f79081..00000000
--- a/src/ac/ed/lurg/ModelContext.java
+++ /dev/null
@@ -1,16 +0,0 @@
-package ac.ed.lurg;
-
-import ac.ed.lurg.demand.DemandManager;
-
-public class ModelContext {
-	private DemandManager demandManager;
-
-	public ModelContext(DemandManager demandManager) {
-		super();
-		this.demandManager = demandManager;
-	}
-
-	public DemandManager getDemandManager() {
-		return demandManager;
-	}
-}
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index e411c8cf..bdaa29e8 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -13,6 +13,7 @@ import ac.ed.lurg.country.CountryAgent;
 import ac.ed.lurg.country.CountryBoundaryItem;
 import ac.ed.lurg.country.CountryBoundaryReader;
 import ac.ed.lurg.demand.DemandManager;
+import ac.ed.lurg.landuse.CommodityData;
 import ac.ed.lurg.landuse.LandCoverItem;
 import ac.ed.lurg.landuse.LandCoverReader;
 import ac.ed.lurg.types.CropType;
@@ -28,8 +29,8 @@ import ac.sac.raster.RasterSet;
 public class ModelMain {
 
 	private Collection<CountryAgent> countryAgents;
-	private ModelContext modelContext;
 	private Map<Country, List<RasterKey>> countryToKeysMap;
+	private DemandManager demandManager;
 
 	public static void main(String[] args)  {
 		ModelMain theModel = new ModelMain();
@@ -39,10 +40,10 @@ public class ModelMain {
 
 	/* setup models, reading inputs, etc. */
 	private void setup() {
-		DemandManager demandManager = new DemandManager(ModelFitType.LOGISTIC, "SSP2_v9_130325");
-		modelContext = new ModelContext(demandManager);
-		countryToKeysMap = getCountryToKeysMap();
-		countryAgents = getAgents();
+		countryToKeysMap = createCountryToKeysMap();
+		demandManager = new DemandManager(ModelFitType.LOGISTIC, ModelConfig.SSP_SCENARIO);
+
+		countryAgents = createCountryAgents();
 	}
 
 	/* run the model */
@@ -60,50 +61,63 @@ public class ModelMain {
 			Collection<RasterKey> countryKeys = countryToKeysMap.get(ca.getCountry());
 
 			YieldRaster countryYieldSurfaces = yieldSurfaces.getSubsetRasterForKeys(countryKeys);
-			ca.determineProduction(timestep, countryYieldSurfaces);
+			ca.determineProduction(timestep, countryYieldSurfaces, getWorldInputEnergy());
 		}
 
 		// examine global trade balance
 	}
 
-	public Map<Country, List<RasterKey>>  getCountryToKeysMap() {
+	public Map<CropType, Double> getWorldInputEnergy() {
+		Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
+		dummyMap.put(CropType.CEREALS, 9.0);
+		dummyMap.put(CropType.FRUIT, 5.0);
+		dummyMap.put(CropType.OILCROPS, 10.0);
+		dummyMap.put(CropType.STARCHY_ROOTS, 5.0);
+		dummyMap.put(CropType.PULSES, 20.0);
+		dummyMap.put(CropType.VEGETABLES, 20.0);
+		dummyMap.put(CropType.MEAT_OR_PASTURE, 80.0);
+		dummyMap.put(CropType.TREENUTS, 5.0);
+		return dummyMap;
+	}
+
+	public Map<Country, List<RasterKey>> createCountryToKeysMap() {
 		CountryBoundaryReader countryReader = new CountryBoundaryReader();
 		RasterSet<CountryBoundaryItem> countryBoundaries = countryReader.getRasterDataFromFile(ModelConfig.COUNTRY_BOUNDARY_FILE);
 
-		Map<Country, List<RasterKey>> countryToKeysMap = new HashMap<Country, List<RasterKey>>();
+		Map<Country, List<RasterKey>> countryMap = new HashMap<Country, List<RasterKey>>();
 
 		for (Map.Entry<RasterKey, CountryBoundaryItem> entry : countryBoundaries.entrySet()) {
-			List<RasterKey> keys = countryToKeysMap.get(entry.getValue().getCountry());
+			List<RasterKey> keys = countryMap.get(entry.getValue().getCountry());
 
 			if (keys == null) {
 				keys = new ArrayList<RasterKey>();
-				countryToKeysMap.put(entry.getValue().getCountry(), keys);
+				countryMap.put(entry.getValue().getCountry(), keys);
 			}
 			keys.add(entry.getKey());
 		}
-		return countryToKeysMap;
+		return countryMap;
 	}
 
+	public Collection<CountryAgent> createCountryAgents() {
 
-	public Collection<CountryAgent> getAgents() {
 		Collection<CountryAgent> countryAgents = new HashSet<CountryAgent>();
 
 		RasterSet<LandCoverItem> initLC = getInitialLandCover();
+		Map<Country, Map<CropType, CommodityData>> commodityDataMap = CommodityData.readCommodityData();
 
 		for (Map.Entry<Country, List<RasterKey>> entry : countryToKeysMap.entrySet()) {
 			Country country = entry.getKey();
 			List<RasterKey> keys = entry.getValue();
 			RasterSet<LandCoverItem> initCountryLC = initLC.popSubsetForKeys(new RasterSet<LandCoverItem>(), keys);
-
-			CountryAgent ca = new CountryAgent(modelContext, country, initCountryLC);
+			Map<CropType, CommodityData> countryCommodityData = commodityDataMap.get(country);
+					
+			CountryAgent ca = new CountryAgent(demandManager, country, initCountryLC, countryCommodityData);
 			countryAgents.add(ca);
 		}
 
 		return countryAgents;
 	}
 
-
-
 	private RasterSet<LandCoverItem> getInitialLandCover() {
 		String rootDir = ModelConfig.INITAL_LAND_COVER_DIR + File.separator;
 		RasterSet<LandCoverItem> initLC = null;
@@ -116,7 +130,6 @@ public class ModelMain {
 		return initLC;
 	}
 
-
 	private YieldRaster getYieldSurfaces(int timestep) {
 		String rootDir = ModelConfig.YIELD_DIR + File.separator + (timestep + ModelConfig.BASE_YEAR) + File.separator;
 		YieldRaster yieldSurfaces = null;
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 20cdb972..87cdf5a5 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -4,12 +4,13 @@ import java.util.HashMap;
 import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.ModelContext;
 import ac.ed.lurg.country.gams.GamsCountryInput;
 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.DemandManager;
 import ac.ed.lurg.landuse.AreasItem;
+import ac.ed.lurg.landuse.CommodityData;
 import ac.ed.lurg.landuse.LandCoverItem;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandDataType;
@@ -20,7 +21,7 @@ import ac.sac.raster.RasterSet;
 
 public class CountryAgent {
 	
-	private ModelContext modelContext;
+	private DemandManager demandManager;
 	private Country country;
 	
 	private Map<Integer, GamsRasterOutput> resultsTimeseries = new HashMap<Integer, GamsRasterOutput>();
@@ -29,36 +30,35 @@ public class CountryAgent {
 	private Map<CropType, Double> currentProjectedDemand;
 	private YieldRaster countryYieldSurfaces;
 		
-	public CountryAgent(ModelContext modelContext, Country country, RasterSet<LandCoverItem> initialLC) {
-		this.modelContext = modelContext;
+	public CountryAgent(DemandManager demandManager, Country country, RasterSet<LandCoverItem> initialLC, Map<CropType, CommodityData> commoditiesData) {
+		this.demandManager = demandManager;
 		this.country = country;
-		
 		RasterSet<AreasItem> cropAreaRaster = convertInitialLC(initialLC);
 		
+		GamsRasterOutput initialData = new GamsRasterOutput(cropAreaRaster, commoditiesData);
 		
-		Map<CropType, Double> feedAmounts = null;
-		Map<CropType, Double> netImports = null;
-		
-		GamsRasterOutput initialData = new GamsRasterOutput(cropAreaRaster, feedAmounts, netImports);
-		
-		resultsTimeseries.put(-1, initialData);
+		resultsTimeseries.put(0, initialData);
 	}
 
 	private RasterSet<AreasItem> convertInitialLC(RasterSet<LandCoverItem> initialLC) {
 		RasterSet<AreasItem> cropAreaRaster = new RasterSet<AreasItem>();
 		
+		
 		for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
 			AreasItem areasItem = new AreasItem();
+			RasterKey key = entry.getKey();
 			LandCoverItem landCover = entry.getValue();
 			
-			areasItem.setCropArea(CropType.CEREALS, landCover.getLandCover(LandDataType.CROPLAND));  // we allow free substitution between crop types so this is ok-ish
-			areasItem.setCropArea(CropType.MEAT_OR_PASTURE, landCover.getLandCover(LandDataType.PASTURE));
+			double cellArea = key.getHalfDegreeArea();
+
+			areasItem.setCropArea(CropType.CEREALS, landCover.getLandCover(LandDataType.CROPLAND) * cellArea);  // we allow free substitution between crop types so this is ok-ish
+			areasItem.setCropArea(CropType.MEAT_OR_PASTURE, landCover.getLandCover(LandDataType.PASTURE) * cellArea);
 			
-			areasItem.setLandCoverArea(LandDataType.LAND, landCover.getLandCover(LandDataType.LAND));
-			areasItem.setLandCoverArea(LandDataType.FOREST, landCover.getLandCover(LandDataType.FOREST));
-			areasItem.setLandCoverArea(LandDataType.OTHER_NATURAL, landCover.getLandCover(LandDataType.OTHER_NATURAL));
+			areasItem.setLandCoverArea(LandDataType.LAND, landCover.getLandCover(LandDataType.LAND) * cellArea);
+			areasItem.setLandCoverArea(LandDataType.FOREST, landCover.getLandCover(LandDataType.FOREST) * cellArea);
+			areasItem.setLandCoverArea(LandDataType.OTHER_NATURAL, landCover.getLandCover(LandDataType.OTHER_NATURAL) * cellArea);
 
-			cropAreaRaster.put(entry.getKey(), areasItem);
+			cropAreaRaster.put(key, areasItem);
 		}
 		return cropAreaRaster;
 	}
@@ -67,16 +67,16 @@ public class CountryAgent {
 		return country;
 	}
 
-	public void determineProduction(int timestep, YieldRaster countryYieldSurfaces) {
+	public void determineProduction(int timestep, YieldRaster countryYieldSurfaces, Map<CropType, Double> worldInputEnergy) {
 		currentTimestep = timestep;
 		this.countryYieldSurfaces = countryYieldSurfaces;
 		
 		// get projected demand
 		int year = ModelConfig.BASE_YEAR + currentTimestep;
-		currentProjectedDemand = modelContext.getDemandManager().getDemand(country, year);
+		currentProjectedDemand = demandManager.getDemand(country, year);
 		
 		// optimize areas and intensity 
-		GamsRasterInput input = getGamsRasterInput();
+		GamsRasterInput input = getGamsRasterInput(worldInputEnergy);
 		GamsRasterOptimiser opti = new GamsRasterOptimiser(input);
 		LogWriter.println("Running " + country.getCountryName() + ", currentTimestep " + currentTimestep);
 		
@@ -88,53 +88,29 @@ public class CountryAgent {
 		return currentProjectedDemand;
 	}
 
-	public GamsRasterInput getGamsRasterInput() {
+	public GamsRasterInput getGamsRasterInput(Map<CropType, Double> worldInputEnergy) {
 
-		int previousTimestep = currentTimestep-1;
+		int previousTimestep = (ModelConfig.START_TIMESTEP == currentTimestep) ? currentTimestep : currentTimestep - 1;
 		GamsRasterOutput prevOutput = resultsTimeseries.get(previousTimestep);
-				
-		GamsCountryInput countryLevelInputs = new GamsCountryInput(getProjectedDemand(), getWorldInputEnergy(), getMaxNetImport(), getMinNetImport());		
+		
+		Map<CropType, CommodityData> prevCommoditiesData = prevOutput.getCommoditiesData();
+		
+		double allowedImportChange;
+		if (currentTimestep == ModelConfig.START_TIMESTEP)
+			allowedImportChange = 0;
+		else
+			allowedImportChange = ModelConfig.MAX_IMPORT_CHANGE;
+		
+		Map<CropType, Double> maxNetImport = new HashMap<CropType, Double>();
+		Map<CropType, Double> minNetImport = new HashMap<CropType, Double>();
+		for (Map.Entry<CropType, CommodityData> entry : prevCommoditiesData.entrySet()) {
+			maxNetImport.put(entry.getKey(), entry.getValue().getNetImports() * (1 + allowedImportChange));
+			minNetImport.put(entry.getKey(), entry.getValue().getNetImports() * (1 - allowedImportChange));
+		}
+		
+		GamsCountryInput countryLevelInputs = new GamsCountryInput(getProjectedDemand(), worldInputEnergy, maxNetImport, maxNetImport);		
 		GamsRasterInput input = new GamsRasterInput(countryYieldSurfaces, prevOutput.getCropAreaRaster(), countryLevelInputs);
 
 		return input;
 	}
-
-	public Map<CropType, Double> getWorldInputEnergy() {
-		Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
-		dummyMap.put(CropType.CEREALS, 9.0);
-		dummyMap.put(CropType.FRUIT, 5.0);
-		dummyMap.put(CropType.OILCROPS, 10.0);
-		dummyMap.put(CropType.STARCHY_ROOTS, 5.0);
-		dummyMap.put(CropType.PULSES, 20.0);
-		dummyMap.put(CropType.VEGETABLES, 20.0);
-		dummyMap.put(CropType.MEAT_OR_PASTURE, 80.0);
-		dummyMap.put(CropType.TREENUTS, 5.0);
-		return dummyMap;
-	}
-
-	public Map<CropType, Double> getMaxNetImport() {
-		Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>();
-		dummyMaxExports.put(CropType.CEREALS, 24.0);
-		dummyMaxExports.put(CropType.FRUIT, 5.0);
-		dummyMaxExports.put(CropType.OILCROPS, 5.0);
-		dummyMaxExports.put(CropType.STARCHY_ROOTS, 5.0);
-		dummyMaxExports.put(CropType.PULSES, 5.0);
-		dummyMaxExports.put(CropType.VEGETABLES, 5.0);
-		dummyMaxExports.put(CropType.MEAT_OR_PASTURE, 5.0);
-		dummyMaxExports.put(CropType.TREENUTS, 5.0);
-		return dummyMaxExports;
-	}
-
-	public Map<CropType, Double> getMinNetImport() {
-		Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>();
-		dummyMaxExports.put(CropType.CEREALS, -24.0);
-		dummyMaxExports.put(CropType.FRUIT, -5.0);
-		dummyMaxExports.put(CropType.OILCROPS, -5.0);
-		dummyMaxExports.put(CropType.STARCHY_ROOTS, -5.0);
-		dummyMaxExports.put(CropType.PULSES, -5.0);
-		dummyMaxExports.put(CropType.VEGETABLES, -5.0);
-		dummyMaxExports.put(CropType.MEAT_OR_PASTURE, -5.0);
-		dummyMaxExports.put(CropType.TREENUTS, -5.0);
-		return dummyMaxExports;
-	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 2c0111d3..6e7ec3c0 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -10,7 +10,7 @@ public class GamsCountryInput {
 	private Map<CropType, Double> worldInputEnergy;
 	private Map<CropType, Double> maxNetImport;
 	private Map<CropType, Double> minNetImport;
-	
+
 	// limits to areas for each location and crop
 	
 /*	private double maxLandUseChange;
@@ -30,26 +30,26 @@ public class GamsCountryInput {
 		this.minNetImport = minNetImport;
 	}
 
-
 	public Map<CropType, Double> getProjectedDemand() {
 		return projectedDemand;
 	}
 
-
 	public Map<CropType, Double> getWorldInputEnergy() {
 		return worldInputEnergy;
 	}
 
-
 	public Map<CropType, Double> getMaxNetImport() {
 		return maxNetImport;
 	}
 
-
 	public Map<CropType, Double> getMinNetImport() {
 		return minNetImport;
 	}
 
+//	public Map<CropType, Double> getMinFeedAmount() {
+//		return minFeedAmount;
+//	}
+
 	public double getMeatEfficiency() {
 		return 1;  // this is already handled by the feed conversion efficiency for each animal product
 	}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 9e806dea..1c5a2cad 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -8,6 +8,7 @@ import java.util.Vector;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.landuse.AreasItem;
+import ac.ed.lurg.landuse.CommodityData;
 import ac.ed.lurg.landuse.IntensitiesItem;
 import ac.ed.lurg.landuse.Intensity;
 import ac.ed.lurg.types.CropType;
@@ -140,8 +141,7 @@ public class GamsLocationOptimiser {
 
 		Map<Integer, IntensitiesItem> intensities = new HashMap<Integer, IntensitiesItem>();
 		Map<Integer, AreasItem> cropAreas = new HashMap<Integer, AreasItem>();
-		Map<CropType, Double> feedAmounts = new HashMap<CropType, Double>();
-		Map<CropType, Double> netImports = new HashMap<CropType, Double>();
+		Map<CropType, CommodityData> commoditiesData = new HashMap<CropType, CommodityData>();
 		CropType prevCropType = null;
 		
 		for (GAMSVariableRecord rec : varAreas) {
@@ -157,8 +157,8 @@ public class GamsLocationOptimiser {
 			if (!cropType.equals(prevCropType)) {
 				feedAmount = varFeedAmount.findRecord(itemName).getLevel();
 				netImport = varNetImports.findRecord(itemName).getLevel();
-				feedAmounts.put(cropType, feedAmount);
-				netImports.put(cropType, netImport);
+				
+				commoditiesData.put(cropType, new CommodityData(feedAmount, netImport));
 				LogWriter.println(String.format("\n%s:\tfeedAmount= %.1f,\tnetImports= %.3f", itemName, feedAmount, netImport)); 
 			}
 
@@ -186,7 +186,7 @@ public class GamsLocationOptimiser {
 		LogWriter.println(String.format("\nTotal area= %.1f", totalArea));
 		//cleanup(ws.workingDirectory());
 		
-		GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, feedAmounts, netImports);
+		GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, commoditiesData);
 		return results ;
 	}
 
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 2fb834a6..bbb55eb0 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -3,6 +3,7 @@ package ac.ed.lurg.country.gams;
 import java.util.Map;
 
 import ac.ed.lurg.landuse.AreasItem;
+import ac.ed.lurg.landuse.CommodityData;
 import ac.ed.lurg.landuse.IntensitiesItem;
 import ac.ed.lurg.types.CropType;
 
@@ -11,20 +12,17 @@ public class GamsLocationOutput {
 	
 	Map<Integer, IntensitiesItem> intensities;  // data mapped from id (not raster)
 	Map<Integer, AreasItem> cropAreas;
-	Map<CropType, Double> feedAmounts;
-	Map<CropType, Double> netImports;
+	private Map<CropType, CommodityData> commoditiesData;
 	
 	public GamsLocationOutput(int status, 
 			Map<Integer, IntensitiesItem> intensities, 
 			Map<Integer, AreasItem> cropAreas, 
-			Map<CropType, Double> feedAmounts,
-			Map<CropType, Double> netImports) {
+			Map<CropType, CommodityData> commoditiesData) {
 		super();
 		this.status = status;
 		this.intensities = intensities;
 		this.cropAreas = cropAreas;
-		this.feedAmounts = feedAmounts;
-		this.netImports = netImports;
+		this.commoditiesData = commoditiesData;
 	}
 	
 	public int getStatus() {
@@ -36,10 +34,7 @@ public class GamsLocationOutput {
 	public Map<Integer, AreasItem> getCropAreas() {
 		return cropAreas;
 	}
-	public Map<CropType, Double> getFeedAmounts() {
-		return feedAmounts;
-	}
-	public Map<CropType, Double> getNetImports() {
-		return netImports;
+	public Map<CropType, CommodityData> getCommoditiesData() {
+		return commoditiesData;
 	}
 }
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 8afa96fb..15184c50 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -48,7 +48,8 @@ public class GamsRasterOptimiser {
 		RasterSet<AreasItem> newAreaRaster = allocAreas(gamsInput.getPreviousAreas(), gamsOutput);
 		RasterSet<IntensitiesItem> newIntensityRaster = allocIntensities(gamsOutput);
 		
-		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getFeedAmounts(), gamsOutput.getNetImports());
+		
+		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData());
 	}
 
 	private RasterSet<AreasItem> allocAreas(Map<Integer, ? extends AreasItem> prevAreasAgg, GamsLocationOutput gamsOutput) {
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
index 10760cea..0f04f589 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
@@ -3,6 +3,7 @@ package ac.ed.lurg.country.gams;
 import java.util.Map;
 
 import ac.ed.lurg.landuse.AreasItem;
+import ac.ed.lurg.landuse.CommodityData;
 import ac.ed.lurg.landuse.IntensitiesItem;
 import ac.ed.lurg.types.CropType;
 import ac.sac.raster.RasterSet;
@@ -13,19 +14,16 @@ public class GamsRasterOutput {
 	
 	private RasterSet<IntensitiesItem> intensityRaster;
 	private RasterSet<AreasItem> cropAreaRaster;
-	private Map<CropType, Double> feedAmounts;
-	private Map<CropType, Double> netImports;
+	private Map<CropType, CommodityData> commoditiesData;
 	
-	public GamsRasterOutput(RasterSet<AreasItem> cropAreaRaster, Map<CropType, Double> feedAmounts, Map<CropType, Double> netImports) {
+	public GamsRasterOutput(RasterSet<AreasItem> cropAreaRaster, Map<CropType, CommodityData> commoditiesData) {
 		super();
 		this.cropAreaRaster = cropAreaRaster;
-		this.feedAmounts = feedAmounts;
-		this.netImports = netImports;
+		this.commoditiesData = commoditiesData;
 	}
 
-	public GamsRasterOutput(int status, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, 
-			Map<CropType, Double> feedAmounts, Map<CropType, Double> netImports) {
-		this(cropAreaRaster, feedAmounts, netImports);
+	public GamsRasterOutput(int status, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, Map<CropType, CommodityData> commoditiesData) {
+		this(cropAreaRaster, commoditiesData);
 		this.status = status;
 		this.intensityRaster = intensityRaster;
 	}
@@ -39,10 +37,7 @@ public class GamsRasterOutput {
 	public RasterSet<AreasItem> getCropAreaRaster() {
 		return cropAreaRaster;
 	}
-	public Map<CropType, Double> getFeedAmounts() {
-		return feedAmounts;
-	}
-	public Map<CropType, Double> getNetImports() {
-		return netImports;
+	public Map<CropType, CommodityData> getCommoditiesData() {
+		return commoditiesData;
 	}
 }
diff --git a/src/ac/ed/lurg/landuse/CommodityData.java b/src/ac/ed/lurg/landuse/CommodityData.java
new file mode 100644
index 00000000..31b8e9ea
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/CommodityData.java
@@ -0,0 +1,78 @@
+package ac.ed.lurg.landuse;
+
+import java.io.BufferedReader;
+import java.io.FileReader;
+import java.util.HashMap;
+import java.util.Map;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.country.Country;
+import ac.ed.lurg.country.CountryManager;
+import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.LogWriter;
+
+public class CommodityData {
+	private static final int COUNTRY_COL = 0; 
+	private static final int COMMODITY_COL = 2; 
+	private static final int FEED_COL = 4; 
+	private static final int NET_IMPORT_COL = 5; 
+
+	private double feedAmount;
+	private double netImports;
+	
+	public CommodityData(double feedAmount, double netImports) {
+		super();
+		this.feedAmount = feedAmount;
+		this.netImports = netImports;
+	}
+
+	public double getFeedAmount() {
+		return feedAmount;
+	}
+
+	public double getNetImports() {
+		return netImports;
+	}
+	
+	public static Map<Country, Map<CropType, CommodityData>> readCommodityData() {
+		Map<Country, Map<CropType, CommodityData>> commodityMap = new HashMap<Country, Map<CropType, CommodityData>>();
+		String filename = ModelConfig.COMMODITY_DATA_FILE;
+		try {
+			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
+			String line, countryName, commodityName;
+			double feedAmount, netImports;
+			fitReader.readLine(); // read header
+
+			while ((line=fitReader.readLine()) != null) {
+				String[] tokens = line.split(",");
+				
+				if (tokens.length < 6)
+					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
+				
+				countryName = tokens[COUNTRY_COL];
+				commodityName = tokens[COMMODITY_COL];
+				feedAmount = Double.parseDouble(tokens[FEED_COL]);
+				netImports = Double.parseDouble(tokens[NET_IMPORT_COL]);
+
+				Country country = CountryManager.getForName(countryName);
+				
+				CommodityData data = new CommodityData(feedAmount, netImports);
+				Map<CropType, CommodityData> countryData = commodityMap.get(country);
+				if (countryData == null) {
+					countryData = new HashMap<CropType, CommodityData>();
+					commodityMap.put(country, countryData);
+				}
+				
+				countryData.put(CropType.getForFaoName(commodityName), data);
+			} 
+			fitReader.close(); 
+		
+		} catch (Exception e) {
+			LogWriter.printlnError("Failed in reading commodity demand fits");
+			e.printStackTrace();
+		}
+		LogWriter.println("Processed " + filename + ", create " + commodityMap.size() + " country commodity maps values");
+		
+		return commodityMap;
+	}
+}
diff --git a/src/ac/sac/raster/RasterKey.java b/src/ac/sac/raster/RasterKey.java
index 2fc3a4e7..4b0def60 100755
--- a/src/ac/sac/raster/RasterKey.java
+++ b/src/ac/sac/raster/RasterKey.java
@@ -38,6 +38,19 @@ public class RasterKey implements Serializable {
 		return d;
 	}
 	
+	/** Assuming coordinates refer to global half degree in Mha */
+	public double getHalfDegreeArea() {
+		// https://badc.nerc.ac.uk/help/coordinates/cell-surf-area.html
+		double cellArea = Math.PI / 360.0 * (Math.sin(getLatInRadians(row)) - Math.sin(getLatInRadians(row+1))) * 6371.0 * 6371.0 / 10000.0;
+		return cellArea;
+	}
+	
+	private static double getLatInRadians(double rowNum) {
+		double latDegrees = 90.0 - rowNum / 2;
+		double latRadian = latDegrees * Math.PI / 180.0;
+		return latRadian;
+	}
+	
 	@Override
 	public int hashCode() {
 		final int prime = 31;
-- 
GitLab