From 923ec8bdab8d7b9a2ea340a886bab9e86b0788d4 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Mon, 24 Aug 2015 15:23:53 +0100
Subject: [PATCH] Output fertiliser application rates

---
 src/ac/ed/lurg/ModelConfig.java               |  5 ++
 src/ac/ed/lurg/ModelMain.java                 | 28 +++----
 .../country/gams/GamsLocationOptimiser.java   | 16 ++--
 src/ac/ed/lurg/landuse/IntensitiesItem.java   | 50 ++++++++++-
 src/ac/ed/lurg/landuse/Intensity.java         |  6 ++
 src/ac/ed/lurg/output/LpjgOutputer.java       | 84 +++++++++++++++++--
 src/ac/ed/lurg/yield/YieldResponse.java       |  5 +-
 7 files changed, 164 insertions(+), 30 deletions(-)

diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 3283c8da..344e9e74 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -123,6 +123,11 @@ public class ModelConfig {
 	public static final int TIMESTEP_SIZE = getIntProperty("TIMESTEP_SIZE", 5);
 	public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2010);
 
+	// Fertiliser application rates in kg/ha
+	public static final double MIN_FERT_AMOUNT = 5;
+	public static final double MID_FERT_AMOUNT = 50;
+	public static final double MAX_FERT_AMOUNT = 200;
+
 	// Other model parameters
 	public static final boolean KEEP_DEMAND_FIXED = getBooleanProperty("KEEP_DEMAND_FIXED", false);
 	public static final String SSP_SCENARIO = getProperty("SSP_SCENARIO", "SSP1_v9_130325");
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 25764f8f..6af6aa81 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -206,28 +206,26 @@ public class ModelMain {
 		}
 	}
 
