flow_recursivedown.cpp

这是一个GPS相关的程序

///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                    Module Library:                    //
//                     ta_hydrology                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                Flow_RecursiveDown.cpp                 //
//                                                       //
//                 Copyright (C) 2003 by                 //
//                      Olaf Conrad                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// 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.                                                  //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     [e-mail]                   //
//                                                       //
//    contact:    Olaf Conrad                            //
//                Institute of Geography                 //
//                University of Goettingen               //
//                Goldschmidtstr. 5                      //
//                37077 Goettingen                       //
//                Germany                                //
//                                                       //
///////////////////////////////////////////////////////////

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


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

//---------------------------------------------------------
#include "Flow_RecursiveDown.h"

//---------------------------------------------------------
#define GET_OUTLET_DIAG__1(in, angle)		(1.0 - (1.0 - in) * tan(M_PI_090 - angle))
#define GET_OUTLET_CROSS_1(in, angle)		(in + tan(angle))

#define GET_OUTLET_DIAG__2(in, angle)		(in * tan(angle))
#define GET_OUTLET_CROSS_2(in, angle)		(in - tan(M_PI_090 - angle))

#define GET_LENGTH(a, b)					(sqrt((a)*(a) + (b)*(b)))


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

//---------------------------------------------------------
CFlow_RecursiveDown::CFlow_RecursiveDown(void)
{
	Set_Name(_TL("Flow Tracing"));

	Set_Author(_TL("Copyrights (c) 2001 by Olaf Conrad"));

	Set_Description	(_TW(
		"Flow tracing algorithms for calculations of flow accumulation and related parameters. "
		"These algorithms trace the flow of each cell in a DEM separately until it finally leaves the DEM or ends in a sink.\n\n"

		"References:\n\n"

		"Rho 8 (this implementation adopted the original algorithm only for the flow routing and will give quite different results):\n"
		"- Fairfield, J. / Leymarie, P. (1991):\n"
		"    'Drainage networks from grid digital elevation models',\n"
		"    Water Resources Research, 27:709-717\n\n"

		"Kinematic Routing Algorithm:\n"
		"- Lea, N.L. (1992):\n"
		"    'An aspect driven kinematic routing algorithm',\n"
		"    in: Parsons, A.J., Abrahams, A.D. (Eds.), 'Overland Flow: hydraulics and erosion mechanics', London, 147-175\n\n"

		"DEMON:\n"
		"- Costa-Cabral, M. / Burges, S.J. (1994):\n"
		"    'Digital Elevation Model Networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas',\n"
		"    Water Resources Research, 30:1681-1692\n\n")
	);


	//-----------------------------------------------------
	// Method...

	Parameters.Add_Choice(
		NULL	, "Method"		, _TL("Method"),
		_TL(""),
		CSG_String::Format(SG_T("%s|%s|%s|"),
			_TL("Rho 8"),
			_TL("Kinematic Routing Algorithm"),
			_TL("DEMON")
		), 1
	);


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

	Parameters.Add_Value(
		NULL	, "MINDQV"		, _TL("DEMON - Min. DQV"),
		_TL("DEMON - Minium Drainage Quota Volume (DQV) for traced flow tubes"),
		PARAMETER_TYPE_Double	,	0.0, 0.0, true, 1.0, true
	);

	Parameters.Add_Value(
		NULL	, "CORRECT"		, _TL("Flow Correction"),
		_TL(""),
		PARAMETER_TYPE_Bool
	);
}

//---------------------------------------------------------
CFlow_RecursiveDown::~CFlow_RecursiveDown(void)
{}


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

//---------------------------------------------------------
void CFlow_RecursiveDown::On_Initialize(void)
{
	int		x, y;
	double	Slope, Aspect;

	//-----------------------------------------------------
	Method			= Parameters("Method"	)->asInt();
	DEMON_minDQV	= Parameters("MINDQV"	)->asDouble();
	bFlowPathWeight	= Parameters("CORRECT"	)->asBool();

	pLinear			= 1 ? SG_Create_Grid(pDTM, GRID_TYPE_Float) : NULL;

	//-----------------------------------------------------
	Lock_Create();

	//-----------------------------------------------------
	switch( Method )
	{
	case 0: default:	// Rho 8...

		pDir		= NULL;
		pDif		= NULL;

		break;

	//-----------------------------------------------------
	case 1:	case 2:		// KRA, DEMON...

		pDir		= SG_Create_Grid(pDTM, GRID_TYPE_Char);
		pDif		= SG_Create_Grid(pDTM, GRID_TYPE_Float);

		for(y=0; y<Get_NY() && Set_Progress(y); y++)
		{
			for(x=0; x<Get_NX(); x++)
			{
				if( !pDTM->is_NoData(x, y) )
				{
					Get_Gradient(x, y, Slope, Aspect);

					if( Aspect >= 0.0 )
					{
						pDir->Set_Value(x, y, 2 * (((int)(Aspect / M_PI_090)) % 4));
						pDif->Set_Value(x, y, fmod(Aspect, M_PI_090));
					}
				}
			}
		}

		break;
	}
}

