1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // grid_operation.cpp //
15 // //
16 // Copyright (C) 2005 by Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. //
22 // //
23 // This library is free software; you can redistribute //
24 // it and/or modify it under the terms of the GNU Lesser //
25 // General Public License as published by the Free //
26 // Software Foundation, either version 2.1 of the //
27 // License, or (at your option) any later version. //
28 // //
29 // This library is distributed in the hope that it will //
30 // be useful, but WITHOUT ANY WARRANTY; without even the //
31 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
32 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
33 // License for more details. //
34 // //
35 // You should have received a copy of the GNU Lesser //
36 // General Public License along with this program; if //
37 // not, see <http://www.gnu.org/licenses/>. //
38 // //
39 //-------------------------------------------------------//
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Goettingen //
44 // Goldschmidtstr. 5 //
45 // 37077 Goettingen //
46 // Germany //
47 // //
48 // e-mail: oconrad@saga-gis.org //
49 // //
50 ///////////////////////////////////////////////////////////
51
52 //---------------------------------------------------------
53 #include <memory.h>
54
55 #include "grid.h"
56
57
58 ///////////////////////////////////////////////////////////
59 // //
60 // Data Assignments - Value //
61 // //
62 ///////////////////////////////////////////////////////////
63
64 //---------------------------------------------------------
Assign_NoData(void)65 void CSG_Grid::Assign_NoData(void)
66 {
67 double Value = Get_NoData_Value();
68
69 #pragma omp parallel for
70 for(sLong i=0; i<Get_NCells(); i++)
71 {
72 Set_Value(i, Value, false);
73 }
74 }
75
76 //---------------------------------------------------------
Assign(double Value)77 bool CSG_Grid::Assign(double Value)
78 {
79 if( !is_Valid() )
80 {
81 return( false );
82 }
83
84 if( Value == 0. && !is_Cached() )
85 {
86 #pragma omp parallel for
87 for(int y=0; y<Get_NY(); y++)
88 {
89 memset(m_Values[y], 0, Get_nLineBytes());
90 }
91 }
92 else
93 {
94 #pragma omp parallel for
95 for(int y=0; y<Get_NY(); y++)
96 {
97 for(int x=0; x<Get_NX(); x++)
98 {
99 Set_Value(x, y, Value);
100 }
101 }
102 }
103
104 //-----------------------------------------------------
105 Get_History().Destroy();
106
107 m_Statistics.Invalidate();
108
109 Set_Update_Flag(false);
110
111 return( true );
112 }
113
114
115 ///////////////////////////////////////////////////////////
116 // //
117 // Data Assignments - Grid //
118 // //
119 ///////////////////////////////////////////////////////////
120
121 //---------------------------------------------------------
Assign(CSG_Data_Object * pObject)122 bool CSG_Grid::Assign(CSG_Data_Object *pObject)
123 {
124 return( pObject && pObject->is_Valid() && pObject->Get_ObjectType() == Get_ObjectType()
125 && Assign((CSG_Grid *)pObject, GRID_RESAMPLING_Undefined) );
126 }
127
Assign(CSG_Grid * pGrid,TSG_Grid_Resampling Interpolation)128 bool CSG_Grid::Assign(CSG_Grid *pGrid, TSG_Grid_Resampling Interpolation)
129 {
130 //-----------------------------------------------------
131 if( !is_Valid() || !pGrid || !pGrid->is_Valid() || is_Intersecting(pGrid->Get_Extent()) == INTERSECTION_None )
132 {
133 return( false );
134 }
135
136 bool bResult = false;
137
138 //---------------------------------------------------------
139 if( m_System == pGrid->m_System )
140 {
141 for(int y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
142 {
143 #pragma omp parallel for
144 for(int x=0; x<Get_NX(); x++)
145 {
146 if( pGrid->is_NoData(x, y) )
147 {
148 Set_NoData(x, y);
149 }
150 else
151 {
152 Set_Value(x, y, pGrid->asDouble(x, y));
153 }
154 }
155 }
156
157 bResult = true;
158 }
159
160 else if( Get_Cellsize() == pGrid->Get_Cellsize() // No-Scaling...
161 && fmod(Get_XMin() - pGrid->Get_XMin(), Get_Cellsize()) == 0.
162 && fmod(Get_YMin() - pGrid->Get_YMin(), Get_Cellsize()) == 0. )
163 {
164 bResult = _Assign_Interpolated(pGrid, GRID_RESAMPLING_NearestNeighbour);
165 }
166
167 //---------------------------------------------------------
168 else switch( Interpolation )
169 {
170 case GRID_RESAMPLING_NearestNeighbour:
171 case GRID_RESAMPLING_Bilinear:
172 case GRID_RESAMPLING_BicubicSpline:
173 case GRID_RESAMPLING_BSpline:
174 bResult = _Assign_Interpolated(pGrid, Interpolation);
175 break;
176
177 case GRID_RESAMPLING_Mean_Nodes:
178 case GRID_RESAMPLING_Mean_Cells:
179 bResult = _Assign_MeanValue (pGrid, Interpolation != GRID_RESAMPLING_Mean_Nodes);
180 break;
181
182 case GRID_RESAMPLING_Minimum:
183 case GRID_RESAMPLING_Maximum:
184 bResult = _Assign_ExtremeValue(pGrid, Interpolation == GRID_RESAMPLING_Maximum);
185 break;
186
187 case GRID_RESAMPLING_Majority:
188 bResult = _Assign_Majority (pGrid);
189 break;
190
191 default:
192 if( Get_Cellsize() < pGrid->Get_Cellsize() ) // Down-Scaling...
193 {
194 bResult = _Assign_Interpolated(pGrid, GRID_RESAMPLING_BSpline);
195 }
196 else // Up-Scaling...
197 {
198 bResult = _Assign_MeanValue(pGrid, Interpolation != GRID_RESAMPLING_Mean_Nodes);
199 }
200 break;
201 }
202
203 //---------------------------------------------------------
204 if( bResult )
205 {
206 Set_Unit(pGrid->Get_Unit());
207
208 if( pGrid->Get_Projection().is_Okay() )
209 {
210 Get_Projection() = pGrid->Get_Projection();
211 }
212
213 Get_History() = pGrid->Get_History();
214 }
215
216 //---------------------------------------------------------
217 SG_UI_Process_Set_Ready();
218
219 return( bResult );
220 }
221
222 //---------------------------------------------------------
_Assign_Interpolated(CSG_Grid * pGrid,TSG_Grid_Resampling Interpolation)223 bool CSG_Grid::_Assign_Interpolated(CSG_Grid *pGrid, TSG_Grid_Resampling Interpolation)
224 {
225 double py = Get_YMin();
226
227 for(int y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++, py+=Get_Cellsize())
228 {
229 #pragma omp parallel for
230 for(int x=0; x<Get_NX(); x++)
231 {
232 double z;
233
234 if( pGrid->Get_Value(Get_XMin() + x * Get_Cellsize(), py, z, Interpolation) )
235 {
236 Set_Value(x, y, z);
237 }
238 else
239 {
240 Set_NoData(x, y);
241 }
242 }
243 }
244
245 return( true );
246 }
247
248 //---------------------------------------------------------
_Assign_ExtremeValue(CSG_Grid * pGrid,bool bMaximum)249 bool CSG_Grid::_Assign_ExtremeValue(CSG_Grid *pGrid, bool bMaximum)
250 {
251 if( Get_Cellsize() < pGrid->Get_Cellsize() )
252 {
253 return( false );
254 }
255
256 //-----------------------------------------------------
257 Assign_NoData();
258
259 double ax = 0.5 + (pGrid->Get_XMin() - Get_XMin()) / Get_Cellsize();
260 double py = 0.5 + (pGrid->Get_YMin() - Get_YMin()) / Get_Cellsize();
261
262 double d = pGrid->Get_Cellsize() / Get_Cellsize();
263
264 for(int y=0; y<pGrid->Get_NY() && SG_UI_Process_Set_Progress(y, pGrid->Get_NY()); y++, py+=d)
265 {
266 int iy = (int)floor(py);
267
268 if( iy >= 0 && iy < Get_NY() )
269 {
270 #pragma omp parallel for
271 for(int x=0; x<pGrid->Get_NX(); x++)
272 {
273 if( !pGrid->is_NoData(x, y) )
274 {
275 int ix = (int)floor(ax + x * d);
276
277 if( ix >= 0 && ix < Get_NX() )
278 {
279 double z = pGrid->asDouble(x, y);
280
281 if( is_NoData(ix, iy)
282 || (bMaximum == true && z > asDouble(ix, iy))
283 || (bMaximum == false && z < asDouble(ix, iy)) )
284 {
285 Set_Value(ix, iy, z);
286 }
287 }
288 }
289 }
290 }
291 }
292
293 //-----------------------------------------------------
294 return( true );
295 }
296
297 //---------------------------------------------------------
_Assign_MeanValue(CSG_Grid * pGrid,bool bAreaProportional)298 bool CSG_Grid::_Assign_MeanValue(CSG_Grid *pGrid, bool bAreaProportional)
299 {
300 if( Get_Cellsize() < pGrid->Get_Cellsize() )
301 {
302 return( false );
303 }
304
305 //-----------------------------------------------------
306 double d = Get_Cellsize() / pGrid->Get_Cellsize();
307
308 double ox = (Get_XMin(true) - pGrid->Get_XMin()) / pGrid->Get_Cellsize();
309 double py = (Get_YMin(true) - pGrid->Get_YMin()) / pGrid->Get_Cellsize();
310
311 for(int y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++, py+=d)
312 {
313 int ay = (int)(bAreaProportional ? floor(py ) : ceil (py ));
314 int by = (int)(bAreaProportional ? ceil (py + d) : floor(py + d));
315
316 //-------------------------------------------------
317 #pragma omp parallel for
318 for(int x=0; x<Get_NX(); x++)
319 {
320 double px = ox + x * d;
321
322 int ax = (int)(bAreaProportional ? floor(px ) : ceil (px ));
323 int bx = (int)(bAreaProportional ? ceil (px + d) : floor(px + d));
324
325 CSG_Rect rMean(px, py, px + d, py + d);
326
327 CSG_Simple_Statistics s;
328
329 for(int iy=ay; iy<=by; iy++)
330 {
331 if( iy >= 0 && iy < pGrid->Get_NY() )
332 {
333 for(int ix=ax; ix<=bx; ix++)
334 {
335 if( ix >= 0 && ix < pGrid->Get_NX() && !pGrid->is_NoData(ix, iy) )
336 {
337 if( bAreaProportional )
338 {
339 CSG_Rect r(ix - 0.5, iy - 0.5, ix + 0.5, iy + 0.5);
340
341 if( r.Intersect(rMean) )
342 {
343 s.Add_Value(pGrid->asDouble(ix, iy), r.Get_Area());
344 }
345 }
346 else
347 {
348 s.Add_Value(pGrid->asDouble(ix, iy));
349 }
350 }
351 }
352 }
353 }
354
355 //---------------------------------------------
356 if( s.Get_Count() > 0 )
357 {
358 Set_Value(x, y, s.Get_Mean());
359 }
360 else
361 {
362 Set_NoData(x, y);
363 }
364 }
365 }
366
367 //-----------------------------------------------------
368 return( true );
369 }
370
371
372 ///////////////////////////////////////////////////////////
373 // //
374 ///////////////////////////////////////////////////////////
375
376 //---------------------------------------------------------
_Assign_Majority(CSG_Grid * pGrid)377 bool CSG_Grid::_Assign_Majority(CSG_Grid *pGrid)
378 {
379 if( Get_Cellsize() < pGrid->Get_Cellsize() )
380 {
381 return( false );
382 }
383
384 //-----------------------------------------------------
385 Set_NoData_Value(pGrid->Get_NoData_Value());
386
387 Assign_NoData();
388
389 //-----------------------------------------------------
390 int ay, by = (int)(1. + (((0 - 0.5) * Get_Cellsize() + Get_YMin()) - pGrid->Get_YMin()) / pGrid->Get_Cellsize());
391
392 for(int y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
393 {
394 ay = by;
395 by = (int)(1. + (((y + 0.5) * Get_Cellsize() + Get_YMin()) - pGrid->Get_YMin()) / pGrid->Get_Cellsize());
396
397 if( ay < pGrid->Get_NY() && by > 0 )
398 {
399 if( ay < 0 )
400 {
401 ay = 0;
402 }
403
404 if( by > pGrid->Get_NY() )
405 {
406 by = pGrid->Get_NY();
407 }
408
409 int ax, bx = (int)(1. + (((0 - 0.5) * Get_Cellsize() + Get_XMin()) - pGrid->Get_XMin()) / pGrid->Get_Cellsize());
410
411 for(int x=0; x<Get_NX(); x++)
412 {
413 ax = bx;
414 bx = (int)(1. + (((x + 0.5) * Get_Cellsize() + Get_XMin()) - pGrid->Get_XMin()) / pGrid->Get_Cellsize());
415
416 if( ax < pGrid->Get_NX() && bx > 0 )
417 {
418 CSG_Unique_Number_Statistics s;
419
420 if( ax < 0 )
421 {
422 ax = 0;
423 }
424
425 if( bx > pGrid->Get_NX() )
426 {
427 bx = pGrid->Get_NX();
428 }
429
430 for(int iy=ay; iy<by; iy++)
431 {
432 for(int ix=ax; ix<bx; ix++)
433 {
434 if( !pGrid->is_NoData(ix, iy) )
435 {
436 s += pGrid->asDouble(ix, iy);
437 }
438 }
439 }
440
441 int n; double z;
442
443 if( s.Get_Majority(z, n) )//&& n > 1 )
444 {
445 Set_Value(x, y, z);
446 }
447 }
448 }
449 }
450 }
451
452 //-----------------------------------------------------
453 return( true );
454 }
455
456
457 ///////////////////////////////////////////////////////////
458 // //
459 // Operatoren //
460 // //
461 ///////////////////////////////////////////////////////////
462
463 //---------------------------------------------------------
operator =(const CSG_Grid & Grid)464 CSG_Grid & CSG_Grid::operator = (const CSG_Grid &Grid)
465 {
466 Assign((CSG_Grid *)&Grid, GRID_RESAMPLING_Undefined);
467
468 return( *this );
469 }
470
operator =(double Value)471 CSG_Grid & CSG_Grid::operator = (double Value)
472 {
473 Assign(Value);
474
475 return( *this );
476 }
477
478 //---------------------------------------------------------
operator +(const CSG_Grid & Grid) const479 CSG_Grid CSG_Grid::operator + (const CSG_Grid &Grid) const
480 {
481 CSG_Grid g(*this);
482
483 return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Addition) );
484 }
485
operator +(double Value) const486 CSG_Grid CSG_Grid::operator + (double Value) const
487 {
488 CSG_Grid g(*this);
489
490 return( g._Operation_Arithmetic(Value, GRID_OPERATION_Addition) );
491 }
492
operator +=(const CSG_Grid & Grid)493 CSG_Grid & CSG_Grid::operator += (const CSG_Grid &Grid)
494 {
495 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Addition) );
496 }
497
operator +=(double Value)498 CSG_Grid & CSG_Grid::operator += (double Value)
499 {
500 return( _Operation_Arithmetic(Value, GRID_OPERATION_Addition) );
501 }
502
Add(const CSG_Grid & Grid)503 CSG_Grid & CSG_Grid::Add (const CSG_Grid &Grid)
504 {
505 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Addition) );
506 }
507
Add(double Value)508 CSG_Grid & CSG_Grid::Add (double Value)
509 {
510 return( _Operation_Arithmetic(Value, GRID_OPERATION_Addition) );
511 }
512
513 //---------------------------------------------------------
operator -(const CSG_Grid & Grid) const514 CSG_Grid CSG_Grid::operator - (const CSG_Grid &Grid) const
515 {
516 CSG_Grid g(*this);
517
518 return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Subtraction) );
519 }
520
operator -(double Value) const521 CSG_Grid CSG_Grid::operator - (double Value) const
522 {
523 CSG_Grid g(*this);
524
525 return( g._Operation_Arithmetic(Value, GRID_OPERATION_Subtraction) );
526 }
527
operator -=(const CSG_Grid & Grid)528 CSG_Grid & CSG_Grid::operator -= (const CSG_Grid &Grid)
529 {
530 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Subtraction) );
531 }
532
operator -=(double Value)533 CSG_Grid & CSG_Grid::operator -= (double Value)
534 {
535 return( _Operation_Arithmetic(Value, GRID_OPERATION_Subtraction) );
536 }
537
Subtract(const CSG_Grid & Grid)538 CSG_Grid & CSG_Grid::Subtract (const CSG_Grid &Grid)
539 {
540 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Subtraction) );
541 }
542
Subtract(double Value)543 CSG_Grid & CSG_Grid::Subtract (double Value)
544 {
545 return( _Operation_Arithmetic(Value, GRID_OPERATION_Subtraction) );
546 }
547
548 //---------------------------------------------------------
operator *(const CSG_Grid & Grid) const549 CSG_Grid CSG_Grid::operator * (const CSG_Grid &Grid) const
550 {
551 CSG_Grid g(*this);
552
553 return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Multiplication) );
554 }
555
operator *(double Value) const556 CSG_Grid CSG_Grid::operator * (double Value) const
557 {
558 CSG_Grid g(*this);
559
560 return( g._Operation_Arithmetic(Value, GRID_OPERATION_Multiplication) );
561 }
562
operator *=(const CSG_Grid & Grid)563 CSG_Grid & CSG_Grid::operator *= (const CSG_Grid &Grid)
564 {
565 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Multiplication) );
566 }
567
operator *=(double Value)568 CSG_Grid & CSG_Grid::operator *= (double Value)
569 {
570 return( _Operation_Arithmetic(Value, GRID_OPERATION_Multiplication) );
571 }
572
Multiply(const CSG_Grid & Grid)573 CSG_Grid & CSG_Grid::Multiply (const CSG_Grid &Grid)
574 {
575 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Multiplication) );
576 }
577
Multiply(double Value)578 CSG_Grid & CSG_Grid::Multiply (double Value)
579 {
580 return( _Operation_Arithmetic(Value, GRID_OPERATION_Multiplication) );
581 }
582
583 //---------------------------------------------------------
operator /(const CSG_Grid & Grid) const584 CSG_Grid CSG_Grid::operator / (const CSG_Grid &Grid) const
585 {
586 CSG_Grid g(*this);
587
588 return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Division) );
589 }
590
operator /(double Value) const591 CSG_Grid CSG_Grid::operator / (double Value) const
592 {
593 CSG_Grid g(*this);
594
595 return( g._Operation_Arithmetic(Value, GRID_OPERATION_Division) );
596 }
597
operator /=(const CSG_Grid & Grid)598 CSG_Grid & CSG_Grid::operator /= (const CSG_Grid &Grid)
599 {
600 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Division) );
601 }
602
operator /=(double Value)603 CSG_Grid & CSG_Grid::operator /= (double Value)
604 {
605 return( _Operation_Arithmetic(Value, GRID_OPERATION_Division) );
606 }
607
Divide(const CSG_Grid & Grid)608 CSG_Grid & CSG_Grid::Divide (const CSG_Grid &Grid)
609 {
610 return( _Operation_Arithmetic(Grid, GRID_OPERATION_Division) );
611 }
612
Divide(double Value)613 CSG_Grid & CSG_Grid::Divide (double Value)
614 {
615 return( _Operation_Arithmetic(Value, GRID_OPERATION_Division) );
616 }
617
618
619 ///////////////////////////////////////////////////////////
620 // //
621 // Operatoren //
622 // //
623 ///////////////////////////////////////////////////////////
624
625 //---------------------------------------------------------
_Operation_Arithmetic(const CSG_Grid & Grid,TSG_Grid_Operation Operation)626 CSG_Grid & CSG_Grid::_Operation_Arithmetic(const CSG_Grid &Grid, TSG_Grid_Operation Operation)
627 {
628 if( is_Intersecting(Grid.Get_Extent()) )
629 {
630 TSG_Grid_Resampling Interpolation =
631 Get_Cellsize() == Grid.Get_Cellsize() && fmod(Get_XMin() - Grid.Get_XMin(), Get_Cellsize()) == 0.
632 && Get_Cellsize() == Grid.Get_Cellsize() && fmod(Get_YMin() - Grid.Get_YMin(), Get_Cellsize()) == 0.
633 ? GRID_RESAMPLING_NearestNeighbour
634 : GRID_RESAMPLING_BSpline;
635
636 #pragma omp parallel for
637 for(int y=0; y<Get_NY(); y++)
638 {
639 double yWorld = Get_YMin() + y * Get_Cellsize();
640
641 for(int x=0; x<Get_NX(); x++)
642 {
643 if( !is_NoData(x, y) )
644 {
645 double xWorld = Get_XMin() + x * Get_Cellsize(), Value;
646
647 if( Grid.Get_Value(xWorld, yWorld, Value, Interpolation) )
648 {
649 switch( Operation )
650 {
651 case GRID_OPERATION_Addition : Add_Value(x, y, Value); break;
652 case GRID_OPERATION_Subtraction : Add_Value(x, y, -Value); break;
653 case GRID_OPERATION_Multiplication: Mul_Value(x, y, Value); break;
654 case GRID_OPERATION_Division :
655 if( Value != 0. )
656 {
657 Mul_Value(x, y, 1. / Value);
658 }
659 else
660 {
661 Set_NoData(x, y);
662 }
663 break;
664 }
665 }
666 }
667 }
668 }
669 }
670
671 return( *this );
672 }
673
674 //---------------------------------------------------------
_Operation_Arithmetic(double Value,TSG_Grid_Operation Operation)675 CSG_Grid & CSG_Grid::_Operation_Arithmetic(double Value, TSG_Grid_Operation Operation)
676 {
677 switch( Operation )
678 {
679 case GRID_OPERATION_Addition :
680 if( Value == 0. )
681 return( *this );
682 break;
683
684 case GRID_OPERATION_Subtraction :
685 if( Value == 0. )
686 return( *this );
687 Value = -Value;
688 break;
689
690 case GRID_OPERATION_Multiplication:
691 if( Value == 1. )
692 return( *this );
693 break;
694
695 case GRID_OPERATION_Division :
696 if( Value == 0. )
697 return( *this );
698 Value = 1. / Value;
699 break;
700 }
701
702 //-----------------------------------------------------
703 #pragma omp parallel for
704 for(int y=0; y<Get_NY(); y++)
705 {
706 for(int x=0; x<Get_NX(); x++)
707 {
708 if( !is_NoData(x, y) )
709 {
710 switch( Operation )
711 {
712 case GRID_OPERATION_Addition :
713 case GRID_OPERATION_Subtraction : Add_Value(x, y, Value); break;
714
715 case GRID_OPERATION_Multiplication:
716 case GRID_OPERATION_Division : Mul_Value(x, y, Value); break;
717 }
718 }
719 }
720 }
721
722 //-----------------------------------------------------
723 return( *this );
724 }
725
726
727 ///////////////////////////////////////////////////////////
728 // //
729 // Grid-Operations - A //
730 // //
731 ///////////////////////////////////////////////////////////
732
733 //---------------------------------------------------------
Invert(void)734 void CSG_Grid::Invert(void)
735 {
736 if( is_Valid() && Get_Range() > 0. )
737 {
738 double Min = Get_Min();
739 double Max = Get_Max();
740
741 #pragma omp parallel for
742 for(int y=0; y<Get_NY(); y++)
743 {
744 for(int x=0; x<Get_NX(); x++)
745 {
746 if( !is_NoData(x, y) )
747 {
748 Set_Value(x, y, Max - (asDouble(x, y) - Min));
749 }
750 }
751 }
752 }
753 }
754
755 //---------------------------------------------------------
Flip(void)756 void CSG_Grid::Flip(void)
757 {
758 if( is_Valid() )
759 {
760 #pragma omp parallel for
761 for(int x=0; x<Get_NX(); x++)
762 {
763 for(int yA=0, yB=Get_NY()-1; yA<yB; yA++, yB--)
764 {
765 double d = asDouble(x, yA);
766 Set_Value(x, yA, asDouble(x, yB));
767 Set_Value(x, yB, d);
768 }
769 }
770 }
771 }
772
773 //---------------------------------------------------------
Mirror(void)774 void CSG_Grid::Mirror(void)
775 {
776 if( is_Valid() )
777 {
778 #pragma omp parallel for
779 for(int y=0; y<Get_NY(); y++)
780 {
781 for(int xA=0, xB=Get_NX()-1; xA<xB; xA++, xB--)
782 {
783 double d = asDouble(xA, y);
784 Set_Value(xA, y, asDouble(xB, y));
785 Set_Value(xB, y, d);
786 }
787 }
788 }
789 }
790
791
792 ///////////////////////////////////////////////////////////
793 // //
794 // Grid-Operations - B //
795 // //
796 ///////////////////////////////////////////////////////////
797
798 //---------------------------------------------------------
Normalise(void)799 bool CSG_Grid::Normalise(void)
800 {
801 if( is_Valid() && Get_Range() > 0. )
802 {
803 double Min = Get_Min ();
804 double Range = Get_Range();
805
806 #pragma omp parallel for
807 for(int y=0; y<Get_NY(); y++)
808 {
809 for(int x=0; x<Get_NX(); x++)
810 {
811 if( !is_NoData(x, y) )
812 {
813 Set_Value(x, y, (asDouble(x, y) - Min) / Range);
814 }
815 }
816 }
817
818 return( true );
819 }
820
821 return( false );
822 }
823
824 //---------------------------------------------------------
DeNormalise(double Minimum,double Maximum)825 bool CSG_Grid::DeNormalise(double Minimum, double Maximum)
826 {
827 if( is_Valid() && Minimum < Maximum )
828 {
829 #pragma omp parallel for
830 for(int y=0; y<Get_NY(); y++)
831 {
832 for(int x=0; x<Get_NX(); x++)
833 {
834 if( !is_NoData(x, y) )
835 {
836 Set_Value(x, y, Minimum + asDouble(x, y) * (Maximum - Minimum));
837 }
838 }
839 }
840
841 return( true );
842 }
843
844 return( false );
845 }
846
847 //---------------------------------------------------------
Standardise(void)848 bool CSG_Grid::Standardise(void)
849 {
850 if( is_Valid() && Get_StdDev() > 0. )
851 {
852 double Mean = Get_Mean ();
853 double StdDev = Get_StdDev();
854
855 #pragma omp parallel for
856 for(int y=0; y<Get_NY(); y++)
857 {
858 for(int x=0; x<Get_NX(); x++)
859 {
860 if( !is_NoData(x, y) )
861 {
862 Set_Value(x, y, (asDouble(x, y) - Mean) / StdDev);
863 }
864 }
865 }
866
867 return( true );
868 }
869
870 return( false );
871 }
872
873 //---------------------------------------------------------
DeStandardise(double Mean,double StdDev)874 bool CSG_Grid::DeStandardise(double Mean, double StdDev)
875 {
876 if( is_Valid() && StdDev > 0. )
877 {
878 #pragma omp parallel for
879 for(int y=0; y<Get_NY(); y++)
880 {
881 for(int x=0; x<Get_NX(); x++)
882 {
883 if( !is_NoData(x, y) )
884 {
885 Set_Value(x, y, Mean + asDouble(x, y) * StdDev);
886 }
887 }
888 }
889
890 return( true );
891 }
892
893 return( false );
894 }
895
896
897 ///////////////////////////////////////////////////////////
898 // //
899 // //
900 // //
901 ///////////////////////////////////////////////////////////
902
903 //---------------------------------------------------------
904 /**
905 * Returns the direction to which the downward gradient is
906 * steepest. This implements the Deterministic 8 (D8) algorithm
907 * for identifying a single flow direction based on elevation
908 * data. Direction is numbered clock-wise beginning with 0 for
909 * the North. If no direction can be identified result will be a -1.
910 * If 'bDown' is not true the cell that the direction is pointing
911 * to might have a higher value than the center cell. If 'bNoEdges'
912 * is true a -1 will be returned for all cells that are at the
913 * edge of the data area.
914 */
915 //---------------------------------------------------------
Get_Gradient_NeighborDir(int x,int y,bool bDown,bool bNoEdges) const916 int CSG_Grid::Get_Gradient_NeighborDir(int x, int y, bool bDown, bool bNoEdges) const
917 {
918 int Direction = -1;
919
920 if( is_InGrid(x, y) )
921 {
922 double z = asDouble(x, y), dzMax = 0.;
923
924 for(int i=0; i<8; i++)
925 {
926 int ix = m_System.Get_xTo(i, x);
927 int iy = m_System.Get_yTo(i, y);
928
929 if( is_InGrid(ix, iy) )
930 {
931 double dz = (z - asDouble(ix, iy)) / m_System.Get_Length(i);
932
933 if( (!bDown || dz > 0.) && (Direction < 0 || dzMax < dz) )
934 {
935 dzMax = dz;
936 Direction = i;
937 }
938 }
939 else if( bNoEdges )
940 {
941 return( -1 );
942 }
943 }
944 }
945
946 return( Direction );
947 }
948
949 //---------------------------------------------------------
950 /**
951 * Calculates the gradient of a cell interpreting all grid cell values
952 * as elevation surface. Calculation uses the formulas proposed
953 * by Zevenbergen & Thorne (1986). Slope and aspect values are calculated
954 * in radians. Aspect is zero for the North direction and increases
955 * clockwise.
956 */
957 //---------------------------------------------------------
Get_Gradient(int x,int y,double & Slope,double & Aspect) const958 bool CSG_Grid::Get_Gradient(int x, int y, double &Slope, double &Aspect) const
959 {
960 if( is_InGrid(x, y) )
961 {
962 double z = asDouble(x, y), dz[4];
963
964 for(int i=0, iDir=0, ix, iy; i<4; i++, iDir+=2)
965 {
966 if( is_InGrid(
967 ix = m_System.Get_xTo (iDir, x),
968 iy = m_System.Get_yTo (iDir, y)) )
969 {
970 dz[i] = asDouble(ix, iy) - z;
971 }
972 else if( is_InGrid(
973 ix = m_System.Get_xFrom(iDir, x),
974 iy = m_System.Get_yFrom(iDir, y)) )
975 {
976 dz[i] = z - asDouble(ix, iy);
977 }
978 else
979 {
980 dz[i] = 0.;
981 }
982 }
983
984 double G = (dz[0] - dz[2]) / (2. * Get_Cellsize());
985 double H = (dz[1] - dz[3]) / (2. * Get_Cellsize());
986
987 Slope = atan(sqrt(G*G + H*H));
988 Aspect = G != 0. ? M_PI_180 + atan2(H, G) : H > 0. ? M_PI_270 : H < 0. ? M_PI_090 : -1.;
989
990 return( true );
991 }
992
993 Slope = 0.;
994 Aspect = -1.;
995
996 return( false );
997 }
998
999 //---------------------------------------------------------
1000 /**
1001 * Calculates the gradient for the given world coordinate.
1002 * Calculation uses the formulas proposed by Zevenbergen & Thorne (1986).
1003 * Slope and aspect values are calulated in radians.
1004 * Aspect is zero for the North direction and increases clockwise.
1005 */
1006 //---------------------------------------------------------
Get_Gradient(double x,double y,double & Slope,double & Aspect,TSG_Grid_Resampling Interpolation) const1007 bool CSG_Grid::Get_Gradient(double x, double y, double &Slope, double &Aspect, TSG_Grid_Resampling Interpolation) const
1008 {
1009 double z, iz, dz[4];
1010
1011 if( Get_Value(x, y, z, Interpolation) )
1012 {
1013 for(int i=0, iDir=0; i<4; i++, iDir+=2)
1014 {
1015 if( Get_Value(
1016 x + Get_Cellsize() * m_System.Get_xTo (iDir),
1017 y + Get_Cellsize() * m_System.Get_yTo (iDir), iz, Interpolation) )
1018 {
1019 dz[i] = iz - z;
1020 }
1021 else if( Get_Value(
1022 x + Get_Cellsize() * m_System.Get_xFrom(iDir),
1023 y + Get_Cellsize() * m_System.Get_yFrom(iDir), iz, Interpolation) )
1024 {
1025 dz[i] = z - iz;
1026 }
1027 else
1028 {
1029 dz[i] = 0.;
1030 }
1031 }
1032
1033 double G = (dz[0] - dz[2]) / (2. * Get_Cellsize());
1034 double H = (dz[1] - dz[3]) / (2. * Get_Cellsize());
1035
1036 Slope = atan(sqrt(G*G + H*H));
1037 Aspect = G != 0. ? M_PI_180 + atan2(H, G) : H > 0. ? M_PI_270 : H < 0. ? M_PI_090 : -1.;
1038
1039 return( true );
1040 }
1041
1042 Slope = 0.;
1043 Aspect = -1.;
1044
1045 return( false );
1046 }
1047
1048 //---------------------------------------------------------
1049 /**
1050 * Calculates the gradient for the given world coordinate.
1051 * Calculation uses the formulas proposed by Zevenbergen & Thorne (1986).
1052 * Slope and aspect values are calculated in radians.
1053 * Aspect is zero for the North direction and increases clockwise.
1054 */
1055 //---------------------------------------------------------
Get_Gradient(const TSG_Point & p,double & Incline,double & Azimuth,TSG_Grid_Resampling Interpolation) const1056 bool CSG_Grid::Get_Gradient(const TSG_Point &p, double &Incline, double &Azimuth, TSG_Grid_Resampling Interpolation) const
1057 {
1058 return( Get_Gradient(p.x, p.y, Incline, Azimuth, Interpolation) );
1059 }
1060
1061
1062 ///////////////////////////////////////////////////////////
1063 // //
1064 // //
1065 // //
1066 ///////////////////////////////////////////////////////////
1067
1068 //---------------------------------------------------------
1069