
///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                     Tool Library                      //
//                     ta_hydrology                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//               Index_of_Connectivity.cpp               //
//                                                       //
//                 Copyright (C) 2026 by                 //
//                    Volker Wichmann                    //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// 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, either version 2 of the     //
// License, or (at your option) any later version.       //
//                                                       //
// 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, see   //
// <http://www.gnu.org/licenses/>.                       //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     wichmann@laserdata.at                  //
//                                                       //
//    contact:    Volker Wichmann                        //
//                LASERDATA GmbH                         //
//                Management and analysis of             //
//                laserscanning data                     //
//                Innsbruck, Austria                     //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "Index_of_Connectivity.h"


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

//---------------------------------------------------------
CIndex_of_Connectivity::CIndex_of_Connectivity(void)
{
	Set_Name		(_TL("Index of Connectivity"));

	Set_Author		("V. Wichmann (c) 2026");

	Set_Description	(_TW(
		"The tool allows one to calculate the 'Index of Connectivity' (IC) as first proposed by "
		"Borselli et al. (2008) and later adapted by Cavalli et al. (2012). The IC can be used to "
		"assess the spatial sediment connectivity in a catchment, i.e. how efficiently hillslopes "
		"and channels (or other sinks) are coupled. The index represents the potential connectivity "
		"between the different parts of a watershed. For each cell the upslope (D_up) and downslope "
		"(D_dn) components of connectivity are evaluated:\n"
		"IC = log10(D_up / D_dn)\n"
		"Please refer to the original publications for further details and formulas.\n"
		"When no weights and slope grids are provided, the tool calculates "
		"the IC as proposed by Cavalli et al. (2012): weights are estimated from "
		"local surface roughness and are limited to the range [0.001;1], local slope is estimated "
		"as the maximum slope (D8) and limited to the range [0.005 m/m;1 m/m].\n"
		"The provided flow routing options are those of the 'Flow Accumulation (Top-Down)' tool "
		"(for references see there).\n"
	));

	Add_Reference("Borselli, L., Cassi, P. & Torri, D.", "2008",
		"Prolegomena to sediment and flow connectivity in the landscape: A GIS and field numerical assessment",
		"Catena, 75, 3: 268-277.",
		SG_T("https://doi.org/10.1016/j.catena.2008.07.006"), SG_T("doi:10.1016/j.catena.2008.07.006")
	);

	Add_Reference("Cavalli, M., Trevisani, S., Comiti, F.. & Marchi, L.", "2012",
		"Geomorphometric assessment of spatial sediment connectivity in small Alpine catchments",
		"Geomorphology, 188: 31-41.",
		SG_T("https://doi.org/10.1016/j.geomorph.2012.05.007"), SG_T("doi:10.1016/j.geomorph.2012.05.007")
	);

	//-----------------------------------------------------
	Parameters.Add_Grid("",
		"ELEVATION"    , _TL("Elevation"),
		_TL("The digital terrain model (DTM). Usually a hydologically sound DTM is used."),
		PARAMETER_INPUT
	);
	
	Parameters.Add_Grid("",
		"SINKS"	, _TL("Sinks"),
		_TL("A grid indicating sink (terminal) cells, e.g. the channel network. Cells with valid values (i.e. not NoData) act as sinks that capture all incoming flow and sediment and do not route anything downslope."),
		PARAMETER_INPUT_OPTIONAL
	);

	Parameters.Add_Grid("",
		"WEIGHTS"      , _TL("Weights"),
		_TL("Pre-computed weights [-], greater than zero. If not given, weights are estimated from local surface roughness within the specified kernel size."),
		PARAMETER_INPUT_OPTIONAL
	);

	Parameters.Add_Grid("",
		"SLOPE", _TL("Slope"),
		_TL("Pre-computed slope [m/m], greater than zero. If not given, slope is estimated as the maximum slope (D8)."),
		PARAMETER_INPUT_OPTIONAL
	);

	Parameters.Add_Grid("",
		"IC", _TL("Index of Connectivity"),
		_TL("The index of connectivity."),
		PARAMETER_OUTPUT
	);

	Parameters.Add_Grid("",
		"UPSLOPE", _TL("Upslope Component"),
		_TL("The upslope component."),
		PARAMETER_OUTPUT_OPTIONAL
	);

	Parameters.Add_Grid("",
		"DOWNSLOPE", _TL("Downslope Component"),
		_TL("The downslope component."),
		PARAMETER_OUTPUT_OPTIONAL
	);

	Parameters.Add_Int("",
		"KERNEL_RADIUS", _TL("Radius Smoothing"),
		_TW("The kernel radius [cells] used for DTM smoothing in local surface roughness / weighting computation."),
		2, 1, true
	);

	Parameters.Add_Choice("",
		"FLOW_ROUTING"		, _TL("Flow Routing"),
		_TL("Choose the flow routing method."),
		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 Direction"),
			_TL("Multiple Maximum Downslope Gradient Based Flow Direction")
		), 4
	);

	Parameters.Add_Bool("",
		"LINEAR_FLOW"		, _TL("Thresholded Linear Flow"),
		_TL("Apply linear flow routing (D8) to all cells that have a flow accumulation greater than the specified threshold."),
		false
	);

	Parameters.Add_Int("LINEAR_FLOW",
		"LINEAR_MIN"	, _TL("Linear Flow Threshold"),
		_TL("The flow accumulation threshold [cells] for linear flow routing."),
		500, 0, true
	);

	Parameters.Add_Double("",
		"CONVERGENCE"	, _TL("Convergence"),
		_TL("The convergence factor for the 'Multiple Flow Direction' and 'Multiple Triangular Flow Direction' algorithms."),
		1.1, 0.001, true
	);

	Parameters.Add_Bool("",
		"MFD_CONTOUR"	, _TL("Contour Length"),
		_TL("Include a (pseudo) contour length as additional weighting factor in multiple flow direction routing, reduces flow to diagonal neighbour cells by a factor of 0.71 (see Quinn et al. 1991 for details)."),
		false
	);

}


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

