From a902de69ed4f8e1830dc075fd95c9c3b1f24386f Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Mon, 30 Aug 2021 19:16:04 +0100
Subject: [PATCH] Changes to how wood yield data is read in; reduces memory
 usage.

---
 debug_config.properties                       |   2 +
 src/ac/ed/lurg/ModelConfig.java               |   9 +-
 src/ac/ed/lurg/ModelMain.java                 |  88 +++--------
 src/ac/ed/lurg/Timestep.java                  |  13 ++
 src/ac/ed/lurg/country/CountryAgent.java      |  66 +++-----
 .../country/gams/GamsLocationOptimiser.java   |  26 +++-
 .../ed/lurg/country/gams/GamsRasterInput.java |   8 +-
 .../country/gams/GamsRasterOptimiser.java     |  61 +++++---
 src/ac/ed/lurg/forestry/ForestManager.java    |   4 +-
 src/ac/ed/lurg/forestry/WoodYieldData.java    |  27 +++-
 src/ac/ed/lurg/forestry/WoodYieldItem.java    | 142 ++++++++++--------
 src/ac/ed/lurg/forestry/WoodYieldReader.java  |  94 ++++++++++--
 src/ac/ed/lurg/landuse/LandUseItem.java       |  89 +++++------
 src/ac/ed/lurg/landuse/LandUseTile.java       |   2 +-
 src/ac/ed/lurg/types/LandCoverType.java       |  12 ++
 15 files changed, 365 insertions(+), 278 deletions(-)

diff --git a/debug_config.properties b/debug_config.properties
index 86fd3726..ff1a0558 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -4,6 +4,8 @@ YIELD_DIR=C:/Users/Bart/Documents/PhD/LURG/rcp60
 
 OUTPUT_DIR = C:/Users/Bart/git/plumv2/output
 
+WOOD_AND_CARBON_DIR = C:/Users/Bart/Documents/PhD/Carbon and wood yields/forestry_dummy_4
+
 BASE_YEAR=2010
 START_TIMESTEP=0
 TIMESTEP_SIZE=1
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index bee3b092..1be1061b 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -248,12 +248,15 @@ public class ModelConfig {
 	
 	// Wood/carbon data
 	public static final String FOREST_DIR = SPATIAL_DATA_DIR + File.separator + "forestry";
+	public static final String WOOD_AND_CARBON_DIR = getProperty("WOOD_AND_CARBON_DIR");
 	//public static final String CARBON_FLUX_FILE = FOREST_DIR + File.separator + "carbon_flux_";
-	public static final String WOOD_YIELD_FILENAME = "wood_yield.out";
-	public static final String TIMBER_FOREST_GROWTH_FILENAME = "wood_yield.out";
-	public static final String CARBON_FOREST_GROWTH_FILENAME = "growth_carbon.out";
+	public static final String WOOD_YIELD_AGRI_FILENAME = "wood_yield_to_agri.out";
+	public static final String WOOD_YIELD_FRST_FILENAME = "wood_yield_to_frst.out";
+	public static final String CARBON_LUC_FILENAME = "carbon_luc.out";
+	public static final String CARBON_NEE_FILENAME = "carbon_nee.out";
 	public static final String NATURAL_FOREST_GROWTH_FILENAME = "growth_natural.out";
 	public static final String LAND_COVER_AGE_DIST_FILENAME = SPATIAL_DATA_DIR + File.separator + "land_cover_age_dist.txt";
+	public static final int WOOD_AND_CARBON_TIMESTEP_SIZE = getIntProperty("WOOD_AND_CARBON_TIMESTEP_SIZE", 20); // years
 
 	// Output
 	public static final String LAND_COVER_OUTPUT_FILE = OUTPUT_DIR + File.separator + "lc.txt";
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 6a1f870e..112f2cdd 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -85,6 +85,8 @@ public class ModelMain {
 	private IrrigationRasterSet currentIrrigationData;
 	private RasterSet<LandUseItem> globalLandUseRaster;
 	private RasterSet<IntegerRasterItem> clusterIdRaster;
+	WoodYieldReader woodYieldReader;
+	
 
 	public static void main(String[] args) {
 		ModelMain theModel = new ModelMain();
@@ -119,6 +121,8 @@ public class ModelMain {
 		globalLandUseRaster = new RasterSet<LandUseItem>(desiredProjection);
 		internationalMarket = new InternationalMarket();
 		
+		woodYieldReader = new WoodYieldReader(desiredProjection);
+		
 		createCountryAgents(compositeCountryManager.getAll());
 	}
 
@@ -136,10 +140,12 @@ 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
+		WoodYieldRasterSet woodYieldData = getWoodYieldData(timestep);
 
 		// When running half earth we can to alter protected areas data at a point in time
 		if(ModelConfig.HALFEARTH && ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR == timestep.getYear() && !ModelConfig.IS_CALIBRATION_RUN) {
@@ -161,7 +167,7 @@ public class ModelMain {
 		Map<LccKey, Double> conversionCosts;
 		conversionCosts = (ModelConfig.CONVERSION_COSTS_FROM_FILE) ? new ConversionCostReader().read() : new ConversionCostReader().calcFromConfig();
 		
-		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, currentWoodYieldData,
+		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, woodYieldData,
 				conversionCosts, carbonDemandIncrease);
 		internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep);
 		
@@ -476,7 +482,7 @@ public class ModelMain {
 		}
 
 		if (ModelConfig.IS_CALIBRATION_RUN) {
-			serializeLandUse(landUseRaster);
+			//serializeLandUse(landUseRaster);
 			countryAgents.serializeCropUsageForAll();
 			countryAgents.serializeWoodUsageForAll();
 			internationalMarket.serializeGlobalPrices();
@@ -551,6 +557,7 @@ public class ModelMain {
 		
 		for (CompositeCountry cc : countryGrouping) {
 			countryAgents.addForCountry(cc, cropUsageDataMap, initLU, woodAndCarbonProdMap);
+			globalLandUseRaster.putAll(initLU);
 		}
 	}
 
@@ -684,67 +691,12 @@ public class ModelMain {
 		RasterSet<LandUseItem> landUseRaster = new RasterSet<LandUseItem>(initialLC.getHeaderDetails());
 		
 		for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
-			landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue()));
-		}
-		
-		LogWriter.println("Reading forest growth data");
-		int initYear = new Timestep(0).getYear();
-		WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
-		for (int yearLower = WoodYieldItem.MIN_YEAR; yearLower < initYear; yearLower += WoodYieldItem.YEAR_INTERVAL) {
-			int yearUpper = yearLower + WoodYieldItem.YEAR_INTERVAL;
-			// Get wood yield data for current interval
-			if (yearLower ==  WoodYieldItem.MIN_YEAR) {
-				getWoodYieldData(woodYieldData, yearLower);
-				getWoodYieldData(woodYieldData, yearUpper);
-			} else {
-				// Remove data no longer needed
-				for (WoodYieldItem w : woodYieldData.values()) {
-					w.removeDataForYear(yearLower - WoodYieldItem.YEAR_INTERVAL);
-				}
-				getWoodYieldData(woodYieldData, yearUpper);
-			}
-			// Iterate over land use items
-			for (Map.Entry<RasterKey, LandUseItem > luEntry: landUseRaster.entrySet()) {
-				LandUseItem luItem = luEntry.getValue();
-				// Which data is needed
-				Map<LandCoverType, Set<Integer>> yearsNeeded = luItem.getTilesAllYearsEstablished();
-				// Wood yield data for grid cell
-				WoodYieldItem wYieldItem = woodYieldData.get(luEntry.getKey());
-				if (wYieldItem == null) {
-					wYieldItem = new WoodYieldItem();
-					wYieldItem.setDefaultForMissingData();
-				}
-				
-				/*
-				LogWriter.println("X=" + landUseRaster.getXCoordin(luEntry.getKey()) + " Y=" + landUseRaster.getYCoordin(luEntry.getKey()));
-				
-				for (RasterKey k : woodYieldData.keySet()) {
-					LogWriter.println(""+woodYieldData.getXCoordin(k)+","+woodYieldData.getYCoordin(k)+","+(woodYieldData.get(k)==null));
-				}
-				*/
-				
-				// For each forest type
-				for (Map.Entry<LandCoverType, Set<Integer>> ynEntry : yearsNeeded.entrySet()) {
-					LandCoverType lcType = ynEntry.getKey();
-					if (!lcType.isNatural()) {
-						continue;
-					}
-					// For each year needed
-					for (int y : ynEntry.getValue()) {
-						// Skip if year not in current interval
-						if (y < yearLower || y > yearUpper)
-							continue;
-						int minAge = initYear - y;
-						// For each age needed
-						for (int a = minAge; a <= WoodYieldItem.MAX_AGE; a++) {
-							double yield = wYieldItem.getInterpolatedYield(lcType, y, a);
-							luItem.setWoodYield(lcType, y, a, yield);
-						}
-					}
-				}
+			if (entry.getKey().getCol() == 177 && entry.getKey().getRow() == 145) {
+				int foo = 1;
 			}
+			landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue()));
 		}
-		
+
 		return landUseRaster;
 	}
 
