From 40e2d1fdbceba257e201fb37bbdf25114678edc3 Mon Sep 17 00:00:00 2001 From: Bart Arendarczyk <s1924442@ed.ac.uk> Date: Tue, 21 Mar 2023 15:50:55 +0000 Subject: [PATCH] Revised implementation of forestry and carbon. --- GAMS/IntExtOpt.gms | 83 +++++--- src/ac/ed/lurg/InternationalMarket.java | 2 +- src/ac/ed/lurg/ModelConfig.java | 17 +- src/ac/ed/lurg/ModelMain.java | 19 +- src/ac/ed/lurg/carbon/CarbonDataOutputer.java | 119 ------------ src/ac/ed/lurg/carbon/CarbonFluxItem.java | 91 +-------- src/ac/ed/lurg/carbon/CarbonFluxReader.java | 6 +- src/ac/ed/lurg/country/CountryAgent.java | 14 +- .../country/gams/GamsLocationOptimiser.java | 66 ++++--- .../lurg/country/gams/GamsLocationOutput.java | 12 +- .../country/gams/GamsRasterOptimiser.java | 80 ++------ .../ed/lurg/demand/AbstractDemandManager.java | 4 +- .../ed/lurg/demand/ElasticDemandManager.java | 2 +- src/ac/ed/lurg/demand/WoodDemandManager.java | 4 +- src/ac/ed/lurg/forestry/ForestManager.java | 178 ----------------- .../lurg/forestry/ForestryDataOutputer.java | 26 +-- src/ac/ed/lurg/forestry/WoodHarvestItem.java | 34 ---- src/ac/ed/lurg/forestry/WoodYieldItem.java | 110 +---------- src/ac/ed/lurg/forestry/WoodYieldReader.java | 42 +--- src/ac/ed/lurg/landuse/LandCoverItem.java | 70 +------ src/ac/ed/lurg/landuse/LandCoverTile.java | 109 +++-------- src/ac/ed/lurg/landuse/LandKey.java | 12 +- src/ac/ed/lurg/landuse/LandTileReader.java | 26 --- .../lurg/landuse/LandUseBinarySerializer.java | 180 ------------------ src/ac/ed/lurg/landuse/LandUseItem.java | 99 +++------- src/ac/ed/lurg/landuse/WoodUsageData.java | 31 ++- src/ac/ed/lurg/landuse/WoodUsageReader.java | 16 +- src/ac/ed/lurg/types/LandCoverType.java | 4 +- 28 files changed, 250 insertions(+), 1206 deletions(-) delete mode 100644 src/ac/ed/lurg/carbon/CarbonDataOutputer.java delete mode 100644 src/ac/ed/lurg/forestry/ForestManager.java delete mode 100644 src/ac/ed/lurg/forestry/WoodHarvestItem.java delete mode 100644 src/ac/ed/lurg/landuse/LandTileReader.java delete mode 100644 src/ac/ed/lurg/landuse/LandUseBinarySerializer.java diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms index a076089a..f487da0a 100644 --- a/GAMS/IntExtOpt.gms +++ b/GAMS/IntExtOpt.gms @@ -58,6 +58,7 @@ PARAMETER woodYieldParam(location); PARAMETER maxRotationIntensity(wood_type); PARAMETER forestBaseCost(managed_forest); + PARAMETER previousRotationIntensity(location); SCALAR meatEfficency efficiency of converting feed and pasture into animal products; SCALAR fertiliserUnitCost fert cost at max fert rate; @@ -69,6 +70,7 @@ SCALAR maxGrossLccRate max rate of gross land cover change; SCALAR maxFertChange max increase in total fertilizer use; SCALAR maxIrrigChange max increase in total irrigation use; + SCALAR tradeAdjustmentCostRate SCALAR forestManagementCost cost $1000 per ha; SCALAR vegClearingCostRate cost of clearing vegetation $1000 per tC; @@ -82,8 +84,8 @@ $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, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate -$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam, maxRotationIntensity, forestBaseCost -$load forestManagementCost, vegClearingCostRate, carbonForestMaxProportion, maxFertChange, maxIrrigChange +$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam, maxRotationIntensity, forestBaseCost, tradeAdjustmentCostRate +$load forestManagementCost, vegClearingCostRate, carbonForestMaxProportion, maxFertChange, maxIrrigChange, previousRotationIntensity $gdxin SCALAR delta "use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf" / 0.00000000001 /; @@ -109,8 +111,8 @@ $gdxin oilcrops 0.2 pulses 0.31 starchyRoots 3.14 - fruitveg 4.5 - sugar 3.8 + fruitveg 4.0 + sugar 3.0 energyCrops 0.34 / ; PARAMETER baseCost(crop); @@ -123,6 +125,7 @@ $gdxin PARAMETER animalProdCost(animal); animalProdCost('ruminants') = 0.28; animalProdCost('monogastrics') = 0.24; + VARIABLES unitCost(crop, location) cost per area for each crop - cost @@ -137,22 +140,25 @@ $gdxin importAmount(all_types) imports of crops and meat - Mt exportAmount(all_types) exports of crops and meat - Mt animalProd(animal) animal production + tradeAdjustmentCost(import_types) totalFeedDM total feed dry matter landCoverArea(land_cover, location) land cover area in Mha landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha totalConversionCost(location) land cover conversion cost - $1000 per ha - timberForestArea(wood_type, location) - woodYieldRota(wood_type, location) - rotationIntensity(wood_type, location) - forestryCost(location) + woodYieldRota(location) wood yield tC per ha + rotationIntensity(location) equivalent to 1 divided by rotation period + woodSupply(location) wood supply MtC + forestryCost(location) total cost in 1000$ + carbonCredits(location) + totalLandValue * A "artificial variable for debugging https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/" total_cost total cost of domestic supply including net imports; POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount, totalFeedDM, - landCoverArea, landCoverChange, totalConversionCost, animalProd, woodYieldRota, rotationIntensity, timberForestArea; + landCoverArea, landCoverChange, totalConversionCost, animalProd, woodYieldRota, rotationIntensity, woodSupply; * POSITIVE VARIABLE A; @@ -182,6 +188,7 @@ $gdxin MIN_NET_IMPORT_CONSTRAINT(import_types) constraint on min net imports PASTURE_IMPORT_CONSTRAINT constraint to not import pasture PASTURE_EXPORT_CONSTRAINT + TRADE_ADJUSTMENT_COST_CALC(import_types) SETASIDE_AREA_CALC(location) CROPLAND_LAND_COVER_CALC(location) cropland area equals sum of all crop areas @@ -190,12 +197,13 @@ $gdxin LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area CONVERSION_COST_CALC(location) cost of land cover conversion - WOOD_YIELD_CALC(wood_type, location) - WOOD_DEMAND_CONSTRAINT(wood_type) - ROTA_INTENSITY_MAX_CONSTRAINT(wood_type, location) - ROTA_INTENSITY_MIN_CONSTRAINT(wood_type, location) - TIMBER_FOREST_AREA_CONSTRAINT(location) - FORESTRY_COST_CALC(location) + WOOD_YIELD_CALC(location) wood yield at rotation - tC per ha + ROTA_INTENSITY_MAX_CONSTRAINT(location) constraint on maximum rotation intensity (= minimum rotation period) + ROTA_INTENSITY_MIN_CONSTRAINT(location) constraint on minimum rotation intensity (= maximum rotation period) + WOOD_SUPPLY_CALC(location) wood supply + WOOD_DEMAND_CONSTRAINT satisfy wood demand and trade + ROUNDWOOD_DEMAND_CONSTRAINT satisfy roundwood demand (can only by obtained from older trees) + FORESTRY_COST_CALC(location) total forestry cost CARBON_CREDIT_CALC(location) CARBON_CREDIT_CONSTRAINT @@ -279,19 +287,24 @@ $gdxin ************* Forestry *********************************** - WOOD_YIELD_CALC(wood_type, location) .. woodYieldRota(wood_type, location) =E= woodYieldMax(location) * (1 - exp(-woodYieldParam(location) / rotationIntensity(wood_type, location))) * rotationIntensity(wood_type, location); - - ROTA_INTENSITY_MAX_CONSTRAINT(wood_type, location) .. rotationIntensity(wood_type, location) =L= maxRotationIntensity(wood_type); + WOOD_YIELD_CALC(location) .. woodYieldRota(location) =E= woodYieldMax(location) * (1 - exp(-woodYieldParam(location) / rotationIntensity(location))); - ROTA_INTENSITY_MIN_CONSTRAINT(wood_type, location) .. rotationIntensity(wood_type, location) =G= 1 / 160; + ROTA_INTENSITY_MAX_CONSTRAINT(location) .. rotationIntensity(location) =L= 1 / 10; - TIMBER_FOREST_AREA_CONSTRAINT(location) .. sum(wood_type, timberForestArea(wood_type, location)) =E= landCoverArea('timberForest', location); + ROTA_INTENSITY_MIN_CONSTRAINT(location) .. rotationIntensity(location) =G= 1 / 150; - WOOD_DEMAND_CONSTRAINT(wood_type) .. demand(wood_type) + exportAmount(wood_type) - importAmount(wood_type) =L= sum(location, woodYieldRota(wood_type, location) * timberForestArea(wood_type, location)); - - FORESTRY_COST_CALC(location) .. forestryCost(location) =E= sum(managed_forest, landCoverArea(managed_forest, location) * forestBaseCost(managed_forest)) + - sum(wood_type, forestManagementCost * rotationIntensity(wood_type, location)); + WOOD_SUPPLY_CALC(location) .. woodSupply(location) =E= woodYieldRota(location) * rotationIntensity(location) * landCoverArea('timberForest', location); + + WOOD_DEMAND_CONSTRAINT .. sum(location, woodSupply(location)) =G= sum(wood_type, demand(wood_type) + exportAmount(wood_type) - importAmount(wood_type)); +* The sigmoid function here calculates the proportion of wood yield that is suitable for roundwood (maximum 85% with halfway point at 25 years). + ROUNDWOOD_DEMAND_CONSTRAINT .. sum(location, woodSupply(location) * 0.85 * sigmoid(0.2 * (1 / rotationIntensity(location) - 25))) =G= demand('roundwood') + exportAmount('roundwood') - importAmount('roundwood'); + + FORESTRY_COST_CALC(location) .. forestryCost(location) =E= sum(managed_forest, landCoverArea(managed_forest, location) * forestBaseCost(managed_forest)) + + forestManagementCost * rotationIntensity(location) * landCoverArea('timberForest', location) + + sum(land_cover$[not sameAs(land_cover, 'carbonForest')], landCoverChange(land_cover, 'carbonForest', location) * forestManagementCost); +* + sum(wood_type, 0.005 * (exp((exportAmount(wood_type) - previousExportAmount(wood_type))) - 1) * exportAmount(wood_type)); + *********** Carbon *********************************** 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)); @@ -304,8 +317,11 @@ $gdxin MAX_NET_IMPORT_CONSTRAINT(import_types) .. importAmount(import_types) - exportAmount(import_types) =L= maxNetImport(import_types); MIN_NET_IMPORT_CONSTRAINT(import_types) .. importAmount(import_types) - exportAmount(import_types) =G= minNetImport(import_types); + PASTURE_IMPORT_CONSTRAINT .. importAmount('pasture') =E= 0; PASTURE_EXPORT_CONSTRAINT .. exportAmount('pasture') =E= 0; + + TRADE_ADJUSTMENT_COST_CALC(import_types) .. tradeAdjustmentCost(import_types) =E= tradeAdjustmentCostRate * power((importAmount(import_types) - exportAmount(import_types)) - (previousImportAmount(import_types) - previousExportAmount(import_types)), 2); ************ Total cost ****************************** @@ -314,9 +330,13 @@ $gdxin sum(location, totalConversionCost(location)) + sum(location, forestryCost(location)) + sum(import_types, importPrices(import_types) * importAmount(import_types)) - - sum(import_types, exportPrices(import_types) * exportAmount(import_types)); + sum(import_types, exportPrices(import_types) * exportAmount(import_types)) + + sum(import_types, tradeAdjustmentCost(import_types)); + + - MODEL LAND_USE /ALL/ ; + MODEL LAND_USE / ALL / ; + fertI.L(crop, location) = previousFertIntensity(crop, location); irrigI.L(crop, location) = previousIrrigIntensity(crop, location); otherIntensity.L(crop, location) = previousOtherIntensity(crop, location); @@ -327,7 +347,8 @@ $gdxin monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop); importAmount.L(all_types) = previousImportAmount(all_types); exportAmount.L(all_types) = previousExportAmount(all_types); - rotationIntensity.L(wood_type, location) = 0.02; + rotationIntensity.L(location) = previousRotationIntensity(location); + * LAND_USE.OptFile = 1; @@ -352,7 +373,8 @@ $gdxin parameter productionShock(all_types); parameter carbonFlux(location); scalar netCarbonCredits; - parameter woodSupplyRota(wood_type); + scalar totalWoodNetDemand; + parameter woodSupplyP(wood_type); * Production quantities based on smaller area (before unhandledCropArea adjustment applied) totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location)); @@ -375,7 +397,8 @@ $gdxin netCarbonCredits = SUM(location, carbonCredits.L(location)); - woodSupplyRota(wood_type) = sum(location, woodYieldRota.L(wood_type, location) * timberForestArea.L(wood_type, location)); + totalWoodNetDemand = sum(wood_type, demand(wood_type) + exportAmount.L(wood_type) - importAmount.L(wood_type)); + woodSupplyP(wood_type) = sum(location, woodSupply.L(location)) * (demand(wood_type) + exportAmount.L(wood_type) - importAmount.L(wood_type)) / (totalWoodNetDemand + 1E-6); Scalar totalCostsLessLU; @@ -390,6 +413,8 @@ $gdxin parameter landCoverDiff(land_cover); parameter totalCropArea(crop); parameter totalCropAreaRatio(crop); + parameter rotationPeriod(location); totalLandCoverChange(land_cover, land_cover_after) = sum(location, landCoverChange.L(land_cover, land_cover_after, location)); landCoverDiff(land_cover) = sum(location, landCoverArea.L(land_cover, location)) - sum(location, previousLandCoverArea(land_cover, location)); totalCropArea(crop) = sum(location, cropArea.L(crop, location)); + rotationPeriod(location) = 1 / rotationIntensity.L(location); diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java index e91aa2f6..69788430 100644 --- a/src/ac/ed/lurg/InternationalMarket.java +++ b/src/ac/ed/lurg/InternationalMarket.java @@ -183,7 +183,7 @@ public class InternationalMarket { for (AbstractCountryAgent ca : countryAgents) { WoodUsageData woodUsage = ca.getWoodUsageData().get(woodType); - totalWoodProduction += woodUsage.getHarvest(); + totalWoodProduction += woodUsage.getProduction(); double netImport = woodUsage.getNetImport(); if (netImport >= 0) { totalWoodImport += netImport; diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java index b53c18e5..f84bc034 100755 --- a/src/ac/ed/lurg/ModelConfig.java +++ b/src/ac/ed/lurg/ModelConfig.java @@ -366,6 +366,7 @@ public class ModelConfig { public static final double ANNUAL_MAX_IMPORT_CHANGE = getDoubleProperty("ANNUAL_MAX_IMPORT_CHANGE", 0.05); public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", ANNUAL_MAX_IMPORT_CHANGE*TIMESTEP_SIZE); public static final double IMPORT_CHANGE_ELASTICITY = getDoubleProperty("IMPORT_CHANGE_ELASTICITY", 3.0); // higher value = steeper trade imbalance adjustment + public static final double TRADE_ADJUSTMENT_COST_RATE = getDoubleProperty("TRADE_ADJUSTMENT_COST_RATE", 0.05); // Price caps public static final String PRICE_CAP_FILE = DATA_DIR + File.separator + "price_caps.csv"; @@ -426,7 +427,7 @@ public class ModelConfig { // public static final double BIOENERGY_HEATING_VALUE_GJ_PER_T = getDoubleProperty("BIOENERGY_HEATING_VALUE_GJ_PER_T", 17.5); // GJ per t DM public static final double MARKET_LAMBDA = getDoubleProperty("MARKET_LAMBDA", 0.04); // controls international market price adjustment rate - public static final double MARKET_DAMPING = getDoubleProperty("MARKET_DAMPING", 0.4); // controls international market price adjustment rate + public static final double MARKET_DAMPING = getDoubleProperty("MARKET_DAMPING", 1.0); // controls international market price adjustment rate public static final String PRICE_UPDATE_METHOD = getProperty("PRICE_UPDATE_METHOD", "StockUseRatio"); // Options: StockUseRatio, MarketImbalance, StockChange public static final double MAX_PRICE_INCREASE = getDoubleProperty("MAX_PRICE_INCREASE", 1.5); @@ -517,17 +518,17 @@ public class ModelConfig { public static final boolean IS_FORESTRY_ON = getBooleanProperty("IS_FORESTRY_ON", false); public static final String WOOD_DEMAND_FILENAME = getProperty("WOOD_DEMAND_FILENAME", "wood_demand.csv"); public static final String WOOD_DEMAND_FILE = getProperty("WOOD_DEMAND_FILE", DATA_DIR + File.separator + WOOD_DEMAND_FILENAME); - public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 1500.0); // MtC-eq + public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 420.0); // MtC-eq public static final double IND_ROUNDWOOD_DEMAND_ELASTICITY = getDoubleProperty("IND_ROUNDWOOD_DEMAND_ELASTICITY", 0.3123881); // fitted to FAO data public static final double FUELWOOD_DEMAND_ELASTICITY = getDoubleProperty("FUELDWOOD_DEMAND_ELASTICITY", -0.3598551); public static final double WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 0.3); // m3 to tC-eq p.16 [https://doi.org/10.5194/gmd-13-5425-2020] - public static final double VEGETATION_CLEARING_COST = getDoubleProperty("VEGETATION_CLEARING_COST", 0.001); //$1000/tC - public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_ESTABLISHMENT_COST", 5.0) : 0.0; // establishment, management etc. $1000/ha - public static final double FOREST_BASE_COST = getDoubleProperty("FOREST_BASE_COST", 0.05); + public static final double VEGETATION_CLEARING_COST = getDoubleProperty("VEGETATION_CLEARING_COST", 0.015); //$1000/tC + public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_MANAGEMENT_COST", 2.0) : 0.0; // establishment, management etc. $1000/ha + public static final double FOREST_BASE_COST = getDoubleProperty("FOREST_BASE_COST", 0.02); // $1000/ha public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.2); //$1000/tC - public static final double INIT_ROUNDWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_ROUNDWOOD_PRICE", 0.1) : 0.0; // $1000/tC-eq - public static final double INIT_FUELWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_FUELWOOD_PRICE", 0.05) : 0.0; // $1000/tC-eq - public static final int MAX_FOREST_ROTATION_PERIOD = getIntProperty("CARBON_WOOD_MAX_TIME", 160); // upper data limit, years + public static final double INIT_ROUNDWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_ROUNDWOOD_PRICE", 0.03) : 0.0; // $1000/tC-eq + public static final double INIT_FUELWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_FUELWOOD_PRICE", 0.02) : 0.0; // $1000/tC-eq + public static final int CARBON_WOOD_MAX_TIME = getIntProperty("CARBON_WOOD_MAX_TIME", 160); // upper data limit, years public static final int LAND_AGE_GROUP_SIZE = getIntProperty("LAND_AGE_GROUP_SIZE", 5); // grouping size for land cover age // Carbon diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java index b17122b7..d275921a 100644 --- a/src/ac/ed/lurg/ModelMain.java +++ b/src/ac/ed/lurg/ModelMain.java @@ -1,6 +1,5 @@ package ac.ed.lurg; -import ac.ed.lurg.carbon.CarbonDataOutputer; import ac.ed.lurg.carbon.CarbonFluxItem; import ac.ed.lurg.carbon.CarbonFluxRasterSet; import ac.ed.lurg.carbon.CarbonFluxReader; @@ -369,7 +368,7 @@ public class ModelMain { if (woodUsage == null) continue; - double prod = woodUsage.getHarvest(); + double prod = woodUsage.getProduction(); double netImports = woodUsage.getNetImport(); CountryPrice px = country.getCurrentCountryWoodPrices().get(woodType); @@ -523,10 +522,7 @@ public class ModelMain { LogWriter.println("Outputing forestry data Year: " + timestep.getYear()); ForestryDataOutputer forestryOutputer = new ForestryDataOutputer(timestep.getYear(), globalLandUseRaster, woodYieldData); forestryOutputer.writeOutput(); - - LogWriter.println("Outputing carbon data Year: " + timestep.getYear()); - CarbonDataOutputer carbonOutputer = new CarbonDataOutputer(timestep.getYear(), globalLandUseRaster, carbonFluxData); - carbonOutputer.writeOutput(); + // don't really need this a LPJ outputs have same data, although in a slightly different format // outputLandCover(timestep.getYear(), landUseRaster, LandCoverType.CROPLAND); @@ -735,14 +731,11 @@ public class ModelMain { private void serializeLandUse(RasterSet<LandUseItem> landUseRaster) { RasterSet<LandUseItem> rasterToSerialise = new RasterSet<LandUseItem>(desiredProjection); String landUseFileStr = ModelConfig.IS_CALIBRATION_RUN ? ModelConfig.SERIALIZED_LAND_USE_FILE : ModelConfig.CHECKPOINT_LAND_USE_FILE; - String landCoverFileStr = ModelConfig.IS_CALIBRATION_RUN ? ModelConfig.SERIALIZED_LAND_COVER_FILE : ModelConfig.CHECKPOINT_LAND_COVER_FILE; for (Map.Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) { LandUseItem newLuItem = new LandUseItem(entry.getValue()); - newLuItem.overwriteLandCoverAreas(new LandCoverTile()); // Clear data - Land cover will be serialized separately due to size rasterToSerialise.put(entry.getKey(), newLuItem); } try { - LogWriter.println("Starting serializing LandUse to " + landUseFileStr); FileOutputStream fileOut = new FileOutputStream(landUseFileStr); ObjectOutputStream out = new ObjectOutputStream(fileOut); @@ -753,11 +746,6 @@ public class ModelMain { } catch (IOException i) { i.printStackTrace(); } - - LogWriter.println("Starting serializing LandCover to " + landCoverFileStr); - LandUseBinarySerializer luSer = new LandUseBinarySerializer(); - luSer.serializeLandUse(globalLandUseRaster, landCoverFileStr); // Land cover - LogWriter.println("LandCover data is saved"); } @SuppressWarnings("unchecked") @@ -770,8 +758,6 @@ public class ModelMain { in.close(); fileIn.close(); LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_LAND_USE_FILE); - LandUseBinarySerializer luSer = new LandUseBinarySerializer(); // Land cover - luSer.deserializeLandUse(initLU, ModelConfig.SERIALIZED_LAND_COVER_FILE); return initLU; } catch (IOException i) { LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_LAND_USE_FILE); @@ -819,7 +805,6 @@ public class ModelMain { new MaxCropAreaReader(initialLC).getRasterDataFromFile(ModelConfig.HIGH_SLOPE_AREAS_FILE); // Fraction unavailable for conversion new LandCoverReader(initialLC).getRasterDataFromFile(ModelConfig.INITAL_LAND_COVER_FILE); // Land cover fractions - new LandTileReader(initialLC).getRasterDataFromFile(ModelConfig.LAND_COVER_AGE_DIST_FILENAME); // Age distribution of land cover new InitProtectedAreasReader(initialLC).getRasterDataFromFile(ModelConfig.PROTECTED_AREAS_FILE); // Protected fraction new CropFractionReader(initialLC).getRasterDataFromFile(ModelConfig.CROP_FRACTIONS_FILE); diff --git a/src/ac/ed/lurg/carbon/CarbonDataOutputer.java b/src/ac/ed/lurg/carbon/CarbonDataOutputer.java deleted file mode 100644 index 2f75c344..00000000 --- a/src/ac/ed/lurg/carbon/CarbonDataOutputer.java +++ /dev/null @@ -1,119 +0,0 @@ -package ac.ed.lurg.carbon; - -import java.io.BufferedWriter; -import java.io.File; -import java.io.FileWriter; -import java.io.IOException; -import java.util.Map; -import java.util.Map.Entry; - -import ac.ed.lurg.landuse.LandUseItem; -import ac.ed.lurg.landuse.LccKey; -import ac.ed.lurg.output.AbstractLandUseOutputer; -import ac.ed.lurg.types.LandCoverType; -import ac.ed.lurg.utils.LogWriter; -import ac.sac.raster.RasterKey; -import ac.sac.raster.RasterSet; - -public class CarbonDataOutputer extends AbstractLandUseOutputer { - private CarbonFluxRasterSet carbonFluxData; - - public CarbonDataOutputer(int year, RasterSet<LandUseItem> landUseRaster, CarbonFluxRasterSet carbonFluxData) { - super(year, landUseRaster); - this.carbonFluxData = carbonFluxData; - } - - @Override - public void writeOutput() { - File outputDir = getOutputDir(year); - String outputFileName = outputDir.getPath() + File.separator + "Carbon.txt"; - BufferedWriter outputWriter = null; - try { - outputWriter = new BufferedWriter(new FileWriter(outputFileName, false)); - outputWriter.write("Lon Lat lccCropToPast lccCropToForT lccCropToForC lccCropToNat lccPastToCrop lccPastToForT lccPastToForC lccPastToNat lccForTToCrop lccForTToPast lccForTToForC lccForTToNat lccForCToCrop lccForCToPast lccForCToForT lccForCToNat lccNatToCrop lccNatToPast lccNatToForT lccNatToForC neeCrop neePast neeForT neeForC neeNat"); - outputWriter.newLine(); - - for (Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) { - RasterKey key = entry.getKey(); - LandUseItem luItem = entry.getValue(); - if (luItem == null) - continue; - CarbonFluxItem cfItem = carbonFluxData.get(key); - if (carbonFluxData == null) - continue; - - double lat = landUseRaster.getXCoordin(key); - double lon = landUseRaster.getYCoordin(key); - - double lccCropToPast = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.NATURAL); - double lccCropToForT = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST); - double lccCropToForC = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.CROPLAND); - double lccCropToNat = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.NATURAL); - double lccPastToCrop = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.CROPLAND); - double lccPastToForT = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST); - double lccPastToForC = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.CROPLAND); - double lccPastToNat = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.NATURAL); - double lccForTToCrop = getCarbonFlux(key, LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND); - double lccForTToPast = getCarbonFlux(key, LandCoverType.TIMBER_FOREST, LandCoverType.NATURAL); - double lccForTToForC = getCarbonFlux(key, LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND); - double lccForTToNat = getCarbonFlux(key, LandCoverType.TIMBER_FOREST, LandCoverType.NATURAL); - double lccForCToCrop = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.CROPLAND); - double lccForCToPast = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.NATURAL); - double lccForCToForT = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST); - double lccForCToNat = getCarbonFlux(key, LandCoverType.CROPLAND, LandCoverType.NATURAL); - double lccNatToCrop = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.CROPLAND); - double lccNatToPast = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.NATURAL); - double lccNatToForT = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST); - double lccNatToForC = getCarbonFlux(key, LandCoverType.NATURAL, LandCoverType.CROPLAND); - double neeCrop = getNeeFlux(key, LandCoverType.CROPLAND); - double neePast = getNeeFlux(key, LandCoverType.PASTURE); - double neeForT = getNeeFlux(key, LandCoverType.TIMBER_FOREST); - double neeForC = getNeeFlux(key, LandCoverType.CARBON_FOREST); - double neeNat = getNeeFlux(key, LandCoverType.NATURAL); - - outputWriter.write(String.format("%.2f %.2f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f %.14f", lat, lon, - lccCropToPast, lccCropToForT, lccCropToForC, lccCropToNat, lccPastToCrop, lccPastToForT, lccPastToForC, - lccPastToNat, lccForTToCrop, lccForTToPast, lccForTToForC, lccForTToNat, lccForCToCrop, lccForCToPast, - lccForCToForT, lccForCToNat, lccNatToCrop, lccNatToPast, lccNatToForT, lccNatToForC, - neeCrop, neePast, neeForT, neeForC, neeNat)); - outputWriter.newLine(); - } - } - catch (IOException e) { - LogWriter.print(e); - } - finally { - if (outputWriter != null) { - try { - outputWriter.close(); - } - catch (IOException e) { - LogWriter.print(e); - } - } - } - } - - private double getCarbonFlux(RasterKey key, LandCoverType fromLc, LandCoverType toLc) { - LandUseItem luItem = landUseRaster.get(key); - CarbonFluxItem cfItem = carbonFluxData.get(key); - if (cfItem == null) - return 0.0; - Map<LccKey, Double> lccMap = luItem.getLandCoverChanges(); - LccKey lccKey = new LccKey(fromLc, toLc); - double area = lccMap.getOrDefault(lccKey, 0.0); - double cFlux = cfItem.getCarbonFluxRate(lccKey); - return area * cFlux; - } - - private double getNeeFlux(RasterKey key, LandCoverType lc) { - LandUseItem luItem = landUseRaster.get(key); - CarbonFluxItem cfItem = carbonFluxData.get(key); - if (cfItem == null) - return 0.0; - LccKey lccKey = new LccKey(lc, lc); - double area = luItem.getLandCoverArea(lc); - double cFlux = cfItem.getCarbonFluxRate(lccKey); - return area * cFlux; - } -} diff --git a/src/ac/ed/lurg/carbon/CarbonFluxItem.java b/src/ac/ed/lurg/carbon/CarbonFluxItem.java index cc164a0e..f8602297 100644 --- a/src/ac/ed/lurg/carbon/CarbonFluxItem.java +++ b/src/ac/ed/lurg/carbon/CarbonFluxItem.java @@ -14,85 +14,14 @@ import ac.ed.lurg.types.LandProtectionType; public class CarbonFluxItem implements RasterItem, Serializable { private static final long serialVersionUID = 440720456140537815L; - - // X to Y conversion (new land cover) = conversion flux + ecosystem flux. X to X conversion (existing land cover) = ecosystem flux - Map<LccKey, Double> carbonFluxRate = new HashMap<LccKey, Double>(); Map<LccKey, Double> carbonCreditRate = new HashMap<LccKey, Double>(); - - public void calcConversionCarbonFlux(LccKey key, Double[] cFluxes, LandCoverTile landTile, int maxAge) { - double totalArea = landTile.getTotalArea(key.getFromLc(), LandProtectionType.CONVERTIBLE); - - if (totalArea == 0) { - carbonFluxRate.put(key, 0.0); - carbonCreditRate.put(key, 0.0); - } else { - double totalFlux = 0; - for (int age : landTile.getAgeKeys()) { - int ageCapped = Math.min(age, maxAge - 1); - totalFlux += cFluxes[ageCapped] * landTile.getArea(key.getFromLc(), LandProtectionType.CONVERTIBLE, age); - } - double meanFlux = totalFlux / totalArea; - carbonFluxRate.put(key, meanFlux); - carbonCreditRate.put(key, -meanFlux); - } - } - - public void calcEcosystemCarbonFlux(LandCoverType lcType, Double[] cFluxes, LandCoverTile landTile, int maxAge) { - LccKey key = new LccKey(lcType, lcType); - - double totalArea = landTile.getTotalArea(lcType, LandProtectionType.CONVERTIBLE); - - if (totalArea == 0) { - carbonFluxRate.put(key, 0.0); - } else { - // Flux rate - double totalFlux = 0; - for (int age : landTile.getAgeKeys()) { - int ageCapped = Math.min(age, maxAge - 1); - totalFlux += cFluxes[ageCapped] * landTile.getArea(lcType, LandProtectionType.CONVERTIBLE, age); - } - double meanFlux = totalFlux / totalArea; - carbonFluxRate.put(key, meanFlux); - } - // Potential carbon flux from this existing land cover - double fromLcCumulFlux = 0; - for (int age : landTile.getAgeKeys()) { - for (int i = age; i <= age + ModelConfig.CARBON_HORIZON; i++) { - fromLcCumulFlux += cFluxes[i]; - } - } + + public void calcCarbonCreditRate(Double[] cFluxes) { - // Potential carbon flux if converting to this land cover - double toLcCumulFlux = 0; - for (int age = 0; age <= ModelConfig.CARBON_HORIZON; age++) { - toLcCumulFlux += cFluxes[age]; - } - - // flip signs as we want positive to mean net offset - fromLcCumulFlux *= -1; - toLcCumulFlux *= -1; - - for (LandCoverType lc : LandCoverType.getConvertibleTypes()) { - LccKey fromKey = new LccKey(lcType, lc); - LccKey toKey = new LccKey(lc, lcType); - if (lcType.equals(lc)) { // no net difference if land cover stays same - carbonCreditRate.put(fromKey, 0.0); - } else { - carbonCreditRate.merge(fromKey, -fromLcCumulFlux, Double::sum); - carbonCreditRate.merge(toKey, toLcCumulFlux, Double::sum); - } - } } - - public void setCarbonFluxRate(LccKey key, double cFlux) { - carbonFluxRate.put(key, cFlux); - } - - public void setCarbonFluxRate(LandCoverType fromLc, LandCoverType toLc, double cFlux) { - carbonFluxRate.put(LccKey.createKey(fromLc, toLc), cFlux); - } + public void setCarbonCreditRate(LccKey lcType, double cFlux) { carbonCreditRate.put(lcType, cFlux); @@ -101,23 +30,12 @@ public class CarbonFluxItem implements RasterItem, Serializable { public void setCarbonCreditRate(LandCoverType fromLc, LandCoverType toLc, double cFlux) { carbonCreditRate.put(LccKey.createKey(fromLc, toLc), cFlux); } - - public double getCarbonFluxRate(LccKey key) { - return carbonFluxRate.getOrDefault(key, Double.NaN); - } - - public double getCarbonFluxRate(LandCoverType fromLc, LandCoverType toLc) { - return carbonFluxRate.getOrDefault(LccKey.createKey(fromLc, toLc), Double.NaN); - } + public double getCarbonCreditRate(LandCoverType fromLc, LandCoverType toLc) { return carbonCreditRate.getOrDefault(LccKey.createKey(fromLc, toLc), Double.NaN); } - public Map<LccKey, Double> getCarbonFluxRateMap() { - return carbonFluxRate; - } - public Map<LccKey, Double> getCarbonCreditRateMap() { return carbonCreditRate; } @@ -127,7 +45,6 @@ public class CarbonFluxItem implements RasterItem, Serializable { for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) { for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) { LccKey key = new LccKey(fromLc, toLc); - item.setCarbonFluxRate(key, 0.0); item.setCarbonCreditRate(key, 0.0); } } diff --git a/src/ac/ed/lurg/carbon/CarbonFluxReader.java b/src/ac/ed/lurg/carbon/CarbonFluxReader.java index 8db05960..ac13d382 100644 --- a/src/ac/ed/lurg/carbon/CarbonFluxReader.java +++ b/src/ac/ed/lurg/carbon/CarbonFluxReader.java @@ -108,11 +108,7 @@ public class CarbonFluxReader { if (landUseRaster.containsKey(key)) { CarbonFluxItem cfItem = carbonFluxData.get(key.getCol(), key.getRow()); LandUseItem luItem = landUseRaster.get(key); - if (neeCalc) { - cfItem.calcEcosystemCarbonFlux(lccKey.getFromLc(), map.get(key), luItem.getLandCoverTiles(), maxAge); - } else { - cfItem.calcConversionCarbonFlux(lccKey, map.get(key), luItem.getLandCoverTiles(), maxAge); - } + cfItem.calcCarbonCreditRate(map.get(key)); } } LogWriter.println("Reading " + fileName + " took " + (System.currentTimeMillis() - startTime) + " ms"); diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java index ded725a0..2f66cef1 100644 --- a/src/ac/ed/lurg/country/CountryAgent.java +++ b/src/ac/ed/lurg/country/CountryAgent.java @@ -141,9 +141,6 @@ public class CountryAgent extends AbstractCountryAgent { GamsRasterOutput result = opti.run(); previousGamsRasterOutput = result; - - if (!ModelConfig.IS_CALIBRATION_RUN) - incrementLandCoverAge(); return result; } @@ -256,12 +253,15 @@ public class CountryAgent extends AbstractCountryAgent { double changeUp; double changeDown; - double maxOfProdOrSupply = woodUsage.getHarvest() + Math.max(baseTrade, 0); + double demand = currentWoodDemand.get(woodType); + double maxOfProdOrSupply = woodUsage.getProduction() + Math.max(baseTrade, 0); + //double baseTrade = demand * woodUsage.getBaseImportDemandRatio(); if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) { changeUp = changeDown = 0.0; } else { changeUp = changeDown = maxOfProdOrSupply * ModelConfig.MAX_IMPORT_CHANGE; + //changeUp = changeDown = demand * 0.04; } woodTradeConstraints.put(woodType, new TradeConstraint(baseTrade - changeDown, baseTrade + changeUp)); @@ -331,12 +331,6 @@ public class CountryAgent extends AbstractCountryAgent { } } - - private void incrementLandCoverAge() { - for (LandUseItem luItem : previousGamsRasterOutput.getLandUses().values()) { - luItem.incrementTilesAge(); - } - } public CarbonUsageData getCarbonUsageData() { return previousGamsRasterOutput.getCarbonUsageData(); diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java index 3573ab1e..24311649 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 prevRotaIntP = inDB.addParameter("previousRotationIntensity", 1); double totalAgriLand = 0; double totalSuitable = 0; @@ -194,7 +195,10 @@ public class GamsLocationOptimiser { if (area == 0) continue; setGamsParamValueTruncate(prevLandCoverP.addRecord(v), area, 6); - } + } + + double foo = landUseItem.getForestRotationIntensity(); + setGamsParamValue(prevRotaIntP.addRecord(locString), landUseItem.getForestRotationIntensity(), 6); } if (DEBUG) LogWriter.println(String.format(" Total agricultural %.1f,\t suitable %.1f", totalAgriLand, totalSuitable)); @@ -355,6 +359,8 @@ public class GamsLocationOptimiser { && inputData.getTimestep().isInitialTimestep()) ? 1 : ModelConfig.MAX_FERT_CHANGE, 5); addScalar(inDB, "maxIrrigChange", (ModelConfig.IS_CALIBRATION_RUN && inputData.getTimestep().isInitialTimestep()) ? 1 : ModelConfig.MAX_IRRIG_CHANGE, 5); + + addScalar(inDB, "tradeAdjustmentCostRate", ModelConfig.TRADE_ADJUSTMENT_COST_RATE, 5); // Forestry GAMSParameter maxRotaIntensityP = inDB.addParameter("maxRotationIntensity", 1); @@ -377,16 +383,24 @@ public class GamsLocationOptimiser { setGamsParamValue(maxRotaIntensityP.addRecord(woodType.getName()), 1.0 / woodType.getMinRotation(), 6); } - addScalar(inDB, "forestManagementCost", ModelConfig.FOREST_MANAGEMENT_COST, 5); - addScalar(inDB, "vegClearingCostRate", ModelConfig.VEGETATION_CLEARING_COST, 5); + addScalar(inDB, "forestManagementCost", ModelConfig.FOREST_MANAGEMENT_COST, -1); + addScalar(inDB, "vegClearingCostRate", ModelConfig.VEGETATION_CLEARING_COST, -1); GAMSParameter forestBaseCostP = inDB.addParameter("forestBaseCost", 1); // No base cost if forestry off to stop unnecessary conversion. setGamsParamValue(forestBaseCostP.addRecord(LandCoverType.TIMBER_FOREST.getName()), - ModelConfig.IS_FORESTRY_ON ? ModelConfig.FOREST_BASE_COST : 0, 1); + ModelConfig.IS_FORESTRY_ON ? ModelConfig.FOREST_BASE_COST : 0, -1); // Keeping base cost even if forestry & carbon off to stop conversion to carbon forest setGamsParamValue(forestBaseCostP.addRecord(LandCoverType.CARBON_FOREST.getName()), - ModelConfig.FOREST_BASE_COST, 1); + ModelConfig.FOREST_BASE_COST, -1); + + for (WoodType wType : WoodType.values()) { + double netImport = countryInput.getPreviousWoodUsageData().get(wType).getNetImport(); + double imports = netImport > 0 ? netImport : 0; + double exports = netImport < 0 ? -netImport : 0; + setGamsParamValue(previousExportAmountP.addRecord(wType.getName()), exports, -1); + setGamsParamValue(previousImportAmountP.addRecord(wType.getName()), imports, -1); + } // Yield from timber forest rotation Map<Integer, ? extends WoodYieldData> woodYieldData = inputData.getWoodYields(); @@ -436,7 +450,6 @@ public class GamsLocationOptimiser { v.add(fromLc.getName()); v.add(toLc.getName()); v.add(locString); - setGamsParamValue(cFluxRateP.addRecord(v), cFlux.getCarbonFluxRate(fromLc, toLc), 5); if (ModelConfig.getCarbonOption(fromLc, toLc)) { setGamsParamValue(cFluxCreditRateP.addRecord(v), cFlux.getCarbonCreditRate(fromLc, toLc), 5); @@ -610,32 +623,33 @@ public class GamsLocationOptimiser { // Wood production Map<WoodType, WoodUsageData> newWoodUsageMap = new HashMap<>(); - GAMSParameter parWoodSupplyRota = outDB.getParameter("woodSupplyRota"); - for (WoodType w : WoodType.values()) { - double supply = getParmValue(parWoodSupplyRota, w.getName()); - double netImports = getParmValue(parmNetImports, w.getName()); - newWoodUsageMap.put(w, new WoodUsageData(supply, netImports)); + /* + GAMSVariable varWoodSupplyRota = outDB.getVariable("woodSupply"); + for (GAMSVariableRecord rec : varWoodSupplyRota) { + WoodType wType = WoodType.getForName(rec.getKey(0)); + double supply = rec.getLevel(); + double netImports = getParmValue(parmNetImports, wType.getName()); + newWoodUsageMap.put(wType, new WoodUsageData(supply, netImports)); + } + */ + GAMSParameter woodSupplyP = outDB.getParameter("woodSupplyP"); + for (WoodType wType : WoodType.values()) { + double netImports = getParmValue(parmNetImports, wType.getName()); + double supply = getParmValue(woodSupplyP, wType.getName()); + WoodUsageData woodUsageData = inputData.getCountryInput().getPreviousWoodUsageData().get(wType); + woodUsageData.setProduction(supply); + woodUsageData.setNetImport(netImports); + newWoodUsageMap.put(wType, woodUsageData); } - Map<Integer, Map<WoodType, Double>> rotationPeriods = new HashMap<>(); + Map<Integer, Double> rotationPeriods = new HashMap<>(); GAMSVariable varRota = outDB.getVariable("rotationIntensity"); for (GAMSVariableRecord rec : varRota) { - WoodType woodType = WoodType.getForName(rec.getKey(0)); - int loc = Integer.parseInt(rec.getKey(1)); + int loc = Integer.parseInt(rec.getKey(0)); double intensity = rec.getLevel(); - Map<WoodType, Double> map = rotationPeriods.computeIfAbsent(loc, k -> new HashMap<>()); - map.put(woodType, 1 / intensity); + rotationPeriods.put(loc, intensity); } - Map<Integer, Map<WoodType, Double>> forestAreaMap = new HashMap<>(); - GAMSVariable varForestArea = outDB.getVariable("timberForestArea"); - for (GAMSVariableRecord rec : varForestArea) { - WoodType woodType = WoodType.getForName(rec.getKey(0)); - int loc = Integer.parseInt(rec.getKey(1)); - double forestArea = rec.getLevel(); - Map<WoodType, Double> map = forestAreaMap.computeIfAbsent(loc, k -> new HashMap<>()); - map.put(woodType, forestArea); - } // Carbon double carbonSequestered = outDB.getParameter("netCarbonCredits").getFirstRecord().getValue(); @@ -643,7 +657,7 @@ public class GamsLocationOptimiser { CarbonUsageData carbonUsageData = new CarbonUsageData(carbonSequestered, netCarbonImport, 0.0); return new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, carbonUsageData, - newWoodUsageMap, rotationPeriods, forestAreaMap); + newWoodUsageMap, rotationPeriods); } private double getParmValue(GAMSParameter aParm, String itemName) { diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java index 4b3bce17..a6d8b4b8 100644 --- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java +++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java @@ -22,16 +22,14 @@ public class GamsLocationOutput { Map<Integer, List<LandCoverChangeItem>> landCoverChanges; private CarbonUsageData carbonUsageData; private Map<WoodType, WoodUsageData> woodUsageData; - private Map<Integer, Map<WoodType, Double>> rotationPeriods; - private Map<Integer, Map<WoodType, Double>> forestAreas; + private Map<Integer, Double> rotationPeriods; public GamsLocationOutput(ModelStat status, Map<Integer, LandUseItem> landUses, Map<CropType, CropUsageData> cropUsageData, Map<Integer, List<LandCoverChangeItem>> landCoverChanges, CarbonUsageData carbonUsageData, Map<WoodType, WoodUsageData> woodUsageData, - Map<Integer, Map<WoodType, Double>> rotationPeriods, - Map<Integer, Map<WoodType, Double>> forestAreas) { + Map<Integer, Double> rotationPeriods) { super(); this.status = status; this.landUses = landUses; @@ -40,7 +38,6 @@ public class GamsLocationOutput { this.carbonUsageData = carbonUsageData; this.woodUsageData = woodUsageData; this.rotationPeriods = rotationPeriods; - this.forestAreas = forestAreas; } public ModelStat getStatus() { @@ -66,11 +63,8 @@ public class GamsLocationOutput { return woodUsageData; } - public Map<Integer, Map<WoodType, Double>> getRotationPeriods() { + public Map<Integer, Double> getRotationPeriods() { return rotationPeriods; } - public Map<Integer, Map<WoodType, Double>> getForestAreas() { - return forestAreas; - } } diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java index fa726076..f671b3cd 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java @@ -6,8 +6,6 @@ import java.util.Map.Entry; import ac.ed.lurg.ModelConfig; import ac.ed.lurg.carbon.CarbonFluxItem; import ac.ed.lurg.country.LandCoverChangeItem; -import ac.ed.lurg.forestry.ForestManager; -import ac.ed.lurg.forestry.WoodHarvestItem; import ac.ed.lurg.forestry.WoodYieldData; import ac.ed.lurg.forestry.WoodYieldItem; import ac.ed.lurg.landuse.*; @@ -104,61 +102,22 @@ public class GamsRasterOptimiser { Map<RasterKey, List<LandCoverChangeItem>> lcChanges = disaggregateLandCoverChanges(gamsOutput.getLandCoverChanges(), newLandUseRaster); - // Calculate amount of wood produced from land cover change - Map<WoodType, Double> lccWoodAmount = ForestManager.handleLccWood(rasterInputData.getWoodYields(), - lcChanges, rasterInputData.getPreviousLandUses()); - // Allocate land cover changes allocLandCoverChanges(newLandUseRaster, lcChanges); - Map<LandCoverType, Double> foo = new HashMap<>(); - for (LandUseItem i : newLandUseRaster.values()) { - for (LandCoverType c : LandCoverType.getConvertibleTypes()) { - foo.merge(c, i.getLandCoverArea(c, LandProtectionType.CONVERTIBLE), Double::sum); - } - } - - // Set wood type forest fractions - Map<Integer, Map<WoodType, Double>> woodAreas = gamsOutput.getForestAreas(); - for (int locId : woodAreas.keySet()) { - Set<RasterKey> keys = new HashSet<RasterKey>(); - for (Entry<RasterKey, IntegerRasterItem> mapEntry : mapping.entrySet()) { - IntegerRasterItem iri = mapEntry.getValue(); - if (locId == iri.getInt()) - keys.add(mapEntry.getKey()); - } - Map<WoodType, Double> areaMap = woodAreas.get(locId); - double totalArea = areaMap.values().stream().reduce(0.0, Double::sum); - for (WoodType wType : WoodType.values()) { - double fraction = totalArea > 0 ? areaMap.get(wType) / totalArea : 0; - for (RasterKey key : keys) { - LandUseItem luItem = newLandUseRaster.get(key); - luItem.setForestFraction(wType, fraction); - } - } - } - - // Calculate land cover changes for timber forest rotation and clear-cut. Resets harvested stands to age 0 - Map<RasterKey, WoodHarvestItem> rotaHarvestMap = ForestManager.handleRotaWood(newLandUseRaster, - rasterInputData.getWoodYields(), gamsOutput, mapping); - - // Set forest rotations - for (RasterKey key : rotaHarvestMap.keySet()) { - WoodHarvestItem whItem = rotaHarvestMap.get(key); + // Forest rotation intensity + for (RasterKey key : newLandUseRaster.keySet()) { LandUseItem luItem = newLandUseRaster.get(key); - luItem.setWoodHarvests(whItem); - } - - - - // Add wood produced from LCC to output - Map<WoodType, WoodUsageData> woodUsageMap = gamsOutput.getWoodUsageData(); - for (WoodType wType : lccWoodAmount.keySet()) { - double amount = lccWoodAmount.get(wType); - WoodUsageData woodUsage = woodUsageMap.get(wType); - woodUsage.setLucHarvest(amount); + Integer locId = mapping.get(key).getInt(); + Double intensity = gamsOutput.getRotationPeriods().get(locId); + if (intensity == null) { + luItem.setForestRotationIntensity(0); + } else { + luItem.setForestRotationIntensity(intensity); + } } + // Crop areas and intensity for (Map.Entry<Integer, LandUseItem> entry : gamsOutput.getLandUses().entrySet()) { Integer locId = entry.getKey(); LandUseItem newLandUseAggItem = entry.getValue(); @@ -318,16 +277,6 @@ public class GamsRasterOptimiser { // Aggregate carbon fluxes and wood yields if (ModelConfig.IS_FORESTRY_ON) { - // wood from LUC - for (LccKey lccKey : woodYieldItem.getYieldLucMap().keySet()) { - double areaThisTime = landUseItem.getLandCoverArea(lccKey.getFromLc(), LandProtectionType.CONVERTIBLE); - double areaSoFar = aggLandUse.getLandCoverArea(lccKey.getFromLc(), LandProtectionType.CONVERTIBLE); - double wYieldThisTime = woodYieldItem.getAggregatedYieldLuc(lccKey, landUseItem); - double wYieldSoFar = aggWYield.getYieldLuc(lccKey); - double wYieldAgg = aggregateMean(wYieldSoFar, areaSoFar, wYieldThisTime, areaThisTime); - aggWYield.setYieldLuc(lccKey, wYieldAgg); - } - for (Map.Entry<Integer, Double> entry : woodYieldItem.getYieldResponseMap().entrySet()) { double areaThisTime = landUseItem.getLandCoverArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE); double areaSoFar = aggLandUse.getLandCoverArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE); @@ -345,11 +294,6 @@ public class GamsRasterOptimiser { double areaThisTime = landUseItem.getLandCoverArea(fromLc); double areaSoFar = aggLandUse.getLandCoverArea(toLc); - double cFluxThisTime = carbonFluxItem.getCarbonFluxRate(fromLc, toLc); - double cFluxSoFar = aggCFlux.getCarbonFluxRate(fromLc, toLc); - double cFluxAgg = aggregateMean(cFluxSoFar, areaSoFar, cFluxThisTime, areaThisTime); - aggCFlux.setCarbonFluxRate(fromLc, toLc, cFluxAgg); - double cCreditThisTime = carbonFluxItem.getCarbonCreditRate(fromLc, toLc); double cCreditSoFar = aggCFlux.getCarbonCreditRate(fromLc, toLc); double cCreditAgg = aggregateMean(cCreditSoFar, areaSoFar, cCreditThisTime, areaThisTime); @@ -408,7 +352,9 @@ public class GamsRasterOptimiser { double protAreaThisTime = landUseItem.getLandCoverArea(lcType, LandProtectionType.PROTECTED); double protAreaSoFar = aggLandUse.getLandCoverArea(lcType, LandProtectionType.PROTECTED); aggLandUse.setLandCoverArea(lcType, LandProtectionType.PROTECTED, protAreaSoFar + protAreaThisTime); - } + } + + aggLandUse.setForestRotationIntensity(landUseItem.getForestRotationIntensity()); } checkedTotalAreas(landUseRaster, "before"); diff --git a/src/ac/ed/lurg/demand/AbstractDemandManager.java b/src/ac/ed/lurg/demand/AbstractDemandManager.java index 8abebcc9..dc7d6d27 100644 --- a/src/ac/ed/lurg/demand/AbstractDemandManager.java +++ b/src/ac/ed/lurg/demand/AbstractDemandManager.java @@ -19,8 +19,8 @@ public abstract class AbstractDemandManager implements Serializable { protected transient CalorieManager calorieManager; protected transient BioenergyDemandManager bioenergyDemandManager; protected transient CerealFractionsManager cerealFractionsManager; - protected CarbonDemandManager carbonDemandManager; - protected WoodDemandManager woodDemandManager; + protected transient CarbonDemandManager carbonDemandManager; + protected transient WoodDemandManager woodDemandManager; public AbstractDemandManager(CalorieManager calorieManager) { setup(calorieManager); diff --git a/src/ac/ed/lurg/demand/ElasticDemandManager.java b/src/ac/ed/lurg/demand/ElasticDemandManager.java index aa3e634e..0ca70b21 100755 --- a/src/ac/ed/lurg/demand/ElasticDemandManager.java +++ b/src/ac/ed/lurg/demand/ElasticDemandManager.java @@ -23,7 +23,7 @@ public class ElasticDemandManager extends AbstractSSPDemandManager { private transient Map<SingleCountry, GamsDemandOutput> previousGamsDemands; private transient BufferedWriter outputFile; private PriceMarkUps priceMarkups; - private int baseYearToRebaseConsumption; + private Integer baseYearToRebaseConsumption; /** Constructor used in calibration only */ public ElasticDemandManager(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager) { diff --git a/src/ac/ed/lurg/demand/WoodDemandManager.java b/src/ac/ed/lurg/demand/WoodDemandManager.java index 2ef60933..6017b1e4 100644 --- a/src/ac/ed/lurg/demand/WoodDemandManager.java +++ b/src/ac/ed/lurg/demand/WoodDemandManager.java @@ -14,7 +14,7 @@ import ac.ed.lurg.utils.LogWriter; import ac.ed.lurg.utils.StringTabularReader; public class WoodDemandManager implements Serializable { - @Serial + private static final long serialVersionUID = 171124405614624397L; private Map<SingleCountry, Map<WoodType, Double>> baseWoodDemand = new HashMap<SingleCountry, Map<WoodType, Double>>(); @@ -30,7 +30,7 @@ public class WoodDemandManager implements Serializable { for (Map<String, String> row : rowList) { String countryCode = row.get("Iso3"); String woodTypeStr = row.get("Type"); - double demand = Double.parseDouble(row.get("Demand")) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; + double demand = Double.parseDouble(row.get("Demand")); SingleCountry sc = CountryManager.getInstance().getForCode(countryCode); WoodType woodType = WoodType.getForName(woodTypeStr); Map<WoodType, Double> wMap = baseWoodDemand.computeIfAbsent(sc, k -> new HashMap<WoodType, Double>()); diff --git a/src/ac/ed/lurg/forestry/ForestManager.java b/src/ac/ed/lurg/forestry/ForestManager.java deleted file mode 100644 index 4fbb77c4..00000000 --- a/src/ac/ed/lurg/forestry/ForestManager.java +++ /dev/null @@ -1,178 +0,0 @@ -package ac.ed.lurg.forestry; - -import ac.ed.lurg.ModelConfig; -import ac.ed.lurg.country.LandCoverChangeItem; -import ac.ed.lurg.country.gams.GamsLocationOutput; -import ac.ed.lurg.landuse.LandCoverTile; -import ac.ed.lurg.landuse.LandUseItem; -import ac.ed.lurg.types.LandCoverType; -import ac.ed.lurg.types.LandProtectionType; -import ac.ed.lurg.types.WoodType; -import ac.ed.lurg.utils.LogWriter; -import ac.sac.raster.IntegerRasterItem; -import ac.sac.raster.RasterKey; -import ac.sac.raster.RasterSet; - -import java.util.*; - -public class ForestManager { - - // Calculates amount of wood harvested from land cover change - public static Map<WoodType, Double> handleLccWood(RasterSet<WoodYieldItem> woodYieldRaster, - Map<RasterKey, List<LandCoverChangeItem>> landCoverChanges, - RasterSet<LandUseItem> inputLandUseRaster) { - - Map<WoodType, Double> woodAmount = new HashMap<>(); - for (WoodType w : WoodType.values()) { - woodAmount.put(w, 0.0); - } - - if (!ModelConfig.IS_FORESTRY_ON) { - return woodAmount; - } - - for (Map.Entry<RasterKey, List<LandCoverChangeItem>> entry : landCoverChanges.entrySet()) { - RasterKey key = entry.getKey(); - List<LandCoverChangeItem> changes = entry.getValue(); - WoodYieldItem wyItem = woodYieldRaster.get(key); - LandCoverTile tiles = inputLandUseRaster.get(key).getLandCoverTiles(); - for (LandCoverChangeItem lccItem : changes) { - LandCoverType fromLc = lccItem.getFromLandCover(); - LandCoverType toLc = lccItem.getToLandCover(); - double changedArea = lccItem.getArea(); - double totalArea = tiles.getTotalArea(fromLc, LandProtectionType.CONVERTIBLE); - if (changedArea > totalArea && Math.abs(changedArea - totalArea) > 1e-10) { - LogWriter.printlnError("Error: ForestManager changed area exceeds total area."); - } - changedArea = Math.min(changedArea, totalArea); // deal with floating point errors - double fraction = changedArea / totalArea; - // Assume stands are converted proportionally by area. - // Type of wood harvested depends on stand age. - for (int age : tiles.getAgeKeys()) { - double standArea = tiles.getArea(fromLc, LandProtectionType.CONVERTIBLE, age); - double standYield = wyItem.getYieldLuc(fromLc, toLc, age); - double standAmount = standYield * standArea * fraction; - // Assuming roundwood if above min rotation age otherwise fuelwood - if (age > WoodType.IND_ROUNDWOOD.getMinRotation()) { - woodAmount.merge(WoodType.IND_ROUNDWOOD, standAmount, Double::sum); - } else { - woodAmount.merge(WoodType.FUELWOOD, standAmount, Double::sum); - } - } - } - } - - return woodAmount; - } - - // Calculates area of timber forest cleared to meet demand - public static Map<RasterKey, WoodHarvestItem> handleRotaWood(RasterSet<LandUseItem> outputLandUseRaster, - RasterSet<WoodYieldItem> woodYieldRaster, - GamsLocationOutput locationOutput, - RasterSet<IntegerRasterItem> mapping) { - - // Since some trees could be used for either roundwood or fuelwood, we harvest the more valuable roundwood first - final List<WoodType> harvestOrder = new ArrayList<>(Arrays.asList(WoodType.IND_ROUNDWOOD, WoodType.FUELWOOD)); - - Map<RasterKey, WoodHarvestItem> harvestMap = new HashMap<>(); - - if (!ModelConfig.IS_FORESTRY_ON) { - return harvestMap; - } - - for (WoodType wType : harvestOrder) { - handleRotaWoodForType(outputLandUseRaster, woodYieldRaster, locationOutput, mapping, wType, harvestMap); - } - - return harvestMap; - } - - private static void handleRotaWoodForType(RasterSet<LandUseItem> outputLandUseRaster, - RasterSet<WoodYieldItem> woodYieldRaster, - GamsLocationOutput locationOutput, - RasterSet<IntegerRasterItem> mapping, - WoodType wType, - Map<RasterKey, WoodHarvestItem> harvestMap) { - - double demand = locationOutput.getWoodUsageData().get(wType).getHarvest() - - locationOutput.getWoodUsageData().get(wType).getLucHarvest(); - if (demand <= 0) { - return; - } - - - // First calculate how much wood is available - Map<RasterKey, Double> stockMap = new HashMap<>(); - for (RasterKey key : outputLandUseRaster.keySet()) { - LandUseItem luItem = outputLandUseRaster.get(key); - LandCoverTile tiles = luItem.getLandCoverTiles(); - WoodYieldItem wyItem = woodYieldRaster.get(key); - Integer locId = mapping.get(key).getInt(); - double stock = 0; - double rotationAge = wType.getMinRotation(); - for (int age : tiles.getAgeKeys()) { - if (age >= rotationAge) { - double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age); - if (area == 0) - continue; - double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age); - stock += yield * area; - } - } - stockMap.put(key, stock); - } - - double totalStock = stockMap.values().stream().reduce(0.0, Double::sum); - LogWriter.println("Total stock for " + wType + " = " + totalStock); - if (demand > totalStock) { - LogWriter.printlnWarning("Insufficient " + wType.getName() + " to meet demand. Shortfall: " + (demand - totalStock)); - } - if (totalStock < 1e-10) { - return; - } - double harvestFraction = Math.min(demand / totalStock, 1); - // Now harvest wood to meet demand. - for (RasterKey key : outputLandUseRaster.keySet()) { - LandUseItem luItem = outputLandUseRaster.get(key); - LandCoverTile tiles = luItem.getLandCoverTiles(); - WoodYieldItem wyItem = woodYieldRaster.get(key); - Integer locId = mapping.get(key).getInt(); - double rotationAge = wType.getMinRotation(); - double totalAmountToRemove = stockMap.get(key) * harvestFraction; - if (totalAmountToRemove < 1e-10) { - continue; - } - // Assuming that oldest trees get harvested first - List<Integer> ageList = new ArrayList<>(tiles.getAgeKeys()); - List<Integer> rotationAgeList = new ArrayList<>(); - for (int i : ageList) { - if (i >= rotationAge) { - rotationAgeList.add(i); - } - } - // Harvest in descending age order - rotationAgeList.sort(Collections.reverseOrder()); - - for (int age : rotationAgeList) { - double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age); - double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age); - if (area == 0 || yield == 0) - continue; - double standStock = area * yield; - double amountToRemove = Math.min(totalAmountToRemove, standStock); - double convertedArea = Math.min(area, amountToRemove / yield); - LandCoverChangeItem lccItem = new LandCoverChangeItem(LandCoverType.TIMBER_FOREST, - LandCoverType.TIMBER_FOREST, age, convertedArea); - WoodHarvestItem whItem = harvestMap.computeIfAbsent(key, k -> new WoodHarvestItem()); - whItem.addAreasHarvested(lccItem); - whItem.addWoodAmount(wType, amountToRemove); - luItem.doForestRotation(age, convertedArea); // WARNING! Side effect - totalAmountToRemove -= amountToRemove; - } - - if (totalAmountToRemove > 0) - LogWriter.printlnWarning("ForestManager: Could not meet wood demand."); - } - } - -} diff --git a/src/ac/ed/lurg/forestry/ForestryDataOutputer.java b/src/ac/ed/lurg/forestry/ForestryDataOutputer.java index 41760185..1434f861 100644 --- a/src/ac/ed/lurg/forestry/ForestryDataOutputer.java +++ b/src/ac/ed/lurg/forestry/ForestryDataOutputer.java @@ -2,8 +2,16 @@ package ac.ed.lurg.forestry; import ac.ed.lurg.landuse.LandUseItem; import ac.ed.lurg.output.AbstractLandUseOutputer; +import ac.ed.lurg.utils.LogWriter; +import ac.sac.raster.RasterKey; import ac.sac.raster.RasterSet; +import java.io.BufferedWriter; +import java.io.File; +import java.io.FileWriter; +import java.io.IOException; +import java.util.Map; + public class ForestryDataOutputer extends AbstractLandUseOutputer { private WoodYieldRasterSet woodYieldData; @@ -14,16 +22,16 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer { @Override public void writeOutput() { - /* + File outputDir = getOutputDir(year); String outputFileName = outputDir.getPath() + File.separator + "Forestry.txt"; BufferedWriter outputWriter = null; try { outputWriter = new BufferedWriter(new FileWriter(outputFileName, false)); - outputWriter.write("Lon Lat RotationAge YieldAtRotation RotationCutArea StandingStockTimber StandingStockCarbon StandingStockNatural"); + outputWriter.write("Lon Lat RotaInt"); outputWriter.newLine(); - for (Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) { + for (Map.Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) { RasterKey key = entry.getKey(); LandUseItem luItem = entry.getValue(); if (luItem == null) @@ -34,15 +42,9 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer { double lat = landUseRaster.getXCoordin(key); double lon = landUseRaster.getYCoordin(key); - int rotaAge = wyItem.getOptimalRotation(); - double rotaYield = wyItem.getYieldAtRotation(); - double rotaArea = wyItem.getRotationArea(); - double stockTimber = wyItem.getYield(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND); - double stockCarbon = wyItem.getYield(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND); - double stockNatural = wyItem.getYield(LandCoverType.NATURAL, LandCoverType.CROPLAND); + double rotaInt = luItem.getForestRotationIntensity(); - outputWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f %.14f", lat, lon, - rotaAge, rotaYield, rotaArea, stockTimber, stockCarbon, stockNatural)); + outputWriter.write(String.format("%.2f %.2f %.14f", lat, lon, rotaInt)); outputWriter.newLine(); } @@ -61,6 +63,6 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer { } } - */ + } } diff --git a/src/ac/ed/lurg/forestry/WoodHarvestItem.java b/src/ac/ed/lurg/forestry/WoodHarvestItem.java deleted file mode 100644 index 228c118a..00000000 --- a/src/ac/ed/lurg/forestry/WoodHarvestItem.java +++ /dev/null @@ -1,34 +0,0 @@ -package ac.ed.lurg.forestry; - -import ac.ed.lurg.country.LandCoverChangeItem; -import ac.ed.lurg.types.WoodType; - -import java.io.Serial; -import java.io.Serializable; -import java.util.ArrayList; -import java.util.HashMap; -import java.util.List; -import java.util.Map; - -public class WoodHarvestItem implements Serializable { - @Serial - private static final long serialVersionUID = -79494635617251642L; - private List<LandCoverChangeItem> areasHarvested = new ArrayList<>(); - private Map<WoodType, Double> woodAmount = new HashMap<>(); - - public List<LandCoverChangeItem> getAreasHarvested() { - return areasHarvested; - } - - public void addAreasHarvested(LandCoverChangeItem lccItem) { - areasHarvested.add(lccItem); - } - - public Map<WoodType, Double> getWoodAmount() { - return woodAmount; - } - - public void addWoodAmount(WoodType woodType, double amount) { - woodAmount.merge(woodType, amount, Double::sum); - } -} diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java index 55266f5b..584e4e1f 100644 --- a/src/ac/ed/lurg/forestry/WoodYieldItem.java +++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java @@ -1,5 +1,6 @@ package ac.ed.lurg.forestry; +import java.util.Arrays; import java.util.HashMap; import java.util.Map; import ac.ed.lurg.ModelConfig; @@ -11,134 +12,31 @@ import ac.ed.lurg.types.LandProtectionType; import ac.ed.lurg.types.WoodType; import ac.sac.raster.RasterItem; -public class WoodYieldItem implements RasterItem { - private Map<LccKey, Map<Integer, Double>> yieldLuc = new HashMap<LccKey, Map<Integer, Double>>(); - private Map<LandCoverType, Map<Integer, Double>> yieldRota = new HashMap<LandCoverType, Map<Integer, Double>>(); +public class WoodYieldItem implements RasterItem { private WoodYieldResponse yieldResponse = new WoodYieldResponse(); private double maxYield = 0; // Yield at maximum age. Used for clustering public WoodYieldItem() {} - public void setYieldLuc(LccKey key, Double[] yields, LandCoverTile landUseTiles, int maxAge) { - Map<Integer, Double> yieldMap = yieldLuc.computeIfAbsent(key, k -> new HashMap<Integer, Double>()); - for (Integer age : landUseTiles.getAgeKeys()) { - int ageCapped = Math.min(age, maxAge); - yieldMap.put(ageCapped, yields[ageCapped]); - } - } - - public void setYieldRota(LandCoverType lcType, Double[] yields, LandCoverTile landUseTiles, int maxAge) { - Map<Integer, Double> yieldMap = yieldRota.computeIfAbsent(lcType, k -> new HashMap<Integer, Double>()); - for (Integer age : landUseTiles.getAgeKeys()) { - int ageCapped = Math.min(age, maxAge); - yieldMap.put(ageCapped, yields[ageCapped]); - } - - // For clustering - if (lcType.equals(LandCoverType.TIMBER_FOREST)) { - maxYield = yields[ModelConfig.MAX_FOREST_ROTATION_PERIOD]; - } - } - public void setYieldResponse(Double[] yields, int modelYear) { int dataYear = Math.floorDiv(modelYear - 10, ModelConfig.CARBON_DATA_TIMESTEP_SIZE) * ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 10; for (int age = modelYear - dataYear; age <= 160; age += ModelConfig.CARBON_DATA_TIMESTEP_SIZE) { yieldResponse.setYield(age, yields[age]); } yieldResponse.setYield(0, 0); // in case missing + + maxYield = Math.max(maxYield, Arrays.stream(yields).reduce(Double::max).get()); } public WoodYieldResponse getYieldResponse() { return yieldResponse; } - - public double getYieldLuc(LandCoverType fromLc, LandCoverType toLc, int age) { - return getYieldLuc(new LccKey(fromLc, toLc), age); - } - - public double getYieldLuc(LccKey key, int age) { - if (yieldLuc.containsKey(key)) { - return yieldLuc.get(key).getOrDefault(age, 0.0); - } else { - return 0.0; - } - } - - public double getYieldRota(LandCoverType lcType, int age) { - if (yieldRota.containsKey(lcType)) { - return yieldRota.get(lcType).getOrDefault(age, 0.0); - } else { - return 0.0; - } - } public double getMaxYield() { return maxYield; } - - public Map<LccKey, Map<Integer, Double>> getYieldLucMap() { - return yieldLuc; - } - - // Average yield assuming proportional conversion of forest stands - public double getAggregatedYieldLuc(LccKey key, LandUseItem luItem) { - LandCoverType fromLc = key.getFromLc(); - LandCoverTile tile = luItem.getLandCoverTiles(); - double area = 0; - double stock = 0; - for (int age : tile.getAgeKeys()) { - double a = tile.getArea(fromLc, LandProtectionType.CONVERTIBLE, age); - double yield = getYieldLuc(key, age); - stock += a * yield; - area += a; - } - return area > 0 ? stock / area : 0; - } - // Yield for each wood type - public Map<LccKey, Map<WoodType, Double>> getYieldLucByType(LandUseItem luItem) { - Map<LccKey, Map<WoodType, Double>> yieldByType = new HashMap<>(); - Map<LccKey, Map<WoodType, Double>> stockByType = new HashMap<>(); - Map<LccKey, Map<WoodType, Double>> areaByType = new HashMap<>(); - - for (LandCoverType fromLc : LandCoverType.values()) { - for (LandCoverType toLc : LandCoverType.values()) { - LccKey key = new LccKey(fromLc, toLc); - LandCoverTile tile = luItem.getLandCoverTiles(); - for (int age : tile.getAgeKeys()) { - double area = tile.getArea(fromLc, LandProtectionType.CONVERTIBLE, age); - double yield = getYieldLuc(key, age); - Map<WoodType, Double> stockMap = stockByType.computeIfAbsent(key, k -> new HashMap<>()); - Map<WoodType, Double> areaMap = stockByType.computeIfAbsent(key, k -> new HashMap<>()); - if (age >= WoodType.IND_ROUNDWOOD.getMinRotation()) { - areaMap.merge(WoodType.IND_ROUNDWOOD, area, Double::sum); - stockMap.merge(WoodType.IND_ROUNDWOOD, yield * area, Double::sum); - } else { - areaMap.merge(WoodType.FUELWOOD, area, Double::sum); - stockMap.merge(WoodType.FUELWOOD, yield * area, Double::sum); - } - } - } - } - - for (Map.Entry<LccKey, Map<WoodType, Double>> entry : stockByType.entrySet()) { - LccKey key = entry.getKey(); - Map<WoodType, Double> stockMap = entry.getValue(); - for (WoodType w : stockMap.keySet()) { - double area = areaByType.get(key).get(w); - double yield = area > 0 ? stockMap.get(w) / area : 0; - yieldByType.computeIfAbsent(key, k -> new HashMap<>()).put(w, yield); - } - } - - return yieldByType; - } - - public Map<LandCoverType, Map<Integer, Double>> getYieldRotaMap() { - return yieldRota; - } - public Map<Integer, Double> getYieldResponseMap() { return yieldResponse.getYieldMap(); } diff --git a/src/ac/ed/lurg/forestry/WoodYieldReader.java b/src/ac/ed/lurg/forestry/WoodYieldReader.java index 53dee70a..f6d7c7cf 100644 --- a/src/ac/ed/lurg/forestry/WoodYieldReader.java +++ b/src/ac/ed/lurg/forestry/WoodYieldReader.java @@ -23,7 +23,7 @@ import ac.sac.raster.RasterKey; import ac.sac.raster.RasterSet; public class WoodYieldReader { - private static final double CONVERSION_FACTOR = 10.0 * ModelConfig.WOOD_YIELD_CALIB_FACTOR; // convert kgC/m2 to tC/ha with calib factor + private static final double CONVERSION_FACTOR = 10.0 * ModelConfig.WOOD_YIELD_CALIB_FACTOR / ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; // convert kgC/m2 to m3/ha with calib factor private final RasterHeaderDetails rasterProj; private int modelYear; private Map<Integer, RasterKey> keyMap; @@ -41,32 +41,8 @@ public class WoodYieldReader { long startTime = System.currentTimeMillis(); // Wood yield from rotation readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_ROTA_FORS_FILENAME, - List.of(new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST))); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_ROTA_NTRL_FILENAME, - List.of(new LccKey(LandCoverType.NATURAL, LandCoverType.NATURAL))); - - // Wood yield from LCC - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_PAST_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), List.of(LandCoverType.PASTURE))); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_CROP_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), List.of(LandCoverType.CROPLAND))); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_NTRL_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), List.of(LandCoverType.NATURAL))); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_FORS_FILENAME, - getLccMappings(LandCoverType.getManagedForestTypes(), LandCoverType.getManagedForestTypes())); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_NTRL_TO_PAST_FILENAME, - getLccMappings(List.of(LandCoverType.NATURAL), List.of(LandCoverType.PASTURE))); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_NTRL_TO_CROP_FILENAME, - getLccMappings(List.of(LandCoverType.NATURAL), List.of(LandCoverType.CROPLAND))); - - readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_NTRL_TO_FORS_FILENAME, - getLccMappings(List.of(LandCoverType.NATURAL), LandCoverType.getManagedForestTypes())); + Arrays.asList(new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST))); + LogWriter.println("Reading forestry data took " + (System.currentTimeMillis() - startTime) + " ms"); return woodYieldData; } @@ -113,17 +89,7 @@ public class WoodYieldReader { for (RasterKey key : map.keySet()) { if (landUseRaster.containsKey(key)) { WoodYieldItem wyItem = woodYieldData.get(key.getCol(), key.getRow()); - LandUseItem luItem = landUseRaster.get(key); - for (LccKey lccKey : lccMappings) { - if (lccKey.getFromLc().equals(lccKey.getToLc())) { - wyItem.setYieldRota(lccKey.getFromLc(), map.get(key), luItem.getLandCoverTiles(), maxAge); - if (lccKey.getFromLc().equals(LandCoverType.TIMBER_FOREST)) { - wyItem.setYieldResponse(map.get(key), modelYear); - } - } else { - wyItem.setYieldLuc(lccKey, map.get(key), luItem.getLandCoverTiles(), maxAge); - } - } + wyItem.setYieldResponse(map.get(key), modelYear); } } LogWriter.println("Reading " + fileName + " took " + (System.currentTimeMillis() - startTime) + " ms"); diff --git a/src/ac/ed/lurg/landuse/LandCoverItem.java b/src/ac/ed/lurg/landuse/LandCoverItem.java index 532bc7aa..861d8d93 100644 --- a/src/ac/ed/lurg/landuse/LandCoverItem.java +++ b/src/ac/ed/lurg/landuse/LandCoverItem.java @@ -6,7 +6,6 @@ import java.util.HashMap; import java.util.List; import java.util.Map; -import ac.ed.lurg.ModelConfig; import ac.ed.lurg.types.CropType; import ac.ed.lurg.types.LandCoverType; import ac.ed.lurg.types.LandProtectionType; @@ -18,7 +17,6 @@ import ac.sac.raster.RasterItem; public class LandCoverItem implements RasterItem { private Map<LandCoverType, Double> landCoverFract = new HashMap<LandCoverType, Double>(); // Land cover fraction - private Map<LandCoverType, Map<Integer, Double>> landUseAgeDist = new HashMap<LandCoverType, Map<Integer, Double>>(); private LandCoverTile landCoverAreas = new LandCoverTile(); private Map<CropType, Double> cropFractions = new HashMap<>(); private double totalArea; @@ -62,15 +60,6 @@ public class LandCoverItem implements RasterItem { public void setMaxCropFraction(double maxCropFraction){ this.unavailableFract=1.0-maxCropFraction; } - - public void setLandUseAgeDist(LandCoverType lcType, double prop, int ageGroup) { - Map<Integer, Double> groupMap = landUseAgeDist.computeIfAbsent(lcType, k -> new HashMap<Integer, Double>()); - groupMap.put(ageGroup, prop); - } - - public Map<LandCoverType, Map<Integer, Double>> getLandUseAgeDist() { - return landUseAgeDist; - } public double getProtectedFraction() { return protectedFraction; @@ -97,7 +86,7 @@ public class LandCoverItem implements RasterItem { public void createLandCoverTiles() { double protectableFraction = 0; - for (LandCoverType lcType : LandCoverType.getProtectibleTypes()) { + for (LandCoverType lcType : LandCoverType.getProtectableTypes()) { protectableFraction += getLandCoverFract(lcType); } // need to scale protectedFraction by protectableFraction and cap at 1.0. Gives fraction of protectable LC types to protect @@ -110,7 +99,7 @@ public class LandCoverItem implements RasterItem { // static land covers if (!lcType.isConvertible()) { - landCoverAreas.addArea(lcType, LandProtectionType.UNAVAILABLE, 0, totalArea); // Assumes land cover age is 0 + landCoverAreas.addArea(lcType, LandProtectionType.UNAVAILABLE, totalArea); continue; } // convertible land covers @@ -124,49 +113,8 @@ public class LandCoverItem implements RasterItem { unprotectedArea = totalArea; } - Map<Integer, Double> ageMap = getLandUseAgeDist().get(lcType); - if (ageMap == null) { - ageMap = new HashMap<Integer, Double>(); - ageMap.put(15, 1.0); - } - double totalProp = ageMap.values().stream().reduce(0.0, Double::sum); - if (totalProp == 0.0) { // If no distribution given then assume land is in oldest age group - int maxAgeGroup = ageMap.keySet().stream().max(Integer::compare).get(); - ageMap.put(maxAgeGroup, 1.0); - } - if (totalProp > 1.0) { - // normalise - for (Map.Entry<Integer, Double> entry : ageMap.entrySet()) { - ageMap.put(entry.getKey(), entry.getValue() / totalProp); - } - } - - for (Map.Entry<Integer, Double> ageEntry : ageMap.entrySet()) { - double prop = ageEntry.getValue(); - if (prop == 0.0) - continue; - int ageGroup = ageEntry.getKey(); - int minAge = (int) (ageGroup - 1) * ModelConfig.LAND_COVER_INIT_AGE_GROUP_SIZE + 1; - int maxAge = minAge + ModelConfig.LAND_COVER_INIT_AGE_GROUP_SIZE - 1; - double groupProp = 0.2; - // Assuming even distribution within each age group -// for (int a = minAge; a <= maxAge; a+=2) { -// double unprotectedA = unprotectedArea * prop * groupProp; -// if (unprotectedA > 0) -// landCoverAreas.setArea(lcType ,LandProtectionType.CONVERTIBLE, a, unprotectedA); -// double protectedA = protectedArea * prop * groupProp; -// if (protectedA > 0) -// landCoverAreas.setArea(lcType, LandProtectionType.PROTECTED, a, protectedA); -// } - // double unprotectedA = unprotectedArea * prop * groupProp; - int a = ageGroup * ModelConfig.LAND_COVER_INIT_AGE_GROUP_SIZE - 5; - double unprotectedA = unprotectedArea * prop; - if (unprotectedA > 0) - landCoverAreas.setArea(lcType ,LandProtectionType.CONVERTIBLE, a, unprotectedA); - double protectedA = protectedArea * prop; - if (protectedA > 0) - landCoverAreas.setArea(lcType, LandProtectionType.PROTECTED, a, protectedA); - } + landCoverAreas.setArea(lcType, LandProtectionType.CONVERTIBLE, unprotectedArea); + landCoverAreas.setArea(lcType, LandProtectionType.PROTECTED, protectedArea); } // Check that areas have been allocated correctly @@ -181,7 +129,7 @@ public class LandCoverItem implements RasterItem { private void setUnavailableArea() { double alreadyUnavail = 0; for (LandCoverType lcType : LandCoverType.values()) { - alreadyUnavail += landCoverAreas.getTotalArea(lcType, LandProtectionType.UNAVAILABLE); + alreadyUnavail += landCoverAreas.getArea(lcType, LandProtectionType.UNAVAILABLE); } double unavailableArea = unavailableFract * totalArea; double shortfall = Math.max(0, unavailableArea - alreadyUnavail); @@ -189,15 +137,15 @@ public class LandCoverItem implements RasterItem { return; for (LandCoverType lcType : unavailPriorityList) { - double cArea = landCoverAreas.getTotalArea(lcType, LandProtectionType.CONVERTIBLE); - double pArea = landCoverAreas.getTotalArea(lcType, LandProtectionType.PROTECTED); + 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, 0, additionalUnavail); + landCoverAreas.addArea(LandCoverType.NATURAL, LandProtectionType.UNAVAILABLE, additionalUnavail); landCoverAreas.removeArea(lcType, LandProtectionType.CONVERTIBLE, additionalUnavail * cProp); landCoverAreas.removeArea(lcType, LandProtectionType.PROTECTED, additionalUnavail * (1 - cProp)); shortfall -= additionalUnavail; @@ -207,7 +155,7 @@ public class LandCoverItem implements RasterItem { break; } - if (shortfall > 1e-10) { + if (shortfall >= 1e-10) { LogWriter.printlnWarning("LandUseItem.setUnavailableArea insufficient area to make unavailable"); } } diff --git a/src/ac/ed/lurg/landuse/LandCoverTile.java b/src/ac/ed/lurg/landuse/LandCoverTile.java index 3fb40566..f19c07a9 100644 --- a/src/ac/ed/lurg/landuse/LandCoverTile.java +++ b/src/ac/ed/lurg/landuse/LandCoverTile.java @@ -19,17 +19,25 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land public LandCoverTile() {} - public double getArea(LandCoverType lcType, LandProtectionType lpType, int age) { - LandKey key = new LandKey(lcType, lpType, age); + public double getArea(LandCoverType lcType, LandProtectionType lpType) { + LandKey key = new LandKey(lcType, lpType); return getArea(key); } public double getArea(LandKey key) { return areas.getOrDefault(key, 0.0); } + + public double getTotalArea(LandCoverType lcType) { + double total = 0; + for (LandProtectionType lpType : LandProtectionType.values()) { + total += getArea(lcType, lpType); + } + return total; + } - public void setArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) { - LandKey key = new LandKey(lcType, lpType, age); + public void setArea(LandCoverType lcType, LandProtectionType lpType, double a) { + LandKey key = new LandKey(lcType, lpType); setArea(key, a); } @@ -37,8 +45,8 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land areas.put(key, a); } - public void addArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) { - LandKey key = new LandKey(lcType, lpType, age); + public void addArea(LandCoverType lcType, LandProtectionType lpType, double a) { + LandKey key = new LandKey(lcType, lpType); addArea(key, a); } @@ -46,94 +54,33 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land areas.merge(key, a, Double::sum); } - // This is used during calibration to preserve Land cover age distribution - public void addArea(LandCoverType lcType, LandProtectionType lpType, double a) { - double totalArea = getTotalArea(lcType, lpType); - if (totalArea == 0) { - addArea(lcType, lpType, 0, a); // assuming all age=0 if no previous land cover - - } else { - double factor = (a + totalArea) / totalArea; // proportional allocation - for (LandKey key : areas.keySet()) { - if (key.getLpType().equals(lpType) && key.getLcType().equals(lcType)) { - double prevArea = getArea(key); - setArea(key, prevArea * factor); - } - } - } - } - - public double getTotalArea(LandCoverType lcType, LandProtectionType lpType) { - double a = 0; - for (LandKey key : areas.keySet()) { - if (key.getLpType().equals(lpType) && key.getLcType().equals(lcType)) { - a += areas.get(key); - } - } - return a; - } - - public void removeArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) { + public void removeArea(LandCoverType lcType, LandProtectionType lpType, double a) { if (a == 0) return; - LandKey key = new LandKey(lcType, lpType, age); + LandKey key = new LandKey(lcType, lpType); double prevArea = getArea(key); - if (a > prevArea) { + if (a - prevArea > 1e-10) { LogWriter.printlnError("LandCoverTile: attempting to remove too much area."); - a = Math.min(a, prevArea); } - setArea(key, prevArea - a); - } + a = Math.min(a, prevArea); - public void removeArea(LandCoverType lcType, LandProtectionType lpType, double a) { // Subtracts area proportionally - if (a == 0) - return; - double totalArea = getTotalArea(lcType, lpType); - double areaToRemove = Math.min(a, totalArea); - if (Math.abs(a - areaToRemove) > 1e-10) { - LogWriter.printlnError("LandCoverTile: attempting to remove too much area."); - } - double factor = (totalArea - areaToRemove) / totalArea; - for (LandKey key : areas.keySet()) { - if (key.getLcType().equals(lcType) && key.getLpType().equals(lpType)) { - double prevArea = getArea(key); - setArea(key, prevArea * factor); - } - } - } - - public double getTotalArea(LandCoverType lcType) { - double area = 0; - for (LandProtectionType lpType : LandProtectionType.values()) { - area += getTotalArea(lcType, lpType); - } - return area; + setArea(key, prevArea - a); } public void removeAllArea(LandCoverType lcType, LandProtectionType lpType) { Map<LandKey, Double> newAreas = new HashMap<LandKey, Double>(); for (LandKey key : areas.keySet()) { if (!(key.getLpType().equals(lpType) && key.getLcType().equals(lcType))) { - newAreas.put(key, getArea(key));; + newAreas.put(key, getArea(key)); } } areas.clear(); areas.putAll(newAreas); } - - public void incrementAge() { - Map<LandKey, Double> newAreas = new HashMap<LandKey, Double>(); - for (LandKey key : areas.keySet()) { - LandKey newKey = new LandKey(key.getLcType(), key.getLpType(), key.getAge() + ModelConfig.TIMESTEP_SIZE); - newAreas.put(newKey, getArea(key)); - } - areas.clear(); - areas.putAll(newAreas); - } // Moves area between land protection types public void moveArea(LandCoverType lcType, LandProtectionType fromPt, LandProtectionType toPt, double area) { - double fromA = getTotalArea(lcType, fromPt); + double fromA = getArea(lcType, fromPt); if (area > fromA) { LogWriter.printlnWarning("LandCoverTile.moveArea cannot protect all required area"); area = Math.min(fromA, area); @@ -142,20 +89,12 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land for (LandKey key : areas.keySet()) { if (key.getLpType().equals(fromPt) & key.getLcType().equals(lcType)) { double areaToMove = getArea(key) * fractToMove; - addArea(lcType, toPt, key.getAge(), areaToMove); - removeArea(lcType, fromPt, key.getAge(), areaToMove); + addArea(lcType, toPt, areaToMove); + removeArea(lcType, fromPt, areaToMove); } } } - - public Set<Integer> getAgeKeys() { // TODO can be more efficient - Set<Integer> ageKeys = new HashSet<Integer>(); - for (LandKey key : areas.keySet()) { - ageKeys.add(key.getAge()); - } - return ageKeys; - } - + // remove 0 areas public void cleanUp() { Map<LandKey, Double> newAreas = new HashMap<LandKey, Double>(); diff --git a/src/ac/ed/lurg/landuse/LandKey.java b/src/ac/ed/lurg/landuse/LandKey.java index 287b3c54..2e48e997 100644 --- a/src/ac/ed/lurg/landuse/LandKey.java +++ b/src/ac/ed/lurg/landuse/LandKey.java @@ -11,13 +11,11 @@ public class LandKey implements Serializable { private LandCoverType lcType; private LandProtectionType lpType; - private int age; // years since conversion - public LandKey(LandCoverType lcType, LandProtectionType lpType, int age) { + public LandKey(LandCoverType lcType, LandProtectionType lpType) { super(); this.lcType = lcType; this.lpType = lpType; - this.age = age; } public LandCoverType getLcType() { @@ -27,20 +25,16 @@ public class LandKey implements Serializable { return lpType; } - public int getAge() { - return age; - } - @Override public boolean equals(Object o) { if (this == o) return true; if (o == null || getClass() != o.getClass()) return false; LandKey landKey = (LandKey) o; - return age == landKey.age && lcType == landKey.lcType && lpType == landKey.lpType; + return lcType == landKey.lcType && lpType == landKey.lpType; } @Override public int hashCode() { - return Objects.hash(lcType, lpType, age); + return Objects.hash(lcType, lpType); } } diff --git a/src/ac/ed/lurg/landuse/LandTileReader.java b/src/ac/ed/lurg/landuse/LandTileReader.java deleted file mode 100644 index 867b1f57..00000000 --- a/src/ac/ed/lurg/landuse/LandTileReader.java +++ /dev/null @@ -1,26 +0,0 @@ -package ac.ed.lurg.landuse; - -import java.util.Map; - -import ac.ed.lurg.types.LandCoverType; -import ac.sac.raster.AbstractTabularRasterReader; -import ac.sac.raster.RasterKey; -import ac.sac.raster.RasterSet; - -public class LandTileReader extends AbstractTabularRasterReader<LandCoverItem> { - private static final int MIN_COLS = 7; - - public LandTileReader(RasterSet<LandCoverItem> data) { - super("[ |\t]+", MIN_COLS, data); - } - - @Override - protected void setData(RasterKey key, LandCoverItem item, Map<String, Double> rowValues) { - int ageGroup = (int) getValueForCol(rowValues, "AgeGroup"); - item.setLandUseAgeDist(LandCoverType.CROPLAND, getValueForCol(rowValues, "CROPLAND"), ageGroup); - item.setLandUseAgeDist(LandCoverType.PASTURE, getValueForCol(rowValues, "PASTURE"), ageGroup); - item.setLandUseAgeDist(LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "MANAGED_FOREST"),ageGroup); - item.setLandUseAgeDist(LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "MANAGED_FOREST"),ageGroup); - item.setLandUseAgeDist(LandCoverType.NATURAL, getValueForCol(rowValues, "NATURAL"), ageGroup); - } -} diff --git a/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java b/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java deleted file mode 100644 index b28fe54c..00000000 --- a/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java +++ /dev/null @@ -1,180 +0,0 @@ -package ac.ed.lurg.landuse; - -import java.io.BufferedInputStream; -import java.io.BufferedOutputStream; -import java.io.DataInputStream; -import java.io.DataOutputStream; -import java.io.FileInputStream; -import java.io.FileOutputStream; -import java.io.IOException; -import java.util.HashMap; -import java.util.Map; - -import ac.ed.lurg.ModelConfig; -import ac.ed.lurg.types.LandCoverType; -import ac.ed.lurg.types.LandProtectionType; -import ac.ed.lurg.utils.LogWriter; -import ac.sac.raster.RasterKey; -import ac.sac.raster.RasterSet; - -public class LandUseBinarySerializer { - private static final int NUM_PROTECTION_CLASSES = 3; // convertible, protected, unavailable -// private static final int NUM_AGE_CLASSES = LandCoverTile.getMaxAgeBin() + 1; - private static final int NUM_AGE_CLASSES = 250 + 1; - private static final int CONVERTIBLE_IDX = 0; - private static final int PROTECTED_IDX = 1; - private static final int UNVAVAILABLE_IDX = 2; - - public LandUseBinarySerializer() {} - - public void serializeLandUse(RasterSet<LandUseItem> luRaster, String filePath) { - long startTime = System.currentTimeMillis(); - - FileOutputStream outputStream = null; - DataOutputStream writer = null; - - try { - outputStream = new FileOutputStream(filePath); - writer = new DataOutputStream(new BufferedOutputStream(outputStream)); - - // TODO file header and checks: num lc types, order, num protection classes, vector size - for (Map.Entry<RasterKey, LandUseItem> luEntry : luRaster.entrySet()) { - RasterKey key = luEntry.getKey(); - LandUseItem luItem = luEntry.getValue(); - - // Write coordinates - writer.writeDouble(key.getCol()); - writer.writeDouble(key.getRow()); - - LandCoverTile tiles = luItem.getLandCoverTiles(); - - // Convertible - for (LandCoverType lcType : LandCoverType.values()) { - for (int a = 0; a < NUM_AGE_CLASSES; a++) { - writer.writeDouble(tiles.getArea(lcType, LandProtectionType.CONVERTIBLE, a)); - } - } - - // Protected - for (LandCoverType lcType : LandCoverType.values()) { - for (int a = 0; a < NUM_AGE_CLASSES; a++) { - writer.writeDouble(tiles.getArea(lcType, LandProtectionType.PROTECTED, a)); - } - } - - // Unavailable - for (LandCoverType lcType : LandCoverType.values()) {; - writer.writeDouble(tiles.getTotalArea(lcType, LandProtectionType.UNAVAILABLE)); - } - } - writer.close(); - - - } - catch (Exception e) { - LogWriter.printlnError("Problem reading data file " + filePath); - LogWriter.print(e); - throw new RuntimeException(e); - } - finally { - if (outputStream != null) { - try { - outputStream.close(); - } catch (IOException e) { - LogWriter.print(e); - } - } - if (writer != null) { - try { - writer.close(); - } catch (IOException e) { - LogWriter.print(e); - } - } - } - LogWriter.println("Serialized land cover " + ModelConfig.SERIALIZED_LAND_COVER_FILE + ", took " + (System.currentTimeMillis() - startTime) + " ms"); - } - - public void deserializeLandUse(RasterSet<LandUseItem> luRaster, String filePath) { - long startTime = System.currentTimeMillis(); - - FileInputStream inputStream = null; - DataInputStream reader = null; - - try { - inputStream = new FileInputStream(filePath); - BufferedInputStream buffInputStream = new BufferedInputStream(inputStream); - buffInputStream.mark(8); - reader = new DataInputStream(buffInputStream); - - while(reader.available() > 0) { - // Coordinates loop - - int col = (int) reader.readDouble(); - int row = (int) reader.readDouble(); - - LandCoverTile landCoverAreas = new LandCoverTile(); - - LandUseItem luItem = luRaster.get(new RasterKey(col, row)); - - // Protection class loop - for (int protIdx = 0; protIdx < NUM_PROTECTION_CLASSES; protIdx++) { - - switch(protIdx) { - case CONVERTIBLE_IDX: - for (LandCoverType lcType : LandCoverType.values()) { - for (int age = 0; age < NUM_AGE_CLASSES; age++) { - double area = reader.readDouble(); - if (area > 0) - landCoverAreas.setArea(lcType, LandProtectionType.CONVERTIBLE, age, area); - } - } - break; - case PROTECTED_IDX: - for (LandCoverType lcType : LandCoverType.values()) { - for (int age = 0; age < NUM_AGE_CLASSES; age++) { - double area = reader.readDouble(); - if (area > 0) - landCoverAreas.setArea(lcType, LandProtectionType.PROTECTED, age, area); - } - } - break; - case UNVAVAILABLE_IDX: - for (LandCoverType lcType : LandCoverType.values()) { - double area = reader.readDouble(); - if (area > 0) - landCoverAreas.setArea(lcType, LandProtectionType.UNAVAILABLE, 0, area); // TODO Do we need to keep track of age? - } - } - } - - luItem.overwriteLandCoverAreas(landCoverAreas); - } - - reader.close(); - } - catch (Exception e) { - LogWriter.printlnError("Problem reading data file " + ModelConfig.SERIALIZED_LAND_COVER_FILE); - LogWriter.print(e); - throw new RuntimeException(e); - } - finally { - if (inputStream != null) { - try { - inputStream.close(); - } catch (IOException e) { - LogWriter.print(e); - } - } - if (reader != null) { - try { - reader.close(); - } catch (IOException e) { - LogWriter.print(e); - } - } - } - LogWriter.println("Deserialized land cover " + ModelConfig.SERIALIZED_LAND_COVER_FILE + ", took " + (System.currentTimeMillis() - startTime) + " ms"); - } - -} diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java index d343e173..fbe118bf 100644 --- a/src/ac/ed/lurg/landuse/LandUseItem.java +++ b/src/ac/ed/lurg/landuse/LandUseItem.java @@ -5,7 +5,6 @@ import java.util.Collection; import java.util.HashMap; import java.util.Map; import ac.ed.lurg.ModelConfig; -import ac.ed.lurg.forestry.WoodHarvestItem; import ac.ed.lurg.types.*; import ac.ed.lurg.utils.Interpolator; import ac.sac.raster.InterpolatingRasterItem; @@ -16,10 +15,9 @@ 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 LandCoverTile landCoverAreas = new LandCoverTile(); - private Map<WoodType, Double> forestFractions = new HashMap<>(); private Map<LccKey, Double> landCoverChange = new HashMap<LccKey, Double>(); // land cover change for previous timestep - private WoodHarvestItem woodHarvests; // forest rotations private double protectedFraction; // protected fraction of total area + private double forestRotationIntensity; public LandUseItem() {} @@ -33,6 +31,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial setProtectedFraction(landCover.getProtectedFraction()); this.landCoverAreas = (landCover.getLandCoverTiles()); } + forestRotationIntensity = 0.02; } public LandUseItem(LandUseItem luItemToCopy) { @@ -41,6 +40,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial cropFractions.putAll(luItemToCopy.cropFractions); landCoverAreas = (luItemToCopy.landCoverAreas); protectedFraction = (luItemToCopy.protectedFraction); + forestRotationIntensity = (luItemToCopy.forestRotationIntensity); } public void overwriteLandCoverAreas(LandCoverTile landCoverAreas) { @@ -68,16 +68,15 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial } public double getLandCoverArea(LandCoverType lcType, LandProtectionType lpType) { - return landCoverAreas.getTotalArea(lcType, lpType); - } - - public double getLandCoverArea(LandCoverType lcType, LandProtectionType lpType, int age) { - return landCoverAreas.getArea(lcType, lpType, age); + return landCoverAreas.getArea(lcType, lpType); } public double getLandCoverArea(LandCoverType lcType) { - Double d = landCoverAreas.getTotalArea(lcType); - return d == null ? 0.0 : d; + double d = 0; + for (LandProtectionType lpType : LandProtectionType.values()) { + d += landCoverAreas.getArea(lcType, lpType); + } + return d; } public double getTotalLandArea() { @@ -123,33 +122,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial throw new RuntimeException("Move unprotected area shortfall. Should not happen due to GAMS constraint."); } - if (ModelConfig.IS_CALIBRATION_RUN) { - // Preserve age distribution of source land cover. - // This is mainly so that forest <> non-forest transitions don't lose trees. - if (fromType.equals(LandCoverType.NATURAL) && toType.equals(LandCoverType.TIMBER_FOREST)) { - double moveFraction = areaMoved / prevFromArea; - for (int age : landCoverAreas.getAgeKeys()) { - double fromArea = landCoverAreas.getArea(fromType, LandProtectionType.CONVERTIBLE, age); - if (fromArea == 0) - continue; - double partialArea = fromArea * moveFraction; - landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, age, partialArea); - landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, age, partialArea); - } - } else { - - double areaMovedCalib = areaMoved; - if (fromType.isNatural()) { - areaMovedCalib *= 0.95; // keep some area during calibration to preserve age distribution - } - landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMovedCalib); - landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, areaMovedCalib); - } - - } else { - landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMoved); - landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, 0, areaMoved); // new area starts at age=0 - } + landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMoved); + landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, areaMoved); // new area starts at age=0 LccKey lccKey = new LccKey(fromType, toType); landCoverChange.merge(lccKey, areaMoved, (a, b) -> a + b); // add to previous, set 0 otherwise @@ -157,27 +131,12 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial return area - areaMoved; } - public void doForestRotation(int age, double area) { - if (ModelConfig.IS_CALIBRATION_RUN) - return; - landCoverAreas.removeArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age, area); - landCoverAreas.addArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, 0, area); - } - - public void setForestFraction(WoodType woodType, double fraction) { - forestFractions.put(woodType, fraction); - } - - public double getForestFraction(WoodType woodType) { - return forestFractions.get(woodType); - } - - public void setWoodHarvests(WoodHarvestItem woodHarvests) { - this.woodHarvests = woodHarvests; + public void setForestRotationIntensity(double i) { + forestRotationIntensity = i; } - public WoodHarvestItem getWoodHarvests() { - return woodHarvests; + public double getForestRotationIntensity() { + return forestRotationIntensity; } public void resetLccAreas() { @@ -195,7 +154,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial double landCover = (d < 0.0) ? 0.0 : d; landCoverAreas.removeAllArea(c, p); - landCoverAreas.addArea(c, p, 0, landCover); + landCoverAreas.addArea(c, p, landCover); } public void updateProtectedFraction(int year) { @@ -215,7 +174,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial double previousProtectedTotal = 0; for (LandCoverType lc : LandCoverType.values()) { - previousProtectedTotal += landCoverAreas.getTotalArea(lc, LandProtectionType.PROTECTED); + previousProtectedTotal += landCoverAreas.getArea(lc, LandProtectionType.PROTECTED); } double newProtectedTotal = protectedFraction * getTotalLandArea(); double factor = newProtectedTotal / previousProtectedTotal; @@ -395,10 +354,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial return changes; } - public void incrementTilesAge() { - landCoverAreas.incrementAge(); - } - private boolean isZeroOrNull(Double d) { return d == null || d == 0; } @@ -408,8 +363,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial cropFractions = new HashMap<CropType, Double>(); landCoverAreas = new LandCoverTile(); - Double fromCropCover = fromItem.landCoverAreas.getTotalArea(LandCoverType.CROPLAND); - Double toCropCover = toItem.landCoverAreas.getTotalArea(LandCoverType.CROPLAND); + Double fromCropCover = fromItem.getLandCoverArea(LandCoverType.CROPLAND); + Double toCropCover = toItem.getLandCoverArea(LandCoverType.CROPLAND); if (!isZeroOrNull(fromCropCover) && isZeroOrNull(toCropCover)) { // if start with crop but end with none, take starting crop fractions cropFractions.putAll(fromItem.cropFractions); @@ -438,15 +393,13 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial intensityMap.put(crop, interpolateIntensity); } } - + for (LandCoverType lcType : LandCoverType.values()) { for (LandProtectionType lpType : LandProtectionType.values()) { - for (int age : landCoverAreas.getAgeKeys()) { - Double from = fromItem.landCoverAreas.getArea(lcType, lpType, age); - Double to = toItem.landCoverAreas.getArea(lcType, lpType, age); - Double d = Interpolator.interpolate(from, to, factor); - landCoverAreas.setArea(lcType, lpType, age, d); - } + Double from = fromItem.landCoverAreas.getArea(lcType, lpType); + Double to = toItem.landCoverAreas.getArea(lcType, lpType); + Double d = Interpolator.interpolate(from, to, factor); + landCoverAreas.setArea(lcType, lpType, d); } } } @@ -476,8 +429,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial Map<LandCoverType, Double> convertibleAreas = new HashMap<LandCoverType, Double>(); Map<LandCoverType, Double> protectedAreas = new HashMap<LandCoverType, Double>(); for (LandCoverType lcType : LandCoverType.values()) { - convertibleAreas.put(lcType, landCoverAreas.getTotalArea(lcType, LandProtectionType.CONVERTIBLE)); - protectedAreas.put(lcType, landCoverAreas.getTotalArea(lcType, LandProtectionType.PROTECTED)); + convertibleAreas.put(lcType, landCoverAreas.getArea(lcType, LandProtectionType.CONVERTIBLE)); + protectedAreas.put(lcType, landCoverAreas.getArea(lcType, LandProtectionType.PROTECTED)); } return "LandUseItem: [landCoverAreas=" + convertibleAreas + ", protectedArea=" + protectedAreas + "]"; } diff --git a/src/ac/ed/lurg/landuse/WoodUsageData.java b/src/ac/ed/lurg/landuse/WoodUsageData.java index 8dedd1d9..3771e670 100644 --- a/src/ac/ed/lurg/landuse/WoodUsageData.java +++ b/src/ac/ed/lurg/landuse/WoodUsageData.java @@ -6,28 +6,39 @@ public class WoodUsageData implements Serializable { private static final long serialVersionUID = -3329080279189782763L; - private double rotaHarvest; + private double expectedProduction; private double netImport; - private double lucHarvest = 0; + private final double initialDemand; + private final double baseImportDemandRatio; - public WoodUsageData(double harvest, double netImport) { - this.rotaHarvest = harvest; + public WoodUsageData(double harvest, double netImport, double initialDemand) { + this.expectedProduction = harvest; this.netImport = netImport; + this.initialDemand = initialDemand; + this.baseImportDemandRatio = initialDemand > 0 ? netImport / initialDemand : 0; } - public double getHarvest() { - return rotaHarvest; + public double getProduction() { + return expectedProduction; + } + + public void setProduction(double production) { + this.expectedProduction = production; } public double getNetImport() { return netImport; } - public double getLucHarvest() { - return lucHarvest; + public void setNetImport(double netImport) { + this.netImport = netImport; + } + + public double getInitialDemand() { + return initialDemand; } - public void setLucHarvest(double harvest) { - lucHarvest = harvest; + public double getBaseImportDemandRatio() { + return baseImportDemandRatio; } } diff --git a/src/ac/ed/lurg/landuse/WoodUsageReader.java b/src/ac/ed/lurg/landuse/WoodUsageReader.java index a9d0862d..0eff756a 100644 --- a/src/ac/ed/lurg/landuse/WoodUsageReader.java +++ b/src/ac/ed/lurg/landuse/WoodUsageReader.java @@ -1,8 +1,5 @@ package ac.ed.lurg.landuse; -import java.io.BufferedReader; -import java.io.FileReader; -import java.io.IOException; import java.util.HashMap; import java.util.List; import java.util.Map; @@ -13,7 +10,6 @@ import ac.ed.lurg.ModelConfig; import ac.ed.lurg.country.CompositeCountry; import ac.ed.lurg.types.WoodType; import ac.ed.lurg.utils.LazyHashMap; -import ac.ed.lurg.utils.LogWriter; import ac.ed.lurg.utils.StringTabularReader; public class WoodUsageReader { @@ -35,7 +31,7 @@ public class WoodUsageReader { for (CompositeCountry cc : CountryManager.getInstance().getAllCompositeCountries()) { Map<WoodType, WoodUsageData> countryData = usageMap.lazyGet(cc); for (WoodType w : WoodType.values()) { - WoodUsageData newData = new WoodUsageData(0, 0); + WoodUsageData newData = new WoodUsageData(0, 0, 0); countryData.put(w, newData); } } @@ -50,8 +46,9 @@ public class WoodUsageReader { for (Map<String, String> row : rowList) { String countryCode = row.get("Iso3"); String woodTypeName = row.get("Type"); - double production = Double.parseDouble(row.get("Production")) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; - double netImport = Double.parseDouble(row.get("NetImports")) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; + double production = Double.parseDouble(row.get("Production")); + double netImport = Double.parseDouble(row.get("NetImports")); + double demand = Double.parseDouble(row.get("Demand")); SingleCountry country = CountryManager.getInstance().getForCode(countryCode); WoodType woodType = WoodType.getForName(woodTypeName); CompositeCountry cc = CountryManager.getInstance().getForSingleCountry(country); @@ -62,10 +59,11 @@ public class WoodUsageReader { // aggregate if multiple countries for cc if (oldData != null) { netImport += oldData.getNetImport(); - production += oldData.getHarvest(); + production += oldData.getProduction(); + demand += oldData.getInitialDemand(); } - WoodUsageData newData = new WoodUsageData(production, netImport); + WoodUsageData newData = new WoodUsageData(production, netImport, demand); countryData.put(woodType, newData); } return usageMap; diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java index 25de106e..9eae0f51 100644 --- a/src/ac/ed/lurg/types/LandCoverType.java +++ b/src/ac/ed/lurg/types/LandCoverType.java @@ -15,7 +15,7 @@ public enum LandCoverType { CROPLAND ("cropland", false, true, false, false, false), PASTURE ("pasture", false, true, false, false, false), BARREN ("barren", false, false, false, false, false), - URBAN("urban", true, false, false, false, false); + URBAN("urban", false, false, false, false, false); private String name; private boolean isProtectable; @@ -84,7 +84,7 @@ public enum LandCoverType { } - public static Collection<LandCoverType> getProtectibleTypes() { + public static Collection<LandCoverType> getProtectableTypes() { Collection<LandCoverType> protectibleTypes = new HashSet<LandCoverType>(); -- GitLab