📄 pit_router.cpp
字号:
///////////////////////////////////////////////////////////
// //
// SAGA //
// //
// System for Automated Geoscientific Analyses //
// //
// Module Library: //
// ta_preprocessor //
// //
//-------------------------------------------------------//
// //
// Pit_Router.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: oconrad@saga-gis.org //
// //
// contact: Olaf Conrad //
// Institute of Geography //
// University of Goettingen //
// Goldschmidtstr. 5 //
// 37077 Goettingen //
// Germany //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#include <memory.h>
#include "Pit_Router.h"
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#define IS_Flat(a,b) (a==b)
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
CPit_Router::CPit_Router(void)
{
Set_Name(_TL("Sink Drainage Route Detection"));
Set_Author(_TL("Copyrights (c) 2001 by Olaf Conrad"));
Set_Description(
_TL("")
);
Parameters.Add_Grid(
NULL , "ELEVATION" , _TL("Elevation"),
_TL(""),
PARAMETER_INPUT
);
Parameters.Add_Grid(
NULL , "SINKROUTE" , _TL("Sink Route"),
_TL(""),
PARAMETER_OUTPUT
);
Parameters.Add_Value(
NULL , "THRESHOLD" , _TL("Threshold"),
_TL(""),
PARAMETER_TYPE_Bool
);
Parameters.Add_Value(
NULL , "THRSHEIGHT" , _TL("Threshold Height"),
_TL(""),
PARAMETER_TYPE_Double , 100
);
}
//---------------------------------------------------------
CPit_Router::~CPit_Router(void)
{}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CPit_Router::On_Execute(void)
{
double Threshold;
CSG_Grid *pDEM, *pRoute;
pDEM = Parameters("ELEVATION")->asGrid();
pRoute = Parameters("SINKROUTE")->asGrid();
Threshold = Parameters("THRESHOLD")->asBool()
? Parameters("THRSHEIGHT")->asDouble() : -1.0;
Get_Routes(pDEM, pRoute, Threshold);
return( true );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
int CPit_Router::Get_Routes(CSG_Grid *pDEM, CSG_Grid *pRoute, double Threshold)
{
int iPit, nPits, n;
TPit_Outlet *pOutlet, *pNext;
//-----------------------------------------------------
m_pDEM = pDEM;
m_pRoute = pRoute;
m_Threshold = Threshold > 0.0 ? Threshold : -1.0;
//-----------------------------------------------------
m_pPits = NULL;
m_Pit = NULL;
m_pFlats = NULL;
m_Flat = NULL;
m_Outlets = NULL;
//-----------------------------------------------------
m_System.Assign(m_pDEM->Get_System());
//-----------------------------------------------------
if( Initialize() )
{
//-------------------------------------------------
// 1. Pits/Flats finden...
SG_UI_Process_Set_Text(_TL("Find Pits"));
nPits = Find_Pits();
if( nPits > 0 )
{
//---------------------------------------------
// 2. Pit/Flat-Zugehoerigkeiten u. pot. m_Outlets finden...
SG_UI_Process_Set_Text(_TL("Find Outlets"));
Find_Outlets(nPits);
//---------------------------------------------
// 3. Routing vornehmen...
SG_UI_Process_Set_Text(_TL("Routing"));
iPit = 0;
do
{
pOutlet = m_Outlets;
while( pOutlet && SG_UI_Process_Get_Okay(false) )
{
pNext = pOutlet->Next;
n = Find_Route(pOutlet);
if( n > 0 )
{
pOutlet = m_Outlets;
iPit += n;
SG_UI_Process_Set_Progress(iPit, nPits);
}
else
{
pOutlet = pNext;
}
}
if( iPit < nPits ) // Thresholding may have prevented total removal of pits...
{
for(n=0; n<nPits; n++)
{
if( !m_Pit[n].bDrained )
{
m_Pit[n].bDrained = true;
iPit++;
break;
}
}
}
}
while( iPit < nPits && SG_UI_Process_Set_Progress(iPit, nPits) );
}
}
//-----------------------------------------------------
Process_Set_Text(_TL("Finalize"));
Finalize();
if( Process_Get_Okay(false) )
{
if( nPits > 0 )
{
Message_Add(CSG_String::Format(_TL("%d sinks have been processed."), nPits));
return( nPits );
}
else
{
Message_Add(_TL("No sinks have been detected."));
}
}
return( 0 );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CPit_Router::Initialize(void)
{
if( m_pDEM && m_pDEM->is_Valid()
&& m_pRoute && m_pRoute->is_Valid()
&& m_pDEM->Get_System() == m_pRoute->Get_System() )
// && m_pDEM->is_Compatible(m_pRoute->Get_System()) )
{
m_pRoute->Assign();
m_pPits = SG_Create_Grid(m_pDEM, GRID_TYPE_Int);
m_pPits->Assign();
m_Pit = NULL;
m_pFlats = NULL;
m_Flat = NULL;
m_Outlets = NULL;
return( true );
}
return( false );
}
//---------------------------------------------------------
void CPit_Router::Finalize(void)
{
TPit_Outlet *pOutlet;
if( m_pPits )
{
delete( m_pPits );
m_pPits = NULL;
}
if( m_Pit )
{
SG_Free(m_Pit);
m_Pit = NULL;
}
if( m_pFlats )
{
delete( m_pFlats );
m_pFlats = NULL;
}
if( m_Flat )
{
SG_Free(m_Flat);
m_Flat = NULL;
}
while( m_Outlets )
{
pOutlet = m_Outlets->Next;
SG_Free(m_Outlets);
m_Outlets = pOutlet;
}
m_Outlets = NULL;
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
int CPit_Router::Find_Pits(void)
{
bool bLower, bFlat;
int x, y, i, ix, iy, nFlats, nPits;
long n;
double z;
TPit *pPit;
//-----------------------------------------------------
nFlats = 0;
nPits = 0;
for(n=0; n<m_System.Get_NCells() && SG_UI_Process_Set_Progress(n, m_System.Get_NCells()); n++)
{
m_pDEM->Get_Sorted(n,x,y,false); // von tief nach hoch...
if( x > 0 && x < m_System.Get_NX() - 1 && y > 0 && y < m_System.Get_NY() - 1 // Randzellen und Missing Values sind
&& !m_pDEM->is_NoData(x, y) // per Definition drainiert (:= 0)...
&& m_pPits->asInt(x, y) == 0 ) // ...oder schon als m_Flat markiert sein...
{
z = m_pDEM->asDouble(x,y);
bLower = false;
bFlat = false;
for(i=0; i<8 && !bLower; i++)
{
ix = m_System.Get_xTo(i,x);
iy = m_System.Get_yTo(i,y);
if( !m_pDEM->is_InGrid(ix, iy) || z > m_pDEM->asDouble(ix, iy) )
{
bLower = true;
}
else if( IS_Flat(z, m_pDEM->asDouble(ix,iy)) )
{
bFlat = true;
}
}
//-----------------------------------------
if( !bLower ) // Pit or Flat...
{
nPits++;
m_pPits->Set_Value(x,y, nPits);
m_Pit = (TPit *)SG_Realloc(m_Pit, nPits * sizeof(TPit));
pPit = m_Pit + nPits - 1;
pPit->bDrained = false;
pPit->z = z;
if( bFlat )
{
nFlats++;
m_Flat = (TGEO_iRect *)SG_Realloc(m_Flat, nFlats * sizeof(TGEO_iRect));
Mark_Flat(x, y, m_Flat + nFlats - 1, nFlats, nPits);
}
}
}
}
return( nPits );
}
//---------------------------------------------------------
int CPit_Router::Find_Outlets(int nPits)
{
bool bOutlet, bExArea, bGoExArea;
int x, y, i, ix, iy, iMin,
iID, j, jID, Pit_ID[8];
long n;
double z, dz, dzMin;
TPit_Outlet *pOutlet;
//-----------------------------------------------------
if( nPits > 0 && SG_UI_Process_Get_Okay(false) )
{
pOutlet = NULL;
m_nJunctions = (int *)SG_Calloc(nPits, sizeof(int ));
m_Junction = (int **)SG_Calloc(nPits, sizeof(int *));
//-------------------------------------------------
for(n=0; n<m_System.Get_NCells() && SG_UI_Process_Set_Progress(n, m_System.Get_NCells()); n++)
{
m_pDEM->Get_Sorted(n, x, y, false); // von tief nach hoch...
if( m_pPits->asInt(x,y) == 0 )
{
z = m_pDEM->asDouble(x,y);
iMin = -1;
bOutlet = false;
bGoExArea = false;
//-----------------------------------------
for(i=0; i<8; i++)
{
ix = m_System.Get_xTo(i,x);
iy = m_System.Get_yTo(i,y);
bExArea = !m_pDEM->is_InGrid(ix, iy);
if( bExArea || z > m_pDEM->asDouble(ix,iy) )
{
Pit_ID[i] = iID = bExArea ? 0 : m_pPits->asInt(ix,iy);
if( iID >= 0 )
{
for(j=0; j<i && !bOutlet; j++)
{
jID = Pit_ID[j];
if( jID >= 0 && !Get_Junction(iID, jID) )
{
bOutlet = true;
}
}
}
//---------------------------------
if( !bGoExArea )
{
if( bExArea )
{
bGoExArea = true;
iMin = i;
}
else
{
dz = (z - m_pDEM->asDouble(ix,iy)) / m_System.Get_Length(i);
if( iMin < 0 || dzMin < dz )
{
iMin = i;
dzMin = dz;
}
}
}
}
else
{
Pit_ID[i] = -1;
}
}
//-----------------------------------------
if( bOutlet )
{
if( pOutlet )
{
pOutlet->Next = (TPit_Outlet *)SG_Malloc(sizeof(TPit_Outlet));
pOutlet->Next->Prev = pOutlet;
pOutlet = pOutlet->Next;
}
else
{
m_Outlets = pOutlet = (TPit_Outlet *)SG_Malloc(sizeof(TPit_Outlet));
m_Outlets->Prev = NULL;
}
pOutlet->Next = NULL;
pOutlet->x = x;
pOutlet->y = y;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -