From dd1a2063d6c1f207cecf091a2f91aec1a7ba235c Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Tue, 7 Jul 2015 12:28:10 +0100
Subject: [PATCH] Changes to fix reading of LPJG file
---
src/ac/ed/lurg/ModelMain.java | 4 +
.../country/gams/GamsRasterOptimiser.java | 107 +++++++++++-------
.../lurg/yield/LPJYieldResponseMapReader.java | 31 +++--
src/ac/ed/lurg/yield/YieldResponse.java | 4 +-
src/ac/sac/raster/RasterKey.java | 2 +-
src/ac/sac/raster/RasterSet.java | 5 +-
src/ac/sac/raster/RasterSetTest.java | 19 ++++
7 files changed, 118 insertions(+), 54 deletions(-)
create mode 100644 src/ac/sac/raster/RasterSetTest.java
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 6c28770b..c937c928 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -67,6 +67,10 @@ public class ModelMain {
private void doTimestep(int timestep) {
YieldRaster yieldSurfaces = getYieldSurfaces(timestep);
+
+// YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-50.0, -4.0);
+// LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.CEREALS) + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+
LogWriter.println("Timestep " + timestep);
CropToDoubleMap totalArea = new CropToDoubleMap();
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 1113aa8f..55dce778 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -25,7 +25,7 @@ import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterSet;
public class GamsRasterOptimiser {
- public static final boolean DEBUG = false;
+ public static final boolean DEBUG = true;
private GamsRasterInput rasterInputData;
private LazyHashMap<Integer, Set<RasterKey>> mapping;
@@ -49,28 +49,28 @@ public class GamsRasterOptimiser {
private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) {
RasterSet<AreasItem> newAreaRaster = allocAreas(gamsInput.getPreviousAreas(), gamsOutput);
RasterSet<IntensitiesItem> newIntensityRaster = allocIntensities(gamsOutput);
-
+
return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData());
}
private RasterSet<AreasItem> allocAreas(Map<Integer, ? extends AreasItem> prevAreasAgg, GamsLocationOutput gamsOutput) {
RasterSet<AreasItem> newAreaRaster = new RasterSet<AreasItem>();
RasterSet<AreasItem> prevAreaRaster = rasterInputData.getPreviousAreas();
-
+
for (Map.Entry<Integer, AreasItem> entry : gamsOutput.getCropAreas().entrySet()) {
Integer locId = entry.getKey();
AreasItem newAreaAggItem = entry.getValue();
AreasItem prevAreaAggItem = prevAreasAgg.get(locId);
Set<RasterKey> keys = mapping.get(locId);
-
+
CropToDouble cropChangeTotals = newAreaAggItem.getCropChanges(prevAreaAggItem);
- if (DEBUG) LogWriter.println("Processing location id " + locId);
+ // if (DEBUG) LogWriter.println("Processing location id " + locId);
double shortfall = allocAllAreasForAgg(newAreaRaster, prevAreaRaster, keys, cropChangeTotals);
if (shortfall > 0)
LogWriter.printlnError("This should never happen, due to GAMS constraint. Not able to incorporate all changes, as not enough forest or natural areas left: " + shortfall);
}
-
+
return newAreaRaster;
}
@@ -79,9 +79,9 @@ public class GamsRasterOptimiser {
double totalShortfall = 0;
Set<RasterKey> keysWithSpace = new HashSet<RasterKey>();
CropToDouble avgCropChange = cropChange.scale(1.0/keys.size());
-
+
for (RasterKey key : keys) {
- if (DEBUG) LogWriter.println(" Processing raster key " + key);
+ //if (DEBUG) LogWriter.println(" Processing raster key " + key);
AreasItem prevAreasRasterItem = prevAreaRaster.get(key);
AreasItem newAreasRasterItem = new AreasItem();
@@ -97,13 +97,13 @@ public class GamsRasterOptimiser {
if (totalShortfall > 0 & keysWithSpace.size() > 0) { // more to allocate and some free areas to allocate into
if (DEBUG) LogWriter.println("A total area of " + totalShortfall + " is unallocated, so trying again with areas with space");
double factor = totalShortfall / cropChange.getTotal();
-
+
allocAllAreasForAgg(newAreaRaster, newAreaRaster, keysWithSpace, cropChange.scale(factor));
}
return totalShortfall;
}
-
+
/** returns the net area of crop change that not been able to be allocated */
private double allocRasterArea(AreasItem newAreasRasterItem, AreasItem prevAreasRasterItem, CropToDouble cropChange) {
@@ -114,7 +114,7 @@ public class GamsRasterOptimiser {
double prevNatural = prevAreasRasterItem.getLandCoverArea(LandDataType.OTHER_NATURAL);
double ratioOfChanges = 1;
-
+
if (netAddedCrop > prevForest + prevNatural) {
ratioOfChanges = (prevForest + prevNatural) / netAddedCrop;
forestRemoved = prevForest;
@@ -139,9 +139,9 @@ public class GamsRasterOptimiser {
for (CropType crop : CropType.values()) {
double prevCropArea = prevAreasRasterItem.getCropArea(crop);
Double d = cropChange.get(crop);
-
+
double additionalArea = (d == null) ? 0 : d * ratioOfChanges;
- // LogWriter.println(String.format("%s, %.1f", crop, additionalArea));
+ // LogWriter.println(String.format("%s, %.1f", crop, additionalArea));
newAreasRasterItem.setCropArea(crop, prevCropArea + additionalArea);
}
@@ -171,11 +171,32 @@ public class GamsRasterOptimiser {
RasterSet<AreasItem> cropAreaRaster = rasterInputData.getPreviousAreas();
RasterSet<IrrigationCostItem> irrigCostRaster = rasterInputData.getIrrigationCost();
- List<Double> cerealDivisions = getDivisions(yieldRaster, CropType.CEREALS, 2);
- List<Double> pastureDivisions = getDivisions(yieldRaster, CropType.MEAT_OR_PASTURE, 2);
+ {
+ // YieldResponsesItem yresp = yieldRaster.getFromCoordinates(-118.0, -33.0);
+ // LogWriter.printlnError("Test key2: " + yresp.getYieldMax(CropType.CEREALS) + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+ }
+
+ // look for inconsistencies
+ for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
+
+ YieldResponsesItem yresp = entry.getValue();
+ RasterKey key = entry.getKey();
+
+ if (yresp != null && (yresp.getYieldMax(CropType.CEREALS) < yresp.getYieldFertOnly(CropType.CEREALS) || yresp.getYieldMax(CropType.CEREALS) < yresp.getYieldIrrigOnly(CropType.CEREALS))) {
+ LogWriter.printlnError(key.getCol() + ", " + key.getRow() + ": " + yieldRaster.getXCoordin(key) + ", " + yieldRaster.getYCoordin(key) + ": " + yresp.getYieldMax(CropType.CEREALS) + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+ }
+ // else if (yresp != null)
+ // LogWriter.println(key.getCol() + ", " + key.getRow() + ": " + yieldRaster.getXCoordin(key) + ", " + yieldRaster.getYCoordin(key) + ": " + yresp.getYieldMax(CropType.CEREALS) + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));
+
+ }
- int numCerealCats = (1+cerealDivisions.size());
- if (DEBUG) LogWriter.println("Making " + numCerealCats * (1+pastureDivisions.size()) + " categories");
+ int numCerealCats = 4;
+ int numPastureCats = 4;
+
+ List<Double> cerealDivisions = getDivisions(yieldRaster, CropType.CEREALS, numCerealCats);
+ List<Double> pastureDivisions = getDivisions(yieldRaster, CropType.MEAT_OR_PASTURE, numPastureCats);
+
+ if (DEBUG) LogWriter.println("Making " + numCerealCats * numPastureCats + " categories");
LazyHashMap<Integer, AveragingYieldResponsesItem> aggregatedYields = new LazyHashMap<Integer, AveragingYieldResponsesItem>() {
protected AveragingYieldResponsesItem createValue() { return new AveragingYieldResponsesItem(); }
@@ -190,14 +211,18 @@ public class GamsRasterOptimiser {
protected Set<RasterKey> createValue() { return new HashSet<RasterKey>(); }
};
+ int countFound = 0, countMissing = 0;
+
for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
YieldResponsesItem yresp = entry.getValue();
RasterKey key = entry.getKey();
-
+
if (yresp == null) {
+ countMissing++;
LogWriter.printlnError("No YieldResponsesItem for " + key + ", x=" + yieldRaster.getXCoordin(key) + ", y=" + yieldRaster.getYCoordin(key));
}
else {
+ countFound++;
//LogWriter.println("YieldResponsesItem found for " + key + ", x=" + yieldRaster.getXCoordin(key) + ", y=" + yieldRaster.getYCoordin(key));
AreasItem cropAreas = cropAreaRaster.get(key);
@@ -205,18 +230,18 @@ public class GamsRasterOptimiser {
LogWriter.printlnError("GamsRasterOptimiser: Can't find cropAreas for " + key);
continue;
}
-
+
IrrigationCostItem irrigCost = irrigCostRaster.get(key);
-
+
int cerealCat = findCategory(cerealDivisions, yresp.getYieldMax(CropType.CEREALS));
int pastureCat = findCategory(pastureDivisions, yresp.getYieldMax(CropType.MEAT_OR_PASTURE));
Integer id = cerealCat + pastureCat * numCerealCats;
-
+
AveragingYieldResponsesItem avgYResp = aggregatedYields.lazyGet(id);
AreasItem aggAreas = aggregatedAreas.lazyGet(id);
IrrigationCostItem avgIrigCost = aggregatedIrrigCosts.lazyGet(id);
mapping.lazyGet(id).add(entry.getKey());
-
+
// Irrigation cost
double suitableAreaThisTime = cropAreas.getSuitableLand();
if (irrigCost!= null && suitableAreaThisTime > 0) {
@@ -230,12 +255,12 @@ public class GamsRasterOptimiser {
for (YieldType yieldType : YieldType.values()) {
avgYResp.setYield(yieldType, crop, yresp.getYield(yieldType, crop));
}
-
+
double areaSoFar = aggAreas.getCropArea(crop);
double areaThisTime = cropAreas.getCropArea(crop);
aggAreas.setCropArea(crop, areaSoFar + areaThisTime);
}
-
+
// Land cover, we only need forest and other natural
for (LandDataType landType : Arrays.asList(LandDataType.FOREST, LandDataType.OTHER_NATURAL)) {
double areaThisTime = cropAreas.getLandCoverArea(landType);
@@ -244,13 +269,14 @@ public class GamsRasterOptimiser {
}
}
}
-
- if (DEBUG) {
- for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) {
- LogWriter.println(e.getKey() + " category has " + e.getValue().size() + " raster areas");
- for (RasterKey key : e.getValue()) {
- LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.MEAT_OR_PASTURE), yieldRaster.get(key).getYieldMax(CropType.CEREALS)));
- }
+
+ LogWriter.println("YieldResponsesItem: " + rasterInputData.getCountryInput().getCountry() + ", countFound=" + countFound + ", countMissing=" + countMissing);
+
+
+ for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) {
+ LogWriter.println(e.getKey() + " category has " + e.getValue().size() + " raster areas");
+ if (DEBUG) for (RasterKey key : e.getValue()) {
+ LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.MEAT_OR_PASTURE), yieldRaster.get(key).getYieldMax(CropType.CEREALS)));
}
}
@@ -262,7 +288,7 @@ public class GamsRasterOptimiser {
int numDivisions = divisions.size();
for (category = 0; category<numDivisions; category++) {
- if (yield > divisions.get(category)) {
+ if (yield < divisions.get(category)) {
break;
}
}
@@ -274,28 +300,31 @@ public class GamsRasterOptimiser {
for (YieldResponsesItem yresp : yieldRaster.values()) {
if (yresp == null) {
- //LogWriter.println("GamsRasterOptimiser: Can't get value for crop " + crop);
+ if (DEBUG) LogWriter.println("GamsRasterOptimiser: Can't get value for crop " + crop);
}
else {
double d = yresp.getYieldMax(crop);
-
- if (Double.isNaN(d))
- LogWriter.println("GamsRasterOptimiser: Got NaN value for crop " + crop);
- else
+ //LogWriter.println("GamsRasterOptimiser: Got value for crop " + crop + " of " + d);
+
+ if (Double.isNaN(d) || d == 0.0) {
+ if (DEBUG) LogWriter.println("GamsRasterOptimiser: Got NaN or zero value for crop " + crop);
+ }
+ else {
yieldValues.add(d);
+ }
}
}
if (yieldValues.size() == 0) {
throw new RuntimeException("No yield values for country, crop = " + crop);
}
-
+
Collections.sort(yieldValues);
List<Double> divisions = new ArrayList<Double>();
for (int i=1; i < numCategories; i++) {
- double d = yieldValues.get(yieldValues.size()/numCategories);
+ double d = yieldValues.get(yieldValues.size()/numCategories*i);
divisions.add(d);
}
diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
index 1f21ede6..343b7abe 100644
--- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
+++ b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
@@ -55,21 +55,29 @@ public class LPJYieldResponseMapReader {
record.y = d;
break;
case 2:
- record.teWW = d;
+ record.teSW = d;
break;
case 3:
- record.trRi = d;
+ record.teSWirr = d;
break;
case 4:
- record.teCo = d;
+ record.teWW = d;
break;
case 5:
- record.teSo = d;
+ record.teWWirr = d;
break;
case 6:
- record.teWWirr = d;
+ record.teCo = d;
+ break;
+ case 7:
+ record.teCoirr = d;
+ break;
+ case 8:
+ record.trRi = d;
+ break;
+ case 9:
+ record.trRiirr = d;
break;
-
}
col++;
@@ -99,18 +107,21 @@ public class LPJYieldResponseMapReader {
//YieldRecord values in kg DM / m2, we want t/ha, so 10 times larger
for (CropType crop : CropType.getAllItems()) {
- item.setYield(noIrrigYieldType, crop, record.teWW * 10);
- item.setYield(maxIrrigYieldType, crop, record.teWWirr * 10);
+ item.setYield(noIrrigYieldType, crop, record.teSW * 10);
+ item.setYield(maxIrrigYieldType, crop, record.teSWirr * 10);
}
}
class YieldRecord {
double x;
double y;
+ double teSW;
+ double teSWirr;
double teWW;
double teWWirr;
- double trRi;
double teCo;
- double teSo;
+ double teCoirr;
+ double trRi;
+ double trRiirr;
}
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index 872cbbac..5aeb8d98 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/src/ac/ed/lurg/yield/YieldResponse.java
@@ -7,7 +7,7 @@ import ac.ed.lurg.types.YieldType;
import ac.ed.lurg.utils.LogWriter;
public class YieldResponse {
- private Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
+ Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
private double fertParm;
private double irrigParm;
@@ -37,7 +37,7 @@ public class YieldResponse {
public double getIrrigParam() {
if (irrigParm == 0)
- irrigParm = calcParam (0, 0.6, 1); // we don't have a mide irrigation figure, so lets assume 60% at mid point
+ irrigParm = calcParam (0, 0.6, 1); // we don't have a mid irrigation figure, so lets assume 60% at mid point
return irrigParm;
}
diff --git a/src/ac/sac/raster/RasterKey.java b/src/ac/sac/raster/RasterKey.java
index 4b0def60..678e1a12 100755
--- a/src/ac/sac/raster/RasterKey.java
+++ b/src/ac/sac/raster/RasterKey.java
@@ -77,6 +77,6 @@ public class RasterKey implements Serializable {
}
public String toString() {
- return "row:" + row + ", col:" + col;
+ return "col:" + col + ", row:" + row;
}
}
diff --git a/src/ac/sac/raster/RasterSet.java b/src/ac/sac/raster/RasterSet.java
index 9e2240aa..9a39a023 100755
--- a/src/ac/sac/raster/RasterSet.java
+++ b/src/ac/sac/raster/RasterSet.java
@@ -75,8 +75,8 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
}
public double getYCoordin(RasterKey key) {
- double x = header.getCellSize() * (key.getRow() + 1) + header.getYllCorner();
- return x;
+ double y = header.getCellSize() * (header.getNrows() - 1 - key.getRow()) + header.getYllCorner();
+ return y;
}
/** Method to really get the data, or create it. Assumes the col and row are in the internal header format */
@@ -123,6 +123,7 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
public RasterSet<D> popSubsetForKeys(RasterSet<D> subset, Collection<RasterKey> keys) {
for (RasterKey key : keys) {
+ //LogWriter.println("popSubsetForKeys: " + key.getCol() + ", " + key.getRow() + ": " + getXCoordin(key) + ", " + getYCoordin(key));
subset.put(key, get(key));
}
return subset;
diff --git a/src/ac/sac/raster/RasterSetTest.java b/src/ac/sac/raster/RasterSetTest.java
new file mode 100644
index 00000000..55ba8547
--- /dev/null
+++ b/src/ac/sac/raster/RasterSetTest.java
@@ -0,0 +1,19 @@
+package ac.sac.raster;
+
+import org.junit.Test;
+
+import ac.ed.lurg.yield.YieldRaster;
+import ac.ed.lurg.yield.YieldResponsesItem;
+
+public class RasterSetTest {
+
+ @Test
+ public void test() {
+ YieldRaster dataset = new YieldRaster(new RasterHeaderDetails(720,360,-180,-90,0.5,""));
+ double y = dataset.getYCoordin(new RasterKey(0,179));
+ System.out.println("y=" +y);
+
+
+ }
+
+}
--
GitLab