From b1a6b46091c2ba6465141970e477dd3788901b64 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Mon, 14 Sep 2015 10:30:32 +0100
Subject: [PATCH] Change to exponential fertiliser response
---
GAMS/IntExtOpt.gms | 12 +++++++---
src/ac/ed/lurg/ModelConfig.java | 2 +-
src/ac/ed/lurg/ModelMain.java | 7 +++---
.../country/gams/GamsRasterOptimiser.java | 6 ++++-
.../lurg/yield/LPJYieldResponseMapReader.java | 8 ++++---
src/ac/ed/lurg/yield/YieldResponse.java | 22 +++++++++++++------
src/ac/ed/lurg/yield/YieldResponsesItem.java | 4 ----
.../raster/AbstractTabularRasterReader.java | 2 +-
8 files changed, 40 insertions(+), 23 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 0ac433f0..de122d72 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -119,10 +119,16 @@ $gdxin
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yieldNone(crop, location) +
- (yieldFertOnly(crop, location) - yieldNone(crop, location)) * ((fertI(crop, location) + delta) ** fertParam(crop, location)) +
- (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) +
+* (yieldFertOnly(crop, location) - yieldNone(crop, location)) * ((fertI(crop, location) + delta) ** fertParam(crop, location)) +
+* (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) +
+* (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
+* ((fertI(crop, location) + delta) ** fertParam(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location))
+
+ (yieldFertOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-fertParam(crop, location)*fertI(crop, location))) +
+ (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location))) +
(yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
- ((fertI(crop, location) + delta) ** fertParam(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location))
+ (1 - exp(-fertParam(crop, location)*fertI(crop, location))) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location)))
+
) * cropAdj(crop) * otherIntensity(crop, location);
NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + netImportAmount(crop);
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 2e1a43b0..1e4673a8 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -128,7 +128,7 @@ public class ModelConfig {
// 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;
+ public static final double MAX_FERT_AMOUNT = 500;
// Other model parameters
public static final boolean KEEP_DEMAND_FIXED = getBooleanProperty("KEEP_DEMAND_FIXED", false);
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index a6827925..b683a6dd 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -20,13 +20,13 @@ import ac.ed.lurg.demand.BaseConsumpManager;
import ac.ed.lurg.demand.DemandManager;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.CropUsageReader;
-import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.IrrigationConstraintReader;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.IrrigationMaxAmountReader;
import ac.ed.lurg.landuse.IrrigiationCostReader;
import ac.ed.lurg.landuse.LandCoverItem;
import ac.ed.lurg.landuse.LandCoverReader;
+import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.output.LpjgOutputer;
import ac.ed.lurg.types.CropToDoubleMap;
import ac.ed.lurg.types.CropType;
@@ -35,6 +35,7 @@ import ac.ed.lurg.types.ModelFitType;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.LPJYieldResponseMapReader;
import ac.ed.lurg.yield.YieldRaster;
+import ac.ed.lurg.yield.YieldResponsesItem;
import ac.sac.raster.IntegerRasterItem;
import ac.sac.raster.InterpolatingRasterSet;
import ac.sac.raster.RasterHeaderDetails;
@@ -97,8 +98,8 @@ public class ModelMain {
YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
- // YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-50.0, -4.0);
- // LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.CEREALS) + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+ YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-90.5, 45.5);
+ LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.MAIZE) + ", " + yresp.getYieldFertOnly(CropType.MAIZE) + ", " + yresp.getYieldIrrigOnly(CropType.MAIZE) + ", " + yresp.getYieldNone(CropType.MAIZE));
CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
CropToDoubleMap totalExportCommodities = new CropToDoubleMap();
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 6df3f1b8..05d9d376 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -357,7 +357,11 @@ public class GamsRasterOptimiser {
// Crops yields and area fractions
for (CropType crop : CropType.getNonMeatTypes()) {
for (YieldType yieldType : YieldType.values()) {
- avgYResp.setYield(yieldType, crop, yresp.getYield(yieldType, crop));
+ double y = yresp.getYield(yieldType, crop);
+ if (Double.isNaN(y))
+ LogWriter.printlnError("Problem getting yield at, col:" + key.getCol() + ", row:" + key.getRow() + ", x:" + yieldRaster.getXCoordin(key) + ", y:" + yieldRaster.getYCoordin(key) + ", type:" + yieldType + ", crop:" + crop);
+ else
+ avgYResp.setYield(yieldType, crop, y);
}
double areaSoFar = aggLandUse.getCropArea(crop);
diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
index 74231fca..3eb218e7 100644
--- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
+++ b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
@@ -81,10 +81,12 @@ public class LPJYieldResponseMapReader extends AbstractTabularRasterReader<Yield
for (IrrigationRate irrig : IrrigationRate.values()) {
YieldType yieldType = YieldType.getYieldType(fert, irrig);
- if (!fert.equals(FertiliserRate.NO_FERT))
- adjFactor = 10; // 10 for kg/m2 to t/ha
+ if (fert.equals(FertiliserRate.NO_FERT))
+ adjFactor = 10 * 0.4; // adjust no fert yield to 40% of figure
+ else if (fert.equals(FertiliserRate.MID_FERT))
+ adjFactor = 10 * 0.9; // adjust mid fert yield to 90% of figure
else
- adjFactor = 5;
+ adjFactor = 10; // 10 for kg/m2 to t/ha
String fertIrrigString = irrig.getId() + fert.getId();
double ww = getValueForCol(rowValues, "TeWW" + fertIrrigString) * adjFactor;
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index d71b2d7d..590031f7 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/src/ac/ed/lurg/yield/YieldResponse.java
@@ -5,13 +5,19 @@ import java.util.Map;
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.types.YieldType;
-import ac.ed.lurg.utils.LogWriter;
public class YieldResponse {
Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
private double fertParm;
private double irrigParm;
+ private static double defaultParam;
+
+ static {
+ defaultParam = calcParam(0, 0.8, 1, 0, 0.5, 1.0); // default assumes 80% at mid point
+ }
+
+
public void setYield(YieldType type, double yield) {
yields.put(type, yield);
fertParm = 0;
@@ -35,23 +41,25 @@ public class YieldResponse {
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)));
+ //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)));
}
return fertParm;
}
public double getIrrigParam() {
if (irrigParm == 0) {
- irrigParm = calcParam (0, 0.7, 1, 0, 0.5, 1.0); // we don't have a mid irrigation figure, so lets assume 60% at mid point
+ irrigParm = defaultParam;
}
+
return irrigParm;
}
- private double calcParam(double yMin, double yMid, double yMax, double xMin, double xMid, double xMax) {
+ private static double calcParam(double yMin, double yMid, double yMax, double xMin, double xMid, double xMax) {
if (yMid <= yMin || yMax <= yMid)
- return 1; // default to linear
+ return defaultParam;
- return Math.log((yMid - yMin)/(yMax - yMin)) / Math.log((xMid - xMin)/(xMax - xMin)) / xMax;
+ return -Math.log(1 - ((yMid-yMin)/(yMax-yMin))) * (xMax-xMin) / (xMid-xMin); // 1-exp(-x) function
+ // return Math.log((yMid - yMin)/(yMax - yMin)) / Math.log((xMid - xMin)/(xMax - xMin)) / xMax; // Cobb-Douglas
}
+
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/yield/YieldResponsesItem.java b/src/ac/ed/lurg/yield/YieldResponsesItem.java
index 159d5554..d90f6ade 100644
--- a/src/ac/ed/lurg/yield/YieldResponsesItem.java
+++ b/src/ac/ed/lurg/yield/YieldResponsesItem.java
@@ -5,7 +5,6 @@ import java.util.Map;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.YieldType;
-import ac.ed.lurg.utils.LogWriter;
import ac.sac.raster.RasterItem;
public class YieldResponsesItem implements RasterItem {
@@ -28,9 +27,6 @@ public class YieldResponsesItem implements RasterItem {
public double getYield(YieldType yieldType, CropType crop) {
YieldResponse yr = getYieldResponseForCrop(crop);
double d = yr.getYield(yieldType);
- if (Double.isNaN(d)) {
- LogWriter.printlnError("YieldResponses: No yield found for type " + yieldType + ", crop " + crop);
- }
return d;
}
diff --git a/src/ac/sac/raster/AbstractTabularRasterReader.java b/src/ac/sac/raster/AbstractTabularRasterReader.java
index dd3fdb47..1cabc85d 100644
--- a/src/ac/sac/raster/AbstractTabularRasterReader.java
+++ b/src/ac/sac/raster/AbstractTabularRasterReader.java
@@ -40,7 +40,7 @@ public abstract class AbstractTabularRasterReader<D extends RasterItem> {
// List<String> colNames = new ArrayList<String>(Arrays.asList(headertokens));
for (int i=0; i<headertokens.length; i++) {
- headertokens[i] = headertokens[i].toLowerCase();
+ headertokens[i] = headertokens[i].toLowerCase().replace("\"", "");
}
return headertokens;
--
GitLab