diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java index ec2acceef7dd96c6ed6c7f07dbbefa24ed731de8..5d1694871c967437a292393de22fe4a5d12c08be 100644 --- a/src/ac/ed/lurg/ModelMain.java +++ b/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); diff --git a/src/ac/ed/lurg/Timestep.java b/src/ac/ed/lurg/Timestep.java index aa783df4ad6b64e00c2bc8f759f51c78be2f0f1d..b72173612081137760cf74714fcb8cac3b321fbf 100644 --- a/src/ac/ed/lurg/Timestep.java +++ b/src/ac/ed/lurg/Timestep.java @@ -1,5 +1,8 @@ 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; + } } diff --git a/src/ac/ed/lurg/landuse/AreasItem.java b/src/ac/ed/lurg/landuse/AreasItem.java index daffeac04ecadcdca553ba4f35207b6b25075011..01b0181deed0542085fb6c0d0cfe0a8b8738b15c 100644 --- a/src/ac/ed/lurg/landuse/AreasItem.java +++ b/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; + } } diff --git a/src/ac/ed/lurg/output/LpjgOutputer.java b/src/ac/ed/lurg/output/LpjgOutputer.java index 3f208c455f3f01b7b4beed674712333733e91311..ff67cce61d4679cdc838e627eceaa13b85a8a2c7 100644 --- a/src/ac/ed/lurg/output/LpjgOutputer.java +++ b/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(); } } diff --git a/src/ac/sac/raster/InterpolatingRasterItem.java b/src/ac/sac/raster/InterpolatingRasterItem.java new file mode 100644 index 0000000000000000000000000000000000000000..9a2a89b225eda3186b5b3ee4d9aeecfa19f07ba3 --- /dev/null +++ b/src/ac/sac/raster/InterpolatingRasterItem.java @@ -0,0 +1,6 @@ +package ac.sac.raster; + +public interface InterpolatingRasterItem<I> extends RasterItem { + + void interpolateAll(I fromItem, I toItem, double factor); +} diff --git a/src/ac/sac/raster/InterpolatingRasterSet.java b/src/ac/sac/raster/InterpolatingRasterSet.java new file mode 100644 index 0000000000000000000000000000000000000000..6302230fd8d1a6f9cb052fdeeb69713514cf5ed3 --- /dev/null +++ b/src/ac/sac/raster/InterpolatingRasterSet.java @@ -0,0 +1,39 @@ +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); + } + } + } + +}