
///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                     Tool Library:                     //
//                        SCIMAP                         //
//  Hydrological Connectivity Indices Collection (HCIC)  //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                      SCIMAP.cpp                       //
//                                                       //
//                 Copyright (C) 2025 by                 //
//                      Sim Reaney                       //
//               sim.reaney@durham.ac.uk                 //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// This file is part of 'SAGA - System for Automated     //
// Geoscientific Analyses'. SAGA 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; version 2 of the License.   //
//                                                       //
// SAGA 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 this program; if not,       //
// write to the Free Software Foundation, Inc.,          //
// 59 Temple Place - Suite 330, Boston, MA 02111-1307,   //
// USA.                                                  //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "SCIMAP.h"

#include <vector>
#include <algorithm>

using namespace std;


///////////////////////////////////////////////////////////
//                                                       //
//                                                       //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CSCIMAP::CSCIMAP(void)
{
	Set_Name		(_TL("SCIMAP - Fine Sediment"));

	Set_Author		("Sim M. Reaney, Durham University 2005-2025");

	Set_Description	(_TW(
		"SCIMAP Risk Maps Generation. SCIMAP considers the spatial pattern of connectivity and source risks to create risk maps of the diffuse pollution. SCIMAP works in relative terms across the landscape to identify the most probable sources of the observed problem. "
	));

	Add_Reference("Reaney S., Lane S., Heathwaite A. and Dugdale L.", "2011",
		 "Risk-based modelling of diffuse land use impacts from rural landscapes upon salmonid fry abundance",
		"Ecological Modelling, 22, (4), 1016-1029.",
		SG_T("https://doi.org/10.1016/j.ecolmodel.2010.08.022")
	);

	//-----------------------------------------------------
	Parameters.Add_Grid         ("", "DEM"              , _TL("Digital Elevation Model"  ), _TL(""), PARAMETER_INPUT );
	Parameters.Add_Grid_or_Const("", "ERODIBILITY"      , _TL("Erodibility"              ), _TL(""), 1., 0., true);
	Parameters.Add_Grid ("", "RAINFALL"         , _TL("Rainfall Pattern"         ), _TL(""), PARAMETER_INPUT);

	Parameters.Add_Grid("", "SLOPE"            , _TL("Slope"                    ), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
	Parameters.Add_Grid("", "CAREA"            , _TL("Catchment Area"           ), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
	Parameters.Add_Grid("", "CHANNELS"         , _TL("River Channels"           ), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
	Parameters.Add_Grid("", "EROSRISK"         , _TL("Erosion Risk"             ),
		_TL("This grid is the potential for material to be exported and connected to the river channels. It is the physical mobilisation risk mutipled by the connectivity"),
		PARAMETER_OUTPUT
	);

	Parameters.Add_Grid("", "WETNESS"          , _TL("Hydrological Connectivity"),
		_TL("This grid is relative strength of the connection between each point in the landscape and the river or lake that it drains to. Scaled between zero and one."),
		PARAMETER_OUTPUT
	);

	Parameters.Add_Grid("", "EROSRISKFACC"     , _TL("Accumulated Erosion Risk" ), _TL(""), PARAMETER_OUTPUT);
	Parameters.Add_Grid("", "EROSRISK_TOP" , _TL("Top 10% Erosion Risk"     ), _TL(""), PARAMETER_OUTPUT);
	// Parameters.Add_Grid("", "CHANNELSRISK"     , _TL("Erosion Risk in Channels" ), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
	Parameters.Add_Grid("", "CHANNELSRISKCONCN", _TL("Erosion Risk in Channels Concn."), _TL(""), PARAMETER_OUTPUT);

	//-----------------------------------------------------
	// Options...

	Parameters.Add_Double("", "CRITAREA", _TL("Critical Area for Channels (m²)"), _TL(""),
		5., 800000., true
	);

//	Parameters.Add_Choice("", "SLOPE_METHOD", _TL("Slope Calculation Method"), _TL(""), CSG_String::Format("%s|%s|%s|%s|%s|%s|%s",
//		_TL("Maximum Slope (Travis et al. 1975)"),
//		_TL("Maximum Triangle Slope (Tarboton 1997)"),
//		_TL("Least Squares Fit Plane (Costa-Cabral & Burgess 1996)"),
//		_TL("Fit 2.Degree Polynom (Bauer, Rohdenburg, Bork 1985)"),
//		_TL("Fit 2.Degree Polynom (Heerdegen & Beran 1982)"),
//		_TL("Fit 2.Degree Polynom (Zevenbergen & Thorne 1987)"),
//		_TL("Fit 3.Degree Polynom (Haralick 1983)")), 5
//	);

	Parameters.Add_Choice("", "PREPROCESS"   , _TL("DEM Pre-Processing"), _TL(""), CSG_String::Format("%s|%s|%s|%s|%s",
		_TL("Deepen Drainage Routes"),
		_TL("Fill Sinks"),
		_TL("Fill Sinks XXL (Wang & Liu)"),
		_TL("Fill Sinks (Planchon/Darboux, 2001)"),
		_TL("none")), 0
	);

	Parameters.Add_Choice("", "FLOWDIRMETHOD", _TL("Flow Direction Method"), _TL(""), CSG_String::Format("%s|%s|%s|%s|%s|%s|%s",
		_TL("Deterministic 8"),
		_TL("Rho 8"),
		_TL("Braunschweiger Reliefmodell"),
		_TL("Deterministic Infinity"),
		_TL("Multiple Flow Direction"),
		_TL("Multiple Triangular Flow Directon"),
		_TL("Multiple Maximum Downslope Gradient Based Flow Directon")), 4
	);

	Parameters.Add_Double("FLOWDIRMETHOD", "CONVERGENCE", _TL("Convergence"), _TL("Convergence factor for Multiple Flow Direction Algorithm (Freeman 1991)"), 1.1, 0., true);
//	Parameters.Add_Bool  ("FLOWDIRMETHOD", "LINEAR_DO"  , _TL("Linear Flow"), _TL("Use D8 if catchment area becomes higher than specified threshold.", false);
//	Parameters.Add_Double("LINEAR_DO"    , "LINEAR_MIN" , _TL("Threshold"  ), _TL("Use D8 if catchment area becomes higher than specified threshold (Cells).", 500.);

	Parameters.Add_Choice("", "WETNESSINDEXMETHOD", _TL("Wetness Index Method"), _TL(""), CSG_String::Format("%s|%s",
		_TL("Classic TWI (Kirkby 1976)"),
		_TL("SAGA Wetness Index")), 0
	);

	Parameters.Add_Choice("", "CONNINDEX", _TL("Connectivity Index"), _TL(""), CSG_String::Format("%s|%s",
		_TL("Network Index of Connectivity"),
		_TL("Percentage Downslope Saturated Length")), 0
	);

	Parameters.Add_Choice("", "STREAMPOWERMETHOD", _TL("Stream Power"), _TL("Decide how to include stream power as part of erosivity"), CSG_String::Format("%s [%s]|%s [%s]|%s [%s]",
		_TL("Stream Power"                       ), SG_T("tan(slope) * area"),
		_TL("Simplified Stream power"            ), SG_T("tan(slope) * (area / cellsize)^0.5"),
		_TL("Stream Power Specifc Catchment Area"), SG_T("tan(slope) * (area / cellsize)")),  0
	);

}


///////////////////////////////////////////////////////////
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
int CSCIMAP::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
{
	if( pParameter->Cmp_Identifier("FLOWDIRMETHOD") ) // flow direction method changed!
	{
		pParameters->Set_Enabled("CONVERGENCE", pParameter->asInt() == 4 || pParameter->asInt() == 5);
	}

	return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
}


///////////////////////////////////////////////////////////
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSCIMAP::On_Execute(void)
{
	CSG_Grid DEM, *pDEM = Parameters("DEM")->asGrid();

	if( Parameters("PREPROCESS")->asInt() < 4 ) // do pre-processing!
	{
		Process_Set_Text(_TL("DEM Pre-Processing"));

		pDEM = &DEM; DEM.Create(Get_System());
	}

	switch( Parameters("PREPROCESS")->asInt() )
	{
	case  0: // Sink Removal - Deepen Drainage Routes
		SG_RUN_TOOL_ExitOnError("ta_preprocessor" , 2,
		   SG_TOOL_PARAMETER_SET("DEM"        , Parameters("DEM"))
		&& SG_TOOL_PARAMETER_SET("DEM_PREPROC", pDEM)
		&& SG_TOOL_PARAMETER_SET("METHOD"     , 0)
		);
		break;

	case  1: // Sink Removal - Fill Sinks
		SG_RUN_TOOL_ExitOnError("ta_preprocessor" , 2,
		   SG_TOOL_PARAMETER_SET("DEM"        , Parameters("DEM"))
		&& SG_TOOL_PARAMETER_SET("DEM_PREPROC", pDEM)
		&& SG_TOOL_PARAMETER_SET("METHOD"     , 1)
		);
		break;

	case  2: // Fill Sinks XXL (Wang & Liu)
		SG_RUN_TOOL_ExitOnError("ta_preprocessor" , 5,
		   SG_TOOL_PARAMETER_SET("ELEV"       , Parameters("DEM"))
		&& SG_TOOL_PARAMETER_SET("FILLED"     , pDEM)
		);
		break;

	case  3: // Fill Sinks (Planchon/Darboux, 2001)
		SG_RUN_TOOL_ExitOnError("ta_preprocessor" , 3,
		   SG_TOOL_PARAMETER_SET("DEM"        , Parameters("DEM"))
		&& SG_TOOL_PARAMETER_SET("RESULT"     , pDEM)
		);
		break;

	default: // keep on going with the supplied DEM as it is!
		break;
	}

	//-----------------------------------------------------
	Process_Set_Text(_TL("Slope"));

	CSG_Grid Slope, *pSlope = Parameters("SLOPE")->asGrid();

	if( !pSlope )
	{
		pSlope = &Slope; Slope.Create(Get_System());
	}

	SG_RUN_TOOL_ExitOnError("ta_morphometry"  , 0, // Slope, Aspect, Curvatures
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , Parameters("DEM"))
		&& SG_TOOL_PARAMETER_SET("SLOPE"      , pSlope)
	);

	//-----------------------------------------------------
	Process_Set_Text(_TL("Catchment Area"));

	CSG_Grid Area, *pArea = Parameters("CAREA")->asGrid();

	if( !pArea )
	{
		pArea = &Area; Area.Create(Get_System());
	}

	SG_RUN_TOOL_ExitOnError("ta_hydrology"    , 0,
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , pDEM)
		&& SG_TOOL_PARAMETER_SET("FLOW"       , pArea)
		&& SG_TOOL_PARAMETER_SET("METHOD"     , Parameters("FLOWDIRMETHOD"))
		&& SG_TOOL_PARAMETER_SET("CONVERGENCE", Parameters("CONVERGENCE"))
	);

	//-----------------------------------------------------
	Process_Set_Text(_TL("Channel Network"));

	CSG_Grid Channels, *pChannels = Parameters("CHANNELS")->asGrid(), Route(Get_System()); CSG_Shapes Shapes;

	if( !pChannels )
	{
		pChannels = &Channels; Channels.Create(Get_System());
	}

	SG_RUN_TOOL_ExitOnError("ta_channels"     , 0,
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , pDEM)
		&& SG_TOOL_PARAMETER_SET("CHNLNTWRK"  , pChannels)
		&& SG_TOOL_PARAMETER_SET("CHNLROUTE"  , &Route)
		&& SG_TOOL_PARAMETER_SET("SHAPES"     , &Shapes)
		&& SG_TOOL_PARAMETER_SET("INIT_METHOD", 2)
		&& SG_TOOL_PARAMETER_SET("INIT_GRID"  , pArea)
		&& SG_TOOL_PARAMETER_SET("INIT_VALUE" , Parameters("CRITAREA")->asDouble())
	);

	//-----------------------------------------------------
	Process_Set_Text(_TL("Stream Power"));

	CSG_Grid *erosRisk = Parameters("EROSRISK")->asGrid();

	int Method = Parameters("STREAMPOWERMETHOD")->asInt(); 
	double Threshold = Parameters("CRITAREA")->asDouble();

	for(int y=0; y<Get_NY() && Set_Progress(y); y++)
	{
		#pragma omp parallel for schedule(dynamic)
		for(int x=0; x<Get_NX(); x++)
		{
			if( pSlope->is_NoData(x, y) || pArea->is_NoData(x, y) )
			{
				erosRisk->Set_NoData(x, y);
			}
			else if( Method > 3 ) // do not use stream power
			{
				erosRisk->Set_Value(x, y, 1.);
			}
			else
			{
				double a = pArea->asDouble(x, y);

				if( a > Threshold ) // the critical area for channel init
				{
					erosRisk->Set_Value(x, y, 0.);
				}
				else switch( Method )
				{
				case  0: // Stream Power, tan(slope) * area
					erosRisk->Set_Value(x, y, tan(pSlope->asDouble(x, y)) * a);
					break;

				case  1: // Simplified Stream power, tan(slope) * (area / cellsize)^0.5
					erosRisk->Set_Value(x, y, tan(pSlope->asDouble(x, y)) * sqrt(a / Get_Cellsize()));
					break;

				case  2: // Stream Power Specifc Catchment Area, tan(slope) * area / cellsize
					erosRisk->Set_Value(x, y, tan(pSlope->asDouble(x, y)) * a / Get_Cellsize());
					break;
				}
			}
		}
	}

	if( Method <= 3 )
	{
		rescale_to_percentiles(erosRisk);
	}

	//-----------------------------------------------------
	Process_Set_Text(_TL("Rescale rainfall data"));

	CSG_Grid Rain(*Parameters("RAINFALL")->asGrid()); // copy construct! we should to keep orignal rain fall data untouched!

	if( Rain.Get_Mean() > 0. )
	{
		Rain.Multiply(1. / Rain.Get_Mean());
	}
	else
	{
		Error_Set(_TL("invalid rainfall data"));

		return( false );
	}

	//-----------------------------------------------------
	Process_Set_Text(_TL("Rainfall Weighted Catchment Area"));

	SG_RUN_TOOL_ExitOnError("ta_hydrology"    , 0,
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , pDEM)
		&& SG_TOOL_PARAMETER_SET("FLOW"       , pArea) // now become rain fall pattern weighted!
		&& SG_TOOL_PARAMETER_SET("WEIGHTS"    , &Rain)
		&& SG_TOOL_PARAMETER_SET("METHOD"     , Parameters("FLOWDIRMETHOD"))
		&& SG_TOOL_PARAMETER_SET("CONVERGENCE", Parameters("CONVERGENCE"))
	);

	//-----------------------------------------------------
	Process_Set_Text(_TL("Topographic Wetness Index"));

	CSG_Grid TWI, *pTWI = Parameters("WETNESS")->asGrid();

	if( !pTWI )
	{
		pTWI = &TWI; TWI.Create(Get_System());
	}

	SG_RUN_TOOL_ExitOnError("ta_hydrology", Parameters("WETNESSINDEXMETHOD")->asInt() == 0 ? 20 : 15, // 'Classic TWI' or 'SAGA TWI'
		   SG_TOOL_PARAMETER_SET("SLOPE" , pSlope)
		&& SG_TOOL_PARAMETER_SET("AREA"  , pArea) // << specific catchment area
		&& SG_TOOL_PARAMETER_SET("TWI"   , pTWI)
		&& SG_TOOL_PARAMETER_SET("CONV"  , 0)
	);

	//-----------------------------------------------------
	Process_Set_Text(_TL("Connectivity Index")); //progress message

	switch( Parameters("CONNINDEX")->asInt() )
	{
	default: // Network Index of Connectivity
		netwet(pDEM, pTWI, pArea);
		rescale_to_percentiles(pTWI); //default is 5th and 95th percentiles
		break;

	case  1: // Percentage Downslope Saturated Length
		rescale_to_percentiles(pTWI); //default is 5th and 95th percentiles
		downslopeSaturatedLength(pDEM, pTWI, pChannels);
		break;
	}	

	//-----------------------------------------------------
	Process_Set_Text("Fine sediment risk per cell");

	CSG_Grid *pErodibility = Parameters("ERODIBILITY")->asGrid  ();
	double     Erodibility = Parameters("ERODIBILITY")->asDouble();
	//CSG_Grid *connRisk = pTWI; // Parameters("WETNESS")->asGrid(); //connectivity information is now stored in the WETNESS grid to save RAM. 

	// Calculate erosion risk
	#pragma omp parallel for schedule(dynamic)
	for(int x=0; x<Get_NX(); x++)
	{
		for(int y=0; y<Get_NY(); y++)
		{
			if( pDEM->is_NoData(x, y) ) // clip erosion risk and connectivity to DEM valid data
			{
				erosRisk->Set_NoData(x, y);
				pTWI    ->Set_NoData(x, y);
			}
			else if( erosRisk->is_NoData(x, y) == false && pChannels->is_NoData(x, y) == true )
			{
				const double stpow = erosRisk->asDouble(x, y);
				const double conn  = pTWI    ->asDouble(x, y);
				const double erod  = pErodibility ? pErodibility->asDouble(x, y) : Erodibility;

				if( stpow < 0 || stpow > 1 || erod < 0 || erod > 1 || conn < 0 || conn > 1 )
				{
					erosRisk->Set_Value(x, y, 0.);
				}
				else
				{
					erosRisk->Set_Value(x, y, stpow * erod * conn);
				}
			}
		}
	}

	//-----------------------------------------------------
	Process_Set_Text(_TL("Accumulate risk through network"));

	CSG_Grid *erosRiskFacc = Parameters("EROSRISKFACC")->asGrid();

	SG_RUN_TOOL_ExitOnError("ta_hydrology"      , 0,
		   SG_TOOL_PARAMETER_SET("ELEVATION"    , pDEM)	// << preprocessed DEM
		&& SG_TOOL_PARAMETER_SET("FLOW"         , erosRiskFacc)
		&& SG_TOOL_PARAMETER_SET("WEIGHTS"      , erosRisk)
		&& SG_TOOL_PARAMETER_SET("METHOD"       , Parameters("FLOWDIRMETHOD"))
	);

	//-----------------------------------------------------
	Process_Set_Text(_TL("Calculating risk in the channel network"));

	//give results for channel cells
	//divide by catchment area

	
	// CSG_Grid *channelRisk      = Parameters("CHANNELSRISK"     )->asGrid();
	CSG_Grid *channelRiskConcn = Parameters("CHANNELSRISKCONCN")->asGrid();

	for (int y = 0; y < Get_NY() && Set_Progress(y); y++)
	{
		for (int x = 0; x < Get_NX(); x++)
		{
			if( pSlope->is_NoData(x, y) == false && pArea->asDouble(x, y) > 0)
			{
				double channelID = pChannels->asDouble(x, y);

				int channel = channelID > 0 || channelID == -1 ? 1 : 0;

				double er = erosRiskFacc->asDouble(x, y);// accmumlated erosion risk
				double a = pArea->asDouble(x, y); // area m2
				// if (channel == 1)
				// 	channelRisk->Set_Value(x, y, er);
				// else
				// 	channelRisk->Set_NoData(x, y);

				double v2 = (er / a);

				if (channel == 1)
					channelRiskConcn->Set_Value(x, y, v2);
				else
					channelRiskConcn->Set_NoData(x, y);
			}
			else
			{
				// channelRisk     ->Set_NoData(x, y);
				channelRiskConcn->Set_NoData(x, y);
			}
		}
	}

	top_percent(erosRisk, Parameters("EROSRISK_TOP")->asGrid());

	return( true );
}


