Landscape.cpp

/*----------------------------------------------------------------------------
 *
 *	Copyright (C) 2026 Greta Bocedi, Stephen C.F. Palmer, Justin M.J. Travis, Anne-Kathleen Malchow, Roslyn Henry, Théo Pannetier, Jette Wolff, Damaris Zurell
 *
 *	This file is part of RangeShifter.
 *
 *	RangeShifter is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation, either version 3 of the License, or
 *	(at your option) any later version.
 *
 *	RangeShifter is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 *	GNU General Public License for more details.
 *
 *	You should have received a copy of the GNU General Public License
 *	along with RangeShifter. If not, see <https://www.gnu.org/licenses/>.
 *
 --------------------------------------------------------------------------*/


 //---------------------------------------------------------------------------

#include "Landscape.h"
//---------------------------------------------------------------------------

ifstream landscape;

ofstream outConnMat;
ofstream outvisits;
#if RS_RCPP
ofstream outMovePaths;
#endif // RS_RCPP

//---------------------------------------------------------------------------

// Initial species distribution functions

InitDist::InitDist(Species* pSp)
{
	pSpecies = pSp;
	resol = 1;
	maxX = 0;
	maxY = 0;
	minEast = 0.0;
	minNorth = 0.0;
}

InitDist::~InitDist() {
	int ncells = (int)cells.size();
	for (int i = 0; i < ncells; i++)
		if (cells[i] != NULL) delete cells[i];
	cells.clear();
}

void InitDist::setDistribution(int nInit) {
	int rr = 0;
	int ncells = (int)cells.size();
	if (nInit == 0) { // set all cells to be initialised
		for (int i = 0; i < ncells; i++) {
			cells[i]->setCell(true);
		}
	}
	else { // set specified number of cells at random to be initialised
		if (nInit > ncells / 2) { // use backwards selection method
			for (int i = 0; i < ncells; i++) cells[i]->setCell(true);
			for (int i = 0; i < (ncells - nInit); i++) {
				do {
					rr = pRandom->IRandom(0, ncells - 1);
				} while (!cells[rr]->selected());
				cells[rr]->setCell(false);
			}
		}
		else { // use forwards selection method
			for (int i = 0; i < ncells; i++) cells[i]->setCell(false);
			for (int i = 0; i < nInit; i++) {
				do {
					rr = pRandom->IRandom(0, ncells - 1);
				} while (cells[rr]->selected());
				cells[rr]->setCell(true);
			}
		}
	}
}

// Set a specified cell (by position in cells vector)
void InitDist::setDistCell(int ix, bool init) {
	cells[ix]->setCell(init);
}

// Set a specified cell (by co-ordinates)
void InitDist::setDistCell(locn loc, bool init) {
	locn cellloc;
	int ncells = (int)cells.size();
	for (int i = 0; i < ncells; i++) {
		cellloc = cells[i]->getLocn();
		if (cellloc.x == loc.x && cellloc.y == loc.y) {
			cells[i]->setCell(init);
			i = ncells;
		}
	}
}

// Specified location is within the initial distribution?
bool InitDist::inInitialDist(locn loc) {
	int ncells = (int)cells.size();
	for (int i = 0; i < ncells; i++) {
		if (cells[i]->toInitialise(loc)) { // cell is to be initialised
			return true;
		}
	}
	return false;
}

int InitDist::cellCount(void) {
	return (int)cells.size();
}

// Return the co-ordinates of a specified initial distribution cell
locn InitDist::getCell(int ix) {
	locn loc;
	if (ix >= 0 && ix < (int)cells.size()) {
		loc = cells[ix]->getLocn();
	}
	else {
		loc.x = loc.y = -666; // indicates invalid index specified
	}
	return loc;
}

// Return the co-ordinates of a specified initial distribution cell if it has been
// selected - otherwise return negative co-ordinates
locn InitDist::getSelectedCell(int ix) {
	locn loc; loc.x = loc.y = -666;
	if (ix < (int)cells.size()) {
		if (cells[ix]->selected()) {
			loc = cells[ix]->getLocn();
		}
	}
	return loc;
}

locn InitDist::getDimensions(void) {
	locn d; d.x = maxX; d.y = maxY; return d;
}

void InitDist::resetDistribution(void) {
	int ncells = (int)cells.size();
	for (int i = 0; i < ncells; i++) {
		cells[i]->setCell(false);
	}
}

//---------------------------------------------------------------------------

// Read species initial distribution file

int InitDist::readDistribution(string distfile) {
#if RS_RCPP
	wstring header;
#else
	string header;
#endif
	int p, nodata;
	int ncols, nrows;
#if RS_RCPP
	wifstream dfile; // species distribution file input stream
#else
	ifstream dfile; // species distribution file input stream
#endif

	// open distribution file
	dfile.open(distfile.c_str());
	if (!dfile.is_open()) return 21;

// read landscape data from header records of distribution file
// NB headers of all files have already been compared
double tmpresol;
dfile >> header >> ncols >> header >> nrows >> header >> minEast >> header >> minNorth
	>> header >> tmpresol >> header >> nodata;
resol = (int) tmpresol;
#if RS_RCPP
	if (!dfile.good()) {
		// corrupt file stream
		StreamErrorR(distfile);
		dfile.close();
		dfile.clear();
		return 144;
	}
#endif

	maxX = ncols - 1; maxY = nrows - 1;

	// set up bad integer value to ensure that valid values are read
	int badvalue = -9; if (nodata == -9) badvalue = -99;

	for (int y = nrows - 1; y >= 0; y--) {
		for (int x = 0; x < ncols; x++) {
			p = badvalue;
#if RS_RCPP
			if (dfile >> p) {
#else
			dfile >> p;
#endif
			if (p == nodata || p == 0 || p == 1) { // only valid values
				if (p == 1) { // species present
					cells.push_back(new DistCell(x, y));
				}
			}
			else { // error in file
				dfile.close(); dfile.clear();
				return 22;
			}
#if RS_RCPP
			}
	else {
		// corrupt file stream
#if RS_RCPP && !R_CMD
		Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
		StreamErrorR(distfile);
		dfile.close();
		dfile.clear();
		return 144;
	}
#endif
		}
	}
#if RS_RCPP
dfile >> p;
if (!dfile.eof()) EOFerrorR(distfile);
#endif

	dfile.close(); dfile.clear();
	return 0;
}

// Read species initial distribution file for threadsafe option
#if RS_RCPP
int InitDist::readDistribution(Rcpp::NumericMatrix distfile, landOrigin habfile_origin, int spResol) {

    int d=0;
    double dfloat=0;
    int ncols,nrows;

    ncols = distfile.ncol();
    nrows = distfile.nrow();

    minEast = habfile_origin.minEast;
    minNorth = habfile_origin.minNorth;
    resol = spResol;
    maxX = ncols-1;
    maxY = nrows-1;

    for (int y = nrows-1; y >= 0; y--) {
        for (int x = 0; x < ncols; x++) {

            dfloat = distfile(nrows-1-y,x);

            if ( !R_IsNA(dfloat) ){ // check for NA
                d = (int)dfloat;
                if ( d == 0 || d == 1) { // only valid values
                    if (d == 1) { // species present
                        cells.push_back(new DistCell(x,y));
                    }
                }
                else { // error in file
#if RS_RCPP && !R_CMD
                    Rcpp::Rcout << "Found invalid value in species distribution raster." <<  std::endl;
#endif
                    return 22;
                }
            }
        }
    }
    return 0;
}

#endif
//---------------------------------------------------------------------------

// Landscape functions

Landscape::Landscape(void) {
	patchModel = false; spDist = false; generated = false; fractal = false; continuous = false;
	dynamic = false; habIndexed = false;
	spatialdemog = false;
	resol = spResol = 1; landNum = 0;
	rasterType = 0;
	nHab = nHabMax = 0;
	dimX = dimY = 100;
	minX = minY = 0;
	maxX = maxY = 99;
	minPct = maxPct = propSuit = hurst = 0.0;
	maxCells = 100;
	gpix = 1.0;
	pix = (int)gpix;
	minEast = minNorth = 0.0;
	cells = 0;
	connectMatrix = 0;
	epsGlobal = 0;
	patchChgMatrix = 0;
	costsChgMatrix = 0;
}

Landscape::~Landscape() {

	if (cells != 0) {
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				if (cells[y][x] != 0) delete cells[y][x];
			}
			if (cells[y] != 0) {
				delete[] cells[y];
			}
		}
		delete[] cells;
		cells = 0;
	}

	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++)
		if (patches[i] != NULL) delete patches[i];
	patches.clear();

	int ndistns = (int)distns.size();
	for (int i = 0; i < ndistns; i++)
		if (distns[i] != NULL) delete distns[i];
	distns.clear();

	int ninitcells = (int)initcells.size();
	for (int i = 0; i < ninitcells; i++)
		if (initcells[i] != NULL) delete initcells[i];
	initcells.clear();

	habCodes.clear();
	landchanges.clear();
	patchchanges.clear();

	deleteConnectMatrix();
	deletePatchChgMatrix();
	if (epsGlobal != 0) delete[] epsGlobal;
}

// Remove all patches and cells
// Used for replicating artificial landscape without deleting the landscape itself
void Landscape::resetLand(void) {
	resetLandLimits();
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) if (patches[i] != NULL) delete patches[i];
	patches.clear();
	if (cells != 0) {
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				if (cells[y][x] != 0) delete cells[y][x];
			}
			if (cells[y] != 0) {
				delete[] cells[y];
			}
		}
		delete[] cells;
		cells = 0;
	}
}

void Landscape::setLandParams(landParams ppp, bool batchMode)
{
	generated = ppp.generated;
	patchModel = ppp.patchModel;
	spDist = ppp.spDist;
	dynamic = ppp.dynamic;
	landNum = ppp.landNum;
	spatialdemog = ppp.spatialdemog;
	if (ppp.resol > 0) resol = ppp.resol;
	if (ppp.spResol > 0 && ppp.spResol % ppp.resol == 0) spResol = ppp.spResol;
	if ((ppp.rasterType >= 0 && ppp.rasterType <= 2) || ppp.rasterType == 9)
		rasterType = ppp.rasterType;
	if (ppp.nHab >= 1) nHab = ppp.nHab;
	if (ppp.nHabMax >= 1) nHabMax = ppp.nHabMax;
	if (ppp.dimX > 0) dimX = ppp.dimX;
	if (ppp.dimY > 0) dimY = ppp.dimY;
	if (ppp.minX >= 0 && ppp.maxX >= 0 && ppp.minX <= ppp.maxX && ppp.maxX < dimX) {
		minX = ppp.minX; maxX = ppp.maxX;
	}
	else {
		minX = 0; maxX = dimX - 1;
	}
	if (ppp.minY >= 0 && ppp.maxY >= 0 && ppp.minY <= ppp.maxY && ppp.maxY < dimY) {
		minY = ppp.minY; maxY = ppp.maxY;
	}
	else {
		minY = 0; maxY = dimY - 1;
	}
	if (batchMode && rasterType == 0) {
		// in batch mode, set up sequential habitat codes if not already present
		if (habCodes.size() == 0) {
			for (int i = 0; i < nHabMax; i++) {
				habCodes.push_back(i + 1);
			}
		}
	}
}

landParams Landscape::getLandParams(void)
{
	landParams ppp;
	ppp.generated = generated; ppp.patchModel = patchModel; ppp.spDist = spDist;
	ppp.dynamic = dynamic;
	ppp.spatialdemog = spatialdemog;
	ppp.landNum = landNum;
	ppp.resol = resol; ppp.spResol = spResol;
	ppp.rasterType = rasterType;
	ppp.nHab = nHab; ppp.nHabMax = nHabMax;
	ppp.dimX = dimX; ppp.dimY = dimY;
	ppp.minX = minX; ppp.minY = minY;
	ppp.maxX = maxX; ppp.maxY = maxY;
	return ppp;
}

landData Landscape::getLandData(void) {
	landData dd;
	dd.resol = resol;
	dd.dimX = dimX; dd.dimY = dimY;
	dd.minX = minX; dd.minY = minY;
	dd.maxX = maxX; dd.maxY = maxY;
	return dd;
}

void Landscape::setGenLandParams(genLandParams ppp)
{
	fractal = ppp.fractal;
	continuous = ppp.continuous;
	if (ppp.minPct > 0.0 && ppp.minPct < 100.0) minPct = ppp.minPct;
	if (ppp.maxPct > 0.0 && ppp.maxPct <= 100.0) maxPct = ppp.maxPct;
	if (ppp.propSuit >= 0.0 && ppp.propSuit <= 1.0) propSuit = ppp.propSuit;
	if (ppp.hurst > 0.0 && ppp.hurst < 1.0) hurst = ppp.hurst;
	if (ppp.maxCells > 0) maxCells = ppp.maxCells;
}

