
///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                     Tool Library:                     //
//                        SCIMAP                         //
//  Hydrological Connectivity Indices Collection (HCIC)  //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                     netIndex.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 "netIndex.h"

#include <vector>
#include <algorithm>


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

//---------------------------------------------------------
CHCIC::CHCIC(void)
{
	Set_Name		("HCIC - Hydrological Connectivity Indices Collection");

	Set_Author		("Sim M. Reaney (c) 2024-25");

	Set_Description	(_TW(
		"This is a collection of different hydrological conenctivity indicies which repersent how different parts of the landscape connect to the rivers or lakes."
		"<ul>"
		"<li>Network Index (see Lane et al. 2004, Lane et al. 2009)"
		"<li>Percentage Downslope Saturated Length" //, see XXX"
		"</ul>"
	));

	Add_Reference("Lane, S.N. and Brookes, C.J. and Kirkby, M.J. and Holden, J.", "2004",
		"A network-index based version of TOPMODEL for use with high-resolution digital topographic data",
		"Hydrological Processes, 18 (1). pp. 191-201."
	);


	Add_Reference("Lane S., Reaney S., Heathwaite A.", "2009",
		"Representation of landscape hydrological connectivity using a topographically driven surface flow index",
		"WATER RESOURCES RESEARCH, 2009. Volume 45, Issue 8.",
		SG_T("https://doi.org/10.1029/2008WR007336")
	);

	//-----------------------------------------------------
	// Data Input/Output...

	Parameters.Add_Grid("", "DEM"     , _TL("Digital Elevation Model"  ), _TL(""), PARAMETER_INPUT);
	Parameters.Add_Grid("", "RAINFALL", _TL("Rainfall Pattern"         ), _TL(""), PARAMETER_INPUT_OPTIONAL);

	Parameters.Add_Grid("", "SLOPE"   , _TL("Slope"                    ), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
	Parameters.Add_Grid("", "CAREA"   , _TL("Catchment Area"           ), _TL(""), PARAMETER_OUTPUT_OPTIONAL);
	Parameters.Add_Grid("", "WETNESS" , _TL("Connectivity"             ), _TL(""), PARAMETER_OUTPUT);

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

	Parameters.Add_Double("", "CRITAREA", _TL("Critical Area for Channels (m2)"), _TL(""),
		5., 80000., 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")), 2
	);

	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_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_Bool("",
		"RESCALE_WET", _TL("Rescale wetness index?"),
		_TL("Rescale wetness index?"),
		1
	);
}


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

//---------------------------------------------------------
int CHCIC::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);
	}

	if( pParameter->Cmp_Identifier("CONNINDEX") )
	{
		pParameters->Set_Enabled("RESCALE_WET", pParameter->asInt() == 0);
	}

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


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