///////////////////////////////////////////////////////////
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSCIMAP::netwet(CSG_Grid *dem, CSG_Grid *wet, CSG_Grid *channels)
{
	//Flow direction (D8)
	Process_Set_Text("Network Index: calc D8");

	CSG_Grid netwet(dem);
	CSG_Grid D8(dem);

	double nodata = dem->Get_NoData_Value();
	// int channelsNoData = (int)channels->Get_NoData_Value();
	Process_Set_Text("Pre-processing for calculating the Network Index");

	#pragma omp parallel for schedule(dynamic)
	for (int y = 0; y < Get_NY(); y++) {
		for (int x = 0; x < Get_NX(); x++) {
			if (channels->is_NoData(x, y) == false) {
				int d = dem->Get_Gradient_NeighborDir(x, y);//get the flow direction 0-7 for the cell x,y using Get_Grad_Neighb
				D8.Set_Value(x, y, d);
			}
			else
				D8.Set_NoData(x, y); //otherwise set the grid value to nodata
		}
	}

	Process_Set_Text("Network Index");
	int n = 0;
	#pragma omp parallel for schedule(dynamic)
	for (int i = 0; i < Get_NX(); i++){
//		Set_Progress(i);
		for (int j = 0; j < Get_NY(); j++)
			if (dem->asDouble(i, j) != nodata) {
				n++;
				double minWet = wet->asDouble(i, j);
				int cX = i; //start x
				int cY = j; //start y
				bool stop = false;
				int pathLength = 0;
				int maxPathLength = Get_NX() * Get_NY(); //ie it goes through every cell once
				int fdir_nodata = (int)D8.Get_NoData_Value();

				while (stop == false) {
					int fdir = D8.asInt(cX, cY);
					int toX = cX;
					int toY = cY;

					if (fdir == 7 || fdir == 6 || fdir == 5)
						toX--;
					else if (fdir == 1 || fdir == 2 || fdir == 3)
						toX++;

					if (fdir == 7 || fdir == 0 || fdir == 1)
						toY++;
					else if (fdir == 5 || fdir == 4 || fdir == 3)
						toY--;

					//chk we are routing within the grid area
					if (toX >= 0 && toY >= 0 && toX < Get_NX() && toY < Get_NY() && fdir != -1 && pathLength < maxPathLength && fdir != fdir_nodata) {
						if (wet->asDouble(toX, toY) < minWet)
							minWet = wet->asDouble(toX, toY);
						pathLength++;
					}
					else {
						stop = true;
						netwet.Set_Value(i, j, minWet);
					}
					cX = toX;
					cY = toY;
				} //close while loop
			} //close for for if loop
		}

	//Overwrite the wetness index grid with the connectivity results
	#pragma omp parallel for schedule(dynamic)
	for (int i = 0; i < Get_NX(); i++){
		for (int j = 0; j < Get_NY(); j++){
			// if (netwet.asDouble(i, j) == 9999999)
			// 	netwet.Set_Value(i, j, nodata);
			if(dem->is_NoData(i,j))
				wet->Set_NoData(i,j);
			else
				wet->Set_Value(i, j, netwet.asDouble(i, j));
		}
	}

	wet->Set_Description("Network Index (connectivity)");
	wet->Set_Name("Network Index of Connectivity");

return (true);
}