genLandParams Landscape::getGenLandParams(void)
{
	genLandParams ppp;
	ppp.fractal = fractal; ppp.continuous = continuous;
	ppp.minPct = minPct; ppp.maxPct = maxPct; ppp.propSuit = propSuit; ppp.hurst = hurst;
	ppp.maxCells = maxCells;
	return ppp;
}

void Landscape::setLandLimits(int x0, int y0, int x1, int y1) {
	if (x0 >= 0 && x1 >= 0 && x0 <= x1 && x1 < dimX
		&& y0 >= 0 && y1 >= 0 && y0 <= y1 && y1 < dimY) {
		minX = x0; maxX = x1; minY = y0; maxY = y1;
	}
}

void Landscape::resetLandLimits(void) {
	minX = minY = 0; maxX = dimX - 1; maxY = dimY - 1;
}

//---------------------------------------------------------------------------

void Landscape::setLandPix(landPix p) {
	if (p.pix > 0) pix = p.pix;
	if (p.gpix > 0.0) gpix = p.gpix;
}

landPix Landscape::getLandPix(void) {
	landPix p;
	p.pix = pix; p.gpix = gpix;
	return p;
}

void Landscape::setOrigin(landOrigin origin) {
	minEast = origin.minEast; minNorth = origin.minNorth;
}

landOrigin Landscape::getOrigin(void) {
	landOrigin origin;
	origin.minEast = minEast; origin.minNorth = minNorth;
	return origin;
}

//---------------------------------------------------------------------------

// Functions to handle habitat codes

bool Landscape::habitatsIndexed(void) { return habIndexed; }

void Landscape::listHabCodes(void) {
	int nhab = (int)habCodes.size();
#if RS_RCPP && !R_CMD
	Rcpp::Rcout << endl;
	for (int i = 0; i < nhab; i++) {
		Rcpp::Rcout << "Habitat code[ " << i << "] = " << habCodes[i] << endl;
	}
	Rcpp::Rcout << endl;
#else
	cout << endl;
	for (int i = 0; i < nhab; i++) {
		cout << "Habitat code[ " << i << "] = " << habCodes[i] << endl;
	}
	cout << endl;
#endif
}

void Landscape::addHabCode(int hab) {
	int nhab = (int)habCodes.size();
	bool addCode = true;
	for (int i = 0; i < nhab; i++) {
		if (hab == habCodes[i]) {
			addCode = false; i = nhab + 1;
		}
	}
	if (addCode) { habCodes.push_back(hab); nHab++; }
}

// Get the index number of the specified habitat in the habitats vector
int Landscape::findHabCode(int hab) {
	int nhab = (int)habCodes.size();
	for (int i = 0; i < nhab; i++) {
		if (hab == habCodes[i]) return i;
	}
	return -999;
}

// Get the specified habitat code
int Landscape::getHabCode(int ixhab) {
	if (ixhab < (int)habCodes.size()) return habCodes[ixhab];
	else return -999;
}

void Landscape::clearHabitats(void) {
	habCodes.clear();
}

//---------------------------------------------------------------------------
void Landscape::setCellArray(void) {
	if (cells != 0) resetLand();
	//cells = new Cell **[maxY+1];
	cells = new Cell * *[dimY];
	for (int y = dimY - 1; y >= 0; y--) {
		cells[y] = new Cell * [dimX];
		for (int x = 0; x < dimX; x++) {
			cells[y][x] = 0;
		}
	}
}

//---------------------------------------------------------------------------
/* Create an artificial landscape (random or fractal), which can be
either binary (habitat index 0 is the matrix, 1 is suitable habitat)
or continuous (0 is the matrix, >0 is suitable habitat) */
void Landscape::generatePatches()
{
	int x, y, ncells;
	double p;
	Patch* pPatch;
	Cell* pCell;

	vector <land> ArtLandscape;
	setCellArray();

	int patchnum = 0;  // initial patch number for cell-based landscape
	// create patch 0 - the matrix patch (even if there is no matrix)
	newPatch(patchnum++);

	// as landscape generator returns cells in a random sequence, first set up all cells
	// in the landscape in the correct sequence, then update them and create patches for
	// habitat cells
	for (int yy = dimY - 1; yy >= 0; yy--) {
		for (int xx = 0; xx < dimX; xx++) {
			addNewCellToLand(xx, yy, 0);
		}
	}

	if (continuous) rasterType = 2;
	else rasterType = 0;
	if (fractal) {
		p = 1.0 - propSuit;
		// fractal_landscape() requires Max_prop > 1 (but does not check it!)
		// as in turn it calls runif(1.0,Max_prop)
		double maxpct;
		if (maxPct < 1.0) maxpct = 100.0; else maxpct = maxPct;

		ArtLandscape = fractal_landscape(dimY, dimX, hurst, p, maxpct, minPct);

		vector<land>::iterator iter = ArtLandscape.begin();
		while (iter != ArtLandscape.end()) {
			x = iter->y_coord; y = iter->x_coord;
			pCell = findCell(x, y);
			if (continuous) {
				if (iter->value > 0.0) { // habitat
					pPatch = newPatch(patchnum++);
					addCellToPatch(pCell, pPatch, iter->value);
				}
				else { // matrix
					addCellToPatch(pCell, patches[0], iter->value);
				}
			}
			else { // discrete
				if (iter->avail == 0) { // matrix
					addCellToPatch(pCell, patches[0]);
				}
				else { // habitat
					pPatch = newPatch(patchnum++);
					addCellToPatch(pCell, pPatch);
					pCell->changeHabIndex(0, 1);
				}
			}
			iter++;
		}
	}
	else { // random landscape
		int hab = 0;
		ncells = (int)((float)(dimX) * (float)(dimY)*propSuit + 0.00001); // no. of cells to initialise
		int i = 0;
		do {
			do {
				x = pRandom->IRandom(0, dimX - 1); y = pRandom->IRandom(0, dimY - 1);
				pCell = findCell(x, y);
				hab = pCell->getHabIndex(0);
			} while (hab > 0);
			pPatch = newPatch(patchnum++);
			pCell = findCell(x, y);
			addCellToPatch(pCell, pPatch);
			pCell->changeHabIndex(0, 1);
			if (continuous) {
				pCell->setHabitat((float)(minPct + pRandom->Random() * (maxPct - minPct)));
			}
			i++;
		} while (i < ncells);
		// remaining cells need to be added to the matrix patch
		p = 0.0;
		x = 0;
		for (int yy = dimY - 1; yy >= 0; yy--) {
			for (int xx = 0; xx < dimX; xx++) {
				pCell = findCell(xx, yy);
				if (continuous) {
					if (pCell->getHabitat(0) <= 0.0)
					{
						addCellToPatch(pCell, patches[0], (float)p);
					}
				}
				else { // discrete
					if (pCell->getHabIndex(0) == 0) {
						addCellToPatch(pCell, patches[0], x);
					}
				}
			}
		}
	}
}

//---------------------------------------------------------------------------

// Landscape patch-management functions

//---------------------------------------------------------------------------
/* Create a patch for each suitable cell of a cell-based landscape (all other
habitat cells are added to the matrix patch) */
void Landscape::allocatePatches(Species* pSpecies)
{
	float habK;
	Patch* pPatch;
	Cell* pCell;

	// delete all existing patches
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) {
		if (patches[i] != NULL) delete patches[i];
	}
	patches.clear();
	// create the matrix patch
	patches.push_back(new Patch(0, 0));
	Patch* matrixPatch = patches[0];
	int patchnum = 1;

	switch (rasterType) {

	case 0: // habitat codes
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				if (cells[y][x] != 0) { // not no-data cell
					pCell = cells[y][x];
					habK = 0.0;
					int nhab = pCell->nHabitats();
					for (int i = 0; i < nhab; i++) {
						habK += pSpecies->getHabK(pCell->getHabIndex(i));
					}
					if (habK > 0.0) { // cell is suitable - create a patch for it
						pPatch = newPatch(patchnum++);
						addCellToPatch(pCell, pPatch);
					}
					else { // cell is not suitable - add to the matrix patch
						addCellToPatch(pCell, matrixPatch);
						pPatch = 0;
					}
				}
			}
		}
		break;
	case 1: // habitat cover
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				if (cells[y][x] != 0) { // not no-data cell
					pCell = cells[y][x];
					habK = 0.0;
					int nhab = pCell->nHabitats();
					for (int i = 0; i < nhab; i++)
					{
						habK += pSpecies->getHabK(i) * pCell->getHabitat(i) / 100.0f;
					}
					if (habK > 0.0) { // cell is suitable - create a patch for it
						pPatch = newPatch(patchnum++);
						addCellToPatch(pCell, pPatch);
					}
					else { // cell is not suitable - add to the matrix patch
						addCellToPatch(pCell, matrixPatch);
						pPatch = 0;
					}
				}
			}
		}
		break;
	case 2: // habitat quality
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				if (cells[y][x] != 0) { // not no-data cell
					pCell = cells[y][x];
					habK = 0.0;
					int nhab = pCell->nHabitats();
					for (int i = 0; i < nhab; i++)
					{
						habK += pSpecies->getHabK(0) * pCell->getHabitat(i) / 100.0f;
					}
					if (habK > 0.0) { // cell is suitable (at some time) - create a patch for it
						pPatch = newPatch(patchnum++);
						addCellToPatch(pCell, pPatch);
					}
					else { // cell is never suitable - add to the matrix patch
						addCellToPatch(pCell, matrixPatch);
						pPatch = 0;
					}
				}
			}
		}
		break;

	} // end of switch (rasterType)
}

Patch* Landscape::newPatch(int num)
{
	int npatches = (int)patches.size();
	patches.push_back(new Patch(num, num));
	return patches[npatches];
}

Patch* Landscape::newPatch(int seqnum, int num)
{
	int npatches = (int)patches.size();
	patches.push_back(new Patch(seqnum, num));
	return patches[npatches];
}

void Landscape::resetPatches(void) {
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) {
		patches[i]->resetLimits();
	}
}

void Landscape::addNewCellToLand(int x, int y, float q) {
	if (q < 0.0) // no-data cell - no Cell created
		cells[y][x] = 0;
	else
		cells[y][x] = new Cell(x, y, 0, q);
}

void Landscape::addNewCellToLand(int x, int y, int hab) {
	if (hab < 0) // no-data cell - no Cell created
		cells[y][x] = 0;
	else
		cells[y][x] = new Cell(x, y, 0, hab);
}

void Landscape::addCellToLand(Cell* c) {
	if (cells == 0) throw runtime_error("Landscape cells member is uninitialised.");
	if (c->getHabIndex(0) < 0.0)
		throw logic_error("Can't add no-data cell to landscape.");
	locn l = c->getLocn();
	cells[l.y][l.x] = c;
}

void Landscape::addNewCellToPatch(Patch* pPatch, int x, int y, float q) {
	if (q < 0.0) { // no-data cell - no Cell created
		cells[y][x] = 0;
	}
	else { // create the new cell
		cells[y][x] = new Cell(x, y, pPatch, q);
		if (pPatch != 0) { // not the matrix patch
			// add the cell to the patch
			pPatch->addCell(cells[y][x], x, y);
		}
	}
}

void Landscape::addNewCellToPatch(Patch* pPatch, int x, int y, int hab) {
	if (hab < 0) // no-data cell - no Cell created
		cells[y][x] = 0;
	else { // create the new cell
		cells[y][x] = new Cell(x, y, pPatch, hab);
		if (pPatch != 0) { // not the matrix patch
			// add the cell to the patch
			pPatch->addCell(cells[y][x], x, y);
		}
	}
}

void Landscape::addCellToPatch(Cell* pCell, Patch* pPatch) {
	pCell->setPatch(pPatch);
	locn loc = pCell->getLocn();
	// add the cell to the patch
	pPatch->addCell(pCell, loc.x, loc.y);
}

void Landscape::addCellToPatch(Cell* pCell, Patch* pPatch, float q) {
	pCell->setPatch(pPatch);
	// update the habitat type of the cell
	pCell->setHabitat(q);
	locn loc = pCell->getLocn();
	// add the cell to the patch
	pPatch->addCell(pCell, loc.x, loc.y);
}

void Landscape::addCellToPatch(Cell* pCell, Patch* pPatch, int hab) {
	pCell->setPatch(pPatch);
	// update the habitat type of the cell
	pCell->setHabIndex(hab);
	locn loc = pCell->getLocn();
	// add the cell to the patch
	pPatch->addCell(pCell, loc.x, loc.y);
}

