Changes to fix reading of LPJG file

src/ac/ed/lurg/ModelMain.java

@@ -67,6 +67,10 @@ public class ModelMain {
 
 	private void doTimestep(int timestep) {
 		YieldRaster yieldSurfaces = getYieldSurfaces(timestep);
+		
+//		YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-50.0, -4.0);
+//		LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.CEREALS)  + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+
 		LogWriter.println("Timestep " + timestep);
 		
 		CropToDoubleMap totalArea = new CropToDoubleMap();

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

@@ -25,7 +25,7 @@ import ac.sac.raster.RasterKey;
 import ac.sac.raster.RasterSet;
 
 public class GamsRasterOptimiser {
-	public static final boolean DEBUG = false;
+	public static final boolean DEBUG = true;
 
 	private GamsRasterInput rasterInputData;
 	private LazyHashMap<Integer, Set<RasterKey>> mapping;
@@ -49,28 +49,28 @@ public class GamsRasterOptimiser {
 	private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) {		
 		RasterSet<AreasItem> newAreaRaster = allocAreas(gamsInput.getPreviousAreas(), gamsOutput);
 		RasterSet<IntensitiesItem> newIntensityRaster = allocIntensities(gamsOutput);
-		
+
 		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData());
 	}
 
 	private RasterSet<AreasItem> allocAreas(Map<Integer, ? extends AreasItem> prevAreasAgg, GamsLocationOutput gamsOutput) {
 		RasterSet<AreasItem> newAreaRaster = new RasterSet<AreasItem>();
 		RasterSet<AreasItem> prevAreaRaster = rasterInputData.getPreviousAreas();
-		
+
 		for (Map.Entry<Integer, AreasItem> entry : gamsOutput.getCropAreas().entrySet()) {
 			Integer locId = entry.getKey();
 			AreasItem newAreaAggItem = entry.getValue();
 			AreasItem prevAreaAggItem = prevAreasAgg.get(locId);
 			Set<RasterKey> keys = mapping.get(locId);
-			
+
 			CropToDouble cropChangeTotals = newAreaAggItem.getCropChanges(prevAreaAggItem);
-			if (DEBUG) LogWriter.println("Processing location id " + locId);
+			// if (DEBUG) LogWriter.println("Processing location id " + locId);
 
 			double shortfall = allocAllAreasForAgg(newAreaRaster, prevAreaRaster, keys, cropChangeTotals);
 			if (shortfall > 0)
 				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);
 		}
-		
+
 		return newAreaRaster;
 	}
 
@@ -79,9 +79,9 @@ public class GamsRasterOptimiser {
 		double totalShortfall = 0;
 		Set<RasterKey> keysWithSpace = new HashSet<RasterKey>();
 		CropToDouble avgCropChange = cropChange.scale(1.0/keys.size());
-		
+
 		for (RasterKey key : keys) {
-			if (DEBUG) LogWriter.println("  Processing raster key " + key);
+			//if (DEBUG) LogWriter.println("  Processing raster key " + key);
 			AreasItem prevAreasRasterItem = prevAreaRaster.get(key);
 
 			AreasItem newAreasRasterItem = new AreasItem();
@@ -97,13 +97,13 @@ public class GamsRasterOptimiser {
 		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));
 		}
 
 		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) {
 
@@ -114,7 +114,7 @@ public class GamsRasterOptimiser {
 		double prevNatural = prevAreasRasterItem.getLandCoverArea(LandDataType.OTHER_NATURAL);
 
 		double ratioOfChanges = 1;
-		
+
 		if (netAddedCrop > prevForest + prevNatural) {
 			ratioOfChanges = (prevForest + prevNatural) / netAddedCrop;
 			forestRemoved = prevForest;
@@ -139,9 +139,9 @@ public class GamsRasterOptimiser {
 		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));
+			//	LogWriter.println(String.format("%s, %.1f", crop, additionalArea));
 
 			newAreasRasterItem.setCropArea(crop, prevCropArea + additionalArea);
 		}
@@ -171,11 +171,32 @@ public class GamsRasterOptimiser {
 		RasterSet<AreasItem> cropAreaRaster = rasterInputData.getPreviousAreas();
 		RasterSet<IrrigationCostItem> irrigCostRaster = rasterInputData.getIrrigationCost();
 
-		List<Double> cerealDivisions = getDivisions(yieldRaster, CropType.CEREALS, 2);
-		List<Double> pastureDivisions = getDivisions(yieldRaster, CropType.MEAT_OR_PASTURE, 2);
+		{
+			//		YieldResponsesItem yresp = yieldRaster.getFromCoordinates(-118.0, -33.0);
+			//		LogWriter.printlnError("Test key2: " + yresp.getYieldMax(CropType.CEREALS)  + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+		}
+
+		// look for inconsistencies
+		for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
+
+			YieldResponsesItem yresp = entry.getValue();
+			RasterKey key = entry.getKey();
+
+			if (yresp != null && (yresp.getYieldMax(CropType.CEREALS) < yresp.getYieldFertOnly(CropType.CEREALS) || yresp.getYieldMax(CropType.CEREALS) < yresp.getYieldIrrigOnly(CropType.CEREALS))) {
+				LogWriter.printlnError(key.getCol() + ", "  + key.getRow() + ": " + yieldRaster.getXCoordin(key) + ", " + yieldRaster.getYCoordin(key) + ": " + yresp.getYieldMax(CropType.CEREALS)  + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+			}
+			//	else if (yresp != null)
+			//		LogWriter.println(key.getCol() + ", "  + key.getRow() + ": " + yieldRaster.getXCoordin(key) + ", " + yieldRaster.getYCoordin(key) + ": " + yresp.getYieldMax(CropType.CEREALS)  + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+
+		}
 
-		int numCerealCats = (1+cerealDivisions.size());
-		if (DEBUG) LogWriter.println("Making " + numCerealCats * (1+pastureDivisions.size()) + " categories");
+		int numCerealCats = 4;
+		int numPastureCats = 4;
+
+		List<Double> cerealDivisions = getDivisions(yieldRaster, CropType.CEREALS, numCerealCats);
+		List<Double> pastureDivisions = getDivisions(yieldRaster, CropType.MEAT_OR_PASTURE, numPastureCats);
+
+		if (DEBUG) LogWriter.println("Making " + numCerealCats * numPastureCats + " categories");
 
 		LazyHashMap<Integer, AveragingYieldResponsesItem> aggregatedYields = new LazyHashMap<Integer, AveragingYieldResponsesItem>() { 
 			protected AveragingYieldResponsesItem createValue() { return new AveragingYieldResponsesItem(); }
@@ -190,14 +211,18 @@ public class GamsRasterOptimiser {
 			protected Set<RasterKey> createValue() { return new HashSet<RasterKey>(); }
 		};
 
+		int countFound = 0, countMissing = 0;
+
 		for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
 			YieldResponsesItem yresp = entry.getValue();
 			RasterKey key = entry.getKey();
-			
+
 			if (yresp == null) {
+				countMissing++;
 				LogWriter.printlnError("No YieldResponsesItem for " + key + ", x=" + yieldRaster.getXCoordin(key) + ", y=" + yieldRaster.getYCoordin(key));
 			}
 			else {
+				countFound++;
 				//LogWriter.println("YieldResponsesItem found for " + key + ", x=" + yieldRaster.getXCoordin(key) + ", y=" + yieldRaster.getYCoordin(key));
 
 				AreasItem cropAreas  = cropAreaRaster.get(key);
@@ -205,18 +230,18 @@ public class GamsRasterOptimiser {
 					LogWriter.printlnError("GamsRasterOptimiser: Can't find cropAreas for " + key);
 					continue;
 				}
-				
+
 				IrrigationCostItem irrigCost = irrigCostRaster.get(key);
-						
+
 				int cerealCat = findCategory(cerealDivisions, yresp.getYieldMax(CropType.CEREALS));
 				int pastureCat = findCategory(pastureDivisions, yresp.getYieldMax(CropType.MEAT_OR_PASTURE));
 				Integer id = cerealCat + pastureCat * numCerealCats;
-	
+
 				AveragingYieldResponsesItem avgYResp = aggregatedYields.lazyGet(id);
 				AreasItem aggAreas = aggregatedAreas.lazyGet(id);
 				IrrigationCostItem avgIrigCost = aggregatedIrrigCosts.lazyGet(id);
 				mapping.lazyGet(id).add(entry.getKey()); 
-	
+
 				// Irrigation cost
 				double suitableAreaThisTime  = cropAreas.getSuitableLand();
 				if (irrigCost!= null && suitableAreaThisTime > 0) {
@@ -230,12 +255,12 @@ public class GamsRasterOptimiser {
 					for (YieldType yieldType : YieldType.values()) {
 						avgYResp.setYield(yieldType, crop, yresp.getYield(yieldType, crop));
 					}
-					
+
 					double areaSoFar = aggAreas.getCropArea(crop);
 					double areaThisTime = cropAreas.getCropArea(crop);
 					aggAreas.setCropArea(crop, areaSoFar + areaThisTime);
 				}
-				
+
 				// Land cover, we only need forest and other natural
 				for (LandDataType landType : Arrays.asList(LandDataType.FOREST, LandDataType.OTHER_NATURAL)) {
 					double areaThisTime = cropAreas.getLandCoverArea(landType);
@@ -244,13 +269,14 @@ public class GamsRasterOptimiser {
 				}
 			}
 		}
-		
-		if (DEBUG) {
-			for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) {
-				LogWriter.println(e.getKey() + " category has " + e.getValue().size() + " raster areas");
-				for (RasterKey key : e.getValue()) {
-					LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.MEAT_OR_PASTURE), yieldRaster.get(key).getYieldMax(CropType.CEREALS)));
-				}
+
+		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()) {
+				LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.MEAT_OR_PASTURE), yieldRaster.get(key).getYieldMax(CropType.CEREALS)));
 			}
 		}
 
