From 331f3aa424c60a93ee64577bb674b81057d4e227 Mon Sep 17 00:00:00 2001
From: Bart Arendarczyk <s1924442@ed.ac.uk>
Date: Wed, 12 Jun 2024 18:40:45 +0100
Subject: [PATCH] Clustering confined to FPUs.
---
src/ac/ed/lurg/ModelConfig.java | 8 +-
src/ac/ed/lurg/country/CountryAgent.java | 73 +++++++++++--------
.../ed/lurg/country/CountryAgentManager.java | 3 +-
src/ac/ed/lurg/landuse/FPUManager.java | 9 +++
src/ac/ed/lurg/landuse/Fpu.java | 4 +-
.../ed/lurg/landuse/IrrigationRasterSet.java | 4 +
src/ac/ed/lurg/yield/YieldClusterPoint.java | 29 ++++----
7 files changed, 79 insertions(+), 51 deletions(-)
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index c65c609f..ce1e7c6c 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -284,7 +284,7 @@ public class ModelConfig {
public static final double PROD_COST_MAIZE = getDoubleProperty("PROD_COST_MAIZE", 0.3);
public static final double PROD_COST_RICE = getDoubleProperty("PROD_COST_RICE", 0.8);
public static final double PROD_COST_OILCROPS_NFIX = getDoubleProperty("PROD_COST_OILCROPS_NFIX", 0.4);
- public static final double PROD_COST_OILCROPS_OTHER = getDoubleProperty("PROD_COST_OILCROPS_OTHER", 0.8);
+ public static final double PROD_COST_OILCROPS_OTHER = getDoubleProperty("PROD_COST_OILCROPS_OTHER", 0.25);
public static final double PROD_COST_PULSES = getDoubleProperty("PROD_COST_PULSES", 0.6);
public static final double PROD_COST_STARCHYROOTS = getDoubleProperty("PROD_COST_STARCHYROOTS", 3.0);
public static final double PROD_COST_MONOGASTRICS = getDoubleProperty("PROD_COST_MONOGASTRICS", 0.05);
@@ -515,7 +515,7 @@ public class ModelConfig {
public static final int LPJG_TIMESTEP_SIZE = getIntProperty("LPJG_TIMESTEP_SIZE", 5);
public static final int LPJ_YEAR_OFFSET = getIntProperty("LPJ_YEAR_OFFSET", -1);
- public static final int NUM_YIELD_CLUSTERS = getIntProperty("NUM_YIELD_CLUSTERS", 2000);
+ public static final int NUM_YIELD_CLUSTERS = getIntProperty("NUM_YIELD_CLUSTERS", 50);
public static final long RANDOM_SEED = getIntProperty("RANDOM_SEED", 1974329); // any number will do
public static final String CHECKPOINT_YEARS = getProperty("CHECKPOINT_YEARS", null); // 2020,2050
@@ -566,8 +566,8 @@ public class ModelConfig {
public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 1200.0); // million m3
public static final double WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 0.3); // m3 to tC-eq p.16 [https://doi.org/10.5194/gmd-13-5425-2020]
public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_MANAGEMENT_COST", 2.0) : 0.0; // establishment, management etc. $1000/ha
- public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.1); //$1000/m3
- public static final double INIT_WOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_WOOD_PRICE", 0.04) : 0.0; // $1000/m3
+ public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.2); //$1000/m3
+ public static final double INIT_WOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_WOOD_PRICE", 0.03) : 0.0; // $1000/m3
public static final int CARBON_WOOD_MAX_TIME = getIntProperty("CARBON_WOOD_MAX_TIME", 160); // upper data limit, years
// Carbon
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 8df574a9..1fc62f44 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -62,63 +62,76 @@ public class CountryAgent extends AbstractCountryAgent {
}
private RasterSet<IntegerRasterItem> calcYieldClusters(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces,
- RasterSet<WoodYieldItem> woodYieldData, RasterSet<SolarPotentialItem> solarData) {
+ RasterSet<WoodYieldItem> woodYieldData, RasterSet<SolarPotentialItem> solarData,
+ RasterSet<Fpu> fpuRasterSet) {
- LogWriter.println("calcYieldClusters: for " + ModelConfig.NUM_YIELD_CLUSTERS + " clusters", 3);
-
- // create collection of ClusteringPoints from countryYieldSurfaces, these have the RasterKey and data for yield (or access to them)
- HashSet<YieldClusterPoint> clusteringPoints = new HashSet<YieldClusterPoint>();
- for (Entry<RasterKey, YieldResponsesItem> entry : countryYieldSurfaces.entrySet()) {
- YieldResponsesItem yieldresp = entry.getValue();
- if (!yieldresp.isMissing()) {
- RasterKey key = entry.getKey();
+ LogWriter.println("calcYieldClusters: for " + ModelConfig.NUM_YIELD_CLUSTERS + " clusters", 3);
+
+ Map<Fpu, Set<RasterKey>> fpuMap = new HashMap<>();
+ for (Map.Entry<RasterKey, Fpu> entry : fpuRasterSet.entrySet()) {
+ Set<RasterKey> keySet = fpuMap.computeIfAbsent(entry.getValue(), k -> new HashSet<>());
+ keySet.add(entry.getKey());
+ }
+
+ RasterSet<IntegerRasterItem> mapping = new RasterSet<IntegerRasterItem>(countryYieldSurfaces.getHeaderDetails());
+
+ int clusterId = 1;
+ // Calculate clusters within each FPU
+ for (Fpu fpu : fpuMap.keySet()) {
+ Set<RasterKey> keySet = fpuMap.get(fpu);
+ HashSet<YieldClusterPoint> clusteringPoints = new HashSet<YieldClusterPoint>();
+
+ // create clustering points
+ for (RasterKey key : keySet) {
+ if (countryYieldSurfaces.get(key).isMissing()) {
+ continue; // missing yields, skipping. Later we put all cells with missing yield into one cluster
+ }
+
+ YieldResponsesItem yieldresp = countryYieldSurfaces.get(key);
IrrigationItem irrigItem = irrigData.get(key);
-
LandUseItem luItem = previousGamsRasterOutput.getLandUses().get(key);
double totalLand = luItem.getTotalLandArea()-luItem.getLandCoverArea(LandCoverType.URBAN);
double protectedAreaFrac = (totalLand <= 0) ? 0.0 : luItem.getLandCoverArea(LandProtectionType.PROTECTED) / totalLand;
-
WoodYieldItem wyItem = woodYieldData.get(key);
SolarPotentialItem solarItem = solarData.get(key);
-
+
clusteringPoints.add(new YieldClusterPoint(key, yieldresp, irrigItem, protectedAreaFrac, wyItem, solarItem));
}
- }
- // do the clustering
- KMeans<String, YieldClusterPoint> kmeans = new KMeans<String, YieldClusterPoint>(clusteringPoints, ModelConfig.NUM_YIELD_CLUSTERS);
- kmeans.calculateClusters(100, 0.1);
- kmeans.printClusters();
+ // do the clustering
+ KMeans<String, YieldClusterPoint> kmeans = new KMeans<String, YieldClusterPoint>(clusteringPoints, ModelConfig.NUM_YIELD_CLUSTERS);
+ kmeans.calculateClusters(100, 0.1);
- // reformat output
- List<Cluster<String, YieldClusterPoint>> yieldClusters = kmeans.getClusters();
- RasterSet<IntegerRasterItem> mapping = new RasterSet<IntegerRasterItem>(countryYieldSurfaces.getHeaderDetails());
-
- int id = 1;
- for (Cluster<String, YieldClusterPoint> c : yieldClusters) {
- if (!c.getPoints().isEmpty()) {
- for (YieldClusterPoint p : c.getPoints()) {
- mapping.put(p.getRasterKey(), new IntegerRasterItem(id));
+ // reformat output
+ List<Cluster<String, YieldClusterPoint>> yieldClusters = kmeans.getClusters();
+
+ for (Cluster<String, YieldClusterPoint> c : yieldClusters) {
+ if (!c.getPoints().isEmpty()) {
+ for (YieldClusterPoint p : c.getPoints()) {
+ mapping.put(p.getRasterKey(), new IntegerRasterItem(clusterId));
+ }
+ clusterId++;
}
- id++;
}
}
// Add cells with missing yields to one cluster
for (RasterKey key : countryYieldSurfaces.keySet()) {
if (countryYieldSurfaces.get(key).isMissing()) {
- IntegerRasterItem item = new IntegerRasterItem(id);
+ IntegerRasterItem item = new IntegerRasterItem(clusterId);
mapping.put(key, item);
}
}
+ LogWriter.println(String.format("%s calculated %d clusters", country, clusterId));
+
return mapping;
}
public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData,
Map<CropType, GlobalPrice> worldPrices, RasterSet<CarbonFluxItem> carbonFluxData,
RasterSet<WoodYieldItem> woodYieldData, RasterSet<SolarPotentialItem> solarPotentialData,
- GlobalPrice carbonPrice, GlobalPrice woodPrice) {
+ GlobalPrice carbonPrice, GlobalPrice woodPrice, RasterSet<Fpu> fpuRasterSet) {
// project demand
calculateCountryPricesAndDemand(worldPrices, carbonPrice, woodPrice, true);
@@ -140,7 +153,7 @@ public class CountryAgent extends AbstractCountryAgent {
}
else {
if (yieldClusters==null) {
- yieldClusters = calcYieldClusters(irrigData, countryYieldSurfaces, woodYieldData, solarPotentialData); // this should only be on the first timestep
+ yieldClusters = calcYieldClusters(irrigData, countryYieldSurfaces, woodYieldData, solarPotentialData, fpuRasterSet); // this should only be on the first timestep
}
// optimize areas and intensity
GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, woodYieldData, carbonFluxData, solarPotentialData);
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index 2aa9af4a..47848214 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -145,6 +145,7 @@ public class CountryAgentManager {
RasterSet<IrrigationItem> irrigData = currentIrrigationData.createSubsetForKeys(countryKeys);
RasterSet<CarbonFluxItem> carbonFluxData = currentCarbonFluxData.createSubsetForKeys(countryKeys);
RasterSet<WoodYieldItem> woodYieldData = currentWoodYieldData.createSubsetForKeys(countryKeys);
+ RasterSet<Fpu> fpuRasterSet = currentIrrigationData.getFpuRasterSet().createSubsetForKeys(countryKeys);
// Need to instantiate Runnable object for multithreading
execService.execute(new Runnable() {
@@ -153,7 +154,7 @@ public class CountryAgentManager {
try {
ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), carbonFluxData, woodYieldData,
- solarPotentialData, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice());
+ solarPotentialData, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice(), fpuRasterSet);
// update global rasters
globalLandUseRaster.putAll(ca.getLandUses());
diff --git a/src/ac/ed/lurg/landuse/FPUManager.java b/src/ac/ed/lurg/landuse/FPUManager.java
index 12620562..3838c027 100644
--- a/src/ac/ed/lurg/landuse/FPUManager.java
+++ b/src/ac/ed/lurg/landuse/FPUManager.java
@@ -25,11 +25,14 @@ public class FPUManager {
private final StringTabularReader fpuGroups;
private final Map<WaterBasin, List<RasterKey>> basinToRasterKeysMap = new HashMap<>();
private final Map<WaterBasin, Map<Integer, Double>> basinToOtherUses = new HashMap<>();
+ private final RasterSet<Fpu> fpuMap; // for clustering
public FPUManager(RasterHeaderDetails desiredProjection) {
fpuGroups = new StringTabularReader(",", new String[]{"fpu", "basin"});
fpuGroups.read(ModelConfig.FPU_GROUPING_FILE);
+ fpuMap = new RasterSet<>(desiredProjection);
+
// Handle rasterkeys
Map<Fpu, List<RasterKey>> fpuToRasterKeys = getMapFpuToKeys(desiredProjection);
for (Map.Entry<Fpu, List<RasterKey>> entry : fpuToRasterKeys.entrySet()) {
@@ -90,6 +93,8 @@ public class FPUManager {
List<RasterKey> rasterKeys = invertedBoundaries.computeIfAbsent(fpu, k -> new ArrayList<>());
rasterKeys.add(entry.getKey());
+
+ fpuMap.put(entry.getKey(), fpu);
}
return invertedBoundaries;
@@ -147,4 +152,8 @@ public class FPUManager {
public Collection<WaterBasin> getWaterBasins() {
return basinToRasterKeysMap.keySet();
}
+
+ public RasterSet<Fpu> getFpuRasterSet() {
+ return fpuMap;
+ }
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/landuse/Fpu.java b/src/ac/ed/lurg/landuse/Fpu.java
index 5164b04a..0c3b031a 100644
--- a/src/ac/ed/lurg/landuse/Fpu.java
+++ b/src/ac/ed/lurg/landuse/Fpu.java
@@ -1,8 +1,10 @@
package ac.ed.lurg.landuse;
+import ac.sac.raster.RasterItem;
+
/** This is currently really just a marker class to aid type safety, using Integer directly would have worked, but been less clear*/
-public class Fpu {
+public class Fpu implements RasterItem {
private Integer id;
Fpu(Integer id) {
diff --git a/src/ac/ed/lurg/landuse/IrrigationRasterSet.java b/src/ac/ed/lurg/landuse/IrrigationRasterSet.java
index 621a579e..fe896bb4 100644
--- a/src/ac/ed/lurg/landuse/IrrigationRasterSet.java
+++ b/src/ac/ed/lurg/landuse/IrrigationRasterSet.java
@@ -88,4 +88,8 @@ public class IrrigationRasterSet extends RasterSet<IrrigationItem> {
}
return lpjIrrigConstraint;
}
+
+ public RasterSet<Fpu> getFpuRasterSet() {
+ return fpuManager.getFpuRasterSet();
+ }
}
diff --git a/src/ac/ed/lurg/yield/YieldClusterPoint.java b/src/ac/ed/lurg/yield/YieldClusterPoint.java
index c3560663..083df7a4 100644
--- a/src/ac/ed/lurg/yield/YieldClusterPoint.java
+++ b/src/ac/ed/lurg/yield/YieldClusterPoint.java
@@ -1,17 +1,16 @@
package ac.ed.lurg.yield;
-import java.util.Arrays;
-import java.util.Collection;
-
import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.solar.SolarPotentialItem;
import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.YieldType;
import ac.ed.lurg.utils.cluster.ClusteringPoint;
import ac.sac.raster.RasterKey;
+import java.util.Arrays;
+import java.util.Collection;
+
public class YieldClusterPoint implements ClusteringPoint<String> {
private final static String PASTURE = "pas";
@@ -45,17 +44,17 @@ public class YieldClusterPoint implements ClusteringPoint<String> {
// not sure if we be better to get a reference to the YieldResponsesItem, rather than caching these values?
// multiplying each variable by scaling factor so all are on a similar scale
- this.wheatMin = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.WHEAT) * 0.1;
- this.riceMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.RICE) * 0.1;
- this.maizeMin = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.MAIZE) * 0.1;
- this.pasture = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.PASTURE) * 0.1;
- this.wheatMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.WHEAT) * 0.1;
- this.maizeMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.MAIZE) * 0.1;
- this.irrigWaterAval = (irrigItem==null) ? 0.0 : irrigItem.getIrrigConstraint() * 0.1;
- this.protectedArea = protectedArea;
- this.woodMax = (woodYields==null) ? 0.0 : woodYields.getYieldResponse().getMaxYield() * 0.01;
- this.woodSlope = (woodYields==null) ? 0.0 : -woodYields.getYieldResponse().getSlope() * 10;
- this.solarPotential = solarPotentialItem == null ? 0.0 : solarPotentialItem.getPvPotential() * 0.001;
+ this.wheatMin = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.WHEAT);
+ this.riceMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.RICE);
+ this.maizeMin = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.MAIZE);
+ this.pasture = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.PASTURE);
+ this.wheatMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.WHEAT);
+ this.maizeMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.MAIZE);
+ this.irrigWaterAval = irrigItem == null ? 0.0 : irrigItem.getIrrigConstraint();
+ this.protectedArea = protectedArea * 10;
+ this.woodMax = woodYields == null ? 0.0 : woodYields.getYieldResponse().getMaxYield() * 0.1;
+ this.woodSlope = woodYields == null ? 0.0 : -woodYields.getYieldResponse().getSlope();
+ this.solarPotential = solarPotentialItem == null ? 0.0 : solarPotentialItem.getPvPotential() * 0.01;
}
public RasterKey getRasterKey() {
--
GitLab