patchData Landscape::getPatchData(int ix) {
	patchData ppp;
	ppp.pPatch = patches[ix]; ppp.patchNum = patches[ix]->getPatchNum();
	ppp.nCells = patches[ix]->getNCells();
	locn randloc; randloc.x = -666; randloc.y = -666;
	Cell* pCell = patches[ix]->getRandomCell();
	if (pCell != 0) {
		randloc = pCell->getLocn();
	}
	ppp.x = randloc.x; ppp.y = randloc.y;
	return ppp;
}

bool Landscape::existsPatch(int num) {
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) {
		if (num == patches[i]->getPatchNum()) return true;
	}
	return false;
}

Patch* Landscape::findPatch(int num) {
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) {
		if (num == patches[i]->getPatchNum()) return patches[i];
	}
	return 0;
}


set<int> Landscape::getPatchNbs() const {
	set<int> patchNbs;
	for (auto& p : patches)
		patchNbs.emplace(p->getPatchNum());
	return patchNbs;
}

set<int> Landscape::samplePatches(const string& samplingOption, int nbToSample, Species* pSpecies) {

	vector<int> sampledPatches;
	vector<int> eligiblePatches;

	// Get list of viable patches where the species is present
	for (auto p : patches) {
		if (p->getPatchNum() == 0) continue; // skip patch 0, the matrix
		if (samplingOption == "random" // then all patches are eligible
			|| p->speciesIsPresent(pSpecies)) // otherwise only patches with at least 1 ind
			eligiblePatches.push_back(p->getPatchNum());
	}

	if (samplingOption == "all") {
		sampledPatches = eligiblePatches;
	}
	else if (samplingOption == "random_occupied" || samplingOption == "random") {
		if (nbToSample > eligiblePatches.size())
			nbToSample = eligiblePatches.size();
		auto rng = pRandom->getRNG();
		sample(eligiblePatches.begin(), eligiblePatches.end(), std::back_inserter(sampledPatches),
			nbToSample, rng);
	}
	else {
		throw logic_error("Sampling option should be random, random_occupied or all when sampling patches.");
	}

	set<int> patchIds;
	copy(sampledPatches.begin(), sampledPatches.end(), inserter(patchIds, patchIds.end()));
	return patchIds;
}

void Landscape::resetPatchPopns(void) {
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) {
		patches[i]->resetPopn();
	}
}

void Landscape::updateCarryingCapacity(Species* pSpecies, int yr, short landIx) {
	envGradParams grad = paramsGrad->getGradient();
	bool gradK = false;
	if (grad.gradient && grad.gradType == 1) gradK = true; // gradient in carrying capacity
	patchLimits landlimits;
	landlimits.xMin = minX; landlimits.xMax = maxX;
	landlimits.yMin = minY; landlimits.yMax = maxY;
	int npatches = (int)patches.size();
	for (int i = 0; i < npatches; i++) {
		if (patches[i]->getPatchNum() != 0) { // not matrix patch
			patches[i]->setCarryingCapacity(pSpecies, landlimits,
				getGlobalStoch(yr), nHab, rasterType, landIx, gradK);
		}
	}

}

void Landscape::updateDemoScalings(short landIx) {

	patchLimits landlimits;
	landlimits.xMin = minX; landlimits.xMax = maxX;
	landlimits.yMin = minY; landlimits.yMax = maxY;

	if(spatialdemog && rasterType == 2) {// demographic scaling only implemented for habitat quality maps
		int npatches = (int)patches.size(); // new: for (auto& p : patches)
		for (int i = 0; i < npatches; i++) {
			if (patches[i]->getPatchNum() != 0) { // not matrix patch
				// calculate local scaling for each patch from its constituent cells
				patches[i]->setPatchDemoScaling(landIx, landlimits);
			}
		}
	}
}

Cell* Landscape::findCell(int x, int y) {
	if (x >= 0 && x < dimX && y >= 0 && y < dimY) return cells[y][x];
	else return 0;
}

bool Landscape::checkDataCell(int x, int y) {
    Cell* pCell;
    pCell = findCell(x, y);
    return true;
}


int Landscape::patchCount(void) {
	return (int)patches.size();
}

void Landscape::listPatches(void) {
	patchLimits p;
	int npatches = (int)patches.size();
#if RS_RCPP && !R_CMD
	Rcpp::Rcout << endl;
	for (int i = 0; i < npatches; i++) {
		p = patches[i]->getLimits();
		Rcpp::Rcout << "Patch " << patches[i]->getPatchNum()
			<< " xMin = " << p.xMin << " xMax = " << p.xMax
			<< " \tyMin = " << p.yMin << " yMax = " << p.yMax
			<< endl;
	}
	Rcpp::Rcout << endl;
#else
	cout << endl;
	for (int i = 0; i < npatches; i++) {
		p = patches[i]->getLimits();
		cout << "Patch " << patches[i]->getPatchNum()
			<< " xMin = " << p.xMin << " xMax = " << p.xMax
			<< " \tyMin = " << p.yMin << " yMax = " << p.yMax
			<< endl;
	}
	cout << endl;
#endif
}

// Check that total cover of any cell does not exceed 100%
// and identify matrix cells
int Landscape::checkTotalCover(void) {
	if (rasterType != 1) return 0; // not appropriate test
	int nCells = 0;
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] != 0)
			{ // not a no-data cell
				float sumCover = 0.0;
				for (int i = 0; i < nHab; i++) {
					sumCover += cells[y][x]->getHabitat(i);
				}
				if (sumCover > 100.00001) nCells++; // decimal part to allow for floating point error
				if (sumCover <= 0.0) // cell is a matrix cell
					cells[y][x]->setHabIndex(0);
				else
					cells[y][x]->setHabIndex(1);
			}
		}
	}
	return nCells;
}

// Convert habitat codes stored on loading habitat codes landscape to
// sequential sorted index numbers
void Landscape::updateHabitatIndices(void) {
	// sort codes
	sort(habCodes.begin(), habCodes.end());
	nHab = (int)habCodes.size();
	// convert codes in landscape
	int h;
	int changes = (int)landchanges.size();
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] != 0) { // not a no-data cell
				for (int c = 0; c <= changes; c++) {
					h = cells[y][x]->getHabIndex(c);
					if (h >= 0) {
						h = findHabCode(h);
						cells[y][x]->changeHabIndex(c, h);
					}
				}
			}
		}
	}
	habIndexed = true;
}

void Landscape::setEnvGradient(Species* pSpecies, bool initial)
{
	float dist_from_opt, dev;
	float habK;
	//int hab;
	double envval;
	// gradient parameters
	envGradParams grad = paramsGrad->getGradient();
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			// NB: gradient lies in range 0-1 for all types, and is applied when necessary...
			// ... implies gradient increment will be dimensionless in range 0-1 (but << 1)
			if (cells[y][x] != 0) { // not no-data cell
				habK = 0.0;
				int nhab = cells[y][x]->nHabitats();
				for (int i = 0; i < nhab; i++) {
					switch (rasterType) {
					case 0:
						habK += pSpecies->getHabK(cells[y][x]->getHabIndex(i));
						break;
					case 1:
						habK += pSpecies->getHabK(i) * cells[y][x]->getHabitat(i) / 100.0f;
						break;
					case 2:
						habK += pSpecies->getHabK(0) * cells[y][x]->getHabitat(i) / 100.0f;
						break;
					}
				}
				if (habK > 0.0) { // suitable cell
					if (initial) { // set local environmental deviation
						cells[y][x]->setEnvDev((float)pRandom->Random() * (2.0f) - 1.0f);
					}
					dist_from_opt = (float)(fabs((double)grad.opt_y - (double)y));
					dev = cells[y][x]->getEnvDev();
					envval = 1.0 - dist_from_opt * grad.grad_inc + dev * grad.factor;
					if (envval < 0.000001) envval = 0.0;
					if (envval > 1.0) envval = 1.0;
				}
				else envval = 0.0;
				cells[y][x]->setEnvVal((float)envval);
			}
		}
	}

}

void Landscape::setGlobalStoch(int nyears) {
	envStochParams env = paramsStoch->getStoch();
	if (epsGlobal != 0) delete[] epsGlobal;
	epsGlobal = new float[nyears];
	epsGlobal[0] = (float)(pRandom->Normal(0.0, env.std) * sqrt(1.0 - (env.ac * env.ac)));
	for (int i = 1; i < nyears; i++) {
		epsGlobal[i] = (float)(env.ac * epsGlobal[i - 1] + pRandom->Normal(0.0, env.std) * sqrt(1.0 - (env.ac * env.ac)));
	}
}

float Landscape::getGlobalStoch(int yr) {
	if (epsGlobal != 0 && yr >= 0) {
		return epsGlobal[yr];
	}
	else return 0.0;
}

void Landscape::updateLocalStoch(void) {
	envStochParams env = paramsStoch->getStoch();
	float randpart;
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] != 0) { // not a no-data cell
				randpart = (float)(pRandom->Normal(0.0, env.std) * sqrt(1.0 - (env.ac * env.ac)));
				cells[y][x]->updateEps((float)env.ac, randpart);
			}
		}
	}

}

void Landscape::resetCosts(void) {
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] != 0) { // not a no-data cell
				cells[y][x]->resetCost();
			}
		}
	}
}

void Landscape::resetEffCosts(void) {
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] != 0) { // not a no-data cell
				cells[y][x]->resetEffCosts();
			}
		}
	}
}

//---------------------------------------------------------------------------

// Dynamic landscape functions

void Landscape::setDynamicLand(bool dyn) { dynamic = dyn; }

void Landscape::addLandChange(landChange c) {
	landchanges.push_back(c);
}

int Landscape::numLandChanges(void) { return (int)landchanges.size(); }

landChange Landscape::getLandChange(short ix) {
	landChange c; c.chgNb = c.chgYear = 0;
	c.pathHabFile = c.pathPatchFile = c.pathCostFile = "none";
	int nchanges = (int)landchanges.size();
	if (ix < nchanges) c = landchanges[ix];
	return c;
}