bool CSCIMAP::downslopeSaturatedLength(CSG_Grid *dem, CSG_Grid *wet, CSG_Grid *channels)
{
	//Flow direction (D8)
	Process_Set_Text("Downslope Saturated Length");

	CSG_Grid pdsl(Get_System());
	CSG_Grid D8(Get_System());

	const double nodata = dem->Get_NoData_Value();
	// const int channelsNoData = (int)channels->Get_NoData_Value();

	#pragma omp parallel for schedule(dynamic)
	for (int y = 0; y < Get_NY(); y++) {
		for (int x = 0; x < Get_NX(); x++) { //iterate through x and y
			if (channels->is_NoData(x, y) == false) { //if the cell is a channel cell
				D8.Set_NoData(x, y); //otherwise set the grid value to nodata	
			}
			else {
				int d = dem->Get_Gradient_NeighborDir(x, y);//get the flow direction 0-7 for the cell x,y using Get_Grad_Neighb
				D8.Set_Value(x, y, d);
			}
		}
	}

	Process_Set_Text("Percent Downslope Saturated Length");
	int n = 0;
	const int maxPathLength = Get_NX() * Get_NY(); //ie it goes through every cell once
	#pragma omp parallel for schedule(dynamic)
	for (int i = 0; i < Get_NX(); i++){
		Set_Progress(i);
		for (int j = 0; j < Get_NY(); j++)
			if (dem->asDouble(i, j) != nodata) {
				n++;
				int dsl_nocells = 0;
				const double minWet = wet->asDouble(i, j);
				int cX = i; //start x
				int cY = j; //start y
				bool stop = false;
				int pathLength = 0;

				const int fdir_nodata = static_cast<int>(D8.Get_NoData_Value());

				while (stop == false) {
					const int fdir = D8.asInt(cX, cY);
					int toX = cX;
					int toY = cY;

					switch (fdir) {
						case 7: case 6: case 5: toX--; break;
						case 1: case 2: case 3: toX++; break;
					}

					switch (fdir) {
						case 7: case 0: case 1: toY++; break;
						case 5: case 4: case 3: toY--; break;
					}

					//chk we are routing within the grid area
					if (toX >= 0 && toY >= 0 && toX < Get_NX() && toY < Get_NY() && fdir != -1 && pathLength < maxPathLength && fdir != fdir_nodata) {
						pathLength++;
						if (wet->asDouble(toX, toY) >= minWet)
							dsl_nocells++;
					}
					else {
						stop = true;
						if (pathLength == 0)
							pdsl.Set_Value(i, j, 1);
						else
							pdsl.Set_Value(i, j, static_cast<double>(dsl_nocells) / static_cast<double>(pathLength));
					}
					cX = toX;
					cY = toY;
				} //close while loop
			} //close for for if loop
		}

	#pragma omp parallel for schedule(dynamic)
	for (int i = 0; i < Get_NX(); i++)
		for (int j = 0; j < Get_NY(); j++){
			wet->Set_Value(i, j, pdsl.asDouble(i, j));
		}
			
	wet->Set_Description("Percentage Downslope Saturated Length (connectivity)");
	wet->Set_Name("Downslope Saturated Length");


	return(true);
}


