/*
 * AbstractGrid2DSquareCellDouble.java
 * Created on 14 December 2001, 17:00
 */
package uk.ac.leeds.ccg.cluster.grids;
import uk.ac.leeds.ccg.raster.*;
import uk.ac.leeds.ccg.geotools.*;
import java.io.*;
import java.lang.*;
import java.util.*;
/**
 * Last updated on 3 Jan 2002
 * @author  andyt
 * @version 1.0
 *
 * Description:
 * An abstract class for representing a type of spatial data. The spatial data
 * are in a Euclidean two-dimensional geometry given by an X axis and a Y axis.
 * The spatial extent of the data is given by a rectangle with its sides aligned
 * with the axes. The spatial data are quanta represented by doubles and
 * referenced to square shaped partitions of the extent, known as cells, which
 * are aligned in rows and columns.
 */

public abstract class AbstractGrid2DSquareCellDouble {
    
    /**
     * debugging level
     * DEBUG=0 produces no debugging information
     * DEBUG=1 provides information about which method is running
     * DEBUG>2 provides information about the running method
     */
    protected final int DEBUG=2;
    
    /**
     * Reference to the instance of the GridStatistics object
     */
    protected AbstractGridStatistics gridStatistics;

    /**
     * Grid descriptor fields
     */
    
    /**
     * For storing the number of columns in the grid.
     */
    protected int ncols;
    
    /**
     * For storing the number of rows in the grid.
     */
    protected int nrows;
    
    /**
     * If the Y axis is verticle and the X axis is horizontal then xllcorner is 
     * for storing the left x-axis cordinate of the grid.
     */
    protected double xllcorner=0.0d;

    /**
     * If the Y axis is verticle and the X axis is horizontal then yllcorner is 
     * for storing the bottom y-axis cordinate of the grid.
     */
    protected double yllcorner=0.0d;

    /**
     * For storing the cellsize - the width (or height) of each cell.
     */
    protected double cellsize=1.0d;

    /**
     * For storing the noDataValue of the grid
     */
    protected double noDataValue=Double.NaN;
    
    
    /**
     * Basic methods
     */

    /**
     * Returns a string description of the grid.
     */
    public abstract String toString();
    
    /**
     * Returns the dataType of the grid.
     * dataType=0 means the data is held in memory
     * dataType=1 means the data is held in a file
     * dataType=2 means the data is held in a database 
     */
    public abstract int getDataType();
    
    /**
     * Returns the reference to the instance of the GridStatistics object
     */
    public final AbstractGridStatistics getGridStatistics() {return this.gridStatistics;}

    /**
     * Initialises the reference to the instance of the GridStatistics object
     */
    protected final void initGridStatistics(AbstractGridStatistics gridStatistics) {this.gridStatistics=gridStatistics;}

    /**
     * Returns the number of columns in the grid.
     */
    public final int getNcols() {return this.ncols;}

    /**
     * Initialises ncols.
     */
    protected final void initNcols(int ncols) {this.ncols=ncols;}

    /**
     * Returns the number of rows in the grid.
     */
    public final int getNrows() {return this.nrows;}
    
    /**
     * Initialises nrows.
     */
    protected final void initNrows(int nrows) {this.nrows=nrows;}

    /**
     * Returns the x-coordinate of the lower left corner of the grid.
     */
    public final double getXllcorner() {return this.xllcorner;}
    
    /**
     * Sets the x-coordinate of the lower left corner of the grid.
     * Returns what was the x-coordinate of the left corner of the grid.
     */
    public final double setXllcorner(double d1) {
        double wasXllcorner=xllcorner;
        this.xllcorner=d1;
        return wasXllcorner;
    }
    
    /**
     * Returns the y-coordinate of the lower left corner of the grid.
     */
    public final double getYllcorner() {return this.yllcorner;}
    
    /**
     * Sets the y-coordinate of the lower left corner of the grid.
     * Returns what was the y-coordinate of the lower left corner of the grid.
     */
    public final double setYllcorner(double d1) {
        double wasYllcorner=yllcorner;
        this.yllcorner=d1;
        return wasYllcorner;
    }
    
    /**
     * Returns the cell width (cell height) of the grid cells.
     */
    public final double getCellsize() {return this.cellsize;}
    
    /**
     * Sets the cell width (cell height) of the grid cells.
     * Returns what was the cellsize of the grid
     * This effectively scales the grid keeping the origin fixed.
     */
    public final double setCellsize(double d1) {
        double wasCellsize=cellsize;
        this.cellsize=d1;
        return wasCellsize;
    }
    
    /**
     * Returns the noDataValue for the grid cells.
     */
    public final double getNoDataValue() {return this.noDataValue;}

    /**
     * Initialises noDataValue
     */
    protected final void initNoDataValue(double noDataValue) {this.noDataValue=noDataValue;}

    /**
     * Sets the noDataValue for the grid cells to d1.
     */
    public abstract void setNoDataValue(double d1);
    
    /**
     * Sets numbers in the range [d1,d2] to noDataValue. All existing noDataValues
     * remain as noDataValues.
     */
    public void setToNoDataValue(double d1,double d2) {
        double d3;
        int ncols = getNcols();
        int nrows = getNrows();
        double noDataValue = getNoDataValue();
        for (int i = 0; i < nrows; i ++) {
            for (int j = 0; j < ncols; j ++) {
                d3 = getCell(i,j);
                if (d1 >= d3 && d3 <= d2) {
                    setCell(i,j,noDataValue);
                }
            }
        }
    }
    
    /**
     * Returns the height of the grid
     */
    public final double getHeight() {return getNrows() * getCellsize();}
    
    /**
     * Returns the width of the grid
     */
    public final double getWidth() {return getNcols() * getCellsize();}
    
    /**
     * Returns true iff x-coordinate x, y-coordinate y is in the Grid.
     */
    public final boolean inGrid(double x,double y) {
        double xllcorner = getXllcorner();
        double yllcorner = getYllcorner();
        double cellsize = getCellsize();
        return (x >= xllcorner && x < (xllcorner + (getNcols() * cellsize)) && y >= yllcorner && y < (yllcorner + (getNrows() * cellsize)));
    }
    
    /**
     * Returns true iff row i,column j is in the Grid.
     */
    public final boolean inGrid(int i,int j) {return (i >= 0 && i < getNrows() && j >= 0 && j < getNcols());}
    
    /**
     * Returns true iff cellID is a valid cellID.
     */
    public final boolean inGrid(int cellID) {return (cellID >= 0 && cellID < getNrows()*getNcols());}
    
    /**
     * Returns the cellID if row and col are respectively row and column indexes
     * that reference a cell in the grid, else returns Integer.MIN_VALUE.
     * Note:
     * Cells are numbered from 0 to getNrows()*getNcols()-1 in row major order and
     * only have an ID if they are in the grid.
     */
    public final int getCellID(int row,int col) {
        if (inGrid(row,col)) {return (row * getNcols()) + col;}
        return Integer.MIN_VALUE;
    }
    public final int getCellID(double x,double y) {return getCellID(getRowIndex(y),getColIndex(x));}
    
