Merging changes from master

(Bart's calibration and land abandonment work)

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);