///////////////////////////////////////////////////////////
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CSCIMAP::top_percent(CSG_Grid *pGrid, CSG_Grid *pTop, double Percentile)
{
	pTop->Set_NoData_Value(0);
	pTop->Fmt_Name("Top %.1f Percentile Pathways", 100. - Percentile);
	pTop->Set_Description(CSG_String::Format("Top %.1f Percentile Pathways", 100. - Percentile));
	
	Percentile = pGrid->Get_Percentile(Percentile);

	#pragma omp parallel for schedule(dynamic)
	for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
	{
		if( pGrid->is_NoData(x, y) || pGrid->asDouble(x, y) <= Percentile )
		{
			pTop->Set_NoData(x, y);
		}
		else
		{
			pTop->Set_Value(x, y, 1);
		}
	}
}


///////////////////////////////////////////////////////////
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSCIMAP::rescale_to_range(CSG_Grid *pGrid)
{
	Process_Set_Text("Rescale");

	if( pGrid->Get_Range() > 0. )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( pGrid->is_NoData(x, y) == false )
			{
				double z = pGrid->asDouble(x, y);

				pGrid->Set_Value(x, y, (z - pGrid->Get_Min()) / pGrid->Get_Range());
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CSCIMAP::rescale_to_percentiles(CSG_Grid *pGrid, double min, double max)
{
	Process_Set_Text("Rescale");

	min = pGrid->Get_Percentile(min);
	max = pGrid->Get_Percentile(max);

	if( min < max )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( pGrid->is_NoData(x, y) == false )
			{
				double z = pGrid->asDouble(x, y);

				pGrid->Set_Value(x, y, z <= min ? 0. : z >= max ? 1. : (z - min) / (max - min));
			}
		}

		return( true );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//                                                       //
//                                                       //
//                                                       //
///////////////////////////////////////////////////////////

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