ModelMain.java

package ac.ed.lurg;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.FileWriter;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.nio.file.FileSystems;
import java.nio.file.Files;
import java.nio.file.StandardCopyOption;
import java.util.Collection;
import java.util.HashSet;
import java.util.List;
import java.util.Map;

import ac.ed.lurg.country.CompositeCountry;
import ac.ed.lurg.country.CompositeCountryManager;
import ac.ed.lurg.country.CountryAgent;
import ac.ed.lurg.country.CountryBoundaryRaster;
import ac.ed.lurg.country.CountryBoundaryReader;
import ac.ed.lurg.country.CountryPrice;
import ac.ed.lurg.country.GlobalPrice;
import ac.ed.lurg.country.SubsidyRateManager;
import ac.ed.lurg.country.TradeManager;
import ac.ed.lurg.demand.AbstractDemandManager;
import ac.ed.lurg.demand.BaseConsumpManager;
import ac.ed.lurg.demand.DemandManagerFromFile;
import ac.ed.lurg.demand.DemandManagerSSP;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.CropUsageReader;
import ac.ed.lurg.landuse.FPUManager;
import ac.ed.lurg.landuse.IrrigationConstraintReader;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.IrrigationMaxAmountReader;
import ac.ed.lurg.landuse.IrrigationRasterSet;
import ac.ed.lurg.landuse.IrrigiationCostReader;
import ac.ed.lurg.landuse.LandCoverItem;
import ac.ed.lurg.landuse.LandCoverReader;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.MaxCropAreaReader;
import ac.ed.lurg.landuse.ProtectedAreasReader;
import ac.ed.lurg.landuse.RunOffReader;
import ac.ed.lurg.output.LandUseOutputer;
import ac.ed.lurg.output.LpjgOutputer;
import ac.ed.lurg.shock.YieldShockReader;
import ac.ed.lurg.shock.YieldShocks;
import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.LPJYieldResponseMapReader;
import ac.ed.lurg.yield.YieldRaster;
import ac.ed.lurg.yield.YieldResponsesItem;
import ac.sac.raster.IntegerRasterItem;
import ac.sac.raster.IntegerRasterReader;
import ac.sac.raster.InterpolatingRasterSet;
import ac.sac.raster.RasterHeaderDetails;
import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterOutputer;
import ac.sac.raster.RasterSet;

public class ModelMain {

	private Collection<CountryAgent> countryAgents;
	private CountryBoundaryRaster countryBoundaryRaster;
	private AbstractDemandManager demandManager;
	private SubsidyRateManager subsidyRateManager;
	private TradeManager tradeManager;
	private CompositeCountryManager compositeCountryManager;
	private RasterHeaderDetails desiredProjection;

	private InternationalMarket internationalMarket;
	private IrrigationRasterSet currentIrrigationData;
	private RasterSet<LandUseItem> globalLandUseRaster;
	private RasterSet<IntegerRasterItem> clusterIdRaster;

	public static void main(String[] args) {
		ModelMain theModel = new ModelMain();
	    	System.out.println("Working Directory = " + System.getProperty("user.dir"));
		theModel.setup();
		theModel.run();
	}

	/* setup models, reading inputs, etc. */
	private void setup() {
		desiredProjection = RasterHeaderDetails.getGlobalHeaderFromCellSize(ModelConfig.CELL_SIZE_X, ModelConfig.CELL_SIZE_Y, "999");

		BaseConsumpManager baseConsumpManager = new BaseConsumpManager();
		compositeCountryManager = new CompositeCountryManager(baseConsumpManager);
		subsidyRateManager = new SubsidyRateManager(compositeCountryManager);
		
		if (ModelConfig.DEMAND_FROM_FILE)
			demandManager = new DemandManagerFromFile(compositeCountryManager);
		else
			demandManager = new DemandManagerSSP(ModelConfig.SSP_SCENARIO, baseConsumpManager, compositeCountryManager);
		
		if (ModelConfig.SHOCKS_POSSIBLE) ModelConfig.readInShocksFile();
			
		tradeManager = new TradeManager(compositeCountryManager);
		currentIrrigationData = getFixedIrrigationData();
		countryBoundaryRaster = getCountryBoundaryRaster();
		clusterIdRaster = ModelConfig.GENERATE_NEW_YIELD_CLUSTERS ? new RasterSet<IntegerRasterItem>(desiredProjection) : getClusterRaster();

		countryAgents = createCountryAgents(compositeCountryManager.getAll());
		globalLandUseRaster = new RasterSet<LandUseItem>(desiredProjection);
		internationalMarket = new InternationalMarket();
	}

