diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index d3cf16ca35792b711aed03b9ba6f87a7d8d1c09a..db5f3d5ff0cc7231f464108beee7e938bc358959 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -52,6 +52,7 @@
PARAMETER subsidyRate(all_types) rates of subsidy compared to costs;
PARAMETER previousLandCoverArea(land_cover, location) land cover area in Mha;
+ PARAMETER maxCroplandArea(location);
PARAMETER carbonCreditRate(land_cover_before, land_cover_after, location) potential carbon credits from LULUCF - tC-eq per ha;
PARAMETER conversionCost(land_cover_before, land_cover_after, location) cost of converting from one land cover to another - $1000 per ha;
PARAMETER woodYieldMax(location);
@@ -84,7 +85,7 @@ $load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousO
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxGrossLccRate, subsidyRate
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandFraction, seedAndWasteRate
-$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam, maxRotationIntensity, forestBaseCost, tradeAdjustmentCostRate
+$load previousLandCoverArea, maxCroplandArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam, maxRotationIntensity, forestBaseCost, tradeAdjustmentCostRate
$load forestManagementCost, vegClearingCostRate, carbonForestMaxProportion, maxFertChange, maxIrrigChange, previousRotationIntensity
$gdxin
@@ -286,7 +287,10 @@ $gdxin
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= previousLandCoverArea(land_cover, location);
CONVERSION_COST_CALC(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after),
- landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after, location));
+ landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after, location));
+
+* Cropland area restriction due to slope
+ landCoverArea.UP('cropland', location) = maxCroplandArea(location);
* MIN_LAND_COVER_AREA_CONSTRAINT(land_cover, location) .. landCoverArea(land_cover, location) =G= minLandCoverArea(land_cover, location);
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 86775fcc41f1bcf603f3e93a767b7e4c04e683e6..05f69af971860b8fc675044489b141d1b00a1d05 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -129,6 +129,7 @@ public class GamsLocationOptimiser {
GAMSParameter seedAndWasteRateP = inDB.addParameter("seedAndWasteRate", 1);
GAMSParameter prevLandCoverP = inDB.addParameter("previousLandCoverArea", 2);
+ GAMSParameter maxCroplandAreaP = inDB.addParameter("maxCroplandArea", 1);
GAMSParameter prevRotaIntP = inDB.addParameter("previousRotationIntensity", 1);
double totalAgriLand = 0;
@@ -197,7 +198,7 @@ public class GamsLocationOptimiser {
setGamsParamValueTruncate(prevLandCoverP.addRecord(v), area, 6);
}
- double foo = landUseItem.getForestRotationIntensity();
+ setGamsParamValue(maxCroplandAreaP.addRecord(locString), landUseItem.getMaxCroplandArea(), 6);
setGamsParamValue(prevRotaIntP.addRecord(locString), landUseItem.getForestRotationIntensity(), 6);
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index f671b3cddb380c01115aa508099d27890d79491c..019d50d9803c061f69949fa3f214b908ef4d85ab 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -161,6 +161,7 @@ public class GamsRasterOptimiser {
keys.add(key);
}
}
+
for (LandCoverChangeItem item : changes) {
LandCoverType fromLc = item.getFromLandCover();
LandCoverType toLc = item.getToLandCover();
@@ -170,12 +171,14 @@ public class GamsRasterOptimiser {
double totalArea = 0;
for (RasterKey key : keys) {
- totalArea += luRaster.get(key).getLandCoverArea(fromLc, LandProtectionType.CONVERTIBLE);
+ double area = luRaster.get(key).getLandCoverArea(fromLc, LandProtectionType.CONVERTIBLE);
+ totalArea += area;
}
for (RasterKey key : keys) {
double area = luRaster.get(key).getLandCoverArea(fromLc, LandProtectionType.CONVERTIBLE);
double fraction = area / totalArea;
+
// Splitting up the area changes in clusters proportionally by areas in raster cells
double changedArea = fraction * totalChangeArea;
if (changedArea < 2 * Double.MIN_VALUE) {
@@ -354,6 +357,9 @@ public class GamsRasterOptimiser {
aggLandUse.setLandCoverArea(lcType, LandProtectionType.PROTECTED, protAreaSoFar + protAreaThisTime);
}
+ // Maximum cropland area
+ aggLandUse.setMaxCroplandArea(aggLandUse.getMaxCroplandArea() + landUseItem.getMaxCroplandArea());
+
aggLandUse.setForestRotationIntensity(landUseItem.getForestRotationIntensity());
}
diff --git a/src/ac/ed/lurg/landuse/LandCoverItem.java b/src/ac/ed/lurg/landuse/LandCoverItem.java
index 861d8d93b53e64a0c644a717c3cc7dfffe230363..6638bd7f68ef603b93727d976be8f3be3e432a7e 100644
--- a/src/ac/ed/lurg/landuse/LandCoverItem.java
+++ b/src/ac/ed/lurg/landuse/LandCoverItem.java
@@ -20,10 +20,8 @@ public class LandCoverItem implements RasterItem {
private LandCoverTile landCoverAreas = new LandCoverTile();
private Map<CropType, Double> cropFractions = new HashMap<>();
private double totalArea;
- private double unavailableFract; // due to slope
+ private double maxCroplandFraction; // due to slope
private double protectedFraction;
- private List<LandCoverType> unavailPriorityList = new ArrayList<LandCoverType>(
- Arrays.asList(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE, LandCoverType.CROPLAND));
/** Area in Mha */
public Double getLandCoverArea(LandCoverType landType) {
@@ -53,12 +51,14 @@ public class LandCoverItem implements RasterItem {
landCoverFract.put(landType, d);
}
- public double getUnavailableFract(){
- return unavailableFract;
+ public double getMaxCroplandArea(){
+ double maxArea = maxCroplandFraction * getTotalArea();
+ double cropArea = landCoverAreas.getTotalArea(LandCoverType.CROPLAND);
+ return Math.max(maxArea, cropArea); // In case of disagreement, trust land cover data
}
public void setMaxCropFraction(double maxCropFraction){
- this.unavailableFract=1.0-maxCropFraction;
+ this.maxCroplandFraction = maxCropFraction;
}
public double getProtectedFraction() {
@@ -71,7 +71,6 @@ public class LandCoverItem implements RasterItem {
public LandCoverTile getLandCoverTiles() {
createLandCoverTiles();
- setUnavailableArea();
cleanUp();
return landCoverAreas;
}
@@ -126,40 +125,6 @@ public class LandCoverItem implements RasterItem {
runAreaCheck(initAreas, landCoverAreas);
}
- private void setUnavailableArea() {
- double alreadyUnavail = 0;
- for (LandCoverType lcType : LandCoverType.values()) {
- alreadyUnavail += landCoverAreas.getArea(lcType, LandProtectionType.UNAVAILABLE);
- }
- double unavailableArea = unavailableFract * totalArea;
- double shortfall = Math.max(0, unavailableArea - alreadyUnavail);
- if (shortfall == 0)
- return;
-
- for (LandCoverType lcType : unavailPriorityList) {
- double cArea = landCoverAreas.getArea(lcType, LandProtectionType.CONVERTIBLE);
- double pArea = landCoverAreas.getArea(lcType, LandProtectionType.PROTECTED);
- double totalArea = cArea + pArea;
- double cProp = cArea / totalArea;
-
- double additionalUnavail = Math.min(shortfall, totalArea);
- if (additionalUnavail > 0) {
- // all additional unavailable area moved to NATURAL land cover
- landCoverAreas.addArea(LandCoverType.NATURAL, LandProtectionType.UNAVAILABLE, additionalUnavail);
- landCoverAreas.removeArea(lcType, LandProtectionType.CONVERTIBLE, additionalUnavail * cProp);
- landCoverAreas.removeArea(lcType, LandProtectionType.PROTECTED, additionalUnavail * (1 - cProp));
- shortfall -= additionalUnavail;
- }
-
- if (shortfall < 1e-10)
- break;
- }
-
- if (shortfall >= 1e-10) {
- LogWriter.printlnWarning("LandUseItem.setUnavailableArea insufficient area to make unavailable");
- }
- }
-
public void runAreaCheck(Map<LandCoverType, Double> lcAreas, LandCoverTile landCoverAreas) {
final double THRESHOLD = 1e-7;
for (LandCoverType lcType : LandCoverType.values()) {
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index fbe118bf6b229a827b1f23b1fddaa76e195c4f5a..4dda40832835f95184fd787d2fc4bc5acaa5e73b 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -7,6 +7,7 @@ import java.util.Map;
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.types.*;
import ac.ed.lurg.utils.Interpolator;
+import ac.ed.lurg.utils.LogWriter;
import ac.sac.raster.InterpolatingRasterItem;
public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serializable {
@@ -18,6 +19,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
private Map<LccKey, Double> landCoverChange = new HashMap<LccKey, Double>(); // land cover change for previous timestep
private double protectedFraction; // protected fraction of total area
private double forestRotationIntensity;
+ private double maxCroplandArea = 0;
public LandUseItem() {}
@@ -27,9 +29,11 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
for (CropType crop : CropType.getCropsLessPasture()) {
setCropFraction(crop, landCover.getCropFraction(crop));
}
-
+ setMaxCroplandArea(landCover.getMaxCroplandArea());
setProtectedFraction(landCover.getProtectedFraction());
this.landCoverAreas = (landCover.getLandCoverTiles());
+ } else {
+ maxCroplandArea = getTotalLandArea();
}
forestRotationIntensity = 0.02;
}
@@ -41,6 +45,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
landCoverAreas = (luItemToCopy.landCoverAreas);
protectedFraction = (luItemToCopy.protectedFraction);
forestRotationIntensity = (luItemToCopy.forestRotationIntensity);
+ maxCroplandArea = (luItemToCopy.maxCroplandArea);
}
public void overwriteLandCoverAreas(LandCoverTile landCoverAreas) {
@@ -100,6 +105,21 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return totalArea==0 ? 0.0 : getLandCoverArea(c) / totalArea;
}
+ public void setMaxCroplandArea(double maxCroplandArea) {
+ this.maxCroplandArea = maxCroplandArea;
+ }
+
+ public double getMaxCroplandArea() {
+ return maxCroplandArea;
+ }
+
+ public double getAreaSuitableForCropland() { // how much land can be converted to cropland
+ if (getLandCoverArea(LandCoverType.CROPLAND) > getMaxCroplandArea()) {
+ LogWriter.printlnWarning("Cropland exceeds maximum area!");
+ }
+ return Math.max(0, getMaxCroplandArea() - getLandCoverArea(LandCoverType.CROPLAND));
+ }
+
public double getSuitableArea() {
double suitableArea = 0;
for (LandCoverType lcType : LandCoverType.getConvertibleTypes()) {