[GRASS-SVN] r43784 - grass-addons/imagery/i.landsat.acca
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Oct 4 16:05:59 EDT 2010
Author: ejtizado
Date: 2010-10-04 20:05:59 +0000 (Mon, 04 Oct 2010)
New Revision: 43784
Modified:
grass-addons/imagery/i.landsat.acca/algorithm.c
grass-addons/imagery/i.landsat.acca/description.html
grass-addons/imagery/i.landsat.acca/main.c
grass-addons/imagery/i.landsat.acca/tools.c
Log:
shift corrected
Modified: grass-addons/imagery/i.landsat.acca/algorithm.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/algorithm.c 2010-10-04 19:09:42 UTC (rev 43783)
+++ grass-addons/imagery/i.landsat.acca/algorithm.c 2010-10-04 20:05:59 UTC (rev 43784)
@@ -1,3 +1,13 @@
+/* File: algorithm.c
+ *
+ * AUTHOR: E. Jorge Tizado, Spain 2010
+ *
+ * COPYRIGHT: (c) 2010 E. Jorge Tizado
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ */
+
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
@@ -9,7 +19,9 @@
#include "local_proto.h"
#define SCALE 200.
+#define K_BASE 230.
+
/* value and count */
#define TOTAL 0
#define WARM 1
@@ -32,18 +44,7 @@
#define SKEW 3
#define DSTD 4
-/*
- Con Landsat-7 funciona bien pero con Landsat-5 falla
- la primera pasada metiendo terreno desertico como nube fría
- y luego se equivoca en resolver los ambiguos en áreas
- con mucha humedad ambiental.
- Habría que reajustar las constantes para Landsat-5
- */
-
-
-
-
/**********************************************************
*
* Automatic Cloud Cover Assessment (ACCA): Irish 2000
@@ -58,7 +59,7 @@
double th_1 = 0.08; /* Band 3 Brightness Threshold */
double th_1_b = 0.07;
-double th_2[] = { -0.25, 0.70 }; /* Normalized Snow Difference Index */
+double th_2[] = { -0.25, 0.70 }; /* Normalized Snow Difference Index */
double th_2_b = 0.8;
double th_3 = 300.; /* Band 6 Temperature Threshold */
@@ -66,13 +67,11 @@
double th_4_b = 0.08;
double th_5 = 2.35; /* Band 4/3 Ratio */
double th_6 = 2.16248; /* Band 4/2 Ratio */
-double th_7 = 1.0; /* Band 4/5 Ratio */ ;
+double th_7 = 1.0; /* Band 4/5 Ratio */ ;
double th_8 = 210.; /* Band 5/6 Composite */
extern int hist_n;
-#define K_BASE 240.
-
void acca_algorithm(int verbose, Gfile * out, Gfile band[],
int single_pass, int with_shadow, int cloud_signature)
{
@@ -99,8 +98,7 @@
value[SOIL] = (double)count[SOIL] / (double)count[TOTAL];
value[0] = (double)(count[WARM] + count[COLD]);
- idesert =
- (value[0] == 0. ? 0. : 1. / (1. + ((double)count[SOIL]) / value[0]));
+ idesert = (value[0] == 0. ? 0. : value[0] / ((double)count[SOIL]));
/* BAND-6 CLOUD SIGNATURE DEVELOPMENT */
if (idesert <= .5 || value[SNOW] > 0.01) {
@@ -119,58 +117,66 @@
hist_cold[i] += hist_warm[i];
}
- signa[KMEAN] = (signa[SUM_COLD] / (double)count[COLD]) * SCALE;
- signa[COVER] = (double)count[COLD] / (double)count[TOTAL];
+ signa[KMEAN] = SCALE * signa[SUM_COLD] / ((double)count[COLD]);
+ signa[COVER] = ((double)count[COLD]) / ((double)count[TOTAL]);
- G_message("> PRELIMINARY SCENE ANALYSIS");
- G_message("> Desert index: %.3lf", idesert);
- G_message("> Snow cover: %.3lf %%", 100. * value[SNOW]);
- G_message("> Cloud cover: %.3lf %%", 100. * signa[COVER]);
- G_message("> Temperature of clouds");
- G_message("> Maximum: %.2lf K", signa[KMAX]);
- G_message("> Mean (%s cloud) : %.2lf K",
- (review_warm ? "cold" : "all"), signa[KMEAN]);
- G_message("> Minimum: %.2lf K", signa[KMIN]);
+ fprintf(stdout, " PRELIMINARY SCENE ANALYSIS\n");
+ fprintf(stdout, " Desert index: %.2lf\n", idesert);
+ fprintf(stdout, " Snow cover: %.2lf %%\n", 100. * value[SNOW]);
+ fprintf(stdout, " Cloud cover: %.2lf %%\n", 100. * signa[COVER]);
+ fprintf(stdout, " Temperature of clouds\n");
+ fprintf(stdout, " Maximum: %.2lf K\n", signa[KMAX]);
+ fprintf(stdout, " Mean (%s cloud) : %.2lf K\n",
+ (review_warm ? "cold" : "all"), signa[KMEAN]);
+ fprintf(stdout, " Minimum: %.2lf K\n", signa[KMIN]);
/* WARNING: re-use of the variable 'value' with new meaning */
/* step 14 */
if (cloud_signature ||
(idesert <= .5 && signa[COVER] > 0.004 && signa[KMEAN] < 295.)) {
+ fprintf(stdout, " HISTOGRAM CLOUD SIGNATURE\n");
+
value[MEAN] = quantile(0.5, hist_cold) + K_BASE;
value[DSTD] = sqrt(moment(2, hist_cold, 1));
value[SKEW] = moment(3, hist_cold, 3) / pow(value[DSTD], 3);
+ fprintf(stdout, " Mean temperature: %.2lf K\n", value[MEAN]);
+ fprintf(stdout, " Standard deviation: %.2lf\n", value[DSTD]);
+ fprintf(stdout, " Skewness: %.2lf\n", value[SKEW]);
+ fprintf(stdout, " Histogram classes: %d\n", hist_n);
+
shift = value[SKEW];
if (shift > 1.)
shift = 1.;
else if (shift < 0.)
shift = 0.;
- shift *= value[DSTD];
-
max = quantile(0.9875, hist_cold) + K_BASE;
value[KUPPER] = quantile(0.975, hist_cold) + K_BASE;
value[KLOWER] = quantile(0.835, hist_cold) + K_BASE;
+
+ fprintf(stdout, " 98.75 percentile: %.2lf K\n", max);
+ fprintf(stdout, " 97.50 percentile: %.2lf K\n", value[KUPPER]);
+ fprintf(stdout, " 83.50 percentile: %.2lf K\n", value[KLOWER]);
+
/* step 17 & 18 */
if (shift > 0.) {
- if ((value[KUPPER] * shift) > max) {
- value[KLOWER] += (value[KUPPER] * (shift - 1.));
+ shift *= value[DSTD];
+
+ if ((value[KUPPER] + shift) > max) {
+ value[KLOWER] += (max - value[KUPPER]);
value[KUPPER] = max;
}
else {
- value[KLOWER] *= shift;
- value[KUPPER] *= shift;
+ value[KLOWER] += shift;
+ value[KUPPER] += shift;
}
}
- G_message("> HISTOGRAM CLOUD SIGNATURE");
- G_message("> Histogram classes: %d", hist_n);
- G_message("> Mean temperature: %.2lf K", value[MEAN]);
- G_message("> Standard deviation: %.2lf", value[DSTD]);
- G_message("> Skewness: %.2lf", value[SKEW]);
- G_message("> Cold cloud: maximun %.2lf K", value[KUPPER]);
- G_message("> Warn cloud: maximun %.2lf K", value[KLOWER]);
+ fprintf(stdout, " Maximum temperature\n");
+ fprintf(stdout, " Cold cloud: %.2lf K\n", value[KUPPER]);
+ fprintf(stdout, " Warn cloud: %.2lf K\n", value[KLOWER]);
}
else {
if (signa[KMEAN] < 295.) {
@@ -256,45 +262,50 @@
nsdi = (pixel[BAND2] - pixel[BAND5]) /
(pixel[BAND2] + pixel[BAND5]);
/* ----------------------------------------------------- */
- /* Brightness Threshold: Eliminates dark images */
+ /* step 1. Brightness Threshold: Eliminates dark images */
if (pixel[BAND3] > th_1) {
- /* Normalized Snow Difference Index: Eliminates many types of snow */
+ /* step 3. Normalized Snow Difference Index: Eliminates many types of snow */
if (nsdi > th_2[0] && nsdi < th_2[1]) {
- /* Temperature Threshold: Eliminates warm image features */
+ /* step 5. Temperature Threshold: Eliminates warm image features */
if (pixel[BAND6] < th_3) {
rat56 = (1. - pixel[BAND5]) * pixel[BAND6];
- /* Band 5/6 Composite: Eliminates numerous categories including ice */
+ /* step 6. Band 5/6 Composite: Eliminates numerous categories including ice */
if (rat56 < th_4) {
- /* Eliminates growing vegetation */
- /* Eliminates senescing vegetation */
- if (pixel[BAND4] / pixel[BAND3] < th_5 &&
- pixel[BAND4] / pixel[BAND2] < th_6) {
- rat45 = pixel[BAND4] / pixel[BAND5];
- /* Eliminates rocks and desert */
- if (rat45 > th_7) {
- /* Distinguishes warm clouds from cold clouds */
- if (rat56 < th_8) {
- code = COLD_CLOUD;
- count[COLD]++;
- /* for statistic */
- stats[SUM_COLD] +=
- (pixel[BAND6] / SCALE);
- hist_put(pixel[BAND6] - K_BASE,
- cold);
+ /* step 8. Eliminates growing vegetation */
+ if ((pixel[BAND4] / pixel[BAND3]) < th_5) {
+ /* step 9. Eliminates senescing vegetation */
+ if ((pixel[BAND4] / pixel[BAND2]) < th_6) {
+ /* step 10. Eliminates rocks and desert */
+ count[SOIL]++;
+ if ((pixel[BAND4] / pixel[BAND5]) >
+ th_7) {
+ /* step 11. Distinguishes warm clouds from cold clouds */
+ if (rat56 < th_8) {
+ code = COLD_CLOUD;
+ count[COLD]++;
+ /* for statistic */
+ stats[SUM_COLD] +=
+ (pixel[BAND6] / SCALE);
+ hist_put(pixel[BAND6] -
+ K_BASE, cold);
+ }
+ else {
+ code = WARM_CLOUD;
+ count[WARM]++;
+ /* for statistic */
+ stats[SUM_WARM] +=
+ (pixel[BAND6] / SCALE);
+ hist_put(pixel[BAND6] -
+ K_BASE, warm);
+ }
+ if (pixel[BAND6] > stats[KMAX])
+ stats[KMAX] = pixel[BAND6];
+ if (pixel[BAND6] < stats[KMIN])
+ stats[KMIN] = pixel[BAND6];
}
else {
- code = WARM_CLOUD;
- count[WARM]++;
- /* for statistic */
- stats[SUM_WARM] +=
- (pixel[BAND6] / SCALE);
- hist_put(pixel[BAND6] - K_BASE,
- warm);
+ code = NO_DEFINED;
}
- if (pixel[BAND6] > stats[KMAX])
- stats[KMAX] = pixel[BAND6];
- if (pixel[BAND6] < stats[KMIN])
- stats[KMIN] = pixel[BAND6];
}
else {
code = NO_DEFINED;
@@ -306,6 +317,7 @@
}
}
else {
+ /* step 7 */
code =
(pixel[BAND5] <
th_4_b) ? NO_CLOUD : NO_DEFINED;
@@ -316,12 +328,14 @@
}
}
else {
+ /* step 3 */
code = NO_CLOUD;
if (nsdi > th_2_b)
count[SNOW]++;
}
}
else {
+ /* step 2 */
code = (pixel[BAND3] < th_1_b) ? NO_CLOUD : NO_DEFINED;
}
/* ----------------------------------------------------- */
@@ -459,12 +473,12 @@
/* I think this filter is better but not in any paper */
if (pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. &&
pixel[BAND4] / pixel[BAND2] > 1. &&
- (pixel[BAND3] - pixel[BAND5]) / (pixel[BAND3] + pixel[BAND5]) <
- 0.10) {
+ (pixel[BAND3] - pixel[BAND5]) / (pixel[BAND3] + pixel[BAND5]) < 0.10)
/*
if (pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. &&
- pixel[BAND4] / pixel[BAND2] > 1.) {
+ pixel[BAND4] / pixel[BAND2] > 1.)
*/
+ {
return IS_SHADOW;
}
Modified: grass-addons/imagery/i.landsat.acca/description.html
===================================================================
--- grass-addons/imagery/i.landsat.acca/description.html 2010-10-04 19:09:42 UTC (rev 43783)
+++ grass-addons/imagery/i.landsat.acca/description.html 2010-10-04 20:05:59 UTC (rev 43784)
@@ -23,13 +23,12 @@
Run the standard ACCA algorithm with filling of small cloud holes (the <b>-f</b> flag):
<p>
-<!-- maybe include i.landsat.toar step as well here to reinforce the point? -->
-
With per-band reflectance raster maps named <tt>226_62.toar.1,
226_62.toar.2, </tt> [...] and LANDSAT-7 thermal band <tt>226_62.toar.61</tt>,
outputing to a new raster map named <tt>226_62.acca</tt>:
<div class="code"><pre>
+i.landsat.toar sensor=7 gain=HHHLHLHHL date=2003-04-07 product_date=2008-11-27 band_prefix=226_62 solar_elevation=49.51654
i.landsat.acca -f band_prefix=226_62.toar output=226_62.acca
</pre></div>
Modified: grass-addons/imagery/i.landsat.acca/main.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/main.c 2010-10-04 19:09:42 UTC (rev 43783)
+++ grass-addons/imagery/i.landsat.acca/main.c 2010-10-04 20:05:59 UTC (rev 43784)
@@ -94,7 +94,7 @@
struct GModule *module;
int i, verbose = 1;
- struct Option *input, *output, *hist, *b56c;
+ struct Option *input, *output, *hist, *b56c, *b45r;
struct Flag *shadow, *filter, *sat5, *pass2, *csig;
char *in_name, *out_name;
struct Categories cats;
@@ -124,9 +124,16 @@
b56c->key = "b56composite";
b56c->type = TYPE_DOUBLE;
b56c->required = NO;
- b56c->description = _("Value for step 6: B56composite");
+ b56c->description = _("B56composite (step 6)");
b56c->answer = "225.";
+ b45r = G_define_option();
+ b45r->key = "b45ratio";
+ b45r->type = TYPE_DOUBLE;
+ b45r->required = NO;
+ b45r->description = _("B45ratio: Desert detection (step 10)");
+ b45r->answer = "1.";
+
hist = G_define_option();
hist->key = "histogram";
hist->type = TYPE_INTEGER;
@@ -161,7 +168,6 @@
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
-
/* stores OPTIONS and FLAGS to variables */
hist_n = atoi(hist->answer);
@@ -187,6 +193,7 @@
/* --------------------------------------- */
th_4 = atof(b56c->answer);
+ th_7 = atof(b45r->answer);
acca_algorithm(verbose, &out, band, pass2->answer, shadow->answer,
csig->answer);
Modified: grass-addons/imagery/i.landsat.acca/tools.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/tools.c 2010-10-04 19:09:42 UTC (rev 43783)
+++ grass-addons/imagery/i.landsat.acca/tools.c 2010-10-04 20:05:59 UTC (rev 43784)
@@ -42,7 +42,7 @@
double moment(int n, int hist[], int k)
{
int i, total;
- double value, mean, cte;
+ double value, mean;
k = 0;
@@ -52,15 +52,15 @@
total += hist[i];
mean += (double)(i * hist[i]);
}
- mean /= (double)total; /* histogram mean */
+ mean /= ((double)total); /* histogram mean */
value = 0.;
- cte = 100. / ((double)hist_n * (double)(total - k));
for (i = 0; i < hist_n; i++) {
- value += (pow((i - mean), n) * (double)hist[i] * cte);
+ value += (pow((i - mean), n) * ((double)hist[i]));
}
+ value /= (double)(total - k);
- return value;
+ return (value / pow((double)hist_n / 100., n) );
}
/* Real data quantile */
More information about the grass-commit
mailing list