From ca481664bab1f91e96baf52951f48c98a1dc15a4 Mon Sep 17 00:00:00 2001 From: Peter Alexander <peter@blackhillock.co.uk> Date: Fri, 12 Dec 2014 10:46:26 +0000 Subject: [PATCH] no message --- .../country/gams/GamsRasterOptimiser.java | 79 ++++++++++++++++--- src/ac/ed/lurg/country/gams/GamsTest.java | 2 +- src/ac/ed/lurg/landuse/AreasItem.java | 20 +++++ src/ac/ed/lurg/yield/YieldResponse.java | 20 +++-- 4 files changed, 103 insertions(+), 18 deletions(-) diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java index 34730bfa..67d1ab7d 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java @@ -12,6 +12,7 @@ import java.util.Set; import ac.ed.lurg.landuse.AreasItem; import ac.ed.lurg.landuse.IntensitiesItem; import ac.ed.lurg.types.CropType; +import ac.ed.lurg.types.LandDataType; import ac.ed.lurg.types.YieldType; import ac.ed.lurg.utils.LogWriter; import ac.ed.lurg.yield.AveragingYieldResponsesItem; @@ -42,20 +43,78 @@ public class GamsRasterOptimiser { } private GamsRasterOutput convertToRaster(GamsInput gamsInput, GamsOutput gamsOutput) { - // intensities - constant over single category? - RasterSet<IntensitiesItem> intensityRaster = new RasterSet<IntensitiesItem>(); + + RasterSet<AreasItem> newAreaRaster = new RasterSet<AreasItem>(); + RasterSet<IntensitiesItem> newIntensityRaster = new RasterSet<IntensitiesItem>(); - for (Map.Entry<Integer, IntensitiesItem> entry : gamsOutput.getIntensities().entrySet()) { - for (RasterKey key : mapping.get(entry.getKey())) { - intensityRaster.put(key, entry.getValue()); + for (Map.Entry<Integer, AreasItem> entry : gamsOutput.getCropAreas().entrySet()) { + Integer locId = entry.getKey(); + AreasItem newAreaAggItem = entry.getValue(); + + Set<RasterKey> keys = mapping.get(locId); + AreasItem prevAreaAggItem = gamsInput.getPreviousAreas().get(locId); + IntensitiesItem newIntensitiesItem = gamsOutput.getIntensities().get(locId); + Map<CropType, Double> avgCropChangeTotals = newAreaAggItem.getScaledCropChanges(prevAreaAggItem, 1.0d/keys.size()); + + double totalShortfall = 0; + Set<RasterKey> keysWithSpace = new HashSet<RasterKey>(); + + for (RasterKey key : keys) { + AreasItem prevAreasRasterItem = rasterInputData.getPreviousAreas().get(key); + + AreasItem newAreasRasterItem = new AreasItem(); + double shortfall = allocationChangeArea(prevAreasRasterItem, avgCropChangeTotals, newAreasRasterItem); + if (shortfall == 0) + keysWithSpace.add(key); + else + totalShortfall += shortfall; + + newAreaRaster.put(key, newAreasRasterItem); + newIntensityRaster.put(key, newIntensitiesItem); // intensities constant over single aggregated land category } } - RasterSet<AreasItem> areaRaster = new RasterSet<AreasItem>(); + return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getFeedAmounts(), gamsOutput.getNetImports()); + } - // crop areas - need some allocation process - - return null; + private double allocationChangeArea(AreasItem prevAreasRasterItem, + Map<CropType, Double> cropChanges, + AreasItem newAreasRasterItem) { + + double netAddedCrop = 0; + for (double d : cropChanges.values()) + netAddedCrop += d; + + double forestRemoved = netAddedCrop/2; + double naturalRemoved = netAddedCrop/2; + double prevForest = prevAreasRasterItem.getLandCoverArea(LandDataType.FOREST); + double prevNatural = prevAreasRasterItem.getLandCoverArea(LandDataType.OTHER_NATURAL); + + if (prevForest < forestRemoved) { + forestRemoved = prevForest; + naturalRemoved = netAddedCrop - forestRemoved; + } + + if (prevNatural < naturalRemoved) { + naturalRemoved = prevNatural; + } + + double ratioOfChanges = (forestRemoved + naturalRemoved)/netAddedCrop; + double shortfall = netAddedCrop - forestRemoved - naturalRemoved; + if (shortfall > 0) + LogWriter.printlnError("Not able to incorporate all changes, as not enough forest or natural areas left: " + ratioOfChanges); + + newAreasRasterItem.setLandCoverArea(LandDataType.LAND, prevAreasRasterItem.getLandCoverArea(LandDataType.LAND)); + newAreasRasterItem.setLandCoverArea(LandDataType.FOREST, prevForest - forestRemoved); + newAreasRasterItem.setLandCoverArea(LandDataType.OTHER_NATURAL, prevNatural - naturalRemoved); + + for (CropType crop : CropType.values()) { + double prevCropArea = prevAreasRasterItem.getCropArea(crop); + double additionalArea = cropChanges.get(crop) * ratioOfChanges; + newAreasRasterItem.setCropArea(crop, prevCropArea + additionalArea); + } + + return shortfall; } @@ -92,6 +151,8 @@ public class GamsRasterOptimiser { } aggAreas.setCropArea(crop, aggAreas.getCropArea(crop) + cropAreas.getCropArea(crop)); } + + aggAreas.setLandCoverArea(LandDataType.LAND, aggAreas.getLandCoverArea(LandDataType.LAND) + cropAreas.getLandCoverArea(LandDataType.LAND)); } return new GamsInput(aggregatedYield, aggregatedAreas, rasterInputData.getCountryInput()); diff --git a/src/ac/ed/lurg/country/gams/GamsTest.java b/src/ac/ed/lurg/country/gams/GamsTest.java index befa8f06..1b491ff0 100644 --- a/src/ac/ed/lurg/country/gams/GamsTest.java +++ b/src/ac/ed/lurg/country/gams/GamsTest.java @@ -22,7 +22,7 @@ public class GamsTest { GamsInput gamsInput = new GamsInput(getYields(), getPreviousArea(), countryLevelInputs); GamsLandUseOptimiser opti = new GamsLandUseOptimiser(gamsInput); - GamsOutput gamsOutput = opti.run(); + opti.run(); } public Map<CropType, Double> getProjectedDemand() { diff --git a/src/ac/ed/lurg/landuse/AreasItem.java b/src/ac/ed/lurg/landuse/AreasItem.java index e4cdd724..f4d6bea6 100644 --- a/src/ac/ed/lurg/landuse/AreasItem.java +++ b/src/ac/ed/lurg/landuse/AreasItem.java @@ -2,6 +2,7 @@ package ac.ed.lurg.landuse; import java.util.HashMap; import java.util.Map; +import java.util.Map.Entry; import ac.ed.lurg.types.CropType; import ac.ed.lurg.types.LandDataType; @@ -33,4 +34,23 @@ public class AreasItem implements RasterItem { public void setLandCoverArea(LandDataType c, double d) { landCoverAreas.put(c, d); } + + public double getTotalCropIncPastureArea() { + double total = 0; + for (double d : cropAreas.values()) { + total += d; + } + return total; + } + + public Map<CropType, Double> getScaledCropChanges(AreasItem prevAreaAggItem, double factor) { + Map<CropType, Double> changes = new HashMap<CropType, Double>(); + + for (Entry<CropType, Double> entry : cropAreas.entrySet()) { + double change = entry.getValue() - prevAreaAggItem.getCropArea(entry.getKey()); + changes.put(entry.getKey(), change*factor); + } + + return changes; + } } diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java index bec1c20f..e11d380d 100644 --- a/src/ac/ed/lurg/yield/YieldResponse.java +++ b/src/ac/ed/lurg/yield/YieldResponse.java @@ -17,15 +17,19 @@ public class YieldResponse { } public double getFertParam() { - double fertMin = 0, fertMid = 0.5, fertMax = 1; - - return Math.log((yields.get(YieldType.FERT_MID_ONLY)-getYield(YieldType.NONE))/(yields.get(YieldType.FERT_MAX_ONLY)-getYield(YieldType.NONE))) / - Math.log((fertMid - fertMin)/(fertMax - fertMin))/fertMax; + return calcParam(getYield(YieldType.NONE), getYield(YieldType.FERT_MID_ONLY), getYield(YieldType.FERT_MAX_ONLY)); } public double getIrrigParam() { - double irrigMin = 0, irrigMid = 0.5, irrigMax = 1; - return Math.log((yields.get(YieldType.IRRIG_MID_ONLY)-getYield(YieldType.NONE))/(yields.get(YieldType.IRRIG_MAX_ONLY)-getYield(YieldType.NONE))) / - Math.log((irrigMid - irrigMin)/(irrigMax - irrigMin))/irrigMax; + return calcParam(getYield(YieldType.NONE), getYield(YieldType.IRRIG_MID_ONLY), getYield(YieldType.IRRIG_MAX_ONLY)); } -} + + private double calcParam(double yMin, double yMid, double yMax) { + double xMin = 0, xMid = 0.5, xMax = 1; + + if (yMid <= yMin || yMax <= yMid) + return 1; // default to linear + + return Math.log((yMid - yMin)/(yMax - yMin)) / Math.log((xMid - xMin)/(xMax - xMin)) / xMax; + } +} \ No newline at end of file -- GitLab