Support non-square grid cells

src/ac/ed/lurg/ModelConfig.java

@@ -101,8 +101,8 @@ public class ModelConfig {
 	// Spatial (gridded) data
 	public static final String SPATIAL_DATA_DIR = DATA_DIR; // + File.separator + "tinyTest";
 	public static final String INITAL_LAND_COVER_DIR = SPATIAL_DATA_DIR + File.separator + "initLandUse";
-	public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries.asc";
-	public static final String IRRIGATION_COST_FILE = SPATIAL_DATA_DIR + File.separator + "irrigation_cost.asc";;
+	public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries_72_96.asc";
+	public static final String IRRIGATION_COST_FILE = SPATIAL_DATA_DIR + File.separator + "irrigation_cost_72_96.asc";;
 	public static final double IRRIG_COST_SCALE_FACTOR = 3.0;
 
 	

src/ac/ed/lurg/ModelMain.java

@@ -54,7 +54,7 @@ public class ModelMain {
 
 	/* setup models, reading inputs, etc. */
 	private void setup() {
-		desiredProjection = new RasterHeaderDetails(720, 360, -180, -90, 0.5, "999");
+		desiredProjection = new RasterHeaderDetails(96, 72, -180, -90, 3.75, 2.5, "999");
 		countryToKeysMap = createCountryToKeysMap();
 		demandManager = new DemandManager(ModelFitType.LOGISTIC, ModelConfig.SSP_SCENARIO);
 		countryAgents = createCountryAgents();
@@ -285,6 +285,12 @@ public class ModelMain {
 			LandCoverReader lcReader = new LandCoverReader(initLC, landType);
 			initLC = lcReader.getRasterDataFromFile(rootDir + landType.getFileName());
 		}
+		
+		// scale all fractional areas read in by the actual area
+		for (Entry<RasterKey, LandCoverItem> entry : initLC.entrySet()) {
+			double cellArea = initLC.getAreaMha(entry.getKey());
+			entry.getValue().scaleAll(cellArea);
+		}
 
 		return initLC;
 	}
@@ -302,7 +308,7 @@ public class ModelMain {
 			}
 		} */
 
-		LPJYieldResponseMapReader yieldReader = new LPJYieldResponseMapReader(ModelConfig.YIELD_DIR);
+		LPJYieldResponseMapReader yieldReader = new LPJYieldResponseMapReader(ModelConfig.YIELD_DIR, desiredProjection);
 
 		for (FertiliserRate fr : FertiliserRate.values()) {
 			yieldReader.getRasterDataFromFile(fr); 

src/ac/ed/lurg/country/gams/GamsCountryInput.java

@@ -97,7 +97,7 @@ public class GamsCountryInput {
 //	}
 
 	public double getMeatEfficiency() {
-		return 0.5;  // this is already handled by the feed conversion efficiency for each animal product
+		return 1.0;  // this is already handled by the feed conversion efficiency for each animal product
 	}
 
 	public double getLandChangeEnergy() {

src/ac/ed/lurg/country/gams/GamsRasterTest.java

@@ -47,7 +47,7 @@ public class GamsRasterTest extends GamsLocationTest {
 	}
 	
 	public RasterSet<IrrigationCostItem> getIrrigationCost() {
-		RasterSet<IrrigationCostItem> testIrrigationCostRaster= new RasterSet<IrrigationCostItem>();
+		RasterSet<IrrigationCostItem> testIrrigationCostRaster= new RasterSet<IrrigationCostItem>(null);
 		
 		for (int i = 0; i<X_RANGE; i++) {
 			for (int j = 0; j<Y_RANGE; j++) {
@@ -63,7 +63,7 @@ public class GamsRasterTest extends GamsLocationTest {
 
 	
 	public RasterSet<AreasItem> getPreviousAreaRaster() {
-		RasterSet<AreasItem> testAreaRaster= new RasterSet<AreasItem>();
+		RasterSet<AreasItem> testAreaRaster= new RasterSet<AreasItem>(null);
 		
 		for (int i = 0; i<X_RANGE; i++) {
 			for (int j = 0; j<Y_RANGE; j++) {

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

@@ -12,7 +12,9 @@ public class IrrigiationCostReader extends AbstractRasterReader<IrrigationCostIt
 
 	@Override
 	public void setData(IrrigationCostItem item, String token) {
-		double irrigCost = Double.parseDouble(token);
-		item.setIrrigCost(irrigCost * ModelConfig.IRRIG_COST_SCALE_FACTOR);
+		if (!"nan".equals(token)) {	
+			double irrigCost = Double.parseDouble(token);
+			item.setIrrigCost(irrigCost * ModelConfig.IRRIG_COST_SCALE_FACTOR);
+		}
 	}
 }

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

@@ -2,6 +2,7 @@ package ac.ed.lurg.landuse;
 
 import java.util.HashMap;
 import java.util.Map;
+import java.util.Map.Entry;
 
 import ac.ed.lurg.types.LandCoverType;
 import ac.sac.raster.RasterItem;
@@ -20,6 +21,13 @@ public class LandCoverItem implements RasterItem {
 	public void setLandCover(LandCoverType landType, double d) {
 		landcover.put(landType, d);
 	}
+	
+	public void scaleAll(double factor) {
+		for (Entry<LandCoverType, Double> entry : landcover.entrySet()) {
+			landcover.put(entry.getKey(), entry.getValue() * factor);
+		}
+
+	}
 		
 	/*public double getTotal() {
 		double total = 0;

src/ac/ed/lurg/types/LandCoverType.java

@@ -18,7 +18,7 @@ public enum LandCoverType {
 	}
 	
 	public String getFileName() {
-		return name + ".asc";
+		return name + "_72_96.asc";
 	}
 }
 

src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java

@@ -16,9 +16,9 @@ public class LPJYieldResponseMapReader {
 	private YieldRaster dataset;
 	private String rootDir;
 
-	public LPJYieldResponseMapReader(String rootDir) {
+	public LPJYieldResponseMapReader(String rootDir, RasterHeaderDetails rasterProj) {
 		this.rootDir = rootDir;
-		dataset = new YieldRaster(new RasterHeaderDetails(720,360,-180,-90,0.5,"999"));
+		dataset = new YieldRaster(rasterProj);
 
 	}
 	

src/ac/sac/raster/AbstractRasterReader.java

@@ -10,6 +10,7 @@ import ac.ed.lurg.utils.LogWriter;
 public abstract class AbstractRasterReader<D extends RasterItem> {
 	protected static final boolean DEBUG = false;
 	private static final int HEADER_LENGTH = 6;
+	private static final int HEADER_TITLE_COL = 0;
 	private static final int HEADER_DATA_COL = 1;
 	private static final int NCOLUMN_ROW = 0;
 	private static final int NROW_ROW = 1;
@@ -37,7 +38,8 @@ public abstract class AbstractRasterReader<D extends RasterItem> {
 		int nrows=0;
 		double xllCorner=0.0;
 		double yllCorner = 0.0;
-		double cellSize = 1000;
+		double xcellSize = 0;
+		double ycellSize = 0;
 		String nodataString=null;
 		
 		for (int i=0; i<HEADER_LENGTH; i++) {
@@ -65,17 +67,26 @@ public abstract class AbstractRasterReader<D extends RasterItem> {
 				yllCorner = Double.parseDouble(tokens[HEADER_DATA_COL]);
 				break;
 			case CELL_SIZE:
-				cellSize =  Double.parseDouble(tokens[HEADER_DATA_COL]);
+				xcellSize =  Double.parseDouble(tokens[HEADER_DATA_COL]);
+				ycellSize = xcellSize;
 				break;
 			case NODATA_VALUE_ROW:
-				nodataString = tokens[HEADER_DATA_COL];
+				// deal with the case of files with x and y cell size is different.  There are still 6 rows in the header as the nodata row is missing
+				if ("dy".equals(tokens[HEADER_TITLE_COL])) {
+					ycellSize = Double.parseDouble(tokens[HEADER_DATA_COL]);
+					nodataString = "0";
+				}
+				else {
+					nodataString = tokens[HEADER_DATA_COL];
+				}
+				
 				break;
 			}
 		}
 		
 		if (DEBUG) LogWriter.println("Creating RasterDataset col:"+ncolumns + ", rows:" + nrows);
 		
-		RasterHeaderDetails headerDetails = new RasterHeaderDetails(ncolumns, nrows, xllCorner, yllCorner, cellSize, nodataString);
+		RasterHeaderDetails headerDetails = new RasterHeaderDetails(ncolumns, nrows, xllCorner, yllCorner, xcellSize, ycellSize, nodataString);
 		if (dataset == null) createDataSet(headerDetails);
 		return headerDetails;
 	}
@@ -96,7 +107,7 @@ public abstract class AbstractRasterReader<D extends RasterItem> {
 			String line;
 			
 			while ((line=in.readLine()) != null) {
-				String[] tokens = parseLine(line);
+				String[] tokens = parseLine(line.trim());
 
 				for (String token : tokens) {
 					

src/ac/sac/raster/RasterHeaderDetails.java

@@ -5,24 +5,23 @@ public class RasterHeaderDetails {
 	private int nrows;
 	private double xllCorner;
 	private double yllCorner;
-	private double cellSize;
+	private double xcellSize;
+	private double ycellSize;
 	private String nodataString;
 	
-	public RasterHeaderDetails(int ncolumns, int nrows, double xllCorner, double yllCorner, double cellSize, String nodataString) {
+	public RasterHeaderDetails(int ncolumns, int nrows, double xllCorner, double yllCorner, double xcellSize, double ycellSize, String nodataString) {
 		super();
 		this.ncolumns = ncolumns;
 		this.nrows = nrows;
 		this.xllCorner = xllCorner;
 		this.yllCorner = yllCorner;
-		this.cellSize = cellSize;
+		this.xcellSize = xcellSize;
+		this.ycellSize = ycellSize;
 		this.nodataString = nodataString;
 	}
-
-	public RasterHeaderDetails() {
-	}
 	
 	public RasterHeaderDetails copy() {
-		return new RasterHeaderDetails(ncolumns, nrows, xllCorner, yllCorner, cellSize, nodataString);
+		return new RasterHeaderDetails(ncolumns, nrows, xllCorner, yllCorner, xcellSize, ycellSize, nodataString);
 	}
 
 	public double getOriginX() {
@@ -30,7 +29,7 @@ public class RasterHeaderDetails {
 	}
 	
 	public double getOriginY() {
-		return yllCorner + cellSize * nrows;
+		return yllCorner + ycellSize * nrows;
 	}
 	
 	public void incrementNcolumns() {
@@ -39,7 +38,7 @@ public class RasterHeaderDetails {
 
 	public void incrementNrows() {
 		nrows++;
-		yllCorner = yllCorner - cellSize;
+		yllCorner = yllCorner - ycellSize;
 	}
 
 	public int getNcolumns() {
@@ -58,10 +57,15 @@ public class RasterHeaderDetails {
 		return yllCorner;
 	}
 
-	public double getCellSize() {
-		return cellSize;
+	public double getXCellSize() {
+		return xcellSize;
+	}
+	
+	public double getYCellSize() {
+		return ycellSize;
 	}
 
+
 	public String getNodataString() {
 		return nodataString;
 	}
@@ -71,15 +75,17 @@ public class RasterHeaderDetails {
 	public int hashCode() {
 		final int prime = 31;
 		int result = 1;
-		long temp;
-		temp = Double.doubleToLongBits(cellSize);
-		result = prime * result + (int) (temp ^ (temp >>> 32));
 		result = prime * result + ncolumns;
 		result = prime * result
 				+ ((nodataString == null) ? 0 : nodataString.hashCode());
 		result = prime * result + nrows;
+		long temp;
+		temp = Double.doubleToLongBits(xcellSize);
+		result = prime * result + (int) (temp ^ (temp >>> 32));
 		temp = Double.doubleToLongBits(xllCorner);
 		result = prime * result + (int) (temp ^ (temp >>> 32));
+		temp = Double.doubleToLongBits(ycellSize);
+		result = prime * result + (int) (temp ^ (temp >>> 32));
 		temp = Double.doubleToLongBits(yllCorner);
 		result = prime * result + (int) (temp ^ (temp >>> 32));
 		return result;
@@ -94,9 +100,6 @@ public class RasterHeaderDetails {
 		if (getClass() != obj.getClass())
 			return false;
 		RasterHeaderDetails other = (RasterHeaderDetails) obj;
-		if (Double.doubleToLongBits(cellSize) != Double
-				.doubleToLongBits(other.cellSize))
-			return false;
 		if (ncolumns != other.ncolumns)
 			return false;
 		if (nodataString == null) {
@@ -106,9 +109,15 @@ public class RasterHeaderDetails {
 			return false;
 		if (nrows != other.nrows)
 			return false;
+		if (Double.doubleToLongBits(xcellSize) != Double
+				.doubleToLongBits(other.xcellSize))
+			return false;
 		if (Double.doubleToLongBits(xllCorner) != Double
 				.doubleToLongBits(other.xllCorner))
 			return false;
+		if (Double.doubleToLongBits(ycellSize) != Double
+				.doubleToLongBits(other.ycellSize))
+			return false;
 		if (Double.doubleToLongBits(yllCorner) != Double
 				.doubleToLongBits(other.yllCorner))
 			return false;
@@ -121,7 +130,8 @@ public class RasterHeaderDetails {
 				", nrows=" + nrows + 
 				", xllCorner=" + xllCorner + 
 				", yllCorner=" + yllCorner + 
-				", cellSize" + cellSize + 
+				", xcellSize" + xcellSize + 
+				", ycellSize" + ycellSize + 
 				", nodataString=" + nodataString;
 	}
 

src/ac/sac/raster/RasterKey.java

@@ -37,20 +37,7 @@ public class RasterKey implements Serializable {
 		double d = Math.sqrt(Math.pow(col - x.col,2) + Math.pow(row - x.row,2));
 		return d;
 	}
-	
-	/** Assuming coordinates refer to global half degree in Mha */
-	public double getHalfDegreeArea() {
-		// https://badc.nerc.ac.uk/help/coordinates/cell-surf-area.html
-		double cellArea = Math.PI / 360.0 * (Math.sin(getLatInRadians(row)) - Math.sin(getLatInRadians(row+1))) * 6371.0 * 6371.0 / 10000.0;
-		return cellArea;
-	}
-	
-	private static double getLatInRadians(double rowNum) {
-		double latDegrees = 90.0 - rowNum / 2;
-		double latRadian = latDegrees * Math.PI / 180.0;
-		return latRadian;
-	}
-	
+		
 	@Override
 	public int hashCode() {
 		final int prime = 31;

src/ac/sac/raster/RasterOutputer.java

@@ -19,6 +19,7 @@ import java.io.IOException;
 import javax.imageio.ImageIO;
 
 import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.utils.LogWriter;
 
 public abstract class RasterOutputer<D extends RasterItem> {
 	protected RasterSet<D> results;
@@ -76,6 +77,9 @@ public abstract class RasterOutputer<D extends RasterItem> {
 	}
 
 	private void writeHeader(BufferedWriter outputFile) throws IOException {
+		if (results.getHeaderDetails().getXCellSize() != results.getHeaderDetails().getYCellSize())
+			LogWriter.printlnError("Headers x and y cell sizes differ");
+			
 		outputFile.write("ncols         " + results.getNcolumns());
     	outputFile.newLine();
 		outputFile.write("nrows         " + + results.getNrows());
@@ -84,7 +88,7 @@ public abstract class RasterOutputer<D extends RasterItem> {
     	outputFile.newLine();
 		outputFile.write("yllcorner     " + results.getHeaderDetails().getYllCorner());
     	outputFile.newLine();
-		outputFile.write("cellsize      " + results.getHeaderDetails().getCellSize());
+		outputFile.write("cellsize      " + results.getHeaderDetails().getXCellSize());
     	outputFile.newLine();
 		outputFile.write("NODATA_value  " + results.getHeaderDetails().getNodataString());
     	outputFile.newLine();

src/ac/sac/raster/RasterSet.java

@@ -13,11 +13,7 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 	
 	private static final long serialVersionUID = 4180188703734898215L;
 	private RasterHeaderDetails header;
-	
-	public RasterSet() {
-		header = new RasterHeaderDetails();
-	}
-	
+		
 	public RasterSet(RasterHeaderDetails header) {
 		this.header = header;
 	}
@@ -29,18 +25,19 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 		int c, r;
 		
 		if (source!=null && !source.equals(header)) {
-			int xOffset = (int) Math.round((source.getXllCorner() - header.getXllCorner())/header.getCellSize());
-			int yOffset = (int) Math.round((-source.getOriginY() + header.getOriginY())/header.getCellSize());  // down from top
+			int xOffset = (int) Math.round((source.getXllCorner() - header.getXllCorner())/header.getXCellSize());
+			int yOffset = (int) Math.round((-source.getOriginY() + header.getOriginY())/header.getYCellSize());  // down from top
 
-			if (source.getCellSize() % header.getCellSize() != 0)
+			if (source.getXCellSize() % header.getXCellSize() != 0 && source.getYCellSize() % header.getYCellSize() != 0)
 				throw new RuntimeException("Cell sizes must be the same or souce an integer multiple of destination");
 			
-			double cellRatio = source.getCellSize() / header.getCellSize();
-			
-			for (int x = 0; x<cellRatio; x++) {
-				for (int y = 0; y<cellRatio; y++) {
-					c = (int)(col*cellRatio) + x + xOffset;
-					r = (int)(row*cellRatio) + y + yOffset;
+			double xcellRatio = source.getXCellSize() / header.getXCellSize();
+			double ycellRatio = source.getYCellSize() / header.getYCellSize();
+
+			for (int x = 0; x<xcellRatio; x++) {
+				for (int y = 0; y<ycellRatio; y++) {
+					c = (int)(col*xcellRatio) + x + xOffset;
+					r = (int)(row*ycellRatio) + y + yOffset;
 
 					D d = get(c, r);
 					if (d != null)
@@ -59,10 +56,10 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 
 	/** Method to get raster item for coordinates */
 	public D getFromCoordinates(double source_x, double source_y) {
-		int col = (int) Math.round((source_x - header.getXllCorner())/header.getCellSize());
+		int col = (int) Math.round((source_x - header.getXllCorner())/header.getXCellSize());
 		
 		// header.getNrows() - 1, lower left basis.  Minus 1 as nrows is number, but indexed from zero
-		int row = header.getNrows() - 1 - (int) Math.round((source_y - header.getYllCorner())/header.getCellSize());
+		int row = header.getNrows() - 1 - (int) Math.round((source_y - header.getYllCorner())/header.getYCellSize());
 	
 		if (row < 0 || col < 0) {
 			LogWriter.println("Got negative values");
@@ -71,12 +68,17 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 	}
 	
 	public double getXCoordin(RasterKey key) {
-		double x = header.getCellSize() * key.getCol() + header.getXllCorner();
+		double x = header.getXCellSize() * key.getCol() + header.getXllCorner();
 		return x;
 	}
 	
 	public double getYCoordin(RasterKey key) {
-		double y = header.getCellSize() * (header.getNrows() - 1 - key.getRow()) + header.getYllCorner();
+		double y = getYCoordin(key.getRow());
+		return y;
+	}
+	
+	private double getYCoordin(int row) {
+		double y = header.getYCellSize() * (header.getNrows() - 1 - row) + header.getYllCorner();
 		return y;
 	}
 
@@ -98,6 +100,22 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 		return data;
 	}
 	
+	// Assuming coordinates refer to global half degree in Mha 
+	public double getAreaMha(RasterKey key) {
+		// https://badc.nerc.ac.uk/help/coordinates/cell-surf-area.html
+		// http://mathforum.org/library/drmath/view/63767.html
+		int aRow = key.getRow();
+		double a = Math.sin(getRadiansFromDeg(getYCoordin(aRow)));
+		double b = Math.sin(getRadiansFromDeg(getYCoordin(aRow+1)));
+		double cellArea = header.getXCellSize() * Math.PI / 180.0 * (a-b) * 6371.0 * 6371.0 / 10000.0;
+		return cellArea;
+	}
+	
+	private static double getRadiansFromDeg(double deg) {
+		double latRadian = deg * Math.PI / 180.0;
+		return latRadian;
+	}
+
 	
 	/** Some classes do not create, those that do will need to override */
 	protected D createRasterData() {
@@ -119,7 +137,7 @@ public class RasterSet<D extends RasterItem> extends HashMap<RasterKey, D> {
 	/** Check passed in header is consistent with one from this set.
 	 *   Actually only checks grid size as can shift rasters*/
 	public boolean isConstistentWithHeader(RasterHeaderDetails h) {
-		return h.getCellSize() == header.getCellSize();
+		return h.getXCellSize() == header.getXCellSize() && h.getYCellSize() == header.getYCellSize();
 	}
 	
 	public RasterSet<D> popSubsetForKeys(RasterSet<D> subset, Collection<RasterKey> keys) {

src/ac/sac/raster/RasterSetTest.java

@@ -8,7 +8,7 @@ public class RasterSetTest {
 
 	@Test
 	public void test() {
-		YieldRaster dataset = new YieldRaster(new RasterHeaderDetails(720,360,-180,-90,0.5,""));
+		YieldRaster dataset = new YieldRaster(new RasterHeaderDetails(720,360,-180,-90,0.5,0.5,""));
 		double y = dataset.getYCoordin(new RasterKey(0,179));
 		System.out.println("y=" +y);