1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // imagery_tools //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // landsat_acca.cpp //
14 // //
15 // Copyright (C) 2013 by //
16 // Benjamin Bechtel & Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "landsat_acca.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 #define NO_DEFINED 1
60 #define IS_SHADOW 2
61 #define IS_COLD_CLOUD 6
62 #define IS_WARM_CLOUD 9
63
64 //---------------------------------------------------------
65 #define LUT_SET_CLASS(id, name, color) { CSG_Table_Record *pR = pLUT->asTable()->Add_Record(); pR->Set_Value(0, color); pR->Set_Value(1, name); pR->Set_Value(3, id); pR->Set_Value(3, id); }
66
67 //---------------------------------------------------------
68 void acca_algorithm (CSG_Grid *pCloud, CSG_Grid *band[], int single_pass, int with_shadow, int cloud_signature);
69 void filter_holes (CSG_Grid *pGrid);
70
71
72 ///////////////////////////////////////////////////////////
73 // //
74 // //
75 // //
76 ///////////////////////////////////////////////////////////
77
78 //---------------------------------------------------------
CLandsat_ACCA(void)79 CLandsat_ACCA::CLandsat_ACCA(void)
80 {
81 //-----------------------------------------------------
82 Set_Name (_TL("Automated Cloud Cover Assessment"));
83
84 Set_Author ("B.Bechtel, O.Conrad (c) 2013");
85
86 Set_Description (_TW(
87 "Automated Cloud-Cover Assessment (ACCA) for Landsat TM/ETM+ imagery as proposed by Irish (2000). "
88 "This tool incorporates E.J. Tizado's GRASS GIS implementation (i.landsat.acca)."
89 ));
90
91 Add_Reference("Irish, R.R.", "2000",
92 "Landsat 7 Automatic Cloud Cover Assessment",
93 "In: Shen, S.S., Descour, M.R. [Eds.]: Algorithms for Multispectral, Hyperspectral, and Ultraspectral Imagery VI. Proceedings of SPIE, 4049: 348-355.",
94 SG_T("http://landsathandbook.gsfc.nasa.gov/pdfs/ACCA_Special_Issue_Final.pdf")
95 );
96
97 Add_Reference("Irish, R.R., Barker, J.L., Goward, S.N., Arvidson, T.", "2006",
98 "Characterization of the Landsat-7 ETM+ Automated Cloud-Cover Assessment (ACCA) Algorithm.",
99 "Photogrammetric Engineering and Remote Sensing vol. 72(10): 1179-1188.",
100 SG_T("https://pdfs.semanticscholar.org/3d20/f685ce83a82b294f0773768624d4166d105d.pdf")
101 );
102
103 Add_Reference("https://grass.osgeo.org/grass72/manuals/i.landsat.acca.html",
104 SG_T("GRASS i.landsat.acca")
105 );
106
107 //-----------------------------------------------------
108 Parameters.Add_Grid("", "BAND2", _TL("Green" ), _TL("Landsat band 2 reflectance"), PARAMETER_INPUT);
109 Parameters.Add_Grid("", "BAND3", _TL("Red" ), _TL("Landsat band 3 reflectance"), PARAMETER_INPUT);
110 Parameters.Add_Grid("", "BAND4", _TL("NIR" ), _TL("Landsat band 4 reflectance"), PARAMETER_INPUT);
111 Parameters.Add_Grid("", "BAND5", _TL("SWIR" ), _TL("Landsat band 5 reflectance"), PARAMETER_INPUT);
112 Parameters.Add_Grid("", "BAND6", _TL("Thermal"), _TL("Landsat band 6 temperature [Kelvin]"), PARAMETER_INPUT, false);
113
114 Parameters.Add_Grid("",
115 "CLOUD" , _TL("Cloud Cover"),
116 _TL(""),
117 PARAMETER_OUTPUT, true, SG_DATATYPE_Char
118 );
119
120 Parameters.Add_Bool("",
121 "FILTER" , _TL("Apply post-processing filter to remove small holes"),
122 _TL(""),
123 true
124 );
125
126 //-----------------------------------------------------
127 Parameters.Add_Node("",
128 "NODE_THRS" , _TL("Thresholds"),
129 _TL("")
130 );
131
132 Parameters.Add_Double("NODE_THRS",
133 "B56C" , _TL("B56 Composite (step 6)"),
134 _TL(""),
135 225.0
136 );
137
138 Parameters.Add_Double("NODE_THRS",
139 "B45R" , _TL("B45 Ratio: Desert detection (step 10)"),
140 _TL(""),
141 1.0
142 );
143
144 //-----------------------------------------------------
145 Parameters.Add_Node("",
146 "NODE_CLOUD", _TL("Cloud Settings"),
147 _TL("")
148 );
149
150 // Parameters.Add_Int("NODE_CLOUD",
151 // "HIST_N" , _TL("Number of classes in the cloud temperature histogram"),
152 // _TL(""),
153 // 100, 10, true
154 // );
155
156 Parameters.Add_Bool("NODE_CLOUD",
157 "CSIG" , _TL("Always use cloud signature (step 14)"),
158 _TL(""),
159 true
160 );
161
162 Parameters.Add_Bool("NODE_CLOUD",
163 "PASS2" , _TL("Bypass second-pass processing, and merge warm (not ambiguous) and cold clouds"),
164 _TL(""),
165 true
166 );
167
168 Parameters.Add_Bool("NODE_CLOUD",
169 "SHADOW" , _TL("Include a category for cloud shadows"),
170 _TL(""),
171 true
172 );
173 }
174
175
176 ///////////////////////////////////////////////////////////
177 // //
178 ///////////////////////////////////////////////////////////
179
180 //---------------------------------------------------------
On_Execute(void)181 bool CLandsat_ACCA::On_Execute(void)
182 {
183 CSG_Grid *pCloud, *pBand[5];
184
185 //-----------------------------------------------------
186 pBand[0] = Parameters("BAND2")->asGrid();
187 pBand[1] = Parameters("BAND3")->asGrid();
188 pBand[2] = Parameters("BAND4")->asGrid();
189 pBand[3] = Parameters("BAND5")->asGrid();
190 pBand[4] = Parameters("BAND6")->asGrid();
191
192 pCloud = Parameters("CLOUD")->asGrid();
193
194 //-----------------------------------------------------
195 CSG_Parameter *pLUT = DataObject_Get_Parameter(pCloud, "LUT");
196
197 if( pLUT && pLUT->asTable() )
198 {
199 pLUT->asTable()->Del_Records();
200
201 LUT_SET_CLASS(IS_SHADOW , _TL("Shadow" ), SG_COLOR_RED);
202 LUT_SET_CLASS(IS_COLD_CLOUD, _TL("Cold Cloud"), SG_COLOR_BLUE);
203 LUT_SET_CLASS(IS_WARM_CLOUD, _TL("Warm Cloud"), SG_COLOR_BLUE_LIGHT);
204
205 DataObject_Set_Parameter(pCloud, pLUT);
206
207 DataObject_Set_Parameter(pCloud, "COLORS_TYPE", 1); // Color Classification Type: Lookup Table
208 }
209
210 pCloud->Set_NoData_Value(0.);
211
212 //-----------------------------------------------------
213 // int hist_n = Parameters("HIST_N")->asInt();
214
215 //-----------------------------------------------------
216 acca_algorithm(pCloud, pBand,
217 Parameters("PASS2" )->asBool(),
218 Parameters("SHADOW")->asBool(),
219 Parameters("CSIG" )->asBool()
220 );
221
222 if( Parameters("FILTER")->asBool() )
223 {
224 filter_holes(pCloud);
225 }
226
227 //-----------------------------------------------------
228 return( true );
229 }
230
231
232 ///////////////////////////////////////////////////////////
233 // //
234 // //
235 // //
236 ///////////////////////////////////////////////////////////
237
238 /****************************************************************************
239 *
240 * TOOL: i.landsat.acca
241 *
242 * AUTHOR(S): E. Jorge Tizado - ej.tizado@unileon.es
243 *
244 * PURPOSE: Landsat TM/ETM+ Automatic Cloud Cover Assessment
245 * Adopted for GRASS 7 by Martin Landa <landa.martin gmail.com>
246 *
247 * COPYRIGHT: (C) 2008, 2010 by the GRASS Development Team
248 *
249 * This program is free software under the GNU General Public
250 * License (>=v2). Read the file COPYING that comes with GRASS
251 * for details.
252 *
253 *****************************************************************************/
254
255 #define SCALE 200.
256 #define K_BASE 230.
257
258 /* value and count */
259 #define TOTAL 0
260 #define WARM 1
261 #define COLD 2
262 #define SNOW 3
263 #define SOIL 4
264
265 /* signa */
266 #define COVER 1
267 #define SUM_COLD 0
268 #define SUM_WARM 1
269 #define KMEAN 2
270 #define KMAX 3
271 #define KMIN 4
272
273 /* re-use value */
274 #define KLOWER 0
275 #define KUPPER 1
276 #define MEAN 2
277 #define SKEW 3
278 #define DSTD 4
279
280
281 /**********************************************************
282 *
283 * Automatic Cloud Cover Assessment (ACCA): Irish 2000
284 *
285 **********************************************************/
286
287 /*--------------------------------------------------------
288 CONSTANTS
289 Usar esta forma para que via extern puedan modificarse
290 como opciones desde el programa main.
291 ---------------------------------------------------------*/
292
293 double th_1 = 0.08; /* Band 3 Brightness Threshold */
294 double th_1_b = 0.07;
295 double th_2[2] = { -0.25, 0.70 }; /* Normalized Snow Difference Index */
296 double th_2_b = 0.8;
297 double th_3 = 300.; /* Band 6 Temperature Threshold */
298 double th_4 = 225.; /* Band 5/6 Composite */
299 double th_4_b = 0.08;
300 double th_5 = 2.35; /* Band 4/3 Ratio */
301 double th_6 = 2.16248; /* Band 4/2 Ratio */
302 double th_7 = 1.0; /* Band 4/5 Ratio */ ;
303 double th_8 = 210.; /* Band 5/6 Composite */
304
305 //---------------------------------------------------------
306 const int hist_n = 100; /* interval of real data 100/hist_n */
307
308 //---------------------------------------------------------
309 #define G_message(s) SG_UI_Msg_Add(s, false)
310
311 void acca_first(CSG_Grid *pCloud, CSG_Grid *band[], int with_shadow, int count[], int cold[], int warm[], double stats[]);
312 void acca_second(CSG_Grid *pCloud, CSG_Grid *band, int review_warm, double upper, double lower);
313 int shadow_algorithm(double pixel[]);
314
315 void hist_put(double t, int hist[]);
316 double quantile(double q, int hist[]);
317 double moment(int n, int hist[], int k);
318
319
320 #define BAND2 0
321 #define BAND3 1
322 #define BAND4 2
323 #define BAND5 3
324 #define BAND6 4
325
326 #define NO_CLOUD 0
327 #define IS_CLOUD 1
328 #define COLD_CLOUD 30
329 #define WARM_CLOUD 50
330
331 //---------------------------------------------------------
acca_algorithm(CSG_Grid * pCloud,CSG_Grid * band[],int single_pass,int with_shadow,int cloud_signature)332 void acca_algorithm(CSG_Grid *pCloud, CSG_Grid *band[], int single_pass, int with_shadow, int cloud_signature)
333 {
334 int i, count[5], hist_cold[hist_n], hist_warm[hist_n], review_warm;
335 double max, value[5], signa[5], idesert, shift;
336
337 /* Reset variables ... */
338 for (i = 0; i < 5; i++) {
339 count[i] = 0;
340 value[i] = 0.;
341 }
342
343 for (i = 0; i < hist_n; i++) {
344 hist_cold[i] = hist_warm[i] = 0;
345 }
346
347 /* FIRST FILTER ... */
348 acca_first(pCloud, band, with_shadow, count, hist_cold, hist_warm, signa);
349 /* CATEGORIES: NO_DEFINED, WARM_CLOUD, COLD_CLOUD, NULL (= NO_CLOUD) */
350
351 value[WARM] = (double)count[WARM] / (double)count[TOTAL];
352 value[COLD] = (double)count[COLD] / (double)count[TOTAL];
353 value[SNOW] = (double)count[SNOW] / (double)count[TOTAL];
354 value[SOIL] = (double)count[SOIL] / (double)count[TOTAL];
355
356 value[0] = (double)(count[WARM] + count[COLD]);
357 idesert = (value[0] == 0. ? 0. : value[0] / ((double)count[SOIL]));
358
359 //-----------------------------------------------------
360 // BAND-6 CLOUD SIGNATURE DEVELOPMENT
361 if( idesert <= .5 || value[SNOW] > 0.01 )
362 { // Only the cold clouds are used if snow or desert soil is present
363 review_warm = 1;
364 }
365 else
366 { // The cold and warm clouds are combined and treated as a single population
367 review_warm = 0;
368
369 count[COLD] += count[WARM];
370 value[COLD] += value[WARM];
371 signa[SUM_COLD] += signa[SUM_WARM];
372
373 for(i=0; i<hist_n; i++)
374 hist_cold[i] += hist_warm[i];
375 }
376
377 signa[KMEAN] = SCALE * signa[SUM_COLD] / ((double)count[COLD]);
378 signa[COVER] = ((double)count[COLD]) / ((double)count[TOTAL]);
379
380 /* G_message(_TL("Preliminary scene analysis:"));
381 G_message(("* Desert index: %.2lf"), idesert);
382 G_message(("* Snow cover: %.2lf %%"), 100. * value[SNOW]);
383 G_message(("* Cloud cover: %.2lf %%"), 100. * signa[COVER]);
384 G_message(("* Temperature of clouds:"));
385 G_message(("** Maximum: %.2lf K"), signa[KMAX]);
386 G_message(("** Mean (%s cloud): %.2lf K"), (review_warm ? "cold" : "all"), signa[KMEAN]);
387 G_message(("** Minimum: %.2lf K"), signa[KMIN]);
388 /**/
389 /* WARNING: re-use of the variable 'value' with new meaning */
390
391 /* step 14 */
392
393 /* To correct Irish2006: idesert has to be bigger than 0.5 to start pass 2 processing (see Irish2000)
394 because then we have no desert condition (thanks to Matthias Eder, Germany) */
395 if( cloud_signature || (idesert > .5 && signa[COVER] > 0.004 && signa[KMEAN] < 295.) )
396 {
397 G_message(_TL("Histogram cloud signature:"));
398
399 value[MEAN] = quantile(0.5, hist_cold) + K_BASE;
400 value[DSTD] = sqrt(moment(2, hist_cold, 1));
401 value[SKEW] = moment(3, hist_cold, 3) / pow(value[DSTD], 3);
402
403 /* G_message(("* Mean temperature: %.2lf K"), value[MEAN]);
404 G_message(("* Standard deviation: %.2lf"), value[DSTD]);
405 G_message(("* Skewness: %.2lf"), value[SKEW]);
406 G_message(("* Histogram classes: %d"), hist_n);
407 /**/
408 shift = value[SKEW];
409 if (shift > 1.)
410 shift = 1.;
411 else if (shift < 0.)
412 shift = 0.;
413
414 max = quantile(0.9875, hist_cold) + K_BASE;
415 value[KUPPER] = quantile(0.975, hist_cold) + K_BASE;
416 value[KLOWER] = quantile(0.835, hist_cold) + K_BASE;
417
418 /* G_message(("* 98.75 percentile: %.2lf K"), max);
419 G_message(("* 97.50 percentile: %.2lf K"), value[KUPPER]);
420 G_message(("* 83.50 percentile: %.2lf K"), value[KLOWER]);
421 /**/
422 /* step 17 & 18 */
423 if (shift > 0.)
424 {
425 shift *= value[DSTD];
426
427 if ((value[KUPPER] + shift) > max)
428 {
429 if ((value[KLOWER] + shift) > max)
430 {
431 value[KLOWER] += (max - value[KUPPER]);
432 }
433 else
434 {
435 value[KLOWER] += shift;
436 }
437
438 value[KUPPER] = max;
439 }
440 else
441 {
442 value[KLOWER] += shift;
443 value[KUPPER] += shift;
444 }
445 }
446
447 /* G_message(("Maximum temperature:"));
448 G_message(("* Cold cloud: %.2lf K"), value[KUPPER]);
449 G_message(("* Warm cloud: %.2lf K"), value[KLOWER]);
450 /**/
451 }
452 else if( signa[KMEAN] < 295. )
453 { // Retained warm and cold clouds
454 G_message(_TL("Result: Scene with clouds"));
455 review_warm = 0;
456 value[KUPPER] = 0.;
457 value[KLOWER] = 0.;
458 }
459 else
460 { // Retained cold clouds
461 G_message(_TL("Result: Scene cloud free"));
462 review_warm = 1;
463 value[KUPPER] = 0.;
464 value[KLOWER] = 0.;
465 }
466
467 //-----------------------------------------------------
468 // SECOND FILTER ...
469
470 // By-pass two processing but it retains warm and cold clouds
471 if( single_pass != 0 )
472 {
473 review_warm = -1;
474 value[KUPPER] = 0.;
475 value[KLOWER] = 0.;
476 }
477
478 // CATEGORIES: IS_WARM_CLOUD, IS_COLD_CLOUD, IS_SHADOW, NULL (= NO_CLOUD)
479 acca_second(pCloud, band[BAND6], review_warm, value[KUPPER], value[KLOWER]);
480
481 //-----------------------------------------------------
482 return;
483 }
484
485 //---------------------------------------------------------
acca_first(CSG_Grid * pCloud,CSG_Grid * band[],int with_shadow,int count[],int cold[],int warm[],double stats[])486 void acca_first(CSG_Grid *pCloud, CSG_Grid *band[], int with_shadow, int count[], int cold[], int warm[], double stats[])
487 {
488 double nsdi, rat56;
489
490 /* Creation of output file */
491 /* ----- ----- */
492 SG_UI_Msg_Add_Execution(_TL("Processing first pass..."), true);
493
494 stats[SUM_COLD] = 0.;
495 stats[SUM_WARM] = 0.;
496 stats[KMAX] = 0.;
497 stats[KMIN] = 10000.;
498
499 for(int y=0; y<pCloud->Get_NY() && SG_UI_Process_Set_Progress(y, pCloud->Get_NY()); y++)
500 {
501 for(int x=0; x<pCloud->Get_NX(); x++)
502 {
503 char code = NO_DEFINED;
504 double pixel[5];
505
506 for(int i=BAND2; i<=BAND6; i++) // Null when null pixel in any band
507 {
508 if( pCloud->Get_System() == band[i]->Get_System() )
509 {
510 if( band[i]->is_NoData(x, y) )
511 {
512 code = NO_CLOUD;
513 break;
514 }
515
516 pixel[i] = band[i]->asDouble(x, y);
517 }
518 else if( !band[i]->Get_Value(pCloud->Get_System().Get_Grid_to_World(x, y), pixel[i]) )
519 {
520 code = NO_CLOUD;
521 break;
522 }
523 }
524
525 /* Determina los pixeles de sombras */
526 if( code == NO_DEFINED && with_shadow )
527 {
528 code = shadow_algorithm(pixel);
529 }
530
531 /* Analiza el valor de los pixeles no definidos */
532 if (code == NO_DEFINED)
533 {
534 code = NO_CLOUD;
535 count[TOTAL]++;
536 nsdi = (pixel[BAND2] - pixel[BAND5]) / (pixel[BAND2] + pixel[BAND5]);
537
538 /* ----------------------------------------------------- */
539 /* step 1. Brightness Threshold: Eliminates dark images */
540 if (pixel[BAND3] > th_1)
541 {
542 /* step 3. Normalized Snow Difference Index: Eliminates many types of snow */
543 if (nsdi > th_2[0] && nsdi < th_2[1])
544 {
545 /* step 5. Temperature Threshold: Eliminates warm image features */
546 if (pixel[BAND6] < th_3)
547 {
548 rat56 = (1. - pixel[BAND5]) * pixel[BAND6];
549 /* step 6. Band 5/6 Composite: Eliminates numerous categories including ice */
550 if (rat56 < th_4)
551 {
552 /* step 8. Eliminates growing vegetation */
553 if ((pixel[BAND4] / pixel[BAND3]) < th_5)
554 {
555 /* step 9. Eliminates senescing vegetation */
556 if ((pixel[BAND4] / pixel[BAND2]) < th_6)
557 {
558 /* step 10. Eliminates rocks and desert */
559 count[SOIL]++;
560
561 if ((pixel[BAND4] / pixel[BAND5]) > th_7)
562 {
563 /* step 11. Distinguishes warm clouds from cold clouds */
564 if (rat56 < th_8)
565 {
566 code = COLD_CLOUD;
567 count[COLD]++;
568 /* for statistic */
569 stats[SUM_COLD] += (pixel[BAND6] / SCALE);
570 hist_put(pixel[BAND6] - K_BASE, cold);
571 }
572 else
573 {
574 code = WARM_CLOUD;
575 count[WARM]++;
576 /* for statistic */
577 stats[SUM_WARM] += (pixel[BAND6] / SCALE);
578 hist_put(pixel[BAND6] - K_BASE, warm);
579 }
580
581 if (pixel[BAND6] > stats[KMAX]) stats[KMAX] = pixel[BAND6];
582 if (pixel[BAND6] < stats[KMIN]) stats[KMIN] = pixel[BAND6];
583 }
584 else
585 {
586 code = NO_DEFINED;
587 }
588 }
589 else
590 {
591 code = NO_DEFINED;
592 count[SOIL]++;
593 }
594 }
595 else
596 {
597 code = NO_DEFINED;
598 }
599 }
600 else
601 {
602 /* step 7 */
603 code = (pixel[BAND5] < th_4_b) ? NO_CLOUD : NO_DEFINED;
604 }
605 }
606 else
607 {
608 code = NO_CLOUD;
609 }
610 }
611 else
612 {
613 /* step 3 */
614 code = NO_CLOUD;
615
616 if (nsdi > th_2_b)
617 count[SNOW]++;
618 }
619 }
620 else
621 {
622 /* step 2 */
623 code = (pixel[BAND3] < th_1_b) ? NO_CLOUD : NO_DEFINED;
624 }
625 }
626
627 //---------------------------------------------
628 if (code == NO_CLOUD)
629 {
630 pCloud->Set_Value(x, y, 1);
631 }
632 else
633 {
634 pCloud->Set_Value(x, y, code);
635 }
636 }
637 }
638
639 return;
640 }
641
642 //---------------------------------------------------------
acca_second(CSG_Grid * pCloud,CSG_Grid * pThermal,int review_warm,double upper,double lower)643 void acca_second(CSG_Grid *pCloud, CSG_Grid *pThermal, int review_warm, double upper, double lower)
644 {
645 SG_UI_Process_Set_Text(upper == 0.0
646 ? _TL("Removing ambiguous pixels...")
647 : _TL("Pass two processing...")
648 );
649
650 //-----------------------------------------------------
651 for(int y=0; y<pCloud->Get_NY() && SG_UI_Process_Set_Progress(y, pCloud->Get_NY()); y++)
652 {
653 double p_y = pCloud->Get_YMin() + y * pCloud->Get_Cellsize();
654
655 #pragma omp parallel for
656 for(int x=0; x<pCloud->Get_NX(); x++)
657 {
658 if( !pCloud->is_NoData(x, y) )
659 {
660 int code = pCloud->asInt(x, y);
661
662 if( code == NO_DEFINED || (code == WARM_CLOUD && review_warm == 1) ) // Resolve ambiguous pixels
663 {
664 double t, p_x = pCloud->Get_XMin() + x * pCloud->Get_Cellsize();
665
666 if( !pThermal->Get_Value(p_x, p_y, t) || t > upper )
667 {
668 pCloud->Set_NoData(x, y);
669 }
670 else
671 {
672 pCloud->Set_Value(x, y, t < lower ? IS_WARM_CLOUD : IS_COLD_CLOUD);
673 }
674 }
675 else if( code == COLD_CLOUD || code == WARM_CLOUD ) // Join warm (not ambiguous) and cold clouds
676 {
677 pCloud->Set_Value(x, y, (code == WARM_CLOUD && review_warm == 0) ? IS_WARM_CLOUD : IS_COLD_CLOUD);
678 }
679 else
680 {
681 pCloud->Set_Value(x, y, IS_SHADOW);
682 }
683 }
684 }
685 }
686
687 //-----------------------------------------------------
688 return;
689 }
690
691
692 ///////////////////////////////////////////////////////////
693 // //
694 // Cloud shadows //
695 // //
696 ///////////////////////////////////////////////////////////
697
698 //---------------------------------------------------------
shadow_algorithm(double pixel[])699 int shadow_algorithm(double pixel[])
700 {
701 // I think this filter is better but not in any paper
702 if( pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. && pixel[BAND4] / pixel[BAND2] > 1.
703 && (pixel[BAND3] - pixel[BAND5]) / (pixel[BAND3] + pixel[BAND5]) < 0.10 )
704 // if( pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. && pixel[BAND4] / pixel[BAND2] > 1. )
705 {
706 return IS_SHADOW;
707 }
708
709 return NO_DEFINED;
710 }
711
712
713 ///////////////////////////////////////////////////////////
714 // //
715 // HISTOGRAM ANALYSIS //
716 // //
717 ///////////////////////////////////////////////////////////
718
719 //---------------------------------------------------------
720 // Define un factor de escala = hist_n/100 con objeto
721 // de dividir el entero 1 por 100/hist_n partes y
722 // aumentar la precision.
723 //
724 // Afecta al almacenamiento en el histograma pero
725 // modifica el calculo de quantiles y momentos.
726 //---------------------------------------------------------
727
728 //---------------------------------------------------------
hist_put(double t,int hist[])729 void hist_put(double t, int hist[])
730 {
731 int i = (int)(t * ((double)hist_n / 100.)); // scale factor
732
733 if( i < 1 )
734 i = 1;
735 else if (i > hist_n)
736 i = hist_n;
737
738 hist[i - 1]++;
739 }
740
741 //---------------------------------------------------------
742 /* histogram moment */
moment(int n,int hist[],int k)743 double moment(int n, int hist[], int k)
744 {
745 int i, total;
746 double value, mean;
747
748 for(i=0, total=0, mean=0; i<hist_n; i++)
749 {
750 total += hist[i];
751 mean += (double)(i * hist[i]);
752 }
753
754 mean /= ((double)total); /* histogram mean */
755
756 for(i=0, value=0.; i<hist_n; i++)
757 {
758 value += (pow((i - mean), n) * ((double)hist[i]));
759 }
760
761 value /= (double)(total); // k = 0; // ???!!!
762 // value /= (double)(total - k);
763
764 return( value / pow((double)hist_n / 100., n) );
765 }
766
767 //---------------------------------------------------------
768 /* Real data quantile */
quantile(double q,int hist[])769 double quantile(double q, int hist[])
770 {
771 int i, total;
772 double value, qmax, qmin;
773
774 for(i=0, total=0; i<hist_n; i++)
775 {
776 total += hist[i];
777 }
778
779 for(i=hist_n-1, value=0, qmax=1.; i>=0; i--)
780 {
781 qmin = qmax - (double)hist[i] / (double)total;
782
783 if( q >= qmin )
784 {
785 value = (q - qmin) / (qmax - qmin) + (i - 1);
786 break;
787 }
788
789 qmax = qmin;
790 }
791
792 /* remove scale factor */
793 return (value / ((double)hist_n / 100.));
794 }
795
796
797 ///////////////////////////////////////////////////////////
798 // //
799 // FILTER HOLES OF CLOUDS //
800 // //
801 ///////////////////////////////////////////////////////////
802
803 //---------------------------------------------------------
804 // This a >=50% filter of 3x3
805 // if >= 50% vecinos cloud then pixel set to cloud
806 //---------------------------------------------------------
filter_holes(CSG_Grid * pCloud)807 void filter_holes(CSG_Grid *pCloud)
808 {
809 if( pCloud->Get_NY() < 3 || pCloud->Get_NX() < 3 )
810 return;
811
812 SG_UI_Process_Set_Text(_TL("Filling small holes in clouds..."));
813
814 CSG_Grid Cloud(*pCloud);
815
816 //-----------------------------------------------------
817 for(int y=0; y<pCloud->Get_NY() && SG_UI_Process_Set_Progress(y, pCloud->Get_NY()); y++)
818 {
819 #pragma omp parallel for
820 for(int x=0; x<pCloud->Get_NX(); x++)
821 {
822 int z = Cloud.asInt(x, y);
823
824 if( z == 0 )
825 {
826 int cold, warm, shadow, nulo;
827
828 cold = warm = shadow = nulo = 0;
829
830 for(int i=0; i<8; i++)
831 {
832 int ix = pCloud->Get_System().Get_xTo(i, x);
833 int iy = pCloud->Get_System().Get_yTo(i, y);
834
835 switch( Cloud.is_InGrid(ix, iy) ? Cloud.asInt(ix, iy) : -1 )
836 {
837 case IS_COLD_CLOUD: cold ++; break;
838 case IS_WARM_CLOUD: warm ++; break;
839 case IS_SHADOW: shadow++; break;
840 default: nulo ++; break;
841 }
842 }
843
844 int lim = (cold + warm + shadow + nulo) / 2;
845
846 // Entra pixel[0] = 0
847 if( nulo < lim )
848 {
849 if( shadow >= (cold + warm) )
850 z = IS_SHADOW;
851 else
852 z = (warm > cold) ? IS_WARM_CLOUD : IS_COLD_CLOUD;
853 }
854 }
855
856 if( z != 0 )
857 {
858 pCloud->Set_Value(x, y, z);
859 }
860 else
861 {
862 pCloud->Set_NoData(x, y);
863 }
864 }
865 }
866
867 //-----------------------------------------------------
868 return;
869 }
870
871
872 ///////////////////////////////////////////////////////////
873 // //
874 // //
875 // //
876 ///////////////////////////////////////////////////////////
877
878 //---------------------------------------------------------
879