-	private double getTotalArea(LandCoverType lcType, RasterSet<AreasItem> globalCropAreaRaster) {
-		double totalArea = 0;
-		for (AreasItem item : globalCropAreaRaster.values()) {
-			totalArea += item.getLandCoverArea(lcType);
-		}
-		return totalArea;
-	}
-
-	private void writeMarketFile(Timestep timestep, RasterSet<AreasItem> cropAreaRaster) {
+	private void writeMarketFile(Timestep timestep, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster) {
 		try {
-			StringBuffer sbHeadings = new StringBuffer("Year, Cropland (Mha), Pasture (Mha), Natural (Mha)");
+			StringBuffer sbHeadings = new StringBuffer("Year, Cropland (Mha), Pasture (Mha), Natural (Mha), Fert crop (Gt), Fert pasture (Gt)");
 			for (CropType crop : CropType.getImportedTypes())
 				sbHeadings.append(",Px_" + crop.getGamsName());
 
 			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.TOTAL_LAND_COVER_FILE, sbHeadings.toString());
 
 			StringBuffer sbData = new StringBuffer();
-			sbData.append(String.format("%d,%.1f,%.1f,%.1f", 
-					timestep.getYear(), getTotalArea(LandCoverType.CROPLAND, cropAreaRaster), getTotalArea(LandCoverType.PASTURE, cropAreaRaster), getTotalArea(LandCoverType.OTHER_NATURAL, cropAreaRaster)));
-
+			sbData.append(String.format("%d, %.1f, %.1f, %.1f", 
+					timestep.getYear(), 
+					AreasItem.getTotalLandCover(cropAreaRaster, LandCoverType.CROPLAND),
+					AreasItem.getTotalLandCover(cropAreaRaster, LandCoverType.PASTURE),
+					AreasItem.getTotalLandCover(cropAreaRaster, LandCoverType.OTHER_NATURAL)));
+
+			sbData.append(String.format(", %.1f", IntensitiesItem.getFertiliserTotal(intensityRaster.values(), CropType.getCropsLessPasture())/1000));
+			sbData.append(String.format(", %.1f", IntensitiesItem.getFertiliserTotal(intensityRaster.values(), CropType.PASTURE)/1000));
+			
 			for (CropType crop : CropType.getImportedTypes() )
-				sbData.append(String.format(",%.3f", prevWorldPrices.get(crop)));
+				sbData.append(String.format(", %.3f", prevWorldPrices.get(crop)));
 
 			outputFile.write(sbData.toString());
 			outputFile.newLine();
@@ -241,7 +239,7 @@ public class ModelMain {
 	private void outputTimestepResults(Timestep timestep, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, 
 			RasterSet<IntegerRasterItem> locationIdRaster, YieldRaster yieldSurfaces) {
 
-		writeMarketFile(timestep, cropAreaRaster);
+		writeMarketFile(timestep, intensityRaster, cropAreaRaster);
 
 		if (ModelConfig.OUTPUT_FOR_LPJG) {
 			for (int outputYear : timestep.getYearsFromLast()) {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index fa044bb5..9e83aea1 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -219,13 +219,13 @@ public class GamsLocationOptimiser {
 		double totalPastureArea = 0;
 		double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, feedAmount, netImport, yield, unitEnergy, cropAdj, prod, prodCost;
 
-		LazyHashMap<Integer, IntensitiesItem> intensities = new LazyHashMap<Integer, IntensitiesItem>() {
-			protected IntensitiesItem createValue() { return new IntensitiesItem(); }
-		};
-		LazyHashMap<Integer, AreasItem> cropAreas = new LazyHashMap<Integer, AreasItem>() { 
+		final LazyHashMap<Integer, AreasItem> cropAreas = new LazyHashMap<Integer, AreasItem>() { 
 			protected AreasItem createValue() { return new AreasItem(); }
 		};
 		
+		HashMap<Integer, IntensitiesItem> intensities = new HashMap<Integer, IntensitiesItem>();
+		
+		
 		Map<CropType, CropUsageData> cropUsageData = new HashMap<CropType, CropUsageData>();
 		Map<CropType, Double> cropAdjs = new HashMap<CropType, Double>();
 
@@ -256,7 +256,13 @@ public class GamsLocationOptimiser {
 
 			if (area > 0) { 
 				if (DEBUG) LogWriter.println(String.format("\t location %s, %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f", locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity));
-				IntensitiesItem intensityItem = intensities.lazyGet(locId);
+				
+				IntensitiesItem intensityItem = intensities.get(locId);
+				if (intensityItem == null) {
+					intensityItem = new IntensitiesItem(cropAreas.lazyGet(locId));
+					intensities.put(locId, intensityItem);
+				}
+				
 				intensityItem.setIntensity(cropType, new Intensity(fertIntensity, irrigIntensity, otherIntensity, yield, unitEnergy));
 			}
 
diff --git a/src/ac/ed/lurg/landuse/IntensitiesItem.java b/src/ac/ed/lurg/landuse/IntensitiesItem.java
index 90c21543..2ae1b174 100644
--- a/src/ac/ed/lurg/landuse/IntensitiesItem.java
+++ b/src/ac/ed/lurg/landuse/IntensitiesItem.java
@@ -1,5 +1,6 @@
 package ac.ed.lurg.landuse;
 
+import java.util.Collection;
 import java.util.HashMap;
 import java.util.Map;
 
@@ -8,7 +9,12 @@ import ac.sac.raster.RasterItem;
 
 public class IntensitiesItem implements RasterItem {
 
-	Map<CropType, Intensity> intensityMap = new HashMap<CropType, Intensity>();
+	private Map<CropType, Intensity> intensityMap = new HashMap<CropType, Intensity>();
+	private AreasItem areaItem;
+	
+	public IntensitiesItem(AreasItem areaItem) {
+		this.areaItem = areaItem;
+	}
 	
 	public Intensity getIntensity(CropType crop) {
 		return intensityMap.get(crop);
@@ -17,4 +23,46 @@ public class IntensitiesItem implements RasterItem {
 	public void setIntensity(CropType crop, Intensity intensityData) {
 		intensityMap.put(crop, intensityData);
 	}
+	
+	
+	public double getFertiliserRate(CropType crop) {
+		Intensity i = getIntensity(crop);
+		return (i == null) ? 0 : i.getFertiliserAmount();
+	}
+		
+	public double getFertiliserAmount(CropType c) {
+		double rate = getFertiliserRate(c);
+		double area = areaItem.getCropArea(c);
+		return rate * area;
+	}
+	
+	public double getFertiliserAverageRate(CropType... crops) {
+		double fertTotal = 0;
+		double areaTotal = 0;
+		
+		for (CropType c : crops) {
+			fertTotal += getFertiliserAmount(c);
+			areaTotal = areaItem.getCropArea(c);
+		}
+		
+		return areaTotal > 0 ? fertTotal / areaTotal : 0;
+	}
+
+	
+	public static double getFertiliserTotal(Collection<? extends IntensitiesItem> items, CropType... crops) {			
+		double total = 0;
+		for (IntensitiesItem a : items) {
+			if (a == null)
+				continue;
+			
+			for (CropType c : crops)
+				total += a.getFertiliserAmount(c);
+		}
+			
+		return total;
+	}
+
+	public static double getFertiliserTotal(Collection<? extends IntensitiesItem> items, Collection<CropType> crops) {
+		return getFertiliserTotal(items, crops.toArray(new CropType[crops.size()]));
+	}
 }
diff --git a/src/ac/ed/lurg/landuse/Intensity.java b/src/ac/ed/lurg/landuse/Intensity.java
index 6844daaf..71a9a22a 100644
--- a/src/ac/ed/lurg/landuse/Intensity.java
+++ b/src/ac/ed/lurg/landuse/Intensity.java
@@ -1,5 +1,7 @@
 package ac.ed.lurg.landuse;
 
+import ac.ed.lurg.ModelConfig;
+
 public class Intensity {
 	double fertiliserIntensity;
 	double irrigationIntensity;
@@ -19,6 +21,10 @@ public class Intensity {
 	public double getFertiliserIntensity() {
 		return fertiliserIntensity;
 	}
+	
+	public double getFertiliserAmount() {
+		return ModelConfig.MIN_FERT_AMOUNT + fertiliserIntensity * (ModelConfig.MAX_FERT_AMOUNT - ModelConfig.MIN_FERT_AMOUNT);
+	}
 
 	public double getIrrigationIntensity() {
 		return irrigationIntensity;
diff --git a/src/ac/ed/lurg/output/LpjgOutputer.java b/src/ac/ed/lurg/output/LpjgOutputer.java
index ff67cce6..fbb15be8 100644
--- a/src/ac/ed/lurg/output/LpjgOutputer.java
+++ b/src/ac/ed/lurg/output/LpjgOutputer.java
@@ -48,8 +48,81 @@ public class LpjgOutputer {
 	}
 
 	private void writeIntensity(File outputDir) {
-		// TODO Auto-generated method stub
-		
+		BufferedWriter fertWriter = null;
+		BufferedWriter irrigWriter = null;
+
+		try {
+			String landCoverFileName = outputDir.getPath() + File.separator + "Fert.txt";
+			fertWriter = new BufferedWriter(new FileWriter(landCoverFileName, false));
+			fertWriter.write("Lon Lat TeWWirr TeSWirr TeCoirr TrRiirr Pasture");
+			fertWriter.newLine();
+
+			String cropFractionFileName = outputDir.getPath() + File.separator + "Irrig.txt";
+			irrigWriter = new BufferedWriter(new FileWriter(cropFractionFileName, false));
+			fertWriter.write("Lon Lat TeWWirr TeSWirr TeCoirr TrRiirr Pasture");
+			irrigWriter.newLine();
+
+			for (Entry<RasterKey, IntensitiesItem> entry : intensityRaster.entrySet()) {
+				RasterKey key = entry.getKey();
+				IntensitiesItem item = entry.getValue();
+				
+				if (item == null)
+					continue;
+				
+				double lat = cropAreaRaster.getXCoordin(key);
+				double lon = cropAreaRaster.getYCoordin(key);
+								
+	/*			double crop = item.getLandCoverArea(LandCoverType.CROPLAND);
+				double pasture = item.getLandCoverArea(LandCoverType.PASTURE);
+				double forest = item.getLandCoverArea(LandCoverType.FOREST);
+				double otherNatural = item.getLandCoverArea(LandCoverType.OTHER_NATURAL);
+				double barren = item.getLandCoverArea(LandCoverType.BARREN);
+				irrigWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, year, crop/area, pasture/area, (forest+otherNatural)/area, barren/area));
+				irrigWriter.newLine();
+
+	*/
+				boolean isSpringWheat = isSpringWheat(key);
+
+				double winterWheatFert = item.getFertiliserAverageRate(CropType.OILCROPS, isSpringWheat ? null : CropType.WHEAT);
+				double springWheatFert = item.getFertiliserAverageRate(CropType.PULSES, CropType.STARCHY_ROOTS, isSpringWheat ? CropType.WHEAT : null);
+				double cornFert = item.getFertiliserAverageRate(CropType.MAIZE);
+				double riceFert = item.getFertiliserAverageRate(CropType.RICE);
+				double pastureFert = item.getFertiliserAverageRate(CropType.PASTURE);
+				fertWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, year, winterWheatFert, springWheatFert, cornFert, riceFert, pastureFert));
+				fertWriter.newLine();
+			}
+		}
+		catch (IOException e) {
+			e.printStackTrace();
+		}
+		finally {
+			if (fertWriter != null) {
+				try {
+					fertWriter.close();
+				} 
+				catch (IOException e) {
+					e.printStackTrace();
+				}
+			}
+			if (irrigWriter != null) {
+				try {
+					irrigWriter.close();
+				} 
+				catch (IOException e) {
+					e.printStackTrace();
+				}
+			}
+		}
+	}
+
+	private boolean isSpringWheat(RasterKey key) {
+		boolean isSpringWheat = false;
+		YieldResponsesItem yieldItem = yieldSurfaces.get(key);
+		if (yieldItem != null)
+			isSpringWheat = yieldSurfaces.get(key).isSpringWheat();
+		else 
+			LogWriter.printlnError("LpjgOutputer: Can't find YieldResponsesItem for " + key);
+		return isSpringWheat;
 	}
 
 	public static void writeMarkerFile(int year, boolean errorStatus) {
@@ -103,12 +176,7 @@ public class LpjgOutputer {
 				landCoverWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, year, crop/area, pasture/area, (forest+otherNatural)/area, barren/area));
 				landCoverWriter.newLine();
 
-				boolean isSpringWheat = false;
-				YieldResponsesItem yieldItem = yieldSurfaces.get(key);
-				if (yieldItem != null)
-					isSpringWheat = yieldSurfaces.get(key).isSpringWheat();
-				else 
-					LogWriter.printlnError("LpjgOutputer: Can't find YieldResponsesItem for " + key);
+				boolean isSpringWheat = isSpringWheat(key);
 
 				double lpjWinterWheatFrac = item.getCropFraction(CropType.OILCROPS, isSpringWheat ? null : CropType.WHEAT);
 				double lpjSpringWheatFrac = item.getCropFraction(CropType.PULSES, CropType.STARCHY_ROOTS, isSpringWheat ? CropType.WHEAT : null);
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index fb5d8d40..d71b2d7d 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/src/ac/ed/lurg/yield/YieldResponse.java
@@ -3,6 +3,7 @@ package ac.ed.lurg.yield;
 import java.util.HashMap;
 import java.util.Map;
 
+import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.types.YieldType;
 import ac.ed.lurg.utils.LogWriter;
 
@@ -31,7 +32,9 @@ public class YieldResponse {
 	public double getFertParam() {
 		if (fertParm == 0) {
 			//	fertParm = calcParam (0, 0.5, 1, 5.0/200, 50.0/200, 200.0/200); // we do have MID fert data, but looks wrong 
-			fertParm = calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG), 5.0/200, 50.0/200, 200.0/200);
+			fertParm = calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG), 
+					ModelConfig.MIN_FERT_AMOUNT/ModelConfig.MAX_FERT_AMOUNT, ModelConfig.MID_FERT_AMOUNT/ModelConfig.MAX_FERT_AMOUNT, ModelConfig.MAX_FERT_AMOUNT/ModelConfig.MAX_FERT_AMOUNT);
+			
 			if (fertParm > 1)
 				LogWriter.println(String.format("%s %s %s", getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG)));
 		}
-- 
GitLab