//---------------------------------------------------------
int CIndex_of_Connectivity::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
{
	if (pParameter->Cmp_Identifier("WEIGHTS"))
	{
		pParameters->Set_Enabled("KERNEL_RADIUS", (*pParameters)["WEIGHTS"].asGrid() == NULL);
	}

	if( pParameter->Cmp_Identifier("FLOW_ROUTING") )
	{
		pParameters->Set_Enabled("CONVERGENCE", pParameter->asInt() == 4 || pParameter->asInt() == 5);
		pParameters->Set_Enabled("MFD_CONTOUR", pParameter->asInt() == 4 || pParameter->asInt() == 5);
	}

	if( pParameter->Cmp_Identifier("LINEAR_FLOW") )
	{
		pParameters->Set_Enabled("LINEAR_MIN", pParameter->asBool());
		pParameters->Set_Enabled("LINEAR_VAL", pParameter->asBool());
	}

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


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

//---------------------------------------------------------
bool CIndex_of_Connectivity::On_Execute(void)
{
	CSG_Grid	*pDTM			= Parameters("ELEVATION")->asGrid();
	CSG_Grid	*pSinks			= Parameters("SINKS")->asGrid();
	CSG_Grid	*pWeights		= Parameters("WEIGHTS")->asGrid();
	CSG_Grid	*pSlope			= Parameters("SLOPE")->asGrid();
	CSG_Grid	*pIC			= Parameters("IC")->asGrid();
	CSG_Grid	*pUp			= Parameters("UPSLOPE")->asGrid();
	CSG_Grid	*pDown			= Parameters("DOWNSLOPE")->asGrid();
	int			Radius			= Parameters("KERNEL_RADIUS")->asInt();
	int			FlowRouting		= Parameters("FLOW_ROUTING")->asInt();
	bool		bLinear			= Parameters("LINEAR_FLOW")->asBool();
	double		MinLinear		= Parameters("LINEAR_MIN")->asDouble();
	double		MFD_Converge	= Parameters("CONVERGENCE")->asDouble();
	bool		MFD_bContour	= Parameters("MFD_CONTOUR")->asBool();

	bool		bCalcWeights	= false;
	bool		bCalcSlope		= false;

	SG_UI_Progress_Lock(true);

	CSG_Grid					Residuals, Weights, Slope, D_down;
	CSG_Grid_Cell_Addressor		Kernel;

	Kernel.Set_Radius(Radius, true);

	if( !pWeights )
	{
		bCalcWeights = true;
		pWeights     = &Weights; Weights.Create(Get_System());
		Residuals.Create(Get_System());
	}

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

	CSG_Grid WeightsFL;		WeightsFL.Create(Get_System());
	CSG_Grid FlowAcc;		FlowAcc.Create(Get_System());
	CSG_Grid WeightsMean;	WeightsMean.Create(Get_System());
	CSG_Grid SlopeMean;		SlopeMean.Create(Get_System());

	SG_UI_Progress_Lock(false);
	
	//-----------------------------------------------------
	if( bCalcWeights || bCalcSlope )
	{
		for(int y=0; y<Get_NY() && Process_Get_Okay(); y++)
		{
			#pragma omp parallel for
			for(int x=0; x<Get_NX(); x++)
			{
				if( bCalcWeights )
				{
					double mean;

					if( Get_Stats(pDTM, &Kernel, x, y, mean, false) )
					{
						Residuals.Set_Value(x, y, pDTM->asDouble(x, y) - mean);
					}
					else
					{
						Residuals.Set_NoData(x, y);
					}
				}

				if (bCalcSlope)
				{
					double s;

					if( Get_Slope(pDTM, x, y, s) )
					{
						if( s < 0.005 )		{ s = 0.005; }
						else if( s > 1. )	{ s = 1.; }
					
						pSlope->Set_Value(x, y, s);
					}
					else
					{
						pSlope->Set_NoData(x, y);
					}
				}
			}
		}
	}

	Set_Progress(12.5, 100.);
	
	//-----------------------------------------------------
	if (bCalcWeights)
	{
		for(int y=0; y<Get_NY() && Process_Get_Okay(); y++)
		{
			#pragma omp parallel for
			for(int x=0; x<Get_NX(); x++)
			{
				double stddev;

				if( Get_Stats(&Residuals, &Kernel, x, y, stddev, true) )
				{
					pWeights->Set_Value(x, y, stddev);
				}
				else
				{
					pWeights->Set_NoData(x, y);
				}
			}
		}
		
		Residuals.Destroy();

		Set_Progress(25., 100.);
		
		SG_UI_Progress_Lock(true);

		pWeights->Set_Max_Samples(pWeights->Get_NCells());

		double w_max = pWeights->Get_Max();

		SG_UI_Progress_Lock(false);	

		for(int y=0; y<Get_NY() && Process_Get_Okay(); y++)
		{
			#pragma omp parallel for
			for(int x=0; x<Get_NX(); x++)
			{
				if( !pWeights->is_NoData(x, y) )
				{
					double w = 1 - (pWeights->asDouble(x, y) / w_max);

					if( w < 0.001 ) { w = 0.001; }

					pWeights->Set_Value(x, y, w);
				}
			}
		}

		Set_Progress(37.5, 100.);
	}
	
	//-----------------------------------------------------
	for(int y=0; y<Get_NY() && Process_Get_Okay(); y++)
	{
		#pragma omp parallel for
		for(int x=0; x<Get_NX(); x++)
		{
			if( pWeights->is_NoData(x, y) || pSlope->is_NoData(x, y) )
			{
				WeightsFL.Set_NoData(x, y);
			}
			else
			{
				WeightsFL.Set_Value(x, y, 1 / (pWeights->asDouble(x, y) * pSlope->asDouble(x, y)));
			}
		}
	}

	Set_Progress(50., 100.);
	
	//-----------------------------------------------------
	SG_UI_Progress_Lock(true);

	SG_RUN_TOOL_ExitOnError("ta_hydrology"    , 0,		// Flow Accumulation (Top-Down)
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , pDTM)
		&& SG_TOOL_PARAMETER_SET("SINKS"      , pSinks)
		&& SG_TOOL_PARAMETER_SET("VAL_INPUT"  , pWeights)
		&& SG_TOOL_PARAMETER_SET("FLOW"       , &FlowAcc)
		&& SG_TOOL_PARAMETER_SET("VAL_MEAN"   , &WeightsMean)
		&& SG_TOOL_PARAMETER_SET("METHOD"     , FlowRouting)
		&& SG_TOOL_PARAMETER_SET("CONVERGENCE", MFD_Converge)
		&& SG_TOOL_PARAMETER_SET("LINEAR_DO"  , bLinear) 
		&& SG_TOOL_PARAMETER_SET("LINEAR_MIN" , MinLinear)
		&& SG_TOOL_PARAMETER_SET("MFD_CONTOUR", MFD_bContour)
		&& SG_TOOL_PARAMETER_SET("FLOW_UNIT"  , 1)
	);

	SG_UI_Progress_Lock(false);

	Set_Progress(62.5, 100.);

	//-----------------------------------------------------
	SG_UI_Progress_Lock(true);

	SG_RUN_TOOL_ExitOnError("ta_hydrology"    , 0,		// Flow Accumulation (Top-Down)
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , pDTM)
		&& SG_TOOL_PARAMETER_SET("SINKS"      , pSinks)
		&& SG_TOOL_PARAMETER_SET("VAL_INPUT"  , pSlope)
		&& SG_TOOL_PARAMETER_SET("FLOW"       , &FlowAcc)
		&& SG_TOOL_PARAMETER_SET("VAL_MEAN"   , &SlopeMean)
		&& SG_TOOL_PARAMETER_SET("METHOD"     , FlowRouting)
		&& SG_TOOL_PARAMETER_SET("CONVERGENCE", MFD_Converge)
		&& SG_TOOL_PARAMETER_SET("LINEAR_DO"  , bLinear) 
		&& SG_TOOL_PARAMETER_SET("LINEAR_MIN" , MinLinear)
		&& SG_TOOL_PARAMETER_SET("MFD_CONTOUR", MFD_bContour)
		&& SG_TOOL_PARAMETER_SET("FLOW_UNIT"  , 1)
	);

	SG_UI_Progress_Lock(false);

	Set_Progress(75., 100.);

	//-----------------------------------------------------
	SG_UI_Progress_Lock(true);

	if( !pDown )
	{
		pDown = &D_down; D_down.Create(Get_System());
	}

	SG_RUN_TOOL_ExitOnError("ta_hydrology"    , 27,		// Maximum Flow Path Length
		   SG_TOOL_PARAMETER_SET("ELEVATION"  , pDTM)
		&& SG_TOOL_PARAMETER_SET("SINKS"      , pSinks)
		&& SG_TOOL_PARAMETER_SET("WEIGHTS"    , &WeightsFL)
		&& SG_TOOL_PARAMETER_SET("DISTANCE"   , pDown)
		&& SG_TOOL_PARAMETER_SET("DIRECTION"  , 0)
	);

	SG_UI_Progress_Lock(false);

	Set_Progress(87.5, 100.);

	//-----------------------------------------------------
	for(int y=0; y<Get_NY() && Process_Get_Okay(); y++)
	{
		#pragma omp parallel for
		for(int x=0; x<Get_NX(); x++)
		{
			if( pDTM->is_NoData(x, y) )
			{
				pIC->Set_NoData(x, y);

				if( pUp )	{ pUp->Set_NoData(x, y); }
				if( pDown )	{ pDown->Set_NoData(x, y); }
			}
			else
			{
				double up = WeightsMean.asDouble(x, y) * SlopeMean.asDouble(x, y) * sqrt(FlowAcc.asDouble(x, y));

				if( pUp )	{ pUp->Set_Value(x, y, up); }

				if( pDown->asDouble(x, y) <= 0. )
				{
					pDown->Set_Value(x, y, WeightsFL.asDouble(x, y));
				}
				
				pIC->Set_Value(x, y, log10(up / pDown->asDouble(x, y)));
			}
		}
	}

	pDown->Set_Name(_TL("Downslope Component"));

	Set_Progress(100., 100.);

	return( true );
}


