Change values held by AreasItem

GAMS/IntExtOpt.gms

@@ -215,6 +215,9 @@ $gdxin
  
  totalProdCost('meat') = sum(feed_crop, feedEnergy(feed_crop))
  
+ parameter totalCropland(location);
+ totalCropland(location) = sum(crop_less_pasture, area.l(crop_less_pasture, location));
+ 
  display totalProdCost, totalProd, totalImportEnergy, feedEnergy, totalProdCost
  
          

src/ac/ed/lurg/ModelMain.java

@@ -28,7 +28,7 @@ import ac.ed.lurg.output.RasterOutputer;
 import ac.ed.lurg.types.CropToDoubleMap;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.FertiliserRate;
-import ac.ed.lurg.types.LandDataType;
+import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.ModelFitType;
 import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.yield.LPJYieldResponseMapReader;
@@ -281,7 +281,7 @@ public class ModelMain {
 
 		String rootDir = ModelConfig.INITAL_LAND_COVER_DIR + File.separator;
 
-		for (LandDataType landType : LandDataType.values()) {
+		for (LandCoverType landType : LandCoverType.values()) {
 			LandCoverReader lcReader = new LandCoverReader(initLC, landType);
 			initLC = lcReader.getRasterDataFromFile(rootDir + landType.getFileName());
 		}

src/ac/ed/lurg/country/CountryAgent.java

@@ -15,7 +15,6 @@ import ac.ed.lurg.landuse.IrrigationCostItem;
 import ac.ed.lurg.landuse.LandCoverItem;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandDataType;
 import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.yield.YieldRaster;
 import ac.sac.raster.RasterKey;
@@ -48,21 +47,15 @@ public class CountryAgent {
 	private RasterSet<AreasItem> convertInitialLC(RasterSet<LandCoverItem> initialLC) {
 		RasterSet<AreasItem> cropAreaRaster = new RasterSet<AreasItem>(initialLC.getHeaderDetails());
 		
-		
 		for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
 			AreasItem areasItem = new AreasItem();
 			RasterKey key = entry.getKey();
-			LandCoverItem landCover = entry.getValue();
-			
-			if (landCover != null) {
-				areasItem.setCropArea(CropType.WHEAT, landCover.getLandCover(LandDataType.CROPLAND));  // we allow free substitution between crop types so this is ok-ish
-				areasItem.setCropArea(CropType.PASTURE, landCover.getLandCover(LandDataType.PASTURE));
-				areasItem.setLandCoverArea(LandDataType.FOREST, landCover.getLandCover(LandDataType.FOREST));
-				areasItem.setLandCoverArea(LandDataType.OTHER_NATURAL, landCover.getLandCover(LandDataType.OTHER_NATURAL));
-			}
-			
+			areasItem.setLandCoverAreas(entry.getValue());
+			areasItem.setCropFraction(CropType.WHEAT, 0.5); // random start, better if we could get data, but free substitution between crops so not critical
+			areasItem.setCropFraction(CropType.MAIZE, 0.5);
 			cropAreaRaster.put(key, areasItem);
 		}
+		
 		return cropAreaRaster;
 	}
 	

src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java

@@ -14,6 +14,8 @@ import ac.ed.lurg.landuse.Intensity;
 import ac.ed.lurg.landuse.IrrigationCostItem;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.types.LandCoverType;
+import ac.ed.lurg.utils.LazyHashMap;
 import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.yield.YieldResponsesItem;
 
@@ -92,21 +94,22 @@ public class GamsLocationOptimiser {
 			Integer locationId = entry.getKey();
 			String locString = Integer.toString(locationId);
 			AreasItem areasItem = entry.getValue();
+			
+			double suitableLand = areasItem.getSuitableLand();
+			if (DEBUG) LogWriter.println(String.format("  %d   %15s,\t %.1f", locationId, "suitableLand", suitableLand));
+			landP.addRecord(locString).setValue(suitableLand);
+			
+			double pastureArea = areasItem.getLandCoverArea(LandCoverType.PASTURE);
+			if (DEBUG) LogWriter.println(String.format("  %d   %15s,\t %.1f", locationId, "pastureArea", pastureArea));
+			setPreviousArea(prevCropP, locString, CropType.PASTURE.getGamsName(), pastureArea);
+			totalAgriLand += pastureArea;
 
-			for (CropType cropType : CropType.values()) {
+			for (CropType cropType : CropType.getCropsLessPasture()) {
 				double d = areasItem.getCropArea(cropType);
 				if (DEBUG && d != 0) LogWriter.println(String.format("  %d   %15s,\t %.1f", locationId, cropType.getGamsName(), d));
-				
-				Vector<String> v = new Vector<String>();
-				v.add(cropType.getGamsName());
-				v.add(locString);
-				prevCropP.addRecord(v).setValue(d);
+				setPreviousArea(prevCropP, locString, cropType.getGamsName(), d);
 				totalAgriLand += d;
 			}
-
-			double suitableLand = areasItem.getSuitableLand();
-			if (DEBUG) LogWriter.println(String.format("  %d   %15s,\t %.1f", locationId, "suitableLand", suitableLand));
-			landP.addRecord(locString).setValue(suitableLand);
 		}
 		if (DEBUG) LogWriter.println(String.format("  Total agricultural\t %.1f", totalAgriLand));
 
@@ -183,6 +186,14 @@ public class GamsLocationOptimiser {
 		addItemMapParm(inDB.addParameter("cropAdj", 1), countryInput.getCropAdjustments());
 	}
 
+	private void setPreviousArea(GAMSParameter prevCropP, String locString, String cropTypeString, double d) {
+		Vector<String> v = new Vector<String>();
+		v.add(cropTypeString);
+		v.add(locString);
+		prevCropP.addRecord(v).setValue(d);
+	}
+
+	@SuppressWarnings("serial")
 	private GamsLocationOutput handleResults(GAMSDatabase outDB) {
 		int modelStatus = (int) outDB.getParameter("ms").findRecord().getValue();
 		System.out.println(
@@ -200,15 +211,20 @@ public class GamsLocationOptimiser {
 		GAMSParameter parmCropAdj = outDB.getParameter("cropAdj");
 		GAMSParameter parmProd = outDB.getParameter("totalProd");
 		GAMSParameter parmProdCost = outDB.getParameter("totalProdCost");
+		GAMSParameter parmCroplandArea = outDB.getParameter("totalCropland");
 
 		double totalArea = 0;
 		double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, feedAmount, netImport, yield, unitEnergy, cropAdj, prod, prodCost;
 
-		Map<Integer, IntensitiesItem> intensities = new HashMap<Integer, IntensitiesItem>();
-		Map<Integer, AreasItem> cropAreas = new HashMap<Integer, AreasItem>();
+		LazyHashMap<Integer, IntensitiesItem> intensities = new LazyHashMap<Integer, IntensitiesItem>() {
+			protected IntensitiesItem createValue() { return new IntensitiesItem(); }
+		};
+		LazyHashMap<Integer, AreasItem> cropAreas = new LazyHashMap<Integer, AreasItem>() { 
+			protected AreasItem createValue() { return new AreasItem(); }
+		};
+		
 		Map<CropType, CropUsageData> cropUsageData = new HashMap<CropType, CropUsageData>();
 		Map<CropType, Double> cropAdjs = new HashMap<CropType, Double>();
-		CropType prevCropType = null;
 
 		for (GAMSVariableRecord rec : varAreas) {
 			String itemName = rec.getKeys()[0];
@@ -223,39 +239,33 @@ public class GamsLocationOptimiser {
 			int locId = Integer.parseInt(locationName);
 			CropType cropType = CropType.getForGamsName(itemName);
 
-			if (!cropType.equals(prevCropType)) {
+			if (!cropUsageData.containsKey(cropType)) {  // then we must not have seen this crop type before, so need to do all non location specific stuff
 				feedAmount = varFeedAmount.findRecord(itemName).getLevel();
 				netImport = cropType.isImportedCrop() ? varNetImports.findRecord(itemName).getLevel() : 0;
 				cropAdj = getParmValue(parmCropAdj, itemName);
 				prod =  getParmValue(parmProd, itemName);
 				prodCost = getParmValue(parmProdCost, itemName);
-
+				
 				cropUsageData.put(cropType, new CropUsageData(feedAmount, netImport, prod, prodCost));
 				cropAdjs.put(cropType, cropAdj);
 				if (DEBUG) LogWriter.println(String.format("\n%s:\tfeedAmount= %.1f,\tnetImports= %.3f,\tprod= %.3f,\tprodCost= %.3f,\tcropAdj= %.3f", itemName, feedAmount, netImport, prod, prodCost, cropAdj)); 
 			}
 
 			if (area > 0) { 
-				if (DEBUG) LogWriter.println(String.format("\t location %s, %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f", 
-						locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity));
-				
-				IntensitiesItem intensityItem = intensities.get(locId);
-				if (intensityItem == null) {
-					intensityItem = new IntensitiesItem();
-					intensities.put(locId, intensityItem);
-				}
+				if (DEBUG) LogWriter.println(String.format("\t location %s, %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f", locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity));
+				IntensitiesItem intensityItem = intensities.lazyGet(locId);
 				intensityItem.setIntensity(cropType, new Intensity(fertIntensity, irrigIntensity, otherIntensity, yield, unitEnergy));
 			}
 
-			AreasItem areaItem = cropAreas.get(locId);
-			if (areaItem == null) {
-				areaItem = new AreasItem();
-				cropAreas.put(locId, areaItem);
-			}
-			areaItem.setCropArea(cropType, area);
-			totalArea += area;
+			double croplandArea = getParmValue(parmCroplandArea, locationName);
+			AreasItem areasItem = cropAreas.lazyGet(locId);
+			if (cropType.equals(CropType.PASTURE))
+				areasItem.setLandCoverArea(LandCoverType.PASTURE, area);
+			else
+				areasItem.setLandCoverArea(LandCoverType.CROPLAND, croplandArea);  // will set this multiple times, once for each arable crop, but doesn't really matter
 
-			prevCropType = cropType;
+			areasItem.setCropFraction(cropType, croplandArea > 0 ? area/croplandArea : 0);
+			totalArea += area;
 		}
 
 		netImport = varNetImports.findRecord(CropType.MEAT.getGamsName()).getLevel();

src/ac/ed/lurg/country/gams/GamsLocationTest.java

@@ -8,7 +8,7 @@ import ac.ed.lurg.landuse.AreasItem;
 import ac.ed.lurg.landuse.IrrigationCostItem;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandDataType;
+import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.YieldType;
 import ac.ed.lurg.yield.YieldResponsesItem;
 
@@ -94,18 +94,19 @@ public class GamsLocationTest {
 	AreasItem getAreaItem(double cereal, double pasture) {
 		AreasItem dummyMap = new AreasItem();
 		
-		dummyMap.setCropArea(CropType.WHEAT, 5.0 + cereal);
-		dummyMap.setCropArea(CropType.MAIZE, 5.0);
-		dummyMap.setCropArea(CropType.RICE, 5.0);
-		dummyMap.setCropArea(CropType.TROPICAL_CEREALS, 5.0);
-		dummyMap.setCropArea(CropType.SOYBEAN, 5.0);
-		dummyMap.setCropArea(CropType.OILCROPS, 5.0);
-		dummyMap.setCropArea(CropType.PULSES, 5.0);
-		dummyMap.setCropArea(CropType.STARCHY_ROOTS, 5.0);
-		dummyMap.setCropArea(CropType.PASTURE, 5.0 + pasture);
+		dummyMap.setCropFraction(CropType.WHEAT, 0.2);
+		dummyMap.setCropFraction(CropType.MAIZE, 0.2);
+		dummyMap.setCropFraction(CropType.RICE, 0.1);
+		dummyMap.setCropFraction(CropType.TROPICAL_CEREALS, 0.1);
+		dummyMap.setCropFraction(CropType.SOYBEAN, 0.1);
+		dummyMap.setCropFraction(CropType.OILCROPS, 0.1);
+		dummyMap.setCropFraction(CropType.PULSES, 0.1);
+		dummyMap.setCropFraction(CropType.STARCHY_ROOTS, 0.1);
 		
-		dummyMap.setLandCoverArea(LandDataType.FOREST, 200 - pasture - cereal);
-		dummyMap.setLandCoverArea(LandDataType.OTHER_NATURAL, 30);
+		dummyMap.setLandCoverArea(LandCoverType.PASTURE, pasture);
+		dummyMap.setLandCoverArea(LandCoverType.CROPLAND, cereal);
+		dummyMap.setLandCoverArea(LandCoverType.FOREST, 200 - pasture - cereal);
+		dummyMap.setLandCoverArea(LandCoverType.OTHER_NATURAL, 30);
 		return dummyMap;
 	}
 

src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java

 package ac.ed.lurg.country.gams;
 
 import java.util.ArrayList;
-import java.util.Arrays;
 import java.util.Collections;
 import java.util.HashSet;
 import java.util.List;
@@ -9,12 +8,13 @@ import java.util.Map;
 import java.util.Map.Entry;
 import java.util.Set;
 
+import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.landuse.AreasItem;
 import ac.ed.lurg.landuse.IntensitiesItem;
 import ac.ed.lurg.landuse.IrrigationCostItem;
-import ac.ed.lurg.types.CropToDouble;
+import ac.ed.lurg.output.RasterOutputer;
 import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandDataType;
+import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.YieldType;
 import ac.ed.lurg.utils.LazyHashMap;
 import ac.ed.lurg.utils.LogWriter;
@@ -61,9 +61,35 @@ public class GamsRasterOptimiser {
 		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData(), gamsOutput.getCropAdjustments(), locationIdRaster);
 	}
 
+	private RasterSet<AreasItem> createWithSameLandCovers(RasterSet<AreasItem> toCopy) {
+		RasterSet<AreasItem> theCopy = new RasterSet<AreasItem>(toCopy.getHeaderDetails());
+		
+		for (Entry<RasterKey, AreasItem> entry : toCopy.entrySet()) {
+			AreasItem newAreasItem = new AreasItem();
+			newAreasItem.setLandCoverAreas(entry.getValue());
+			theCopy.put(entry.getKey(), newAreasItem);
+		}
+		
+		return theCopy;
+	}
+	
+	@SuppressWarnings("unused")
+	private void checkedTotalAreas(RasterSet<AreasItem> areaRaster, String comment) {
+		for (LandCoverType l : LandCoverType.values()) {
+			double total = 0;
+			for (AreasItem a : areaRaster.values()) {
+				total += a.getLandCoverArea(l);
+			}
+			
+			LogWriter.printlnError("Total Area " + comment + ": " + l.getName() + ": " + total);
+		}
+	}
+	
 	private RasterSet<AreasItem> allocAreas(Map<Integer, ? extends AreasItem> prevAreasAgg, GamsLocationOutput gamsOutput) {
-		RasterSet<AreasItem> prevAreaRaster = rasterInputData.getPreviousAreas();
-		RasterSet<AreasItem> newAreaRaster = new RasterSet<AreasItem>(prevAreaRaster.getHeaderDetails());
+		RasterSet<AreasItem> newAreaRaster = createWithSameLandCovers(rasterInputData.getPreviousAreas());
+
+		//checkedTotalAreas(rasterInputData.getPreviousAreas(), "old");
+		//checkedTotalAreas(newAreaRaster, "new");
 
 		for (Map.Entry<Integer, AreasItem> entry : gamsOutput.getCropAreas().entrySet()) {
 			Integer locId = entry.getKey();
@@ -71,11 +97,13 @@ public class GamsRasterOptimiser {
 			AreasItem prevAreaAggItem = prevAreasAgg.get(locId);
 			Set<RasterKey> keys = mapping.get(locId);
 
-			CropToDouble cropChangeTotals = newAreaAggItem.getCropChanges(prevAreaAggItem);
+			//checkedTotalAreas(newAreaRaster.popSubsetForKeys(new RasterSet<AreasItem>(), keys), locId + " setset");
+
+			
 			// if (DEBUG) LogWriter.println("Processing location id " + locId);
 			
-			double pastureChange = cropChangeTotals.get(CropType.PASTURE);
-			double croplandChange = cropChangeTotals.getCroplandTotal();
+			double pastureChange = newAreaAggItem.getLandCoverArea(LandCoverType.PASTURE) - prevAreaAggItem.getLandCoverArea(LandCoverType.PASTURE);
+			double croplandChange = newAreaAggItem.getLandCoverArea(LandCoverType.CROPLAND) - prevAreaAggItem.getLandCoverArea(LandCoverType.CROPLAND);
 			double prevForestToNaturalFraction = prevAreaAggItem.getForestToNaturalFraction();
 
 			double pastureFromCrop = 0;
@@ -84,103 +112,74 @@ public class GamsRasterOptimiser {
 
 			if (pastureChange > 0 && croplandChange < 0) {
 				pastureFromCrop = Math.min(pastureChange, -croplandChange);
-				pastureFromNatural = (pastureChange > -croplandChange) ? 0 : pastureChange+croplandChange;
+				pastureFromNatural = (pastureChange > -croplandChange) ? pastureChange+croplandChange : 0;
+				cropFromNatural = croplandChange + pastureFromCrop;
 			}
 			else if (pastureChange < 0 && croplandChange > 0) {
 				pastureFromCrop = -Math.min(-pastureChange, croplandChange);
-				cropFromNatural = (croplandChange > -pastureChange) ? 0 : -(pastureChange+croplandChange);
+				cropFromNatural = (croplandChange > -pastureChange) ? -(pastureChange+croplandChange) : 0;
+				pastureFromNatural = pastureChange - pastureFromCrop;
 			}
 			else {
 				pastureFromNatural = pastureChange;
 				cropFromNatural = croplandChange;
 			}
 			
-	//		allocAllAreasForAgg(newAreaRaster, prevAreaRaster, keys, pastureFromNatural);
+			double shortfall = 0;
+			shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.PASTURE, LandCoverType.FOREST, pastureFromNatural * prevForestToNaturalFraction);
+			shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.PASTURE, LandCoverType.OTHER_NATURAL, pastureFromNatural * (1-prevForestToNaturalFraction));
+			shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.PASTURE, LandCoverType.CROPLAND, pastureFromCrop);
+			shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.CROPLAND, LandCoverType.FOREST, cropFromNatural * prevForestToNaturalFraction);
+			shortfall += allocAllLandCrops(newAreaRaster, keys, LandCoverType.CROPLAND, LandCoverType.OTHER_NATURAL, cropFromNatural * (1-prevForestToNaturalFraction));
 
-			
-			double shortfall = allocAllAreasForAgg(newAreaRaster, prevAreaRaster, keys, cropChangeTotals);
 			if (shortfall > 0.00001)
 				LogWriter.printlnError("This should never happen, due to GAMS constraint. Not able to incorporate all changes, as not enough forest or natural areas left: " + shortfall);
+
+			for (RasterKey key : keys) {
+				AreasItem newAreasRasterItem = newAreaRaster.get(key);
+				for (CropType crop : CropType.values()) {
+					newAreasRasterItem.setCropFraction(crop, newAreaAggItem.getCropFraction(crop));
+				}
+			}
 		}
 
 		return newAreaRaster;
 	}
-
-	private double allocAllAreasForAgg(RasterSet<AreasItem> newAreaRaster, RasterSet<AreasItem> prevAreaRaster, Set<RasterKey> keys, CropToDouble cropChange) {
-
+	
+	private double allocAllLandCrops(RasterSet<AreasItem> newAreaRaster, Set<RasterKey> keys, LandCoverType toLC, LandCoverType fromLC, double change) {
+		LandCoverType toType = toLC;
+		LandCoverType fromType = fromLC;
+		
+		// reverse direction if negative
+		if (change < 0) {
+			change = -change;
+			toType = fromLC;
+			fromType = toLC;
+		}
+		
 		double totalShortfall = 0;
 		Set<RasterKey> keysWithSpace = new HashSet<RasterKey>();
-		CropToDouble avgCropChange = cropChange.scale(1.0/keys.size());
+		double avgChange = change/keys.size();
 
 		for (RasterKey key : keys) {
 			//if (DEBUG) LogWriter.println("  Processing raster key " + key);
-			AreasItem prevAreasRasterItem = prevAreaRaster.get(key);
-
-			AreasItem newAreasRasterItem = new AreasItem();
-			double shortfall = allocRasterArea(newAreasRasterItem, prevAreasRasterItem, avgCropChange);
+			AreasItem newAreasRasterItem = newAreaRaster.get(key);
+			
+			double shortfall = newAreasRasterItem.moveAreas(toType, fromType, avgChange);
 			if (shortfall == 0)
 				keysWithSpace.add(key);
 			else
 				totalShortfall += shortfall;
-
-			newAreaRaster.put(key, newAreasRasterItem);
 		}
 
 		if (totalShortfall > 0 & keysWithSpace.size() > 0) {  // more to allocate and some free areas to allocate into
-			if (DEBUG) LogWriter.println("A total area of " + totalShortfall + " is unallocated, so trying again with areas with space");
-			double factor = totalShortfall / cropChange.getTotal();
-
-			allocAllAreasForAgg(newAreaRaster, newAreaRaster, keysWithSpace, cropChange.scale(factor));
+			if (DEBUG) LogWriter.println("A total land cover of " + totalShortfall + " is unallocated, so trying again with areas with space");
+			totalShortfall = allocAllLandCrops(newAreaRaster, keysWithSpace, toType, fromType, totalShortfall);
 		}
 
 		return totalShortfall;
 	}
 
-	/** returns the net area of crop change that not been able to be allocated */
-	private double allocRasterArea(AreasItem newAreasRasterItem, AreasItem prevAreasRasterItem, CropToDouble cropChange) {
-
-		double netAddedCrop = cropChange.getTotal();
-		double forestRemoved = netAddedCrop/2;
-		double naturalRemoved = netAddedCrop/2;
-		double prevForest = prevAreasRasterItem.getLandCoverArea(LandDataType.FOREST);
-		double prevNatural = prevAreasRasterItem.getLandCoverArea(LandDataType.OTHER_NATURAL);
-
-		double ratioOfChanges = 1;
-
-		if (netAddedCrop > prevForest + prevNatural) {
-			ratioOfChanges = (prevForest + prevNatural) / netAddedCrop;
-			forestRemoved = prevForest;
-			naturalRemoved = prevNatural;
-			if (DEBUG) LogWriter.println("Not able to incorporate all changes for raster location: " + netAddedCrop);
-		}
-		else if (prevForest < forestRemoved) {
-			forestRemoved = prevForest;
-			naturalRemoved = netAddedCrop - forestRemoved;
-		}
-		else if (prevNatural < naturalRemoved) {
-			naturalRemoved = prevNatural;
-			forestRemoved = netAddedCrop - naturalRemoved;
-		}
-		else {
-			// enough space to keep 50/50 allocation
-		}
-
-		newAreasRasterItem.setLandCoverArea(LandDataType.FOREST, prevForest - forestRemoved);
-		newAreasRasterItem.setLandCoverArea(LandDataType.OTHER_NATURAL, prevNatural - naturalRemoved);
-
-		for (CropType crop : CropType.values()) {
-			double prevCropArea = prevAreasRasterItem.getCropArea(crop);
-			Double d = cropChange.get(crop);
-
-			double additionalArea = (d == null) ? 0 : d * ratioOfChanges;	
-			//	LogWriter.println(String.format("%s, %.1f", crop, additionalArea));
-
-			newAreasRasterItem.setCropArea(crop, prevCropArea + additionalArea);
-		}
-
-		return netAddedCrop - prevForest - prevNatural;
-	}
-
 	private RasterSet<IntensitiesItem> allocIntensities(GamsLocationOutput gamsOutput) {
 
 		RasterSet<IntensitiesItem> newIntensityRaster = new RasterSet<IntensitiesItem>(rasterInputData.getPreviousAreas().getHeaderDetails());
@@ -201,6 +200,35 @@ public class GamsRasterOptimiser {
 		// as a first attempt only going to look at pasture and cereal yields, assuming other yields will be correlated to one or the other.
 		YieldRaster yieldRaster = rasterInputData.getYields();
 		RasterSet<AreasItem> cropAreaRaster = rasterInputData.getPreviousAreas();
+		
+	/*	checkedTotalAreas(cropAreaRaster, "before");
+		new RasterOutputer<AreasItem>(cropAreaRaster, "BeforeLandCover") {
+			@Override
+			public Double getValue(RasterKey location) {
+				AreasItem area = results.get(location);
+				if (area == null)
+					return null;
+
+				int count = 0;
+				double maxLC = 0;
+				int maxId = 0;
+				for (LandCoverType l : LandCoverType.values()) {
+					count++;
+					if (maxLC < area.getLandCoverArea(l)) {
+						maxLC = area.getLandCoverArea(l);
+						maxId = count;
+					}
+				}
+				return (double)maxId;
+			}
+
+			@Override
+			public int convertToPixelValue(double value) {
+				return (int) (1 + value * 2);
+			}
+		}.writeOutput(true); */
+
+
 		RasterSet<IrrigationCostItem> irrigCostRaster = rasterInputData.getIrrigationCost();
 
 		{
@@ -247,6 +275,8 @@ public class GamsRasterOptimiser {
 
 		int countFound = 0, countMissing = 0;
 
+		
+		
 		for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
 			YieldResponsesItem yresp = entry.getValue();
 			RasterKey key = entry.getKey();
@@ -264,7 +294,7 @@ public class GamsRasterOptimiser {
 					LogWriter.printlnError("GamsRasterOptimiser: Can't find cropAreas for " + key);
 					continue;
 				}
-
+				
 				IrrigationCostItem irrigCost = irrigCostRaster.get(key);
 
 				int cerealCat = findCategory(wheatlDivisions, yresp.getYieldMax(CropType.WHEAT));
@@ -274,7 +304,7 @@ public class GamsRasterOptimiser {
 				AveragingYieldResponsesItem avgYResp = aggregatedYields.lazyGet(id);
 				AreasItem aggAreas = aggregatedAreas.lazyGet(id);
 				IrrigationCostItem avgIrigCost = aggregatedIrrigCosts.lazyGet(id);
-				mapping.lazyGet(id).add(entry.getKey()); 
+				mapping.lazyGet(id).add(key); 
 
 				// Irrigation cost
 				double suitableAreaThisTime  = cropAreas.getSuitableLand();
@@ -284,7 +314,7 @@ public class GamsRasterOptimiser {
 					//LogWriter.println(id + ", " + key + ", " + avgIrigCost.getIrrigCost() + ", from " + areaSoFar + ", " + suitableAreaThisTime + ", " + irrigCost.getIrrigCost());
 				}
 
-				// Crop yields
+				// Crops yields and area fractions
 				for (CropType crop : CropType.getNonMeatTypes()) {
 					for (YieldType yieldType : YieldType.values()) {
 						avgYResp.setYield(yieldType, crop, yresp.getYield(yieldType, crop));
@@ -292,11 +322,14 @@ public class GamsRasterOptimiser {
 
 					double areaSoFar = aggAreas.getCropArea(crop);
 					double areaThisTime = cropAreas.getCropArea(crop);
-					aggAreas.setCropArea(crop, areaSoFar + areaThisTime);
+					double updateCropland = (cropAreas.getLandCoverArea(LandCoverType.CROPLAND) + aggAreas.getLandCoverArea(LandCoverType.CROPLAND));
+					
+					if (updateCropland > 0)
+						aggAreas.setCropFraction(crop, (areaSoFar + areaThisTime) / updateCropland);
 				}
-
-				// Land cover, we only need forest and other natural
-				for (LandDataType landType : Arrays.asList(LandDataType.FOREST, LandDataType.OTHER_NATURAL)) {
+				
+				// Land covers ares
+				for (LandCoverType landType : LandCoverType.values()) {
 					double areaThisTime = cropAreas.getLandCoverArea(landType);
 					double areaSoFar = aggAreas.getLandCoverArea(landType);
 					aggAreas.setLandCoverArea(landType, areaSoFar + areaThisTime);
@@ -306,7 +339,6 @@ public class GamsRasterOptimiser {
 
 		LogWriter.println("YieldResponsesItem: " + rasterInputData.getCountryInput().getCountry() + ", countFound=" + countFound + ", countMissing=" + countMissing);
 
-
 		for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) {
 			LogWriter.println(e.getKey() + " category has " + e.getValue().size() + " raster areas");
 			if (DEBUG) for (RasterKey key : e.getValue()) {

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

@@ -6,7 +6,7 @@ import java.util.Map.Entry;
 
 import ac.ed.lurg.types.CropToDouble;
 import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandDataType;
+import ac.ed.lurg.types.LandCoverType;
 import ac.sac.raster.RasterItem;
 
 /*
@@ -14,57 +14,84 @@ import ac.sac.raster.RasterItem;
  */
 public class AreasItem implements RasterItem {
 
-	private Map<CropType, Double> cropAreas = new HashMap<CropType, Double>();
-	private Map<LandDataType, Double> landCoverAreas = new HashMap<LandDataType, Double>();
+	private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
+	private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>();
+	
+	public void setLandCoverAreas(LandCoverItem landCover) {
+		if (landCover != null) {
+			for (LandCoverType lcType : LandCoverType.values())
+				landCoverAreas.put(lcType, landCover.getLandCover(lcType));
+		}
+	}
+	
+	public void setLandCoverAreas(AreasItem otherAreasItem) {
+		for (Entry<LandCoverType, Double> entry : otherAreasItem.landCoverAreas.entrySet()) {
+			landCoverAreas.put(entry.getKey(), entry.getValue());
+		}
+	}
+
+	/** move areas from one land cover to another, return an residual not possible */
+	public double moveAreas(LandCoverType toType, LandCoverType fromType, double changeReq) {
+		
+		double prevTo = getLandCoverArea(toType);
+		double prevFrom = getLandCoverArea(fromType);		
+		double actuallyChanged = Math.min(prevFrom, changeReq);
+		
+		setLandCoverArea(toType, prevTo + actuallyChanged);
+		setLandCoverArea(fromType, prevFrom - actuallyChanged);
+
+		return changeReq - actuallyChanged;
+	}
+
 
 	public double getCropArea(CropType c) {
-		Double d = cropAreas.get(c);
+		Double d = getLandCoverArea(LandCoverType.CROPLAND) * getCropFraction(c);
+		return d;
+	}
+	
+	public double getCropFraction(CropType c) {
+		Double d = cropFractions.get(c);
 		return d == null ? 0.0 : d;
 	}
+
 	
-	public void setCropArea(CropType c, double area) {
-		cropAreas.put(c, area);
+	public void setCropFraction(CropType c, double area) {
+		cropFractions.put(c, area);
 	}
 	
-	public double getLandCoverArea(LandDataType c) {
+	public double getLandCoverArea(LandCoverType c) {
 		Double d = landCoverAreas.get(c);
 		return d == null ? 0.0 : d;
 	}
 	
-	public void setLandCoverArea(LandDataType c, double d) {
+	public void setLandCoverArea(LandCoverType c, double d) {
 		landCoverAreas.put(c, d);
 	}
-	
-	public double getTotalCropIncPastureArea() {
-		double total = 0;
-		for (double d : cropAreas.values()) {
-			total += d;
-		}
-		return total;
-	}
-	
+		
 	public double getSuitableLand() {
-		double d = getTotalCropIncPastureArea() + getLandCoverArea(LandDataType.FOREST) + getLandCoverArea(LandDataType.OTHER_NATURAL);
+		double d = 0;
+		for (LandCoverType l : LandCoverType.values())
+			d += getLandCoverArea(l);
+		
 		return d;
 	}
 	
 	public double getForestToNaturalFraction() {
-		double forest = getLandCoverArea(LandDataType.FOREST);
-		double natural = forest + getLandCoverArea(LandDataType.OTHER_NATURAL);
+		double forest = getLandCoverArea(LandCoverType.FOREST);
+		double natural = forest + getLandCoverArea(LandCoverType.OTHER_NATURAL);
 		double d = forest / natural;
 		return d;
 	}
 
-
 	public CropToDouble getCropChanges(AreasItem prevAreaAggItem) {
 		CropToDouble changes = new CropToDouble();
 		
-		for (Entry<CropType, Double> entry : prevAreaAggItem.cropAreas.entrySet()) {
-			CropType crop = entry.getKey();
-			double change = getCropArea(crop) - entry.getValue();
-			changes.put(crop, change);
+		for (CropType c : CropType.getCropsLessPasture()) {
+			double change = getCropArea(c) - prevAreaAggItem.getCropArea(c);
+			changes.put(c, change);
 		}
 		
 		return changes;
 	}
+
 }

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

@@ -3,21 +3,21 @@ package ac.ed.lurg.landuse;
 import java.util.HashMap;
 import java.util.Map;
 
-import ac.ed.lurg.types.LandDataType;
+import ac.ed.lurg.types.LandCoverType;
 import ac.sac.raster.RasterItem;
 
 /** Used to hold less detailed land-cover information
  *  This is used in the initalisation phase, after that land-use is used */
 public class LandCoverItem implements RasterItem {
 	
-	Map<LandDataType, Double> landcover = new HashMap<LandDataType, Double>();
+	Map<LandCoverType, Double> landcover = new HashMap<LandCoverType, Double>();
 
 	/** Area in Mha */ 
-	public Double getLandCover(LandDataType landType) {
+	public Double getLandCover(LandCoverType landType) {
 		return landcover.get(landType);
 	}
 
-	public void setLandCover(LandDataType landType, double d) {
+	public void setLandCover(LandCoverType landType, double d) {
 		landcover.put(landType, d);
 	}
 		

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

 package ac.ed.lurg.landuse;
 
-import ac.ed.lurg.types.LandDataType;
+import ac.ed.lurg.types.LandCoverType;
 import ac.sac.raster.AbstractRasterReader;
 import ac.sac.raster.RasterSet;
 
 public class LandCoverReader extends AbstractRasterReader<LandCoverItem> {	
 
-	private LandDataType landType;
+	private LandCoverType landType;
 	
-	public LandCoverReader (RasterSet<LandCoverItem> dataset, LandDataType landType) {
+	public LandCoverReader (RasterSet<LandCoverItem> dataset, LandCoverType landType) {
 		super(dataset);
 		this.landType = landType;
 	}

src/ac/ed/lurg/types/CropToDouble.java

@@ -22,13 +22,4 @@ public class CropToDouble extends HashMap<CropType, Double> {
 		
 		return total;
 	}
-	
-	public double getCroplandTotal() {
-		double total = 0;
-		for (CropType c : CropType.getCropsLessPasture())
-			total += get(c);
-		
-		return total;
-
-	}
 }

src/ac/ed/lurg/types/CropType.java

@@ -40,7 +40,7 @@ public enum CropType {
 		Collection<CropType> comms = new HashSet<CropType>();
 		
 		for (CropType c : values())
-			if (c.importedCrop || !c.isMeat)
+			if (c.importedCrop && !c.isMeat)
 				comms.add(c);
 		
 		return comms;

src/ac/ed/lurg/types/LandDataType.java → src/ac/ed/lurg/types/LandCoverType.java

 package ac.ed.lurg.types;
 
-public enum LandDataType {
+public enum LandCoverType {
 	
-//	LAND ("land"),
 	FOREST ("forest"),
 	OTHER_NATURAL ("otherNatural"),
 	CROPLAND ("cropland"),
@@ -10,7 +9,7 @@ public enum LandDataType {
 
 	private String name;
 	
-	LandDataType(String name) {
+	LandCoverType(String name) {
 		this.name = name;
 	}
 	

src/ac/sac/raster/RasterSet.java

@@ -5,6 +5,7 @@ import java.util.Collection;
 import java.util.HashMap;
 import java.util.Map;
 
+import ac.ed.lurg.landuse.AreasItem;
 import ac.ed.lurg.utils.LogWriter;
 
 /* Class which holds raster data.  The generics used to defines the type of raster data held.