void Landscape::deleteLandChanges(void) {
	while (landchanges.size() > 0) landchanges.pop_back();
	landchanges.clear();
}
#if RS_RCPP
int Landscape::readLandChange(int filenum, Rcpp::NumericMatrix habfile, Rcpp::NumericMatrix pchfile, Rcpp::NumericMatrix costfile, Rcpp::NumericVector scalinglayers){

	if (filenum < 0) return 19;

	int h = 0, p = 0, c = 0, pchseq = 0;
	double hfloat = 0,pfloat = 0,cfloat = 0;
	bool costs = false;
	if(costfile.nrow()>0 && costfile.ncol()>0) costs = true;

	arma::vec cellDemoScalings;  // vector to store local demog scalings
	Rcpp::IntegerVector DSdim;
	int nrDemogScaleLayers = 0;
	if(scalinglayers.attr("dim")==R_NilValue) DSdim = Rcpp::IntegerVector::create(1,1,1);
	else{
	    DSdim = scalinglayers.attr("dim");
	    if(DSdim.size()>2) nrDemogScaleLayers = DSdim[2]; //nr of slices on cube
	    else nrDemogScaleLayers = 1;
	}
	arma::cube scalingCube(scalinglayers.begin(),DSdim[0],DSdim[1],nrDemogScaleLayers,false); // turn scaling layers into a cube

	simParams sim = paramsSim->getSim();

	if (patchModel) pchseq = patchCount();

	switch (rasterType) {

	case 0: // raster with habitat codes - 100% habitat each cell

		for (int y = dimY-1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {

				// get numerics from each raster for this cell
				hfloat = habfile(dimY-1-y,x);
				if (patchModel) pfloat = pchfile(dimY-1-y,x);
				if (costs) cfloat = costfile(dimY-1-y,x);

				if (cells[y][x] != 0) { // not a no data cell (in initial landscape)
					if ( R_IsNA(hfloat) ){ // invalid no data cell in change map
						return 36;
					}
					else {
						h = (int)hfloat;
						if (h < 0 || (sim.batchMode && (h < 1 || h > nHabMax))) { // invalid habitat code
							return 33;
						}
						else {
							addHabCode(h);
							cells[y][x]->setHabIndex(h);
						}
					}
					if (patchModel) {
						if ( R_IsNA(pfloat) ){
							#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Found patch NA in valid habitat cell." <<  std::endl;
							#endif
							return 34;
						}
						else {
							p = (int)pfloat;
							if (p < 0 ) { // invalid patch code
								#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Found negative patch ID in valid habitat cell." <<  std::endl;
								#endif
								return 34;
							}
							else {
								patchChgMatrix[y][x][2] = p;
								if (p > 0 && !existsPatch(p)) {
									// addPatchNum(p); // was removed - why?
									newPatch(pchseq++,p);
								}
							}
						}
					}
					if (costs) {
						if ( R_IsNA(cfloat) ){ // invalid cost
							return 38;
						}
						else{
							c = (int)cfloat;
							if (c < 1) { // invalid cost
								return 38;
							}
							else {
								costsChgMatrix[y][x][2] = c;
							}
						}
					}
				}
			}
		}
		break;

		case 2: // habitat quality

		for (int y = dimY-1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {

				// get numerics from each raster for this cell
				hfloat = habfile(dimY-1-y,x);
				if (patchModel) pfloat = pchfile(dimY-1-y,x);
				if (costs) cfloat = costfile(dimY-1-y,x);

				if (cells[y][x] != 0) { // not a no data cell (in initial landscape)
					if ( R_IsNA(hfloat) ){ // invalid no data cell in change map
					    Rcpp::Rcout << "Found NA in valid habitat cell. For landscape nb "<< filenum + 2 <<  std::endl;
						return 36;
					}
					else {
						if (hfloat < 0.0 || hfloat > 100.0) { // invalid quality score
						    Rcpp::Rcout << "Found invalid habitat quality value " << hfloat << " in valid habitat cell." <<  std::endl;
							return 37;
						}
						else {
							cells[y][x]->setHabitat((float)hfloat);
						}
					}
					if (patchModel) {
						if ( R_IsNA(pfloat) ){
							#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Found patch NA in valid habitat cell." <<  std::endl;
							#endif
							return 34;
						}
						else {
							p = (int)pfloat;
							if (p < 0 ) { // invalid patch code
								#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Found negative patch ID in valid habitat cell." <<  std::endl;
								#endif
								return 34;
							}
							else {
								patchChgMatrix[y][x][2] = p;
								if (p > 0 && !existsPatch(p)) {
									//addPatchNum(p); // was removed - why?
									newPatch(pchseq++,p);
								}
							}
						}
					}
					if (costs) {
						if ( R_IsNA(cfloat) ){ // invalid cost
							return 38;
						}
						else{
							c = (int)cfloat;
							if (c < 1) { // invalid cost
								return 38;
							}
							else {
								costsChgMatrix[y][x][2] = c;
							}
						}
					}
					//SPATIALDEMOG
					// read demographic scalings
					if(nrDemogScaleLayers>0){
						// get tube at (y/x)
						cellDemoScalings = scalingCube(arma::span(dimY-1-y), arma::span(x), arma::span::all);
						if(cellDemoScalings.n_elem==(unsigned)nDSlayer){
							// set vector percentage values in cell
							cells[y][x]->addchgDemoScaling(arma::conv_to< std::vector<float> >::from(cellDemoScalings));
						}
						else{// invalid patch code
							#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Wrong number of demographic scaling layers in cell " << x << " ," << y << " at dyn land change nr " << filenum << std::endl;
							#endif
							return 39;
						}
					}// SPATIALDEMOG End

				}
			}
		}
		break;

	default:
		break;
	}

	return 0;
}
#endif

#if RS_RCPP && !R_CMD // normal file input
int Landscape::readLandChange(int filenum, bool costs, wifstream& hfile, wifstream& pfile, wifstream& cfile, int habnodata, int pchnodata, int costnodata, vector <string>  scalinglayers)
#else
int Landscape::readLandChange(int filenum, bool costs, vector <string>  scalinglayers)
#endif
{

#if RS_RCPP
	wstring header;
#else
	string header;
	int ncols, nrows, habnodata, costnodata, pchnodata;
	costnodata = 0;
	pchnodata = 0;
#endif
	int h = 0, p = 0, c = 0, pchseq = 0;
	float hfloat, pfloat, cfloat;
	simParams sim = paramsSim->getSim();

	if (filenum < 0) return 19;
	if (patchModel) pchseq = patchCount();

#if !RS_RCPP
	ifstream hfile; // habitat file input stream
	ifstream pfile; // patch file input stream
	ifstream cfile; // costs file input stream
#endif

#if !RS_RCPP || R_CMD
	// open habitat file and optionally also patch and costs files
	hfile.open(landchanges[filenum].pathHabFile.c_str());
	if (!hfile.is_open()) return 30;
	if (patchModel) {
		pfile.open(landchanges[filenum].pathPatchFile.c_str());
		if (!pfile.is_open()) {
			hfile.close(); hfile.clear();
			return 31;
		}
	}
	if (costs) {
		cfile.open(landchanges[filenum].pathCostFile.c_str());
		if (!cfile.is_open()) {
			hfile.close(); hfile.clear();
			if (pfile.is_open()) {
				pfile.close(); pfile.clear();
			}
			return 32;
		}
	}

	// read header records of habitat (and patch) file(s)
	// NB headers of all files have already been compared
	hfile >> header >> ncols >> header >> nrows >> header >> hfloat >> header >> hfloat
		>> header >> hfloat >> header >> habnodata;
	if (patchModel) {
		for (int i = 0; i < 5; i++) pfile >> header >> pfloat;
		pfile >> header >> pchnodata;
	}
	if (costs) {
		for (int i = 0; i < 5; i++) cfile >> header >> cfloat;
		cfile >> header >> costnodata;
	}
#endif

	// set up bad float values to ensure that valid values are read
	float badhfloat = -9.0; if (habnodata == -9) badhfloat = -99.0;
	float badpfloat = -9.0; if (pchnodata == -9) badpfloat = -99.0;
	float badcfloat = -9.0; if (costnodata == -9) badcfloat = -99.0;

	switch (rasterType) {

	case 0: // raster with habitat codes - 100% habitat each cell
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				hfloat = badhfloat;
#if RS_RCPP
				if (hfile >> hfloat) {
#else
				hfile >> hfloat;
#endif
				h = (int)hfloat;
#if RS_RCPP
				}
	else {
		// corrupt file stream
#if RS_RCPP && !R_CMD
		Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
		StreamErrorR("habitatchgfile");
		hfile.close();
		hfile.clear();
		pfile.close();
		pfile.clear();
		return 171;
	}
#endif
	if (patchModel) {
		pfloat = badpfloat;
#if RS_RCPP
		if (pfile >> pfloat) {
#else
		pfile >> pfloat;
#endif
		p = (int)pfloat;
#if RS_RCPP
		}
	else {
		// corrupt file stream
#if RS_RCPP && !R_CMD
		Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
		StreamErrorR("patchchgfile");
		hfile.close();
		hfile.clear();
		pfile.close();
		pfile.clear();
		return 172;
	}
#endif
	}
		if (costs) {
			cfloat = badcfloat;
#if RS_RCPP
			if (cfile >> cfloat) {
#else
			cfile >> cfloat;
#endif
			c = (int)cfloat;
#if RS_RCPP
			}
		else {
			// corrupt file stream
#if RS_RCPP && !R_CMD
			Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
			StreamErrorR("costchgfile");
			hfile.close();
			hfile.clear();
			pfile.close();
			pfile.clear();
			return 173;
		}
#endif
		}
	if (cells[y][x] != 0) { // not a no data cell (in initial landscape)
		if (h == habnodata) { // invalid no data cell in change map
			hfile.close(); hfile.clear();
			return 36;
		}
		else {
			if (h < 0 || (sim.batchMode && (h < 1 || h > nHabMax))) {
				// invalid habitat code
				hfile.close(); hfile.clear();
				if (patchModel) { pfile.close(); pfile.clear(); }
				return 33;
			}
			else {
				addHabCode(h);
				cells[y][x]->setHabIndex(h);
			}
		}
		if (patchModel) {
			if (p < 0 || p == pchnodata) { // invalid patch code
#if RS_RCPP && !R_CMD
				if (p == pchnodata) Rcpp::Rcout << "Found patch NA in valid habitat cell." << std::endl;
				else Rcpp::Rcout << "Found negative patch ID in valid habitat cell." << std::endl;
#endif
				hfile.close(); hfile.clear();
				pfile.close(); pfile.clear();
				return 34;
			}
			else {
				patchChgMatrix[y][x][2] = p;
				if (p > 0 && !existsPatch(p)) {
					// addPatchNum(p); why is it removed?
					newPatch(pchseq++, p);
				}
			}
		}
		if (costs) {
			if (c < 1) { // invalid cost
#if RS_RCPP
			    Rcpp::Rcout << "Found invalid cost value of " << c << " in cell x " << x << " and y  " << y << std::endl;
#endif
				hfile.close(); hfile.clear();
				if (pfile.is_open()) {
					pfile.close(); pfile.clear();
				}
				return 38;
			}
			else {
				costsChgMatrix[y][x][2] = c;
			}
		}
	}
			} // for x
		} // for y
#if RS_RCPP
		hfile >> hfloat;
		if (!hfile.eof()) EOFerrorR("habitatchgfile");
		if (patchModel)
		{
			pfile >> pfloat;
			if (!pfile.eof()) EOFerrorR("patchchgfile");
		}
		if (costs)
		{
			cfile >> cfloat;
			if (!cfile.eof()) EOFerrorR("costchgfile");
		}
#endif
		break;

	case 2: // habitat quality
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				hfloat = badhfloat;
#if RS_RCPP
				if (hfile >> hfloat) {
#else
				hfile >> hfloat;
#endif
				h = (int)hfloat;
#if RS_RCPP
				}
		else {
			// corrupt file stream
#if RS_RCPP && !R_CMD
			Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
			StreamErrorR("habitatchgfile");
			hfile.close();
			hfile.clear();
			pfile.close();
			pfile.clear();
			return 172;
		}
#endif
		if (patchModel) {
			pfloat = badpfloat;
#if RS_RCPP
			if (pfile >> pfloat) {
#else
			pfile >> pfloat;
#endif
			p = (int)pfloat;
#if RS_RCPP
			}
		else {
			// corrupt file stream
#if RS_RCPP && !R_CMD
			Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
			StreamErrorR("patchchgfile");
			hfile.close();
			hfile.clear();
			pfile.close();
			pfile.clear();
			return 175;
		}
#endif
		}
		if (costs) {
			cfloat = badcfloat;
#if RS_RCPP
			if (cfile >> cfloat) {
#else
			cfile >> cfloat;
#endif
			c = (int)cfloat;
#if RS_RCPP
			}
		else {
			// corrupt file stream
#if RS_RCPP && !R_CMD
			Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
			StreamErrorR("costchgfile");
			hfile.close();
			hfile.clear();
			pfile.close();
			pfile.clear();
			return 173;
		}
#endif
		}
		if (cells[y][x] != 0) { // not a no data cell (in initial landscape)
			if (h == habnodata) { // invalid no data cell in change map
				hfile.close(); hfile.clear();
				if (patchModel) { pfile.close(); pfile.clear(); }
				return 36;
			}
			else {
				if (hfloat < 0.0 || hfloat > 100.0) { // invalid quality score
					hfile.close(); hfile.clear();
					if (patchModel) { pfile.close(); pfile.clear(); }
					return 37;
				}
				else {
					cells[y][x]->setHabitat(hfloat);
				}
			}
			if (patchModel) {
				if (p < 0 || p == pchnodata) { // invalid patch code
#if RS_RCPP && !R_CMD
					if (p == pchnodata) Rcpp::Rcout << "Found patch NA in valid habitat cell." << std::endl;
					else Rcpp::Rcout << "Found negative patch ID in valid habitat cell." << std::endl;
#endif
					hfile.close(); hfile.clear();
					pfile.close(); pfile.clear();
					return 34;
				}
				else {
					patchChgMatrix[y][x][2] = p;
					if (p > 0 && !existsPatch(p)) {
						// addPatchNum(p); // was removed - why?
						newPatch(pchseq++, p);
					}
				}
			}
			if (costs) {
				if (c < 1) { // invalid cost
#if RS_RCPP
				    Rcpp::Rcout << "Found invalid cost value of " << c << "in cell x " << x << " and y  " << y << std::endl;
#endif
					hfile.close(); hfile.clear();
					if (pfile.is_open()) {
						pfile.close(); pfile.clear();
					}
					return 38;
				}
				else {
					costsChgMatrix[y][x][2] = c;
				}
			}
		}
			} // end x
		} // end y
#if RS_RCPP
		hfile >> hfloat;
		if (!hfile.eof()) EOFerrorR("habitatchgfile");
		if (patchModel)
		{
			pfile >> pfloat;
			if (!pfile.eof()) EOFerrorR("patchchgfile");
		}
		if (costs)
		{
			cfile >> cfloat;
			if (!cfile.eof()) EOFerrorR("costchgfile");
		}
#endif
		break;

default:
	break;
	}

	if (hfile.is_open()) { hfile.close(); hfile.clear(); }
	if (pfile.is_open()) { pfile.close(); pfile.clear(); }
	if (cfile.is_open()) { cfile.close(); cfile.clear(); }

	// add here the reading of demographic scaling layers
	if(scalinglayers.size()>0){
		int retcode = readDemographicScaling(scalinglayers);
		if (retcode < 0) return 54; //change number
	}

	return 0;

}

