Changes to force protected areas to be respected (over time).

src/ac/ed/lurg/ModelConfig.java

@@ -318,5 +318,8 @@ public class ModelConfig {
 	public static final int NUM_YIELD_CLUSTERS = getIntProperty("NUM_YIELD_CLUSTERS", 8000);
 	public static final long RANDOM_SEED = getIntProperty("RANDOM_SEED", 1974329);  // any number will do
 	
-
+	// Protected areas forcing parameters
+	public static final boolean FORCE_PROTECTED_AREAS = IS_CALIBRATION_RUN ? false : getBooleanProperty("FORCE_PROTECTED_AREAS", false);
+	public static final int FORCE_PROTECTED_AREAS_START_YEAR  = getIntProperty("FORCE_PROTECTED_AREAS_START_YEAR", 2020);
+	public static final int FORCE_PROTECTED_AREAS_END_YEAR  = getIntProperty("FORCE_PROTECTED_AREAS_END_YEAR", 2050);
 }

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

@@ -106,7 +106,7 @@ public class GamsLocationOptimiser {
 			String locString = Integer.toString(locationId);
 			LandUseItem landUseItem = entry.getValue();
 			
-			double suitableLand = landUseItem.getSuitableLand();
+			double suitableLand = landUseItem.getSuitableLand(inputData.getTimestep().getYear());
 			totalSuitable += suitableLand;
 			if (DEBUG) LogWriter.println(String.format("  %d   %15s,\t %.3f", locationId, "suitableLand", suitableLand));
 			setGamsParamValueTruncate(landP.addRecord(locString), suitableLand, 3);

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

@@ -76,17 +76,15 @@ public class GamsRasterOptimiser {
 			LogWriter.printlnError("Total Area " + comment + ": " + l.getName() + ": " + total);
 		}
 		
-		double unprotectedNatural=0, suitableArea=0, protectedArea=0, unprotectedArea=0;
+		double unprotectedNatural=0, suitableArea=0, protectedArea=0;
 		for (LandUseItem a : areaRaster.values()) {
 			if (a != null) {
 				protectedArea += a.getProtectedArea();
-				unprotectedArea += a.getUnprotectedArea();
 				unprotectedNatural += a.getUnprotectedNatural();
-				suitableArea += a.getSuitableLand();
+				suitableArea += a.getSuitableLand(rasterInputData.getTimestep().getYear());
 			}
 		}
 		LogWriter.println("Total protectedArea " + comment + ": " + protectedArea);
-		LogWriter.println("Total unprotectedArea " + comment + ": " + unprotectedArea);
 		LogWriter.println("Total unprotectedNatural " + comment + ": " + unprotectedNatural);
 		LogWriter.println("Total suitableArea " + comment + ": " + suitableArea);
 	}
@@ -120,17 +118,17 @@ public class GamsRasterOptimiser {
 				LogWriter.println("croplandChange" + croplandChange);
 			}
 
-			double prevUnprotectedManagedForest = 0, prevUnprotectedUnmanagedForest = 0, prevUnprotectedNatural = 0;  
+			double prevManagedForest = 0, prevUnmanagedForest = 0, prevNatural = 0;  
 			
 			for (LandUseItem luItem: landUseItemsForLocation.values()) {
-				prevUnprotectedManagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.MANAGED_FOREST);
-				prevUnprotectedUnmanagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST);
-				prevUnprotectedNatural += luItem.getUnprotectedNatural();
+				prevManagedForest += luItem.getLandCoverArea(LandCoverType.MANAGED_FOREST);
+				prevUnmanagedForest += luItem.getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
+				prevNatural += luItem.getTotalNatural();
 			}
 			
-			double prevForestTotal = prevUnprotectedManagedForest + prevUnprotectedUnmanagedForest;
-			double prevForestToNaturalFraction = (prevUnprotectedNatural > 0) ? prevForestTotal / prevUnprotectedNatural : 0;
-			double prevForestManagedFraction = (prevForestTotal > 0) ? prevUnprotectedManagedForest / prevForestTotal : 0;
+			double prevForestTotal = prevManagedForest + prevUnmanagedForest;
+			double prevForestToNaturalFraction = (prevNatural > 0) ? prevForestTotal / prevNatural : 0;
+			double prevForestManagedFraction = (prevForestTotal > 0) ? prevManagedForest / prevForestTotal : 0;
 			
 			double pastureFromCrop = 0;
 			double pastureFromNatural = 0;
@@ -335,8 +333,8 @@ public class GamsRasterOptimiser {
 				IrrigationItem aggIrig = aggregatedIrrigCosts.lazyGet(clusterId);
 
 				// Irrigation cost
-				double suitableAreaThisTime  = landUseItem.getSuitableLand();
-				double suitableAreaSoFar = aggLandUse.getSuitableLand();
+				double suitableAreaThisTime  = landUseItem.getSuitableLand(rasterInputData.getTimestep().getYear());
+				double suitableAreaSoFar = aggLandUse.getSuitableLand(rasterInputData.getTimestep().getYear());
 				
 				if (irrigItem!= null) {
 					aggIrig.setIrrigCost( aggregateMean(aggIrig.getIrrigCost(), suitableAreaSoFar, irrigItem.getIrrigCost(), suitableAreaThisTime));
@@ -379,12 +377,12 @@ public class GamsRasterOptimiser {
 					aggLandUse.setIntensity(crop, intensityThisTime);  // intensity currently is always the same within a location, so can just that any.  If this changes need to average here.
 				}
 
-				// Unavailable Fraction - confusingly, this must be done before protected areas and land covers are set due to getUnprotectedArea() call
-				double unprotectedAreaSoFar = aggLandUse.getUnprotectedArea();     
-				double unprotectedAreaThisTime = landUseItem.getUnprotectedArea();
+				// Unavailable Fraction
+				double totalAreaSoFar = aggLandUse.getTotalLandCoverArea();     
+				double totalAreaThisTime = landUseItem.getTotalLandCoverArea();
 				double unavailFractThisTime = landUseItem.getUnavailableFract();
 				double unavailFractSoFar = aggLandUse.getUnavailableFract();
-				double averagedUnavailFract = aggregateMean(unavailFractSoFar, unprotectedAreaSoFar, unavailFractThisTime, unprotectedAreaThisTime);
+				double averagedUnavailFract = aggregateMean(unavailFractSoFar, totalAreaSoFar, unavailFractThisTime, totalAreaThisTime);
 				aggLandUse.setUnavailableFract(averagedUnavailFract);
 				
 				// Protected areas

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

@@ -50,10 +50,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return unavailableFract;
 	}
 
-	public double getUnprotectedArea() {
-		return getTotalLandCoverArea()-getProtectedArea();	
-	}
-	
 	public Intensity getIntensity(CropType crop) {
 		return intensityMap.get(crop);
 	}
@@ -241,43 +237,60 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return total;
 	}
 	
+	public double getTotalNatural() {
+		double totalNatural = getLandCoverArea(LandCoverType.OTHER_NATURAL) + getLandCoverArea(LandCoverType.MANAGED_FOREST) + getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
+		return totalNatural;
+	}
+	
 	/** averages protected areas across land cover types */
 	public double getUnprotectedLandCoverArea(LandCoverType c) {
 		double lcTypeArea = getLandCoverArea(c);
 		
 		if (c.isProtectable()) {
-			double unprotectedNaturalFract = getUnprotectedNaturalFract();
-			return lcTypeArea * unprotectedNaturalFract;
+			double unprotectedNat = getUnprotectedNatural();
+			double totalNatural = getTotalNatural();
+			double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0;
+			return lcTypeArea * unprotectedFract;
 		}
 		else {
 			return lcTypeArea;
 		}
 	}
 	