//---------------------------------------------------------
bool CIndex_of_Connectivity::Get_Stats(CSG_Grid *pGrid, CSG_Grid_Cell_Addressor *pKernel, int x, int y, double &Value, bool bStdDev)
{
	CSG_Simple_Statistics s;

	if( pGrid->is_InGrid(x, y) )
	{
		for(int i=0; i<pKernel->Get_Count(); i++)
		{
			int ix = pKernel->Get_X(i, x);
			int iy = pKernel->Get_Y(i, y);

			if( pGrid->is_InGrid(ix, iy) )
			{
				s += pGrid->asDouble(ix, iy);
			}
		}
	}

	if( s.Get_Count() > 0 )
	{
		if( bStdDev )
		{
			Value = s.Get_StdDev();
		}
		else
		{
			Value = s.Get_Mean();
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CIndex_of_Connectivity::Get_Slope(CSG_Grid *pGrid, int x, int y, double& Value)
{
	if (!pGrid->is_InGrid(x, y))
	{
		return( false );
	}

	double Z[8];
	double z = pGrid->asDouble(x, y);
	Value = 0.;

	for(int i=0; i<8; i++)
	{
		int ix = Get_xTo(i, x);
		int iy = Get_yTo(i, y);

		if( !pGrid->is_InGrid(ix, iy) )
		{
			Z[i] = 0.;
		}
		else
		{
			Z[i] = (z - pGrid->asDouble(ix, iy)) / Get_Length(i);

			if (Z[i] > Value)
			{
				Value = Z[i];
			}
		}
	}

	return( true );
}


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

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