From 9736907d4f37026fcdba05f804414686f9d9f023 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Wed, 20 Jun 2018 11:18:46 +0100
Subject: [PATCH] Fix to issues with PA not being aggregated correctly

---
 .../country/gams/GamsRasterOptimiser.java     | 46 +++++++++++--------
 src/ac/ed/lurg/landuse/LandUseItem.java       | 43 +++++++++--------
 2 files changed, 51 insertions(+), 38 deletions(-)

diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 0d576a6e..bb18c250 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -65,7 +65,7 @@ public class GamsRasterOptimiser {
 		return theCopy;
 	}
 
-	private void checkedTotalAreas(RasterSet<LandUseItem> areaRaster, String comment) {
+	private void checkedTotalAreas(Map<? extends Object, LandUseItem> areaRaster, String comment) {
 		for (LandCoverType l : LandCoverType.values()) {
 			double total = 0;
 			for (LandUseItem a : areaRaster.values()) {
@@ -75,16 +75,17 @@ public class GamsRasterOptimiser {
 			LogWriter.printlnError("Total Area " + comment + ": " + l.getName() + ": " + total);
 		}
 		
-		if (DEBUG) {
-			double protectedArea = 0;
-			double suitableArea = 0;
-			for (LandUseItem a : areaRaster.values()) {
-				protectedArea += a.getProtectedAreaIncMinNatural();
-				suitableArea += a.getSuitableLand();
-			}
-			LogWriter.println("Total protected " + comment + ": " + protectedArea);
-			LogWriter.println("Total suitableArea " + comment + ": " + suitableArea);
+		double protectedAreaIncMinNatural=0, suitableArea=0, protectedArea=0, unprotectedArea=0;
+		for (LandUseItem a : areaRaster.values()) {
+			protectedArea += a.getProtectedArea();
+			unprotectedArea += a.getUnprotectedArea();
+			protectedAreaIncMinNatural += a.getProtectedAreaIncMinNatural();
+			suitableArea += a.getSuitableLand();
 		}
+		LogWriter.println("Total protectedArea " + comment + ": " + protectedArea);
+		LogWriter.println("Total unprotectedArea " + comment + ": " + unprotectedArea);
+		LogWriter.println("Total protectedAreaIncMinNatural " + comment + ": " + protectedAreaIncMinNatural);
+		LogWriter.println("Total suitableArea " + comment + ": " + suitableArea);
 	}
 
 	private RasterSet<LandUseItem> allocAreas(Map<Integer, ? extends LandUseItem> prevAreasAgg, GamsLocationOutput gamsOutput) {
@@ -243,11 +244,11 @@ 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<LandUseItem> landUseRaster = rasterInputData.getPreviousLandUses();
+		checkedTotalAreas(landUseRaster, "before");
 
 	//	LandUseItem test = landUseRaster.getFromCoordinates(71.25, 35.25);
 
-		/*	checkedTotalAreas(cropAreaRaster, "before");
-		new RasterOutputer<AreasItem>(cropAreaRaster, "BeforeLandCover") {
+			/*	new RasterOutputer<AreasItem>(cropAreaRaster, "BeforeLandCover") {
 			@Override
 			public Double getValue(RasterKey location) {
 				AreasItem area = results.get(location);
@@ -372,11 +373,19 @@ 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.
 				}
 
-				// Protected areas
-				double protectedThisTime = landUseItem.getProtectedAreaIncMinNatural();
-				double protectedSoFar = aggLandUse.getProtectedAreaIncMinNatural();
-				aggLandUse.setProtectedArea(protectedSoFar+ protectedThisTime);
+				// 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();
+				double unavailFractThisTime = landUseItem.getUnavailableFract();
+				double unavailFractSoFar = aggLandUse.getUnavailableFract();
+				double averagedUnavailFract = aggregateMean(unavailFractSoFar, unprotectedAreaSoFar, unavailFractThisTime, unprotectedAreaThisTime);
+				aggLandUse.setUnavailableFract(averagedUnavailFract);
 				
+				// Protected areas
+				double protectedThisTime = landUseItem.getProtectedArea();
+				double protectedSoFar = aggLandUse.getProtectedArea();
+				aggLandUse.setProtectedArea(protectedSoFar + protectedThisTime);
+
 				// Land covers ares
 				for (LandCoverType landType : LandCoverType.values()) {
 					double areaThisTime = landUseItem.getLandCoverArea(landType);
@@ -403,10 +412,7 @@ public class GamsRasterOptimiser {
 			} */
 	//	}
 		
-		double baseCropland = LandUseItem.getTotalLandCover(aggregatedAreas.values(), LandCoverType.CROPLAND);
-		double basePasture = LandUseItem.getTotalLandCover(aggregatedAreas.values(), LandCoverType.PASTURE);
-		LogWriter.println("\nbaseCropland=" + baseCropland + ", basePasture=" + basePasture);
-
+		checkedTotalAreas(aggregatedAreas, "after");
 		return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, rasterInputData.getCountryInput());
 	}
 	
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index b0ed042f..819a4d5d 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -10,7 +10,6 @@ import ac.ed.lurg.types.CropToDouble;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.utils.Interpolator;
-import ac.ed.lurg.utils.LogWriter;
 import ac.sac.raster.InterpolatingRasterItem;
 
 public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serializable {
@@ -36,21 +35,34 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		this.protectedArea = protectedArea;
 	}
 	
-	public double getProtectedAreaIncMinNatural() {
-		
-		double protectedLand = 0;
-		double minOrProtectedA = 0;
+	public double getProtectedArea() {
 		if(ModelConfig.PROTECTED_AREAS_ENABLED)
-			protectedLand = protectedArea;
-		
-		double unprotectedLand = getTotalLandCoverArea()-protectedLand;
-		double unsuitableLand = Math.max(unprotectedLand * ModelConfig.MIN_NATURAL_RATE,unprotectedLand*unavailableFract);
-		
-		minOrProtectedA = protectedLand + unsuitableLand;
-		return minOrProtectedA;
-		
+			return protectedArea;
+		else 
+			return 0.0;	
+	}
+	
+	public void setUnavailableFract(double unavailF) {
+		unavailableFract = unavailF;
+	}
+	
+	public double getUnavailableFract() {
+		return unavailableFract;
+	}
+
+	public double getUnprotectedArea() {
+		return getTotalLandCoverArea()-getProtectedArea();	
 	}
 	
+	public double getProtectedAreaIncMinNatural() {
+		double protectedLand = getProtectedArea();
+		double unprotectedLand = getUnprotectedArea();
+		double unsuitableLand = Math.max(unprotectedLand * ModelConfig.MIN_NATURAL_RATE, unprotectedLand*unavailableFract);
+		
+		double minOrProtectedA = protectedLand + unsuitableLand;
+		return minOrProtectedA;	
+	}
+
 	public Intensity getIntensity(CropType crop) {
 		return intensityMap.get(crop);
 	}
@@ -84,7 +96,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return areaTotal > 0 ? fertTotal / areaTotal : 0;
 	}
 
-	
 	public static double getFertiliserTotal(Collection<? extends LandUseItem> items, CropType... crops) {			
 		double total = 0;
 		for (LandUseItem a : items) {
@@ -365,8 +376,4 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	public String toString() {
 		return "LandUseItem: [landCoverAreas=" + landCoverAreas + ", protectedArea=" + protectedArea + "], suitableLand " +  getSuitableLand();
 	}
-
-	public void setUnavailableFract(double unavailableFract) {
-		this.unavailableFract = unavailableFract;
-	}
 }
\ No newline at end of file
-- 
GitLab