📄 saga_wetness_index.cpp
字号:
///////////////////////////////////////////////////////////
// //
// SAGA //
// //
// System for Automated Geoscientific Analyses //
// //
// Module Library: //
// ta_hydrology //
// //
//-------------------------------------------------------//
// //
// SAGA_Wetness_Index.cpp //
// //
// Copyright (C) 2005 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. //
// //
//-------------------------------------------------------//
// //
// contact: Olaf Conrad //
// Institute of Geography //
// University of Goettingen //
// Goldschmidtstr. 5 //
// 37077 Goettingen //
// Germany //
// //
// e-mail: oconrad@saga-gis.org //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#include "SAGA_Wetness_Index.h"
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
CSAGA_Wetness_Index::CSAGA_Wetness_Index(void)
{
//-----------------------------------------------------
Set_Name (_TL("SAGA Wetness Index"));
Set_Author (_TL("Copyrights (c) 2005 by Juergen Boehner, Olaf Conrad"));
Set_Description (_TW(
"The 'SAGA Wetness Index' is, as the name says, similar to the "
"'Topographic Wetness Index' (TWI), but it is based on a modified "
"catchment area calculation ('Modified Catchment Area'), which does "
"not think of the flow as very thin film. As result it predicts for "
"cells situated in valley floors with a small vertical distance to "
"a channel a more realistic, higher potential soil moisture compared "
"to the standard TWI calculation.\n\n"
"References\n"
"- Boehner, J., Koethe, R. Conrad, O., Gross, J., Ringeler, A., Selige, T. (2002): "
"Soil Regionalisation by Means of Terrain Analysis and Process Parameterisation. "
"In: Micheli, E., Nachtergaele, F., Montanarella, L. [Ed.]: Soil Classification 2001. "
"European Soil Bureau, Research Report No. 7, EUR 20398 EN, Luxembourg. pp.213-222."
));
//-----------------------------------------------------
Parameters.Add_Grid(
NULL , "DEM" , _TL("Elevation"),
_TL(""),
PARAMETER_INPUT
);
Parameters.Add_Grid(
NULL , "C" , _TL("Catchment area"),
_TL(""),
PARAMETER_OUTPUT
);
Parameters.Add_Grid(
NULL , "GN" , _TL("Catchment slope"),
_TL(""),
PARAMETER_OUTPUT
);
Parameters.Add_Grid(
NULL , "CS" , _TL("Modified catchment area"),
_TL(""),
PARAMETER_OUTPUT
);
Parameters.Add_Grid(
NULL , "SB" , _TL("Wetness index"),
_TL(""),
PARAMETER_OUTPUT
);
Parameters.Add_Value(
NULL , "T" , _TL("t"),
_TL(""),
PARAMETER_TYPE_Double , 10.0, 0.0, true
);
}
//---------------------------------------------------------
CSAGA_Wetness_Index::~CSAGA_Wetness_Index(void)
{}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CSAGA_Wetness_Index::On_Execute(void)
{
double a, t;
CSG_Grid *pDEM, *pC, *pGN, *pCS, *pSB;
//-----------------------------------------------------
pDEM = Parameters("DEM") ->asGrid();
pC = Parameters("C") ->asGrid();
pGN = Parameters("GN") ->asGrid();
pCS = Parameters("CS") ->asGrid();
pSB = Parameters("SB") ->asGrid();
t = Parameters("T") ->asDouble();
a = 0.0174532;
//-----------------------------------------------------
Get_Area_Catchment (pDEM, pC, pGN);
Get_Area_Modified (pDEM, pC, pCS, t);
Get_Wetness_Index (pDEM, pCS, pSB, a);
//-----------------------------------------------------
return( true );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
// Im Unterschied zu Freeman's urspruenglichen Verfahren
// wird der Winkel (atan(dz / dx)) und nicht das Gefaelle
// (dz / dx) fuer die Gewichtung der Abfluszanteile benutzt!
//---------------------------------------------------------
bool CSAGA_Wetness_Index::Get_Area_Catchment(CSG_Grid *pDEM, CSG_Grid *pC, CSG_Grid *pS)
{
const double MFD_Converge = 1.1;
int n, x, y, i, ix, iy;
double z, d, dz[8], dzSum, c, s;
//-----------------------------------------------------
Process_Set_Text(_TL("Catchment area calculation..."));
pC->Assign(1.0);
pS->Assign(0.0);
//-----------------------------------------------------
for(n=0; n<Get_NCells() && Set_Progress_NCells(n); n++)
{
pDEM->Get_Sorted(n, x, y);
if( pDEM->is_NoData(x, y) )
{
pC->Set_NoData(x, y);
}
else
{
pDEM->Get_Gradient(x, y, s, c);
c = pC->asDouble(x, y);
s = pS->asDouble(x, y) + s;
pC->Set_Value(x, y, c);
pS->Set_Value(x, y, s / c);
for(i=0, dzSum=0.0, z=pDEM->asDouble(x, y); i<8; i++)
{
ix = Get_xTo(i, x);
iy = Get_yTo(i, y);
if( is_InGrid(ix, iy) && !pDEM->is_NoData(ix, iy) && (d = z - pDEM->asDouble(ix, iy)) > 0.0 )
{
dzSum += (dz[i] = pow(atan(d / Get_Length(i)), MFD_Converge));
}
else
{
dz[i] = 0.0;
}
}
if( dzSum > 0.0 )
{
for(i=0; i<8; i++)
{
if( dz[i] > 0.0 )
{
ix = Get_xTo(i, x);
iy = Get_yTo(i, y);
d = dz[i] / dzSum;
pC->Add_Value(ix, iy, d * c);
pS->Add_Value(ix, iy, d * s);
}
}
}
}
}
//-----------------------------------------------------
*pC *= pC->Get_Cellarea();
return( true );
}
/*/---------------------------------------------------------
bool CSAGA_Wetness_Index::Get_Area_Catchment(CSG_Grid *pDEM, CSG_Grid *pC, CSG_Grid *pGN)
{
int n, x, y, i, ix, iy;
double z, d, dz[8], dzSum, c, w, sw, sn, Slope, Azimuth,
Weight = 0.5,
a = Get_Cellsize() * Get_Cellsize(),
MFD_Converge = 1.1;
CSG_Grid W, SW, SN;
//-----------------------------------------------------
Process_Set_Text("Catchment area/slope calculation...");
pC ->Assign(0.0);
pGN->Assign(0.0);
W .Create(pDEM);
SW.Create(pDEM);
SN.Create(pDEM);
//-----------------------------------------------------
for(n=0; n<Get_NCells() && Set_Progress_NCells(n); n++)
{
pDEM->Get_Sorted(n, x, y);
if( pDEM->is_NoData(x, y) )
{
pC ->Set_NoData(x, y);
pGN->Set_NoData(x, y);
}
else
{
pDEM->Get_Gradient(x, y, Slope, Azimuth);
z = pDEM->asDouble(x, y);
c = pC ->asDouble(x, y) + a;
w = W .asDouble(x, y) + pow(c, Weight); //{ W[p] = C[p]^w;}
sw = SW .asDouble(x, y) + w; //{SW[p] = W[p];}
sn = SN .asDouble(x, y) + w * Slope; //{SN[p] = W[p] * N[p];}
for(i=0, dzSum=0.0; i<8; i++)
{
ix = Get_xTo(i, x);
iy = Get_yTo(i, y);
if( is_InGrid(ix, iy) && !pDEM->is_NoData(ix, iy) && (d = z - pDEM->asDouble(ix, iy)) > 0.0 )
{
dzSum += (dz[i] = pow(atan(d / Get_Length(i)), MFD_Converge));
}
else
{
dz[i] = 0.0;
}
}
if( dzSum > 0.0 )
{
for(i=0; i<8; i++)
{
if( dz[i] > 0.0 )
{
ix = Get_xTo(i, x);
iy = Get_yTo(i, y);
d = dz[i] / dzSum;
pC->Add_Value(ix, iy, d * c);
W .Add_Value(ix, iy, d * w); //{ W[p] = W[p] + W[p+pul] * UL[p] + W[p+pl] * LL[p] + W[p+pol] * OL[p] + W[p+po] * OO[p] + W[p+por] * OR[p] + W[p+pr] * RR[p] + W[p+pur] * UR[p] + W[p+pu] * UU[p]; gefunden = 1;}
SW .Add_Value(ix, iy, d * sw); //{SW[p] = SW[p] + SW[p+pul] * UL[p] + SW[p+pl] * LL[p] + SW[p+pol] * OL[p] + SW[p+po] * OO[p] + SW[p+por] * OR[p] + SW[p+pr] * RR[p] + SW[p+pur] * UR[p] + SW[p+pu] * UU[p]; gefunden = 1;}
SN .Add_Value(ix, iy, d * sn); //{SN[p] = SN[p] + SN[p+pul] * UL[p] + SN[p+pl] * LL[p] + SN[p+pol] * OL[p] + SN[p+po] * OO[p] + SN[p+por] * OR[p] + SN[p+pr] * RR[p] + SN[p+pur] * UR[p] + SN[p+pu] * UU[p]; gefunden = 1;}
}
}
}
}
pC ->Set_Value(x, y, c);
pGN->Set_Value(x, y, sn / sw); //{GN[p] = SN[p] / SW[p];}
}
//-----------------------------------------------------
return( true );
}/**/
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
// in den folgenden drei Schritten wird die gesaugte
// Einzugsgebietsgroe遝 CS ermittelt, wobei der t-Parameter
// die Staerke der Saugung steuert. Werte unter 10 (z.B. 5)
// fuehren zu einer starken Saugung, Werte ueber 10 (z.B. 15)
// zu einer schwachen Saugung. Die gesaugte
// Einzugsgebietsgroe遝 selbst stellt bereits einen Parameter
// dar, der die raeumliche Relief-bedingte Feuchteverteilung
// in guter Weise annaehert
//---------------------------------------------------------
bool CSAGA_Wetness_Index::Get_Area_Modified(CSG_Grid *pDEM, CSG_Grid *pC, CSG_Grid *pCS, double t)
{
bool bRecalculate;
int x, y, i, ix, iy, n;
double z, d;
CSG_Grid C, C_Last, T;
//-----------------------------------------------------
Process_Set_Text(_TL("Modify: pre-processing..."));
T.Create(pC);
for(y=0; y<Get_NY() && Set_Progress(y); y++)
{
for(x=0; x<Get_NX(); x++)
{
if( !pDEM->is_NoData(x, y) && !pC->is_NoData(x, y) )
{
pDEM->Get_Gradient(x, y, z, d);
z = pow(t, z);
z = pow(1.0 / z, exp(z));
T.Set_Value(x, y, z);
}
else
{
T.Set_NoData(x, y);
}
}
}
C .Create(pC);
C .Assign(pC);
C_Last.Create(pC);
C_Last.Assign(pC);
//-----------------------------------------------------
for(i=0, n=1; n>0; i++)
{
for(y=0, n=0; y<Get_NY() && Process_Get_Okay(false); y++)
{
for(x=0; x<Get_NX(); x++)
{
if( !T.is_NoData(x, y) )
{
z = T.asDouble(x, y) * Get_Local_Maximum(&C, x, y);
if( z > C.asDouble(x, y) )
{
n++;
C.Set_Value(x, y, z);
}
}
}
}
if( n > 0 )
{
for(y=0, n=0; y<Get_NY() && Process_Get_Okay(false); y++)
{
for(x=0; x<Get_NX(); x++)
{
if( C.asDouble(x, y) != C_Last.asDouble(x, y) )
{
n++;
C_Last.Set_Value(x, y, C.asDouble(x, y));
}
}
}
}
Process_Set_Text(CSG_String::Format(SG_T("Modify: %d. iteration [%d > 0]"), i, n));
}
//-----------------------------------------------------
Process_Set_Text(_TL("Modify: post-processing..."));
for(y=0; y<Get_NY() && Set_Progress(y); y++)
{
for(x=0; x<Get_NX(); x++)
{
if( pC->is_NoData(x, y) )
{
pCS->Set_NoData(x, y);
}
else
{
for(iy=y-1, bRecalculate=false; iy<=y+1 && !bRecalculate; iy++)
{
for(ix=x-1; ix<=x+1 && !bRecalculate; ix++)
{
if( C.is_InGrid(ix, iy) && pC->is_InGrid(ix, iy) && C.asDouble(ix, iy) > pC->asDouble(ix, iy) )
{
bRecalculate = true;
}
}
}
if( bRecalculate )
{
for(iy=y-1, z=0.0, n=0; iy<=y+1; iy++)
{
for(ix=x-1; ix<=x+1; ix++)
{
if( C.is_InGrid(ix, iy) )
{
n++;
z += C.asDouble(ix, iy);
}
}
}
z = z / (double)n;
}
else
{
z = C.asDouble(x, y);
}
pCS->Set_Value(x, y, z);
}
}
}
//-----------------------------------------------------
return( true );
}
//---------------------------------------------------------
double CSAGA_Wetness_Index::Get_Local_Maximum(CSG_Grid *pGrid, int x, int y)
{
int i, ix, iy;
double z = pGrid->asDouble(x, y);
for(i=0; i<8; i++)
{
ix = Get_xTo(i, x);
iy = Get_yTo(i, y);
if( pGrid->is_InGrid(ix, iy) && pGrid->asDouble(ix, iy) > z )
{
z = pGrid->asDouble(ix, iy);
}
}
return( z );
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -