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.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;

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.GlobalPrice;
import ac.ed.lurg.country.StockReader;
import ac.ed.lurg.country.TradeManager;
import ac.ed.lurg.country.gams.GamsRasterOutput;
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.types.CommodityType;
import ac.ed.lurg.types.CropToDoubleMap;
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 TradeManager tradeManager;
	private FPUManager fpuManager;
	private CompositeCountryManager compositeCountryManager;
	private RasterHeaderDetails desiredProjection;

	private Map<CropType, GlobalPrice> prevWorldPrices;
	private Map<CropType, Double> prevStockLevel;
	private RasterSet<LandUseItem> prevLandUseRaster;
	private IrrigationRasterSet currentIrrigationData;
	private RasterSet<LandUseItem> globalLandUseRaster;
	private RasterSet<IntegerRasterItem> clusterIdRaster;

	public static void main(String[] args)  {
		ModelMain theModel = new ModelMain();
		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);
		
		if (ModelConfig.DEMAND_FROM_FILE)
			demandManager = new DemandManagerFromFile(compositeCountryManager);
		else
			demandManager = new DemandManagerSSP(ModelConfig.SSP_SCENARIO, baseConsumpManager, compositeCountryManager);
		
		tradeManager = new TradeManager(compositeCountryManager);
		fpuManager = new FPUManager(desiredProjection);
		currentIrrigationData = getFixedIrrigationData();
		countryBoundaryRaster = getCountryBoundaryRaster();
		clusterIdRaster = (ModelConfig.IS_CALIBRATION_RUN) ? new RasterSet<IntegerRasterItem>(desiredProjection) : getClusterRaster();

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

		// in first timestep we don't have this info, but ok as constrained to import/export specified amount, values based on http://www.indexmundi.com/commodities/ for Jun 2010
		prevWorldPrices = new HashMap<CropType, GlobalPrice>();  
		prevWorldPrices.put(CropType.WHEAT, GlobalPrice.createInitial(0.157 * ModelConfig.INITIAL_PRICE_SHIFT));
		prevWorldPrices.put(CropType.MAIZE, GlobalPrice.createInitial(0.152 * ModelConfig.INITIAL_PRICE_SHIFT));
		prevWorldPrices.put(CropType.RICE, GlobalPrice.createInitial(0.182 * ModelConfig.INITIAL_PRICE_SHIFT));
		prevWorldPrices.put(CropType.OILCROPS, GlobalPrice.createInitial((0.820 * .4 + 0.314 * .6)*0.5 * ModelConfig.INITIAL_PRICE_SHIFT));

		prevWorldPrices.put(CropType.PULSES, GlobalPrice.createInitial(0.4 * ModelConfig.INITIAL_PRICE_SHIFT));
		prevWorldPrices.put(CropType.STARCHY_ROOTS, GlobalPrice.createInitial(0.1 * ModelConfig.INITIAL_PRICE_SHIFT));
		prevWorldPrices.put(CropType.MONOGASTRICS, GlobalPrice.createInitial(0.4 * 0.5 * ModelConfig.INITIAL_PRICE_SHIFT));  // quantities is in feed equivalent term (0.4 is weighted average price per feed, and 0.5 accounts for mark-up for additional processing)
		prevWorldPrices.put(CropType.RUMINANTS, GlobalPrice.createInitial(0.1 * 0.6 * ModelConfig.INITIAL_PRICE_SHIFT));  // quantities is in feed equivalent term
		prevWorldPrices.put(CropType.ENERGY_CROPS, GlobalPrice.createInitial(0.04 * ModelConfig.INITIAL_PRICE_SHIFT));
		prevStockLevel = getInitialStockLevels();
	}

	/* 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
		RasterSet<IrrigationItem> allIrrigData = getUpdateIrrigationData(timestep, yieldSurfaces);

		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() ? 0 : demandManager.getSecondGenBioenergyDemand(timestep.getPreviousTimestep());
		double gen2EcDDemand = demandManager.getSecondGenBioenergyDemand(timestep);
		CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
		CropToDoubleMap totalExportCommodities = new CropToDoubleMap();

		for (CountryAgent ca : countryAgents) {

			LogWriter.println("Country " + ca.getCountry());
			Collection<RasterKey> countryKeys = countryBoundaryRaster.getKeysFor(ca.getCountry());
			YieldRaster countryYieldSurfaces = yieldSurfaces.createSubsetForKeys(countryKeys);
			RasterSet<IrrigationItem> irrigData = allIrrigData.createSubsetForKeys(countryKeys);

			GamsRasterOutput result = null;

			// do the optimization
			try {
				result = ca.determineProduction(timestep, countryYieldSurfaces, irrigData, prevWorldPrices, (gen2EcDDemand-previousGen2EcDDemand));
			}
			catch (Exception e) {
				LogWriter.printlnError("Exception processing " + ca.getCountry() + " will continue with other countries");	
				LogWriter.print(e);
			}

			if (result == null) {
				LogWriter.printlnError("No results for " + ca.getCountry());
				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(result.getLandUses());

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

			// Get values for world input costs
			Map<CropType, CropUsageData> cropUsage = result.getCropUsageData();

			for (Entry<CropType, CropUsageData> entry : cropUsage.entrySet()) {
				CropType c = entry.getKey();
				double countryNetImports = entry.getValue().getNetImports();

				if (countryNetImports > 0)
					totalImportCommodities.incrementValue(c, countryNetImports);
				else
					totalExportCommodities.incrementValue(c, -countryNetImports);
			}
		}
		
		// energycrops are a special case where demand in global and exogenously specified
		totalImportCommodities.incrementValue(CropType.ENERGY_CROPS, gen2EcDDemand);
		
		// Look at trade balance and adjust appropriately
		for (CropType crop : CropType.getImportedTypes()) {
			GlobalPrice prevPrice = prevWorldPrices.get(crop);
			double imports = totalImportCommodities.containsKey(crop) ? totalImportCommodities.get(crop) : 0.0;
			double exports = totalExportCommodities.containsKey(crop) ? totalExportCommodities.get(crop) * (1.0-ModelConfig.TRANSPORT_LOSSES): 0.0;
			
			double previousStockLevel = prevStockLevel.get(crop);
			GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(imports, exports, (ModelConfig.MARKET_ADJ_PRICE));
			LogWriter.println(String.format("Price for %s updated from %s (imports amount %.0f, exports amount %.0f) to %s ", crop.getGamsName(), prevPrice, imports, exports, adjustedPrice));
			prevWorldPrices.put(crop, adjustedPrice);
			double updatedStockLevel = previousStockLevel + exports-imports;
			prevStockLevel.put(crop, updatedStockLevel);
			LogWriter.println(String.format("Global stocks for %s updated from %s to %s ", crop.getGamsName(),previousStockLevel, updatedStockLevel));
			if(updatedStockLevel < 0) LogWriter.println(String.format("Global stocks for %s is below zero ", crop.getGamsName(),previousStockLevel));
		}


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

		// keep last to allow interpolation
		prevLandUseRaster = 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,EnergyCrop,FertCrop,IrrigCrop"); //Year,Cropland (Mha),Pasture (Mha),Managed_Forest (Mha),Unmanaged_Forest (Mha),Natural (Mha),Energycrops (Mha),Fert crop (Mt),,Irrig crop (km3)");
			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", 
					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)));

			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.getFertiliserTotal(landUseRaster.values(), CropType.PASTURE)/1000));
			sbData.append(String.format(",%.1f", LandUseItem.getIrrigationTotal(landUseRaster.values(), CropType.getCropsLessPasture())));
	//		sbData.append(String.format(",%.1f", LandUseItem.getIrrigationTotal(landUseRaster.values(), CropType.PASTURE)));

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

	private void writeGlobalMarketFile(Timestep timestep) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year,Crop,Imports (Mt),Exports (Mt),New export price, Stock Levels (Mt)");
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.PRICES_OUTPUT_FILE, sbHeadings.toString());

			for (CropType crop : CropType.getImportedTypes() ) {
				StringBuffer sbData = new StringBuffer();
				sbData.append(String.format("%d,%s", timestep.getYear(), crop.getGamsName()));
				sbData.append(String.format(",%.1f,%.1f", prevWorldPrices.get(crop).getImportAmount(), prevWorldPrices.get(crop).getExportAmount()));
				sbData.append(String.format(",%.3f,%.3f", prevWorldPrices.get(crop).getExportPrice(), prevStockLevel.get(crop)));


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

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

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

			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 outputTimestepResults(Timestep timestep, RasterSet<LandUseItem> landUseRaster, YieldRaster yieldSurfaces) {

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

		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(prevLandUseRaster, landUseRaster, timestep.getPreviousTimestep().getYear(), timestep.getYear(), outputYear);
					landUseToOutput = intermediateLandUse;
				}

				if (landUseToOutput != null) {
					LpjgOutputer lpjOutputer = new LpjgOutputer(outputYear, landUseToOutput, yieldSurfaces);
					lpjOutputer.writeOutput();
				}
			}
		}

		if (ModelConfig.IS_CALIBRATION_RUN) {
			serializeLandUse(landUseRaster);

			if (timestep.isInitialTimestep()) {
				outputClusters(clusterIdRaster);
				outputWaterAvailablity(currentIrrigationData);
			}
		}

		// 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(IrrigationRasterSet irrigiationRS) {
		new RasterOutputer<Double, IrrigationItem>(irrigiationRS, ModelConfig.OUTPUT_DIR + File.separator + "irrig_constraint.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(); // read in all land use data saved from previous state, e.g. calibrations run

		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.IS_CALIBRATION_RUN) ? null : clusterIdRaster.createSubsetForKeys(keys);

				if (initCountryLandUse.size() == 0) {
					LogWriter.printlnError("No initial land use for " +cc + ", so skipping");
					continue;
				}

				CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLandUse, countryCommodityData, countryTradeBarriers, yieldClusters);
				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;
		}

		// alternative to deserialization is to use tabular data
		/*	RasterSet<LandUseItem> initLU = new RasterSet<LandUseItem>(desiredProjection) {
			private static final long serialVersionUID = 1317102740312961042L;
			protected LandUseItem createRasterData() {
				return new LandUseItem();
			}
		};
		new LandUseReader(initLU).getRasterDataFromFile(ModelConfig.INITAL_LAND_USE_FILE);
		return initLU; */
	} 

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

		/*	new RasterOutputer<LandCoverItem>(initLC, "InitialCropland") {
			@Override
			public Double getValue(RasterKey location) {
				LandCoverItem item = results.get(location);
				if (item == null)
					return null;
				return item.getLandCoverFract(LandCoverType.CROPLAND);
			}
		}.writeOutput(); */

		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, better if we could get data, but free substitution between crops so not critical
			areasItem.setCropFraction(CropType.MAIZE, 0.5);
			areasItem.setProtectedArea(landCover.getProtectedArea());
			landUseRaster.put(key, areasItem);
		}

		return landUseRaster;
	} 

	private YieldRaster getYieldSurfaces(Timestep timestep) {
		LPJYieldResponseMapReader yieldReader = new LPJYieldResponseMapReader(desiredProjection);
		return yieldReader.getRasterData(timestep); 
	}
	

	private IrrigationRasterSet getFixedIrrigationData() {	
	    IrrigationRasterSet fixedIrrigData = new IrrigationRasterSet(desiredProjection, fpuManager);
		new IrrigiationCostReader(fixedIrrigData).getRasterDataFromFile(ModelConfig.IRRIGATION_COST_FILE);
        
		if (ModelConfig.USE_BLUE_WATER_FILE_IRRIG_CONSTRAINT)
			new IrrigationConstraintReader(fixedIrrigData).getRasterDataFromFile(ModelConfig.IRRIGATION_CONSTRAINT_FILE);
		
        return fixedIrrigData;
	}
	
	private RasterSet<IrrigationItem> getUpdateIrrigationData(Timestep timestep, YieldRaster yieldSurfaces) {		
		String rootDir = timestep.getYearSubDir(ModelConfig.YIELD_DIR);
		new IrrigationMaxAmountReader(currentIrrigationData, yieldSurfaces).getRasterDataFromFile(rootDir + File.separator + ModelConfig.IRRIG_MAX_WATER_FILENAME);
		
		if (!ModelConfig.USE_BLUE_WATER_FILE_IRRIG_CONSTRAINT) {
			new RunOffReader(currentIrrigationData).getRasterDataFromFile(rootDir + File.separator + ModelConfig.IRRIG_RUNOFF_FILE);
			currentIrrigationData.updateConstraintByFPU(timestep);
		}
		return currentIrrigationData;
	}
	
	private Map<CropType,Double> getInitialStockLevels(){
		Map<CropType, Double> initialStockLevels = new StockReader().read();
		initialStockLevels.put(CropType.ENERGY_CROPS, 0.0);
		return initialStockLevels	;
	}
}