Interpolate missing years

src/ac/ed/lurg/ModelMain.java

@@ -35,6 +35,7 @@ import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.yield.LPJYieldResponseMapReader;
 import ac.ed.lurg.yield.YieldRaster;
 import ac.sac.raster.IntegerRasterItem;
+import ac.sac.raster.InterpolatingRasterSet;
 import ac.sac.raster.RasterHeaderDetails;
 import ac.sac.raster.RasterKey;
 import ac.sac.raster.RasterOutputer;
@@ -46,9 +47,13 @@ public class ModelMain {
 	private CountryBoundaryRaster countryBoundaryRaster;
 	private DemandManager demandManager;
 	private CompositeCountryManager compositeCountryManager;
-	private Map<CropType, Double> prevWorldPrices;
 	private RasterHeaderDetails desiredProjection;
 
+	private Map<CropType, Double> prevWorldPrices;
+	private RasterSet<IntensitiesItem> prevIntensityRaster;
+	private RasterSet<AreasItem> prevCropAreaRaster;
+	
+	
 	public static void main(String[] args)  {
 		ModelMain theModel = new ModelMain();
 		theModel.setup();
@@ -82,7 +87,7 @@ public class ModelMain {
 				doTimestep(timestep);
 			}
 			catch (Exception e) {
-				LpjgOutputer.writeMarkerFile(timestep, true);
+				LpjgOutputer.writeMarkerFile(timestep.getYear(), true);
 				throw new RuntimeException(e);
 			}
 		}
@@ -161,8 +166,13 @@ public class ModelMain {
 			prevWorldPrices.put(crop, adjustedPrice);
 		}
 		
+		
 		// output results
 		outputTimestepResults(timestep, globalIntensityRaster, globalCropAreaRaster, globalLocationIdRaster, yieldSurfaces);
+		
+		// keep last to allow interpolation
+		prevIntensityRaster = globalIntensityRaster;
+		prevCropAreaRaster = globalCropAreaRaster;
 	}
 	
 	public double updateMarketPrices(double previousPrice, double demand, double supply) {
@@ -232,10 +242,28 @@ public class ModelMain {
 			RasterSet<IntegerRasterItem> locationIdRaster, YieldRaster yieldSurfaces) {
 		
 		writeMarketFile(timestep, cropAreaRaster);
-
+		
 		if (ModelConfig.OUTPUT_FOR_LPJG) {
-			LpjgOutputer lpjOutputer = new LpjgOutputer(timestep, intensityRaster, cropAreaRaster, yieldSurfaces);
-			lpjOutputer.writeOutput();
+
+			for (int outputYear : timestep.getYearsFromLast()) {
+				LogWriter.printlnError("Year" + outputYear);
+				LpjgOutputer lpjOutputer;
+				if (outputYear == timestep.getYear()) {
+					lpjOutputer = new LpjgOutputer(timestep.getYear(), intensityRaster, cropAreaRaster, yieldSurfaces);
+				}
+				else {
+					InterpolatingRasterSet<AreasItem> intermediateAreas = new InterpolatingRasterSet<AreasItem>(cropAreaRaster.getHeaderDetails()) {
+						private static final long serialVersionUID = 1306045141011047760L;
+						protected AreasItem createRasterData() {
+							return new AreasItem();
+						}
+					};
+					intermediateAreas.setup(prevCropAreaRaster, cropAreaRaster, timestep.getPreviousTimestep().getYear(), timestep.getYear(), outputYear);
+					lpjOutputer = new LpjgOutputer(outputYear, intensityRaster, intermediateAreas, yieldSurfaces);
+				}
+			
+				lpjOutputer.writeOutput();
+			}
 		}
 		
 		outputLandCover(timestep.getYear(), cropAreaRaster, LandCoverType.CROPLAND);
@@ -273,9 +301,9 @@ public class ModelMain {
 		for (CompositeCountry cc : countryGrouping) {
 			
 			// DEBUG code
-		//	if (!(cc.getName().equals("United States of America") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
-		//		continue;
-		//	}
+			if (!(cc.getName().equals("United States of America") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
+				continue;
+			}
 			
 			List<RasterKey> keys = countryBoundaryRaster.getKeysFor(cc);
 			RasterSet<LandCoverItem> initCountryLC = initLC.popSubsetForKeys(new RasterSet<LandCoverItem>(initLC.getHeaderDetails()), keys);

src/ac/ed/lurg/Timestep.java

 package ac.ed.lurg;
 
+import java.util.ArrayList;
+import java.util.Collection;
+
 
 public class Timestep {
 
@@ -63,4 +66,20 @@ public class Timestep {
 			return false;
 		return true;
 	}
+
+	public Collection<Integer> getYearsFromLast() {
+		Collection<Integer> years = new ArrayList<Integer>();
+		int currentYear = getYear();
+		
+		if (isInitialTimestep()) {
+			years.add(currentYear);
+			return years;
+		}
+		
+		for (int i=0; i<ModelConfig.TIMESTEP_SIZE; i++) {
+			years.add(currentYear - i);
+		}
+		
+		return years;
+	}
 }

src/ac/ed/lurg/landuse/AreasItem.java

@@ -7,12 +7,12 @@ import java.util.Map.Entry;
 import ac.ed.lurg.types.CropToDouble;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandCoverType;
-import ac.sac.raster.RasterItem;
+import ac.sac.raster.InterpolatingRasterItem;
 
 /*
  * Data for one year and country, but it does not know which one
  */
-public class AreasItem implements RasterItem {
+public class AreasItem implements InterpolatingRasterItem<AreasItem> {
 
 	private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
 	private Map<LandCoverType, Double> landCoverAreas = new HashMap<LandCoverType, Double>();
@@ -115,4 +115,37 @@ public class AreasItem implements RasterItem {
 		
 		return changes;
 	}
+
+	@Override
+	public void interpolateAll(AreasItem fromItem, AreasItem toItem, double factor) {		
+		cropFractions = new HashMap<CropType, Double>();
+		landCoverAreas = new HashMap<LandCoverType, Double>();
+
+		for (CropType crop : CropType.values()) {
+			Double from = fromItem.cropFractions.get(crop);
+			Double to = toItem.cropFractions.get(crop);
+			Double d = interpolate(from, to, factor);
+			cropFractions.put(crop, d);
+		}
+		
+		for (LandCoverType landCover : LandCoverType.values()) {
+			Double from = fromItem.landCoverAreas.get(landCover);
+			Double to = toItem.landCoverAreas.get(landCover);
+			Double d = interpolate(from, to, factor);
+			landCoverAreas.put(landCover, d);
+		}
+	}
+	
+	private Double interpolate(Double from, Double to, double factor) {
+		if (from == null) {
+			if (to == null) 
+				return 0.0;
+			return to;
+		}
+		else if (to == null)
+			return from;
+			
+		Double res = from + factor * (to - from);
+		return res;
+	}
 }

src/ac/ed/lurg/output/LpjgOutputer.java

@@ -8,7 +8,6 @@ import java.io.IOException;
 import java.util.Map.Entry;
 
 import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.Timestep;
 import ac.ed.lurg.landuse.AreasItem;
 import ac.ed.lurg.landuse.IntensitiesItem;
 import ac.ed.lurg.types.CropType;
@@ -24,25 +23,25 @@ public class LpjgOutputer {
 	private RasterSet<IntensitiesItem> intensityRaster;
 	private RasterSet<AreasItem> cropAreaRaster;
 	private YieldRaster yieldSurfaces;
-	private Timestep timestep;
+	private int year;
 
-	public LpjgOutputer(Timestep timestep, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, YieldRaster yieldSurfaces) {
-		this.timestep = timestep;
+	public LpjgOutputer(int year, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, YieldRaster yieldSurfaces) {
+		this.year = year;
 		this.intensityRaster = intensityRaster;
 		this.cropAreaRaster = cropAreaRaster;
 		this.yieldSurfaces = yieldSurfaces;
 	}
 
 	public void writeOutput() {
-		File outputDir = getOutputDir(timestep);
+		File outputDir = getOutputDir(year);
 
 		writeLandCoverAndCrop(outputDir);
 		writeIntensity(outputDir);
-		writeMarkerFile(timestep, false);
+		writeMarkerFile(year, false);
 	}
 
-	private static File getOutputDir(Timestep timestep) {
-		String outputDirName = ModelConfig.OUTPUT_DIR + File.separator + timestep.getYear();
+	private static File getOutputDir(int year) {
+		String outputDirName = ModelConfig.OUTPUT_DIR + File.separator + year;
 		File outputDir = new File(outputDirName);
 		outputDir.mkdirs();
 		return outputDir;
@@ -53,8 +52,8 @@ public class LpjgOutputer {
 		
 	}
 
-	public static void writeMarkerFile(Timestep timestep, boolean errorStatus) {
-		File outputDir = getOutputDir(timestep);
+	public static void writeMarkerFile(int year, boolean errorStatus) {
+		File outputDir = getOutputDir(year);
 
 		File markerFile = new File(outputDir.getPath() + File.separator + (errorStatus ? "error" : "done"));
 		
@@ -101,7 +100,7 @@ public class LpjgOutputer {
 				double forest = item.getLandCoverArea(LandCoverType.FOREST);
 				double otherNatural = item.getLandCoverArea(LandCoverType.OTHER_NATURAL);
 				double barren = item.getLandCoverArea(LandCoverType.BARREN);
-				landCoverWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, timestep.getYear(), crop/area, pasture/area, (forest+otherNatural)/area, barren/area));
+				landCoverWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, year, crop/area, pasture/area, (forest+otherNatural)/area, barren/area));
 				landCoverWriter.newLine();
 
 				boolean isSpringWheat = false;
@@ -115,7 +114,7 @@ public class LpjgOutputer {
 				double lpjSpringWheatFrac = item.getCropFraction(CropType.PULSES, CropType.STARCHY_ROOTS, isSpringWheat ? CropType.WHEAT : null);
 				double lpjCornFrac = item.getCropFraction(CropType.MAIZE);
 				double lpjRiceFrac = item.getCropFraction(CropType.RICE);
-				cropFractWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, timestep.getYear(), lpjWinterWheatFrac, lpjSpringWheatFrac, lpjCornFrac, lpjRiceFrac));
+				cropFractWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f", lat, lon, year, lpjWinterWheatFrac, lpjSpringWheatFrac, lpjCornFrac, lpjRiceFrac));
 				cropFractWriter.newLine();
 			}
 		}

src/ac/sac/raster/InterpolatingRasterItem.java

+package ac.sac.raster;
+
+public interface InterpolatingRasterItem<I> extends RasterItem {
+
+	void interpolateAll(I fromItem, I toItem, double factor);
+}

src/ac/sac/raster/InterpolatingRasterSet.java

+package ac.sac.raster;
+
+import ac.ed.lurg.utils.LogWriter;
+
+public class InterpolatingRasterSet<D extends InterpolatingRasterItem<D>> extends RasterSet<D> {
+
+	private static final long serialVersionUID = 7372086730280503581L;
+
+	public InterpolatingRasterSet(RasterHeaderDetails header) {
+		super(header);
+	}
+
+	public void setup(RasterSet<D> fromSet, RasterSet<D> toSet, double fromValue, double toValue, double desiredValue) {
+		clear();
+		
+		if ( !(isConstistentWithHeader(toSet.getHeaderDetails()) && isConstistentWithHeader(fromSet.getHeaderDetails())) ) {
+			throw new RuntimeException("Raster headers inconsistent");
+		}
+		
+		double factor = ((double)(desiredValue - fromValue)) / (toValue - fromValue);
+		if (factor > 1 || factor < 0) {
+			LogWriter.printlnError("InterpolatingRasterSet: Not interpolating, but extrapolating " + factor);
+		}
+		
+		for (RasterKey key : fromSet.keySet()) {
+			D fromItem = fromSet.get(key);
+			D toItem = toSet.get(key);
+			D newItem = createRasterData();
+
+			if (fromItem == null || toItem == null)
+				LogWriter.printlnError("InterpolatingRasterSet: Got nulls for key" + key);
+			else {
+				newItem.interpolateAll(fromItem, toItem, factor);
+				put(key, newItem);
+			}
+		}
+	}
+
+}