1 /**********************************************************
2 * Version $Id: Helper.cpp 1016 2011-04-27 18:40:36Z oconrad $
3 *********************************************************/
4 /*******************************************************************************
5 Helper.cpp
6 Copyright (C) Victor Olaya
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin Street, 5th Floor, Boston, MA 02110-1301, USA
21 *******************************************************************************/
22
23 #include "Helper.h"
24
Init_FlowDirectionsD8(CSG_Grid * pDEM,CSG_Grid * pDirection)25 void Init_FlowDirectionsD8(CSG_Grid *pDEM, CSG_Grid *pDirection)
26 {
27 for(int y=0; y<pDEM->Get_NY() && SG_UI_Process_Set_Progress(y, pDEM->Get_NY()); y++)
28 {
29 #pragma omp parallel for
30 for(int x=0; x<pDEM->Get_NX(); x++)
31 {
32 pDirection->Set_Value(x, y, pDEM->Get_Gradient_NeighborDir(x, y), false);
33 }
34 }
35 }
36
37
getNextCell(CSG_Grid * g,int iX,int iY,int & iNextX,int & iNextY)38 void getNextCell(
39 CSG_Grid *g,
40 int iX,
41 int iY,
42 int &iNextX,
43 int &iNextY) {
44
45 float fDist;
46 float fMaxSlope;
47 float fSlope;
48
49 fMaxSlope = 0;
50 fSlope = 0;
51
52 if (iX < 1 || iX >= g->Get_NX() - 1 || iY < 1 || iY >= g->Get_NY() - 1
53 || g->is_NoData(iX,iY)) {
54 iNextX = iX;
55 iNextY = iY;
56 return;
57 }// if
58
59 for (int i = -1; i < 2; i++) {
60 for (int j = -1; j < 2; j++) {
61 if (!g->is_NoData(iX + i, iY + i)){
62 if (i == 0 || j == 0) {
63 fDist = 1.0f;
64 }// if
65 else {
66 fDist = 1.414f;
67 }// else
68 fSlope = (g->asFloat(iX+i,iY+j)
69 - g->asFloat(iX,iY)) / fDist;
70 if (fSlope <= fMaxSlope) {
71 iNextX = iX+i;
72 iNextY = iY+j;
73 fMaxSlope = fSlope;
74 }// if
75 }//if
76 }// for
77 }// for
78
79 }// method
80
getNextCell(CSG_Grid * g,CSG_Grid * g2,int iX,int iY,int & iNextX,int & iNextY)81 void getNextCell(
82 CSG_Grid *g,
83 CSG_Grid *g2,
84 int iX,
85 int iY,
86 int &iNextX,
87 int &iNextY) {
88
89 double fDist;
90 double fMaxSlope;
91 double fSlope;
92
93 fMaxSlope = 0.0000001;
94 fSlope = 0;
95
96 if (iX < 1 || iX >= g->Get_NX() - 1 || iY < 1 || iY >= g->Get_NY() - 1
97 || g->is_NoData(iX,iY)) {
98 iNextX = iX;
99 iNextY = iY;
100 return;
101 }// if
102
103 for (int i = -1; i < 2; i++) {
104 for (int j = -1; j < 2; j++) {
105 if (!g->is_NoData(iX+i,iY+j) &&
106 !g2->is_NoData(iX+i,iY+j)){
107 if (i == 0 || j == 0) {
108 fDist = 1.0f;
109 }// if
110 else {
111 fDist = 1.414f;
112 }// else
113 fSlope = (g->asFloat(iX+i,iY+j)
114 - g->asFloat(iX,iY)) / fDist;
115 if (fSlope < fMaxSlope) {
116 iNextX = iX+i;
117 iNextY = iY+j;
118 fMaxSlope = fSlope;
119 }// if
120 }//if
121 }// for
122 }// for
123
124 }// method
125
FlowDistance(CSG_Grid * pDEM,CSG_Grid * pBasinGrid,int iBasin,int iX,int iY,int iX2,int iY2)126 double FlowDistance(CSG_Grid *pDEM,
127 CSG_Grid *pBasinGrid,
128 int iBasin,
129 int iX,
130 int iY,
131 int iX2,
132 int iY2){// result in m, coords in grid coords. Returns 0 if no distance is calculated.
133
134 bool bIsInBasin;
135 double dDist = 1;
136 int iNextX = iX;
137 int iNextY = iY;
138
139 if (iX2 <= 0 || iX2 >= pDEM->Get_NX() || iY2 <= 0 || iY2 >= pDEM->Get_NY() ||
140 iX <= 0 || iX >= pDEM->Get_NX() || iY <= 0 || iY >= pDEM->Get_NY() ) {
141 return 0;
142 }// if
143 do {
144 iX = iNextX;
145 iY = iNextY;
146 getNextCell(pDEM,iX,iY,iNextX,iNextY);
147 if (fabs((double)(iX - iNextX + iY - iNextY)) == 1.0) {
148 dDist = dDist + pDEM->Get_Cellsize();
149 }// if
150 else {
151 dDist = dDist + 1.414 * pDEM->Get_Cellsize();
152 }// else
153 if (iX == iX2 && iY == iY2) {
154 return dDist;
155 }// if
156 if (iBasin == GLOBAL_BASIN){
157 bIsInBasin = !pBasinGrid->is_NoData(iX,iY);
158 }//if
159 else{
160 bIsInBasin = (pBasinGrid->asInt(iX,iY) == iBasin);
161 }//else
162 }while (bIsInBasin
163 &&(iX!=iNextX || iY!=iNextY));
164
165 return 0;
166
167 }// method
168
169
AccFlow(CSG_Grid * pGrid,CSG_Grid * pDEM,int iX,int iY)170 double AccFlow(CSG_Grid *pGrid, CSG_Grid *pDEM, int iX, int iY){
171
172 int iNextX, iNextY;
173 double dAccFlow = pGrid->Get_Cellsize()*pGrid->Get_Cellsize();
174
175 for (int i = -1; i<2; i++){
176 for (int j = -1; j<2; j++){
177 if (!(i == 0) || !(j == 0)) {
178 getNextCell(pDEM, iX + i, iY + j, iNextX, iNextY);
179 if (iNextY == iY && iNextX == iX) {
180 if (pGrid->asDouble(iX+i,iY+j)!=0){
181 dAccFlow += pGrid->asDouble(iX+i,iY+j);
182 }//if
183 else{
184 dAccFlow += AccFlow (pGrid,pDEM,iX+i,iY+j);
185 }//else
186 }// if
187 }//if
188 }//for
189 }//for
190
191 pGrid->Set_Value(iX,iY,dAccFlow);
192
193 return dAccFlow;
194
195 }// function
196
CalculateBasinGrid(CSG_Grid * pBasinGrid,CSG_Grid * pDEM,int iOutletX,int iOutletY)197 double CalculateBasinGrid(CSG_Grid *pBasinGrid,
198 CSG_Grid *pDEM,
199 int iOutletX,
200 int iOutletY){
201
202 pBasinGrid->Assign((double)0);
203 double dAccFlow = AccFlow(pBasinGrid,pDEM,iOutletX,iOutletY);
204
205 return dAccFlow;
206
207 }//function
208
CalculateFlowAccGrid(CSG_Grid * pFlowAccGrid,CSG_Grid * pDEM)209 void CalculateFlowAccGrid(CSG_Grid *pFlowAccGrid,
210 CSG_Grid *pDEM){
211
212 int x,y;
213
214 pFlowAccGrid->Assign((double)0);
215 for(y=0; y<pDEM->Get_NY(); y++){
216 for(x=0; x<pDEM->Get_NX(); x++){
217 AccFlow(pFlowAccGrid,pDEM,x,y);
218 }//for
219 }//for
220
221 pFlowAccGrid->Set_Description(_TL("Acc. Area"));
222 pFlowAccGrid->Set_Unit(_TL("m2"));
223
224 }//function
225
226
227 /*
228 #include <vector>
229
230 TSG_Point ** RiverProfile(int iX,
231 int iY,
232 CSG_Grid* pDEM,
233 CSG_Grid* pBasinGrid,
234 CSG_Grid* pExtGrid,
235 int &iProfileLength){
236
237 int i;
238 float fLength = 0;
239 int iNextX, iNextY;
240 CSG_Points Profile;
241 CSG_Points Ext;
242 TSG_Point P;
243 TSG_Point **pProfile;
244 iProfileLength = 0;
245
246 if (!pBasinGrid->is_NoData(iX,iY)) {
247 iNextX = iX;
248 iNextY = iY;
249 do {
250 iX = iNextX;
251 iY = iNextY;
252 getNextCell(pDEM, iX, iY, iNextX, iNextY);
253
254 if (fabs(iX - iNextX + iY - iNextY) == 1) {
255 fLength += pDEM->Get_DX();
256 }//if
257 else {
258 fLength += (1.414f * pDEM->Get_DX());
259 }//else
260 P.x = fLength;
261 P.y = pDEM->asFloat(iNextX, iNextY);
262 Profile.push_back(P);
263 P.y = pExtGrid->asFloat(iNextX, iNextY);
264 Ext.push_back(P);
265 }while (!pBasinGrid->is_NoData(iX, iY)
266 && (iX != iNextX || iY != iNextY));
267
268 pProfile = new Pt*[2];
269 for (i = 0; i<2; i++){
270 pProfile[i] = new Pt [Profile.size()];
271 }//for
272 for (i = 0; i<Profile.size(); i++){
273 pProfile[0][i]=Profile.at(i);
274 pProfile[1][i]=Ext.at(i);
275 }//for
276 iProfileLength = Profile.size();
277
278 }// if
279
280 return pProfile;
281
282 }//method
283
284 Pt* RiverCoords(int iX, //the resulting coords are grid coords,
285 int iY,
286 CSG_Grid* pDEM,
287 CSG_Grid* pBasinGrid,
288 int &iProfileLength){
289
290 float fLength = 0;
291 int iNextX, iNextY;
292 CSG_Points Profile;
293 TSG_Point P;
294 TSG_Point *pProfile;
295 iProfileLength = 0;
296
297 if (!pBasinGrid->is_NoData(iX,iY)) {
298 iNextX = iX;
299 iNextY = iY;
300 do {
301 iX = iNextX;
302 iY = iNextY;
303 getNextCell(pDEM, iX, iY, iNextX, iNextY);
304 P.x = iX;
305 P.y = iY;
306 Profile.push_back(P);
307 }while (!pBasinGrid->is_NoData(iX, iY)
308 && (iX != iNextX || iY != iNextY));
309
310 pProfile = new Pt [Profile.size()];
311 for (int i = 0; i<Profile.size(); i++){
312 pProfile[i]=Profile.at(i);
313 }//for
314 iProfileLength = Profile.size();
315
316 }// if
317
318 return pProfile;
319
320 }//method
321
322 float DrainageDensity(CSG_Shapes *pHeaders,
323 CSG_Shapes *pBasins,
324 CSG_Grid *pBasinGrid,
325 CSG_Grid *pDEM){
326
327 CSG_Grid * pChannelsGrid;
328 float fLength = 0;
329 int iX, iY;
330 int iNextX, iNextY;
331
332 pChannelsGrid = new CSG_Grid(pDEM, SG_DATATYPE_Byte);
333 pChannelsGrid->Assign(0);
334
335 for (int i = 0; i < pHeaders->Get_Count(); i++) {
336
337 iX = (pHeaders->Get_Shape(i)->Get_Point(0).x - pBasinGrid->Get_XMin())/pBasinGrid->Get_DX();
338 iY = (pHeaders->Get_Shape(i)->Get_Point(0).y - pBasinGrid->Get_YMin())/pBasinGrid->Get_DX();
339
340 if (!pBasinGrid->is_NoData(iX,iY)) {
341 iNextX = iX;
342 iNextY = iY;
343 do {
344 iX = iNextX;
345 iY = iNextY;
346 getNextCell(pDEM, iX, iY, iNextX, iNextY);
347
348 if (fabs(iX - iNextX + iY - iNextY) == 1) {
349 fLength += pChannelsGrid->Get_DX();
350 }//if
351 else {
352 fLength += (1.414f * pChannelsGrid->Get_DX());
353 }//else
354 if (pChannelsGrid->asFloat(iX,iY) == 0) {
355 pChannelsGrid->Set_Value(iX, iY, 1);
356 }// if
357 else {
358 break;
359 }// else
360 }while (!pBasinGrid->is_NoData(iX, iY)
361 && (iX != iNextX || iY != iNextY));
362
363 }// if
364 }// for
365
366 float fArea = pBasins->Get_Shape(0)->asFloat(4);
367
368 return fLength / fArea / 10000.0;
369
370 }//method
371
372 void ClosingPoint(CSG_Grid* pDEM,
373 CSG_Grid* pBasinGrid,
374 int &iClosingX,
375 int &iClosingY){
376
377 int x,y;
378 int iX,iY;
379 int iNextX, iNextY;
380
381 for (x = 0; x<pBasinGrid->Get_NX(); x++){
382 for (y = 0; y<pBasinGrid->Get_NY(); y++){
383 if (!pBasinGrid->is_NoData(x,y)){
384 iX = x;
385 iY = y;
386 goto out;
387 }//if
388 }//for
389 }//for
390 return;
391 out:
392
393 iNextX = iX;
394 iNextY = iY;
395 do {
396 iX = iNextX;
397 iY = iNextY;
398 getNextCell(pDEM, iX, iY, iNextX, iNextY);
399 }while (!pBasinGrid->is_NoData(iX, iY)
400 && (iX != iNextX || iY != iNextY));
401
402 iClosingX = iX;
403 iClosingY = iY;
404
405 }//method*/