watersheds_ext.cpp

这是一个GPS相关的程序

/*******************************************************************************
    Watersheds.cpp
    Copyright (C) Victor Olaya
    
    This program 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.

    This program 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 "Watersheds_ext.h"
#include "Helper.h"
#include <string>

#define NO_EQUIVALENT_RECTANGLE -1

CWatersheds_ext::CWatersheds_ext(void){

	Parameters.Set_Name(_TL("Watershed Basins (extended)"));
	Parameters.Set_Description(_TW("(c) 2004 by Victor Olaya. Subbasin" ));

	Parameters.Add_Grid(NULL, 
						"DEM", 
						_TL("DEM"), 
						_TL(""), 
						PARAMETER_INPUT);

	Parameters.Add_Grid(NULL, 
						"CHANNELS", 
						_TL("Drainage Network"), 
						_TL(""), 
						PARAMETER_INPUT);

	Parameters.Add_Grid(NULL, 
						"SOILLOSS", 
						_TL("Soil Losses"), 
						_TL(""), 
						PARAMETER_INPUT_OPTIONAL);

	Parameters.Add_Grid(NULL, 
						"CN", 
						_TL("Curve Number"), 
						_TL(""), 
						PARAMETER_INPUT_OPTIONAL);

	Parameters.Add_Grid(NULL, 
						"BASINGRID", 
						_TL("Subbasins"), 
						_TL(""), 
						PARAMETER_OUTPUT);

	Parameters.Add_Shapes(NULL, 
						"BASINS", 
						_TL("Subbasins"), 
						_TL(""), 
						PARAMETER_OUTPUT);

	Parameters.Add_Shapes(NULL, 
						"HEADERS", 
						_TL("River Headers"), 
						_TL(""), 
						PARAMETER_OUTPUT);
	
	Parameters.Add_Choice(NULL, 
						"FRAGMENTATION", 
						_TL("Basin subdivision"), 
						_TL(""), 
						CSG_String::Format(SG_T("%s|%s|"),
							_TL("Only closing points on main stream"),
							_TL("All")
						),0
					);

}//constructor

CWatersheds_ext::~CWatersheds_ext(void){}

bool CWatersheds_ext::On_Execute(void){
	
	m_pDEM = Parameters("DEM")->asGrid(); 
	m_pCNGrid = Parameters("CN")->asGrid();
	m_pBasinGrid = Parameters("BASINGRID")->asGrid();
	m_pSoilLossGrid = Parameters("SOILLOSS")->asGrid();
	m_pChannelsGrid = Parameters("CHANNELS")->asGrid();
	m_pBasins = Parameters("BASINS")->asShapes();
	m_pHeaders = Parameters("HEADERS")->asShapes();
	m_iFragmentationType = Parameters("FRAGMENTATION")->asInt();

	m_pBasinGrid->Assign((double)0);
	m_pBasinGrid->Set_Name(_TL("Subbasins"));
	m_pBasinGrid->Set_Description(_TW("Subbasins"));

	CalculateBasin();
	CreateShapesLayer();

	delete[] m_fMaxDistance;
	delete[] m_fHeightDif;
	delete m_pDistanceGrid;


	return( true );

}//method


void CWatersheds_ext::CalculateBasin() {

	int i,j;
	float fMinHeight = 1000000000;
	int iNextX, iNextY;
	int iX, iY;
    bool bAddPoint;
	TSG_Point		P;
	CSG_Points	Crossings;  

	Process_Set_Text(_TL("Calculating Subbasins..."));
	
	for(int y=0; y<Get_NY() && Set_Progress(y); y++){		
		for(int x=0; x<Get_NX(); x++){
			if(!m_pChannelsGrid->is_NoData(x,y)){
				if (isHeader(x,y)){
					P.x=x;
					P.y=y;
					m_Headers.Add(P);				
				}//if
				if (m_pDEM->asFloat(x,y)<fMinHeight){
					m_iClosingX = x;
					m_iClosingY = y;
					fMinHeight = m_pDEM->asFloat(x,y);
				}//if
			}//if
		}// for
	}// for
				
	for (i = 0; i < m_Headers.Get_Count(); i++) {
	    iNextX = m_Headers.Get_Point(i).x;
		iNextY = m_Headers.Get_Point(i).y;
	    do {
	        iX = iNextX;
	        iY = iNextY;
			getNextCell(m_pDEM, m_pChannelsGrid, iX, iY, iNextX, iNextY);	                
	        if (m_pChannelsGrid->asInt(iNextX,iNextY) != m_pChannelsGrid->asInt(iX,iY)) {
	            bAddPoint = true;
	            for (j = 0; j < Crossings.Get_Count(); j++) {
	                if (Crossings.Get_Point(j).x == iX && Crossings.Get_Point(j).y == iY){
	                    bAddPoint = false;
	                }//if
	            }// for
	            if (bAddPoint) {							
					P.x = iX;
					P.y = iY;
	                Crossings.Add(P);
	            }// if	
	        }// if
	    }while (!(iX == m_iClosingX && iY == m_iClosingY)
	            && (iX != iNextX || iY != iNextY));
	}// for
	
	float fMaxDistance = -1;
	float fDistance = 0;
	TSG_Point	LongestRiverHeader;
	
	for (i = 0; i < m_Headers.Get_Count(); i++) {
	    fDistance = DistanceToClosingPoint(m_Headers.Get_Point(i).x, m_Headers.Get_Point(i).y);
	    if (fDistance > fMaxDistance) {
	        fMaxDistance = fDistance;
	        LongestRiverHeader.x = m_Headers.Get_Point(i).x; 
			LongestRiverHeader.y = m_Headers.Get_Point(i).y; 
	    }// if
	}// for		

	float fHeightDif = m_pDEM->asFloat((int)LongestRiverHeader.x, (int)LongestRiverHeader.y);
		
	if (m_iFragmentationType == FRAGMENTATION_TO_MAIN) {

		CSG_Points	CrossingsNext;
		CSG_Points	CrossingsAtPrincipal;
		
	    for (i = 0; i < Crossings.Get_Count(); i++) {
	        getNextCell(m_pDEM, m_pChannelsGrid, Crossings.Get_Point(i).x,Crossings.Get_Point(i).y, iNextX,iNextY);
			P.x = iNextX;
			P.y = iNextY;
	        CrossingsNext.Add(P);
	    }// for

		iNextX = LongestRiverHeader.x;
		iNextY = LongestRiverHeader.y;
		
	    do {
	        iX = iNextX;
	        iY = iNextY;
	        getNextCell(m_pDEM, m_pChannelsGrid, iX, iY, iNextX, iNextY);
	        for (i = 0; i < CrossingsNext.Get_Count(); i++) {
	            if (iX == CrossingsNext.Get_Point(i).x && iY == CrossingsNext.Get_Point(i).y) {
	                CrossingsAtPrincipal.Add(Crossings.Get_Point(i));
	            }//if
	        }// for
	    }while (!(iX == m_iClosingX && iY == m_iClosingY)
	            && (iX != iNextX || iY != iNextY));
	    Crossings = CrossingsAtPrincipal;			
	}// if
	
	P.x = m_iClosingX;
	P.y = m_iClosingY;
	Crossings.Add(P);
	
	bool *bCalculated = new bool[Crossings.Get_Count()];
	for (i = 0; i < Crossings.Get_Count(); i++){
		bCalculated[i] = false;
	}//for

	m_fMaxDistance = new float[Crossings.Get_Count()+1];
	m_fMaxDistance[0] = fMaxDistance;
	m_fHeightDif = new float[Crossings.Get_Count()+1];
	m_fHeightDif[0] = fHeightDif;

	m_pDistanceGrid = SG_Create_Grid(m_pDEM,GRID_TYPE_Float);
	m_pDistanceGrid->Assign((double)0);

	m_iNumBasins = 1;
	for (i = 0; i < Crossings.Get_Count(); i++) {
	    if (!bCalculated[i]) {
	        if (isTopHeader(Crossings, i, bCalculated)) {
				m_fCells=0;
	            WriteBasin(Crossings.Get_Point(i).x, Crossings.Get_Point(i).y, m_iNumBasins);
	            bCalculated[i] = true;
				if (m_fCells<100){
					DeleteBasin(Crossings.Get_Point(i).x, Crossings.Get_Point(i).y, m_iNumBasins);
					m_iNumBasins--;
				}//if
	            i = 0;
				m_iNumBasins++;
	        }// if
	    }// if
	}// for
		
	m_pBasinGrid->Set_NoData_Value(0);		
	
	//DataObject_Add(m_pDistanceGrid);

	
}// method

void CWatersheds_ext::CreateShapesLayer(){ //first shape (0) is the whole basin. 

	CSG_Table *pTable;
	CSG_Table_Record *pRecord, *pRecord2;
	CSG_Shape *pSubbasin;
	TSG_Point	Point;	
	float fArea=0, fPerim=0;
	float fSide1, fSide2;
	float fConcTime;
	float fMinHeight;
	int i,j,k;	
	int x,y;
	int iNextX, iNextY;
	int iX, iY;
	int iXOrig, iYOrig;
	int iIndex;
	int iUpstreamBasin;
	int iDownstreamBasin;
	int iCode;
	float * fCN;
	float * fSoilLoss;
	float * fCNValidCells;
	float * fSoilLossValidCells;
	CSG_String sSubbasins;
	int iOffsetX[] = {0, 1, 1, 1, 0, -1, -1, -1};
    int iOffsetY[] = {-1, -1, 0, 1, 1, 1, 0, -1};
    int iDirection;
    int iDir = 1;
    int iNewDirection[] = {5, 6, 7, 0, 1, 2, 3, 4};
		
	
	Process_Set_Text(_TL("Vectorizing headers..."));
		
	m_pHeaders->Create(SHAPE_TYPE_Point, _TL("Headers"));
	pTable=&m_pHeaders->Get_Table();
	pTable->Add_Field("Header", TABLE_FIELDTYPE_Int);
	for (i = 0; i < m_Headers.Get_Count(); i++){
		m_pHeaders->Add_Shape()->Add_Point(m_Headers.Get_Point(i).x * m_pDEM->Get_Cellsize() + m_pDEM->Get_XMin(), 
											m_Headers.Get_Point(i).y * m_pDEM->Get_Cellsize() + m_pDEM->Get_YMin());
		m_pHeaders->Get_Shape(i)->Get_Record()->Set_Value(0,i);
	}//for

	Process_Set_Text(_TL("Vectorizing subbasins..."));

	m_pBasins->Create(SHAPE_TYPE_Polygon, _TL("Subbasins"));
	pTable = &m_pBasins->Get_Table();
	pTable->Add_Field(_TL("Basin Code"), TABLE_FIELDTYPE_Int);

	pSubbasin = m_pBasins->Add_Shape();		

	for (x = 0; x < m_pDEM->Get_NX(); x++) {
	    for (y = 0; y < m_pDEM->Get_NY(); y++) {
	        if (m_pBasinGrid->asInt(x,y) !=0) {
	            iX = x;
	            iY = y;
	            goto out;
	        }// if		
	    }// for
	}// for
out:
	iXOrig = iX;
	iYOrig = iY;
	pSubbasin->Add_Point(iX * m_pDEM->Get_Cellsize() + m_pDEM->Get_XMin(),
						iY* m_pDEM->Get_Cellsize() + m_pDEM->Get_YMin());
	iDirection = 1;
	do {
	    if (iDirection > 7) {
	        iDirection =  iDirection % 8;
	    }//if
	    for (i = iDirection; i < iDirection + 8; i++) {
	        if (i > 7) {
	            iDir = i % 8;
	        }//if
	        else {
	            iDir = i;
	        }//else
	        if (m_pBasinGrid->asInt(iX + iOffsetX[iDir],iY + iOffsetY[iDir]) !=0) {
	            iX = iX + iOffsetX[iDir];
				iY = iY + iOffsetY[iDir];	                        
	            iDirection = iNewDirection[iDir];
	            break;
	        }// if
	    }// for
	    pSubbasin->Add_Point(iX * m_pDEM->Get_Cellsize() + m_pDEM->Get_XMin(),
							iY* m_pDEM->Get_Cellsize() + m_pDEM->Get_YMin());
	}while ((iY != iYOrig) || (iX != iXOrig));
	
	pSubbasin->Get_Record()->Set_Value(0,(double)0);

	for (k = 1; k < m_iNumBasins; k++) {

		pSubbasin = m_pBasins->Add_Shape();		

	    for (x = 0; x < m_pDEM->Get_NX(); x++) {
	        for (y = 0; y < m_pDEM->Get_NY(); y++) {
	            if (m_pBasinGrid->asInt(x,y) == k) {
	                iX = x;
	                iY = y;
	                goto out2;
	            }// if		
	        }// for
	    }// for
out2:
	    iXOrig = iX;
	    iYOrig = iY;
	    pSubbasin->Add_Point(iX * m_pDEM->Get_Cellsize() + m_pDEM->Get_XMin(),
							iY* m_pDEM->Get_Cellsize() + m_pDEM->Get_YMin());
	    iDirection = 1;
	    do {
	        if (iDirection > 7) {
	            iDirection =  iDirection % 8;
	        }//if
	        for (i = iDirection; i < iDirection + 8; i++) {
	            if (i > 7) {
	                iDir = i % 8;
	            }//if
	            else {
	                iDir = i;
	            }//else
	            if (m_pBasinGrid->asInt(iX + iOffsetX[iDir],iY + iOffsetY[iDir]) == k) {