// Create & initialise patch change matrix
void Landscape::createPatchChgMatrix(void)
{
	Patch* pPatch;
	Cell* pCell;
	if (patchChgMatrix != 0) deletePatchChgMatrix();
	patchChgMatrix = new int** [dimY];
	for (int y = dimY - 1; y >= 0; y--) {
		patchChgMatrix[y] = new int* [dimX];
		for (int x = 0; x < dimX; x++) {
			patchChgMatrix[y][x] = new int[3];
			pCell = findCell(x, y);
			if (pCell == 0) { // no-data cell
				patchChgMatrix[y][x][0] = patchChgMatrix[y][x][1] = 0;
			}
			else {
				// record initial patch number
				pPatch = pCell->getPatch();
				if (pPatch == nullptr) { // matrix cell
					patchChgMatrix[y][x][0] = patchChgMatrix[y][x][1] = 0;
				}
				else {
					patchChgMatrix[y][x][0] = patchChgMatrix[y][x][1] = pPatch->getPatchNum();
				}
			}
			patchChgMatrix[y][x][2] = 0;
		}
	}
}

void Landscape::recordPatchChanges(int landIx) {
	if (patchChgMatrix == 0) return; // should not occur
	patchChange chg;

	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (landIx == 0) { // reset to original landscape
				if (patchChgMatrix[y][x][0] != patchChgMatrix[y][x][2]) {
					// record change of patch for current cell
					chg.chgnum = 666666; chg.x = x; chg.y = y;
					chg.oldpatch = patchChgMatrix[y][x][2];
					chg.newpatch = patchChgMatrix[y][x][0];
					patchchanges.push_back(chg);
				}
			}
			else { // any other change
				if (patchChgMatrix[y][x][2] != patchChgMatrix[y][x][1]) {
					// record change of patch for current cell
					chg.chgnum = landIx; chg.x = x; chg.y = y;
					chg.oldpatch = patchChgMatrix[y][x][1];
					chg.newpatch = patchChgMatrix[y][x][2];
					patchchanges.push_back(chg);
				}
			}
			// reset cell for next landscape change
			patchChgMatrix[y][x][1] = patchChgMatrix[y][x][2];
		}
	}

}

int Landscape::numPatchChanges(void) { return (int)patchchanges.size(); }

patchChange Landscape::getPatchChange(int i) {
	patchChange c;
	c.chgnum = 99999999;
	c.x = c.y = c.oldpatch = c.newpatch = -1;
	if (i >= 0 && i < (int)patchchanges.size())
		c = patchchanges[i];
	return c;
}

void Landscape::deletePatchChgMatrix(void) {
	if (patchChgMatrix != 0) {
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				delete[] patchChgMatrix[y][x];
			}
			delete[] patchChgMatrix[y];
		}
	}
	patchChgMatrix = 0;
}

// Create & initialise costs change matrix
void Landscape::createCostsChgMatrix(void)
{
	//intptr patch;
	//Patch *pPatch;
	Cell* pCell;
	if (costsChgMatrix != 0) deleteCostsChgMatrix();
	costsChgMatrix = new int** [dimY];
	for (int y = dimY - 1; y >= 0; y--) {
		costsChgMatrix[y] = new int* [dimX];
		for (int x = 0; x < dimX; x++) {
			costsChgMatrix[y][x] = new int[3];
			pCell = findCell(x, y);
			if (pCell == 0) { // no-data cell
				costsChgMatrix[y][x][0] = costsChgMatrix[y][x][1] = 0;
			}
			else {
				// record initial cost
				costsChgMatrix[y][x][0] = costsChgMatrix[y][x][1] = pCell->getCost();
			}
			costsChgMatrix[y][x][2] = 0;
		}
	}
}

void Landscape::recordCostChanges(int landIx) {

	if (costsChgMatrix == 0) return; // should not occur
	costChange chg;

	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (landIx == 0) { // reset to original landscape
				if (costsChgMatrix[y][x][0] != costsChgMatrix[y][x][2]) {
					// record change of cost for current cell
					chg.chgnum = 666666; chg.x = x; chg.y = y;
					chg.oldcost = costsChgMatrix[y][x][2];
					chg.newcost = costsChgMatrix[y][x][0];
					costschanges.push_back(chg);
				}
			}
			else { // any other change
				if (costsChgMatrix[y][x][2] != costsChgMatrix[y][x][1]) {
					// record change of cost for current cell
					chg.chgnum = landIx; chg.x = x; chg.y = y;
					chg.oldcost = costsChgMatrix[y][x][1];
					chg.newcost = costsChgMatrix[y][x][2];
					costschanges.push_back(chg);
				}
			}
			// reset cell for next landscape change
			costsChgMatrix[y][x][1] = costsChgMatrix[y][x][2];
		}
	}

}

int Landscape::numCostChanges(void) { return (int)costschanges.size(); }

costChange Landscape::getCostChange(int i) {
	costChange c; c.chgnum = 99999999; c.x = c.y = c.oldcost = c.newcost = -1;
	if (i >= 0 && i < (int)costschanges.size()) c = costschanges[i];
	return c;
}

void Landscape::deleteCostsChgMatrix(void) {
	if (costsChgMatrix != 0) {
		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
				delete[] costsChgMatrix[y][x];
			}
			delete[] costsChgMatrix[y];
		}
	}
	costsChgMatrix = 0;
}

//---------------------------------------------------------------------------

// Species distribution functions

//If file input as R objects only for #RS_RCPP
#if RS_RCPP
int Landscape::newDistribution(Species *pSpecies, Rcpp::NumericMatrix distname, int spResol) {
	int readcode;
	int ndistns = (int)distns.size();
	distns.push_back(new InitDist(pSpecies));
	landOrigin habfile_origin = this->getOrigin();
	readcode = distns[ndistns]->readDistribution(distname,habfile_origin,spResol);
	if (readcode != 0) { // error encountered
		// delete the distribution created above
		delete distns[ndistns];
		distns.pop_back();
	}
	return readcode;
}
#endif
// Standard file input
int Landscape::newDistribution(Species *pSpecies, string distname) {
	int readcode;
	int ndistns = (int)distns.size();
	distns.push_back(new InitDist(pSpecies));
	readcode = distns[ndistns]->readDistribution(distname);
	if (readcode != 0) { // error encountered
		// delete the distribution created above
		delete distns[ndistns];
		distns.pop_back();
	}
	return readcode;
}


void Landscape::setDistribution(Species* pSpecies, int nInit) {
	// WILL NEED TO SELECT DISTRIBUTION FOR CORRECT SPECIES ...
	// ... CURRENTLY IT IS THE ONLY ONE
	distns[0]->setDistribution(nInit);
}

// Specified cell match one of the distribution cells to be initialised?
bool Landscape::inInitialDist(Species* pSpecies, locn loc) {
	// convert landscape co-ordinates to distribution co-ordinates
	locn initloc;
	initloc.x = loc.x * resol / spResol;
	initloc.y = loc.y * resol / spResol;
	// WILL HAVE TO GET CORRECT SPECIES WHEN THERE ARE MULTIPLE SPECIES ...
	bool initialise = distns[0]->inInitialDist(initloc);
	return initialise;
}

void Landscape::deleteDistribution(Species* pSpecies) {
	if (distns[0] != 0) delete distns[0]; distns.clear();
}

// Return no. of initial distributions
int Landscape::distnCount(void) {
	return (int)distns.size();
}

int Landscape::distCellCount(int dist) {
	return distns[dist]->cellCount();
}

// Set a cell in a specified initial distribution (by position in cells vector)
void Landscape::setDistnCell(int dist, int ix, bool init) {
	distns[dist]->setDistCell(ix, init);
}

// Set a cell in a specified initial distribution (by given co-ordinates)
void Landscape::setDistnCell(int dist, locn loc, bool init) {
	distns[dist]->setDistCell(loc, init);
}

// Get the co-ordinates of a specified cell in a specified initial distribution
locn Landscape::getDistnCell(int dist, int ix) {
	return distns[dist]->getCell(ix);
}

// Get the co-ordinates of a specified cell in a specified initial distribution
// Returns negative co-ordinates if the cell is not selected
locn Landscape::getSelectedDistnCell(int dist, int ix) {
	return distns[dist]->getSelectedCell(ix);
}

// Get the dimensions of a specified initial distribution
locn Landscape::getDistnDimensions(int dist) {
	return distns[dist]->getDimensions();
}

// Reset the distribution for a given species so that all cells are deselected
void Landscape::resetDistribution(Species* pSp) {
	// CURRENTLY WORKS FOR FIRST SPECIES ONLY ...
	distns[0]->resetDistribution();
}

//---------------------------------------------------------------------------

// Initialisation cell functions

int Landscape::initCellCount(void) {
	return (int)initcells.size();
}

void Landscape::addInitCell(int x, int y) {
	initcells.push_back(new DistCell(x, y));
}

locn Landscape::getInitCell(int ix) {
	return initcells[ix]->getLocn();
}

void Landscape::clearInitCells(void) {
	int ncells = (int)initcells.size();
	for (int i = 0; i < ncells; i++) {
		delete initcells[i];
	}
	initcells.clear();
}

//---------------------------------------------------------------------------

// Read landscape file(s)
// Returns error code or zero if read correctly

