From 38a0cc11c066f93fad38f2f399d3dc88fde066aa Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Wed, 1 Sep 2021 13:35:31 +0100
Subject: [PATCH] Changes to carbon flux data read-in.
Wood yields and C flux bug fixes.
---
GAMS/IntExtOpt.gms | 3 +-
src/ac/ed/lurg/ModelMain.java | 34 +----
src/ac/ed/lurg/carbon/CarbonFluxItem.java | 100 +++++++-------
src/ac/ed/lurg/carbon/CarbonFluxReader.java | 122 ++++++++++++------
.../country/gams/GamsLocationOptimiser.java | 28 ++--
.../country/gams/GamsRasterOptimiser.java | 46 +++----
src/ac/ed/lurg/forestry/WoodYieldItem.java | 5 +-
src/ac/ed/lurg/forestry/WoodYieldReader.java | 8 +-
8 files changed, 194 insertions(+), 152 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index f1fbb16c..2d7b2a04 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -52,6 +52,7 @@
PARAMETER previousLandCoverArea(land_cover, location) land cover area in Mha;
PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location) carbon flux - MtC-eq per Mha;
+ PARAMETER carbonFluxRateNEE(land_cover, location);
PARAMETER woodYieldLuc(land_cover_before, land_cover_after, location) wood yield - MtC-eq per Mha;
PARAMETER woodYieldRota(land_cover_before, land_cover_after, location);
@@ -89,7 +90,7 @@ $load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImpo
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
$load previousLandCoverArea, carbonFluxRate, woodYieldLuc, carbonPrice, woodPrice
$load conversionCost, woodYieldRota, infrastructureCost, forestRotationPeriod, woodMaxNetImport, woodMinNetImport, managedForestAreaMaxChangeRate, woodDemand
-$load forestEstablishmentCost, naturalAreaValue, woodFromRotation, vegClearingCost, discountRate
+$load forestEstablishmentCost, naturalAreaValue, woodFromRotation, vegClearingCost, discountRate, carbonFluxRateNEE
$gdxin
SCALAR delta "use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf" / 0.00000000001 /;
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 112f2cdd..f2a1adce 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -86,6 +86,7 @@ public class ModelMain {
private RasterSet<LandUseItem> globalLandUseRaster;
private RasterSet<IntegerRasterItem> clusterIdRaster;
WoodYieldReader woodYieldReader;
+ CarbonFluxReader carbonFluxReader;
public static void main(String[] args) {
@@ -122,6 +123,7 @@ public class ModelMain {
internationalMarket = new InternationalMarket();
woodYieldReader = new WoodYieldReader(desiredProjection);
+ carbonFluxReader = new CarbonFluxReader(desiredProjection);
createCountryAgents(compositeCountryManager.getAll());
}
@@ -145,7 +147,6 @@ public class ModelMain {
YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
getUpdateIrrigationData(timestep, yieldSurfaces); // updating currentIrrigationData
- WoodYieldRasterSet woodYieldData = getWoodYieldData(timestep);
// When running half earth we can to alter protected areas data at a point in time
if(ModelConfig.HALFEARTH && ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR == timestep.getYear() && !ModelConfig.IS_CALIBRATION_RUN) {
@@ -161,13 +162,13 @@ public class ModelMain {
double carbonDemand = demandManager.getCarbonDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
double carbonDemandIncrease = (carbonDemand > previousCarbonDemand) ? carbonDemand - previousCarbonDemand: 0;
- CarbonFluxRasterSet currentCarbonFluxData = getCarbonFluxData(timestep);
- WoodYieldRasterSet currentWoodYieldData = getInterpolatedWoodYieldData(timestep);
+ CarbonFluxRasterSet carbonFluxData = getCarbonFluxData(timestep);
+ WoodYieldRasterSet woodYieldData = getWoodYieldData(timestep);
Map<LccKey, Double> conversionCosts;
conversionCosts = (ModelConfig.CONVERSION_COSTS_FROM_FILE) ? new ConversionCostReader().read() : new ConversionCostReader().calcFromConfig();
- countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, woodYieldData,
+ countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, carbonFluxData, woodYieldData,
conversionCosts, carbonDemandIncrease);
internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep);
@@ -691,9 +692,6 @@ public class ModelMain {
RasterSet<LandUseItem> landUseRaster = new RasterSet<LandUseItem>(initialLC.getHeaderDetails());
for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
- if (entry.getKey().getCol() == 177 && entry.getKey().getRow() == 145) {
- int foo = 1;
- }
landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue()));
}
@@ -719,33 +717,13 @@ public class ModelMain {
/** Get carbon flux data */
private CarbonFluxRasterSet getCarbonFluxData(Timestep timestep) {
- CarbonFluxRasterSet carbonFluxData = new CarbonFluxRasterSet(desiredProjection);
- LogWriter.println("Reading carbon flux data");
- new CarbonFluxReader(carbonFluxData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + "cflux_net.out");
- return carbonFluxData;
+ return carbonFluxReader.getCarbonFluxes(globalLandUseRaster, timestep);
}
/** Get wood yield data */
private WoodYieldRasterSet getWoodYieldData(Timestep timestep) {
return woodYieldReader.getWoodYields(globalLandUseRaster, timestep, internationalMarket.getWoodPrice().getExportPrice());
}
-
-
- private WoodYieldRasterSet getInterpolatedWoodYieldData(Timestep timestep) {
- WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
- /*int yearCurrent = timestep.getYear();
- int yearLower = Math.floorDiv(yearCurrent, WoodYieldItem.YEAR_INTERVAL) * WoodYieldItem.YEAR_INTERVAL + WoodYieldItem.YEAR_OFFSET;
- int yearUpper = yearLower + WoodYieldItem.YEAR_INTERVAL;
- int age =
- getWoodYieldData(woodYieldData, yearLower);
- getWoodYieldData(woodYieldData, yearUpper);
- for (WoodYieldItem wyItem : woodYieldData.values()) {
- wyItem.interpolateAllForYear(yearCurrent);
- wyItem.removeDataForYear(yearLower);
- wyItem.removeDataForYear(yearUpper);
- }*/
- return woodYieldData;
- }
/** Ugly in situ update of currentIrrigationData, better if IrrigationRasterSets were handled more immutably */
private void getUpdateIrrigationData(Timestep timestep, YieldRaster yieldSurfaces) {
diff --git a/src/ac/ed/lurg/carbon/CarbonFluxItem.java b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
index 044782d5..e3ccfabd 100644
--- a/src/ac/ed/lurg/carbon/CarbonFluxItem.java
+++ b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
@@ -3,75 +3,85 @@ package ac.ed.lurg.carbon;
import ac.sac.raster.RasterItem;
import java.io.Serializable;
+import java.util.ArrayList;
import java.util.HashMap;
import java.util.Map;
+import ac.ed.lurg.Timestep;
+import ac.ed.lurg.landuse.LandUseTile;
import ac.ed.lurg.landuse.LccKey;
import ac.ed.lurg.types.LandCoverType;
public class CarbonFluxItem implements RasterItem, Serializable {
private static final long serialVersionUID = 440720456140537815L;
- Map<LandCoverType, Map<LandCoverType, Double>> carbonFluxes = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
- Map<LccKey, AgeFunction> conversionCarbonFlux = new HashMap<LccKey, AgeFunction>();
- Map<LandCoverType, AgeFunction> ecosystemCarbonFlux = new HashMap<LandCoverType, AgeFunction>();
+ Map<LccKey, Double> conversionCarbonFlux = new HashMap<LccKey, Double>();
+ Map<LandCoverType, Double> ecosystemCarbonFlux = new HashMap<LandCoverType, Double>();
- public void setConversionCarbonFlux(LandCoverType fromLc, LandCoverType toLc, int year, int age, double cFlux) {
- LccKey key = new LccKey(fromLc, toLc);
- AgeFunction ageFunc = conversionCarbonFlux.get(key);
- if (ageFunc == null) {
- ageFunc = new AgeFunction();
- conversionCarbonFlux.put(key, ageFunc);
- }
- ageFunc.setData(year, age, cFlux);
+ public void setConversionCarbonFlux(LccKey key, double cFlux) {
+ conversionCarbonFlux.put(key, cFlux);
}
- public void setEcosystemCarbonFlux(LandCoverType lcType, int year, int age, double cFlux) {
- AgeFunction ageFunc = ecosystemCarbonFlux.get(lcType);
- if (ageFunc == null) {
- ageFunc = new AgeFunction();
- ecosystemCarbonFlux.put(lcType, ageFunc);
- }
- ageFunc.setData(year, age, cFlux);
- }
-
- public double getConversionCarbonFlux(LccKey key, int year, int age) {
- return conversionCarbonFlux.get(key).getValue(year, age);
+ public void calcConversionCarbonFlux(LccKey key, Double[] cFluxes, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep) {
+
+ ArrayList<LandUseTile> tiles = landUseTiles.get(key.getFromLc());
+
+ if (tiles.size() == 0) {
+ this.conversionCarbonFlux.put(key, 0.0);
+ } else {
+ double totalFlux = 0;
+ double totalArea = 0;
+ for (LandUseTile tile : tiles) {
+ if (tile.isProtected())
+ continue;
+ int age = tile.getAge(timestep);
+ totalFlux += cFluxes[age] * tile.getArea();
+ totalArea += tile.getArea();
+ }
+ double meanFlux = totalFlux / totalArea;
+ this.conversionCarbonFlux.put(key, meanFlux);
+ }
}
- public double getEcosystemCarbonFlux(LandCoverType lcType, int year, int age) {
- return ecosystemCarbonFlux.get(lcType).getValue(year, age);
+ public void setEcosystemCarbonFlux(LandCoverType lcType, double cFlux) {
+ ecosystemCarbonFlux.put(lcType, cFlux);
}
- public void setCarbonFlux(LandCoverType previousLandCover, LandCoverType newLandCover, double carbonFlux) {
+ public void calcEcosystemCarbonFlux(LandCoverType lcType, Double[] cFluxes, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep) {
- if (carbonFluxes.containsKey(previousLandCover)) {
- carbonFluxes.get(previousLandCover).put(newLandCover, carbonFlux);
- } else {
- Map<LandCoverType, Double> cfMap = new HashMap<LandCoverType, Double>();
- cfMap.put(newLandCover, carbonFlux);
- carbonFluxes.put(previousLandCover, cfMap);
+ ArrayList<LandUseTile> tiles = landUseTiles.get(lcType);
+
+ if (tiles.size() == 0) {
+ this.ecosystemCarbonFlux.put(lcType, 0.0);
+ } else {
+ double totalFlux = 0;
+ double totalArea = 0;
+ for (LandUseTile tile : tiles) {
+ if (tile.isProtected())
+ continue;
+ int age = tile.getAge(timestep);
+ totalFlux += cFluxes[age] * tile.getArea();
+ totalArea += tile.getArea();
+ }
+ double meanFlux = totalFlux / totalArea;
+ this.ecosystemCarbonFlux.put(lcType, meanFlux);
}
}
+
+ public double getConversionCarbonFlux(LccKey key) {
+ return conversionCarbonFlux.getOrDefault(key, Double.NaN);
+ }
- public double getCarbonFlux(LandCoverType previousLandCover, LandCoverType newLandCover) {
- return carbonFluxes.get(previousLandCover).get(newLandCover);
+ public double getEcosystemCarbonFlux(LandCoverType lcType) {
+ return ecosystemCarbonFlux.getOrDefault(lcType, Double.NaN);
}
- public boolean checkForKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
- if (carbonFluxes.containsKey(previousLandCover)) {
- if (carbonFluxes.get(previousLandCover).containsKey(newLandCover)) {
- return true;
- } else {
- return false;
- }
- } else {
- return false;
- }
+ public Map<LccKey, Double> getConversionCarbonFluxMap() {
+ return conversionCarbonFlux;
}
- public Map<LandCoverType, Map<LandCoverType, Double>> getCarbonFluxes() {
- return carbonFluxes;
+ public Map<LandCoverType, Double> getEcosystemCarbonFluxMap() {
+ return ecosystemCarbonFlux;
}
}
diff --git a/src/ac/ed/lurg/carbon/CarbonFluxReader.java b/src/ac/ed/lurg/carbon/CarbonFluxReader.java
index 9b712557..79f40500 100644
--- a/src/ac/ed/lurg/carbon/CarbonFluxReader.java
+++ b/src/ac/ed/lurg/carbon/CarbonFluxReader.java
@@ -1,53 +1,103 @@
package ac.ed.lurg.carbon;
+import java.io.File;
import java.util.Map;
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.Timestep;
+import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.LccKey;
import ac.ed.lurg.types.LandCoverType;
import ac.sac.raster.AbstractTabularRasterReader;
+import ac.sac.raster.RasterHeaderDetails;
import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterSet;
-public class CarbonFluxReader extends AbstractTabularRasterReader<CarbonFluxItem> {
-
- private static final int MIN_COLS = 27;
- private static final double CONVERSION_FACTOR = 10; // convert kgC/m2 to tC/ha
+public class CarbonFluxReader {
+ private static final int MIN_COLS = 263;
+ private static final String DELIMITER = "[ |\t]+";
+ private static final double CONVERSION_FACTOR = 10.0; // convert kgC/m2 to tC/ha
+ private RasterHeaderDetails rasterProj;
- public CarbonFluxReader(RasterSet<CarbonFluxItem> carbonFluxes) {
- super("[ |\t]+", MIN_COLS, carbonFluxes);
+ public CarbonFluxReader(RasterHeaderDetails rasterProj) {
+ this.rasterProj = rasterProj;
}
- @Override
- protected void setData(RasterKey key, CarbonFluxItem item, Map<String, Double> rowValues) {
- item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.NATURAL, getValueForCol(rowValues, "crop_to_ntrl") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "crop_to_forT") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "crop_to_forC") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.PASTURE, getValueForCol(rowValues, "crop_to_past") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CROPLAND, LandCoverType.CROPLAND, getValueForCol(rowValues, "crop_to_crop") * CONVERSION_FACTOR);
-
- item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.NATURAL, getValueForCol(rowValues, "past_to_ntrl") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "past_to_forT") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "past_to_forC") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.CROPLAND, getValueForCol(rowValues, "past_to_crop") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.PASTURE, LandCoverType.PASTURE, getValueForCol(rowValues, "past_to_past") * CONVERSION_FACTOR);
+ public CarbonFluxRasterSet getCarbonFluxes(RasterSet<LandUseItem> landUseRaster, Timestep timestep) {
+ CarbonFluxRasterSet cFluxData = new CarbonFluxRasterSet(rasterProj);
+
+ getCarbonFluxesConversion("cflux_conv_ntrl_to_past.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.NATURAL, LandCoverType.PASTURE));
+ getCarbonFluxesConversion("cflux_conv_ntrl_to_forC.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.NATURAL, LandCoverType.CARBON_FOREST));
+ getCarbonFluxesConversion("cflux_conv_ntrl_to_forT.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST));
+ getCarbonFluxesConversion("cflux_conv_ntrl_to_crop.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.NATURAL, LandCoverType.CROPLAND));
+
+ getCarbonFluxesConversion("cflux_conv_past_to_ntrl.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.PASTURE, LandCoverType.NATURAL));
+ getCarbonFluxesConversion("cflux_conv_past_to_forC.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.PASTURE, LandCoverType.CARBON_FOREST));
+ getCarbonFluxesConversion("cflux_conv_past_to_forT.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST));
+ getCarbonFluxesConversion("cflux_conv_past_to_crop.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.PASTURE, LandCoverType.CROPLAND));
+
+ getCarbonFluxesConversion("cflux_conv_forC_to_ntrl.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CARBON_FOREST, LandCoverType.NATURAL));
+ getCarbonFluxesConversion("cflux_conv_forC_to_past.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE));
+ getCarbonFluxesConversion("cflux_conv_forC_to_forT.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST));
+ getCarbonFluxesConversion("cflux_conv_forC_to_crop.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND));
+
+ getCarbonFluxesConversion("cflux_conv_forT_to_ntrl.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.NATURAL));
+ getCarbonFluxesConversion("cflux_conv_forT_to_past.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE));
+ getCarbonFluxesConversion("cflux_conv_forT_to_forC.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.CARBON_FOREST));
+ getCarbonFluxesConversion("cflux_conv_forT_to_forT.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST));
+ getCarbonFluxesConversion("cflux_conv_forT_to_crop.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND));
- item.setCarbonFlux(LandCoverType.NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "ntrl_to_past") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "ntrl_to_forC") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "ntrl_to_crop") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.NATURAL, LandCoverType.NATURAL, getValueForCol(rowValues, "ntrl_to_ntrl") * CONVERSION_FACTOR);
+ getCarbonFluxesConversion("cflux_conv_crop_to_ntrl.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CROPLAND, LandCoverType.NATURAL));
+ getCarbonFluxesConversion("cflux_conv_crop_to_past.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CROPLAND, LandCoverType.PASTURE));
+ getCarbonFluxesConversion("cflux_conv_crop_to_forC.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CROPLAND, LandCoverType.CARBON_FOREST));
+ getCarbonFluxesConversion("cflux_conv_crop_to_forT.out", cFluxData, landUseRaster, timestep, new LccKey(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST));
- item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "forT_to_past") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.NATURAL, getValueForCol(rowValues, "forT_to_ntrl") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "forT_to_forC") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "forT_to_crop") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forT_to_forT") * CONVERSION_FACTOR);
-
- item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "forC_to_past") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.NATURAL, getValueForCol(rowValues, "forC_to_ntrl") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "forC_to_crop") * CONVERSION_FACTOR);
- item.setCarbonFlux(LandCoverType.CARBON_FOREST, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "forC_to_forC") * CONVERSION_FACTOR);
+ getCarbonFluxesEcosystem("cflux_sts_ntrl.out", cFluxData, landUseRaster, timestep, LandCoverType.NATURAL);
+ getCarbonFluxesEcosystem("cflux_sts_past.out", cFluxData, landUseRaster, timestep, LandCoverType.PASTURE);
+ getCarbonFluxesEcosystem("cflux_sts_forC.out", cFluxData, landUseRaster, timestep, LandCoverType.CARBON_FOREST);
+ getCarbonFluxesEcosystem("cflux_sts_forT.out", cFluxData, landUseRaster, timestep, LandCoverType.TIMBER_FOREST);
+ getCarbonFluxesEcosystem("cflux_sts_crop.out", cFluxData, landUseRaster, timestep, LandCoverType.CROPLAND);
-
+ return cFluxData;
+
+ }
+
+ public void getCarbonFluxesConversion(String filename, CarbonFluxRasterSet cFluxData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, LccKey lccKey) {
+
+ AbstractTabularRasterReader<CarbonFluxItem> cFluxReader = new AbstractTabularRasterReader<CarbonFluxItem>(DELIMITER, MIN_COLS, cFluxData) {
+ protected void setData(RasterKey key, CarbonFluxItem item, Map<String, Double> rowValues) {
+
+ Double[] fluxes = getArrayFromRowValues(rowValues);
+ item.calcConversionCarbonFlux(lccKey, fluxes, landUseRaster.get(key).getLandUseTiles(), timestep);
+ }
+ };
+
+ cFluxReader.getRasterDataFromFile(getDataDir(filename, timestep));
+ }
+
+ public void getCarbonFluxesEcosystem(String filename, CarbonFluxRasterSet cFluxData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, LandCoverType lcType) {
+
+ AbstractTabularRasterReader<CarbonFluxItem> cFluxReader = new AbstractTabularRasterReader<CarbonFluxItem>(DELIMITER, MIN_COLS, cFluxData) {
+ protected void setData(RasterKey key, CarbonFluxItem item, Map<String, Double> rowValues) {
+
+ Double[] fluxes = getArrayFromRowValues(rowValues);
+ item.calcEcosystemCarbonFlux(lcType, fluxes, landUseRaster.get(key).getLandUseTiles(), timestep);
+ }
+ };
+
+ cFluxReader.getRasterDataFromFile(getDataDir(filename, timestep));
+ }
+
+ private Double[] getArrayFromRowValues(Map<String, Double> rowValues) {
+ int arrSize = rowValues.size();
+ Double[] arr = new Double[arrSize];
+ for (int i=0; i<arrSize; i++) {
+ arr[i] = rowValues.get(Integer.toString(i)) * CONVERSION_FACTOR;
+ }
+ return arr;
+ }
+
+ private String getDataDir(String filename, Timestep timestep) {
+ return timestep.getWoodAndCarbonYearSubDir(ModelConfig.WOOD_AND_CARBON_DIR) + File.separator + filename;
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 082106bf..0569bf7d 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -361,24 +361,28 @@ public class GamsLocationOptimiser {
// Carbon fluxes
GAMSParameter cFluxRateP = inDB.addParameter("carbonFluxRate", 3);
+ GAMSParameter cFluxRateNeeP = inDB.addParameter("carbonFluxRateNEE", 2);
for (Entry<Integer, ? extends CarbonFluxItem> entry : inputData.getCarbonFluxes().entrySet()) {
Integer locationId = entry.getKey();
String locString = Integer.toString(locationId);
CarbonFluxItem cFlux = entry.getValue();
-
- for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
- for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
- Vector<String> v = new Vector<String>();
- v.add(prevLC.getName());
- v.add(newLC.getName());
- v.add(locString);
- // TODO Need a better way to deal with missing values
- if (cFlux.checkForKeys(prevLC, newLC))
- setGamsParamValue(cFluxRateP.addRecord(v), cFlux.getCarbonFlux(prevLC, newLC), 3);
- }
- }
+ for (Map.Entry<LccKey, Double> cfEntry : cFlux.getConversionCarbonFluxMap().entrySet()) {
+ LccKey key = cfEntry.getKey();
+ Vector<String> v = new Vector<String>();
+ v.add(key.getFromLc().getName());
+ v.add(key.getToLc().getName());
+ v.add(locString);
+ setGamsParamValue(cFluxRateP.addRecord(v), cfEntry.getValue(), 2);
+ }
+
+ for (Map.Entry<LandCoverType, Double> cfEntry : cFlux.getEcosystemCarbonFluxMap().entrySet()) {
+ Vector<String> v = new Vector<String>();
+ v.add(cfEntry.getKey().getName());
+ v.add(locString);
+ setGamsParamValue(cFluxRateNeeP.addRecord(v), cfEntry.getValue(), 3);
+ }
}
// Wood yields
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 4c80bbc5..1f675e9a 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -249,8 +249,6 @@ public class GamsRasterOptimiser {
protected WoodYieldData createValue() { return new WoodYieldData(); }
};
Map<Integer, Double> aggregatedForestRotation = new HashMap<Integer, Double>();
-
- int currentYear = rasterInputData.getTimestep().getYear();
int countFound = 0, countMissing = 0;
@@ -304,32 +302,36 @@ public class GamsRasterOptimiser {
}
// Aggregate carbon fluxes and wood yields
- for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
- for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
- // Carbon flux
- double carbonFlux;
- if (carbonFluxItem != null) {
- carbonFlux = carbonFluxItem.getCarbonFlux(prevLC, newLC);
- } else {
- carbonFlux = 0.0; // if missing data assume 0 TODO is this right?
- }
-
- if (!aggCFlux.checkForKeys(prevLC, newLC)) { // if not seen yet
- aggCFlux.setCarbonFlux(prevLC, newLC, carbonFlux);
- } else {
- aggCFlux.setCarbonFlux(prevLC, newLC, aggregateMean(aggCFlux.getCarbonFlux(prevLC, newLC), suitableAreaSoFar,
- carbonFlux, suitableAreaThisTime));
- }
+ if (carbonFluxItem != null) {
+ for (Map.Entry<LccKey, Double> cfEntry : carbonFluxItem.getConversionCarbonFluxMap().entrySet()) {
+ LccKey lccKey = cfEntry.getKey();
+ double cFluxThisTime = cfEntry.getValue();
+ double cFluxSoFar = aggCFlux.getConversionCarbonFlux(lccKey);
+ double areaThisTime = landUseItem.getLandCoverArea(lccKey.getFromLc());
+ double areaSoFar = aggLandUse.getLandCoverArea(lccKey.getFromLc());
+ double cFluxAgg = aggregateMean(cFluxSoFar, areaSoFar, cFluxThisTime, areaThisTime);
+ aggCFlux.setConversionCarbonFlux(lccKey, cFluxAgg);
}
- }
+
+ for (Map.Entry<LandCoverType, Double> cfEntry : carbonFluxItem.getEcosystemCarbonFluxMap().entrySet()) {
+ LandCoverType lcType = cfEntry.getKey();
+ double cFluxThisTime = cfEntry.getValue();
+ double cFluxSoFar = aggCFlux.getEcosystemCarbonFlux(lcType);
+ double areaThisTime = landUseItem.getLandCoverArea(lcType);
+ double areaSoFar = aggLandUse.getLandCoverArea(lcType);
+ double cFluxAgg = aggregateMean(cFluxSoFar, areaSoFar, cFluxThisTime, areaThisTime);
+ aggCFlux.setEcosystemCarbonFlux(lcType, cFluxAgg);
+ }
+ }
if (woodYieldItem != null) { // TODO Fix
for (Map.Entry<LccKey, Double> wyEntry : woodYieldItem.getYieldMap().entrySet()) {
LccKey lccKey = wyEntry.getKey();
- double yield = wyEntry.getValue();
- double wYieldThisTime = woodYieldItem.getYield(lccKey);
+ double wYieldThisTime = wyEntry.getValue();
double wYieldSoFar = aggWYield.getYieldLuc(lccKey);
- double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
+ double areaThisTime = landUseItem.getLandCoverArea(lccKey.getFromLc());
+ double areaSoFar = aggLandUse.getLandCoverArea(lccKey.getFromLc());
+ double wYieldAgg = aggregateMean(wYieldSoFar, areaSoFar, wYieldThisTime, areaThisTime);
aggWYield.setYieldLuc(lccKey, wYieldAgg);
}
diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java
index c04c4b23..608ee0aa 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldItem.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java
@@ -37,11 +37,11 @@ public class WoodYieldItem implements RasterItem {
ArrayList<LandUseTile> tiles = landUseTiles.get(key.getFromLc());
- if (tiles.size() == 0) {
+ double totalArea = tiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
+ if (tiles.size() == 0 || totalArea == 0) {
this.yield.put(key, 0.0) ;
} else {
double totalYield = 0;
- double totalArea = 0;
for (LandUseTile tile : tiles) {
if (tile.isProtected())
continue;
@@ -50,6 +50,7 @@ public class WoodYieldItem implements RasterItem {
totalArea += tile.getArea();
}
double meanYield = totalYield / totalArea;
+
this.yield.put(key, meanYield);
}
}
diff --git a/src/ac/ed/lurg/forestry/WoodYieldReader.java b/src/ac/ed/lurg/forestry/WoodYieldReader.java
index d1d1557c..3f7076e4 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldReader.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldReader.java
@@ -31,7 +31,7 @@ public class WoodYieldReader {
return woodYieldData;
}
- public WoodYieldRasterSet getWoodYieldsForest(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
+ public void getWoodYieldsForest(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
AbstractTabularRasterReader<WoodYieldItem> yieldReader = new AbstractTabularRasterReader<WoodYieldItem>(DELIMITER, MIN_COLS, woodYieldData) {
protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
@@ -53,11 +53,9 @@ public class WoodYieldReader {
String filename = getDataDir(ModelConfig.WOOD_YIELD_FRST_FILENAME, timestep);
yieldReader.getRasterDataFromFile(filename);
-
- return woodYieldData;
}
- public WoodYieldRasterSet getWoodYieldsAgri(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
+ public void getWoodYieldsAgri(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
AbstractTabularRasterReader<WoodYieldItem> yieldReader = new AbstractTabularRasterReader<WoodYieldItem>(DELIMITER, MIN_COLS, woodYieldData) {
protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
@@ -78,8 +76,6 @@ public class WoodYieldReader {
String filename = getDataDir(ModelConfig.WOOD_YIELD_AGRI_FILENAME, timestep);
yieldReader.getRasterDataFromFile(filename);
-
- return woodYieldData;
}
private String getDataDir(String filename, Timestep timestep) {
--
GitLab