@@ -774,24 +726,24 @@ public class ModelMain {
 	}
 	
 	/** Get wood yield data */
-	private WoodYieldRasterSet getWoodYieldData(WoodYieldRasterSet woodYieldData, int year) {
-		String rootDir = "C:\\Users\\Bart\\Documents\\PhD\\Carbon and wood yields\\forestry_dummy";
-		new WoodYieldReader(woodYieldData, year).getRasterDataFromFile(rootDir + File.separator + year + File.separator + "wood_yield.out");
-		return woodYieldData;		
+	private WoodYieldRasterSet getWoodYieldData(Timestep timestep) {
+		return woodYieldReader.getWoodYields(globalLandUseRaster, timestep, internationalMarket.getWoodPrice().getExportPrice());		
 	}
 	
+	
 	private WoodYieldRasterSet getInterpolatedWoodYieldData(Timestep timestep) {
 		WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
-		int yearCurrent = timestep.getYear();
+		/*int yearCurrent = timestep.getYear();
 		int yearLower = Math.floorDiv(yearCurrent, WoodYieldItem.YEAR_INTERVAL) * WoodYieldItem.YEAR_INTERVAL + WoodYieldItem.YEAR_OFFSET;
 		int yearUpper = yearLower + WoodYieldItem.YEAR_INTERVAL;
+		int age = 
 		getWoodYieldData(woodYieldData, yearLower);
 		getWoodYieldData(woodYieldData, yearUpper);
 		for (WoodYieldItem wyItem : woodYieldData.values()) {
 			wyItem.interpolateAllForYear(yearCurrent);
 			wyItem.removeDataForYear(yearLower);
 			wyItem.removeDataForYear(yearUpper);
-		}
+		}*/
 		return woodYieldData;
 	}
 
diff --git a/src/ac/ed/lurg/Timestep.java b/src/ac/ed/lurg/Timestep.java
index 48b7c7b2..c92f0b95 100644
--- a/src/ac/ed/lurg/Timestep.java
+++ b/src/ac/ed/lurg/Timestep.java
@@ -105,6 +105,19 @@ public class Timestep {
 		return getYearSubDir(rootDir, endYear);
 	}
 	
+	public String getWoodAndCarbonYearSubDir(String rootDir) {
+		int startYear;
+		
+		if (ModelConfig.CHANGE_YIELD_DATA_YEAR) {
+			startYear = (getYear()/ModelConfig.WOOD_AND_CARBON_TIMESTEP_SIZE) * ModelConfig.WOOD_AND_CARBON_TIMESTEP_SIZE;  // truncation in division as int
+		}
+		else {
+			startYear = ModelConfig.BASE_YEAR;
+		}
+		
+		return rootDir + File.separator + startYear;
+	}
+	
 	public static String getYearSubDir(String rootDir, int endYear) {
 		return rootDir + File.separator + (endYear-ModelConfig.LPJG_TIMESTEP_SIZE+1) + "-" + endYear;
 	}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 0bbff8c2..8a4136c7 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -115,11 +115,6 @@ public class CountryAgent extends AbstractCountryAgent {
 		// project demand
 		calculateCountryPricesAndDemand(worldPrices, carbonPrice, woodPrice, false);
 		
-		if (!currentTimestep.isInitialTimestep()) {
-			updateWoodYieldData(woodYieldData);
-			calculateForestRotation(woodPrice.getExportPrice(), ModelConfig.DISCOUNT_RATE);
-		}
-		
 		if (saveGamsGdxFiles && ModelConfig.PRICE_ELASTIC_DEMAND)
 			saveGDXFile("demand");
 		
@@ -137,7 +132,7 @@ public class CountryAgent extends AbstractCountryAgent {
 			}
 			
 			// optimize areas and intensity 
-			GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, carbonFluxData, conversionCosts, 
+			GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, woodYieldData, carbonFluxData, conversionCosts, 
 					carbonDemandIncrease);
 			GamsRasterOptimiser opti = new GamsRasterOptimiser(input, yieldClusters);
 			LogWriter.println("Running " + country.getName() + ", currentTimestep " + currentTimestep);
@@ -176,7 +171,7 @@ public class CountryAgent extends AbstractCountryAgent {
 	}
 
 	private GamsRasterInput getGamsRasterInput(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces, double gen2EcIncrease,
-			RasterSet<CarbonFluxItem> carbonFluxData, Map<LccKey, Double> conversionCosts,
+			RasterSet<WoodYieldItem> woodYieldData, RasterSet<CarbonFluxItem> carbonFluxData, Map<LccKey, Double> conversionCosts,
 			double carbonDemandIncrease) {
 		double allowedImportChange;
 
@@ -264,9 +259,9 @@ public class CountryAgent extends AbstractCountryAgent {
 
 		GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, 
 				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates, currentCarbonPrice, carbonTradeConstraint, currentWoodDemand, 
-				currentTimberPrice, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData(), getTotalRotaWoodHarvest());	
+				currentTimberPrice, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData(), harvestWoodFromRotation(woodYieldData));	
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, 
-				carbonFluxData, conversionCosts, countryLevelInputs);
+				woodYieldData, carbonFluxData, conversionCosts, countryLevelInputs);
 
 		return input;
 	}
