Fix for land allocation error in allocAreas method. Unprotected natural area...

Fix for land allocation error in allocAreas method. Unprotected natural area was being calculated on more that one occasion within a timestep resulting in changing proportions of 'other natural' and 'forest' that should be protected. Area(s) of natural land that should be protected now updated at the start of each timestep and remain constant throughout the timestep.

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

@@ -110,11 +110,11 @@ public class GamsRasterOptimiser {
 			double gamsPastureChange = newLandUseAggItem.getLandCoverArea(LandCoverType.PASTURE) - prevLandUseAggItem.getLandCoverArea(LandCoverType.PASTURE);
 			double gamsCroplandChange = newLandUseAggItem.getLandCoverArea(LandCoverType.CROPLAND) - prevLandUseAggItem.getLandCoverArea(LandCoverType.CROPLAND);
 			RasterSet<LandUseItem> landUseItemsForLocation = newLandUseRaster.createSubsetForKeys(keys);
-
+			
 			if (DEBUG) {
 				checkedTotalAreas(landUseItemsForLocation, locId + " before");
-				LogWriter.println("pastureChange" + gamsPastureChange);
-				LogWriter.println("croplandChange" + gamsCroplandChange);
+				LogWriter.println("pastureChange " + gamsPastureChange);
+				LogWriter.println("croplandChange " + gamsCroplandChange);
 			}
 
 			double prevManagedForest = 0, prevUnmanagedForest = 0, prevNatural = 0, protectedAreasPastureChange = 0, protectedAreasCroplandChange = 0;  
@@ -122,13 +122,13 @@ public class GamsRasterOptimiser {
 			for (LandUseItem luItem: landUseItemsForLocation.values()) {
 				prevManagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.MANAGED_FOREST);
 				prevUnmanagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST);
-				prevNatural += luItem.getUnprotectedNatural();
+				prevNatural += luItem.getTotalUnprotectedNatural();
 			}
 
 			double prevForestTotal = prevManagedForest + prevUnmanagedForest;
 			double prevForestToNaturalFraction = (prevNatural > 0) ? prevForestTotal / prevNatural : 0;
 			double prevForestManagedFraction = (prevForestTotal > 0) ? prevManagedForest / prevForestTotal : 0;
-
+			
 			//this tries to honour spatial distribution of protected areas, force removal of pasture and cropland from protected areas 
 			if (ModelConfig.FORCE_PROTECTED_AREAS && year >= ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR) {
 				double protectedLandShortfall = 0;
@@ -206,9 +206,9 @@ public class GamsRasterOptimiser {
 			}
 
 			if (DEBUG) {
-				LogWriter.println("pastureFromCrop" + pastureFromCrop);
-				LogWriter.println("pastureFromNatural" + pastureFromNatural);
-				LogWriter.println("cropFromNatural" + cropFromNatural);
+				LogWriter.println("pastureFromCrop " + pastureFromCrop);
+				LogWriter.println("pastureFromNatural " + pastureFromNatural);
+				LogWriter.println("cropFromNatural " + cropFromNatural);
 			}
 
 			allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.PASTURE, LandCoverType.MANAGED_FOREST, pastureFromNatural * prevForestToNaturalFraction * prevForestManagedFraction, locId);
@@ -232,7 +232,7 @@ public class GamsRasterOptimiser {
 				}
 			}
 		}
-
+		
 		return newLandUseRaster;
 	}
 
@@ -262,7 +262,7 @@ public class GamsRasterOptimiser {
 		double totalShortfall = 0;
 		Set<RasterKey> keysWithSpace = new HashSet<RasterKey>();
 		double totalUnprotectedLC = 0;
-
+		
 		for (RasterKey key : keys) {
 			LandUseItem newLandUseItem = newLandUseRaster.get(key);
 			totalUnprotectedLC += newLandUseItem.getUnprotectedLandCoverArea(fromType);

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

@@ -18,6 +18,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	private Map<CropType, Intensity> intensityMap = new HashMap<CropType, Intensity>();
 	private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
 	private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>();
+	private Map<LandCoverType, Double> unprotectedAreas = new HashMap<LandCoverType, Double>();
 	private double protectedArea; //protected area in Mha
 	private double unavailableArea; //area unavailable due to altitude etc 
 	private double suitableArea;
@@ -42,6 +43,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		intensityMap.putAll(luItemToCopy.intensityMap);
 		cropFractions.putAll(luItemToCopy.cropFractions);
 		landCoverAreas.putAll(luItemToCopy.landCoverAreas);
+		unprotectedAreas.putAll(luItemToCopy.unprotectedAreas);
 		protectedArea = (luItemToCopy.protectedArea);
 		suitableArea = (luItemToCopy.suitableArea);
 	}
@@ -193,7 +195,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 
         double prevTo = getLandCoverArea(toType);
         double prevFrom = getLandCoverArea(fromType);
-        double availLC = fromType.isProtectable() ? Math.min(prevFrom, getUnprotectedNatural()) : prevFrom;  // if protected can not go past protected area
+        double availLC = fromType.isProtectable() ? Math.min(prevFrom, unprotectedAreas.get(fromType)) : prevFrom;  // if protected can not go past protected area
         
         double actuallyChanged = Math.min(availLC, changeReq);
 
@@ -237,7 +239,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	}
 	
 	public void setLandCoverArea(LandCoverType c, double d) {
-		if (Double.isNaN(d) || Double.isInfinite(d))
+		if (Double.isNaN(d) || Double.isInfinite(d) || d < 0)
 			throw new RuntimeException("AreasItem for " + c + " is " + d);
 		
 		landCoverAreas.put(c, d);
@@ -280,19 +282,18 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return totalNatural;
 	}
 	
-	/** averages protected areas across land cover types */
-	public double getUnprotectedLandCoverArea(LandCoverType c) {
-		double lcTypeArea = getLandCoverArea(c);
-		
-		if (c.isProtectable()) {
-			double unprotectedNat = getUnprotectedNatural();
-			double totalNatural = getTotalNatural();
-			double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0;
-			return lcTypeArea * unprotectedFract;
-		}
-		else {
-			return lcTypeArea;
+	public double getTotalUnprotectedNatural() {	
+		double unprotectedNatural = 0;
+		for (LandCoverType landType : LandCoverType.values()) {
+			if (landType.isProtectable()) {
+				unprotectedNatural += unprotectedAreas.get(landType);
+				}
 		}
+		return unprotectedNatural;
+	}
+	
+	public double getUnprotectedLandCoverArea(LandCoverType landType) {
+		return unprotectedAreas.get(landType);
 	}
 
 	private double getProtectedandUnavailable() {
@@ -303,16 +304,23 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return Math.min(totalSuitableLandCover, minAndProtectedA);
 	}
 	
-	
-	public double getUnprotectedNatural() {
+	/** averages protected areas across land cover types */
+	private void updateUnprotectedAreas() {
 		double desiredProtected = getProtectedandUnavailable();
 		double totalNatural = getTotalNatural();
 		double unprotectedNat = Math.max(0, totalNatural - desiredProtected);  // if we are already using more than is protected then the unprotectedArea is 0
-		//LogWriter.println(" totalNatural - desiredProtected " + (totalNatural - desiredProtected) + "unprotectedNat " + unprotectedNat);
-		return unprotectedNat;
+		double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0;
+		
+		for (LandCoverType landType : LandCoverType.values()) {
+			if (landType.isProtectable()) {
+		      unprotectedAreas.put(landType, getLandCoverArea(landType)*unprotectedFract);
+			}
+			else unprotectedAreas.put(landType, getLandCoverArea(landType));
+			}
 	}
 	
 	public void updateSuitableArea(int year) {			
+		updateUnprotectedAreas(); //update unprotected areas 
 		double currentAgri = getLandCoverArea(LandCoverType.PASTURE) + getLandCoverArea(LandCoverType.CROPLAND);
 		double totalNatural = getTotalNatural();
 		double natAvailForAgri =  totalNatural - getProtectedandUnavailable();  // could be negative, i.e. excess agricultural area already used
@@ -327,7 +335,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 			suitableArea = Math.max(0, currentAgri + natAvailForAgri * proportion); // netNatAvailForAgri is negative, but suitable area < 0 is not sensible (seems to happen with high barren areas)	
 		}
 		else //not honouring protected areas if agriculture > natural that should be protected
-			suitableArea = currentAgri + getUnprotectedNatural();  
+			suitableArea = currentAgri + getTotalUnprotectedNatural();  
 	}
 	
 	public double getForestManagedFraction() {