@@ -262,7 +288,7 @@ public class GamsRasterOptimiser {
 		int numDivisions = divisions.size();
 
 		for (category = 0; category<numDivisions; category++) {
-			if (yield > divisions.get(category)) {
+			if (yield < divisions.get(category)) {
 				break;
 			}
 		}
@@ -274,28 +300,31 @@ public class GamsRasterOptimiser {
 
 		for (YieldResponsesItem yresp : yieldRaster.values()) {
 			if (yresp == null) {
-				//LogWriter.println("GamsRasterOptimiser: Can't get value for crop " + crop);
+				if (DEBUG) LogWriter.println("GamsRasterOptimiser: Can't get value for crop " + crop);
 			}
 			else {
 				double d = yresp.getYieldMax(crop);
-				
-				if (Double.isNaN(d))
-					LogWriter.println("GamsRasterOptimiser: Got NaN value for crop " + crop);
-				else
+				//LogWriter.println("GamsRasterOptimiser: Got value for crop " + crop + " of " + d);
+
+				if (Double.isNaN(d) || d == 0.0) {
+					if (DEBUG) LogWriter.println("GamsRasterOptimiser: Got NaN or zero value for crop " + crop);
+				}
+				else {
 					yieldValues.add(d);
+				}
 			}
 		}
 
 		if (yieldValues.size() == 0) {
 			throw new RuntimeException("No yield values for country, crop = " + crop);
 		}
-		
+
 		Collections.sort(yieldValues);
 
 		List<Double> divisions = new ArrayList<Double>();
 
 		for (int i=1; i < numCategories; i++) {
-			double d = yieldValues.get(yieldValues.size()/numCategories);
+			double d = yieldValues.get(yieldValues.size()/numCategories*i);
 			divisions.add(d);
 		}	
 

src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java

@@ -55,21 +55,29 @@ public class LPJYieldResponseMapReader {
 						record.y = d;
 						break;
 					case 2:
-						record.teWW = d;
+						record.teSW = d;
 						break;
 					case 3:
-						record.trRi = d;
+						record.teSWirr = d;
 						break;
 					case 4:
-						record.teCo = d;
+						record.teWW = d;
 						break;
 					case 5:
-						record.teSo = d;
+						record.teWWirr = d;
 						break;
 					case 6:
-						record.teWWirr = d;
+						record.teCo = d;
+						break;
+					case 7:
+						record.teCoirr = d;
+						break;
+					case 8:
+						record.trRi = d;
+						break;
+					case 9:
+						record.trRiirr = d;
 						break;
-
 					}
 
 					col++;
@@ -99,18 +107,21 @@ public class LPJYieldResponseMapReader {
 		
 		//YieldRecord values in kg DM / m2, we want t/ha, so 10 times larger
 		for (CropType crop : CropType.getAllItems()) {
-			item.setYield(noIrrigYieldType, crop, record.teWW * 10);
-			item.setYield(maxIrrigYieldType, crop, record.teWWirr * 10);
+			item.setYield(noIrrigYieldType, crop, record.teSW * 10);
+			item.setYield(maxIrrigYieldType, crop, record.teSWirr * 10);
 		}
 	}
 	
 	class YieldRecord {
 		double x;
 		double y;
+		double teSW;
+		double teSWirr;
 		double teWW;
 		double teWWirr;
-		double trRi;
 		double teCo;
-		double teSo;
+		double teCoirr;
+		double trRi;
+		double trRiirr;
 	}
 }
\ No newline at end of file

src/ac/ed/lurg/yield/YieldResponse.java

@@ -7,7 +7,7 @@ import ac.ed.lurg.types.YieldType;
 import ac.ed.lurg.utils.LogWriter;
 
 public class YieldResponse {
-	private Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
+	Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
 	private double fertParm;
 	private double irrigParm;
 
@@ -37,7 +37,7 @@ public class YieldResponse {
 	
 	public double getIrrigParam() {
 		if (irrigParm == 0)
-			irrigParm = calcParam (0, 0.6, 1);  // we don't have a mide irrigation figure, so lets assume 60% at mid point
+			irrigParm = calcParam (0, 0.6, 1);  // we don't have a mid irrigation figure, so lets assume 60% at mid point
 		
 		return irrigParm;
 	}

src/ac/sac/raster/RasterKey.java

@@ -77,6 +77,6 @@ public class RasterKey implements Serializable {
 	}	
 	
 	public String toString() {
-		return "row:" + row + ", col:" + col;
+		return "col:" + col + ", row:" + row;
 	}
 }

src/ac/sac/raster/RasterSet.java

@@ -75,8 +75,8 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 	}
 	
 	public double getYCoordin(RasterKey key) {
-		double x = header.getCellSize() * (key.getRow() + 1) + header.getYllCorner();
-		return x;
+		double y = header.getCellSize() * (header.getNrows() - 1 - key.getRow()) + header.getYllCorner();
+		return y;
 	}
 
 	/** Method to really get the data, or create it. Assumes the col and row are in the internal header format */
@@ -123,6 +123,7 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 	
 	public RasterSet<D> popSubsetForKeys(RasterSet<D> subset, Collection<RasterKey> keys) {
 		for (RasterKey key : keys) {
+			//LogWriter.println("popSubsetForKeys: " + key.getCol() + ", "  + key.getRow() + ": " + getXCoordin(key) + ", " + getYCoordin(key));
 			subset.put(key, get(key));
 		}
 		return subset;

src/ac/sac/raster/RasterSetTest.java

+package ac.sac.raster;
+
+import org.junit.Test;
+
+import ac.ed.lurg.yield.YieldRaster;
+import ac.ed.lurg.yield.YieldResponsesItem;
+
+public class RasterSetTest {
+
+	@Test
+	public void test() {
+		YieldRaster dataset = new YieldRaster(new RasterHeaderDetails(720,360,-180,-90,0.5,""));
+		double y = dataset.getYCoordin(new RasterKey(0,179));
+		System.out.println("y=" +y);
+		
+		
+	}
+
+}