//---------------------------------------------------------
void CFlow_RecursiveDown::On_Finalize(void)
{
	int		x, y, dir;
	long	n;
	double	qFlow;

	//-----------------------------------------------------
	if( pDir )
	{
		delete(pDir);
	}

	if( pDif )
	{
		delete(pDif);
	}

	Lock_Destroy();

	//-----------------------------------------------------
	if( pLinear )
	{
		for(n=0; n<Get_NCells() && Set_Progress_NCells(n); n++)
		{
			pDTM->Get_Sorted(n, x, y);

			if( !pDTM->is_NoData(x, y) && (qFlow = pLinear->asDouble(x, y)) > 0.0 )
			{
				Add_Flow(x, y, qFlow);

				if( (dir = pDTM->Get_Gradient_NeighborDir(x, y)) >= 0 )
				{
					x	= Get_xTo(dir, x);
					y	= Get_yTo(dir, y);

					if( pDTM->is_InGrid(x, y) )
					{
						pLinear->Add_Value(x, y, qFlow);
					}
				}
			}
		}

		delete(pLinear);

		pLinear	= NULL;
	}
}


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

//---------------------------------------------------------
void CFlow_RecursiveDown::Calculate(void)
{
	int		x, y;

	for(y=0; y<Get_NY() && Set_Progress(y); y+=Step)
	{
		//if( !(y%2) )DataObject_Update(pFlow);

		for(x=0; x<Get_NX(); x+=Step)
		{
			Calculate(x, y);
		}
	}
}

//---------------------------------------------------------
void CFlow_RecursiveDown::Calculate(int x, int y)
{
	double 	Slope, Aspect, qFlow;

	if( !pDTM->is_NoData(x, y) && (qFlow = pWeight ? pWeight->asDouble(x, y) : 1.0) > 0.0 )
	{
		Get_Gradient(x, y, Slope, Aspect);

		Src_Height	= pDTM->asDouble(x,y);
		Src_Slope	= Slope;

		Add_Flow(x, y, qFlow);
		Lock_Set(x, y, 1);

		//-------------------------------------------------
		switch( Method )
		{
		case 0:
			Rho8_Start(		x, y, qFlow );
			break;

		case 1:
			KRA_Start(		x, y, qFlow );
			break;

		case 2:
			DEMON_Start(	x, y, qFlow );
			break;
		}

		Lock_Set(x, y, 0);
	}
}


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

//---------------------------------------------------------
void CFlow_RecursiveDown::Add_Flow(int x, int y, double Fraction)
{
	if( pCatch )
	{
		pCatch			->Add_Value(x, y, Fraction );
	}

	if( pCatch_Height )
	{
		pCatch_Height	->Add_Value(x, y, Fraction * Src_Height );
	}

	if( pCatch_Slope )
	{
		pCatch_Slope	->Add_Value(x, y, Fraction * Src_Slope );
	}
}


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

//---------------------------------------------------------
void CFlow_RecursiveDown::Rho8_Start(int x, int y, double qFlow)
{
	int		dir, ix, iy;
	double 	Slope, Aspect;

	//-----------------------------------------------------
	Get_Gradient(x, y, Slope, Aspect);

	if( Aspect >= 0.0 )
	{
		dir		= (int)(Aspect / M_PI_045);

		if( fmod(Aspect, M_PI_045) / M_PI_045 > rand() / (double)RAND_MAX )
		{
			dir++;
		}

		dir	%= 8;

		ix		= Get_xTo(dir, x);
		iy		= Get_yTo(dir, y);

		//-------------------------------------------------
		if( is_InGrid(ix, iy) )