From f726ca41642bf691f50f5f17a53ae4259e7dcf8e Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Fri, 28 Aug 2015 11:12:41 +0100
Subject: [PATCH] Interpolate intensities

---
 src/ac/ed/lurg/ModelMain.java                 | 10 ++--
 .../country/gams/GamsRasterOptimiser.java     |  8 ++--
 src/ac/ed/lurg/landuse/Intensity.java         | 11 +++++
 src/ac/ed/lurg/landuse/LandUseItem.java       | 48 ++++++++++---------
 src/ac/ed/lurg/output/LpjgOutputer.java       | 20 ++++----
 src/ac/ed/lurg/utils/Interpolator.java        | 21 ++++++++
 6 files changed, 79 insertions(+), 39 deletions(-)
 create mode 100644 src/ac/ed/lurg/utils/Interpolator.java

diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index d3b9f236..a6827925 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -213,11 +213,15 @@ public class ModelMain {
 			StringBuffer sbData = new StringBuffer();
 			sbData.append(String.format("%d, %.1f, %.1f, %.1f", 
 					timestep.getYear(), 
-					LandUseItem.getTotalLandCover(landUseRaster, LandCoverType.CROPLAND),
-					LandUseItem.getTotalLandCover(landUseRaster, LandCoverType.PASTURE),
-					LandUseItem.getTotalLandCover(landUseRaster, LandCoverType.OTHER_NATURAL)));
+					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.CROPLAND),
+					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.PASTURE),
+					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.OTHER_NATURAL)));
 
 			sbData.append(String.format(", %.1f", LandUseItem.getFertiliserTotal(landUseRaster.values(), CropType.getCropsLessPasture())/1000));
+			
+//			double f = LandUseItem.getFertiliserTotal(landUseRaster.values(), CropType.MAIZE)/1000;
+//			double c = LandUseItem.getTotalCropArea(landUseRaster.values(), CropType.MAIZE);
+
 			sbData.append(String.format(", %.1f", LandUseItem.getFertiliserTotal(landUseRaster.values(), CropType.PASTURE)/1000));
 			
 			sbData.append(String.format(", %.1f", LandUseItem.getIrrigationTotal(landUseRaster.values(), CropType.getCropsLessPasture())));
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 5494bd8e..6df3f1b8 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -61,8 +61,8 @@ public class GamsRasterOptimiser {
 				cropAdjustments.put(c, 1.0);
 			}
 
-			double baseCropland = LandUseItem.getTotalLandCover(gamsInput.getPreviousLandUse(), LandCoverType.CROPLAND);
-			double basePasture = LandUseItem.getTotalLandCover(gamsInput.getPreviousLandUse(), LandCoverType.PASTURE);
+			double baseCropland = LandUseItem.getTotalLandCover(gamsInput.getPreviousLandUse().values(), LandCoverType.CROPLAND);
+			double basePasture = LandUseItem.getTotalLandCover(gamsInput.getPreviousLandUse().values(), LandCoverType.PASTURE);
 
 			LogWriter.println("baseCropland=" + baseCropland + ", basePasture=" + basePasture);
 