-	private double getUnprotectedNaturalFract() {
-		double unprotectedNat = getUnprotectedNatural();
-		double totalNatural = getLandCoverArea(LandCoverType.OTHER_NATURAL) + getLandCoverArea(LandCoverType.MANAGED_FOREST) + getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
-		double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0;
-		return unprotectedFract;
-	}
-	
-	public double getUnprotectedNatural() {
+	private double getDesiredProtected() {
 		double totalLand = getTotalLandCoverArea();
 		if (totalLand <=0)
 			return 0;
-		double totalNatural = getLandCoverArea(LandCoverType.OTHER_NATURAL) + getLandCoverArea(LandCoverType.MANAGED_FOREST) + getLandCoverArea(LandCoverType.UNMANAGED_FOREST);
-
 		double protectedFract = getProtectedArea() / totalLand;
 		double minAndProtectedA = totalLand * (1 - (1-protectedFract) * (1-unavailableFract) * (1-ModelConfig.MIN_NATURAL_RATE));
-
-		double unprotectedNat = Math.max(0, totalNatural - minAndProtectedA);  // if we are already using more than is protected then the unprotectedArea is 0
+		return minAndProtectedA;
+	}
+	
+	public double getUnprotectedNatural() {
+		double desiredProtected = getDesiredProtected();
+		double totalNatural = getTotalNatural();
+		double unprotectedNat = Math.max(0, totalNatural - desiredProtected);  // if we are already using more than is protected then the unprotectedArea is 0
 		return unprotectedNat;
 	}
 	
-	public double getSuitableLand() {
-		double totalUnprotectedNatural  = getUnprotectedNatural();
+	public double getSuitableLand(int year) {
+		double suitable;			
 		double currentAgri = getLandCoverArea(LandCoverType.PASTURE) + getLandCoverArea(LandCoverType.CROPLAND);
-		double suitable = currentAgri + totalUnprotectedNatural;
+		double desiredProtected  = getDesiredProtected();
+		double totalNatural = getTotalNatural();
+		double natAvailForAgri =  totalNatural - desiredProtected;  // could be negative, i.e. excess agricultural area already used
+		
+		if (ModelConfig.FORCE_PROTECTED_AREAS && natAvailForAgri < 0 && year >= ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR) {
+			double proportion = 1;
+			
+			if (year < ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR)
+				proportion = 1.0 / (ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR - year);
+			
+			suitable = Math.max(0, currentAgri + natAvailForAgri * proportion); // netNatAvailForAgri is negative, but suitable area < 0 is not sensible (seems to happen with high barren areas)
+		}
+		else
+			suitable = currentAgri + Math.max(0, natAvailForAgri);  // should be identical to currentAgri + getUnprotectedNatural() in this case
+		
 		return suitable;
 	}
 	
@@ -369,6 +382,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 
 	@Override
 	public String toString() {
-		return "LandUseItem: [landCoverAreas=" + landCoverAreas + ", protectedArea=" + protectedArea + "], suitableLand " +  getSuitableLand();
+		return "LandUseItem: [landCoverAreas=" + landCoverAreas + ", protectedArea=" + protectedArea + ", unavailableFract=" + unavailableFract + "]";
 	}
 }
\ No newline at end of file

src/ac/ed/lurg/output/LandUseOutputer.java

@@ -57,7 +57,7 @@ public class LandUseOutputer extends AbstractLandUseOutputer {
 				StringBuffer sbData = new StringBuffer(String.format("%.2f %.2f", lat, lon));
 
 				sbData.append(String.format(" %.8f", item.getTotalLandCoverArea()));
-				sbData.append(String.format(" %.8f", item.getSuitableLand()));
+				sbData.append(String.format(" %.8f", item.getSuitableLand(year)));
 
 				for (LandCoverType cover : LandCoverType.values()) {
 					sbData.append(String.format(" %.8f", item.getLandCoverArea(cover)));