diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms index 87cdb46b6d4eb88284092adb4a19abff4cede7d8..2abce6ab8cb4590cabc29133200aa5886abbb6a5 100644 --- a/GAMS/IntExtOpt.gms +++ b/GAMS/IntExtOpt.gms @@ -270,7 +270,7 @@ $gdxin CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location)); - CARBON_CREDIT_CALC(location) .. carbonCredits(location) =E= -sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonCreditRate(land_cover_before, land_cover_after, location)); + CARBON_CREDIT_CALC(location) .. carbonCredits(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonCreditRate(land_cover_before, land_cover_after, location)); CARBON_CREDIT_CONSTRAINT .. sum(location, carbonCredits(location)) =G= demand('carbonCredits') + exportAmount('carbonCredits') - importAmount('carbonCredits'); diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java index 537f4feb18cb995aa84ca54dfc918a8c7d8ff7ce..62149474472e7300dc2cde1d998da4a2d01ff68f 100755 --- a/src/ac/ed/lurg/ModelConfig.java +++ b/src/ac/ed/lurg/ModelConfig.java @@ -516,6 +516,7 @@ public class ModelConfig { public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.2); //$1000/tC public static final double INIT_WOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_WOOD_PRICE", 0.38) : 0.0; // $1000/tC-eq public static final boolean ENABLE_VEGETATION_CLEARANCE_COST = getBooleanProperty("ENABLE_VEGETATION_CLEARANCE_COST", false); // calculates cost of clearing vegetation even if forestry off + public static final int MIN_FOREST_ROTATION_PERIOD = getIntProperty("MIN_FOREST_ROTATION_PERIOD", 5); // Minimum forest rotation period; Bauhus et al., 2010. Ecosystem goods and services from plantation forests. // Carbon public static final boolean IS_CARBON_ON = getBooleanProperty("IS_CARBON_ON", false); diff --git a/src/ac/ed/lurg/carbon/CarbonFluxItem.java b/src/ac/ed/lurg/carbon/CarbonFluxItem.java index e62718d6372b1bb10d0cac2dd0f150ee91e1444c..e50123495f6554d3630c0e1f0a354da59133a4e2 100644 --- a/src/ac/ed/lurg/carbon/CarbonFluxItem.java +++ b/src/ac/ed/lurg/carbon/CarbonFluxItem.java @@ -35,7 +35,7 @@ public class CarbonFluxItem implements RasterItem, Serializable { } double meanFlux = totalFlux / totalArea; carbonFluxRate.put(key, meanFlux); - carbonCreditRate.put(key, meanFlux); + carbonCreditRate.put(key, -meanFlux); } } @@ -58,24 +58,24 @@ public class CarbonFluxItem implements RasterItem, Serializable { double meanFlux = totalFlux / totalArea; carbonFluxRate.put(key, meanFlux); - // Carbon credit rate + // Carbon credit rate for unchanged land cover for (int age : tiles.getAgeKeys()) { for (int i = age; i < ModelConfig.CARBON_HORIZON + age; i++) { int ageCapped = Math.min(i, ModelConfig.CARBON_WOOD_MAX_TIME - 1); totalFlux += cFluxes[ageCapped] * tiles.getArea(LandProtectionType.CONVERTIBLE, age); } } - double annualisedFlux = (totalFlux / totalArea) / ModelConfig.CARBON_HORIZON; - carbonCreditRate.put(key, annualisedFlux); + double annualisedFlux = (totalFlux / totalArea) / ModelConfig.CARBON_HORIZON; + carbonCreditRate.put(key, -annualisedFlux); // switching sign as we want +ve to mean C offset } - // Credit rate + // Credit rate for destination land cover { double totalFlux = 0; for (int i = 0; i < ModelConfig.CARBON_HORIZON; i++) { totalFlux += cFluxes[i]; } - double annualisedFlux = totalFlux / ModelConfig.CARBON_HORIZON; + double annualisedFlux = -totalFlux / ModelConfig.CARBON_HORIZON; // switching sign as we want +ve to mean C offset for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) { if (!fromLc.equals(lcType)) { LccKey lccKey = new LccKey(fromLc, lcType); diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java index 30cd5e2fe1cc172b1edc3f9e7200bfd8926b889a..f1f9ffab8f493d318a32a72df050c211c5024615 100644 --- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java @@ -395,7 +395,7 @@ public class GamsLocationOptimiser { v.add(locString); setGamsParamValue(woodYieldRotaP.addRecord(v), wYield.getYieldRota(), 5); - for (LandCoverType fromLc : LandCoverType.getNaturalTypes()) { + for (LandCoverType fromLc : LandCoverType.getForestedTypes()) { for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) { if (fromLc.equals(toLc)) continue; @@ -614,8 +614,8 @@ public class GamsLocationOptimiser { // Timber harvest - // Ignore harvest from LUC during initial calib to avoid flooding the market. - double woodSupplyLuc = (inputData.getTimestep().getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) ? 0.0 : outDB.getVariable("woodSupplyLuc").getFirstRecord().getLevel(); + // Ignore harvest from LUC during calibration as we want all wood demand to be met by rotation harvest + double woodSupplyLuc = (ModelConfig.IS_CALIBRATION_RUN) ? 0.0 : outDB.getVariable("woodSupplyLuc").getFirstRecord().getLevel(); double woodSupplyRota = outDB.getVariable("woodSupplyRota").getFirstRecord().getLevel(); double totalWoodHarvest = woodSupplyLuc + woodSupplyRota; double netWoodImport = getParmValue(parmNetImports, "wood"); @@ -638,7 +638,7 @@ public class GamsLocationOptimiser { } // Carbon - double carbonSequestered = -outDB.getParameter("netCarbonCredits").getFirstRecord().getValue(); + double carbonSequestered = outDB.getParameter("netCarbonCredits").getFirstRecord().getValue(); double netCarbonImport = getParmValue(parmNetImports, "carbonCredits"); double netCarbonFlux = outDB.getParameter("netCarbonFlux").getFirstRecord().getValue(); CarbonUsageData carbonUsageData = new CarbonUsageData(carbonSequestered, netCarbonImport, netCarbonFlux); diff --git a/src/ac/ed/lurg/demand/SspManager.java b/src/ac/ed/lurg/demand/SspManager.java index e284d666d17eaf1487aff5370a646021be4f6d97..6f22b54b567c1d5939736c40eddff7af750deb59 100644 --- a/src/ac/ed/lurg/demand/SspManager.java +++ b/src/ac/ed/lurg/demand/SspManager.java @@ -95,7 +95,7 @@ public class SspManager { return new SspData(pop, gdp); } else { - LogWriter.printlnWarning("Can't find straddling year for " + year + ", " + scenario + ", " + country.getCountryName()); + LogWriter.println("Warning: Can't find straddling year for " + year + ", " + scenario + ", " + country.getCountryName(), 3); return null; } } diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java index 7b8becb406d4050d202aede51f1e4197dd9b48b5..222d4994b6190ab5878311cce79b8cb39815c8a9 100644 --- a/src/ac/ed/lurg/forestry/WoodYieldItem.java +++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java @@ -30,7 +30,7 @@ public class WoodYieldItem implements RasterItem { LandCoverTile tiles = landUseTiles.get(key.getFromLc()); double totalArea = tiles.getTotalArea(LandProtectionType.CONVERTIBLE); // Assuming no harvest from protected areas - if (totalArea < 1e-6) { + if (totalArea < 1e-6) { // Avoid floating point errors due to small areas this.yield.put(key, 0.0) ; } else { double totalYield = 0; @@ -77,6 +77,7 @@ public class WoodYieldItem implements RasterItem { } optimalRotation = Collections.min(candidates); // Choose shortest rotation + optimalRotation = Math.max(ModelConfig.MIN_FOREST_ROTATION_PERIOD, optimalRotation); // set lower bound yieldAtRotation = yields[optimalRotation] / optimalRotation; } diff --git a/src/ac/ed/lurg/forestry/WoodYieldKey.java b/src/ac/ed/lurg/forestry/WoodYieldKey.java deleted file mode 100644 index 5929d90672188e51d9632337ffa68fb8230559e3..0000000000000000000000000000000000000000 --- a/src/ac/ed/lurg/forestry/WoodYieldKey.java +++ /dev/null @@ -1,56 +0,0 @@ -package ac.ed.lurg.forestry; - -import ac.ed.lurg.types.LandCoverType; - -public class WoodYieldKey { - private LandCoverType forestType; - private int year; - private int age; - - public WoodYieldKey(LandCoverType forestType, int year, int age) { - super(); - this.forestType = forestType; - this.year = year; - this.age = age; - } - - public LandCoverType getForestType() { - return forestType; - } - - public int getYear() { - return year; - } - - public int getAge() { - return age; - } - - @Override - public int hashCode() { - final int prime = 31; - int result = 1; - result = prime * result + age; - result = prime * result + ((forestType == null) ? 0 : forestType.hashCode()); - result = prime * result + year; - return result; - } - - @Override - public boolean equals(Object obj) { - if (this == obj) - return true; - if (obj == null) - return false; - if (getClass() != obj.getClass()) - return false; - WoodYieldKey other = (WoodYieldKey) obj; - if (age != other.age) - return false; - if (forestType != other.forestType) - return false; - if (year != other.year) - return false; - return true; - } -} diff --git a/src/ac/ed/lurg/forestry/WoodYieldReader.java b/src/ac/ed/lurg/forestry/WoodYieldReader.java index 9dc454f54fd9869ebc2f124abb8eae8b95817870..f7cfc47ee1c20647066095b261ec9191b5533181 100644 --- a/src/ac/ed/lurg/forestry/WoodYieldReader.java +++ b/src/ac/ed/lurg/forestry/WoodYieldReader.java @@ -1,119 +1,119 @@ -package ac.ed.lurg.forestry; - -import java.io.File; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Collection; -import java.util.HashMap; -import java.util.List; -import java.util.Map; - -import ac.ed.lurg.ModelConfig; -import ac.ed.lurg.Timestep; -import ac.ed.lurg.landuse.LandUseItem; -import ac.ed.lurg.landuse.LccKey; -import ac.ed.lurg.types.LandCoverType; -import ac.sac.raster.AbstractBinaryRasterReader; -import ac.sac.raster.RasterHeaderDetails; -import ac.sac.raster.RasterKey; -import ac.sac.raster.RasterSet; - -public class WoodYieldReader { - private static final int MIN_COLS = ModelConfig.CARBON_WOOD_MAX_TIME + 2 + 1; // number of years + 2 for coords + 1 for age=0 - private static final double CONVERSION_FACTOR = 10.0 * ModelConfig.WOOD_YIELD_CALIB_FACTOR; // convert kgC/m2 to tC/ha with calib factor - private RasterHeaderDetails rasterProj; - private String[] header = new String[MIN_COLS]; - - public WoodYieldReader(RasterHeaderDetails rasterProj) { - this.rasterProj = rasterProj; - header[0] = "Lon"; - header[1] = "Lat"; - for (int i = 0; i < MIN_COLS - 1; i++) { - header[i] = Integer.toString(i); - } - } - - public WoodYieldRasterSet getWoodYields(RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) { - WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(rasterProj); - - // Wood yield from rotation - List<LccKey> rotaMapping = new ArrayList<LccKey>(); - rotaMapping.add(new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST)); // Not using getLccMappings as that excludes X to X conversions - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_ROTA_FILENAME, rotaMapping , true); - - // Wood yield from LCC - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_PAST_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.PASTURE)), false); - - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_CROP_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.CROPLAND)), false); - - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_NTRL_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.NATURAL)), false); - - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_FORS_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), LandCoverType.getManagedForestTypes()), false); - - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_NTRL_TO_PAST_FILENAME, - getLccMappings(Arrays.asList(LandCoverType.NATURAL), Arrays.asList(LandCoverType.PASTURE)), false); - - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_NTRL_TO_CROP_FILENAME, - getLccMappings(Arrays.asList(LandCoverType.NATURAL), Arrays.asList(LandCoverType.CROPLAND)), false); - - getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_NTRL_TO_FORS_FILENAME, - getLccMappings(Arrays.asList(LandCoverType.NATURAL), LandCoverType.getManagedForestTypes()), false); - - return woodYieldData; - } - - public void getWoodYields(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice, - String filePath, List<LccKey> lccMappings, boolean calcRotation) { - - AbstractBinaryRasterReader<WoodYieldItem> yieldReader = new AbstractBinaryRasterReader<WoodYieldItem>(header, MIN_COLS - 2, woodYieldData) { - protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) { - - if (!landUseRaster.containsKey(key)) { - return; - } - - Double[] yields = getYieldsFromRowValues(rowValues); - Map<LccKey, Double[]> yieldMap = new HashMap<LccKey, Double[]>(); - for (LccKey lccKey : lccMappings) { - yieldMap.put(lccKey, yields); - } - - item.calcYieldData(yieldMap, landUseRaster.get(key).getLandCoverTiles(), timestep); - - if (calcRotation) // only need to do this for forest-to-forest data - item.calcRotationData(yieldMap, landUseRaster.get(key).getLandCoverTiles(), timestep, woodPrice); - } - }; - - String filename = getDataDir(filePath, timestep); - yieldReader.getRasterDataFromFile(filename); - } - - private ArrayList<LccKey> getLccMappings(Collection<LandCoverType> fromTypes, Collection<LandCoverType> toTypes) { - ArrayList<LccKey> lccMappings = new ArrayList<LccKey>(); - for (LandCoverType fromLc : fromTypes) { - for (LandCoverType toLc : toTypes) { - if (!fromLc.equals(toLc)) // exclude no change as yield should be 0 - lccMappings.add(new LccKey(fromLc, toLc)); - } - } - return lccMappings; - } - - private String getDataDir(String filename, Timestep timestep) { - return timestep.getWoodAndCarbonYearSubDir(ModelConfig.FORESTRY_DIR) + File.separator + filename; - } - - private Double[] getYieldsFromRowValues(Map<String, Double> rowValues) { - Double[] yields = new Double[ModelConfig.CARBON_WOOD_MAX_TIME]; - for (int i=0; i<ModelConfig.CARBON_WOOD_MAX_TIME; i++) { - yields[i] = rowValues.get(Integer.toString(i)) * CONVERSION_FACTOR; - } - return yields; - } - -} +package ac.ed.lurg.forestry; + +import java.io.File; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collection; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + +import ac.ed.lurg.ModelConfig; +import ac.ed.lurg.Timestep; +import ac.ed.lurg.landuse.LandUseItem; +import ac.ed.lurg.landuse.LccKey; +import ac.ed.lurg.types.LandCoverType; +import ac.sac.raster.AbstractBinaryRasterReader; +import ac.sac.raster.RasterHeaderDetails; +import ac.sac.raster.RasterKey; +import ac.sac.raster.RasterSet; + +public class WoodYieldReader { + private static final int MIN_COLS = ModelConfig.CARBON_WOOD_MAX_TIME + 2 + 1; // number of years + 2 for coords + 1 for age=0 + private static final double CONVERSION_FACTOR = 10.0 * ModelConfig.WOOD_YIELD_CALIB_FACTOR; // convert kgC/m2 to tC/ha with calib factor + private RasterHeaderDetails rasterProj; + private String[] header = new String[MIN_COLS]; + + public WoodYieldReader(RasterHeaderDetails rasterProj) { + this.rasterProj = rasterProj; + header[0] = "Lon"; + header[1] = "Lat"; + for (int i = 0; i < MIN_COLS - 1; i++) { + header[i] = Integer.toString(i); + } + } + + public WoodYieldRasterSet getWoodYields(RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) { + WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(rasterProj); + + // Wood yield from rotation + List<LccKey> rotaMapping = new ArrayList<LccKey>(); + rotaMapping.add(new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST)); // Not using getLccMappings as that excludes X to X conversions + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_ROTA_FILENAME, rotaMapping , true); + + // Wood yield from LCC + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_PAST_FILENAME, + getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.PASTURE)), false); + + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_CROP_FILENAME, + getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.CROPLAND)), false); + + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_NTRL_FILENAME, + getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.NATURAL)), false); + + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_FORS_TO_FORS_FILENAME, + getLccMappings(LandCoverType.getManagedForestTypes(), LandCoverType.getManagedForestTypes()), false); + + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_NTRL_TO_PAST_FILENAME, + getLccMappings(Arrays.asList(LandCoverType.NATURAL), Arrays.asList(LandCoverType.PASTURE)), false); + + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_NTRL_TO_CROP_FILENAME, + getLccMappings(Arrays.asList(LandCoverType.NATURAL), Arrays.asList(LandCoverType.CROPLAND)), false); + + getWoodYields(woodYieldData, landUseRaster, timestep, woodPrice, ModelConfig.WOOD_YIELD_NTRL_TO_FORS_FILENAME, + getLccMappings(Arrays.asList(LandCoverType.NATURAL), LandCoverType.getManagedForestTypes()), false); + + return woodYieldData; + } + + public void getWoodYields(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice, + String filePath, List<LccKey> lccMappings, boolean calcRotation) { + + AbstractBinaryRasterReader<WoodYieldItem> yieldReader = new AbstractBinaryRasterReader<WoodYieldItem>(header, MIN_COLS - 2, woodYieldData) { + protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) { + + if (!landUseRaster.containsKey(key)) { + return; + } + + Double[] yields = getYieldsFromRowValues(rowValues); + Map<LccKey, Double[]> yieldMap = new HashMap<LccKey, Double[]>(); + for (LccKey lccKey : lccMappings) { + yieldMap.put(lccKey, yields); + } + + item.calcYieldData(yieldMap, landUseRaster.get(key).getLandCoverTiles(), timestep); + + if (calcRotation) // only need to do this for forest-to-forest data + item.calcRotationData(yieldMap, landUseRaster.get(key).getLandCoverTiles(), timestep, woodPrice); + } + }; + + String filename = getDataDir(filePath, timestep); + yieldReader.getRasterDataFromFile(filename); + } + + private ArrayList<LccKey> getLccMappings(Collection<LandCoverType> fromTypes, Collection<LandCoverType> toTypes) { + ArrayList<LccKey> lccMappings = new ArrayList<LccKey>(); + for (LandCoverType fromLc : fromTypes) { + for (LandCoverType toLc : toTypes) { + if (!fromLc.equals(toLc)) // exclude no change as yield should be 0 + lccMappings.add(new LccKey(fromLc, toLc)); + } + } + return lccMappings; + } + + private String getDataDir(String filename, Timestep timestep) { + return timestep.getWoodAndCarbonYearSubDir(ModelConfig.FORESTRY_DIR) + File.separator + filename; + } + + private Double[] getYieldsFromRowValues(Map<String, Double> rowValues) { + Double[] yields = new Double[ModelConfig.CARBON_WOOD_MAX_TIME]; + for (int i=0; i<ModelConfig.CARBON_WOOD_MAX_TIME; i++) { + yields[i] = rowValues.get(Integer.toString(i)) * CONVERSION_FACTOR; + } + return yields; + } + +} diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java index d9aac67874a670c815ea53063c0039a2ef01a919..dad537656bd284596278af324c4787dd6c8bc0df 100644 --- a/src/ac/ed/lurg/types/LandCoverType.java +++ b/src/ac/ed/lurg/types/LandCoverType.java @@ -105,7 +105,7 @@ public enum LandCoverType { } - public static Collection<LandCoverType> getNaturalTypes() { + public static Collection<LandCoverType> getForestedTypes() { Collection<LandCoverType> naturalTypes = new HashSet<LandCoverType>();