@@ -354,49 +349,28 @@ public class CountryAgent extends AbstractCountryAgent {
 		
 	}
 	
-	public double getNetCarbonFlux() {
-		return previousGamsRasterOutput.getNetCarbonFlux();
-	}
-	
-	public Map<WoodType, WoodUsageData> getWoodUsageData() {
-		return previousGamsRasterOutput.getWoodUsageData();
-	}
-	
-	private void updateWoodYieldData(RasterSet<WoodYieldItem> woodYieldData) {
-		RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
-		for (Map.Entry<RasterKey, LandUseItem> entry : previousLandUses.entrySet()) {
+	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();
-			if (luItem == null)
-				continue;
-			
 			WoodYieldItem wyItem = woodYieldData.get(key);
-			if (wyItem == null) 
-				continue;
-			
-			luItem.copyWoodYields(wyItem);
+			if (wyItem == null)
+				continue; // TODO Deal with this properly
+			totalHarvest += wyItem.getCurrentRotationHarvest();
+			luItem.doForestRotation(currentTimestep, wyItem.getOptimalRotation()); // WARNING: Side effect
 		}
+		return totalHarvest;
 	}
-
-	private void calculateForestRotation(double woodPrice, double discountRate) {
-		RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
-		for (LandUseItem luItem : previousLandUses.values()) {
-			if (luItem == null)
-				continue;
-			
-			luItem.calculateOptimalRotation(woodPrice, discountRate, currentTimestep);;
-		}
+	
+	public double getNetCarbonFlux() {
+		return previousGamsRasterOutput.getNetCarbonFlux();
 	}
 	
-	private double getTotalRotaWoodHarvest() {
-		RasterSet<LandUseItem> previousLandUses = previousGamsRasterOutput.getLandUses();
-		double totalHarvest = 0.0;
-		for (LandUseItem luItem : previousLandUses.values()) {
-			if (luItem == null)
-				continue;
-			
-			totalHarvest += luItem.getWoodHarvestFromRotation(currentTimestep);
-		}
-		return totalHarvest;
+	public Map<WoodType, WoodUsageData> getWoodUsageData() {
+		return previousGamsRasterOutput.getWoodUsageData();
 	}
+
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index e5c4bcd8..082106bf 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -392,6 +392,7 @@ public class GamsLocationOptimiser {
 			String locString = Integer.toString(locationId);
 			WoodYieldData wYield = entry.getValue();
 
+			/*
 			for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
 				for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
 					Vector<String> v = new Vector<String>();
@@ -405,7 +406,30 @@ public class GamsLocationOptimiser {
 					double wYieldRota = (toLc.equals(LandCoverType.TIMBER_FOREST)) ? wYield.getYieldRota() : 0.0;
 					setGamsParamValue(woodYieldRotaP.addRecord(v), wYieldRota, 5);
 				}
-			}			
+			}
+			*/
+			
+			for (Map.Entry<LccKey, Double> wyEntry : wYield.getYieldLucMap().entrySet()) {
+				LandCoverType fromLc = wyEntry.getKey().getFromLc();
+				LandCoverType toLc = wyEntry.getKey().getToLc();
+				double yield = wyEntry.getValue();
+				
+				Vector<String> v = new Vector<String>();
+				v.add(fromLc.getName());
+				v.add(toLc.getName());
+				v.add(locString);
+			
+				setGamsParamValue(woodYieldLucP.addRecord(v), yield, 5);
+			}
+			
+			// Yield from timber forest rotation
+			for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
+				Vector<String> v = new Vector<String>();
+				v.add(LandCoverType.TIMBER_FOREST.getName());
+				v.add(toLc.getName());
+				v.add(locString);
+				setGamsParamValue(woodYieldRotaP.addRecord(v), wYield.getYieldRota(), 5);
+			}
 		}	
 		
 		// Land cover conversion cost
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index 32a1b30b..71f255e3 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -19,18 +19,20 @@ public class GamsRasterInput {
 	private YieldRaster yields;
 	private RasterSet<LandUseItem> previousLandsUses;
 	private RasterSet<IrrigationItem> irrigationCost;
+	private RasterSet<WoodYieldItem> woodYields;
 	private RasterSet<CarbonFluxItem> carbonFluxes;
 	private Map<LccKey, Double> conversionCosts;
 	private GamsCountryInput countryInput;
 
 	public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses, RasterSet<IrrigationItem> irrigationCost, 
-			RasterSet<CarbonFluxItem> carbonFluxes,
+			RasterSet<WoodYieldItem> woodYields, RasterSet<CarbonFluxItem> carbonFluxes,
 			Map<LccKey, Double> conversionCosts, GamsCountryInput countryInput) {
 		super();
 		this.timestep = timestep;
 		this.yields = yields;
 		this.previousLandsUses = previousLandsUses;
 		this.irrigationCost = irrigationCost;
+		this.woodYields = woodYields;
 		this.carbonFluxes = carbonFluxes;
 		this.conversionCosts = conversionCosts;
 		this.countryInput = countryInput;
@@ -48,6 +50,10 @@ public class GamsRasterInput {
 		return irrigationCost;
 	}
 	
+	public RasterSet<WoodYieldItem> getWoodYields() {
+		return woodYields;
+	}
+
 	public RasterSet<CarbonFluxItem> getCarbonFluxes() {
 		return carbonFluxes;
 	}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 4182c182..4c80bbc5 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -10,9 +10,12 @@ import java.util.Set;
 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;
+import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.YieldType;
@@ -191,7 +194,7 @@ public class GamsRasterOptimiser {
 			if (newLandUseItem!=null) {
 				double cellChange = (totalUnprotectedLC > 0) ? change * newLandUseItem.getUnprotectedLandCoverArea(fromType)/totalUnprotectedLC : change / keys.size();
 				double shortfall = newLandUseItem.moveAreas(toType, fromType, cellChange);
-				shortfall += newLandUseItem.moveUnprotectedArea(toType, fromType, cellChange);
+				shortfall += newLandUseItem.moveUnprotectedArea(toType, fromType, cellChange, rasterInputData.getTimestep());
 				if (shortfall == 0)
 					keysWithSpace.add(key);
 				else
@@ -215,7 +218,7 @@ public class GamsRasterOptimiser {
 		
 		RasterSet<CarbonFluxItem> carbonFluxRaster = rasterInputData.getCarbonFluxes();
 		
-//		RasterSet<WoodYieldItem> woodYieldRaster = rasterInputData.getWoodYields();
+		RasterSet<WoodYieldItem> woodYieldRaster = rasterInputData.getWoodYields();
 
 		// look for inconsistencies
 		for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
@@ -270,7 +273,7 @@ public class GamsRasterOptimiser {
 				}
 
 				IrrigationItem irrigItem = irrigRaster.get(key);
-				
+				WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
 				CarbonFluxItem carbonFluxItem = carbonFluxRaster.get(key);
 //				WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
 
@@ -320,29 +323,32 @@ public class GamsRasterOptimiser {
 					}
 				}
 				
-				for (LandCoverType forestType : LandCoverType.getNaturalTypes()) {
-					double wYieldThisTime = landUseItem.getWoodYield(forestType, rasterInputData.getTimestep());
-					double wYieldSoFar = aggWYield.getYieldLuc(forestType);
-					double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
-					aggWYield.setYieldLuc(forestType, wYieldAgg);		
-				}
-				
-				{
-					// Rotation wood yield
-					double wYieldThisTime = landUseItem.getWoodYieldAtRotation(currentYear);
-					double wYieldSoFar = aggWYield.getYieldRota();
-					double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
-					aggWYield.setYieldRota(wYieldAgg);	
+				if (woodYieldItem != null) { // TODO Fix
+					for (Map.Entry<LccKey, Double> wyEntry : woodYieldItem.getYieldMap().entrySet()) {
+						LccKey lccKey = wyEntry.getKey();
+						double yield = wyEntry.getValue();
+						double wYieldThisTime = woodYieldItem.getYield(lccKey);
+						double wYieldSoFar = aggWYield.getYieldLuc(lccKey);
+						double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
+						aggWYield.setYieldLuc(lccKey, wYieldAgg);		
+					}
+					
+					{
+						// Rotation wood yield
+						double wYieldThisTime = woodYieldItem.getYieldAtRotation();
+						double wYieldSoFar = aggWYield.getYieldRota();
+						double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
+						aggWYield.setYieldRota(wYieldAgg);	
+					}
+					
+					
+					// Rotation period
+					double forestRotaThisTime = woodYieldItem.getOptimalRotation();
+					double forestRotaSoFar = (aggregatedForestRotation.containsKey(clusterId)) ? aggregatedForestRotation.get(clusterId) : forestRotaThisTime;
+					double forestRotaAgg = aggregateMean(forestRotaSoFar, suitableAreaSoFar, forestRotaThisTime, suitableAreaThisTime);
+					aggregatedForestRotation.put(clusterId, forestRotaAgg);				
 				}
 				
-				
-				// Rotation period
-				double forestRotaThisTime = landUseItem.getForestRotation();
-				double forestRotaSoFar = (aggregatedForestRotation.containsKey(clusterId)) ? aggregatedForestRotation.get(clusterId) : forestRotaThisTime;
-				double forestRotaAgg = aggregateMean(forestRotaSoFar, suitableAreaSoFar, forestRotaThisTime, suitableAreaThisTime);
-				aggregatedForestRotation.put(clusterId, forestRotaAgg);
-				
-				
 				// Crops yields and area fractions
 				for (CropType crop : CropType.getNonMeatTypes()) {
 					if (!crop.equals(CropType.SETASIDE)) {
@@ -418,6 +424,13 @@ public class GamsRasterOptimiser {
 		checkedTotalAreas(landUseRaster, "before");
 		checkedTotalAreas(aggregatedAreas, "after");
 		
+		// 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))
+				aggregatedForestRotation.put(locId, meanRotation);
+		}
+		
 		return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, 
 				aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput(), aggregatedForestRotation);
 	}
diff --git a/src/ac/ed/lurg/forestry/ForestManager.java b/src/ac/ed/lurg/forestry/ForestManager.java
index 54608da4..3027c3f8 100644
--- a/src/ac/ed/lurg/forestry/ForestManager.java
+++ b/src/ac/ed/lurg/forestry/ForestManager.java
@@ -22,7 +22,7 @@ public class ForestManager implements Serializable {
 	public ForestManager() {
 		this.woodYields = new WoodYieldItem();
 	}
-	
+	/*
 	public void setWoodYield(LandCoverType forestType, int year, int age, double yield) {
 		woodYields.setYield(forestType, year, age, yield);
 	}
@@ -72,6 +72,6 @@ public class ForestManager implements Serializable {
 	public double getYieldAtRotation(int year) {
 		return woodYields.getYield(LandCoverType.TIMBER_FOREST, year, optimalRotation);
 	}
-	
+	*/
 
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldData.java b/src/ac/ed/lurg/forestry/WoodYieldData.java
index a3c878b3..ab8b8b00 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldData.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldData.java
@@ -1,26 +1,43 @@
 package ac.ed.lurg.forestry;
 
+import java.util.HashMap;
 import java.util.Map;
 
+import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
 
 public class WoodYieldData {
 	private double yieldRota; // wood yield at current rotation length
-	private Map<LandCoverType, Double> yieldLuc; // wood yield from clearing current forest
+	private Map<LccKey, Double> yieldLuc = new HashMap<LccKey, Double>(); // wood yield from clearing current forest
 	
 	public void setYieldRota(double yield) {
 		yieldRota = yield;
 	}
 	
-	public void setYieldLuc(LandCoverType forestType, double yield) {
-		yieldLuc.put(forestType, yield);
+	public void setYieldLuc(LccKey key, double yield) {
+		yieldLuc.put(key, yield);
 	}
 	
 	public double getYieldRota() {
 		return yieldRota;
 	}
 	
-	public double getYieldLuc(LandCoverType forestType) {
-		return yieldLuc.get(forestType);
+	public double getYieldLuc(LandCoverType fromLc, LandCoverType toLc) {
+		LccKey key = new LccKey(fromLc, toLc);
+		if (yieldLuc.containsKey(key)) 
+			return yieldLuc.get(key);
+		
+		return Double.NaN;
+	}
+	
+	public double getYieldLuc(LccKey key) {
+		if (yieldLuc.containsKey(key)) 
+			return yieldLuc.get(key);
+		
+		return Double.NaN;
+	}
+	
+	public Map<LccKey, Double> getYieldLucMap() {
+		return yieldLuc;
 	}
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java
index 25ae9ed6..c04c4b23 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldItem.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java
@@ -1,97 +1,115 @@
 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 java.util.stream.Collectors;
 
+import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.Timestep;
-import ac.ed.lurg.carbon.AgeFunction;
+import ac.ed.lurg.landuse.LandUseTile;
+import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.BilinearInterpolator;
-import ac.ed.lurg.utils.DoubleMap;
-import ac.ed.lurg.utils.Interpolator;
-import ac.ed.lurg.utils.LogWriter;
 import ac.sac.raster.RasterItem;
 
-public class WoodYieldItem implements RasterItem, Serializable {
-	private static final long serialVersionUID = -4822328490070047362L;
-	public static final int MAX_AGE = 150;
+public class WoodYieldItem implements RasterItem {
+	public static final int MAX_AGE = 250;
 	public static final int YEAR_INTERVAL = 20;
 	public static final int YEAR_OFFSET = 10;
 	public static final int AGE_INTERVAL = 1;
 	public static final int MIN_YEAR = 1850;
 	public static final int MAX_YEAR = 2100;
 	
-	private Map<WoodYieldKey, Double> data = new HashMap<WoodYieldKey, Double>();
+	private Map<LccKey, Double> yield = new HashMap<LccKey, Double>();
+	private int optimalRotation; // Only applies to TIMBER_FOREST
+	private double yieldAtRotation; // Only applies to TIMBER_FOREST
+	private double currentRotationHarvest; // timber harvested from current rotation
 	
-	public WoodYieldItem() {}
-	
-	public WoodYieldItem(WoodYieldItem toCopy) {
-		data.putAll(toCopy.data);
-	}
 	
-	public void setYield(LandCoverType forestType, int year, int age, double yield) {
-		WoodYieldKey key = new WoodYieldKey(forestType, year, age);
-		data.put(key, yield);
-	}
-	
-	public void copyAll(WoodYieldItem woodYieldItem) {
-		data.putAll(woodYieldItem.data);
-	}
+	public WoodYieldItem() {}
 	
-	public void setDefaultForMissingData() {
-		for (LandCoverType forestType : LandCoverType.getNaturalTypes()) {
-			for (int y = MIN_YEAR; y <= MAX_YEAR; y+=YEAR_INTERVAL) {
-				for (int t = 0; t <= MAX_AGE; t+=AGE_INTERVAL) {
-					setYield(forestType, y, t, 0);
+	public void calcYieldData(Map<LccKey, Double[]> woodYields, Map<LandCoverType, ArrayList<LandUseTile>> 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());
+			
+			if (tiles.size() == 0) {
+				this.yield.put(key, 0.0) ;
+			} else {				
+				double totalYield = 0;
+				double totalArea = 0;
+				for (LandUseTile tile : tiles) {
+					if (tile.isProtected())
+						continue;
+					int age = tile.getAge(timestep);
+					totalYield += yields[age] * tile.getArea();
+					totalArea += tile.getArea();
 				}
-			}	
+				double meanYield = totalYield / totalArea;
+				this.yield.put(key, meanYield);
+			}
 		}
 	}
 	
-	public double getYield(LandCoverType forestType, int year, int age) {
-		age = Math.min(age, MAX_AGE); // assuming no further growth
-		WoodYieldKey key = new WoodYieldKey(forestType, year, age);
-		Double d = data.get(key);
-		if (d == null) {
-			for (Map.Entry<WoodYieldKey, Double> entry : data.entrySet()) {
-				LogWriter.println(entry.getKey().getYear()+","+entry.getKey().getAge()+","+entry.getValue());
-			}
-			int foo = 1;
+	public void calcRotationData(Map<LccKey, Double[]> woodYields, Map<LandCoverType, ArrayList<LandUseTile>> landUseTiles, Timestep timestep, double woodPrice) {
+		// Optimal rotation
+		Map<Double, Integer> levMap = new HashMap<Double, Integer>(); // Land Expected Value
+		
+		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);
+		}
+		// 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) {
+			currentRotationHarvest = 0.0;
+		} else {
+			double areaHarvested = timberTiles.stream()
+					.filter(o -> !o.isProtected() && o.getAge(timestep) >= optimalRotation)
+					.mapToDouble(o -> o.getArea())
+					.sum();
+			currentRotationHarvest = areaHarvested * yieldAtRotation;
 		}
-		return d;
 	}
 	
-	public double getInterpolatedYield(LandCoverType forestType, int year, int age) {
-		int yearLower = Math.floorDiv(year - YEAR_OFFSET, YEAR_INTERVAL) * YEAR_INTERVAL + YEAR_OFFSET;
-		int yearUpper = yearLower + YEAR_INTERVAL;
-		if (yearUpper > MAX_YEAR)
-			yearUpper = yearLower;
-		double q1 = getYield(forestType, yearLower, age);
-		double q2 = getYield(forestType, yearUpper, age);
-		double d = q1 + (year - yearLower) * ((q2 - q1) / (yearUpper - yearLower));;
-		return d;		
+	public void setDefaultForMissingData() {
+		
 	}
 	
-	public void interpolateAllForYear(int year) {
-		for (LandCoverType lcType: LandCoverType.getNaturalTypes()) {
-			for (int age = 0; age <= MAX_AGE; age++) {
-				double yield = getInterpolatedYield(lcType, year, age);
-				setYield(lcType, year, age, yield);
-			}
-		}
+	public double getYield(LandCoverType fromLc, LandCoverType toLc) {
+		return yield.get(new LccKey(fromLc, toLc));
 	}
 	
-	public void removeDataForYear(int year) {
-		for (WoodYieldKey key : data.keySet()) {
-			if (key.getYear() == year)
-				data.remove(key);
-		}
+	public double getYield(LccKey key) {
+		return yield.get(key);
 	}
 	
-	public void removePastData(Timestep timestep) {
-		// TODO
+	public Map<LccKey, Double> getYieldMap() {
+		return yield;
 	}
+
+	public int getOptimalRotation() {
+		return optimalRotation;
+	}
+
+	public double getYieldAtRotation() {
+		return yieldAtRotation;
+	}
+
+	public double getCurrentRotationHarvest() {
+		return currentRotationHarvest;
+	}
+	
+	
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldReader.java b/src/ac/ed/lurg/forestry/WoodYieldReader.java
index 10a1f56f..d1d1557c 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldReader.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldReader.java
@@ -1,27 +1,97 @@
 package ac.ed.lurg.forestry;
 
+import java.io.File;
+import java.util.HashMap;
 import java.util.Map;
 
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.Timestep;
+import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
 import ac.sac.raster.AbstractTabularRasterReader;
+import ac.sac.raster.RasterHeaderDetails;
 import ac.sac.raster.RasterKey;
 import ac.sac.raster.RasterSet;
 
-public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem>{
-	private static final int MIN_COLS = 4;
-	private static final double conversionFactor = 10.0; // convert kgC/m2 to tC/ha
-	private int year;
+public class WoodYieldReader {
+	private static final int MIN_COLS = 253;
+	private static final String DELIMITER = "[ |\t]+";
+	private static final double CONVERSION_FACTOR = 10.0; // convert kgC/m2 to tC/ha
+	private RasterHeaderDetails rasterProj;
 	
-	public WoodYieldReader(WoodYieldRasterSet woodYieldData, int year) {
-		super("[ |\t]+", MIN_COLS, woodYieldData);
-		this.year = year;
+	public WoodYieldReader(RasterHeaderDetails rasterProj) {
+		this.rasterProj = rasterProj;
 	}
+	
+	public WoodYieldRasterSet getWoodYields(RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
+		WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(rasterProj);
+		getWoodYieldsForest(woodYieldData, landUseRaster, timestep, woodPrice);
+		getWoodYieldsAgri(woodYieldData, landUseRaster, timestep, woodPrice);
+		return woodYieldData;
+	}
+
+	public WoodYieldRasterSet getWoodYieldsForest(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
+
+		AbstractTabularRasterReader<WoodYieldItem> yieldReader = new AbstractTabularRasterReader<WoodYieldItem>(DELIMITER, MIN_COLS, woodYieldData) {
+			protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
 
-	@Override
-	protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
-		item.setYield(LandCoverType.TIMBER_FOREST, year, (int) getValueForCol(rowValues, "Time"), getValueForCol(rowValues, "yield") * conversionFactor);
-		item.setYield(LandCoverType.CARBON_FOREST, year, (int) getValueForCol(rowValues, "Time"), getValueForCol(rowValues, "yield") * conversionFactor);
-		item.setYield(LandCoverType.NATURAL,       year, (int) getValueForCol(rowValues, "Time"), getValueForCol(rowValues, "yield") * conversionFactor);                                                                        
+				Double[] yields = getYieldsFromRowValues(rowValues);
+				Map<LccKey, Double[]> yieldMap = new HashMap<LccKey, Double[]>();
+				// For now, same yields for each land cover.
+				for (LandCoverType fromLc : LandCoverType.getNaturalTypes()) {
+					for (LandCoverType toLc : LandCoverType.getNaturalTypes()) {
+						LccKey lccKey = new LccKey(fromLc, toLc);
+						yieldMap.put(lccKey, yields);
+					}	
+				}
+				
+				item.calcYieldData(yieldMap, landUseRaster.get(key).getLandUseTiles(), timestep);
+				item.calcRotationData(yieldMap, landUseRaster.get(key).getLandUseTiles(), timestep, woodPrice);
+			}			
+		};
+	
+		String filename = getDataDir(ModelConfig.WOOD_YIELD_FRST_FILENAME, timestep);
+		yieldReader.getRasterDataFromFile(filename);
+		
+		return woodYieldData;
+	}
+	
+	public WoodYieldRasterSet getWoodYieldsAgri(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, Timestep timestep, double woodPrice) {
+
+		AbstractTabularRasterReader<WoodYieldItem> yieldReader = new AbstractTabularRasterReader<WoodYieldItem>(DELIMITER, MIN_COLS, woodYieldData) {
+			protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
+
+				Double[] yields = getYieldsFromRowValues(rowValues);
+				Map<LccKey, Double[]> yieldMap = new HashMap<LccKey, Double[]>();
+				// For now, same yields for each land cover.
+				for (LandCoverType fromLc : LandCoverType.getNaturalTypes()) {
+					for (LandCoverType toLc : LandCoverType.getAgriculturalTypes()) {
+						LccKey lccKey = new LccKey(fromLc, toLc);
+						yieldMap.put(lccKey, yields);
+					}	
+				}
+				
+				item.calcYieldData(yieldMap, landUseRaster.get(key).getLandUseTiles(), timestep);
+			}			
+		};
+		
+		String filename = getDataDir(ModelConfig.WOOD_YIELD_AGRI_FILENAME, timestep);
+		yieldReader.getRasterDataFromFile(filename);
+		
+		return woodYieldData;
+	}
+	
+	private String getDataDir(String filename, Timestep timestep) {
+		return timestep.getWoodAndCarbonYearSubDir(ModelConfig.WOOD_AND_CARBON_DIR) + File.separator + filename;
+	}
+	
+	private Double[] getYieldsFromRowValues(Map<String, Double> rowValues) {
+		Double[] yields = new Double[251];
+		for (int i=0; i<=250; i++) {
+			yields[i] = rowValues.get(Integer.toString(i)) * CONVERSION_FACTOR;
+		}
+		return yields;
 	}
 
 }
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index 368b1ffa..a378194e 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -30,7 +30,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	private double protectedArea; //protected area in Mha
 	private double unavailableArea; //area unavailable due to altitude etc 
 	private double suitableArea; //Mha
-	private ForestManager forestManager;
 	
 	public LandUseItem() {}
 	
@@ -45,7 +44,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 			setUnavailableArea(landCover.getUnavailableFract());
 			updateSuitableArea(ModelConfig.BASE_YEAR);
 			addLandUseTiles(landCover.getLandUseAgeDist()); // Should only be run after unprotected areas calculated
-			forestManager = new ForestManager();
 			runAreaCheck(); // check that tile areas add up to total areas
 		}
 	}
@@ -59,39 +57,38 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		landUseTiles.putAll(luItemToCopy.landUseTiles);
 		protectedArea = (luItemToCopy.protectedArea);
 		suitableArea = (luItemToCopy.suitableArea);
-		forestManager = (luItemToCopy.forestManager);
 	}
 	
 	
-	public void moveTilesArea(LandCoverType toType, LandCoverType fromType, double area) {
+	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;
 		}
-	
+
 		// Remove area from tiles
-			double totalCurrentArea = fromTiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
-			for (LandUseTile s: fromTiles) {
-				if (s.isProtected())
-					continue;
+		double totalCurrentArea = fromTiles.stream().filter(o -> !o.isProtected()).mapToDouble(o -> o.getArea()).sum();
+		for (LandUseTile s: fromTiles) {
+			if (s.isProtected())
+				continue;
+
+			double areaToRemove = s.getArea() / totalCurrentArea * area; 
+			s.removeArea(areaToRemove);
+		}		
 
-				double areaToRemove = s.getArea() / totalCurrentArea * area; 
-				s.removeArea(areaToRemove);
-			}		
-		
 		// Add new tiles	
-			LandUseTile newStand = new LandUseTile(toType, area, 0, false);
-			toTiles.add(newStand);	
-			
-			runAreaCheck();
+		LandUseTile newStand = new LandUseTile(toType, area, timestep.getYear(), false);
+		toTiles.add(newStand);	
+
+		runAreaCheck();
 	}
 	
 	public void runAreaCheck() {
@@ -133,9 +130,16 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 				int maxAgeGroup = ageMap.keySet().stream().max(Integer::compare).get();
 				ageMap.put(maxAgeGroup, 1.0);
 			}
-
+			if (totalProp > 1.0) {
+				// normalise
+				for (Map.Entry<Integer, Double> entry : ageMap.entrySet()) {
+					ageMap.put(entry.getKey(), entry.getValue() / totalProp);
+				}
+			}
+			
+			// Assuming land use distributed evenly in each age group
 			for (Map.Entry<Integer, Double> ageEntry : ageMap.entrySet()) {
-				double prop = ageEntry.getValue() / 15;
+				double prop = ageEntry.getValue() / 10;
 				if (prop == 0.0)
 					continue;
 				int ageGroup = ageEntry.getKey();
@@ -170,12 +174,9 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return reqData;
 	}
 	
-	public void setWoodYield(LandCoverType forestType, int year, int age, double yield) {
-		forestManager.setWoodYield(forestType, year, age, yield);
-	}
 	
-	public void copyWoodYields(WoodYieldItem woodYieldItem) {
-		forestManager.copyWoodYields(woodYieldItem);
+	public Map<LandCoverType, ArrayList<LandUseTile>> getLandUseTiles() {
+		return landUseTiles;
 	}
 	
 	public void setProtectedFract(double protFrac) {
@@ -343,7 +344,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
         return changeReq - actuallyChanged;
     }
     
-	public double moveUnprotectedArea(LandCoverType toType, LandCoverType fromType, double area) {
+	public double moveUnprotectedArea(LandCoverType toType, LandCoverType fromType, double area, Timestep timestep) {
 		
         double prevTo = getUnprotectedLandCoverArea(toType);
         double prevFrom = getUnprotectedLandCoverArea(fromType);
@@ -362,7 +363,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
         }
         */
         
-        moveTilesArea(toType, fromType, areaMoved);
+        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;
@@ -526,37 +530,16 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		
 		return changes;
 	}
-
-	public void calculateOptimalRotation(double woodPrice, double discountRate, Timestep timestep) {
-		forestManager.calculateOptimalRotation(woodPrice, discountRate, timestep);
-	}
 	
-	public double getWoodYield(LandCoverType forestType, Timestep timestep) {
-		return forestManager.getAverageYield(landUseTiles.get(forestType), timestep);
-	}
-	
-	public double getWoodYieldAtRotation(int year) {
-		return forestManager.getYieldAtRotation(year);
-	}
-	
-	public double getWoodHarvestFromRotation(Timestep timestep) {
+	public void doForestRotation(Timestep timestep, double rotationAge) {
 		ArrayList<LandUseTile> tiles = landUseTiles.get(LandCoverType.TIMBER_FOREST);
 		if (tiles == null) 
-			return 0;
-		double totalHarvest = 0;
-		int rotationAge = forestManager.getOptimalRotation();
+			return;
 		for (LandUseTile t : tiles) {
 			if (t.getAge(timestep) < rotationAge)
 				continue;
-			double yield = forestManager.getStandYield(t, timestep);
-			totalHarvest += yield * t.getArea();
-			t.resetAge(timestep); // Stand age to 0
+			t.resetAge(timestep); // Set stand age to 0
 		}
-		return totalHarvest;
-	}
-	
-	public int getForestRotation() {
-		return forestManager.getOptimalRotation();
 	}
 
 	private boolean isZeroOrNull(Double d) {
diff --git a/src/ac/ed/lurg/landuse/LandUseTile.java b/src/ac/ed/lurg/landuse/LandUseTile.java
index 9bde2d77..7af6fb13 100644
--- a/src/ac/ed/lurg/landuse/LandUseTile.java
+++ b/src/ac/ed/lurg/landuse/LandUseTile.java
@@ -49,7 +49,7 @@ public class LandUseTile implements Serializable, InterpolatingRasterItem<LandUs
 
 	public void removeArea(double a) {
 		if (a > area) {
-			LogWriter.printlnError("LandUseItem: attempting to remove too much area: " + a + ", available: " + area);
+			LogWriter.printlnError("LandUseTile.removeArea: attempting to remove too much area: " + a + ", available: " + area);
 			area -= Math.min(area, a);
 		} else {
 			area -= a;
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index fc9346e8..311704e9 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -112,5 +112,17 @@ public enum LandCoverType {
 		return naturalTypes;
 
 	}
+	
+	public static Collection<LandCoverType> getAgriculturalTypes() {
+
+		Collection<LandCoverType> agriTypes = new HashSet<LandCoverType>();
+
+		for (LandCoverType c : values())
+			if (!c.isNatural && c.isConvertible)
+				agriTypes.add(c);
+
+		return agriTypes;
+
+	}
 }
 
-- 
GitLab