// for new landscape input using R objects AND spatial demography: only for RS_RCPP
// RS_THREADSAFE and SPATIALDEMOG
#if RS_RCPP
int Landscape::readLandscape(int fileNum, Rcpp::NumericMatrix habfile, Rcpp::NumericMatrix pchfile, Rcpp::NumericMatrix costfile, Rcpp::NumericVector scalinglayers) {
	if (fileNum < 0) return 19;

    int habCode,seq,patchCode,ncols,nrows,hc,maxcost = 0;
    double habFloat,patchFloat,costFloat;
	Patch *pPatch;
	Cell *pCell;
	simParams sim = paramsSim->getSim();
	initParams init = paramsInit->getInit();

	// initialise landscape size
	ncols = habfile.ncol();
	nrows = habfile.nrow();
	dimX = ncols; dimY = nrows;
	minX = maxY = 0;
	maxX = dimX-1; maxY = dimY-1;
	if (fileNum == 0) {
		// set initialisation limits to landscape limits
		init.minSeedX = init.minSeedY = 0;
		init.maxSeedX = maxX; init.maxSeedY = maxY;
		paramsInit->setInit(init);
		setCellArray();
	}

	seq = 0; 	// initial sequential patch landscape
    patchCode= 0; 		// initial patch number for cell-based landscape
	// create patch 0 - the matrix patch (even if there is no matrix)
    if (fileNum == 0) {
        newPatch(seq,patchCode);
        seq++;
        patchCode++;
    }

	switch (rasterType) {

	case 0: // raster with habitat codes - 100% habitat each cell
		if (fileNum > 0) return 19; // error condition - should not occur

		for (int y = dimY-1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {

				// read value from raster cell
                habFloat = habfile(dimY-1-y,x);
				// check for NA
                if ( R_IsNA(habFloat) )
					addNewCellToLand(x,y,-1); // add cell only to landscape
				else {
                    habCode= static_cast<int>(habFloat);
                    if (habCode< 0 || (sim.batchMode && (habCode< 1 || habCode> nHabMax))) {
						// invalid habitat code
						#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Found invalid habitat code." <<  std::endl;
						#endif
						return 13;
					}
                    else { // valid habitat code
                        addHabCode(habCode);
						if (patchModel) {
                            patchFloat = pchfile(dimY-1-y,x);
                            if ( R_IsNA(patchFloat) ) { // invalid patch code
								#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Found patch NA in valid habitat cell." <<  std::endl;
								#endif
								return 14;
							}
                            patchCode= static_cast<int>(patchFloat);
                            if (patchCode< 0 ) { // invalid patch code
								#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Found negative patch ID in valid habitat cell." <<  std::endl;
								#endif
								return 14;
							}

                            // Does the patch already exists?
                            pPatch = patchCode == 0 ? nullptr : ( // matrix cell
                                existsPatch(patchCode) ?
                            findPatch(patchCode) :
                                newPatch(seq++, patchCode)
                            );
                            addNewCellToPatch(pPatch, x, y, habCode);
						}
						else { // cell-based model
							// add cell to landscape (patches created later)
                            addNewCellToLand(x,y,habCode);
						}
					}
                } // end of habFloat is NA
            } // end x
        } // end y
		break;

	case 1: // multiple % cover

		for (int y = dimY-1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {

                habFloat = habfile(dimY-1-y,x);
				if (fileNum == 0) { // first habitat cover layer
                    if ( R_IsNA(habFloat) ) { // check for NA
						addNewCellToLand(x,y,-1); // add cell only to landscape
					}
					else {
                        if (habFloat < 0.0 || habFloat > 100.0) { // invalid cover score
							#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Found invalid habitat cover score." <<  std::endl;
							#endif
							return 17;
						}
						else {
							if (patchModel) {
                                patchFloat = pchfile(dimY-1-y,x);
                                if ( R_IsNA(patchFloat) ) { // invalid patch code
									#if RS_RCPP && !R_CMD
									Rcpp::Rcout << "Found patch NA in valid habitat cell." <<  std::endl;
									#endif
									return 14;
								}

                                patchCode= static_cast<int>(patchFloat);

                                if (patchCode< 0 ) { // invalid patch code
									#if RS_RCPP && !R_CMD
									Rcpp::Rcout << "Found negative patch ID in valid habitat cell." <<  std::endl;
									#endif
									return 14;
								}

                                pPatch = patchCode == 0 ? nullptr : ( // matrix cell
                                    existsPatch(patchCode) ?
                                findPatch(patchCode) :
                                    newPatch(seq++, patchCode)
                                );
                                addNewCellToPatch(pPatch, x, y, (float)habFloat);
							}
							else { // cell-based model
								// add cell to landscape (patches created later)
                                addNewCellToLand(x,y,(float)habFloat);
							}
						}
					}
				}
				else { // additional habitat cover layers
                    if ( !R_IsNA(habFloat) ) {
                        if (habFloat < 0.0 || habFloat > 100.0) { // invalid cover score
							#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Found invalid habitat cover score." <<  std::endl;
							#endif
							return 17;
						}
						else {
                            cells[dimY-1-y][x]->setHabitat((float)habFloat);
						}
                    } // end of habFloat is not NA
                } // end additional habitat cover layers
            } // end x
        } // end y
		habIndexed = true; // habitats are already numbered 1...n in correct order

		break;

	case 2: // habitat quality
		if (fileNum > 0) return 19; // error condition - should not occur

		for (int y = dimY-1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {
                habFloat = habfile(dimY-1-y,x);
                if ( R_IsNA(habFloat) ) { // check for NA
					addNewCellToLand(x,y,-1); // add cell only to landscape
				}
				else {
                    if (habFloat < 0.0 || habFloat > 100.0) { // invalid quality score
						#if RS_RCPP && !R_CMD
						Rcpp::Rcout << "Found invalid habitat quality score." <<  std::endl;
						#endif
						return 17;
					}
                    else { // valid quality score
						if (patchModel) {
                            patchFloat = pchfile(dimY-1-y,x);
                            if ( R_IsNA(patchFloat) ) { // invalid patch code
								#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Found patch NA in valid habitat cell." <<  std::endl;
								#endif
								return 14;
							}
                            patchCode= static_cast<int>(patchFloat);
                            if (patchCode< 0 ) { // invalid patch code
								#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Found negative patch ID in valid habitat cell." <<  std::endl;
								#endif
								return 14;
							}

                            pPatch = patchCode == 0 ? nullptr : ( // matrix cell
                                existsPatch(patchCode) ?
                            findPatch(patchCode) :
                                newPatch(seq++, patchCode)
                            );
                            addNewCellToPatch(pPatch, x, y, (float)habFloat);
						}
						else { // cell-based model
							// add cell to landscape (patches created later)
                            addNewCellToLand(x,y,(float)habFloat);
						}
                    } // end of valid quality score
                } // end of habFloat is not NA
            } // end x
        } // end y
		break;

	default:
		break;
	} // end switch(rasterType)

//#if SPATIALDEMOG
	int SMScosts = costfile.nrow()*costfile.ncol();

	arma::vec cellDemoScalings;  // vector to store local demog scalings
	Rcpp::IntegerVector DSdim;
	int nrDemogScaleLayers = 0;
	if(scalinglayers.attr("dim")==R_NilValue) DSdim = Rcpp::IntegerVector::create(1,1,1);
	else{
		DSdim = scalinglayers.attr("dim");
		if(DSdim.size()>2) nrDemogScaleLayers = DSdim[2]; //nr of slices on cube
		else nrDemogScaleLayers = 1;
	}
	arma::cube scalingCube(scalinglayers.begin(),DSdim[0],DSdim[1],nrDemogScaleLayers,false); // turn scaling layers into a cube, what happens if it is not a spatial demog?

	if(SMScosts || nrDemogScaleLayers){ // are there SMS costs or demographic scaling layers to read?

		for (int y = dimY-1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {

				// find the cell
				pCell = findCell(x,y);
				if (pCell != 0) { // not no-data cell

					// read cost raster
					if(SMScosts) {
                        costFloat = costfile(dimY-1-y,x);
                        if ( !R_IsNA(costFloat) ) {
                            hc = (int)costFloat;
							if ( hc < 1 ) {
							#if RS_RCPP && !R_CMD
								Rcpp::Rcout << "Cost map may only contain values of 1 or higher, but found " << hc << "." << endl;
							#endif
								return 54;
							}
							// set cost value
							pCell->setCost(hc);
							if (hc > maxcost) maxcost = hc;
						}
					}

					// read demographic scalings
					if(nrDemogScaleLayers){
						// get tube at (y/x)
						cellDemoScalings = scalingCube(arma::span(dimY-1-y), arma::span(x), arma::span::all);
						if(cellDemoScalings.n_elem==(unsigned)nDSlayer){
							// set vector percentage values in cell
							pCell->addchgDemoScaling(arma::conv_to< std::vector<float> >::from(cellDemoScalings));
						}
						else{// invalid patch code
							#if RS_RCPP && !R_CMD
							Rcpp::Rcout << "Wrong number of demographic scaling layers in cell " << x << " ," << y << " of first layer array." << std::endl;
							#endif
							return 64;
						}
					}
				}
			}
		}
	}
	return 0;
}
#endif