@@ -73,7 +73,7 @@ public class GamsRasterOptimiser {
 				result = opti.run();
 
 				if (baseCropland > 0) {
-					double cropAdjSingle = LandUseItem.getTotalLandCover(result.getLandUses(), LandCoverType.CROPLAND) / baseCropland;
+					double cropAdjSingle = LandUseItem.getTotalLandCover(result.getLandUses().values(), LandCoverType.CROPLAND) / baseCropland;
 					cropAdjSingle = Math.min(Math.max(cropAdjSingle, 0.5), 2);
 
 					for (CropType c : CropType.getCropsLessPasture()) {
@@ -82,7 +82,7 @@ public class GamsRasterOptimiser {
 				}
 
 				if (basePasture > 0) {
-					double pastureAdjSingle = LandUseItem.getTotalLandCover(result.getLandUses(), LandCoverType.PASTURE) / basePasture;
+					double pastureAdjSingle = LandUseItem.getTotalLandCover(result.getLandUses().values(), LandCoverType.PASTURE) / basePasture;
 					pastureAdjSingle = Math.min(Math.max(pastureAdjSingle, 0.5), 2);
 					cropAdjustments.put(CropType.PASTURE, Math.min(Math.max(cropAdjustments.get(CropType.PASTURE) * pastureAdjSingle, 0.2), 5));
 				}
diff --git a/src/ac/ed/lurg/landuse/Intensity.java b/src/ac/ed/lurg/landuse/Intensity.java
index 9022ec21..23f62525 100644
--- a/src/ac/ed/lurg/landuse/Intensity.java
+++ b/src/ac/ed/lurg/landuse/Intensity.java
@@ -1,6 +1,7 @@
 package ac.ed.lurg.landuse;
 
 import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.utils.Interpolator;
 
 public class Intensity {
 	private double fertiliserIntensity;
@@ -47,4 +48,14 @@ public class Intensity {
 	public double getIrrigationRate() {
 		return maxIrrigRate * irrigationIntensity;
 	}
+
+	public Intensity(Intensity from, Intensity to, double factor) {
+		this(
+			Interpolator.interpolate(from.fertiliserIntensity, to.fertiliserIntensity, factor), 
+			Interpolator.interpolate(from.irrigationIntensity, to.irrigationIntensity, factor), 
+			Interpolator.interpolate(from.otherIntensity, to.otherIntensity, factor), 
+			Interpolator.interpolate(from.yield, to.yield, factor), 
+			Interpolator.interpolate(from.unitEnergy, to.unitEnergy, factor), 
+			Interpolator.interpolate(from.maxIrrigRate, to.maxIrrigRate, factor));
+	}
 }
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index c625f59a..7838ff59 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -8,6 +8,7 @@ import java.util.Map.Entry;
 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.sac.raster.InterpolatingRasterItem;
 
 public class LandUseItem implements InterpolatingRasterItem<LandUseItem> {
@@ -177,10 +178,9 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem> {
 	
 	public double getTotalLandCoverArea() {
 		double d = 0;
-		for (LandCoverType l : LandCoverType.values()) {
+		for (LandCoverType l : LandCoverType.values()) 
 			d += getLandCoverArea(l);
-		}
-		
+	
 		return d;
 	}
 		
@@ -220,38 +220,42 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem> {
 		for (CropType crop : CropType.values()) {
 			Double from = fromItem.cropFractions.get(crop);
 			Double to = toItem.cropFractions.get(crop);
-			Double d = interpolate(from, to, factor);
+			Double d = Interpolator.interpolate(from, to, factor);
 			cropFractions.put(crop, d);
+			
+			Intensity fromIntensity = fromItem.intensityMap.get(crop);
+			Intensity toIntensity = toItem.intensityMap.get(crop);
+			Intensity interpolateIntensity = toIntensity;  // might still be null
+			
+			if (fromIntensity != null && toIntensity != null)
+				interpolateIntensity = new Intensity(fromIntensity, toIntensity, factor);  // both non-null really interpolate
+			else if (fromIntensity != null)
+				interpolateIntensity = fromIntensity; // just fromIntensity non-null
+			
+			intensityMap.put(crop, interpolateIntensity);
 		}
 		
 		for (LandCoverType landCover : LandCoverType.values()) {
 			Double from = fromItem.landCoverAreas.get(landCover);
 			Double to = toItem.landCoverAreas.get(landCover);
-			Double d = interpolate(from, to, factor);
+			Double d = Interpolator.interpolate(from, to, factor);
 			landCoverAreas.put(landCover, d);
 		}
 	}
 	
-	private Double interpolate(Double from, Double to, double factor) {
-		if (from == null) {
-			if (to == null) 
-				return 0.0;
-			return to;
-		}
-		else if (to == null)
-			return from;
-			
-		Double res = from + factor * (to - from);
-		return res;
-	}
-	
-	public static double getTotalLandCover(Map<?, ? extends LandUseItem> areaRaster, LandCoverType landCover) {
+	public static double getTotalLandCover(Collection<? extends LandUseItem> landUses, LandCoverType landCover) {
 		double total = 0;
-		for (LandUseItem a : areaRaster.values()) {
+		for (LandUseItem a : landUses)
 			total += a.getLandCoverArea(landCover);
-		}
-			
+		
 		return total;
 	}
 
+	public static double getTotalCropArea(Collection<LandUseItem> landUses, CropType crop) {
+		double total = 0;
+		for (LandUseItem a : landUses)
+			total += a.getCropArea(crop);
+		
+		return total;
+	}
 }
diff --git a/src/ac/ed/lurg/output/LpjgOutputer.java b/src/ac/ed/lurg/output/LpjgOutputer.java
index 219c8901..fd9b2713 100644
--- a/src/ac/ed/lurg/output/LpjgOutputer.java
+++ b/src/ac/ed/lurg/output/LpjgOutputer.java
@@ -19,13 +19,13 @@ import ac.sac.raster.RasterSet;
 
 public class LpjgOutputer {
 
-	private RasterSet<LandUseItem> intensityRaster;
+	private RasterSet<LandUseItem> landUseRaster;
 	private YieldRaster yieldSurfaces;
 	private int year;
 
-	public LpjgOutputer(int year, RasterSet<LandUseItem> intensityRaster, YieldRaster yieldSurfaces) {
+	public LpjgOutputer(int year, RasterSet<LandUseItem> landUseRaster, YieldRaster yieldSurfaces) {
 		this.year = year;
-		this.intensityRaster = intensityRaster;
+		this.landUseRaster = landUseRaster;
 		this.yieldSurfaces = yieldSurfaces;
 	}
 
@@ -59,15 +59,15 @@ public class LpjgOutputer {
 			fertWriter.write("Lon Lat TeWWirr TeSWirr TeCoirr TrRiirr Pasture");
 			irrigWriter.newLine();
 
-			for (Entry<RasterKey, LandUseItem> entry : intensityRaster.entrySet()) {
+			for (Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) {
 				RasterKey key = entry.getKey();
 				LandUseItem item = entry.getValue();
 				
 				if (item == null)
 					continue;
 				
-				double lat = intensityRaster.getXCoordin(key);
-				double lon = intensityRaster.getYCoordin(key);
+				double lat = landUseRaster.getXCoordin(key);
+				double lon = landUseRaster.getYCoordin(key);
 				boolean isSpringWheat = isSpringWheat(key);
 
 				double winterWheatIrrig = item.getIrrigationAverageRate(CropType.OILCROPS, isSpringWheat ? null : CropType.WHEAT);
@@ -150,13 +150,13 @@ public class LpjgOutputer {
 			cropFractWriter.write("Lon Lat TeWWirr TeSWirr TeCoirr TrRiirr");
 			cropFractWriter.newLine();
 
-			for (Entry<RasterKey, LandUseItem> entry : intensityRaster.entrySet()) {
+			for (Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) {
 				RasterKey key = entry.getKey();
 				LandUseItem item = entry.getValue();
-				double lat = intensityRaster.getXCoordin(key);
-				double lon = intensityRaster.getYCoordin(key);
+				double lat = landUseRaster.getXCoordin(key);
+				double lon = landUseRaster.getYCoordin(key);
 
-				double expectedArea = intensityRaster.getAreaMha(key);
+				double expectedArea = landUseRaster.getAreaMha(key);
 				double area = item.getTotalLandCoverArea();
 				
 				if (area > 0 && Math.abs((expectedArea-area)/expectedArea) > 0.01) {  // zero area, due to data not being found, already reported so ignored here
diff --git a/src/ac/ed/lurg/utils/Interpolator.java b/src/ac/ed/lurg/utils/Interpolator.java
new file mode 100644
index 00000000..d49da626
--- /dev/null
+++ b/src/ac/ed/lurg/utils/Interpolator.java
@@ -0,0 +1,21 @@
+package ac.ed.lurg.utils;
+
+public class Interpolator {
+
+	public static Double interpolate(Double from, Double to, double factor) {
+		if (from == null) {
+			if (to == null) 
+				return 0.0;
+			return to;
+		}
+		else if (to == null)
+			return from;
+			
+		Double res = interpolate(from.doubleValue(), to.doubleValue(), factor);
+		return res;
+	}
+	
+	public static double interpolate(double from, double to, double factor) {
+		return from + factor * (to - from);
+	}
+}
-- 
GitLab