//---------------------------------------------------------
bool CHCIC::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("Rainfall Weighted 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("WEIGHTS"    , Parameters("RAINFALL"))
		&& 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());
	}

	switch( Parameters("WETNESSINDEXMETHOD")->asInt() )
	{
	default: // 'Classic TWI'
		SG_RUN_TOOL_ExitOnError("ta_hydrology", 20,
			   SG_TOOL_PARAMETER_SET("SLOPE"  , pSlope)
			&& SG_TOOL_PARAMETER_SET("AREA"   , pArea)
			&& SG_TOOL_PARAMETER_SET("TWI"    , pTWI)
			&& SG_TOOL_PARAMETER_SET("CONV"   , 0)
		);
		break;

	case  1: // 'SAGA TWI'
		SG_RUN_TOOL_ExitOnError("ta_hydrology", 15,
			   SG_TOOL_PARAMETER_SET("DEM"    , pDEM)
			&& SG_TOOL_PARAMETER_SET("TWI"    , pTWI)
		);
		break;
	}


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

	switch( Parameters("CONNINDEX")->asInt() )
	{
	default: // Network Index of Connectivity
		netwet(pDEM, pTWI, pArea);

		if( Parameters("RESCALE_WET")->asBool() ) // reclass wetness with 5th and 95th percentiles
		{
			rescale_to_percentiles(pTWI);
		}
		break;

	case  1: // Percentage Downslope Saturated Length
		rescale_to_percentiles(pTWI);

		downslopeSaturatedLength(pDEM, pTWI, pArea);
		break;

	case  2:
		Error_Set("missing function requested => 'connectivityIndex()'!"); return( false );
		//connectivityIndex();
		//if( rescale )
		//{
		//	rescale_to_percentiles(pTWI);
		//}
		//break;
	}	

	//-----------------------------------------------------
	// clip erosion risk and netwet to DEM valid data
	#pragma omp parallel
	for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
	{
		if( pDEM->is_NoData(x, y) )
		{
			pTWI->Set_NoData(x, y);
		}
	}

	return( true );
}


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

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

	CSG_Grid netwet(Get_System()), D8(dem);
	double critAreaForChannels = Parameters("CRITAREA")->asDouble();

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

	//setting catchment area to 1 / nodata to represent channels
	for (int y = 0; y < Get_NY() && Set_Progress(y); y++) {
		for (int x = 0; x < Get_NX(); x++) { //iterate through x and y
			if (channels->asDouble(x, y) > critAreaForChannels)
				channels->Set_Value(x, y, 1);
			else
				channels->Set_Value(x, y, channelsNoData);
		}
	}

	#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)) {
				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);
		}
	}

	Process_Set_Text("Network Index"); //progress message
	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
		}


	//set cells that still == 9999999 to nodata
	// 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);

	//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++)
			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 CHCIC::downslopeSaturatedLength(CSG_Grid *dem, CSG_Grid *wet, CSG_Grid *channels)
{
	//Flow direction (D8)
	Process_Set_Text("Downslope Saturated Length: calc D8"); //progress message

	CSG_Grid netwet(Get_System()), D8(Get_System());

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

	Process_Set_Text("Pre-processing for calculating the Network Index"); //progress message
	double critAreaForChannels = Parameters("CRITAREA")->asDouble();

//setting catchment area to 1 / nodata to repersent channels
	for (int y = 0; y < Get_NY() && Set_Progress(y); y++) {
		for (int x = 0; x < Get_NX(); x++) { //iterate through x and y
			if (channels->asDouble(x, y) > critAreaForChannels)
				channels->Set_Value(x, y, 1);
			else
				channels->Set_Value(x, y, channelsNoData);
		}
	}

	for (int y = 0; y < Get_NY() && Set_Progress(y); y++) {
		for (int x = 0; x < Get_NX(); x++) { //iterate through x and y

			if (channels->asInt(x, y) == channelsNoData) {
				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); //set the value of the grid D8 at location x,y to d
			}
			else
				D8.Set_NoData(x, y); //otherwise set the grid value to nodata
		}
	}

	Process_Set_Text("Percent Downslope Saturated Length"); //progress message
	// int n = 0;
	const int maxPathLength = Get_NX() * Get_NY(); //ie it goes through every cell once
	for (int j = 0; j < Get_NY() && Set_Progress(j); j++)
	{
		#pragma omp parallel for schedule(dynamic)
		for (int i = 0; i < Get_NX(); i++)
		{
			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)
							netwet.Set_Value(i, j, 1);
						else
							netwet.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
		}
	}


	//set cells that still == 9999999 to nodata
	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);

	for (int i = 0; i < Get_NX(); i++)
		for (int j = 0; j < Get_NY(); j++)
			wet->Set_Value(i, j, netwet.asDouble(i, j));

	wet->Set_Description("Percentage Downslope Saturated Length index of hydrological connectivity.");
	wet->Set_Name("Percentage Downslope Saturated Length");

	return(true);
}


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

//---------------------------------------------------------
bool CHCIC::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 CHCIC::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 );
}


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

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