int Landscape::readLandscape(int fileNum, string habfile, string pchfile, string costfile, vector <string> scalinglayers) // vector <string> scalinglayers is the vector of filenames for spatial demography of year 1
{
	// fileNum == 0 for (first) habitat file and optional patch file
	// fileNum > 0  for subsequent habitat files under the %cover option

#if RS_RCPP
	wstring header;
#else
	string header;
#endif
	int habCode, seq, patchCode, noDataHabCode;
	int noDataPatch = 0;
	int ncols, nrows;
	float habFloat, patchFloat;
	Patch* pPatch;
	simParams sim = paramsSim->getSim();

	if (fileNum < 0) return 19;

#if RS_RCPP
	wifstream ifsHabMap; // habitat file input stream
	wifstream ifsPatchMap; // patch file input stream
#else
	ifstream ifsHabMap; // habitat file input stream
	ifstream ifsPatchMap; // patch file input stream
#endif
	initParams init = paramsInit->getInit();

	// Open habitat file and optionally also patch file
// #if RS_RCPP
// 	hfile.open(habfile, std::ios::binary);
// 	if (landraster.utf) {
// 		// apply BOM-sensitive UTF-16 facet
// 		hfile.imbue(std::locale(hfile.getloc(), new std::codecvt_utf16<wchar_t, 0x10ffff, std::consume_header>));
// 	}
// #else
	ifsHabMap.open(habfile.c_str());
// #endif
// #endif

	if (!ifsHabMap.is_open()) return 11;

	if (fileNum == 0) {
		if (patchModel) {
// #if RS_RCPP
// 			pfile.open(pchfile, std::ios::binary);
// 			if (patchraster.utf) {
// 				// apply BOM-sensitive UTF-16 facet
// 				pfile.imbue(std::locale(pfile.getloc(), new std::codecvt_utf16<wchar_t, 0x10ffff, std::consume_header>));
// 			}
// #else
			ifsPatchMap.open(pchfile.c_str());
// #endif
// #endif
			if (!ifsPatchMap.is_open()) {
				ifsHabMap.close();
				ifsHabMap.clear();
				return 12;
			}
		}
	}

// read landscape data from header records of habitat file
// NB headers of all files have already been compared
// !! DO NOT REMOVE THIS WORKAROUND !! NEEDED FOR LANDSCAPES GENERATED BY TERRA (R PACKAGE)
double tmpresol;
	ifsHabMap >> header >> ncols
		>> header >> nrows
		>> header >> minEast
		>> header >> minNorth
		>> header >> tmpresol
		>> header >> noDataHabCode;
resol = (int) tmpresol;
#if RS_RCPP
	if (!ifsHabMap.good()) {
		// corrupt file stream
		StreamErrorR(habfile);
		ifsHabMap.close();
		ifsHabMap.clear();
		if (patchModel) {
			ifsPatchMap.close();
			ifsPatchMap.clear();
		}
		return 131;
	}
#endif

	dimX = ncols;
	dimY = nrows;
	minX = maxY = 0; // bottom-left / south-west corner
	maxX = dimX - 1;
	maxY = dimY - 1;

	if (fileNum == 0) { // First map layer

		// Set initialisation limits to landscape limits
		init.minSeedX = init.minSeedY = 0;
		init.maxSeedX = maxX;
		init.maxSeedY = maxY;
		paramsInit->setInit(init);

		if (patchModel) {
			for (int i = 0; i < 5; i++)
				ifsPatchMap >> header >> patchFloat;
			ifsPatchMap >> header >> noDataPatch;
		}

#if RS_RCPP
		if (!ifsPatchMap.good()) {
			// corrupt file stream
			StreamErrorR(pchfile);
			ifsHabMap.close();
			ifsHabMap.clear();
			ifsPatchMap.close();
			ifsPatchMap.clear();
			return 135;
		}
#endif
		setCellArray();
	}

	// set up bad float values to ensure that valid values are read
	float badHabFloat = -9.0;
	if (noDataHabCode == -9) badHabFloat = -99.0;
	float badPatchFloat = -9.0;
	if (noDataPatch == -9) badPatchFloat = -99.0;

	seq = 0; 	// initial sequential patch landscape
	patchCode = 0; 		// initial patch number for cell-based landscape

	// Create the matrix patch (even if there is no matrix)
	if (fileNum == 0) {
		newPatch(seq, patchCode);
		seq++;
		patchCode++;
	}

	switch (rasterType) {

	case 0: // Raster with habitat codes - 100% habitat each cell
		if (fileNum > 0) return 19; // error condition - should not occur

		for (int y = dimY - 1; y >= 0; y--) {
			for (int x = 0; x < dimX; x++) {

				habFloat = badHabFloat;
#if RS_RCPP
				if (ifsHabMap >> habFloat) {
					habCode = (int)habFloat;
					if (patchModel) {
						patchFloat = badPatchFloat;
						if (ifsPatchMap >> patchFloat) {
							patchCode = (int)patchFloat;
						}
						else { // corrupt file stream
#if !R_CMD
							Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
							StreamErrorR(pchfile);
							ifsHabMap.close();
							ifsHabMap.clear();
							ifsPatchMap.close();
							ifsPatchMap.clear();
							return 132;
						}
				}
				}
				else { // corrupt file stream
#if !R_CMD
			Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
					StreamErrorR(habfile);
					ifsHabMap.close();
					ifsHabMap.clear();
			if (patchModel) {
						ifsPatchMap.close();
						ifsPatchMap.clear();
					}
			return 135;
		}
#else
				// Read habitat code in this cell
				ifsHabMap >> habFloat;
				habCode = static_cast<int>(habFloat);
				if (patchModel) {
					// Read patch code in this cell
					patchFloat = badPatchFloat;
					ifsPatchMap >> patchFloat;
					patchCode = static_cast<int>(patchFloat);
				}
#endif
				if (habCode == noDataHabCode) {
					// No-data cell
					addNewCellToLand(x, y, -1); // x, y is a null cell
				}
				else if (habCode < 0 || (sim.batchMode && (habCode < 1 || habCode > nHabMax))) {
					// Invalid habitat code
#if RS_RCPP && !R_CMD
				Rcpp::Rcout << "Found invalid habitat code." << std::endl;
#endif
					ifsHabMap.close();
					ifsHabMap.clear();
					if (patchModel) {
						ifsPatchMap.close();
						ifsPatchMap.clear();
					}
					return 13;
				}
				else { // Valid habitat code
					addHabCode(habCode);
					if (patchModel) {
						if (patchCode < 0 || patchCode == noDataPatch) { // invalid patch code
#if RS_RCPP && !R_CMD
							if (patchCode == noDataPatch) Rcpp::Rcout << "Found patch NA in valid habitat cell." << std::endl;
							else Rcpp::Rcout << "Found negative patch ID in valid habitat cell." << std::endl;
#endif
							ifsHabMap.close();
							ifsHabMap.clear();
							ifsPatchMap.close();
							ifsPatchMap.clear();
							return 14;
						}
						// Does the patch already exists?
						pPatch = patchCode == 0 ? nullptr : ( // matrix cell
							existsPatch(patchCode) ?
							findPatch(patchCode) :
							newPatch(seq++, patchCode)
							);
						addNewCellToPatch(pPatch, x, y, habCode);
					}
					else { // cell-based model
						// add cell to landscape (patches created later)
						addNewCellToLand(x, y, habCode);
					}
				}

			} // for x
		} // for y

#if RS_RCPP
		ifsHabMap >> habFloat;
		if (!ifsHabMap.eof()) EOFerrorR(habfile);
		if (patchModel) {
			ifsPatchMap >> patchFloat;
			if (!ifsPatchMap.eof()) EOFerrorR(pchfile);
		}
#endif
break;

case 1: // multiple % cover
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {

				habFloat = badHabFloat;
#if RS_RCPP
				if (ifsHabMap >> habFloat) {
					habCode = static_cast<int>(habFloat);
					if (fileNum == 0) { // first habitat cover layer
						if (patchModel) {
							patchFloat = badPatchFloat;
							if (ifsPatchMap >> patchFloat) {
								patchCode = static_cast<int>(patchFloat);
					}
				else {
					// corrupt file stream
#if !R_CMD
					Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
								StreamErrorR(pchfile);
								ifsHabMap.close();
								ifsHabMap.clear();
								ifsPatchMap.close();
								ifsPatchMap.clear();
					return 135;
				}
                                  } // end if patchModel
#else
				ifsHabMap >> habFloat;
				habCode = static_cast<int>(habFloat);

				if (fileNum == 0) { // first habitat cover layer
					if (patchModel) {
						patchFloat = badPatchFloat;
						ifsPatchMap >> patchFloat;
						patchCode = static_cast<int>(patchFloat);
                                      } // end if patchModel
#endif
					if (habCode == noDataHabCode) {
				addNewCellToLand(x, y, -1); // add cell only to landscape
			}
					else if (habFloat < 0.0 || habFloat > 100.0) { // invalid cover score
#if RS_RCPP && !R_CMD
					Rcpp::Rcout << "Found invalid habitat cover score." << std::endl;
#endif
						ifsHabMap.close();
						ifsHabMap.clear();
						if (patchModel) {
							ifsPatchMap.close();
							ifsPatchMap.clear();
						}
						return 17;
					}
					else {
						if (patchModel) {
							if (patchCode < 0 || patchCode == noDataPatch) { // invalid patch code
#if RS_RCPP && !R_CMD
								if (patchCode == noDataPatch) Rcpp::Rcout << "Found patch NA in valid habitat cell." << std::endl;
							else Rcpp::Rcout << "Found negative patch ID in valid habitat cell." << std::endl;
#endif
								ifsHabMap.close();
								ifsHabMap.clear();
								ifsPatchMap.close();
								ifsPatchMap.clear();
								return 14;
							}

							pPatch = patchCode == 0 ? nullptr : ( // matrix cell
								existsPatch(patchCode) ?
								findPatch(patchCode) :
								newPatch(seq++, patchCode)
								);
							addNewCellToPatch(pPatch, x, y, habFloat);
                                          } // end patchModel
						else { // cell-based model
							// add cell to landscape (patches created later)
							addNewCellToLand(x, y, habFloat);
						}
                                      } // if valid habFloat
				}
else { // additional habitat cover layers
					if (habCode != noDataHabCode) {
						if (habFloat < 0.0 || habFloat > 100.0) { // invalid cover score
#if RS_RCPP && !R_CMD
			Rcpp::Rcout << "Found invalid habitat cover score." << std::endl;
#endif
							ifsHabMap.close();
							ifsHabMap.clear();
			if (patchModel) {
								ifsPatchMap.close();
								ifsPatchMap.clear();
			}
			return 17;
		}
		else {
							cells[y][x]->setHabitat(habFloat);
						}
					}
				}
#if RS_RCPP
                              }
		else { // not ifsHabMap >> habFloat
	// corrupt file stream
#if RS_RCPP && !R_CMD
	Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
			StreamErrorR(habfile);
			ifsHabMap.close();
			ifsHabMap.clear();
			if (patchModel) {
				ifsPatchMap.close();
				ifsPatchMap.clear();
	}
	return 133;
                              }// end not ifsHabMap >> habFloat
#endif

			} // for x
		} // for y

	habIndexed = true; // habitats are already numbered 1...n in correct order

#if RS_RCPP
		ifsHabMap >> habFloat;
		if (!ifsHabMap.eof()) EOFerrorR(habfile);
		if (patchModel) {
			ifsPatchMap >> patchFloat;
			if (!ifsPatchMap.eof()) EOFerrorR(pchfile);
	}
#endif
	break;

case 2: // habitat quality
	if (fileNum > 0) return 19; // error condition - should not occur
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {

				habFloat = badHabFloat;
#if RS_RCPP
				if (ifsHabMap >> habFloat) {
					habCode = static_cast<int>(habFloat);

			}
	else {
		// corrupt file stream
#if !R_CMD
		Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
					StreamErrorR(habfile);
					ifsHabMap.close();
					ifsHabMap.clear();
					if (patchModel) {
						ifsPatchMap.close();
						ifsPatchMap.clear();
		}
		return 134;
	}
	if (patchModel) {
					patchFloat = badPatchFloat;
					if (ifsPatchMap >> patchFloat) {
						patchCode = static_cast<int>(patchFloat);
		}
	else {
		// corrupt file stream
#if RS_RCPP && !R_CMD
		Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
						StreamErrorR(pchfile);
						ifsHabMap.close();
						ifsHabMap.clear();
						ifsPatchMap.close();
						ifsPatchMap.clear();
		return 135;
	}
				}
#else
				ifsHabMap >> habFloat;
				habCode = static_cast<int>(habFloat);
				if (patchModel) {
					patchFloat = badPatchFloat;
					ifsPatchMap >> patchFloat;
					patchCode = static_cast<int>(patchFloat);
				}
#endif
				if (habCode == noDataHabCode) {
		addNewCellToLand(x, y, -1); // add cell only to landscape
	}
				else if (habFloat < 0.0 || habFloat > 100.0) { // invalid quality score
#if RS_RCPP && !R_CMD
			Rcpp::Rcout << "Found invalid habitat quality score." << std::endl;
#endif
					ifsHabMap.close();
					ifsHabMap.clear();
					if (patchModel) {
						ifsPatchMap.close();
						ifsPatchMap.clear();
					}
					return 17;
				}
				else {
					if (patchModel) {
						if (patchCode < 0 || patchCode == noDataPatch) { // invalid patch code
#if RS_RCPP && !R_CMD
							if (patchCode == noDataPatch) Rcpp::Rcout << "Found patch NA in valid habitat cell." << std::endl;
					else Rcpp::Rcout << "Found negative patch ID in valid habitat cell." << std::endl;
#endif
							ifsHabMap.close();
							ifsHabMap.clear();
							ifsPatchMap.close();
							ifsPatchMap.clear();
							return 14;
						}
							if (patchCode == 0) { // cell is in the matrix
								addNewCellToPatch(0, x, y, habFloat);
				}
				else {
								if (existsPatch(patchCode)) {
									pPatch = findPatch(patchCode);
					}
					else {
									pPatch = newPatch(seq++, patchCode);
					}
								addNewCellToPatch(pPatch, x, y, habFloat);
				}
					}
					else { // cell-based model
						// add cell to landscape (patches created later)
						addNewCellToLand(x, y, habFloat);
					}
				}

			} // for x
		} // for y

#if RS_RCPP
		ifsHabMap >> habFloat;
		if (!ifsHabMap.eof()) EOFerrorR(habfile);
		if (patchModel)
		{
			ifsPatchMap >> patchFloat;
			if (!ifsPatchMap.eof()) EOFerrorR(pchfile);
	}
#endif
	break;

default:
	break;
	} // end switch(rasterType)

	if (ifsHabMap.is_open()) {
		ifsHabMap.close();
		ifsHabMap.clear();
	}
	if (ifsPatchMap.is_open()) {
		ifsPatchMap.close();
		ifsPatchMap.clear();
	}

	if (sim.batchMode) {
		if (costfile != "NULL") {
			int retcode = readCosts(costfile);
			if (retcode < 0) return 54;
        }
        if (scalinglayers.size() > 0) {
            if (scalinglayers.size() == nDSlayer) {
                int retcode = readDemographicScaling(scalinglayers);
                if (retcode < 0) return 54; //change number
            }
		}
	}

	return 0;
}

//---------------------------------------------------------------------------

//---------------------------------------------------------------------------

//---------------------------------------------------------------------------

