From bba5e9ca24db780537b1220594295792061941d8 Mon Sep 17 00:00:00 2001 From: Roslyn Henry <roslyn.henry@ed.ac.uk> Date: Wed, 13 Jan 2021 15:09:04 +0000 Subject: [PATCH] Fix for land allocation error in allocAreas method. Unprotected natural area was being calculated on more that one occasion within a timestep resulting in changing proportions of 'other natural' and 'forest' that should be protected. Area(s) of natural land that should be protected now updated at the start of each timestep and remain constant throughout the timestep. --- .../country/gams/GamsRasterOptimiser.java | 20 ++++---- src/ac/ed/lurg/landuse/LandUseItem.java | 46 +++++++++++-------- 2 files changed, 37 insertions(+), 29 deletions(-) diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java index 7c516c50..d0504d77 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java @@ -110,11 +110,11 @@ public class GamsRasterOptimiser { double gamsPastureChange = newLandUseAggItem.getLandCoverArea(LandCoverType.PASTURE) - prevLandUseAggItem.getLandCoverArea(LandCoverType.PASTURE); double gamsCroplandChange = newLandUseAggItem.getLandCoverArea(LandCoverType.CROPLAND) - prevLandUseAggItem.getLandCoverArea(LandCoverType.CROPLAND); RasterSet<LandUseItem> landUseItemsForLocation = newLandUseRaster.createSubsetForKeys(keys); - + if (DEBUG) { checkedTotalAreas(landUseItemsForLocation, locId + " before"); - LogWriter.println("pastureChange" + gamsPastureChange); - LogWriter.println("croplandChange" + gamsCroplandChange); + LogWriter.println("pastureChange " + gamsPastureChange); + LogWriter.println("croplandChange " + gamsCroplandChange); } double prevManagedForest = 0, prevUnmanagedForest = 0, prevNatural = 0, protectedAreasPastureChange = 0, protectedAreasCroplandChange = 0; @@ -122,13 +122,13 @@ public class GamsRasterOptimiser { for (LandUseItem luItem: landUseItemsForLocation.values()) { prevManagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.MANAGED_FOREST); prevUnmanagedForest += luItem.getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST); - prevNatural += luItem.getUnprotectedNatural(); + prevNatural += luItem.getTotalUnprotectedNatural(); } double prevForestTotal = prevManagedForest + prevUnmanagedForest; double prevForestToNaturalFraction = (prevNatural > 0) ? prevForestTotal / prevNatural : 0; double prevForestManagedFraction = (prevForestTotal > 0) ? prevManagedForest / prevForestTotal : 0; - + //this tries to honour spatial distribution of protected areas, force removal of pasture and cropland from protected areas if (ModelConfig.FORCE_PROTECTED_AREAS && year >= ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR) { double protectedLandShortfall = 0; @@ -206,9 +206,9 @@ public class GamsRasterOptimiser { } if (DEBUG) { - LogWriter.println("pastureFromCrop" + pastureFromCrop); - LogWriter.println("pastureFromNatural" + pastureFromNatural); - LogWriter.println("cropFromNatural" + cropFromNatural); + LogWriter.println("pastureFromCrop " + pastureFromCrop); + LogWriter.println("pastureFromNatural " + pastureFromNatural); + LogWriter.println("cropFromNatural " + cropFromNatural); } allocAllLandCropsTop(newLandUseRaster, keys, LandCoverType.PASTURE, LandCoverType.MANAGED_FOREST, pastureFromNatural * prevForestToNaturalFraction * prevForestManagedFraction, locId); @@ -232,7 +232,7 @@ public class GamsRasterOptimiser { } } } - + return newLandUseRaster; } @@ -262,7 +262,7 @@ public class GamsRasterOptimiser { double totalShortfall = 0; Set<RasterKey> keysWithSpace = new HashSet<RasterKey>(); double totalUnprotectedLC = 0; - + for (RasterKey key : keys) { LandUseItem newLandUseItem = newLandUseRaster.get(key); totalUnprotectedLC += newLandUseItem.getUnprotectedLandCoverArea(fromType); diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java index 782ff889..32f285bf 100644 --- a/src/ac/ed/lurg/landuse/LandUseItem.java +++ b/src/ac/ed/lurg/landuse/LandUseItem.java @@ -18,6 +18,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial private Map<CropType, Intensity> intensityMap = new HashMap<CropType, Intensity>(); private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>(); private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>(); + private Map<LandCoverType, Double> unprotectedAreas = new HashMap<LandCoverType, Double>(); private double protectedArea; //protected area in Mha private double unavailableArea; //area unavailable due to altitude etc private double suitableArea; @@ -42,6 +43,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial intensityMap.putAll(luItemToCopy.intensityMap); cropFractions.putAll(luItemToCopy.cropFractions); landCoverAreas.putAll(luItemToCopy.landCoverAreas); + unprotectedAreas.putAll(luItemToCopy.unprotectedAreas); protectedArea = (luItemToCopy.protectedArea); suitableArea = (luItemToCopy.suitableArea); } @@ -193,7 +195,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial double prevTo = getLandCoverArea(toType); double prevFrom = getLandCoverArea(fromType); - double availLC = fromType.isProtectable() ? Math.min(prevFrom, getUnprotectedNatural()) : prevFrom; // if protected can not go past protected area + double availLC = fromType.isProtectable() ? Math.min(prevFrom, unprotectedAreas.get(fromType)) : prevFrom; // if protected can not go past protected area double actuallyChanged = Math.min(availLC, changeReq); @@ -237,7 +239,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial } public void setLandCoverArea(LandCoverType c, double d) { - if (Double.isNaN(d) || Double.isInfinite(d)) + if (Double.isNaN(d) || Double.isInfinite(d) || d < 0) throw new RuntimeException("AreasItem for " + c + " is " + d); landCoverAreas.put(c, d); @@ -280,19 +282,18 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial return totalNatural; } - /** averages protected areas across land cover types */ - public double getUnprotectedLandCoverArea(LandCoverType c) { - double lcTypeArea = getLandCoverArea(c); - - if (c.isProtectable()) { - double unprotectedNat = getUnprotectedNatural(); - double totalNatural = getTotalNatural(); - double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0; - return lcTypeArea * unprotectedFract; - } - else { - return lcTypeArea; + public double getTotalUnprotectedNatural() { + double unprotectedNatural = 0; + for (LandCoverType landType : LandCoverType.values()) { + if (landType.isProtectable()) { + unprotectedNatural += unprotectedAreas.get(landType); + } } + return unprotectedNatural; + } + + public double getUnprotectedLandCoverArea(LandCoverType landType) { + return unprotectedAreas.get(landType); } private double getProtectedandUnavailable() { @@ -303,16 +304,23 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial return Math.min(totalSuitableLandCover, minAndProtectedA); } - - public double getUnprotectedNatural() { + /** averages protected areas across land cover types */ + private void updateUnprotectedAreas() { double desiredProtected = getProtectedandUnavailable(); double totalNatural = getTotalNatural(); double unprotectedNat = Math.max(0, totalNatural - desiredProtected); // if we are already using more than is protected then the unprotectedArea is 0 - //LogWriter.println(" totalNatural - desiredProtected " + (totalNatural - desiredProtected) + "unprotectedNat " + unprotectedNat); - return unprotectedNat; + double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0; + + for (LandCoverType landType : LandCoverType.values()) { + if (landType.isProtectable()) { + unprotectedAreas.put(landType, getLandCoverArea(landType)*unprotectedFract); + } + else unprotectedAreas.put(landType, getLandCoverArea(landType)); + } } public void updateSuitableArea(int year) { + updateUnprotectedAreas(); //update unprotected areas double currentAgri = getLandCoverArea(LandCoverType.PASTURE) + getLandCoverArea(LandCoverType.CROPLAND); double totalNatural = getTotalNatural(); double natAvailForAgri = totalNatural - getProtectedandUnavailable(); // could be negative, i.e. excess agricultural area already used @@ -327,7 +335,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial suitableArea = Math.max(0, currentAgri + natAvailForAgri * proportion); // netNatAvailForAgri is negative, but suitable area < 0 is not sensible (seems to happen with high barren areas) } else //not honouring protected areas if agriculture > natural that should be protected - suitableArea = currentAgri + getUnprotectedNatural(); + suitableArea = currentAgri + getTotalUnprotectedNatural(); } public double getForestManagedFraction() { -- GitLab