From 2491d15ce7c7d28c80327a7ea0543d9890229fda Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Tue, 14 Sep 2021 18:18:55 +0100
Subject: [PATCH] Refactored land cover age tracking.
Changes to protected area calculation and update.
---
GAMS/IntExtOpt.gms | 6 +-
debug_config.properties | 16 +-
src/ac/ed/lurg/InternationalMarket.java | 2 +-
src/ac/ed/lurg/ModelMain.java | 38 +-
src/ac/ed/lurg/Timestep.java | 4 +
src/ac/ed/lurg/carbon/CarbonFluxItem.java | 36 +-
src/ac/ed/lurg/country/CountryAgent.java | 24 +-
.../ed/lurg/country/CountryAgentManager.java | 2 -
.../ed/lurg/country/LandCoverChangeItem.java | 8 +-
.../lurg/country/gams/GamsLocationInput.java | 3 -
.../country/gams/GamsLocationOptimiser.java | 19 +-
.../lurg/country/gams/GamsLocationOutput.java | 2 -
.../ed/lurg/country/gams/GamsRasterInput.java | 2 -
.../country/gams/GamsRasterOptimiser.java | 42 +-
src/ac/ed/lurg/forestry/ForestManager.java | 77 ---
src/ac/ed/lurg/forestry/WoodYieldItem.java | 62 +-
.../landuse/InitProtectedAreasReader.java | 18 +
src/ac/ed/lurg/landuse/LandCoverItem.java | 16 +-
src/ac/ed/lurg/landuse/LandCoverTile.java | 156 +++++
src/ac/ed/lurg/landuse/LandUseItem.java | 539 ++++++++----------
src/ac/ed/lurg/landuse/LandUseReader.java | 4 +-
src/ac/ed/lurg/landuse/LandUseTile.java | 64 ---
.../ed/lurg/landuse/ProtectedAreasReader.java | 2 +-
src/ac/ed/lurg/output/LandUseOutputer.java | 4 +-
src/ac/ed/lurg/types/LandCoverType.java | 4 +
25 files changed, 557 insertions(+), 593 deletions(-)
delete mode 100644 src/ac/ed/lurg/forestry/ForestManager.java
create mode 100644 src/ac/ed/lurg/landuse/InitProtectedAreasReader.java
create mode 100644 src/ac/ed/lurg/landuse/LandCoverTile.java
delete mode 100644 src/ac/ed/lurg/landuse/LandUseTile.java
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 2d7b2a04..1fb6a0c0 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -137,6 +137,8 @@ $gdxin
maxNetImport(import_crop) = min(maxNetImport(import_crop), demand(import_crop));
minNetImport(import_crop) = min(minNetImport(import_crop), demand(import_crop));
+
+ woodDemand = 0;
VARIABLES
cropArea(crop, location) total area for crops - Mha
@@ -280,11 +282,11 @@ $gdxin
NEW_FOREST_POTENTIAL_HARVEST(location) .. newForestPotentialHarvest(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location));
- WOOD_HARVEST_NATURAL_CALC(location) .. woodHarvestNat(location) =E= naturalDeforestedArea(location) * woodYieldLuc('natural', 'pasture', location);
+ WOOD_HARVEST_NATURAL_CALC(location) .. woodHarvestNat(location) =E= naturalDeforestedArea(location) * woodYieldLuc('natural', 'natural', location);
WOOD_HARVEST_NATURAL_CONSTRAINT(location) .. naturalDeforestedArea(location) =L= landCoverChange('natural', 'natural', location);
- WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((forested, land_cover_after), landCoverChange(forested, land_cover_after, location) * woodYieldLuc(forested, land_cover_after, location));
+ WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((forested, exc_natural), landCoverChange(forested, exc_natural, location) * woodYieldLuc(forested, exc_natural, location));
WOOD_SUPPLY_CALC .. woodSupply =E= woodFromRotation + sum(location, woodHarvestNat(location) + woodHarvestLuc(location));
diff --git a/debug_config.properties b/debug_config.properties
index ff1a0558..bdbbde9e 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -18,19 +18,19 @@ GENERATE_NEW_YIELD_CLUSTERS=false
YIELD_FILENAME=yield.out
-DEBUG_LIMIT_COUNTRIES=true
-DEBUG_COUNTRY_NAME=Central America
-GAMS_COUNTRY_TO_SAVE=Central America
+DEBUG_LIMIT_COUNTRIES=false
+DEBUG_COUNTRY_NAME=United States of America
+GAMS_COUNTRY_TO_SAVE=United States of America
-INIT_WOOD_PRICE=0.4
+INIT_WOOD_PRICE=0.0
INIT_WOOD_STOCK=1000
INIT_CARBON_PRICE=0.0
FOREST_ROTATION_PERIOD=30
INFRASTRUCTURE_EXPANSION_COST=0.0
MANAGED_FOREST_INCREASE_COST=0.2
MANAGED_FOREST_DECREASE_COST=0.2
-FOREST_MANAGEMENT_COST=0.2
+FOREST_MANAGEMENT_COST=0.25
FOREST_MAX_CHANGE=0.02
-LAND_CHANGE_COST=0.2
-NATURAL_AREA_VALUE=0.01
-VEGETATION_CLEARING_COST=0.25
\ No newline at end of file
+LAND_CHANGE_COST=0.5
+NATURAL_AREA_VALUE=0.0
+VEGETATION_CLEARING_COST=0.05
\ No newline at end of file
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index 06deaee9..cf260102 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -88,7 +88,6 @@ public class InternationalMarket {
double countryNetImports = entry.getValue().getShockedNetImports();
totalProduction.incrementValue(c, (entry.getValue().getProductionExpected()-entry.getValue().getProductionShock()));
- //TODO this doesn't work when both imports and exports > 0
if (countryNetImports > 0)
totalImportCommodities.incrementValue(c, countryNetImports);
else
@@ -138,6 +137,7 @@ public class InternationalMarket {
double totalWoodImport = 0;
double totalWoodExport = 0;
double totalWoodProduction = 0;
+
for (AbstractCountryAgent ca : countryAgents) {
Map<WoodType, WoodUsageData> woodUsageMap = ca.getWoodUsageData();
for (WoodUsageData woodUsage : woodUsageMap.values()) {
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index f2a1adce..f25cddc1 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -7,12 +7,9 @@ import java.io.FileOutputStream;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
-import java.util.ArrayList;
import java.util.Collection;
import java.util.Map;
import java.util.Map.Entry;
-import java.util.Set;
-
import ac.ed.lurg.carbon.CarbonFluxRasterSet;
import ac.ed.lurg.carbon.CarbonFluxReader;
import ac.ed.lurg.country.AbstractCountryAgent;
@@ -30,12 +27,12 @@ import ac.ed.lurg.demand.CalorieManager;
import ac.ed.lurg.demand.DemandManagerFromFile;
import ac.ed.lurg.demand.DemandManagerSSP;
import ac.ed.lurg.demand.ElasticDemandManager;
-import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.forestry.WoodYieldRasterSet;
import ac.ed.lurg.landuse.ConversionCostReader;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.CropUsageReader;
import ac.ed.lurg.landuse.FPUManager;
+import ac.ed.lurg.landuse.InitProtectedAreasReader;
import ac.ed.lurg.landuse.IrrigationConstraintReader;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.IrrigationMaxAmountReader;
@@ -51,7 +48,6 @@ import ac.ed.lurg.landuse.ProtectedAreasReader;
import ac.ed.lurg.landuse.RunOffReader;
import ac.ed.lurg.landuse.WoodUsageData;
import ac.ed.lurg.landuse.WoodUsageReader;
-import ac.ed.lurg.forestry.WoodYieldRasterSet;
import ac.ed.lurg.forestry.WoodYieldReader;
import ac.ed.lurg.output.LandUseOutputer;
import ac.ed.lurg.output.LpjgOutputer;
@@ -85,8 +81,8 @@ public class ModelMain {
private IrrigationRasterSet currentIrrigationData;
private RasterSet<LandUseItem> globalLandUseRaster;
private RasterSet<IntegerRasterItem> clusterIdRaster;
- WoodYieldReader woodYieldReader;
- CarbonFluxReader carbonFluxReader;
+ private WoodYieldReader woodYieldReader;
+ private CarbonFluxReader carbonFluxReader;
public static void main(String[] args) {
@@ -144,7 +140,7 @@ public class ModelMain {
private void doTimestep(Timestep timestep) {
System.gc();
LogWriter.println("Timestep: " + timestep.toString());
-
+
YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
getUpdateIrrigationData(timestep, yieldSurfaces); // updating currentIrrigationData
@@ -162,9 +158,9 @@ public class ModelMain {
double carbonDemand = demandManager.getCarbonDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
double carbonDemandIncrease = (carbonDemand > previousCarbonDemand) ? carbonDemand - previousCarbonDemand: 0;
- CarbonFluxRasterSet carbonFluxData = getCarbonFluxData(timestep);
WoodYieldRasterSet woodYieldData = getWoodYieldData(timestep);
-
+ CarbonFluxRasterSet carbonFluxData = getCarbonFluxData(timestep);
+
Map<LccKey, Double> conversionCosts;
conversionCosts = (ModelConfig.CONVERSION_COSTS_FROM_FILE) ? new ConversionCostReader().read() : new ConversionCostReader().calcFromConfig();
@@ -180,7 +176,7 @@ public class ModelMain {
internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep); // calculate prices
}
internationalMarket.applyPriceShocks(timestep);
-
+
// output results
outputTimestepResults(timestep, globalLandUseRaster);
}
@@ -483,10 +479,12 @@ public class ModelMain {
}
if (ModelConfig.IS_CALIBRATION_RUN) {
- //serializeLandUse(landUseRaster);
- countryAgents.serializeCropUsageForAll();
- countryAgents.serializeWoodUsageForAll();
- internationalMarket.serializeGlobalPrices();
+ if (timestep.isFinalTimestep()) {
+ serializeLandUse(landUseRaster);
+ countryAgents.serializeCropUsageForAll();
+ countryAgents.serializeWoodUsageForAll();
+ internationalMarket.serializeGlobalPrices();
+ }
}
if (timestep.isInitialTimestep() && ModelConfig.GENERATE_NEW_YIELD_CLUSTERS)
@@ -685,13 +683,18 @@ public class ModelMain {
}
};
- new MaxCropAreaReader(initialLC).getRasterDataFromFile(ModelConfig.HIGH_SLOPE_AREAS_FILE);
+ 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
RasterSet<LandUseItem> landUseRaster = new RasterSet<LandUseItem>(initialLC.getHeaderDetails());
for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
+ //LogWriter.println(initialLC.getXCoordin(entry.getKey()) + " " + initialLC.getYCoordin(entry.getKey()));
+ if (entry.getKey().getCol()==38 && entry.getKey().getRow()==60) {
+ int foo = 1;
+ }
landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue()));
}
@@ -717,7 +720,8 @@ public class ModelMain {
/** Get carbon flux data */
private CarbonFluxRasterSet getCarbonFluxData(Timestep timestep) {
- return carbonFluxReader.getCarbonFluxes(globalLandUseRaster, timestep);
+ //return carbonFluxReader.getCarbonFluxes(globalLandUseRaster, timestep);
+ return new CarbonFluxRasterSet(desiredProjection);
}
/** Get wood yield data */
diff --git a/src/ac/ed/lurg/Timestep.java b/src/ac/ed/lurg/Timestep.java
index c92f0b95..0798d353 100644
--- a/src/ac/ed/lurg/Timestep.java
+++ b/src/ac/ed/lurg/Timestep.java
@@ -47,6 +47,10 @@ public class Timestep {
return timestep == ModelConfig.START_TIMESTEP;
}
+ public boolean isFinalTimestep() {
+ return timestep == ModelConfig.END_TIMESTEP;
+ }
+
public Timestep getPreviousTimestep() {
if (timestep == ModelConfig.START_TIMESTEP)
throw new RuntimeException("Can't getPreviousTimestep for " + timestep);
diff --git a/src/ac/ed/lurg/carbon/CarbonFluxItem.java b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
index e3ccfabd..3514e65a 100644
--- a/src/ac/ed/lurg/carbon/CarbonFluxItem.java
+++ b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
@@ -3,12 +3,11 @@ 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.LandCoverTile;
import ac.ed.lurg.landuse.LccKey;
import ac.ed.lurg.types.LandCoverType;
@@ -22,21 +21,16 @@ public class CarbonFluxItem implements RasterItem, Serializable {
conversionCarbonFlux.put(key, cFlux);
}
- public void calcConversionCarbonFlux(LccKey key, Double[] cFluxes, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep) {
+ public void calcConversionCarbonFlux(LccKey key, Double[] cFluxes, Map<LandCoverType, LandCoverTile> landUseTiles, Timestep timestep) {
+ LandCoverTile tiles = landUseTiles.get(key.getFromLc());
- ArrayList<LandUseTile> tiles = landUseTiles.get(key.getFromLc());
-
- if (tiles.size() == 0) {
+ double totalArea = tiles.getTotalConvertibleArea();
+ if (totalArea == 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();
+ for (int age=0; age<=tiles.getMaxAgeBin(); age++) {
+ totalFlux += cFluxes[age] * tiles.getConvertibleArea(age);
}
double meanFlux = totalFlux / totalArea;
this.conversionCarbonFlux.put(key, meanFlux);
@@ -47,21 +41,17 @@ public class CarbonFluxItem implements RasterItem, Serializable {
ecosystemCarbonFlux.put(lcType, cFlux);
}
- public void calcEcosystemCarbonFlux(LandCoverType lcType, Double[] cFluxes, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep) {
+ public void calcEcosystemCarbonFlux(LandCoverType lcType, Double[] cFluxes, Map<LandCoverType, LandCoverTile> landUseTiles, Timestep timestep) {
- ArrayList<LandUseTile> tiles = landUseTiles.get(lcType);
+ LandCoverTile tiles = landUseTiles.get(lcType);
- if (tiles.size() == 0) {
+ double totalArea = tiles.getTotalConvertibleArea();
+ if (totalArea == 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();
+ for (int age=0; age<=tiles.getMaxAgeBin(); age++) {
+ totalFlux += cFluxes[age] * tiles.getConvertibleArea(age);
}
double meanFlux = totalFlux / totalArea;
this.ecosystemCarbonFlux.put(lcType, meanFlux);
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 8a4136c7..a581b02a 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -28,7 +28,6 @@ import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.WoodType;
-import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.utils.cluster.Cluster;
import ac.ed.lurg.utils.cluster.KMeans;
@@ -78,7 +77,7 @@ public class CountryAgent extends AbstractCountryAgent {
LandUseItem luItem = previousGamsRasterOutput.getLandUses().get(key);
double totalLand = luItem.getTotalLandCoverArea()-luItem.getLandCoverArea(LandCoverType.URBAN);
- double protectedAreaFrac = (totalLand <= 0) ? 0.0 : luItem.getProtectedArea() / totalLand;
+ double protectedAreaFrac = (totalLand <= 0) ? 0.0 : luItem.getTotalProtectedArea() / totalLand;
clusteringPoints.add(new YieldClusterPoint(key, yieldresp, irrigItem, protectedAreaFrac));
}
@@ -141,8 +140,10 @@ public class CountryAgent extends AbstractCountryAgent {
if (saveGamsGdxFiles)
saveGDXFile("landuse");
-
+
previousGamsRasterOutput = result;
+
+ incrementLandCoverAge();
return result;
}
@@ -343,28 +344,29 @@ public class CountryAgent extends AbstractCountryAgent {
public void updateProtectedAreas(RasterSet<LandUseItem> newProtectedAreas) {
for (Entry<RasterKey, LandUseItem> entry : previousGamsRasterOutput.getLandUses().entrySet()) {
- double updatedPA = newProtectedAreas.get(entry.getKey()).getProtectedArea();
- entry.getValue().setProtectedArea(updatedPA);
+ double updatedProtectedFract = newProtectedAreas.get(entry.getKey()).getProtectedFraction();
+ entry.getValue().setProtectedFraction(updatedProtectedFract);
}
}
private double harvestWoodFromRotation(RasterSet<WoodYieldItem> woodYieldData) {
- //RasterSet<WoodYieldItem> woodYieldDataSubset = woodYieldData.createSubsetForKeys(previousGamsRasterOutput.getLandUses().keySet());
- //return woodYieldDataSubset.values().stream().mapToDouble(o -> o.getCurrentRotationHarvest()).sum();
double totalHarvest = 0;
- for (Map.Entry<RasterKey, LandUseItem> entry : previousGamsRasterOutput.getLandUses().entrySet()) {
- RasterKey key = entry.getKey();
- LandUseItem luItem = entry.getValue();
+ for (RasterKey key : previousGamsRasterOutput.getLandUses().keySet()) {
WoodYieldItem wyItem = woodYieldData.get(key);
if (wyItem == null)
continue; // TODO Deal with this properly
totalHarvest += wyItem.getCurrentRotationHarvest();
- luItem.doForestRotation(currentTimestep, wyItem.getOptimalRotation()); // WARNING: Side effect
}
return totalHarvest;
}
+ private void incrementLandCoverAge() {
+ for (LandUseItem luItem : previousGamsRasterOutput.getLandUses().values()) {
+ luItem.incrementTilesAge();
+ }
+ }
+
public double getNetCarbonFlux() {
return previousGamsRasterOutput.getNetCarbonFlux();
}
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index b5951f96..a7408a40 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -26,9 +26,7 @@ import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.LccKey;
import ac.ed.lurg.landuse.WoodUsageData;
import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.WoodType;
-import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.YieldRaster;
import ac.sac.raster.IntegerRasterItem;
diff --git a/src/ac/ed/lurg/country/LandCoverChangeItem.java b/src/ac/ed/lurg/country/LandCoverChangeItem.java
index 5e06a36f..348b93e2 100644
--- a/src/ac/ed/lurg/country/LandCoverChangeItem.java
+++ b/src/ac/ed/lurg/country/LandCoverChangeItem.java
@@ -3,24 +3,18 @@ package ac.ed.lurg.country;
import ac.ed.lurg.types.LandCoverType;
public class LandCoverChangeItem {
- private int location;
private LandCoverType fromLandCover;
private LandCoverType toLandCover;
private double area;
public LandCoverChangeItem() {}
- public LandCoverChangeItem(int location, LandCoverType fromLandCover, LandCoverType toLandCover, double area) {
- this.location = location;
+ public LandCoverChangeItem(LandCoverType fromLandCover, LandCoverType toLandCover, double area) {
this.fromLandCover = fromLandCover;
this.toLandCover = toLandCover;
this.area = area;
}
- public int getLocation() {
- return location;
- }
-
public LandCoverType getFromLandCover() {
return fromLandCover;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index a1e398cc..211b4974 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -5,12 +5,9 @@ import java.util.Map;
import ac.ed.lurg.Timestep;
import ac.ed.lurg.carbon.CarbonFluxItem;
import ac.ed.lurg.forestry.WoodYieldData;
-import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.LccKey;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.yield.YieldResponsesItem;
public class GamsLocationInput {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 0569bf7d..7fe7179a 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -172,7 +172,7 @@ public class GamsLocationOptimiser {
Vector<String> v = new Vector<String>();
v.add(lc.getName());
v.add(locString);
- setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(lc), 4);
+ setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getConvertibleLandCoverArea(lc), 4);
}
// Infrastructure cost
@@ -569,15 +569,7 @@ public class GamsLocationOptimiser {
double croplandArea = getParmValue(parmCroplandArea, locationName);
- if (cropType.equals(CropType.PASTURE)) {
- landUseItem.setLandCoverArea(LandCoverType.PASTURE, area);
- totalPastureArea += area;
- }
- else {
- landUseItem.setLandCoverArea(LandCoverType.CROPLAND, croplandArea); // will set this multiple times, once for each arable crop, but doesn't really matter
- totalCropArea += area;
- }
-
+
landUseItem.setCropFraction(cropType, croplandArea > 0 ? area/croplandArea : 0); // TODO gives wrong pasture numbers
}
@@ -606,7 +598,7 @@ public class GamsLocationOptimiser {
int locId = Integer.parseInt(locationName);
double change = rec.getLevel();
- // Skip if no change
+ // Skip if no change or no conversion
if (change == 0 || fromLcStr.equals(toLcStr))
continue;
@@ -614,7 +606,7 @@ public class GamsLocationOptimiser {
LandCoverType toLc = LandCoverType.getForName(toLcStr);
ArrayList<LandCoverChangeItem> changesList = landCoverChanges.computeIfAbsent(locId, k -> new ArrayList<LandCoverChangeItem>());
- changesList.add(new LandCoverChangeItem(locId, fromLc, toLc, change));
+ changesList.add(new LandCoverChangeItem(fromLc, toLc, change));
}
// Carbon flux
@@ -635,7 +627,8 @@ public class GamsLocationOptimiser {
double previousNetImport = previousWoodUsageMap.get(wt).getNetImport();
double newNetImport = previousNetImport + netWoodImportChange / numWoodTypes;
double previousHarvest = previousWoodUsageMap.get(wt).getHarvest();
- double newHarvest = totalWoodHarvest * (previousHarvest / previousWoodHarvestSum);
+ // assuming equal split if no previous harvest
+ double newHarvest = (previousWoodHarvestSum > 0) ? totalWoodHarvest * (previousHarvest / previousWoodHarvestSum) : totalWoodHarvest / numWoodTypes;
newWoodUsageMap.put(wt, new WoodUsageData(newHarvest, newNetImport));
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index edce072f..08bb3f85 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -9,8 +9,6 @@ import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.WoodUsageData;
import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.utils.TripleMap;
-import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.WoodType;
import ac.ed.lurg.country.LandCoverChangeItem;
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index 71f255e3..ddc7ee8d 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -7,8 +7,6 @@ import ac.ed.lurg.carbon.CarbonFluxItem;
import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.LccKey;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.yield.YieldRaster;
import ac.sac.raster.RasterSet;
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 1f675e9a..fc6c36c4 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -11,7 +11,6 @@ import ac.ed.lurg.carbon.CarbonFluxItem;
import ac.ed.lurg.country.LandCoverChangeItem;
import ac.ed.lurg.forestry.WoodYieldData;
import ac.ed.lurg.forestry.WoodYieldItem;
-import ac.ed.lurg.forestry.WoodYieldRasterSet;
import ac.ed.lurg.landuse.Intensity;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandUseItem;
@@ -19,7 +18,6 @@ import ac.ed.lurg.landuse.LccKey;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.YieldType;
-import ac.ed.lurg.utils.DoubleMap;
import ac.ed.lurg.utils.LazyTreeMap;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.YieldRaster;
@@ -63,7 +61,7 @@ public class GamsRasterOptimiser {
}
private RasterSet<LandUseItem> createWithSameLandCovers(RasterSet<LandUseItem> toCopy) {
- RasterSet<LandUseItem> theCopy = new RasterSet<LandUseItem>(toCopy.getHeaderDetails());
+ RasterSet<LandUseItem> theCopy = new RasterSet<LandUseItem>(toCopy.getHeaderDetails()); // TODO why do we make a copy???
for (Entry<RasterKey, LandUseItem> entry : toCopy.entrySet()) {
if (entry.getValue() != null) {
@@ -82,7 +80,7 @@ public class GamsRasterOptimiser {
for (LandUseItem a : areaRaster.values()) {
if (a != null) {
total += a.getLandCoverArea(l);
- unprotected += a.getUnprotectedLandCoverArea(l);
+ unprotected += a.getConvertibleLandCoverArea(l);
}
}
@@ -93,7 +91,7 @@ public class GamsRasterOptimiser {
double suitableArea=0, protectedArea=0;
for (LandUseItem a : areaRaster.values()) {
if (a != null) {
- protectedArea += a.getProtectedArea();
+ protectedArea += a.getTotalProtectedArea();
suitableArea += a.getSuitableArea();
}
}
@@ -184,7 +182,7 @@ public class GamsRasterOptimiser {
for (RasterKey key : keys) {
LandUseItem newLandUseItem = newLandUseRaster.get(key);
- totalUnprotectedLC += newLandUseItem.getUnprotectedLandCoverArea(fromType);
+ totalUnprotectedLC += newLandUseItem.getConvertibleLandCoverArea(fromType);
}
for (RasterKey key : keys) {
@@ -192,9 +190,8 @@ public class GamsRasterOptimiser {
LandUseItem newLandUseItem = newLandUseRaster.get(key);
if (newLandUseItem!=null) {
- double cellChange = (totalUnprotectedLC > 0) ? change * newLandUseItem.getUnprotectedLandCoverArea(fromType)/totalUnprotectedLC : change / keys.size();
+ double cellChange = (totalUnprotectedLC > 0) ? change * newLandUseItem.getConvertibleLandCoverArea(fromType)/totalUnprotectedLC : change / keys.size();
double shortfall = newLandUseItem.moveAreas(toType, fromType, cellChange);
- shortfall += newLandUseItem.moveUnprotectedArea(toType, fromType, cellChange, rasterInputData.getTimestep());
if (shortfall == 0)
keysWithSpace.add(key);
else
@@ -273,7 +270,6 @@ public class GamsRasterOptimiser {
IrrigationItem irrigItem = irrigRaster.get(key);
WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
CarbonFluxItem carbonFluxItem = carbonFluxRaster.get(key);
-// WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
int clusterId = mapping.get(key).getInt();
@@ -392,30 +388,12 @@ public class GamsRasterOptimiser {
Intensity intensityThisTime = landUseItem.getIntensity(crop);
aggLandUse.setIntensity(crop, intensityThisTime); // intensity currently is always the same within a location, so can just that any. If this changes need to average here.
}
-
- // Unavailable Fraction
- double unavailAreaThisTime = landUseItem.getUnavailableArea();
- double unavailAreaSoFar = aggLandUse.getUnavailableArea();
- aggLandUse.setUnavailableArea(unavailAreaSoFar + unavailAreaThisTime);
-
- // Protected areas
- double protectedThisTime = landUseItem.getProtectedArea();
- double protectedSoFar = aggLandUse.getProtectedArea();
- aggLandUse.setProtectedArea(protectedSoFar + protectedThisTime);
-
- //Suitable areas
- aggLandUse.setSuitableArea(suitableAreaSoFar + suitableAreaThisTime);
-
// Land covers areas
for (LandCoverType landType : LandCoverType.values()) {
- double unprotectedAreaThisTime = landUseItem.getUnprotectedLandCoverArea(landType);
- double unprotectedAreaSoFar = aggLandUse.getUnprotectedLandCoverArea(landType);
- aggLandUse.setUnprotectedLandCoverArea(landType, unprotectedAreaSoFar + unprotectedAreaThisTime);
-
- double totalAreaThisTime = landUseItem.getLandCoverArea(landType);
- double totalAreaSoFar = aggLandUse.getLandCoverArea(landType);
- aggLandUse.setLandCoverArea(landType, totalAreaSoFar + totalAreaThisTime);
+ double totalAreaThisTime = landUseItem.getConvertibleLandCoverArea(landType);
+ double totalAreaSoFar = aggLandUse.getConvertibleLandCoverArea(landType);
+ aggLandUse.setConvertibleLandCoverArea(landType, totalAreaSoFar + totalAreaThisTime);
}
}
@@ -429,8 +407,10 @@ public class GamsRasterOptimiser {
// TODO shouldn't have missing data in the first place
double meanRotation = aggregatedForestRotation.values().stream().mapToDouble(o -> o.doubleValue()).sum() / aggregatedForestRotation.keySet().size();
for (Integer locId : aggregatedYields.keySet()) {
- if (!aggregatedForestRotation.containsKey(locId))
+ if (!aggregatedForestRotation.containsKey(locId)) {
aggregatedForestRotation.put(locId, meanRotation);
+ LogWriter.printlnWarning("GamsRasterOptimister.convertFromRaster missing timber rotation; substituting average");
+ }
}
return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts,
diff --git a/src/ac/ed/lurg/forestry/ForestManager.java b/src/ac/ed/lurg/forestry/ForestManager.java
deleted file mode 100644
index 3027c3f8..00000000
--- a/src/ac/ed/lurg/forestry/ForestManager.java
+++ /dev/null
@@ -1,77 +0,0 @@
-package ac.ed.lurg.forestry;
-
-import java.io.Serializable;
-import java.util.ArrayList;
-import java.util.Collections;
-import java.util.HashMap;
-import java.util.Map;
-
-import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.Timestep;
-import ac.ed.lurg.carbon.CarbonFluxItem;
-import ac.ed.lurg.landuse.LandUseTile;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.LogWriter;
-
-public class ForestManager implements Serializable {
- private static final long serialVersionUID = 3093028776147903435L;
-
- private WoodYieldItem woodYields;
- private int optimalRotation;
-
- public ForestManager() {
- this.woodYields = new WoodYieldItem();
- }
- /*
- public void setWoodYield(LandCoverType forestType, int year, int age, double yield) {
- woodYields.setYield(forestType, year, age, yield);
- }
-
- public void copyWoodYields(WoodYieldItem woodYieldItem) {
- woodYields.copyAll(woodYieldItem);
- }
-
- public void calculateOptimalRotation(double woodPrice, double discountRate, Timestep timestep) {
- Map<Double, Integer> levMap = new HashMap<Double, Integer>(); // Land Expected Value
- int maxAge = WoodYieldItem.MAX_AGE;
- int year = timestep.getYear();
-
- for (int age = 0; age <= maxAge; age++) {
- double yield = woodYields.getYield(LandCoverType.TIMBER_FOREST, year, age);
- double lev = (woodPrice * yield * Math.exp(-discountRate * age) - ModelConfig.FOREST_ESTABLISHMENT_COST) / (1 - Math.exp(-discountRate * age));
- levMap.put(lev, age);
- }
- // TODO potential edge case with tied LEV
- double maxLev = Collections.max(levMap.keySet());
- optimalRotation = levMap.get(maxLev);
- }
-
- public int getOptimalRotation() {
- return optimalRotation;
- }
-
- public double getStandYield(LandUseTile tile, Timestep timestep) {
- return woodYields.getYield(tile.getLcType(), timestep.getYear(), tile.getAge(timestep));
- }
-
- public double getAverageYield(ArrayList<LandUseTile> forestTiles, Timestep timestep) {
- if (forestTiles == null)
- return 0;
- Double totalWoodStock = forestTiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea() * getStandYield(o, timestep)).sum();
- if (totalWoodStock > 0) {
- int foo = 1;
- }
- double totalArea = forestTiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
- if (totalWoodStock == null || totalWoodStock == 0) {
- return 0.0;
- } else {
- return totalWoodStock / totalArea;
- }
- }
-
- public double getYieldAtRotation(int year) {
- return woodYields.getYield(LandCoverType.TIMBER_FOREST, year, optimalRotation);
- }
- */
-
-}
diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java
index 608ee0aa..e0de28c4 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldItem.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java
@@ -3,14 +3,14 @@ package ac.ed.lurg.forestry;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
+import java.util.List;
import java.util.Map;
-import java.util.stream.Collectors;
-
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.Timestep;
-import ac.ed.lurg.landuse.LandUseTile;
+import ac.ed.lurg.landuse.LandCoverTile;
import ac.ed.lurg.landuse.LccKey;
import ac.ed.lurg.types.LandCoverType;
+import ac.ed.lurg.utils.LogWriter;
import ac.sac.raster.RasterItem;
public class WoodYieldItem implements RasterItem {
@@ -29,25 +29,21 @@ public class WoodYieldItem implements RasterItem {
public WoodYieldItem() {}
- public void calcYieldData(Map<LccKey, Double[]> woodYields, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep) {
+ public void calcYieldData(Map<LccKey, Double[]> woodYields, Map<LandCoverType, LandCoverTile> landUseTiles, Timestep timestep) {
// Mean wood yield for grid cell
for (Map.Entry<LccKey, Double[]> entry : woodYields.entrySet()) {
LccKey key = entry.getKey();
Double[] yields = entry.getValue();
- ArrayList<LandUseTile> tiles = landUseTiles.get(key.getFromLc());
+ LandCoverTile tiles = landUseTiles.get(key.getFromLc());
- double totalArea = tiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
- if (tiles.size() == 0 || totalArea == 0) {
+ double totalArea = tiles.getTotalConvertibleArea(); // Assuming no harvest from protected areas
+ if (totalArea == 0) {
this.yield.put(key, 0.0) ;
} else {
double totalYield = 0;
- for (LandUseTile tile : tiles) {
- if (tile.isProtected())
- continue;
- int age = tile.getAge(timestep);
- totalYield += yields[age] * tile.getArea();
- totalArea += tile.getArea();
+ for (int age=0; age<=tiles.getMaxAgeBin(); age++) {
+ totalYield += yields[age] * tiles.getConvertibleArea(age);
}
double meanYield = totalYield / totalArea;
@@ -56,31 +52,45 @@ public class WoodYieldItem implements RasterItem {
}
}
- public void calcRotationData(Map<LccKey, Double[]> woodYields, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep, double woodPrice) {
+ public void calcRotationData(Map<LccKey, Double[]> woodYields, Map<LandCoverType, LandCoverTile> landUseTiles, Timestep timestep, double woodPrice) {
// Optimal rotation
- Map<Double, Integer> levMap = new HashMap<Double, Integer>(); // Land Expected Value
+ List<Double> levArr = new ArrayList<Double>();
Double[] yields = woodYields.get(new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST));
for (int age = 0; age <= WoodYieldItem.MAX_AGE; age++) {
double yield = yields[age];
double lev = (woodPrice * yield * Math.exp(-ModelConfig.DISCOUNT_RATE * age) - ModelConfig.FOREST_ESTABLISHMENT_COST) / (1 - Math.exp(-ModelConfig.DISCOUNT_RATE * age));
- levMap.put(lev, age);
+ levArr.add(lev);
+ }
+
+ double maxLev = Collections.max(levArr);
+ List<Integer> candidates = new ArrayList<Integer>(); // There may be several equal maximum values
+ for (int i = 0; i < levArr.size(); i ++) {
+ if (levArr.get(i) == maxLev) {
+ candidates.add(i);
+ }
+ }
+
+ try {
+ optimalRotation = Collections.min(candidates); // Choose shortest rotation
+ } catch(Exception e) {
+ LogWriter.print(e);
}
- // TODO potential edge case with tied LEV
- double maxLev = Collections.max(levMap.keySet());
- optimalRotation = levMap.get(maxLev);
yieldAtRotation = yields[optimalRotation];
- ArrayList<LandUseTile> timberTiles = landUseTiles.get(LandCoverType.TIMBER_FOREST);
- if (timberTiles.size() == 0) {
+ LandCoverTile timberTiles = landUseTiles.get(LandCoverType.TIMBER_FOREST);
+ double timberForestArea = timberTiles.getTotalConvertibleArea();
+ if (timberForestArea == 0) {
currentRotationHarvest = 0.0;
} else {
- double areaHarvested = timberTiles.stream()
- .filter(o -> !o.isProtected() && o.getAge(timestep) >= optimalRotation)
- .mapToDouble(o -> o.getArea())
- .sum();
- currentRotationHarvest = areaHarvested * yieldAtRotation;
+ double areaHarvested = 0;
+ for (int age=optimalRotation; age<=timberTiles.getMaxAgeBin(); age++) {
+ areaHarvested += timberTiles.getConvertibleArea(age);
+ timberTiles.resetAreaForAge(age); // Warning side effect
+ }
+ currentRotationHarvest = areaHarvested * yieldAtRotation; // TODO how to better deal with age > rotationAge at initialisation?
+
}
}
diff --git a/src/ac/ed/lurg/landuse/InitProtectedAreasReader.java b/src/ac/ed/lurg/landuse/InitProtectedAreasReader.java
new file mode 100644
index 00000000..24774e43
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/InitProtectedAreasReader.java
@@ -0,0 +1,18 @@
+package ac.ed.lurg.landuse;
+
+import ac.sac.raster.AbstractRasterReader;
+import ac.sac.raster.RasterSet;
+
+public class InitProtectedAreasReader extends AbstractRasterReader<LandCoverItem> {
+ public InitProtectedAreasReader (RasterSet<LandCoverItem> dataset) {
+ super(dataset);
+ }
+
+ @Override
+ public void setData(LandCoverItem lcData, String token) {
+ if (!"nan".equals(token)) {
+ double protFrac = Double.parseDouble(token);
+ lcData.setProtectedFraction(protFrac);
+ }
+ }
+}
diff --git a/src/ac/ed/lurg/landuse/LandCoverItem.java b/src/ac/ed/lurg/landuse/LandCoverItem.java
index 767624f5..66936e43 100644
--- a/src/ac/ed/lurg/landuse/LandCoverItem.java
+++ b/src/ac/ed/lurg/landuse/LandCoverItem.java
@@ -1,6 +1,5 @@
package ac.ed.lurg.landuse;
-import java.util.ArrayList;
import java.util.HashMap;
import java.util.Map;
@@ -11,11 +10,12 @@ import ac.sac.raster.RasterItem;
* This is used in the initalisation phase, after that land-use is used */
public class LandCoverItem implements RasterItem {
- private Map<LandCoverType, Double> landcover = new HashMap<LandCoverType, Double>();
+ private Map<LandCoverType, Double> landcover = new HashMap<LandCoverType, Double>(); // Land cover fraction
private Map<LandCoverType, Map<Integer, Double>> landUseAgeDist = new HashMap<LandCoverType, Map<Integer, Double>>();
private double totalArea;
private double unavailableFract; // due to slope
private double initialForestedNaturalFract;
+ private double protectedFraction;
/** Area in Mha */
public Double getLandCoverArea(LandCoverType landType) {
@@ -27,6 +27,10 @@ public class LandCoverItem implements RasterItem {
return d==null ? 0 : d.doubleValue();
}
+ public Map<LandCoverType, Double> getLandCoverMap() {
+ return landcover;
+ }
+
/** Area in Mha */
public double getTotalArea() {
return totalArea;
@@ -65,4 +69,12 @@ public class LandCoverItem implements RasterItem {
public Map<LandCoverType, Map<Integer, Double>> getLandUseAgeDist() {
return landUseAgeDist;
}
+
+ public double getProtectedFraction() {
+ return protectedFraction;
+ }
+
+ public void setProtectedFraction(double protectedFraction) {
+ this.protectedFraction = protectedFraction;
+ }
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/landuse/LandCoverTile.java b/src/ac/ed/lurg/landuse/LandCoverTile.java
new file mode 100644
index 00000000..3f0bdb38
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/LandCoverTile.java
@@ -0,0 +1,156 @@
+package ac.ed.lurg.landuse;
+
+import java.io.Serializable;
+import java.util.Arrays;
+
+import ac.ed.lurg.utils.LogWriter;
+import ac.sac.raster.InterpolatingRasterItem;
+
+public class LandCoverTile implements Serializable, InterpolatingRasterItem<LandCoverTile> {
+ private static final long serialVersionUID = 3208516181615267472L;
+ private static final int MAX_AGE = 250;
+
+ private Double[] convertibleAreas; // area Mha; index = age
+ private Double[] protectedAreas; // area Mha; index = age
+ private double unavailableArea; // unavailable for conversion due to slope, altitude TODO age?
+
+ public LandCoverTile() {
+ this.convertibleAreas = new Double[MAX_AGE+1];
+ this.protectedAreas = new Double[MAX_AGE+1];
+ Arrays.fill(convertibleAreas, 0.0);
+ Arrays.fill(protectedAreas, 0.0);
+ unavailableArea = 0.0;
+ }
+
+ public double getConvertibleArea(int age) {
+ return convertibleAreas[age];
+ }
+
+ public void setConvertibleArea(int age, double area) {
+ convertibleAreas[age] = area;
+ }
+
+ public void addConvertibleArea(double a) {
+ convertibleAreas[0] += a;
+ }
+
+ public double getTotalConvertibleArea() {
+ double total = 0;
+ for (double d : convertibleAreas) {
+ total += d;
+ }
+ return total;
+ }
+
+ public void removeConvertibleArea(double a) { // Subtracts area uniformly
+ if (a <= 0)
+ return;
+ double totalArea = getTotalConvertibleArea();
+ double areaToRemove = a;
+ if (a > totalArea) {
+ LogWriter.printlnError("LandUseTile.removeConvertibleArea: attempting to remove too much area: " + a + ", available: " + totalArea);
+ areaToRemove = totalArea;
+ }
+ double factor = (totalArea - areaToRemove) / totalArea;
+ for (int i=0; i<convertibleAreas.length; i++) {
+ convertibleAreas[i] *= factor;
+ }
+ }
+
+ public void setProtectedArea(int age, double area) {
+ protectedAreas[age] = area;
+ }
+
+ public void addProtectedArea(double d) {
+ protectedAreas[0] += d;
+
+ }
+
+ public double getTotalProtectedArea() {
+ double total = 0;
+ for (double d : protectedAreas) {
+ total += d;
+ }
+ return total;
+ }
+
+ public void removeProtectedArea(double a) { // Subtracts area uniformly
+ if (a <= 0)
+ return;
+ double totalArea = getTotalProtectedArea();
+ double areaToRemove = a;
+ if (a > totalArea) {
+ LogWriter.printlnError("LandUseTile.removeProtectedArea: attempting to remove too much area: " + a + ", available: " + totalArea);
+ areaToRemove = totalArea;
+ }
+ double factor = (totalArea - areaToRemove) / totalArea;
+ for (int i=0; i<protectedAreas.length; i++) {
+ protectedAreas[i] *= factor;
+ }
+ }
+
+ public void addUnavailableArea(double d) {
+ unavailableArea += d;
+ }
+
+ public double getTotalArea() {
+ return getTotalConvertibleArea() + getTotalProtectedArea() + unavailableArea;
+ }
+
+ public void resetAreaForAge(int age) {
+ convertibleAreas[0] += convertibleAreas[age];
+ convertibleAreas[age] = 0.0;
+ }
+
+ public void removeAllArea() {
+ Arrays.fill(convertibleAreas, 0.0);
+ Arrays.fill(protectedAreas, 0.0);
+ unavailableArea = 0.0;
+ }
+
+ public void incrementAge() {
+ int n = convertibleAreas.length;
+ convertibleAreas[n-1] += convertibleAreas[n-2]; // final age bin is all absorbing
+ convertibleAreas[n-2] = 0.0;
+ // shift values
+ for (int i=(n-2); i>0; i--) {
+ convertibleAreas[i] = convertibleAreas[i-1];
+ }
+ convertibleAreas[0]= 0.0;
+
+ int m = protectedAreas.length;
+ protectedAreas[m-1] += protectedAreas[m-2]; // final age bin is all absorbing
+ protectedAreas[m-2] = 0.0;
+ // shift values
+ for (int i=(m-2); i>0; i--) {
+ protectedAreas[i] = protectedAreas[i-1];
+ }
+ protectedAreas[0]= 0.0;
+ }
+
+ // Moves area from convertible to protected, maintaining age classes
+ public void moveAreaToProtected(double area) {
+ double convertibleA = getTotalConvertibleArea();
+ if (area > convertibleA) {
+ LogWriter.printlnWarning("LandUseTile.moveAreaToProtected cannot protect all required area");
+ area = Math.min(convertibleA, area);
+ }
+ double fractToMove = area / convertibleA;
+ for (int i = 0; i <= MAX_AGE; i ++) {
+ double areaToMove = convertibleAreas[i] * fractToMove;
+ protectedAreas[i] += areaToMove;
+ convertibleAreas[i] -= areaToMove;
+ }
+ }
+
+ public int getMaxAgeBin() {
+ return MAX_AGE;
+ }
+
+ @Override
+ public void interpolateAll(LandCoverTile fromItem, LandCoverTile toItem, double factor) {
+ // TODO Auto-generated method stub
+
+ }
+
+}
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index a378194e..a00e2419 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -1,17 +1,10 @@
package ac.ed.lurg.landuse;
import java.io.Serializable;
-import java.util.ArrayList;
import java.util.Collection;
import java.util.HashMap;
-import java.util.HashSet;
import java.util.Map;
-import java.util.Set;
-
import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.Timestep;
-import ac.ed.lurg.forestry.WoodYieldItem;
-import ac.ed.lurg.forestry.ForestManager;
import ac.ed.lurg.types.CropToDouble;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
@@ -24,27 +17,24 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
private Map<CropType, Intensity> intensityMap = new HashMap<CropType, Intensity>();
private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
- private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>(); //Mha
- private Map<LandCoverType, Double> unprotectedAreas = new HashMap<LandCoverType, Double>(); //Mha
- private Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles = new HashMap<LandCoverType, ArrayList<LandUseTile>>();
- private double protectedArea; //protected area in Mha
- private double unavailableArea; //area unavailable due to altitude etc
- private double suitableArea; //Mha
+ private Map<LandCoverType, LandCoverTile> landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
+ private double protectedFraction; // protected fraction of total area
- public LandUseItem() {}
+ public LandUseItem() {
+ for (LandCoverType lcType : LandCoverType.values()) {
+ landCoverAreas.put(lcType, new LandCoverTile());
+ }
+ }
public LandUseItem(LandCoverItem landCover) {
this();
if (landCover != null) {
- for (LandCoverType lcType : LandCoverType.values())
- landCoverAreas.put(lcType, landCover.getLandCoverArea(lcType)); // Converts land cover fractions to areas
-
setCropFraction(CropType.WHEAT, 0.5); // random start as don't have better data
setCropFraction(CropType.MAIZE, 0.5);
+ setProtectedFraction(landCover.getProtectedFraction());
+ addLandCoverTiles(landCover);
setUnavailableArea(landCover.getUnavailableFract());
- updateSuitableArea(ModelConfig.BASE_YEAR);
- addLandUseTiles(landCover.getLandUseAgeDist()); // Should only be run after unprotected areas calculated
- runAreaCheck(); // check that tile areas add up to total areas
+
}
}
@@ -53,80 +43,42 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
intensityMap.putAll(luItemToCopy.intensityMap);
cropFractions.putAll(luItemToCopy.cropFractions);
landCoverAreas.putAll(luItemToCopy.landCoverAreas);
- unprotectedAreas.putAll(luItemToCopy.unprotectedAreas);
- landUseTiles.putAll(luItemToCopy.landUseTiles);
- protectedArea = (luItemToCopy.protectedArea);
- suitableArea = (luItemToCopy.suitableArea);
+ protectedFraction = (luItemToCopy.protectedFraction);
}
-
-
- public void moveTilesArea(LandCoverType toType, LandCoverType fromType, double area, Timestep timestep) {
- // Note if toType == fromType, area will be removed from fromType and a new tile will be created with age 0
- if (area == 0)
- return;
-
- ArrayList<LandUseTile> toTiles = landUseTiles.get(toType);
- ArrayList<LandUseTile> fromTiles = landUseTiles.get(fromType);
-
- if (fromTiles == null) {
- LogWriter.printlnError("LandUseItem: cannot move area between tiles as source is null");
- return;
+ public void addLandCoverTiles(LandCoverItem landCover) {
+ double protectableFraction = 0;
+ for (LandCoverType lcType : LandCoverType.getProtectibleTypes()) {
+ protectableFraction += landCover.getLandCoverFract(lcType);
}
-
- // Remove area from tiles
- double totalCurrentArea = fromTiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
- for (LandUseTile s: fromTiles) {
- if (s.isProtected())
+ // need to scale protectedFraction by protectableFraction and cap at 1.0. Gives fraction of protectable LC types to protect
+ double adjProtectedFraction = Math.min(1.0, protectedFraction / protectableFraction);
+
+ for (LandCoverType lcType : LandCoverType.values()) {
+ LandCoverTile tiles = landCoverAreas.get(lcType);
+ double totalArea = landCover.getLandCoverArea(lcType);
+ if (totalArea == 0)
continue;
- double areaToRemove = s.getArea() / totalCurrentArea * area;
- s.removeArea(areaToRemove);
- }
-
- // Add new tiles
- LandUseTile newStand = new LandUseTile(toType, area, timestep.getYear(), false);
- toTiles.add(newStand);
-
- runAreaCheck();
- }
-
- public void runAreaCheck() {
- double threshold = 1e-10;
- for (Map.Entry<LandCoverType, ArrayList<LandUseTile>> entry : landUseTiles.entrySet()) {
- LandCoverType lcType = entry.getKey();
- ArrayList<LandUseTile> tiles = entry.getValue();
- double totalArea = tiles.stream().mapToDouble(o -> o.getArea()).sum();
- double totalUnprotectedArea = tiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
- if (totalArea - landCoverAreas.get(lcType) > threshold) {
- LogWriter.printlnError("LandUseItem: total area and tiles area different: " + lcType.getName() + " " + landCoverAreas.get(lcType) + ", " + totalArea);
+ // static land covers
+ if (!lcType.isConvertible()) {
+ tiles.addConvertibleArea(totalArea); // Assumes land cover age is 0
+ continue;
}
- if (totalUnprotectedArea - unprotectedAreas.get(lcType) > threshold) {
- LogWriter.printlnError("LandUseItem: unprotected area and tiles area different: " + lcType.getName() + " " + unprotectedAreas.get(lcType) + ", " + totalUnprotectedArea);
+ // convertible land covers
+ double protectedArea;
+ double unprotectedArea;
+ if (lcType.isProtectable()) {
+ protectedArea = totalArea * adjProtectedFraction;
+ unprotectedArea = totalArea * (1 - adjProtectedFraction);
+ } else {
+ protectedArea = 0;
+ unprotectedArea = totalArea;
}
- }
- }
-
- public void addLandUseTiles(Map<LandCoverType, Map<Integer, Double>> ageDist) {
-
- // Initialise array lists
- for (LandCoverType lcType : LandCoverType.getConvertibleTypes()) {
- ArrayList<LandUseTile> tiles = new ArrayList<LandUseTile>();
- landUseTiles.put(lcType, tiles);
- }
-
- // Add tiles
- for (Map.Entry<LandCoverType, Map<Integer, Double>> lcEntry : ageDist.entrySet()) {
- LandCoverType lcType = lcEntry.getKey();
- if (landCoverAreas.get(lcType) == 0)
- continue;
-
- double protectedArea = landCoverAreas.get(lcType) - unprotectedAreas.get(lcType);
- double unprotectedArea = unprotectedAreas.get(lcType);
- Map<Integer, Double> ageMap = lcEntry.getValue();
+
+ Map<Integer, Double> ageMap = landCover.getLandUseAgeDist().get(lcType);
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
+ 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);
}
@@ -136,7 +88,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
ageMap.put(entry.getKey(), entry.getValue() / totalProp);
}
}
-
+
// Assuming land use distributed evenly in each age group
for (Map.Entry<Integer, Double> ageEntry : ageMap.entrySet()) {
double prop = ageEntry.getValue() / 10;
@@ -146,80 +98,208 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
int minAge = (int) (ageGroup - 1) * 10 + 1;
int maxAge = minAge + 9;
for (int a = minAge; a <= maxAge; a++) {
- int yearEstablished = ModelConfig.BASE_YEAR - a;
-
- LandUseTile unprotectedTile = new LandUseTile(lcType, unprotectedArea * prop, yearEstablished, false);
- LandUseTile protectedTile = new LandUseTile(lcType, protectedArea * prop, yearEstablished, true);
-
- landUseTiles.get(lcType).add(unprotectedTile);
- landUseTiles.get(lcType).add(protectedTile);
+ landCoverAreas.get(lcType).setConvertibleArea(a, unprotectedArea * prop);
+ landCoverAreas.get(lcType).setProtectedArea(a, protectedArea * prop);
}
}
- }
+ }
+
+ // Check that areas have been allocated correctly
+ Map<LandCoverType, Double> initAreas = new HashMap<LandCoverType, Double>();
+ for (LandCoverType lc : LandCoverType.values()) {
+ initAreas.put(lc, landCover.getLandCoverArea(lc));
+ }
+ runAreaCheck(initAreas);
}
- public Map<LandCoverType, Set<Integer>> getTilesAllYearsEstablished() {
- Map<LandCoverType, Set<Integer>> reqData = new HashMap<LandCoverType, Set<Integer>>();
- for (Map.Entry<LandCoverType, ArrayList<LandUseTile>> entry : landUseTiles.entrySet()) {
+
+ public void runAreaCheck(Map<LandCoverType, Double> lcAreas) {
+ double threshold = 1e-10;
+ for (Map.Entry<LandCoverType, LandCoverTile> entry : landCoverAreas.entrySet()) {
LandCoverType lcType = entry.getKey();
- Set<Integer> yearList = new HashSet<Integer>(); // <yearEstablished>
- ArrayList<LandUseTile> tiles = entry.getValue();
- if (tiles.isEmpty())
- continue;
- for (LandUseTile t : tiles) {
- yearList.add(t.getYearEstablished());
+ LandCoverTile tiles = entry.getValue();
+ double totalArea = tiles.getTotalConvertibleArea();
+ if (totalArea - lcAreas.get(lcType) > threshold) {
+ LogWriter.printlnError("LandUseItem: unprotected area and tiles area different: " + lcType.getName() + " " + lcAreas.get(lcType) + ", " + totalArea);
}
- reqData.put(lcType, yearList);
}
- return reqData;
- }
-
-
- public Map<LandCoverType, ArrayList<LandUseTile>> getLandUseTiles() {
- return landUseTiles;
}
- public void setProtectedFract(double protFrac) {
- if (!Double.isNaN(ModelConfig.CONSTANT_PROTECTED_AREA_RATE))
- protFrac=ModelConfig.CONSTANT_PROTECTED_AREA_RATE;
+ public void setUnavailableArea(double unavailF) {
+ double barrenAndUrban = getLandCoverArea(LandCoverType.BARREN) + getLandCoverArea(LandCoverType.URBAN);
+ double unavailableArea = unavailF * getTotalLandCoverArea();
+ double shortfall = Math.max(0.0, unavailableArea - barrenAndUrban);
+ if (shortfall == 0.0)
+ return;
+
+ double convertibleNatural = getConvertibleLandCoverArea(LandCoverType.NATURAL);
+ double protectedNatural = getProtectedLandCoverArea(LandCoverType.NATURAL);
+ double convertibleTimberF = getConvertibleLandCoverArea(LandCoverType.TIMBER_FOREST);
+ /*
+ Map<LandCoverType, Double> foo = getLandCoverAreaMap();
+ for (Map.Entry<LandCoverType, Double> entry : foo.entrySet()) {
+ LogWriter.println(entry.getKey().getName() + "," + entry.getValue());
+ }
+ */
+ double unavailFromConvertibleNatural = Math.min(shortfall, convertibleNatural);
+ shortfall -= unavailFromConvertibleNatural;
+ if (unavailFromConvertibleNatural > 0) {
+ landCoverAreas.get(LandCoverType.NATURAL).addUnavailableArea(unavailFromConvertibleNatural);
+ landCoverAreas.get(LandCoverType.NATURAL).removeConvertibleArea(unavailFromConvertibleNatural);
+ }
- double totalArea = getTotalLandCoverArea();
- double urbanArea = getLandCoverArea(LandCoverType.URBAN);
- double barrenArea = getLandCoverArea(LandCoverType.BARREN);
+ double unavailFromProtectedNatural = Math.min(shortfall, protectedNatural);
+ shortfall -= unavailFromProtectedNatural;
+ if (unavailFromProtectedNatural > 0) {
+ landCoverAreas.get(LandCoverType.NATURAL).addUnavailableArea(unavailFromProtectedNatural);
+ landCoverAreas.get(LandCoverType.NATURAL).removeProtectedArea(unavailFromProtectedNatural);
+ }
- this.protectedArea = Math.min(totalArea-urbanArea-barrenArea, protFrac*(totalArea-urbanArea));
+ double unavailFromTimberF = Math.min(shortfall, convertibleTimberF);
+ shortfall -= unavailFromTimberF;
+ if (unavailFromTimberF > 0) {
+ landCoverAreas.get(LandCoverType.TIMBER_FOREST).addUnavailableArea(unavailFromTimberF);
+ landCoverAreas.get(LandCoverType.TIMBER_FOREST).removeConvertibleArea(unavailFromTimberF);
+ }
+ if (shortfall > 0) {
+ LogWriter.printlnWarning("LandUseItem.setUnavailableArea insufficient area to make unavailable");
+ }
+ }
+
+ public Map<LandCoverType, LandCoverTile> getLandUseTiles() {
+ return landCoverAreas;
}
+
+ public void moveTilesArea(LandCoverType toType, LandCoverType fromType, double area) {
+ // Note if toType == fromType, area will be removed from fromType and added back in with age 0
+ if (area == 0)
+ return;
- public void setProtectedArea(double protectedArea) {
- this.protectedArea = protectedArea;
+ LandCoverTile toTiles = landCoverAreas.get(toType);
+ LandCoverTile fromTiles = landCoverAreas.get(fromType);
+
+ if (fromTiles == null) {
+ LogWriter.printlnError("LandUseItem: cannot move area between tiles as source is null");
+ return;
+ }
+
+ fromTiles.removeConvertibleArea(area);
+ toTiles.addConvertibleArea(area);
+ }
+
+ public void setProtectedFraction(double protFrac) {
+ if (!ModelConfig.PROTECTED_AREAS_ENABLED) {
+ this.protectedFraction = 0;
+ return;
+ }
+ if (!Double.isNaN(ModelConfig.CONSTANT_PROTECTED_AREA_RATE)) {
+ this.protectedFraction = ModelConfig.CONSTANT_PROTECTED_AREA_RATE; // TODO is this a fraction or area?
+ } else {
+ this.protectedFraction = Math.max(protFrac, ModelConfig.MIN_NATURAL_RATE);
+ }
}
- public double getProtectedArea() {
+ public double getProtectedFraction() {
+ return protectedFraction;
+ }
+
+ public double getTotalProtectedArea() {
if(ModelConfig.PROTECTED_AREAS_ENABLED)
- return protectedArea;
+ return protectedFraction * getTotalLandCoverArea();
else
return 0.0;
}
- public void setSuitableArea(double suitable) {
- this.suitableArea = suitable;
+ public double getProtectedLandCoverArea(LandCoverType lcType) {
+ return landCoverAreas.get(lcType).getTotalProtectedArea();
}
-
+
public double getSuitableArea() {
+ double suitableArea = 0;
+ for (LandCoverType lcType : LandCoverType.getConvertibleTypes()) {
+ suitableArea += getConvertibleLandCoverArea(lcType);
+ }
return suitableArea;
}
+ /** move areas from one land cover to another, return any residual not possible */
+ public double moveAreas(LandCoverType toType, LandCoverType fromType, double area) {
+
+ double prevFrom = getConvertibleLandCoverArea(fromType);
+ double areaMoved;
+
+ areaMoved = Math.min(area, prevFrom);
+
+ if (Math.abs(area - areaMoved) / area > 1e-6) {
+ throw new RuntimeException("Move unprotected area shortfall. Should not happen due to GAMS constraint.");
+ }
+
+ moveTilesArea(toType, fromType, areaMoved);
+
+ return area - areaMoved;
+ }
- public void setUnavailableArea(double unavailF) {
- // subtract BARREN from unavailableArea, set to 0 if negative
- unavailableArea = Math.max(unavailF*getTotalLandCoverArea() - getLandCoverArea(LandCoverType.BARREN), 0.0);
+ public double getLandCoverArea(LandCoverType c) {
+ Double d = landCoverAreas.get(c).getTotalArea();
+ return d == null ? 0.0 : d;
}
- public double getUnavailableArea() {
- return unavailableArea;
+ public double getLandCoverFract(LandCoverType c) {
+ double totalArea = getTotalLandCoverArea();
+ return totalArea==0 ? 0.0 : getLandCoverArea(c) / totalArea;
}
+
+ public void setConvertibleLandCoverArea(LandCoverType c, double d) {
+ if (Double.isNaN(d) || Double.isInfinite(d))
+ throw new RuntimeException("AreasItem for " + c + " is " + d);
+ double landCover = (d < 0.0) ? 0.0 : d;
+
+ landCoverAreas.get(c).removeAllArea();
+ landCoverAreas.get(c).addConvertibleArea(landCover);;
+ }
+
+ public double getTotalLandCoverArea() {
+ double d = 0;
+ for (LandCoverType l : LandCoverType.values())
+ d += getLandCoverArea(l);
+
+ return d;
+ }
+
+ public Map<LandCoverType, Double> getLandCoverAreaMap() {
+ Map<LandCoverType, Double> areaMap = new HashMap<LandCoverType, Double>();
+ for (LandCoverType lcType : LandCoverType.values()) {
+ areaMap.put(lcType, getLandCoverArea(lcType));
+ }
+ return areaMap;
+ }
+
+ public double getTotalConvertibleNaturalArea() {
+ double unprotectedNatural = 0;
+ for (LandCoverType landType : LandCoverType.getNaturalTypes()) {
+ unprotectedNatural += landCoverAreas.get(landType).getTotalConvertibleArea();
+ }
+ return unprotectedNatural;
+ }
+
+ public double getConvertibleLandCoverArea(LandCoverType c) {
+ return landCoverAreas.get(c).getTotalConvertibleArea();
+ }
+
+
+ public void updateProtectedFraction(int year) {
+ if (!ModelConfig.FORCE_PROTECTED_AREAS)
+ return;
+
+ if (year >= ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR && year < ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR)
+ protectedFraction = 1.0 / (ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR - year);
+
+ if (year >= ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR)
+ protectedFraction = 1.0;
+ }
+
public Intensity getIntensity(CropType crop) {
return intensityMap.get(crop);
}
@@ -321,56 +401,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
public static double getIrrigationTotal(Collection<? extends LandUseItem> items, Collection<CropType> crops) {
return getIrrigationTotal(items, crops.toArray(new CropType[crops.size()]));
}
-
- /** move areas from one land cover to another, return any residual not possible */
- public double moveAreas(LandCoverType toType, LandCoverType fromType, double changeReq) {
-
- if (fromType.equals(toType))
- return 0;
-
- double prevTo = getLandCoverArea(toType);
- double prevFrom = getLandCoverArea(fromType);
-
- double actuallyChanged = Math.min(prevFrom, changeReq);
-
- /*
- if (changeReq - actuallyChanged > 0) {
- throw new RuntimeException("Move area shortfall. Should not happen due to GAMS constraint.");
- }*/
-
- setLandCoverArea(toType, prevTo + actuallyChanged);
- setLandCoverArea(fromType, prevFrom - actuallyChanged);
-
- return changeReq - actuallyChanged;
- }
-
- public double moveUnprotectedArea(LandCoverType toType, LandCoverType fromType, double area, Timestep timestep) {
-
- double prevTo = getUnprotectedLandCoverArea(toType);
- double prevFrom = getUnprotectedLandCoverArea(fromType);
- double areaMoved;
-
- areaMoved = Math.min(area, prevFrom);
-
- if (!fromType.equals(toType)) {
- unprotectedAreas.put(toType, prevTo + areaMoved);
- unprotectedAreas.put(fromType, prevFrom - areaMoved);
- }
-
- /*
- if (area - areaMoved > 0) {
- throw new RuntimeException("Move unprotected area shortfall. Should not happen due to GAMS constraint.");
- }
- */
-
- double foo = landUseTiles.get(fromType).stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
- double foo2 = landUseTiles.get(fromType).stream().mapToDouble(o -> o.getArea()).sum();
-
- moveTilesArea(toType, fromType, areaMoved, timestep);
-
-
- return area - areaMoved;
- }
public double getCropArea(CropType c) {
Double d = getLandCoverArea(LandCoverType.CROPLAND) * getCropFraction(c);
@@ -390,38 +420,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return totalFract;
}
-
public void setCropFraction(CropType c, double areaFract) {
cropFractions.put(c, areaFract);
}
-
- public double getLandCoverArea(LandCoverType c) {
- Double d = landCoverAreas.get(c);
- return d == null ? 0.0 : d;
- }
-
- public double getLandCoverFract(LandCoverType c) {
- double totalArea = getTotalLandCoverArea();
- return totalArea==0 ? 0.0 : getLandCoverArea(c) / totalArea;
- }
-
- public void setLandCoverArea(LandCoverType c, double d) {
- if (Double.isNaN(d) || Double.isInfinite(d))
- throw new RuntimeException("AreasItem for " + c + " is " + d);
- double landCover = (d < 0.0) ? 0.0 : d;
-
- landCoverAreas.put(c, landCover);
- }
-
- public double getTotalLandCoverArea() {
- double d = 0;
- for (LandCoverType l : LandCoverType.values())
- d += getLandCoverArea(l);
-
- return d;
- }
-
public static double getAbandonedPasture(Collection<? extends LandUseItem> landUses) {
double total = 0;
for (LandUseItem a : landUses) {
@@ -445,81 +447,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
}
return total;
}
-
- public double getTotalNatural() {
- double totalNatural = getLandCoverArea(LandCoverType.NATURAL);
- return totalNatural;
- }
-
- public double getTotalUnprotectedNatural() {
- double unprotectedNatural = 0;
- for (LandCoverType landType : LandCoverType.values()) {
- if (landType.isNatural()) {
- unprotectedNatural += unprotectedAreas.get(landType);
- }
- }
- return unprotectedNatural;
- }
-
- public double getUnprotectedLandCoverArea(LandCoverType landType) {
- Double d = unprotectedAreas.get(landType);
- return d == null ? 0.0 : d;
- }
-
- public void setUnprotectedLandCoverArea(LandCoverType landType, double d) {
- unprotectedAreas.put(landType, d);
-
- }
-
- private double getProtectedandUnavailable() {
- double totalSuitableLandCover = getTotalLandCoverArea() - getLandCoverArea(LandCoverType.BARREN) - getLandCoverArea(LandCoverType.URBAN);
- if (totalSuitableLandCover < 0)
- throw new RuntimeException("totalSuitableLandCover < 0");
- //return 0;
- // Protected area should be at least MIN_NATURAL_RATE
- double desiredProtectedArea = Math.max(getProtectedArea(), totalSuitableLandCover * ModelConfig.MIN_NATURAL_RATE);
- // if desired protected is more than suitable then protect all suitable
- if (desiredProtectedArea + getUnavailableArea() > totalSuitableLandCover) {
- LogWriter.printlnWarning("Desired protected area exceeds total suitable area");
- }
- return Math.min(totalSuitableLandCover, desiredProtectedArea + getUnavailableArea());
- }
-
- /** averages protected areas across land cover types */
- private void updateUnprotectedAreas() {
- double desiredProtected = getProtectedandUnavailable();
- double totalNatural = getTotalNatural();
- double unprotectedNat = Math.max(0, totalNatural - desiredProtected); // if desiredProtected more than total Natural then protect all
-
- double unprotectedFract = (totalNatural > 0) ? unprotectedNat/totalNatural : 0;
-
- for (LandCoverType landType : LandCoverType.values()) {
- if (landType.isProtectable()) {
- unprotectedAreas.put(landType, getLandCoverArea(landType)*unprotectedFract);
- }
- else unprotectedAreas.put(landType, getLandCoverArea(landType));
- }
- }
-
- public void updateSuitableArea(int year) {
- updateUnprotectedAreas(); //update unprotected areas
- double currentAgri = getLandCoverArea(LandCoverType.PASTURE) + getLandCoverArea(LandCoverType.CROPLAND);
- double totalNatural = getTotalNatural();
- double natAvailForAgri = totalNatural - getProtectedandUnavailable(); // could be negative, i.e. excess agricultural area already used
- //this forces protected areas, min_natural_rate and slope constraints to be obeyed
- if (ModelConfig.FORCE_PROTECTED_AREAS && natAvailForAgri < 0 && year >= ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR) {
- double proportion = 1.0;
-
- if (year < ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR)
- proportion = 1.0 / (ModelConfig.FORCE_PROTECTED_AREAS_END_YEAR - year);
-
- suitableArea = Math.max(0, currentAgri + natAvailForAgri * proportion); // netNatAvailForAgri is negative, but suitable area < 0 is not sensible (seems to happen with high barren areas)
- }
- else //not honouring protected areas if agriculture > natural that should be protected
- suitableArea = currentAgri + getTotalUnprotectedNatural();
- }
-
public CropToDouble getCropChanges(LandUseItem prevAreaAggItem) {
CropToDouble changes = new CropToDouble();
@@ -531,14 +459,16 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return changes;
}
- public void doForestRotation(Timestep timestep, double rotationAge) {
- ArrayList<LandUseTile> tiles = landUseTiles.get(LandCoverType.TIMBER_FOREST);
- if (tiles == null)
- return;
- for (LandUseTile t : tiles) {
- if (t.getAge(timestep) < rotationAge)
- continue;
- t.resetAge(timestep); // Set stand age to 0
+ public void doForestRotation(int rotationAge) {
+ LandCoverTile tiles = landCoverAreas.get(LandCoverType.TIMBER_FOREST);
+ for (int a=rotationAge; a<=tiles.getMaxAgeBin(); a++) {
+ tiles.resetAreaForAge(a);
+ }
+ }
+
+ public void incrementTilesAge() {
+ for (LandCoverTile tile : landCoverAreas.values()) {
+ tile.incrementAge();
}
}
@@ -549,10 +479,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
@Override
public void interpolateAll(LandUseItem fromItem, LandUseItem toItem, double factor) {
cropFractions = new HashMap<CropType, Double>();
- landCoverAreas = new HashMap<LandCoverType, Double>();
+ landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
- Double fromCropCover = fromItem.landCoverAreas.get(LandCoverType.CROPLAND);
- Double toCropCover = toItem.landCoverAreas.get(LandCoverType.CROPLAND);
+ Double fromCropCover = fromItem.landCoverAreas.get(LandCoverType.CROPLAND).getTotalArea();
+ Double toCropCover = toItem.landCoverAreas.get(LandCoverType.CROPLAND).getTotalArea();
if (!isZeroOrNull(fromCropCover) && isZeroOrNull(toCropCover)) { // if start with crop but end with none, take starting crop fractions
cropFractions.putAll(fromItem.cropFractions);
@@ -583,10 +513,19 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
}
for (LandCoverType landCover : LandCoverType.values()) {
- Double from = fromItem.landCoverAreas.get(landCover);
- Double to = toItem.landCoverAreas.get(landCover);
+ Double from = fromItem.landCoverAreas.get(landCover).getTotalConvertibleArea();
+ Double to = toItem.landCoverAreas.get(landCover).getTotalConvertibleArea();
+ Double d = Interpolator.interpolate(from, to, factor);
+ LandCoverTile tile = landCoverAreas.computeIfAbsent(landCover, k -> new LandCoverTile());
+ tile.addConvertibleArea(d);
+ }
+
+ for (LandCoverType landCover : LandCoverType.values()) {
+ Double from = fromItem.landCoverAreas.get(landCover).getTotalProtectedArea();
+ Double to = toItem.landCoverAreas.get(landCover).getTotalProtectedArea();
Double d = Interpolator.interpolate(from, to, factor);
- landCoverAreas.put(landCover, d);
+ LandCoverTile tile = landCoverAreas.computeIfAbsent(landCover, k -> new LandCoverTile());
+ tile.addProtectedArea(d);
}
}
@@ -624,6 +563,12 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
@Override
public String toString() {
- return "LandUseItem: [landCoverAreas=" + landCoverAreas + ", protectedArea=" + protectedArea + ", unavailableArea=" + unavailableArea + "]";
+ 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.get(lcType).getTotalConvertibleArea());
+ protectedAreas.put(lcType, landCoverAreas.get(lcType).getTotalProtectedArea());
+ }
+ return "LandUseItem: [landCoverAreas=" + convertibleAreas + ", protectedArea=" + protectedAreas + "]";
}
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/landuse/LandUseReader.java b/src/ac/ed/lurg/landuse/LandUseReader.java
index a857eff3..bc896330 100644
--- a/src/ac/ed/lurg/landuse/LandUseReader.java
+++ b/src/ac/ed/lurg/landuse/LandUseReader.java
@@ -21,10 +21,10 @@ public class LandUseReader extends AbstractTabularRasterReader<LandUseItem> {
protected void setData(RasterKey key, LandUseItem luData, Map<String, Double> rowValues) {
for (LandCoverType cover : LandCoverType.values()) {
- luData.setLandCoverArea(cover, getValueForCol(rowValues, cover.getName()));
+ luData.setConvertibleLandCoverArea(cover, getValueForCol(rowValues, cover.getName()));
}
- luData.setProtectedArea(getValueForCol(rowValues, "protected"));
+ luData.setProtectedFraction(getValueForCol(rowValues, "protected") / luData.getTotalLandCoverArea());
for (CropType crop : CropType.getAllItems()) {
double cropFract = getValueForCol(rowValues, crop.getGamsName() + "_A");
diff --git a/src/ac/ed/lurg/landuse/LandUseTile.java b/src/ac/ed/lurg/landuse/LandUseTile.java
deleted file mode 100644
index 7af6fb13..00000000
--- a/src/ac/ed/lurg/landuse/LandUseTile.java
+++ /dev/null
@@ -1,64 +0,0 @@
-package ac.ed.lurg.landuse;
-
-import java.io.Serializable;
-
-import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.Timestep;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.LogWriter;
-import ac.sac.raster.InterpolatingRasterItem;
-
-public class LandUseTile implements Serializable, InterpolatingRasterItem<LandUseTile> {
- private static final long serialVersionUID = 3208516181615267472L;
-
- private LandCoverType lcType;
- private double area; // Mha
- private int yearEstablished;
- private boolean isProtected;
-
- public LandUseTile(LandCoverType lcType, double area, int yearEstablished, boolean isProtected) {
- this.lcType = lcType;
- this.area = area;
- this.yearEstablished = yearEstablished;
- this.isProtected = isProtected;
- }
-
- public LandCoverType getLcType() {
- return lcType;
- }
-
- public double getArea() {
- return area;
- }
-
- public int getYearEstablished() {
- return yearEstablished;
- }
-
- public int getAge(Timestep timestep) {
- return timestep.getYear() - yearEstablished;
- }
-
- public boolean isProtected() {
- return isProtected;
- }
-
- public void resetAge(Timestep timestep) {
- yearEstablished = timestep.getYear();
- }
-
- public void removeArea(double a) {
- if (a > area) {
- LogWriter.printlnError("LandUseTile.removeArea: attempting to remove too much area: " + a + ", available: " + area);
- area -= Math.min(area, a);
- } else {
- area -= a;
- }
- }
-
- @Override
- public void interpolateAll(LandUseTile fromItem, LandUseTile toItem, double factor) {
- // TODO Auto-generated method stub
-
- }
-}
diff --git a/src/ac/ed/lurg/landuse/ProtectedAreasReader.java b/src/ac/ed/lurg/landuse/ProtectedAreasReader.java
index dcf8ca6a..5a8b1e4d 100644
--- a/src/ac/ed/lurg/landuse/ProtectedAreasReader.java
+++ b/src/ac/ed/lurg/landuse/ProtectedAreasReader.java
@@ -13,7 +13,7 @@ public class ProtectedAreasReader extends AbstractRasterReader<LandUseItem> {
public void setData(LandUseItem lcData, String token) {
if (!"nan".equals(token)) {
double protFrac = Double.parseDouble(token);
- lcData.setProtectedFract(protFrac);
+ lcData.setProtectedFraction(protFrac);
}
}
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/output/LandUseOutputer.java b/src/ac/ed/lurg/output/LandUseOutputer.java
index 3d341eb2..82eed7ea 100644
--- a/src/ac/ed/lurg/output/LandUseOutputer.java
+++ b/src/ac/ed/lurg/output/LandUseOutputer.java
@@ -58,8 +58,8 @@ public class LandUseOutputer extends AbstractLandUseOutputer {
sbData.append(String.format(" %.8f", item.getTotalLandCoverArea()));
sbData.append(String.format(" %.8f", item.getSuitableArea()));
- sbData.append(String.format(" %.8f", item.getProtectedArea()));
- sbData.append(String.format(" %.8f", item.getProtectedArea()/item.getTotalLandCoverArea()));
+ sbData.append(String.format(" %.8f", item.getTotalProtectedArea()));
+ sbData.append(String.format(" %.8f", item.getTotalProtectedArea()/item.getTotalLandCoverArea()));
for (LandCoverType cover : LandCoverType.values()) {
sbData.append(String.format(" %.8f", item.getLandCoverArea(cover)));
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index 311704e9..7398c8a2 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -39,6 +39,10 @@ public enum LandCoverType {
return isProtectable;
}
+ public boolean isConvertible() {
+ return isConvertible;
+ }
+
public boolean isNatural() {
return isNatural;
}
--
GitLab