    /**
     * Returns the row index for the input y-coordinate. Rows in the grid are
     * indexed from 0 to nrows-1.
     */
    public final int getRowIndex(double y) {
        double cellsize = getCellsize();
        return ((int) Math.floor((((double) getNrows() * cellsize) + getYllcorner() - y) / cellsize));
    }
    /**
     * Returns the row index for the input cellID.
     */
    public final int getRowIndex(int cellID) {
        if (inGrid(cellID)) {return (int) cellID / getNcols();}
        return Integer.MIN_VALUE; // This should not happen given a valid cellID.
    }
    
    /**
     * Returns the column index for the input x-coordinate. Columns in the grid are
     * indexed from 0 to ncols-1.
     */
    public final int getColIndex(double x) {return (int) Math.floor((x - getXllcorner()) / getCellsize());}
    /**
     * Returns the column index for the input cellID.
     */
    public final int getColIndex(int cellID) {
        if (inGrid(cellID)) {return cellID - (getRowIndex(cellID) * getNcols());}
        return Integer.MIN_VALUE; // This should not happen given a valid cellID.
    }
    
    /**
     * Returns the x-coordinate of the cell centroid.
     */
    public final double getCellX(int cellID) {
        double cellsize = getCellsize();
        return getXllcorner() + (cellsize * getColIndex(cellID)) + (cellsize / 2.0d);
    }
    public final double getCellX(int i,int j) {
        //return getCellX(getCellID(i,j));}
        double cellsize=getCellsize();
        return getXllcorner() + (cellsize * i) + (cellsize / 2.0d);
    }

    /**
     * Returns the y-coordinate of the cell centroid.
     */
    public final double getCellY(int cellID) {
        double cellsize = getCellsize();
        return getYllcorner() + (cellsize * (getNrows() - getRowIndex(cellID))) - (cellsize / 2.0d);
    }
    public final double getCellY(int i,int j) {
        //return getCellY(getCellID(i,j));
        double cellsize = getCellsize();
        return getYllcorner() + (cellsize * (getNrows() - j)) - (cellsize / 2.0d);
    }
    
    /**
     * Returns the value at row i, column j.
     */
    public abstract double getCell(int i,int j);
    /** Returns the value of the cell containing x-coordinate x, y-coordinate y. */
    public double getCell(double x,double y) {return getCell(getRowIndex(y),getColIndex(x));}
    /** Returns the value of the cell with cellID. */
    public double getCell(int cellID) {return getCell(getRowIndex(cellID),getColIndex(cellID));}
    
    /**
     * Initilises the value at row i column j to d1.
     */
    protected abstract void initCell(int i,int j,double d1);
    // Initialises the value at x-coordinate x, y-coordinate y to d1.
    protected final void initCell(double x,double y,double d1) {initCell(getRowIndex(y),getColIndex(x),d1);}
    // Initialises the value at cellID to d1.
    protected final void initCell(int cellID,double d1) {initCell(getRowIndex(cellID),getColIndex(cellID),d1);}
    
    /**
     * Sets the value of the cell at row i column j to d1 and returns the original value of the cell.
     */
    public abstract double setCell(int i,int j,double d1);
    /** Sets the value at x-coordinate x, y-coordinate y to d1. */
    public final double setCell(double x,double y,double d1) {return setCell(getRowIndex(y),getColIndex(x),d1);}
    /** Sets the value of the cell cellID to d1. */
    public final double setCell(int cellID,double d1) {return setCell(getRowIndex(cellID),getColIndex(cellID),d1);}
    
    /**
     * Adds d1 to cell at row i column j and returns the cell value that was replaced if cell is in Grid otherwise returns noDataValue.
     */
    public final double addToCell(int i,int j,double d1) {
        double noDataValue=getNoDataValue();
        if (inGrid(i,j)) {
            double d2 = getCell(i,j);
            if (d2 != noDataValue) {
                if (d1 != noDataValue) {
                    setCell(i,j,(d1 + d2));
                } else {
                    setCell(i,j,d2);
                }
            } else {
                if (d1 != noDataValue) {
                    setCell(i,j,d1);
                }
            }
            return d2;
        }
        return noDataValue;
    }
    // Adds d1 to cell containing the x-coordinate x, y-coordinate y.
    public final double addToCell(double x,double y,double d1) {return addToCell(getRowIndex(y),getColIndex(x),d1);}
    // Adds d1 to cell cellID.
    public final double addToCell(int cellID,double d1) {return addToCell(getRowIndex(cellID),getColIndex(cellID),d1);}
    
    /**
     * Returns the bounding GeoRectangle of the grid.
     */
    public final GeoRectangle getBoundingGeoRectangle() {
        double cellsize = getCellsize();
        return new GeoRectangle(getXllcorner(),getYllcorner(),cellsize * getNcols(),cellsize * getNrows());
    }
    
    /**
     * Returns the distance between two cell centroids in map units
     */
    public final double distance(int cellID1, int cellID2) {
        return Math.sqrt(Math.pow(getCellX(cellID1) - getCellX(cellID2),2.0d) + Math.pow(getCellY(cellID1) - getCellY(cellID2),2.0d));
    }
    
    /**
     * Returns the distance between two cell centroids in cell units
     */
    public final double cellDistance(int cellID1, int cellID2) {
        return (Math.sqrt(Math.pow(getCellX(cellID1) - getCellX(cellID2),2.0d) + Math.pow(getCellY(cellID1) - getCellY(cellID2),2.0d))) / cellsize;
    }
    
    /**
     * Grid operations
     */
    
    /**
     * Mask methods
     * Masking a grid involves setting values in it to its noDataValue if values in the mask have the mask noDataValue.
     * ToDos;
     * 1. Ability to handle inputGrids with different spatial frames.
     * @param mask - the AbstractGrid2DSquareCellDouble with which to mask
     */
    public final void mask(AbstractGrid2DSquareCellDouble mask) {
        if (DEBUG>1) {System.out.println("mask(AbstractGrid2DSquareCellDouble "+mask.toString()+") {");}
        double d1;
        double maskNoDataValue = mask.getNoDataValue();
        double noDataValue = getNoDataValue();
        for (int i = 0; i < mask.getNrows(); i ++) {
            for (int j = 0; j < mask.getNcols(); j ++) {
                d1 = mask.getCell(i,j);
                if (d1 == maskNoDataValue) {
                    setCell(i,j,noDataValue);
                }
            }
        }
        if (DEBUG > 1) {System.out.println("}");}
    }
    
    /**
     * Subtracts input grid values and returns the answer as a new grid
     * ToDos:
     * 1. Extend for other measures of difference (or error):
     *      - absolute difference, SSE, RMSE etc...
     * 2. Optional handling of noDataValues.
     * 3. Ability to handle inputGrids with different spatial frames.
     * @param inputGrid - the Abstract2DSquareCellDouble from which the difference will be calculated 
     */
     public final AbstractGrid2DSquareCellDouble diff(AbstractGrid2DSquareCellDouble inputGrid) {
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return diff(inputGrid,gridFactory);
        } catch (OutOfMemoryError oome1) {
            gridFactory = new Grid2DSquareCellDoubleFileFactory();
            return diff(inputGrid,gridFactory);
        }
     }
     /*
      * @param gridFactory - the Abstract2DSquareCellDoubleFactory used to create the output grid
      */
     public final AbstractGrid2DSquareCellDouble diff(AbstractGrid2DSquareCellDouble inputGrid,AbstractGrid2DSquareCellDoubleFactory gridFactory) {
        if (DEBUG > 2) {System.out.println("diff(inputGrid("+inputGrid.toString()+"),AbstractGrid2DSquareCellDoubleFactory("+gridFactory.toString()+")) {");}
        int inputGridNrows = inputGrid.getNrows();
        int inputGridNcols = inputGrid.getNcols();
        double inputGridXllcorner = inputGrid.getXllcorner();
        double inputGridYllcorner = inputGrid.getYllcorner();
        double inputGridCellsize = inputGrid.getCellsize();
        double inputGridNoDataValue = inputGrid.getNoDataValue();
        int nrows = getNrows();
        int ncols = getNcols();
        double xllcorner = getXllcorner();
        double yllcorner = getYllcorner();
        double cellsize = getCellsize();
        double noDataValue = getNoDataValue();
        double d1 = 0.0;
        double d2 = 0.0;
        // Report parameter differences
        if (nrows != inputGridNrows) {System.out.println("nrows("+nrows+") != inputGridNrows("+inputGridNrows+")");}
        if (ncols != inputGridNcols) {System.out.println("ncols("+ncols+") != inputGridNcols("+inputGridNcols+")");}
        if (xllcorner != inputGridXllcorner) {System.out.println("xllcorner("+xllcorner+") != inputGridXllcorner("+inputGridXllcorner+")");}
        if (yllcorner != inputGridYllcorner) {System.out.println("yllcorner("+yllcorner+") != inputGridYllcorner("+inputGridYllcorner+")");}
        if (cellsize != inputGridCellsize) {System.out.println("cellsize("+cellsize+") != inputGridCellsize("+inputGridCellsize+")");}
        if (noDataValue != inputGridNoDataValue) {System.out.println("noDataValue("+noDataValue+") != inputGridNoDataValue("+inputGridNoDataValue+")");}
        // Assess if inputGrid has same spatial frame
        // Calculate the intersection
        int inputGridStartRow = (int) ((inputGridYllcorner - yllcorner) / cellsize);
        int inputGridEndRow = (inputGridStartRow + inputGridNrows - 1);
        int inputGridStartCol = (int) ((inputGridXllcorner - xllcorner) / cellsize);
        int inputGridEndCol = (inputGridStartCol + inputGridNcols - 1);
        int startRow = Math.max(inputGridStartRow,0);
        int endRow = Math.min(inputGridEndRow,(nrows - 1));
        int startCol = Math.max(inputGridStartCol,0);
        int endCol = Math.min(inputGridEndCol,(ncols - 1));
        // Initialise output grid
        //AbstractGrid2DSquareCellDouble outputGrid = gridFactory.createGrid2DSquareCellDouble(inputGridNrows,inputGridNcols,inputGridXllcorner,inputGridYllcorner,inputGridCellsize,inputGridNoDataValue);
        AbstractGrid2DSquareCellDouble outputGrid = gridFactory.createGrid2DSquareCellDouble((endRow - startRow +1),(endCol - startCol +1),Math.max(xllcorner,inputGridXllcorner),Math.max(yllcorner,inputGridYllcorner),cellsize,noDataValue);
        // Perform difference operation
        //for (int i = 0; i < nrows; i ++) {
        //    for (int j = 0; j < ncols; j ++) {
        for (int i = startRow; i <= endRow; i ++) {
            for (int j = startCol; j <= endCol; j ++) {
                d1 = getCell(i,j);
                d2 = inputGrid.getCell(i+inputGridStartRow,j+inputGridStartCol);
                if (d1 != noDataValue && d2 != inputGridNoDataValue) {
                    outputGrid.setCell(i,j,(d1 - d2));
                } else {
                    outputGrid.setCell(i,j,noDataValue);
                }
            }
        }
        if (DEBUG > 1) {System.out.println("}");}
        return outputGrid;
    }
    
    /**
     * Aggregate methods
     * ToDos:
     * 1. Add normalisationLevel parameter
     * 2. Add visualisationLevel parameter
     * 3. Make it so xllcorner and yllcorner can represent any corner of a spatial frame
     * @param xllcorner - the x-coordinate of the aggregate grid lower left corner
     * @param yllcorner - the y-coordinate of the aggregate grid lower left corner
     * @param cellfactor - the number of times wider/heigher the aggregated grid cells are
     */
     public final AbstractGrid2DSquareCellDouble aggregate(int cellFactor) {
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return aggregate(cellFactor,gridFactory);
        } catch (java.lang.OutOfMemoryError oome1) {
            gridFactory = new Grid2DSquareCellDoubleFileFactory();
            return aggregate(cellFactor,gridFactory);
        }
     }
    /**
     * @param factory - the Abstract2DSquareCellDoubleFactory used to create temporary and output grids
     */
    public final AbstractGrid2DSquareCellDouble aggregate(int cellFactor,AbstractGrid2DSquareCellDoubleFactory gridFactory) {
        return aggregate(cellFactor,gridFactory,getXllcorner(),getYllcorner(),getNoDataValue());
    }
    public final AbstractGrid2DSquareCellDouble aggregate(int cellFactor,double outputGridXllcorner,double outputGridYllcorner) {
        double noDataValue = getNoDataValue();
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return aggregate(cellFactor,gridFactory,outputGridXllcorner,outputGridYllcorner,noDataValue);
        } catch (java.lang.OutOfMemoryError oome1) {
            gridFactory = new Grid2DSquareCellDoubleFileFactory();
            return aggregate(cellFactor,gridFactory,outputGridXllcorner,outputGridYllcorner,noDataValue);
        }
    }
    public final AbstractGrid2DSquareCellDouble aggregate(int cellFactor,AbstractGrid2DSquareCellDoubleFactory gridFactory,double outputGridXllcorner,double outputGridYllcorner,double noDataValue) {
        if (DEBUG > 2) {System.out.println("aggregate(cellFactor("+cellFactor+"),AbstractGrid2DSquareCellDoubleFactory("+gridFactory.toString()+"),outputGridXllcorner("+outputGridXllcorner+"),outputGridYllcorner("+outputGridYllcorner+"),noDataValue("+noDataValue+")) {");}
        if (cellFactor == 0) {
            if (DEBUG>0) {System.out.println("That's odd cellFactor is zero returning empty one cell grid!");}
            return gridFactory.createGrid2DSquareCellDouble(1,1,getXllcorner(),getYllcorner(),Math.max(getNrows(),getNcols())*getCellsize(),getNoDataValue());
        }
        AbstractGrid2DSquareCellDouble outputGrid;
        int nrows = getNrows();
        int ncols = getNcols();
        double xllcorner = getXllcorner();
        double yllcorner = getYllcorner();
        double cellsize = getCellsize();
        double outputGridCellsize = cellsize * (double) cellFactor;
        double width = cellsize * ncols;
        double height = cellsize * nrows;
        double d1 = 0.0d;
        double d2 = 0.0d;
        double d3 = 0.0d;
        // Calculate outputGridNrows and outputGridHeight
        int outputGridNrows = 1;
        double outputGridHeight = outputGridCellsize;
        while ((yllcorner + outputGridHeight) < (yllcorner + height)) {
            outputGridNrows ++;
            outputGridHeight += outputGridCellsize;
        }
        // Calculate outputGridNcols and outputGridWidth
        int outputGridNcols = 1;
        double outputGridWidth = outputGridCellsize;
        while ((xllcorner + outputGridWidth) < (xllcorner + width)) {
            outputGridNcols ++;
            outputGridWidth += outputGridCellsize;
        }
        // Initialise outputGrid
        outputGrid = gridFactory.createGrid2DSquareCellDouble(outputGridNrows,outputGridNcols,outputGridXllcorner,outputGridYllcorner,outputGridCellsize,noDataValue);
        // Aggregate
        // NB. The values are normalised by dividing the aggregate Grid sum by the proportion of cells with grid values.
        for (int i = 0; i < outputGridNrows; i ++) {
            for (int j = 0; j < outputGridNcols; j ++) {
                d2 = 0.0d;
                d3 = 0.0d;
                for (int p = (0 - (int) ((yllcorner-yllcorner) / cellsize)); p < (cellFactor - (int) ((yllcorner-yllcorner) / cellsize)); p ++) {
                    for (int q = (0 - (int) ((xllcorner - xllcorner) / cellsize)); q < (cellFactor - (int) ((xllcorner - xllcorner) / cellsize)); q ++) {
                        d1 = getCell((i * cellFactor) + p,(j * cellFactor) + q);
                        if (d1 != noDataValue) {
                            d2 += d1;
                            d3 ++;
                        }
                    }
                }
                if (d3 != 0.0d) {
                    outputGrid.setCell(i,j,(d2 / (Math.pow(cellFactor,2.0d) / d3)));
                }
            }
        }
        System.out.println(outputGrid.toString());
        if (DEBUG > 1) {System.out.println("}");}
        return outputGrid;
    }
    
    /**
     * Smooth methods
     * ToDos:
     * 1. Add normalisationLevel parameter
     * 2. Add visualisationLevel parameter
     * 3. Enhance distance weighting functionality so that various kernels can be used
     * @param cellDistance - the cell distance within which values will be used
     */
    public final AbstractGrid2DSquareCellDouble smooth(int cellDistance) {
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return smooth(cellDistance,1.0,gridFactory);
        } catch (java.lang.OutOfMemoryError oome1) {
            gridFactory = new Grid2DSquareCellDoubleFileFactory();
            return smooth(cellDistance,1.0,gridFactory);
        }
    }
    /**
     * @param distanceWeightFactor - indicates the distance weighting factor applied during smoothing
     */
    public final AbstractGrid2DSquareCellDouble smooth(int cellDistance,double distanceWeightFactor) {
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return smooth(cellDistance,distanceWeightFactor,gridFactory);
        } catch (java.lang.OutOfMemoryError oome1) {
            gridFactory = new Grid2DSquareCellDoubleFileFactory();
            return smooth(cellDistance,distanceWeightFactor,gridFactory);
        }
    }
    /**
     * @param gridFactory - the AbstractFactory2DSquareCellDouble used to create outputGrids
     */
    public final AbstractGrid2DSquareCellDouble smooth(int cellDistance,double distanceWeightFactor,AbstractGrid2DSquareCellDoubleFactory gridFactory) {
        return smooth((double)cellDistance*getCellsize(),distanceWeightFactor,gridFactory);
    }
    /**
     * @param distance - the map distance within which values will be used
     */
    public final AbstractGrid2DSquareCellDouble smooth(double distance,double distanceWeightFactor,AbstractGrid2DSquareCellDoubleFactory gridFactory) {
        if (DEBUG > 1) {System.out.println("smooth(distance("+distance+"),distanceWeightFactor("+distanceWeightFactor+"),AbstractGrid2DSquareCellDoubleFactory("+gridFactory.toString()+")) {");}
        AbstractGrid2DSquareCellDouble outputGrid;
        int ncols = getNcols();
        int nrows = getNrows();
        double xllcorner = getXllcorner();
        double yllcorner = getYllcorner();
        double cellsize = getCellsize();
        double noDataValue = getNoDataValue();
        int cellDistance = (int) Math.ceil(distance / cellsize);
        int cellID1 = 0;
        int cellID2 = 0;
        double thisCellDistance = 0.0d;
        double d1 = 0.0d;
        double d2 = 0.0d;
        double d3 = 0.0d;
        // Initialise outputGrid
        outputGrid = gridFactory.createGrid2DSquareCellDouble(nrows,ncols,xllcorner,yllcorner,cellsize,noDataValue);
        // Smooth
        // NB. The number of noDataValue cells are taken into account
        for (int i = 0; i < nrows; i ++) {
            for (int j = 0; j < ncols; j ++) {
                d1 = 0.0d;
                d2 = 0.0d;
                cellID1 = getCellID(i,j);
                for (int p = -cellDistance; p <= cellDistance; p ++) {
                    for (int q = -cellDistance; q <= cellDistance; q ++) {
                        d3 = getCell((i + p),(j + q));
                        if (d3 != noDataValue) {
                            cellID2 = getCellID((i + p),(j + q));
                            thisCellDistance = cellDistance(cellID1,cellID2);
                            if (thisCellDistance <= cellDistance) {
                                if (thisCellDistance != 0.0d) {
                                    d1 += (d3 / Math.pow(thisCellDistance,distanceWeightFactor));
                                } else {
                                    d1 += d3;
                                }
                                d2 += 1.0d;
                            }
                        }
                    }
                }
                if (d2!=0.0d) {
                    outputGrid.setCell(i,j,(d1 / d2));
                }
            }
        }
        if (DEBUG>1) {System.out.println("}");}
        return outputGrid;
    }
    
    /**
     * geometricDensity methods
     * The algorithm for generating a geometric density surface is described in:
     * Turner A (2000) Density Data Generation for Spatial Data Mining Applications.
     * http://www.geog.leeds.ac.uk/people/a.turner/papers/geocomp00/gc_017.htm
     *
     * Notes:
     * 1. If normalisationLevel=0 then there is no normalisation.
     * 2. If normalisationLevel=1 then the input and the output are forced into the
     *   range [0.0d,1.0d].
     * 3. If normalization level=2 then the input the output and the density output
     *   at each geometric scale (iteration) are forced into the range [0.0d,1.0d].
     * 4. The cellDistance parameter allows control for the number of iterations.
     *   By default it is the maximum of the rows and columns. It is effectively
     *   a limit beyond which cells of the grid are unrelated.
     * 5. The chunkSize parameter controls the size of each chunk processed. By
     *   default it is the cellDistance. If it is not >=cellDistance then it is set
     *   to cellDistance.
     *
     * ToDos:
     * 1. Enable mid-scale results to be kept if possible for further investigation.
     * 2. Enable visualisations of the algorithm via a visualisationLevel.
     * 3. Enable the algortihm to restart on partially completed runs.
     */
    public final AbstractGrid2DSquareCellDouble geometricDensity() {
        return geometricDensity(0,(Math.max(getNrows(),getNcols()) / 2));
    }
    public final AbstractGrid2DSquareCellDouble geometricDensity(int normalisationLevel,int cellDistance) {
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return geometricDensity(normalisationLevel,cellDistance,gridFactory);
        } catch (java.lang.OutOfMemoryError oome1) {
            if (cellDistance < (Math.max(getNrows(),getNcols()) / 2)) {
                int chunkSize = cellDistance;
                return geometricDensity(normalisationLevel,cellDistance,chunkSize);
            } else {
                gridFactory = new Grid2DSquareCellDoubleFileFactory();
                return geometricDensity(normalisationLevel,cellDistance,gridFactory);
            }
        }
    }
    public final AbstractGrid2DSquareCellDouble geometricDensity(int normalisationLevel,int cellDistance,AbstractGrid2DSquareCellDoubleFactory gridFactory) {
        if (DEBUG > 1) {System.out.println("geometricDensity(normalisationLevel("+normalisationLevel+"),cellDistance("+cellDistance+"),gridFactory("+gridFactory+")) {");}
        AbstractGridStatistics gridStatistics = getGridStatistics();
        AbstractGridStatistics outputGridStatistics;
        AbstractGridStatistics tempGridStatistics;
        double sparseness = gridStatistics.getSparseness();
        int nrows = getNrows();
        int ncols = getNcols();
        double xllcorner = getXllcorner();
        double yllcorner = getYllcorner();
        double cellsize = getCellsize();
        double noDataValue = getNoDataValue();
        AbstractGrid2DSquareCellDouble outputGrid,g2,g3,g4,g5,g6,g7;
        double d1 = 0.0d;
        double d2 = 0.0d;
        double d3 = 0.0d;
        int height = nrows;
        int width = ncols;
        int i1;
        int numberOfIterations = 0;
        int doubler = 1;
        int growth = 1;
        // Normalize input
        outputGrid = gridFactory.createGrid2DSquareCellDouble(nrows,ncols,xllcorner,yllcorner,cellsize,noDataValue);
        outputGridStatistics = outputGrid.getGridStatistics();
        if (sparseness == 1.0d) {return outputGrid;} // If all values are noDataValues the answer is g1
        g2 = gridFactory.createGrid2DSquareCellDouble(nrows,ncols,xllcorner,yllcorner,cellsize,0.0d);
        g3 = gridFactory.createGrid2DSquareCellDouble(nrows,ncols,xllcorner,yllcorner,cellsize,0.0d);
        double min;
        double max;
        if (normalisationLevel > 0) {
            min = gridStatistics.getMin();
            max = gridStatistics.getMax();
        } else {
            min = 0.0d;
            max = 1.0d;
        }
        for (int i = 0; i < nrows; i ++) {
            for (int j = 0; j < ncols; j ++) {
                d1 = getCell(i,j);
                if (d1 != noDataValue) {
                    // initcell() here is faster than setCell() which is faster than addToCell().
                    outputGrid.initCell(i,j,(d1 - min) / (max - min));
                    g2.initCell(i,j,(d1 - min) / (max - min));
                    g3.initCell(i,j,1.0d);
                }
            }
        }
        if (DEBUG > 2) {
            System.out.println("initial answer:");
            System.out.println(outputGrid.toString());
            System.out.println("numerator:");
            System.out.println(g2.toString());
            System.out.println("denominator:");
            System.out.println(g3.toString());
        }
        // Calculate number of iterations
        i1 = Math.min(cellDistance,(Math.max(nrows,ncols) / 2));
        for (int i = 0; i < cellDistance; i++) {
            if (i1 > 1) {
                i1 = i1 / 2;
                numberOfIterations ++;
            } else {break;}
        }
        // Densification loop
        for (int iteration = 0; iteration < numberOfIterations; iteration ++) {
            if (DEBUG > 1) {System.out.println("Iteration "+(iteration + 1)+" out of "+numberOfIterations);}
            // Set loop variables
            height += doubler;
            width += doubler;
            growth *= 2;
            // Step 1: Aggregate
            g4 = gridFactory.createGrid2DSquareCellDouble(height,width,0.0d,0.0d,cellsize,0.0d);
            g5 = gridFactory.createGrid2DSquareCellDouble(height,width,0.0d,0.0d,cellsize,0.0d);
            g6 = gridFactory.createGrid2DSquareCellDouble(height,width,0.0d,0.0d,cellsize,0.0d);
            for (int p = 0; p < doubler; p ++) {
                for (int q = 0; q < doubler; q ++) {
                    for (int i = 0; i < height; i += doubler) {
                        for (int j = 0; j < width; j += doubler) {
                            d1 = g2.getCell((i + p),(j + q)) + g2.getCell((i + p),(j + q - doubler)) + g2.getCell((i + p - doubler),(j + q)) + g2.getCell((i + p - doubler),(j + q - doubler));
                            g4.initCell((i + p),(j + q),d1);
                            d2 = g3.getCell((i + p),(j + q)) + g3.getCell((i + p),(j + q - doubler)) + g3.getCell((i + p - doubler),(j + q)) + g3.getCell((i + p - doubler),(j + q - doubler));
                            g5.initCell((i + p),(j + q),d2);
                            if (d2 != 0.0d) {
                                g6.initCell(i,j,(d1 / d2));
                            }
                        }
                    }
                }
            }
            //g2.clear();
            //g3.clear();
            if (DEBUG > 3) {
                System.out.println("aggregated numerator:");
                System.out.println(g4.toString());
                System.out.println("aggregated denominator:");
                System.out.println(g5.toString());
                System.out.println("division:");
                System.out.println(g6.toString());
            }
            // Step 2: Average over output region.
            // 1. This is probably the slowest part of the algorithm and gets slower
            //    with each iteration. Using alternative data structures and
            //    processing strategies this step can probably be speeded up a lot.
            g7 = gridFactory.createGrid2DSquareCellDouble(nrows,ncols,0.0d,0.0d,cellsize,0.0d);
            for (int i = 0; i < nrows; i += doubler) {
                for (int p = 0; p < doubler; p ++) {
                    for (int j = 0; j < ncols; j += doubler) {
                        for (int q = 0; q < doubler; q ++) {
                            d1 = 0.0d;
                            d2 = 0.0d;
                            for (int a = 0; a < growth; a ++) {
                                for (int b = 0; b < growth; b ++) {
                                    if (g6.inGrid((i + p + a),(j + q + b))) {
                                        d1 += g6.getCell((i + p + a),(j + q + b));
                                        d2 += 1.0d;
                                    }
                                }
                            }
                            if (d2 != 0.0d) {
                                g7.initCell((i + p),(j + q),(d1 / d2));
                            }
                        }
                    }
                }
            }
            System.out.println(g7.toString());
            System.out.println(outputGrid.toString());
            //g6.clear();
            // Add g7 to outputGrid
            if (normalisationLevel > 1) {
                tempGridStatistics = g7.getGridStatistics();
                min = tempGridStatistics.getMin();
                max = tempGridStatistics.getMax();
            } else {
                min = 0.0d;
                max = 1.0d;
            }
            for (int i = 0; i < nrows; i ++) {
                for (int j = 0; j < ncols; j ++) {
                    if (g7.inGrid(i,j)) {
                        outputGrid.addToCell(i,j,((g7.getCell(i,j) - min) / (max - min)));
                    }
                }
            }
            //g7.clear();
            if (DEBUG > 3) {
                System.out.println("result at iteration "+(iteration+1)+":");
                System.out.println(outputGrid.toString());
            }
            // Set loop variables for next iteration
            doubler *= 2;
            g2 = g4;
            g3 = g5;
        }
        // Mask
        outputGrid.mask(this);
        // Normalise
        if (normalisationLevel > 0) {
            max = outputGridStatistics.getMax();
            min = outputGridStatistics.getMin();
            for (int i = 0; i < nrows; i ++) {
                for (int j = 0; j < ncols; j ++) {
                    d3 = outputGrid.getCell(i,j);
                    if (d3 != noDataValue) {
                        outputGrid.setCell(i,j,((d3 - min) / (max - min)));
                    }
                }
            }
        }
        if (DEBUG > 2) {
            System.out.println("outputGrid:");
            System.out.println(outputGrid.toString());
        }
        if (DEBUG > 1) {System.out.println("}");}
        return outputGrid;
    }
    public final AbstractGrid2DSquareCellDouble geometricDensity(int normalisationLevel,int cellDistance,int chunkSize) {
        AbstractGrid2DSquareCellDoubleFactory gridFactory;
        try {
            gridFactory = new Grid2DSquareCellDoubleFactory();
            return geometricDensity(normalisationLevel,cellDistance,chunkSize,gridFactory);
        } catch (java.lang.OutOfMemoryError oome1) {
            gridFactory = new Grid2DSquareCellDoubleFileFactory();
            return geometricDensity(normalisationLevel,cellDistance,chunkSize,gridFactory);
        }
    }
    public final AbstractGrid2DSquareCellDouble geometricDensity(int normalisationLevel,int cellDistance,int chunkSize,AbstractGrid2DSquareCellDoubleFactory gridFactory) {
        if (DEBUG > 1) {System.out.println("geometricDensity(normalisationLevel("+normalisationLevel+"),cellDistance("+cellDistance+"),chunkSize("+chunkSize+"),gridFactory("+gridFactory+")) {");}
        String dataDirectory = System.getProperty("java.io.tmpdir");
        int nrows = getNrows();
        int ncols = getNcols();
        double xllcorner = getXllcorner();
        double yllcorner = getYllcorner();
        double cellsize = getCellsize();
        double noDataValue = getNoDataValue();
        int chunks = 0;
        int colChunks = 0;
        int rowChunks = 0;
        int startRow = 0;
        int endRow = 0;
        int startCol = 0;
        int endCol = 0;
        AbstractGrid2DSquareCellDouble outputGrid;
        int i1;
        int numberOfIterations,doubler,growth,lastrows,lastcols,thisrows,thiscols;
        double n1,d1,d2;
        int startRowIndex,endRowIndex,startColIndex,endColIndex;
        boolean alternator;
        // First try doing it by loading chunks into memory.
        if (DEBUG > 2) {System.out.println("First try doing it by loading chunks into memory...");}
        // Create outputGrid
        if (gridFactory instanceof Grid2DSquareCellDoubleFileFactory) {
            Grid2DSquareCellDoubleFileFactory factory = (Grid2DSquareCellDoubleFileFactory) gridFactory;
            if (this instanceof Grid2DSquareCellDoubleFile) {
                Grid2DSquareCellDoubleFile dummy = (Grid2DSquareCellDoubleFile) this;
                dataDirectory = dummy.getDataDirectory();
                factory.setDataDirectory(dataDirectory);
            }
            outputGrid = (Grid2DSquareCellDoubleFile) factory.createGrid2DSquareCellDouble(nrows,ncols,xllcorner,yllcorner,cellsize,noDataValue);
        } else {
            outputGrid = (Grid2DSquareCellDouble) gridFactory.createGrid2DSquareCellDouble(nrows,ncols,xllcorner,yllcorner,cellsize,noDataValue);
        }
        // Initialise temporary file f1
        String fileName1=new String("f"+System.currentTimeMillis());
        File f1=new File(dataDirectory,fileName1);
        try {
            if (f1.exists()) {
                boolean setUpFile=false;
                while (!setUpFile) {
                    fileName1="f"+System.currentTimeMillis();
                    f1=new File(dataDirectory,fileName1);
                    if (!f1.exists()) {
                        f1.createNewFile();
                        setUpFile=true;
                    }
                }
            } else {
                f1.createNewFile();
            }
        } catch (java.io.IOException ioe1) {
            System.out.println(""+ioe1+" creating temporary file in geometricDensity(dataDirectory("+dataDirectory+"),normalisationLevel("+normalisationLevel+"),cellDistance("+cellDistance+"),chunkSize("+chunkSize+"))");
            System.exit(0);
        }
        // Set to delete
        f1.deleteOnExit();
        try {
            RandomAccessFile rf1=new RandomAccessFile(f1,new String("rw"));
            // Split the problem into chunks
            if (chunkSize < (cellDistance * 3)) {
                chunkSize = cellDistance * 3;
            }
            // Calculate the number of the chunks of this size
            if ((ncols % chunkSize) > 0) {
                colChunks = (ncols / chunkSize) + 1;
            } else {
                colChunks = ncols / chunkSize;
            }
            if ((nrows % chunkSize) > 0) {
                rowChunks = (nrows / chunkSize) + 1;
            } else {
                rowChunks = nrows / chunkSize;
            }
            chunks  = rowChunks * colChunks;
            if (DEBUG > 2) {
                System.out.println("chunks to process="+chunks);
                System.out.println("row process chunks="+rowChunks);
                System.out.println("col process chunks="+colChunks);
                System.out.println("chunkSize="+chunkSize);
                System.out.println("cellDistance="+cellDistance);
            }
            // Set mallarchery loop variables
            startRow = (0 - cellDistance);
            endRow = (chunkSize - 1) + cellDistance;
            startCol = (0 - cellDistance - chunkSize);
            endCol = (cellDistance - 1);
            int extraCols = 0;
            int extraRows = 0;
            int tendRow = endRow;
            int tncols = 0;
            int tnrows = 0;
            for (int mallarchery = 0; mallarchery < chunks; mallarchery ++) {
                if (DEBUG > 1) {System.out.println("Chunk "+(mallarchery+1)+" out of "+chunks);}
                // Set mallarchery loop variables
                startCol += chunkSize;
                endCol += chunkSize;
                // Correct for short Col chunk
                if ((endCol - cellDistance) > (ncols - 1)) {
                    extraCols = endCol - ((ncols - 1) + cellDistance);
                    tncols = chunkSize - extraCols;
                    endCol = (ncols - 1) + cellDistance;
                } else {
                    extraCols = 0;
                    tncols = chunkSize;
                }
                // Correct for short Row chunks
                if ((endRow - cellDistance) > (nrows - 1)) {
                    tendRow = endRow;
                    extraRows = tendRow - ((nrows - 1) + cellDistance);
                    tnrows = chunkSize - extraRows;
                    endRow = (nrows - 1) + cellDistance;
                } else {
                    tendRow = endRow;
                    extraRows = 0;
                    tnrows = chunkSize;
                }
                if (DEBUG > 2) {
                    System.out.println("Values used:");
                    System.out.println("(nrows,ncols)=("+nrows+","+ncols+")");
                    System.out.println("(cellDistance,chunkSize)=("+cellDistance+","+chunkSize+")");
                    System.out.println("(startRow,startCol)=("+startRow+","+startCol+")");
                    System.out.println("(endRow,endCol)=("+endRow+","+endCol+")");
                    System.out.println("(extraRows,extraCols)=("+extraRows+","+extraCols+")");
                    System.out.println("(tnrows,tncols)=("+tnrows+","+tncols+")");
                }
                // Try to process chunk using fast access memory.
                //try {
                    if (DEBUG > 2) {System.out.println("Attempting to process chunk in memory...");}
                    // Construct input grid chunk.
                    Grid2DSquareCellDouble chunk = new Grid2DSquareCellDouble(this,startRow,startCol,endRow,endCol);
                    if (DEBUG > 3) {System.out.println(chunk.toString());}
                    // Generate geometricDensity
                    AbstractGrid2DSquareCellDouble chunkDensity = chunk.geometricDensity(normalisationLevel,cellDistance);
                    if (DEBUG > 3) {System.out.println(chunkDensity.toString());}
                    // Write out this chunks result to rf1.
                    chunkDensity.toBinaryFile(rf1,true,cellDistance,cellDistance,(cellDistance + tnrows - 1),(cellDistance + tncols - 1));
                /*} catch (java.lang.OutOfMemoryError oome1) {
                    // If that fails then send this chunk to be processed via the other geometricDensity() method.
                    if (DEBUG > 2) {System.out.println("...chunk too large to process in memory! Processing chunk using filespace...");}
                    File f3 = new File(dataDirectory,new String("f"+System.currentTimeMillis()));
                    Grid2DSquareCellDoubleFile chunk = new Grid2DSquareCellDoubleFile1_0(this,startRow,startCol,endRow,endCol,f3);
                    // Generate geometricDensiyy
                    Grid2DSquareCellDoubleFile chunkDensity = chunk.geometricDensity(dataDirectory,normalisationLevel,cellDistance);
                    // Write out this chunks result to rf1.
                    chunkDensity.toBinaryFile(rf1,true,cellDistance,cellDistance,cellDistance+(tnrows-1),cellDistance+(tncols-1));
                    // Tidy up
                    try {
                        chunk.getRandomAccessFile().close();
                        chunkDensity.getRandomAccessFile().close();
                    } catch (java.io.IOException ioe1) {
                        System.out.println(""+ioe1+" closing chunk.gridRandomAccessFile or chunkDensity.gridRandomAccessFile in geometricDensity(dataDirectory("+dataDirectory+"),cellDistance("+cellDistance+"),chunkSize("+chunkSize+"))");
                        System.exit(0);
                    }
                    chunk.getFile().delete();
                    //chunkDensity.grid.delete(); //one or the other!
                    f3.delete();
                }*/   
                // Set mallarchery loop variables
                // Set startRow and endRow for next iteration
                startCol += chunkSize;
                endCol += chunkSize;
                if ((startCol + cellDistance) > (ncols - 1)) {
                    startRow += chunkSize;
                    endRow += chunkSize;
                    startCol = (0 - cellDistance - chunkSize);
                    endCol = (cellDistance - 1);
                } else {
                    startCol -= chunkSize;
                    endCol -= chunkSize;
                    endRow = tendRow;
                }
            }
            // Read temporary file back in and write out result to g2.
            // This is effectively a row by row grid merge operation
            try {
                for (int p = 0; p < rowChunks; p ++) {
                    //System.out.println("rowChunk "+p);
                    for (int i = 0; i < chunkSize; i ++) {
                        for (int q = 0; q < colChunks; q ++) {
                            //System.out.println("colChunk "+q);
                            // if not beyond last row
                            if (((p * chunkSize) + i) < nrows) {
                                // The last rowChunk and the last colChunk are special cases
                                if ((p + 1) == rowChunks) {
                                    if ((q + 1) == colChunks) {
                                        rf1.seek(((p * (chunkSize * ((chunkSize * (colChunks - 1)) + tncols))) + (q * chunkSize * tnrows) + (i * tncols)) * 8);
                                    } else {
                                        rf1.seek(((p * (chunkSize * ((chunkSize * (colChunks - 1)) + tncols))) + (q * chunkSize * tnrows) + (i * chunkSize)) * 8);
                                    }
                                } else {
                                    if ((q + 1) == colChunks) {
                                        rf1.seek(((p * (chunkSize * ((chunkSize * (colChunks - 1)) + tncols))) + (q * chunkSize * chunkSize) + (i * tncols)) * 8);
                                    } else {
                                        rf1.seek(((p * (chunkSize * ((chunkSize * (colChunks - 1)) + tncols))) + (q * chunkSize * chunkSize) + (i * chunkSize)) * 8);
                                    }
                                }
                                for (int j = 0; j < chunkSize; j ++) {
                                    // if not beyond last column
                                    if (((q * chunkSize) + j) < ncols) {
                                        //g2.initCell(((int)p*chunkSize)+i,(int)(q*chunkSize)+j,rf1.readDouble());
                                        outputGrid.setCell(((p * chunkSize) + i),((q * chunkSize) + j),rf1.readDouble());
                                    }
                                }
                            }
                        }
                    }
                }
            } catch (java.io.IOException ioe1) {
                System.out.println(""+ioe1+" in geometricDensity(dataDirectory("+dataDirectory+"),cellDistance("+cellDistance+"),chunkSize("+chunkSize+"))");
                System.exit(0);
            }
            // Tidy up
            try {
                rf1.close();
            } catch (java.io.IOException ioe1) {
                System.out.println(""+ioe1+" closing rf1 in geometricDensity(dataDirectory("+dataDirectory+"),cellDistance("+cellDistance+"),chunkSize("+chunkSize+"))");
                System.exit(0);
            }
            f1.delete();
        } catch (java.io.FileNotFoundException fnfe1) {
            System.out.println(""+fnfe1+" in geometricDensity(dataDirectory("+dataDirectory+"),cellDistance("+cellDistance+"),chunkSize("+chunkSize+"))");
            System.exit(0);
        }
        if (DEBUG > 2) {System.out.println("}");}
        // Return result
        return outputGrid;
    }

    /**
     * Conversion methods
     */
    
    /**
     * Writes grid out to file in ArcInfo gridascii format and returns the File
     * to which it was written. The asciiFile is written with a 6 line header 
     * followed by getNrows() lines containing the values of the grid written 
     * out in column order. If the grid is already stored as a file then the 
     * asciiFile is written out in the same directory else to the temporary 
     * directory by default.
     */
    public File toAsciiFile() {
        String dataDirectory;
        if (this instanceof Grid2DSquareCellDoubleFile) {
            Grid2DSquareCellDoubleFile dummy = (Grid2DSquareCellDoubleFile) this;
            dataDirectory = dummy.getDataDirectory();
        } else {
            dataDirectory = System.getProperty("java.io.tmp.dir");
        }
        File f1 = new File(dataDirectory + new String("f"+System.currentTimeMillis()));
        try {
            while (!f1.exists()) {
                f1 = new File(dataDirectory + new String("f"+System.currentTimeMillis()));
            }
            f1.createNewFile();
        } catch (java.io.IOException ioe1) {
            System.out.println(""+ioe1+" creating file "+f1+"in toAsciiFile()");
            System.exit(0);
        }
        return toAsciiFile(f1);
    }
    public File toAsciiFile(File f1) {
        if (DEBUG > 2) {System.out.println("toAsciiFile(File("+f1+")) {");}
        PrintWriter pw1;
        try {
            pw1 = new PrintWriter(new BufferedWriter(new OutputStreamWriter(new FileOutputStream(f1))));
            // Write header info
            int nrows = getNrows();
            int ncols = getNcols();
            pw1.println("ncols "+getNcols());
            pw1.println("nrows "+getNrows());
            pw1.println("xllcorner "+getXllcorner());
            pw1.println("yllcorner "+getYllcorner());
            pw1.println("cellsize "+getCellsize());
            pw1.println("noDataValue "+getNoDataValue());
            // Write values
            for (int i = 0; i < nrows; i ++) {
                for (int j = 0 ; j < ncols; j ++) {
                    pw1.print(getCell(i,j)+" ");
                }
                pw1.println("");
            }
            // Close output
            pw1.flush();
            pw1.close();
        } catch (java.io.FileNotFoundException fnfe1) {
            System.out.println(""+fnfe1+" in toAsciiFile(File("+f1+"))");
            System.exit(0);
        }
        if (DEBUG > 2) {System.out.println("}");}
        return f1;
    }
    
    /**
     * Writes grid cell values in row major order as binary doubles to file using a
     * RandomAccessFile and returns the RandomAccessFile. If append is true then the
     * values are appended to the file. If append is false then the values are
     * written from the start of the file after which the number of rows, the number
     * of columns, xllcorner, yllcorner, cellsize and noDataValue are written and
     * the file length is set.
     */
    public RandomAccessFile toBinaryFile() {
        String dataDirectory;
        if (this instanceof Grid2DSquareCellDoubleFile) {
            Grid2DSquareCellDoubleFile dummy = (Grid2DSquareCellDoubleFile) this;
            dataDirectory = dummy.getDataDirectory();
        } else {
            dataDirectory = System.getProperty("java.io.tmp.dir");
        }
        File f1 = new File(dataDirectory + new String("f"+System.currentTimeMillis()));
        try {
            while (!f1.exists()) {
                f1 = new File(dataDirectory + new String("f"+System.currentTimeMillis()));
            }
            f1.createNewFile();
        } catch (java.io.IOException ioe1) {
            System.out.println(""+ioe1+" creating file in toBinaryFile()");
            System.exit(0);
        }
        return toBinaryFile(f1);
    }
    public RandomAccessFile toBinaryFile(File f1) {
        return toBinaryFile(f1,false);
    }
    public RandomAccessFile toBinaryFile(File f1,boolean append) {
        return toBinaryFile(f1,append,0,0,(getNrows() - 1),(getNcols() - 1));
    }
    public RandomAccessFile toBinaryFile(File f1,boolean append,int startRow,int startCol,int endRow,int endCol) {
        RandomAccessFile raf1;
        try {
            raf1 = new RandomAccessFile(f1,"rw");
            return toBinaryFile(raf1,append,0,0,(getNrows() - 1),(getNcols() - 1));
        } catch (java.io.IOException ioe1) {
            System.out.println(""+ioe1+" creating RandomAccessfile in toBinaryFile(File("+f1+")startRow("+startRow+"),startCol("+startCol+"),endRow("+endRow+"),endCol("+endCol+"))");
            System.exit(0);
        }
        return null;
    }
    public RandomAccessFile toBinaryFile(RandomAccessFile raf1) {
        return toBinaryFile(raf1,false);
    }
    public RandomAccessFile toBinaryFile(RandomAccessFile raf1,boolean append) {
        return toBinaryFile(raf1,append,0,0,(getNrows() - 1),(getNcols() - 1));
    }
    public RandomAccessFile toBinaryFile(RandomAccessFile raf1,boolean append,int startRow,int startCol,int endRow,int endCol) {
        if (DEBUG > 2) {System.out.println("toBinaryFile(RandomAccessFile(raf1),append("+append+"),startRow("+startRow+"),startCol("+startCol+"),endRow("+endRow+"),endCol("+endCol+"))");}
        try {
            if (append) {
                raf1.seek(raf1.length());
            } else {
                raf1.seek(0);
            }
            for (int i = startRow; i <= endRow; i ++) {
                for (int j = startCol; j <= endCol; j ++) {
                    raf1.writeDouble(getCell(i,j));
                }
            }
            if (!append) {
                raf1.writeDouble((double) (endRow - startRow + 1));
                raf1.writeDouble((double) (endCol - startCol + 1));
                raf1.writeDouble(getXllcorner());
                raf1.writeDouble(getYllcorner());
                raf1.writeDouble(getCellsize());
                raf1.writeDouble(getNoDataValue());
                raf1.setLength(raf1.getFilePointer());
            }
        } catch (java.io.IOException ioe1) {
            System.out.println(""+ioe1+" in toBinaryFile(RandomAccessFile(raf1),append("+append+"),startRow("+startRow+"),startCol("+startCol+"),endRow("+endRow+"),endCol("+endCol+"))");
            System.exit(0);
        }
        if (DEBUG > 2) {System.out.println("}");}
        return raf1;
    }
    
    /**
     * Returns a uk.ac.leeds.ccg.raster.Raster constructed from an array of grid
     * cell values.
     */
    public Raster toRaster() {
        if (DEBUG > 2) {System.out.println("toRaster() {");}
        int nrows = getNrows();
        int ncols = getNcols();
        double cellsize = getCellsize();
        GeoRectangle gr1 = new GeoRectangle(getXllcorner(),getYllcorner(),(double)(ncols*cellsize),(double)(nrows*cellsize));
        Raster r1 = new Raster(gr1,cellsize);
        try {
            double[] data = new double[nrows*ncols];
            double d1;
            for (int i = 0; i < nrows; i ++) {
                for (int j = 0; j < ncols; j ++) {
                    d1 = getCell(i,j);
                    data[(i * ncols) + j] = d1;
                    //data[getCellID(i,j)] = d1;
                }
            }
            //r1.setNoDataValue(getNoDataValue());
            r1.setMissingValueCode(getNoDataValue());
            r1.setData(data);
        } catch (OutOfMemoryError oome1) {
            System.out.println(""+oome1+" in toRaster()");
            System.exit(0);
        }
        if (DEBUG > 2) {System.out.println("}");}
        return r1;
    }

}