From 0658b14fd0faa538f204e9364fb314ae56a9cb53 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Thu, 17 Sep 2015 13:07:20 +0100
Subject: [PATCH] Calc irrigation usage more correctly
---
GAMS/IntExtOpt.gms | 20 ++++++++++---------
src/ac/ed/lurg/ModelConfig.java | 3 ++-
src/ac/ed/lurg/ModelMain.java | 4 ++--
.../country/gams/GamsRasterOptimiser.java | 9 ++++++++-
src/ac/ed/lurg/landuse/Intensity.java | 1 +
src/ac/ed/lurg/landuse/IrrigationItem.java | 2 +-
src/ac/ed/lurg/landuse/LandUseItem.java | 15 ++++++++------
src/ac/ed/lurg/utils/LogWriter.java | 2 +-
8 files changed, 35 insertions(+), 21 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index db113c5c..4defc99f 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -23,9 +23,9 @@
PARAMETER worldImportPrices(import_crop) prices for imports;
PARAMETER maxNetImport(import_crop) maximum net import for each crop based on world market;
PARAMETER minNetImport(import_crop) minimum net import for each crop based on world market;
- PARAMETER irrigCost(location) irrigation cost;
- PARAMETER irrigMaxRate(crop, location) max water application rate irrigation in mm;
- PARAMETER irrigConstraint(location) max water available for irrigation in mm;
+ PARAMETER irrigCost(location) irrigation cost in energy per 1000 Mlitre or Mha for each litre per m2;
+ PARAMETER irrigMaxRate(crop, location) max water application rate irrigation in litre per m2;
+ PARAMETER irrigConstraint(location) max water available for irrigation in litre per m2;
PARAMETER cropAdj(crop) this is the fudge factor that allows the model to fit the observed data
SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
@@ -112,7 +112,7 @@ $gdxin
UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= ( baseCost(crop) +
fertiliserUnitCost * fertI(crop, location) +
- irrigCost(location) * irrigI(crop, location) +
+ irrigCost(location) * irrigMaxRate(crop, location) * irrigI(crop, location) +
SQR(otherIntensity(crop, location)) * 1
) ;
@@ -166,11 +166,13 @@ $gdxin
ENERGY_EQ .. energy =E=
(
SUM((crop, location), area(crop, location) * unitEnergy(crop, location)) +
- sum(location, agriLandExpansion(location) * landChangeEnergy) +
- sum(location, (SQR(cropIncrease(location) + cropDecrease(location)) / (suitableLandArea(location) * 0.01 + sum(crop_less_pasture, previousArea(crop_less_pasture, location)))
- + (cropIncrease(location) + cropDecrease(location)) * 0.1) * landChangeEnergy) +
- sum(location, (SQR(pastureIncrease(location) + pastureDecrease(location)) / (suitableLandArea(location) * 0.01 + previousArea('pasture', location))
- + (pastureIncrease(location) + pastureDecrease(location)) * 0.1) * landChangeEnergy) +
+ sum(location,
+ agriLandExpansion(location) +
+ 0.45 * SQR(cropIncrease(location) + cropDecrease(location)) / (suitableLandArea(location) * 0.01 + sum(crop_less_pasture, previousArea(crop_less_pasture, location))) +
+ 0.05 * (cropIncrease(location) + cropDecrease(location)) +
+ 0.45 * SQR(pastureIncrease(location) + pastureDecrease(location)) / (suitableLandArea(location) * 0.01 + previousArea('pasture', location)) +
+ 0.05 * (pastureIncrease(location) + pastureDecrease(location))
+ ) * landChangeEnergy +
sum(import_crop, importAmount(import_crop) * worldImportPrices(import_crop) - exportAmount(import_crop) * worldExportPrices(import_crop))
) / 1000000;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index c507777c..59019ffa 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -148,7 +148,7 @@ public class ModelConfig {
public static final double MARKET_LAMBA = getDoubleProperty("MARKET_LAMBA", 0.3); // controls international market price adjustment rate
public static final double POPULATION_AGGREG_LIMIT = getDoubleProperty("POPULATION_AGGREG_LIMIT", 40.0); // in millions, smaller countries are aggregated on a regional basis
- public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 2.0);
+ public static final double IRRIG_COST_SCALE_FACTOR = getDoubleProperty("IRRIG_COST_SCALE_FACTOR", 0.1);
public static final double FERTILISER_MAX_COST = getDoubleProperty("FERTILISER_MAX_COST", 2.5);
public static final double TRANSPORT_LOSSES = getDoubleProperty("TRANSPORT_LOSSES", 0.15); // in international trade
public static final double TRADE_BARRIER_FACTOR = getDoubleProperty("TRADE_BARRIER_FACTOR", 1.2); // price factor in international trade, transport cost and real trade barriers
@@ -157,4 +157,5 @@ public class ModelConfig {
public static final int NUM_PASTURE_CATEGORIES = getIntProperty("NUM_PASTURE_CATEGORIES", 1);
public static final boolean DEBUG_LIMIT_COUNTRIES = getBooleanProperty("DEBUG_LIMIT_COUNTRIES", false);
+ public static final double DEFAULT_MAX_IRRIGATION_RATE = getDoubleProperty("DEFAULT_MAX_IRRIGATION_RATE", 50.0); // should need this but some areas crops don't have a value, but was causing them to be selected
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 5f9982f5..2f25d8d8 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -190,7 +190,7 @@ public class ModelMain {
private void writeLandCoverFile(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {
try {
- StringBuffer sbHeadings = new StringBuffer("Year, Cropland (Mha), Pasture (Mha), Natural (Mha), Fert crop (Mt), Fert pasture (Mt), Irrig crop (M litre), Irrig pasture (M litre)");
+ StringBuffer sbHeadings = new StringBuffer("Year, Cropland (Mha), Pasture (Mha), Natural (Mha), Fert crop (Mt), Fert pasture (Mt), Irrig crop (km3), Irrig pasture (km3)");
BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.LAND_COVER_OUTPUT_FILE, sbHeadings.toString());
StringBuffer sbData = new StringBuffer();
@@ -303,7 +303,7 @@ public class ModelMain {
// DEBUG code
if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {
- if (!(cc.getName().equals("United States of America") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
+ if (!(cc.getName().equals("United States of Americaxx") || cc.getName().equals("India") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
continue;
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 105e2978..fea86a05 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -355,7 +355,14 @@ public class GamsRasterOptimiser {
// Crops yields and area fractions
for (CropType crop : CropType.getNonMeatTypes()) {
if (irrigItem!= null) {
- aggIrig.setMaxIrrigAmount(crop, aggregateMean(aggIrig.getMaxIrrigAmount(crop), suitableAreaSoFar, irrigItem.getMaxIrrigAmount(crop), suitableAreaThisTime));
+ double irrigMax = irrigItem.getMaxIrrigAmount(crop);
+
+ if (irrigMax == 0.0) {
+ LogWriter.printlnError("Can't find irrig max amount for col:" + key.getCol() + ", row:" + key.getRow() + ", x:" + yieldRaster.getXCoordin(key) + ", y:" + yieldRaster.getYCoordin(key) + ", crop:" + crop);
+ irrigMax = ModelConfig.DEFAULT_MAX_IRRIGATION_RATE;
+ }
+
+ aggIrig.setMaxIrrigAmount(crop, aggregateMean(aggIrig.getMaxIrrigAmount(crop), suitableAreaSoFar, irrigMax, suitableAreaThisTime));
}
double areaSoFar = aggLandUse.getCropArea(crop);
diff --git a/src/ac/ed/lurg/landuse/Intensity.java b/src/ac/ed/lurg/landuse/Intensity.java
index 23f62525..79bdedf5 100644
--- a/src/ac/ed/lurg/landuse/Intensity.java
+++ b/src/ac/ed/lurg/landuse/Intensity.java
@@ -45,6 +45,7 @@ public class Intensity {
return unitEnergy;
}
+ /** irrigation rate in litre/m2 */
public double getIrrigationRate() {
return maxIrrigRate * irrigationIntensity;
}
diff --git a/src/ac/ed/lurg/landuse/IrrigationItem.java b/src/ac/ed/lurg/landuse/IrrigationItem.java
index b209ed0f..fa266640 100644
--- a/src/ac/ed/lurg/landuse/IrrigationItem.java
+++ b/src/ac/ed/lurg/landuse/IrrigationItem.java
@@ -11,7 +11,7 @@ public class IrrigationItem implements RasterItem {
// when is the limited to be >= 0 (only Greenland had negative values)
// this gives a value from 1 to 0, with 1 being the most costly and 0 the least.
- private double irrigationCost; // cost per ? of irrigation
+ private double irrigationCost; // cost per 10000 Mlitre or 10^10 l or 10^7 m3 or 10^-2 km3, when multipled by rates from areas in Mha makes sense to have it in this apparently odd unit
private double irrigConstraint; // mm or litre/m2 (they are the same) of water available
private Map<CropType, Double> maxIrrigAmounts = new HashMap<CropType, Double>(); // amount of water (litres/m2) used at maximum irrigation, i.e. to achieve no water stress
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index 7838ff59..94c2d60d 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -69,32 +69,35 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem> {
return getFertiliserTotal(items, crops.toArray(new CropType[crops.size()]));
}
- /** Irrigation rate in litre/ha */
+ /** Irrigation rate in litre/m2 */
public double getIrrigationRate(CropType crop) {
Intensity i = getIntensity(crop);
return (i == null) ? 0 : i.getIrrigationRate();
}
- /** Irrigation amount in M litres, for this location */
- public double getIrrigationAmount(CropType c) {
+ /** Irrigation amount in km3 or 10^9 m3, for this location */
+ private double getIrrigationAmount(CropType c) {
double rate = getIrrigationRate(c);
double area = getCropArea(c);
- return rate * area;
+ return rate * area * 0.01; // rate(10^-3m or mm or l/m2) * area(10^6ha) = 10^3 m.ha = 10^7 m3 = 1/100 km3
}
+ /** Irrigation rate in litre/m2 */
public double getIrrigationAverageRate(CropType... crops) {
double irrigTotal = 0;
double areaTotal = 0;
for (CropType c : crops) {
- irrigTotal += getIrrigationAmount(c);
- areaTotal = getCropArea(c);
+ double area = getCropArea(c);
+ irrigTotal += getIrrigationRate(c) * area;
+ areaTotal += area;
}
return areaTotal > 0 ? irrigTotal / areaTotal : 0;
}
+ /** Irrigation amount in km3, for this location */
public static double getIrrigationTotal(Collection<? extends LandUseItem> items, CropType... crops) {
double total = 0;
for (LandUseItem a : items) {
diff --git a/src/ac/ed/lurg/utils/LogWriter.java b/src/ac/ed/lurg/utils/LogWriter.java
index 27ca8df9..32b5da0c 100644
--- a/src/ac/ed/lurg/utils/LogWriter.java
+++ b/src/ac/ed/lurg/utils/LogWriter.java
@@ -52,7 +52,7 @@ public class LogWriter {
private void print(String s, PrintStream stream, boolean newLine) {
logFile.write(s);
if (newLine) logFile.write(System.lineSeparator());
- // logFile.flush();
+ logFile.flush();
if (!ModelConfig.SUPPRESS_STD_OUTPUT) {
stream.print(s);
--
GitLab