int Landscape::readCosts(string fname)
{

#if RS_RCPP
	wifstream costs; // cost map file input stream
#else
	ifstream costs; // cost map file input stream
#endif

	//int hc,maxYcost,maxXcost,NODATACost,hab;
	int hc, maxYcost, maxXcost, NODATACost;
	float minLongCost, minLatCost; int resolCost;
	float fcost;
#if RS_RCPP
	wstring header;
#else
	string header;
#endif
	Cell* pCell;

	int maxcost = 0;

 // open cost file
// #if RS_RCPP
// 	costs.open(fname, std::ios::binary);
// 	if (costsraster.utf) {
// 		// apply BOM-sensitive UTF-16 facet
// 		costs.imbue(std::locale(costs.getloc(), new std::codecvt_utf16<wchar_t, 0x10ffff, std::consume_header>));
// 	}
// #else
	costs.open(fname.c_str());
// #endif
	// read headers and check that they correspond to the landscape ones
	costs >> header;
#if RS_RCPP
	if (!costs.good()) {
		// corrupt file stream
		StreamErrorR(fname);
		costs.close();
		costs.clear();
		return -181;
	}
	if (header != L"ncols" && header != L"NCOLS") {
#else
	if (header != "ncols" && header != "NCOLS") {
#endif
	costs.close(); costs.clear();
	return -1;
}
double tmpresolCost;
costs >> maxXcost >> header >> maxYcost >> header >> minLongCost;
costs >> header >> minLatCost >> header >> tmpresolCost >> header >> NODATACost;
resolCost = (int) tmpresolCost;

	for (int y = maxYcost - 1; y > -1; y--) {
		for (int x = 0; x < maxXcost; x++) {
#if RS_RCPP
			if (costs >> fcost) {
#else
			costs >> fcost;
#endif
			hc = (int)fcost; // read as float and convert to int
#if RS_RCPP
			}
	else {
		// corrupt file stream
#if RS_RCPP && !R_CMD
		Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
		StreamErrorR(fname);
		costs.close();
		costs.clear();
		return -181;
	}
#endif
	if (hc < 1 && hc != NODATACost) {
#if RS_RCPP && !R_CMD
		Rcpp::Rcout << "Cost map may only contain values of 1 or higher, but found " << fcost << "." << endl;
#endif
		// error - zero / negative cost not allowed
		costs.close(); costs.clear();
		return -999;
	}
	pCell = findCell(x, y);
	if (pCell != 0) { // not no-data cell
	    if (hc > 0){ // only if cost value is  above 0 in a data cell
	        pCell->setCost(hc);
		    if (hc > maxcost) maxcost = hc;
	    } else { // if cost value is below 0
#if RS_RCPP && !R_CMD
	    Rcpp::Rcout << "Cost map may only contain values of 1 or higher in habitat cells, but found " << hc << " in cell x: " << x << " y: " << y << "." << endl;
#endif
		throw runtime_error("Found negative- or zero-cost habitat cell.");
	    }

	} // end not no data cell
		}
	}
#if RS_RCPP
costs >> fcost;
if (costs.eof()) {
#if RS_RCPP && !R_CMD
	Rcpp::Rcout << "Costs map loaded." << endl;
#endif
}
else EOFerrorR(fname);
#endif

costs.close(); costs.clear();

return maxcost;

}

    //---------------------------------------------------------------------------
    int Landscape::readDemographicScaling(vector <string> scalinglayers){
        // Create a temporary landscape to store the vectors for each cell
        // I bet there is a better way to implement it, but I couldn't think of it for now
        // each cell will contain a vector of three floats
        std::vector<std::vector<std::vector<float>>> landscape(dimY,
                                                               std::vector<std::vector<float>>(dimX,
                                                                                               std::vector<float>(scalinglayers.size(), 0.0f))); //length of string is determined by DSlayer

#if RS_RCPP
        wstring header;
#else
        string header;
#endif

        int DSnb = 0; // first position is 0
        int DS;
        int DSnodata;
        float DSfloat;   // float for reading in data from file

#if RS_RCPP
        wifstream DSfile; // DS file input stream
#else
        ifstream DSfile; // DS file input stream
#endif



        // for each element of scalinglayers
        for (const auto& DSlayer : scalinglayers) { //DSlayer is a reference to file string
            // open file
            DSfile.open(DSlayer.c_str());
            // if file couldn't be opened, return a failure message (and the readLandscape function: close all file connections and exit)
            if (!DSfile.is_open()) return -11; // might need to change the code number

            // if it opened
            // skip header, but store no data value
            for (int i = 0; i < 5; i++)
                DSfile >> header >> DSfloat;
            DSfile >> header >> DSnodata; // 6th line

#if RS_RCPP
            if (!DSfile.good()) {
                // corrupt file stream
                StreamErrorR(DSlayer);
                DSfile.close();
                DSfile.clear();
                //do I also need to close the habitat, and potentially patchfile? or is this taken care of?
                return -181; // need to change the code
            }
#endif

            // set badfloat
            float badDSfloat = -9.0; if (DSnodata == -9) DSfloat = -99.0;

            // loop over landscape x+y
            for (int y = dimY - 1; y >= 0; y--) {
                for (int x = 0; x < dimX; x++) {
                    // read in cell value
                    DSfloat = badDSfloat; // set to bad float value
                    // do I need to check this value for any inconsistencies? -> better safe than sorry
#if RS_RCPP
                    if (DSfile >> DSfloat) {
#else
                        DSfile >> DSfloat;
#endif
                        //  DS = (int)DSfloat; I guess non integer values are totally fine?
#if RS_RCPP
                    }
                    else {
                        // corrupt file stream
#if RS_RCPP && !R_CMD
                        Rcpp::Rcout << "At (x,y) = " << x << "," << y << " :" << std::endl;
#endif
                        StreamErrorR(DSlayer);
                        DSfile.close();
                        DSfile.clear();
                        return -134; // need to change number code
                    }
#endif
                    // check if the value is a no data value
                    if (DSfloat == DSnodata) {
                        // if it is, do I need to do anything?
                    }
                    else {
                        if (DSfloat < 0.0 || DSfloat > 100.0) { // check for negative values or values above 100
#if RS_RCPP && !R_CMD
                            Rcpp::Rcout << "Found invalid demographic scaling score." << std::endl;
#endif
                            DSfile.close(); DSfile.clear();
                            return -17; // need to change code
                        }
                        else {
                            // if the value is ok,
                            // set the vector value of this cell at this location to the float value read in:
                            landscape[y][x][DSnb] = DSfloat; // Assignment to specific cell's vector
                        }
                    }
                } // end x loop
            } // end y loop

            // use Rcpp::Rcout to print the landscape at DSnb=0:
            // close file connection
#if RS_RCPP
            DSfile >> DSfloat;
            if (!DSfile.eof()) EOFerrorR(DSlayer);
            else {
#if RS_RCPP && !R_CMD
                Rcpp::Rcout << "Demographic scaling layer " << DSnb + 1 << " loaded." << endl;
#endif
            }
#endif
            if (DSfile.is_open()) { DSfile.close(); DSfile.clear(); }
            DSnb ++; // increase the reference number and read in next file
        }; // end loop over scalinglayers

        // add the cells' vector of DSfloats to the actual landscape using the function addchgDemoScaling
        // loop over all cells
        for (int y = dimY - 1; y >= 0; y--) {
            for (int x = 0; x < dimX; x++) {
                // extract the vector in landscape at x, y
                std::vector<float> localDS = landscape[y][x];
                cells[y][x]->addchgDemoScaling(localDS); // add the vector to the cell
            }
        }
        return 0; // return success code
    };

//---------------------------------------------------------------------------

rasterdata CheckRasterFile(string fname)
{
	rasterdata r;
	string header;
	int inint;
	ifstream infile;

	r.ok = true;
	r.errors = r.ncols = r.nrows = r.cellsize = 0;
	r.xllcorner = r.yllcorner = 0.0;

	infile.open(fname.c_str());
	if (infile.is_open()) {
		infile >> header >> r.ncols;
		if (header != "ncols" && header != "NCOLS") r.errors++;
		infile >> header >> r.nrows;
		if (header != "nrows" && header != "NROWS") r.errors++;
		infile >> header >> r.xllcorner;
		if (header != "xllcorner" && header != "XLLCORNER") r.errors++;
		infile >> header >> r.yllcorner;
		if (header != "yllcorner" && header != "YLLCORNER") r.errors++;
		double tmpcellsize;
		infile >> header >> tmpcellsize;
		r.cellsize = (int) tmpcellsize;
		if (header != "cellsize" && header != "CELLSIZE") r.errors++;
		infile >> header >> inint;
		if (header != "NODATA_value" && header != "NODATA_VALUE") r.errors++;
		infile.close();
		infile.clear();
		if (r.errors > 0) r.ok = false;
	}
	else {
		r.ok = false; r.errors = -111;
	}
	infile.clear();

	return r;
}

//---------------------------------------------------------------------------

// Patch connectivity functions

// Create & initialise connectivity matrix
void Landscape::createConnectMatrix(void)
{
	if (connectMatrix != 0) deleteConnectMatrix();
	int npatches = (int)patches.size();
	connectMatrix = new int* [npatches];
	for (int i = 0; i < npatches; i++) {
		connectMatrix[i] = new int[npatches];
		for (int j = 0; j < npatches; j++) connectMatrix[i][j] = 0;
	}
}

// Re-initialise connectivity matrix
void Landscape::resetConnectMatrix(void)
{
	if (connectMatrix != 0) {
		int npatches = (int)patches.size();
		for (int i = 0; i < npatches; i++) {
			for (int j = 0; j < npatches; j++) connectMatrix[i][j] = 0;
		}
	}
}

// Increment connectivity count between two specified patches
void Landscape::incrConnectMatrix(int p0, int p1) {
	int npatches = (int)patches.size();
	if (connectMatrix == 0 || p0 < 0 || p0 >= npatches || p1 < 0 || p1 >= npatches) return;
	connectMatrix[p0][p1]++;
}

// Delete connectivity matrix
void Landscape::deleteConnectMatrix(void)
{
	if (connectMatrix != 0) {
		int npatches = (int)patches.size();
		for (int j = 0; j < npatches; j++) {
			if (connectMatrix[j] != 0)
				delete connectMatrix[j];
		}
		delete[] connectMatrix;
		connectMatrix = 0;
	}
}

// Close connectivity file
bool Landscape::outConnectFinishLandscape()
{
		if (outConnMat.is_open()) outConnMat.close();
		outConnMat.clear();
		return true;
	}

// Open connectivity file and write header record
bool Landscape::outConnectStartLandscape()
{
	simParams sim = paramsSim->getSim();

	string name = paramsSim->getDir(2);
	if (sim.batchMode) {
		name += "Batch" + to_string(sim.batchNum) + "_";
		name += "Sim" + to_string(sim.simulation) + "_Land" + to_string(landNum);
	}
	else
		name += "Sim" + to_string(sim.simulation);
	name += "_Connect.txt";
	outConnMat.open(name.c_str());

	outConnMat << "Rep\tYear\tStartPatch\tEndPatch\tNinds" << endl;

	return outConnMat.is_open();
}

#if RS_RCPP
// Close movement paths file
void Landscape::outPathsFinishReplicate()
{
		if (outMovePaths.is_open()) outMovePaths.close();
		outMovePaths.clear();
	}

// Open movement paths file and write header record
void Landscape::outPathsStartReplicate(int rep)
{
		simParams sim = paramsSim->getSim();
		string name = paramsSim->getDir(2);
		if (sim.batchMode) {
			name += "Batch" + to_string(sim.batchNum)
				+ "_Sim" + to_string(sim.simulation)
				+ "_Land" + to_string(landNum)
				+ "_Rep" + to_string(rep);
		}
		else {
			name += "Sim" + to_string(sim.simulation)
				+ "_Rep" + to_string(rep);
		}
		name += "_MovePaths.txt";

		outMovePaths.open(name.c_str());
		if (outMovePaths.is_open()) {
			outMovePaths << "Year\tIndID\tStep\tx\ty\tStatus" << endl;
		}
		else {
			outMovePaths.clear();
	}
}
#endif

void Landscape::outConnect(int rep, int yr)
{
	int patchnum0, patchnum1;
	int npatches = (int)patches.size();
	int* emigrants = new int[npatches]; // 1D array to hold emigrants from each patch
	int* immigrants = new int[npatches]; // 1D array to hold immigrants to  each patch

	for (int i = 0; i < npatches; i++) {
		emigrants[i] = immigrants[i] = 0;
	}

	for (int i = 0; i < npatches; i++) {
		patchnum0 = patches[i]->getPatchNum();
		if (patchnum0 != 0) {
			for (int j = 0; j < npatches; j++) {
				patchnum1 = patches[j]->getPatchNum();
				if (patchnum1 != 0) {
					emigrants[i] += connectMatrix[i][j];
					immigrants[j] += connectMatrix[i][j];
					if (connectMatrix[i][j] > 0) {
						outConnMat << rep << "\t" << yr
							<< "\t" << patchnum0 << "\t" << patchnum1
							<< "\t" << connectMatrix[i][j] << endl;
					}
				}
			}
		}
	}

	for (int i = 0; i < npatches; i++) {
		patchnum0 = patches[i]->getPatchNum();
		if (patchnum0 != 0) {
			if (patches[i]->getK() > 0.0)
			{ // suitable patch
				outConnMat << rep << "\t" << yr
					<< "\t" << patchnum0 << "\t-999\t" << emigrants[i] << endl;
				outConnMat << rep << "\t" << yr
					<< "\t-999\t" << patchnum0 << "\t" << immigrants[i] << endl;
			}
		}
	}

	delete[] emigrants;
	delete[] immigrants;

}

//---------------------------------------------------------------------------

void Landscape::resetVisits(void) {
	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] != 0) { // not a no-data cell
				cells[y][x]->resetVisits();
			}
		}
	}
}

// Save SMS path visits map to raster text file
void Landscape::outVisits(int rep, int landNr) {

	string name;
	simParams sim = paramsSim->getSim();

	if (sim.batchMode) {
		name = paramsSim->getDir(3)
#if RS_RCPP
			+ "Batch" + to_string(sim.batchNum) + "_"
			+ "Sim" + to_string(sim.simulation)
			+ "_Land" + to_string(landNr) + "_Rep" + to_string(rep)
#else
			+ "Batch" + to_string(sim.batchNum) + "_"
			+ "Sim" + to_string(sim.simulation)
			+ "_land" + to_string(landNr) + "_rep" + to_string(rep)
#endif
			+ "_Visits.txt";
	}
	else {
		name = paramsSim->getDir(3)
			+ "Sim" + to_string(sim.simulation)
			+ "_land" + to_string(landNr) + "_rep" + to_string(rep)
			+ "_Visits.txt";
	}
	outvisits.open(name.c_str());

	outvisits << "ncols " << dimX << endl;
	outvisits << "nrows " << dimY << endl;
	outvisits << "xllcorner " << minEast << endl;
	outvisits << "yllcorner " << minNorth << endl;
	outvisits << "cellsize " << resol << endl;
	outvisits << "NODATA_value -9" << endl;

	for (int y = dimY - 1; y >= 0; y--) {
		for (int x = 0; x < dimX; x++) {
			if (cells[y][x] == 0) { // no-data cell
				outvisits << "-9 ";
			}
			else {
				outvisits << cells[y][x]->getVisits() << " ";
			}
		}
		outvisits << endl;
	}

	outvisits.close(); 
	outvisits.clear();
}

//---------------------------------------------------------------------------

#ifdef UNIT_TESTS
// Tests only: shortcut setup utilities

Landscape createLandscapeFromCells(vector<Cell*> cells, const landParams& lp, Species sp) {
	// Set up landscape
	Landscape ls;
	ls.setLandParams(lp, true);
	// Add cells
	ls.setCellArray();
	for (auto c : cells) {
		ls.addCellToLand(c);
	}
	ls.allocatePatches(&sp);
	return ls;
}

landParams createDefaultLandParams(const int& dim) {

	landParams ls_params;
	ls_params.dimX = ls_params.dimY = dim;
	ls_params.minX = ls_params.minY = 0;
	ls_params.maxX = ls_params.maxY = ls_params.dimX - 1;
	ls_params.resol = ls_params.spResol = 1;
	ls_params.rasterType = 0; // habitat types

	ls_params.patchModel = false;
	ls_params.spDist = false;
	ls_params.generated = false;
	ls_params.dynamic = false;
	ls_params.landNum = 0;
	ls_params.nHab = ls_params.nHabMax = 0; // irrelevant for habitat codes
	return ls_params;
}

void testLandscape() {
	// test coordinate system...
}
#endif // UNIT_TESTS

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------