diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java index 2f5e043ea678b05e7dd8b15c6f1b8ee26194119c..0113d0721fa5ec8e17ceb5db754d837e2992eee4 100644 --- a/src/ac/ed/lurg/ModelConfig.java +++ b/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; diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java index 90fbf0a41302507480d5308d558f1e2291486f9b..e8f6cbfab99ac6c5a0cc573eb2dd388280483bfe 100644 --- a/src/ac/ed/lurg/ModelMain.java +++ b/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); diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java index 65ad31b9b7d71cf48357686a9d8022942c801d5a..309335a22f2d7f4b87ef0bc04cc998aedc8617cc 100644 --- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java +++ b/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() { diff --git a/src/ac/ed/lurg/country/gams/GamsRasterTest.java b/src/ac/ed/lurg/country/gams/GamsRasterTest.java index 8cc2a1512d16e14cf0c0c90693014f1fbdf44cb7..98a9b4e6b9ada411334962b2b7a906007b3368cc 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterTest.java +++ b/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++) { diff --git a/src/ac/ed/lurg/landuse/IrrigiationCostReader.java b/src/ac/ed/lurg/landuse/IrrigiationCostReader.java index 742bb88886f21ec1f73c8df6a50ac27dc2528a57..d3a0d5888e24291dd55dcd6a16703d653e8ed763 100644 --- a/src/ac/ed/lurg/landuse/IrrigiationCostReader.java +++ b/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); + } } } diff --git a/src/ac/ed/lurg/landuse/LandCoverItem.java b/src/ac/ed/lurg/landuse/LandCoverItem.java index 13f6bc77857db44644dbba4152baea079ac135f6..0ff7f0c40cdb7692c9a37b6036f084121c66a901 100644 --- a/src/ac/ed/lurg/landuse/LandCoverItem.java +++ b/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; diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java index 6908a9d759e9ce9127913b8b2615128842b0ee7f..2acb6998fcf23c24a3b2ed4e48c6ff5f422f97dc 100644 --- a/src/ac/ed/lurg/types/LandCoverType.java +++ b/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"; } } diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java index 2cde5a8b250f7d859d57c9e14906ffed87407955..bcdea941b68e3697f28bdf2faf6d7bd76ff1623f 100644 --- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java +++ b/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); } diff --git a/src/ac/sac/raster/AbstractRasterReader.java b/src/ac/sac/raster/AbstractRasterReader.java index 4812275f50f08d638494022045b4b9332dbcd8d0..2042959e88dba5b8327c5be0da22223d5cb1eee4 100755 --- a/src/ac/sac/raster/AbstractRasterReader.java +++ b/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) { diff --git a/src/ac/sac/raster/RasterHeaderDetails.java b/src/ac/sac/raster/RasterHeaderDetails.java index 1878cab389b501d87a6d547156c54e83720c0ec9..45f77a566f600d2be99a7597198058812da05ad4 100755 --- a/src/ac/sac/raster/RasterHeaderDetails.java +++ b/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; } diff --git a/src/ac/sac/raster/RasterKey.java b/src/ac/sac/raster/RasterKey.java index 678e1a128e0193a6fb3e3872a3cd4b4ebfda79c0..cde885f873724bcb21206ae837718a7cd5891c00 100755 --- a/src/ac/sac/raster/RasterKey.java +++ b/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; diff --git a/src/ac/sac/raster/RasterOutputer.java b/src/ac/sac/raster/RasterOutputer.java index ff4c3b7f06de9b6577ee6218599b2af65ceda2a2..8932b11bea10bcd2d115a8910952bee0ccfe0529 100644 --- a/src/ac/sac/raster/RasterOutputer.java +++ b/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(); diff --git a/src/ac/sac/raster/RasterSet.java b/src/ac/sac/raster/RasterSet.java index 768c30e7e8427f7cbf50e91972d50bf221d2de14..c67d7ebb9f9b04e1227253f22a4ae52ac7e5c4a6 100755 --- a/src/ac/sac/raster/RasterSet.java +++ b/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) { diff --git a/src/ac/sac/raster/RasterSetTest.java b/src/ac/sac/raster/RasterSetTest.java index 857588fe3bc3707c20a921b39a78be705ee469e1..1871791e9a7538aaad48f209df34d01394a63430 100644 --- a/src/ac/sac/raster/RasterSetTest.java +++ b/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);