	/* run the model */
	private void run() {
		for (int i = ModelConfig.START_TIMESTEP; i <= ModelConfig.END_TIMESTEP; i++) {
			Timestep timestep = new Timestep(i);
			try {
				doTimestep(timestep);
			} catch (Exception e) {
				LpjgOutputer.writeMarkerFile(timestep.getYear(), true);
				throw new RuntimeException(e);
			}
		}
	}

	private void doTimestep(Timestep timestep) {
		LogWriter.println(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

		YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-90.5, 45.5);
		LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.MAIZE) + ", " + yresp.getYieldFertOnly(CropType.MAIZE) + ", "
						+ yresp.getYieldIrrigOnly(CropType.MAIZE) + ", " + yresp.getYieldNone(CropType.MAIZE));

		double previousGen2EcDDemand = (timestep.isInitialTimestep() || ModelConfig.IS_CALIBRATION_RUN ) ? 0: demandManager.getSecondGenBioenergyDemand(timestep.getPreviousTimestep());
		double gen2EcDDemand = demandManager.getSecondGenBioenergyDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
		double gen2Increase = (gen2EcDDemand>previousGen2EcDDemand) ? gen2EcDDemand - previousGen2EcDDemand : 0.0;

		for (CountryAgent ca : countryAgents) {

			LogWriter.println("Country " + ca.getCountry());
			Collection<RasterKey> countryKeys = countryBoundaryRaster.getKeysFor(ca.getCountry());
			YieldRaster countryYieldSurfaces = yieldSurfaces.createSubsetForKeys(countryKeys);
			RasterSet<IrrigationItem> irrigData = currentIrrigationData.createSubsetForKeys(countryKeys);
			
			// do the optimization
			try {
				ca.determineProduction(timestep, countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase);
			} catch (Exception e) {
				LogWriter.printlnError("Exception processing " + ca.getCountry() + " will continue with other countries");
				LogWriter.print(e);
				continue;
			}

			// some hacky code for debug purposes that keeps each gams gdx file
			// for one country
			if (ModelConfig.GAMS_COUNTRY_TO_SAVE != null && ca.getCountry().getName().equals(ModelConfig.GAMS_COUNTRY_TO_SAVE)) {
				try {
					Files.copy(
							FileSystems.getDefault().getPath(ModelConfig.TEMP_DIR + File.separator + "_gams_java_gdb1.gdx"),
							FileSystems.getDefault().getPath(ModelConfig.TEMP_DIR + File.separator + ModelConfig.GAMS_COUNTRY_TO_SAVE + timestep.getYear() + ".gdx"),
							StandardCopyOption.REPLACE_EXISTING);
				} catch (IOException e) {
					LogWriter.print(e);
				}
			}

			// update global rasters
			globalLandUseRaster.putAll(ca.getLandUses());

			// if first timestep and calibration get the clustering info, which
			// doesn't change through time
			if (ModelConfig.GENERATE_NEW_YIELD_CLUSTERS && timestep.isInitialTimestep())
				clusterIdRaster.putAll(ca.getYieldClusters());
		}
		
		internationalMarket.determineInternationalTrade(countryAgents, gen2EcDDemand, timestep);

		// output results
		outputTimestepResults(timestep, globalLandUseRaster);
	}

	private BufferedWriter getFileWriter(Timestep timestep, String file, String columnHeadings) throws IOException {
		if (timestep.isInitialTimestep()) {
			FileWriter fstream = new FileWriter(file, false);
			BufferedWriter outputFile = new BufferedWriter(fstream);
			outputFile.write(columnHeadings);
			outputFile.newLine();
			return outputFile;
		} else {
			FileWriter fstream = new FileWriter(file, true);
			return new BufferedWriter(fstream);
		}
	}

	private void writeLandCoverFile(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year,Cropland,Pasture,ManForest,UnmanForest,Natural,AbPasture,EnergyCrop,FertCrop,IrrigCrop");
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.LAND_COVER_OUTPUT_FILE, sbHeadings.toString());

			StringBuffer sbData = new StringBuffer();
			sbData.append(String.format("%d,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f", timestep.getYear(),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.CROPLAND),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.PASTURE),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.MANAGED_FOREST),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.UNMANAGED_FOREST),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.OTHER_NATURAL),
					LandUseItem.getAbandonedPasture(landUseRaster.values()))
					);
						
			sbData.append(String.format(",%.1f", LandUseItem.getTotalCropArea(landUseRaster.values(), CropType.ENERGY_CROPS)));
			sbData.append(String.format(",%.1f", LandUseItem.getFertiliserTotal(landUseRaster.values(), CropType.getCropsLessPasture()) / 1000));
			sbData.append(String.format(",%.1f", LandUseItem.getIrrigationTotal(landUseRaster.values(), CropType.getCropsLessPasture())));

			outputFile.write(sbData.toString());
			outputFile.newLine();
			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}

	private void writeGlobalMarketFile(Timestep timestep) {
		try {
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.PRICES_OUTPUT_FILE, "Year,Crop,Imports (Mt),Exports (Mt),New export price, Stock Levels (Mt)");
			internationalMarket.writeGlobalMarketFile(timestep, outputFile);
			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}
	
	private void writeGlobalFoodBalanceSheet(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {
		try {
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.FOOD_BALANCE_SHEET_FILE, "Year,Crop,Production,Imports,Export,TransportLosses,StockVar,Supply,MonogastricsFeed,RuminantsFeed,SeedAndOtherLosses,FoodAnd1stGen,ProdArea");
			Map<CropType, GlobalPrice> worldPrices = internationalMarket.getWorldPrices();
			double harvestedAndFallowArea = 0;
			
			for (CropType crop : CropType.getCropsLessPasture()) {
				GlobalPrice priceQuantity = worldPrices.get(crop); // some specific logic for import/exports and this has been aggregated already, so best to use it
				double prod=0, prodArea=0, feedMonogastrics=0, feedRuminants=0;
				double exportsBeforeTL=0, imports=0, transportloss=0, stockChange=0;
				
				if (priceQuantity != null) {
					exportsBeforeTL = priceQuantity.getExportsBeforeTransportLoss();
					imports = priceQuantity.getImportAmount();
					transportloss = priceQuantity.getTransportLosses();
					stockChange = priceQuantity.getStockChange();
				}
				
				for (CountryAgent ca : countryAgents) {
					Map<CropType, CropUsageData> allCropUsage = ca.getCropUsageData();
					CropUsageData cropUsage = allCropUsage.get(crop);
					if (cropUsage != null) {
						prod += cropUsage.getProduction();
						prodArea += cropUsage.getArea();
						feedMonogastrics += cropUsage.getMonogastricFeed();
						feedRuminants += cropUsage.getRuminantFeed();
					}
				}

				double seedAndWaste = prod * crop.getSeedAndWasteRate();
				double netSupply = prod - exportsBeforeTL + imports;
				double foodAnd1stGen = netSupply - feedMonogastrics - feedRuminants - seedAndWaste;
				
				if (!crop.equals(CropType.SETASIDE))
					prodArea *= (1-ModelConfig.UNHANDLED_CROP_RATE);  // remove unhandled crop area
				
				harvestedAndFallowArea += prodArea;
				
				StringBuffer sbData = new StringBuffer();
				sbData.append(String.format("%d,%s", timestep.getYear(), crop.getGamsName()));
				sbData.append(String.format(",%.2f", prod));
				sbData.append(String.format(",%.2f,%.2f,%.2f,%.2f", imports, exportsBeforeTL, transportloss, stockChange));
				sbData.append(String.format(",%.2f,%.2f,%.2f,%.2f,%.2f", netSupply, feedMonogastrics, feedRuminants, seedAndWaste, foodAnd1stGen));
				sbData.append(String.format(",%.2f", prodArea));
				outputFile.write(sbData.toString());
				outputFile.newLine();
			}

			double cropLandArea = LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.CROPLAND);
			double unhandledArea = cropLandArea - harvestedAndFallowArea;
			outputFile.write(String.format("%d,%s,,,,,,,,,,,%.2f", timestep.getYear(), "unhandled", unhandledArea));
			outputFile.newLine();

			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}


	private void writeDemandFile(Timestep timestep) {
		try {
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.DEMAND_OUTPUT_FILE, "Year,Commodity,Amount (Mt)");

			for (CommodityType comm : CommodityType.getAllItems()) {
				double demandAmount = 0;

				for (CountryAgent country : countryAgents) {
					Double d = country.getCurrentProjectedDemand().get(comm);
					if (d != null) {
						demandAmount += d.doubleValue();
						LogWriter.println(String.format("%s,%s,%.4f", country.getCountry(), comm.getGamsName(), d));
					}
				}
				StringBuffer sbData = new StringBuffer();
				sbData.append(String.format("%d,%s", timestep.getYear(), comm.getGamsName()));
				sbData.append(String.format(",%.1f", demandAmount));

				LogWriter.println("Global demand " + timestep.getYear() + " " + comm.getGamsName() + " " + demandAmount + "\n");
				outputFile.write(sbData.toString());
				outputFile.newLine();
			}
			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}
	
private void writeDomesticProductionFile(Timestep timestep) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year, Country, Crop, Area, Production, Production_cost, Import_price, Export_price, Net_imports, Rum_feed_amount, Mon_feed_amount");
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.DOMESTIC_OUTPUT_FILE, sbHeadings.toString());

			for (CropType crop : CropType.getAllItems()) {
				for (CountryAgent country : countryAgents) {
					
					Map<CropType, CropUsageData> cropUsageAllCrops = country.getCropUsageData();
					CropUsageData cropUsage = cropUsageAllCrops.get(crop);
					
					if (cropUsage == null) 
						continue;
					
					Double prodCosts = cropUsage.getProdCost();
					Double prod = cropUsage.getProduction();
					Double area = cropUsage.getArea();
					Double rumFeedAmount = cropUsage.getRuminantFeed();
					Double monFeedAmount = cropUsage.getMonogastricFeed();
					
					Double importPrice = null;
					Double exportPrice = null;
					Double netImports = null;

					if (crop.isImportedCrop()) {
						CountryPrice px = country.getCurrentCountryPrices().get(crop);
						importPrice = px.getImportPrice(); 
						exportPrice = px.getExportPrice();
						netImports = cropUsage.getNetImports();	//this isn't accounting for transport losses in exports 
					}
	
					StringBuffer sbData = new StringBuffer();
					sbData.append(String.format("%d,%s,%s", timestep.getYear(), country.getCountry(), crop.getGamsName()));
					sbData.append(String.format(",%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f", area, prod, prodCosts, importPrice, exportPrice, netImports, rumFeedAmount, monFeedAmount));

					outputFile.write(sbData.toString());
					outputFile.newLine();
				}
			}
			outputFile.close();

		} catch (IOException e) {
			LogWriter.print(e);
		}
	}

	private void outputTimestepResults(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {

		writeLandCoverFile(timestep, landUseRaster);
		writeGlobalMarketFile(timestep);
		writeDemandFile(timestep);
		writeDomesticProductionFile(timestep);
		writeGlobalFoodBalanceSheet(timestep, landUseRaster);

		if (ModelConfig.OUTPUT_FOR_LPJG) {
			for (int outputYear : timestep.getYearsFromLast()) {
				LogWriter.printlnError("Outputing Year: " + outputYear);
				RasterSet<LandUseItem> landUseToOutput = null;

				if (outputYear == timestep.getYear()) {
					landUseToOutput = landUseRaster;
				} 
				else if (ModelConfig.INTERPOLATE_OUTPUT_YEARS) {
					InterpolatingRasterSet<LandUseItem> intermediateLandUse = new InterpolatingRasterSet<LandUseItem>( landUseRaster.getHeaderDetails()) {
						private static final long serialVersionUID = 1306045141011047760L;
						protected LandUseItem createRasterData() {
							return new LandUseItem();
						}
					};
					intermediateLandUse.setup(globalLandUseRaster, landUseRaster, timestep.getPreviousTimestep().getYear(), timestep.getYear(), outputYear);
					landUseToOutput = intermediateLandUse;
				}

				if (landUseToOutput != null) {
					LpjgOutputer lpjOutputer = new LpjgOutputer(outputYear, landUseToOutput);
					lpjOutputer.writeOutput();
				}
			}
			outputWaterAvailablity(timestep, currentIrrigationData);  // uses the year directory structure created above
		}

		if (ModelConfig.IS_CALIBRATION_RUN)
			serializeLandUse(landUseRaster);

		if (timestep.isInitialTimestep() && ModelConfig.GENERATE_NEW_YIELD_CLUSTERS)
				outputClusters(clusterIdRaster);

		// Output LandUses to tabular file, for analysis (perhaps)
		LogWriter.println("Outputing land uses Year: " + timestep.getYear());
		LandUseOutputer landuseOutputer = new LandUseOutputer(timestep.getYear(), landUseRaster);
		landuseOutputer.writeOutput();

		// don't really need this a LPJ outputs have same data, although in a slightly different format
		// outputLandCover(timestep.getYear(), landUseRaster, LandCoverType.CROPLAND);
		// outputLandCover(timestep.getYear(), landUseRaster, LandCoverType.PASTURE);
	}

	private void outputWaterAvailablity(Timestep timestep, IrrigationRasterSet irrigiationRS) {
		new RasterOutputer<Double, IrrigationItem>(irrigiationRS, ModelConfig.OUTPUT_DIR + File.separator + timestep.getYear() + File.separator + "IrrigConstraint.asc") {
			@Override
			public Double getValue(RasterKey location) {
				IrrigationItem item = results.get(location);
				if (item == null)
					return null;

				return item.getIrrigConstraint();
			}
		}.writeOutput();
	}

	private void outputClusters(RasterSet<IntegerRasterItem> landUseRaster) {
		new RasterOutputer<Integer, IntegerRasterItem>(landUseRaster, ModelConfig.CLUSTERED_YIELD_FILE) {
			@Override
			public Integer getValue(RasterKey location) {
				IntegerRasterItem item = results.get(location);
				if (item == null)
					return null;

				return item.getInt();
			}
		}.writeOutput();
	}

	public RasterSet<IntegerRasterItem> getClusterRaster() {
		RasterSet<IntegerRasterItem> clusters = new RasterSet<IntegerRasterItem>(desiredProjection) {
			private static final long serialVersionUID = 2467452274591854417L;

			@Override
			protected IntegerRasterItem createRasterData() {
				return new IntegerRasterItem(0);
			}
		};

		IntegerRasterReader clusterReader = new IntegerRasterReader(clusters);
		clusterReader.getRasterDataFromFile(ModelConfig.CLUSTERED_YIELD_FILE);
		return clusters;
	}

	public CountryBoundaryRaster getCountryBoundaryRaster() {
		CountryBoundaryRaster countryBoundaries = new CountryBoundaryRaster(desiredProjection, compositeCountryManager);
		CountryBoundaryReader countryReader = new CountryBoundaryReader(countryBoundaries);
		countryReader.getRasterDataFromFile(ModelConfig.COUNTRY_BOUNDARY_FILE);
		return countryBoundaries;
	}

	public Collection<CountryAgent> createCountryAgents(Collection<CompositeCountry> countryGrouping) {
		Collection<CountryAgent> countryAgents = new HashSet<CountryAgent>();
		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap = new CropUsageReader(compositeCountryManager).getCommodityData();
		RasterSet<LandUseItem> initLU = getInitialLandUse();

		for (CompositeCountry cc : countryGrouping) {

			// DEBUG code
			if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {
				if (!(cc.getName().equals(ModelConfig.DEBUG_COUNTRY_NAME)))
					continue;
			}

			List<RasterKey> keys = countryBoundaryRaster.getKeysFor(cc);
			Map<CropType, CropUsageData> countryCommodityData = cropUsageDataMap.get(cc);
			Map<CropType, Double> countryTradeBarriers = tradeManager.getTradeBarriers(cc);

			if (countryCommodityData == null) {
				LogWriter.printlnError("No commodities data for " + cc + ", so skipping");
			} else {
				RasterSet<LandUseItem> initCountryLandUse = initLU.createSubsetForKeys(keys);
				RasterSet<IntegerRasterItem> yieldClusters = ModelConfig.GENERATE_NEW_YIELD_CLUSTERS ? null : clusterIdRaster.createSubsetForKeys(keys);

				if (initCountryLandUse.size() == 0) {
					LogWriter.printlnError("No initial land use for " + cc + ", so skipping");
					continue;
				}
				
				Map<CropType, Double> subsidyRates = subsidyRateManager.getSubsidyRates(cc);

				CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLandUse, countryCommodityData, countryTradeBarriers, yieldClusters, subsidyRates);
				countryAgents.add(ca);
				LogWriter.println("Creating country agent for: " + cc);
			}
		}

		return countryAgents;
	}

	private RasterSet<LandUseItem> getInitialLandUse() {
		if (ModelConfig.IS_CALIBRATION_RUN)
			return getLandUseFromBaseline();
		else
			return deserializeLandUse();
	}

	private void serializeLandUse(RasterSet<LandUseItem> landUseRaster) {
		try {
			LogWriter.println("Starting serializing LandUse to " + ModelConfig.SERIALIZED_LAND_USE_FILE);
			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_LAND_USE_FILE);
			ObjectOutputStream out = new ObjectOutputStream(fileOut);
			out.writeObject(landUseRaster);
			out.close();
			fileOut.close();
			LogWriter.println("Serialized data is saved");
		} catch (IOException i) {
			i.printStackTrace();
		}
	}

	@SuppressWarnings("unchecked")
	private RasterSet<LandUseItem> deserializeLandUse() {
		try {
			RasterSet<LandUseItem> initLU;
			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_LAND_USE_FILE);
			ObjectInputStream in = new ObjectInputStream(fileIn);
			initLU = (RasterSet<LandUseItem>) in.readObject();
			in.close();
			fileIn.close();
			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_LAND_USE_FILE);
			return initLU;
		} catch (IOException i) {
			LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_LAND_USE_FILE);
			LogWriter.print(i);
			return null;
		} catch (ClassNotFoundException c) {
			LogWriter.printlnError("RasterSet<LandUseItem> class not found");
			c.printStackTrace();
			return null;
		}
	}

	/** this is if we are starting from Hurtt of other initial land covers (so we don't have land uses and intensity data) */
	private RasterSet<LandUseItem> getLandUseFromBaseline() {
		RasterSet<LandCoverItem> initialLC = new RasterSet<LandCoverItem>(desiredProjection) {
			private static final long serialVersionUID = 4642550777741425501L;

			protected LandCoverItem createRasterData() {
				return new LandCoverItem();
			}
		};

		new ProtectedAreasReader(initialLC).getRasterDataFromFile(ModelConfig.PROTECTED_AREAS_FILE);
		new MaxCropAreaReader(initialLC).getRasterDataFromFile(ModelConfig.HIGH_SLOPE_AREAS_FILE);
		new LandCoverReader(initialLC).getRasterDataFromFile(ModelConfig.INITAL_LAND_COVER_FILE);

		RasterSet<LandUseItem> landUseRaster = new RasterSet<LandUseItem>(initialLC.getHeaderDetails());

		for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet()) {
			LandUseItem areasItem = new LandUseItem();
			RasterKey key = entry.getKey();
			LandCoverItem landCover = entry.getValue();
			areasItem.setLandCoverAreas(landCover);
			areasItem.setCropFraction(CropType.WHEAT, 0.5); // random start as don't have better data
			areasItem.setCropFraction(CropType.MAIZE, 0.5);
			areasItem.setProtectedArea(landCover.getProtectedArea());
			areasItem.setUnavailableFract(landCover.getUnavailableFract());
			landUseRaster.put(key, areasItem);
		}

		return landUseRaster;
	}

	private YieldRaster getYieldSurfaces(Timestep timestep) {
		LPJYieldResponseMapReader yieldReader = new LPJYieldResponseMapReader(desiredProjection);
		YieldRaster yieldRaster = yieldReader.getRasterData(timestep);
		
		YieldResponsesItem yresp = yieldRaster.getFromCoordinates(-90.5, 45.5);
		LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.MAIZE) + ", " + yresp.getYieldFertOnly(CropType.MAIZE) + ", "
						+ yresp.getYieldIrrigOnly(CropType.MAIZE) + ", " + yresp.getYieldNone(CropType.MAIZE));

		// read in shocks for timestep
		// for each, get shock map and apply to yieldRaster
		if (false){
			double yieldAdjustment = 0.6;  // these would come from shock descriptor file
			CropType crop = CropType.MAIZE;
			
			YieldShocks yieldShocks = new YieldShocks(desiredProjection);
			YieldShockReader yieldShockReader = new YieldShockReader(yieldShocks);
			yieldShockReader.getRasterDataFromFile(ModelConfig.SPATIAL_DATA_DIR + File.separator + "yieldshockmap.txt");
			yieldShocks.adjustForShocks(yieldRaster, crop, yieldAdjustment);
		}
		
		yresp = yieldRaster.getFromCoordinates(-90.5, 45.5);
		LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.MAIZE) + ", " + yresp.getYieldFertOnly(CropType.MAIZE) + ", "
						+ yresp.getYieldIrrigOnly(CropType.MAIZE) + ", " + yresp.getYieldNone(CropType.MAIZE));

		return yieldRaster;
	}

	/** Get irrigation data that does not change with time, should only be called once */
	private IrrigationRasterSet getFixedIrrigationData() {
		IrrigationRasterSet fixedIrrigData = new IrrigationRasterSet(desiredProjection, new FPUManager(desiredProjection));
		new IrrigiationCostReader(fixedIrrigData).getRasterDataFromFile(ModelConfig.IRRIGATION_COST_FILE);
		new IrrigationConstraintReader(fixedIrrigData).getRasterDataFromFile(ModelConfig.IRRIGATION_CONSTRAINT_FILE);
		
		String baseTimestepRootDir = Timestep.getYearSubDir(ModelConfig.YIELD_DIR, ModelConfig.ELLIOTT_BASEYEAR); // needs to be Elliott base timestep
		new RunOffReader(fixedIrrigData, true).getRasterDataFromFile(baseTimestepRootDir + File.separator + ModelConfig.IRRIG_RUNOFF_FILE);
		
		fixedIrrigData.calcIrrigConstraintOffsets();  // should have everything we need to calc offset between Elliott and LPJ data
		return fixedIrrigData;
	}

	/** Ugly in situ update of currentIrrigationData, better if IrrigationRasterSets were handled more immutably */
	private void getUpdateIrrigationData(Timestep timestep, YieldRaster yieldSurfaces) {
		String rootDir = timestep.getYearSubDir(ModelConfig.YIELD_DIR);
		
		new IrrigationMaxAmountReader(currentIrrigationData).getRasterDataFromFile(rootDir + File.separator + ModelConfig.IRRIG_MAX_WATER_FILENAME);

		if (!ModelConfig.USE_BLUE_WATER_FILE_IRRIG_CONSTRAINT) {
			new RunOffReader(currentIrrigationData, false).getRasterDataFromFile(rootDir + File.separator + ModelConfig.IRRIG_RUNOFF_FILE);
			currentIrrigationData.updateIrrigConstraints(timestep);
		}
	}
}