📄 watersheds_ext.cpp
字号:
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,k);
}// for
pTable->Add_Field(_TL("Closing Pt. X"), TABLE_FIELDTYPE_Double); //1
pTable->Add_Field(_TL("Closing Pt. Y"), TABLE_FIELDTYPE_Double); //2
pTable->Add_Field(_TL("Perimeter (m)"), TABLE_FIELDTYPE_Double); //3
pTable->Add_Field(_TL("Area (ha)"), TABLE_FIELDTYPE_Double); //4
pTable->Add_Field(_TL("Avg. Curve Number"), TABLE_FIELDTYPE_Double); //5
pTable->Add_Field(_TL("Soil losses (t/ha\xc2\xb7year)"), TABLE_FIELDTYPE_Double); //6
pTable->Add_Field(_TL("Concentration time (h)"), TABLE_FIELDTYPE_Double); //7
pTable->Add_Field(_TL("Upslope Basins"), TABLE_FIELDTYPE_String); //8
pTable->Add_Field(_TL("Downslope basins"), TABLE_FIELDTYPE_Int); //9
pTable->Add_Field(_TL("Basin Type (Gravelius)"), TABLE_FIELDTYPE_String); //10
pTable->Add_Field(_TL("Equivalente Rectangle (Side 1)(m)"), TABLE_FIELDTYPE_Double); //11
pTable->Add_Field(_TL("Equivalente Rectangle (Side 2)(m)"), TABLE_FIELDTYPE_Double); //12
pTable->Add_Field(_TL("Orographic coef."), TABLE_FIELDTYPE_Double); //13
pTable->Add_Field(_TL("Massivity coef."), TABLE_FIELDTYPE_Double); //14
pTable->Add_Field(_TL("Centroid X"), TABLE_FIELDTYPE_Double); //15
pTable->Add_Field(_TL("Centroid Y"), TABLE_FIELDTYPE_Double); //16
pTable->Add_Field(_TL("Maximum flow distance(m)"), TABLE_FIELDTYPE_Double); //17
for (i = 0; i < m_pBasins->Get_Count(); i++) {
pSubbasin = m_pBasins->Get_Shape(i);
fArea = ((CSG_Shape_Polygon*)pSubbasin)->Get_Area();
fPerim = ((CSG_Shape_Polygon*)pSubbasin)->Get_Perimeter();
Point = ((CSG_Shape_Polygon*)pSubbasin)->Get_Centroid();
pSubbasin->Get_Record()->Set_Value(15, Point.x);
pSubbasin->Get_Record()->Set_Value(16, Point.y);
EquivalentRectangle(fPerim,fArea, fSide1, fSide2);
pSubbasin->Get_Record()->Set_Value(3, fPerim);
pSubbasin->Get_Record()->Set_Value(4, fArea/10000.0);
pSubbasin->Get_Record()->Set_Value(10, GraveliusType(fPerim, fArea).c_str());
if (fSide1!=NO_EQUIVALENT_RECTANGLE){
pSubbasin->Get_Record()->Set_Value(11, SG_Get_String(fSide1).c_str());
pSubbasin->Get_Record()->Set_Value(12, SG_Get_String(fSide2).c_str());
}//if
else{
pSubbasin->Get_Record()->Set_Value(11, SG_T("---"));
pSubbasin->Get_Record()->Set_Value(12, SG_T("---"));
}//else
pSubbasin->Get_Record()->Set_Value(13, OrographicIndex(fPerim, fArea/10000));
pSubbasin->Get_Record()->Set_Value(14, MassivityIndex(fPerim, fArea/10000));
}//for
//upstream and downstream basins
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);
iDownstreamBasin = m_pBasinGrid->asInt(iNextX,iNextY);
iUpstreamBasin = m_pBasinGrid->asInt(iX,iY);
if (iUpstreamBasin != iDownstreamBasin && iDownstreamBasin!=0 && iUpstreamBasin !=0) {
for (j = 0; j < pTable->Get_Record_Count(); j++) {
pRecord = pTable->Get_Record(j);
iCode = pRecord->asInt(0);
if (iCode == iUpstreamBasin){
pRecord->Set_Value(9, iDownstreamBasin);
pRecord->Set_Value(1, iX * m_pDEM->Get_Cellsize() + m_pDEM->Get_XMin());
pRecord->Set_Value(2, iY * m_pDEM->Get_Cellsize() + m_pDEM->Get_YMin());
}//if
}// for
}// if
}while (!(iX == m_iClosingX && iY == m_iClosingY)
&& (iX != iNextX || iY != iNextY));
}// for
pRecord = pTable->Get_Record(pTable->Get_Record_Count()-1);
pRecord->Set_Value(1, iX * m_pDEM->Get_Cellsize() + m_pDEM->Get_XMin());
pRecord->Set_Value(2, iY * m_pDEM->Get_Cellsize() + m_pDEM->Get_YMin());
for (i = 1; i<pTable->Get_Record_Count(); i++){
pTable->Get_Record(i)->Set_Value(8,SG_T(" --- "));
iCode = pTable->Get_Record(i)->asInt(0);
for (j = 0; j<pTable->Get_Record_Count(); j++){
iDownstreamBasin = pTable->Get_Record(j)->asInt(9);
if (iDownstreamBasin == iCode){
sSubbasins = CSG_String(pTable->Get_Record(i)->asString(8)) + SG_T(" ")
+ SG_Get_String(pTable->Get_Record(j)->asInt(0));
pTable->Get_Record(i)->Set_Value(8, sSubbasins.c_str());
}//if
}//for
}//for*/
fCN = new float[m_pBasins->Get_Count()];
fCNValidCells = new float[m_pBasins->Get_Count()];
fSoilLoss = new float[m_pBasins->Get_Count()];
fSoilLossValidCells = new float[m_pBasins->Get_Count()];
for (i = 0; i < m_pBasins->Get_Count(); i++){
fCN[i] = 0;
fCNValidCells[i] = 0;
fSoilLoss[i] = 0;
fSoilLossValidCells[i] = 0;
}//for
for(y=0; y<Get_NY() && Set_Progress(y); y++){
for(x=0; x<Get_NX(); x++){
if (!m_pBasinGrid->is_NoData(x,y)){
iIndex = m_pBasinGrid->asInt(x,y);
fSoilLossValidCells[iIndex]++;
fSoilLossValidCells[0]++;
if (m_pSoilLossGrid){
fSoilLoss[iIndex] += m_pSoilLossGrid->asFloat(x,y);
fSoilLoss[0] += m_pSoilLossGrid->asFloat(x,y);
}//if
fCNValidCells[iIndex]++;
fCNValidCells[0]++;
if (m_pCNGrid){
fCN[iIndex] += m_pCNGrid->asFloat(x,y);
fCN[0] += m_pCNGrid->asFloat(x,y);
}//if
}//if
}//for
}//for
for (i = 0; i < m_pBasins->Get_Count(); i++){
fCN[i] /= fCNValidCells[i];
fSoilLoss[i] /= fSoilLossValidCells[i];
}//for
for (j = 1; j < pTable->Get_Record_Count(); j++) {
pRecord = pTable->Get_Record(j);
iCode = pRecord->asInt(0);
pRecord->Set_Value(5, fCN[iCode]);
pRecord->Set_Value(6, fSoilLoss[iCode]);
iX = (pRecord->asInt(1) - m_pDEM->Get_XMin()) / m_pDEM->Get_Cellsize();
iY = (pRecord->asInt(2) - m_pDEM->Get_YMin()) / m_pDEM->Get_Cellsize();
fMinHeight = m_pDEM->asFloat(iX,iY);
fConcTime = pow(0.87 * pow(m_fMaxDistance[iCode] / 1000.0, 3.0)
/ (m_fHeightDif[iCode]-fMinHeight), 0.385);
pRecord->Set_Value(7, fConcTime);
pRecord->Set_Value(17, m_fMaxDistance[iCode]);
}// for */
pRecord = pTable->Get_Record(0);
fConcTime = pow(0.87 * pow(m_fMaxDistance[0] / 1000.0, 3.0)
/ (m_fHeightDif[0]-fMinHeight), 0.385);
pRecord->Set_Value(5, fCN[0]);
pRecord->Set_Value(6, fSoilLoss[0]);
pRecord->Set_Value(7, fConcTime);
pRecord2 = pTable->Get_Record(j-1);
pRecord->Set_Value(1, pRecord2->asInt(1));
pRecord->Set_Value(2, pRecord2->asInt(2));
pRecord->Set_Value(17, m_fMaxDistance[0]);
}//method
bool CWatersheds_ext::isHeader(int iX, int iY){
int iNextX, iNextY;
for (int i = -1; i < 2; i++) {
for (int j = -1; j < 2; j++) {
if (is_InGrid(iX+i,iY+j)){
if (!m_pChannelsGrid->is_NoData(iX+i,iY+j) && (i!=0 || j!=0)){
getNextCell(m_pDEM, m_pChannelsGrid, iX + i, iY + j, iNextX, iNextY);
if (iNextX == iX && iNextY == iY) {
return false;
}//if
}//if
}// if
}// for
}// for
return true;
}// method
bool CWatersheds_ext::isLastCell(int iX, int iY){
int iNextX, iNextY;
getNextCell(m_pDEM, m_pChannelsGrid, iX, iY, iNextX, iNextY);
if ((iNextX == iX && iNextY == iY) ||
(m_pChannelsGrid->is_NoData(iNextX,iNextY))) {
return true;
}//if
return false;
}// method
void CWatersheds_ext::WriteBasin(int iX, int iY, int iBasinNumber) {
int iNextX, iNextY;
float fDistance;
if ((iX <= 0 || iX >= Get_NX() || iY <= 0 || iY >= Get_NY())
|| m_pDEM->is_NoData(iX,iY)) {} else {
m_pBasinGrid->Set_Value(iX,iY, iBasinNumber);
m_fCells++;
for (int i = -1; i < 2; i++) {
for (int j = -1; j < 2; j++) {
if (!(i == 0) || !(j == 0)) {
if (m_pBasinGrid->asInt(iX + i, iY + j) == 0) {
getNextCell(m_pDEM, iX + i, iY + j, iNextX, iNextY);
if (iNextX == iX && iNextY == iY) {
fDistance = m_pDistanceGrid->asDouble(iX,iY)
+ M_GET_LENGTH(i,j)*m_pDEM->Get_Cellsize();
m_pDistanceGrid->Set_Value(iX+i, iY+j,fDistance) ;
if (fDistance > m_fMaxDistance[iBasinNumber]){
m_fMaxDistance[iBasinNumber] = fDistance;
m_fHeightDif[iBasinNumber] = m_pDEM->asFloat(iX+i,iY+j);
}//if
WriteBasin(iX + i, iY + j, iBasinNumber);
}//if
}// if
}// if
}// for
}// for
}// else
}// method
void CWatersheds_ext::DeleteBasin(int iX, int iY, int iBasinNumber) {
int iNextX, iNextY;
if ((iX <= 0 || iX >= Get_NX() || iY <= 0 || iY >= Get_NY())
|| m_pDEM->is_NoData(iX,iY)) {} else {
m_pBasinGrid->Set_Value(iX,iY, 0);
for (int i = -1; i < 2; i++) {
for (int j = -1; j < 2; j++) {
if (!(i == 0) || !(j == 0)) {
if (m_pBasinGrid->asInt(iX + i, iY + j) == iBasinNumber) {
getNextCell(m_pDEM, iX + i, iY + j, iNextX, iNextY);
if (iNextX == iX && iNextY == iY) {
DeleteBasin(iX + i, iY + j, iBasinNumber);
}//if
}// if
}// if
}// for
}// for
}// else
}// method
bool CWatersheds_ext::isTopHeader(CSG_Points &Headers,
int iIndex,
bool* bCalculated) {
int iX, iY;
int iNextX, iNextY;
int iTargetX = Headers.Get_Point(iIndex).x;
int iTargetY = Headers.Get_Point(iIndex).y;
for (int i = 0; i < Headers.Get_Count(); i++) {
if (i != iIndex && !bCalculated[i]) {
iNextX = Headers.Get_Point(i).x;
iNextY = Headers.Get_Point(i).y;
do {
iX = iNextX;
iY = iNextY;
getNextCell(m_pDEM, m_pChannelsGrid, iX, iY, iNextX, iNextY);
if (iNextX == iTargetX && iNextY == iTargetY) {
return false;
}//if
}while (!(iX == m_iClosingX && iY == m_iClosingY)
&& !(iNextX == iX && iNextY == iY));
}// if
}// for
return true;
}// method
float CWatersheds_ext::DistanceToClosingPoint(int iX, int iY){ // in m
float fDist = 1;
int iNextX = iX;
int iNextY = iY;
if ((iX <= 0 || iX >= Get_NX() || iY <= 0 || iY >= Get_NY())
|| m_pDEM->is_NoData(iX,iY)) {
return 0;
}// if
do {
iX = iNextX;
iY = iNextY;
getNextCell(m_pDEM, iX, iY, iNextX, iNextY);
if (fabs((double)iX - iNextX + iY - iNextY) == 1) {
fDist = fDist + Get_Cellsize();
}// if
else {
fDist = fDist + 1.414f * Get_Cellsize();
}// else
if (iX == m_iClosingX && iY == m_iClosingY) {
return fDist;
}// if
}while (iX != iNextX || iY != iNextY);
return 0;
}// method
CSG_String CWatersheds_ext::GraveliusType(float fPerimeter, //in m
float fArea){ //in m2
CSG_String sType;
float fGraveliusIndex = 0.28 * fPerimeter / sqrt(fArea);
if (fGraveliusIndex > 1.75) {
sType = "Rectangular";
}//if
else if (fGraveliusIndex > 1.5) {
sType = "Ovalooblonga-rectangularoblonga";
}//else if
else if (fGraveliusIndex > 1.5) {
sType = "Ovaloredonda-ovalooblonga";
}//else if
else {
sType = "Redonda-ovaloredonda";
}//else
return sType;
}// method
void CWatersheds_ext::EquivalentRectangle(float fPerimeter, //in m
float fArea, // in m2
float &fSide1, //in m
float &fSide2){// in m
float fDisc = pow(fPerimeter, 2) - 8.0 * fArea;
if (fDisc > 0) {
fSide1 = (fPerimeter +sqrt(fDisc)) / 4.0;
fSide2 = (fPerimeter - (2 * fSide1)) / 2.0;
}// if
else {
fSide1 = NO_EQUIVALENT_RECTANGLE;
fSide2 = NO_EQUIVALENT_RECTANGLE;
}// if
}// method
float CWatersheds_ext::MassivityIndex(float fMeanHeight,//in m
float fArea){ //in ha
return fMeanHeight / fArea;
}// method
float CWatersheds_ext::OrographicIndex(float fMeanHeight,//in m
float fArea){ //in ha
return fMeanHeight * fMeanHeight / fArea;
}// method
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -