[GRASS-SVN] r43723 - grass-addons/imagery/i.landsat.acca
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 28 08:12:02 EDT 2010
Author: ejtizado
Date: 2010-09-28 12:12:02 +0000 (Tue, 28 Sep 2010)
New Revision: 43723
Modified:
grass-addons/imagery/i.landsat.acca/algorithm.c
grass-addons/imagery/i.landsat.acca/description.html
grass-addons/imagery/i.landsat.acca/local_proto.h
grass-addons/imagery/i.landsat.acca/main.c
grass-addons/imagery/i.landsat.acca/tools.c
Log:
G_message
Modified: grass-addons/imagery/i.landsat.acca/algorithm.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/algorithm.c 2010-09-28 01:15:57 UTC (rev 43722)
+++ grass-addons/imagery/i.landsat.acca/algorithm.c 2010-09-28 12:12:02 UTC (rev 43723)
@@ -33,13 +33,13 @@
#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.
+ 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
-*/
+ Habría que reajustar las constantes para Landsat-5
+ */
@@ -56,407 +56,389 @@
como opciones desde el programa main.
---------------------------------------------------------*/
-double th_1 = 0.08; /* Band 3 Brightness Thershold */
+double th_1 = 0.08; /* Band 3 Brightness Thershold */
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 Thershold */
-double th_4 = 225.; /* Band 5/6 Composite */
+double th_3 = 300.; /* Band 6 Temperature Thershold */
+double th_4 = 225.; /* Band 5/6 Composite */
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_8 = 210.; /* Band 5/6 Composite */
+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_8 = 210.; /* Band 5/6 Composite */
extern int hist_n;
#define K_BASE 240.
void acca_algorithm(int verbose, Gfile * out, Gfile band[],
- int two_pass, int with_shadow)
+ int two_pass, int with_shadow)
{
- int i, count[5], hist_cold[hist_n], hist_warm[hist_n];
- double x, value[5], signa[5], idesert, warm_ambiguous, shift;
+ int i, count[5], hist_cold[hist_n], hist_warm[hist_n];
+ double x, value[5], signa[5], idesert, warm_ambiguous, shift;
- /* Reset variables ... */
- for (i = 0; i < 5; i++)
- {
- count[i] = 0;
- value[i] = 0.;
- }
- for (i = 0; i < hist_n; i++)
- {
- hist_cold[i] = hist_warm[i] = 0;
- }
+ /* Reset variables ... */
+ for (i = 0; i < 5; i++) {
+ count[i] = 0;
+ value[i] = 0.;
+ }
+ for (i = 0; i < hist_n; i++) {
+ hist_cold[i] = hist_warm[i] = 0;
+ }
- /* FIRST FILTER ... */
- acca_first(verbose, out, band, with_shadow,
- count, hist_cold, hist_warm, signa);
- /* CATEGORIES: NO_DEFINED, WARM_CLOUD, COLD_CLOUD, NULL (= NO_CLOUD) */
+ /* FIRST FILTER ... */
+ acca_first(verbose, out, band, with_shadow,
+ count, hist_cold, hist_warm, signa);
+ /* CATEGORIES: NO_DEFINED, WARM_CLOUD, COLD_CLOUD, NULL (= NO_CLOUD) */
- value[WARM] = (double)count[WARM] / (double)count[TOTAL];
- value[COLD] = (double)count[COLD] / (double)count[TOTAL];
- value[SNOW] = (double)count[SNOW] / (double)count[TOTAL];
- value[SOIL] = (double)count[SOIL] / (double)count[TOTAL];
+ value[WARM] = (double)count[WARM] / (double)count[TOTAL];
+ value[COLD] = (double)count[COLD] / (double)count[TOTAL];
+ value[SNOW] = (double)count[SNOW] / (double)count[TOTAL];
+ 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]);
+ value[0] = (double)(count[WARM] + count[COLD]);
+ idesert =
+ (value[0] == 0. ? 0. : 1. / (1. + (double)count[SOIL]) / value[0]);
- /* BAND-6 CLOUD SIGNATURE DEVELOPMENT */
- if (value[SNOW] > 0.01)
- {
- /* Only the cold clouds are used
- if snow or desert soil is present */
- warm_ambiguous = 1;
- }
- else
- {
- /* The cold and warm clouds are combined
- and treated as a single population */
- warm_ambiguous = 0;
- count[COLD] += count[WARM];
- value[COLD] += value[WARM];
- signa[SUM_COLD] += signa[SUM_WARM];
- for (i = 0; i < hist_n; i++) hist_cold[i] += hist_warm[i];
- }
+ /* BAND-6 CLOUD SIGNATURE DEVELOPMENT */
+ if (value[SNOW] > 0.01)
+ // if (idesert <= .5 || value[SNOW] > 0.01)
+ {
+ /* Only the cold clouds are used
+ if snow or desert soil is present */
+ warm_ambiguous = 1;
+ }
+ else {
+ /* The cold and warm clouds are combined
+ and treated as a single population */
+ warm_ambiguous = 0;
+ count[COLD] += count[WARM];
+ value[COLD] += value[WARM];
+ signa[SUM_COLD] += signa[SUM_WARM];
+ for (i = 0; i < hist_n; i++)
+ hist_cold[i] += hist_warm[i];
+ }
- signa[KMEAN] = (signa[SUM_COLD] / (double)count[COLD]) * SCALE;
- /* signa[KMEAN] = quantile( 0.50, hist_cold ) :::::: median */
- signa[COVER] = (double)count[COLD] / (double)count[TOTAL];
+ signa[KMEAN] = (signa[SUM_COLD] / (double)count[COLD]) * SCALE;
+ signa[COVER] = (double)count[COLD] / (double)count[TOTAL];
- fprintf(stdout, " PRELIMINARY SCENE ANALYSIS\n");
- fprintf(stdout, " Desert index: %.3lf\n", idesert);
- fprintf(stdout, " Snow cover : %.3lf %%\n", 100.*value[SNOW]);
- fprintf(stdout, " Cloud cover : %.3lf %%\n", 100.*signa[COVER]);
- fprintf(stdout, " Temperature of %s clouds\n", (warm_ambiguous ? "cold" : "all"));
- fprintf(stdout, " Maximum: %.2lf K\n", signa[KMAX]);
- fprintf(stdout, " Mean : %.2lf K\n", signa[KMEAN]);
- fprintf(stdout, " Minimum: %.2lf K\n", signa[KMIN]);
+ G_message(" PRELIMINARY SCENE ANALYSIS\n");
+ G_message(" Desert index: %.3lf\n", idesert);
+ G_message(" Snow cover : %.3lf %%\n", 100. * value[SNOW]);
+ G_message(" Cloud cover : %.3lf %%\n", 100. * signa[COVER]);
+ G_message(" Temperature of %s clouds\n",
+ (warm_ambiguous ? "cold" : "all"));
+ G_message(" Maximum: %.2lf K\n", signa[KMAX]);
+ G_message(" Mean : %.2lf K\n", signa[KMEAN]);
+ G_message(" Minimum: %.2lf K\n", signa[KMIN]);
- /* WARNING: variable 'value' reutilizada con nuevos valores */
+ /* WARNING: variable 'value' reutilizada con nuevos valores */
- if (idesert > 0.5 && signa[COVER] > 0.004 && signa[KMEAN] < 295.)
- {
- value[KUPPER] = quantile( 0.975, hist_cold ) + K_BASE;
- value[KLOWER] = quantile( 0.835, hist_cold ) + K_BASE;
- value[MEAN] = mean(hist_cold) + K_BASE; /* quantile( 0.5, hist_cold ): */
- value[DSTD] = sqrt( moment(2, hist_cold, 1) );
- value[SKEW] = moment( 3, hist_cold, 3 ) / pow( value[DSTD], 3);
+ if (idesert > 0.5 && signa[COVER] > 0.004 && signa[KMEAN] < 295.) {
+ value[KUPPER] = quantile(0.975, hist_cold) + K_BASE;
+ value[KLOWER] = quantile(0.835, hist_cold) + K_BASE;
+ value[MEAN] = mean(hist_cold) + K_BASE; /* quantile( 0.5, hist_cold ): */
+ value[DSTD] = sqrt(moment(2, hist_cold, 1));
+ value[SKEW] = moment(3, hist_cold, 3) / pow(value[DSTD], 3);
- if (value[SKEW] > 0.)
- {
- shift = value[SKEW];
- if (shift > 1.) shift = 1.;
- shift *= value[DSTD];
+ if (value[SKEW] > 0.) {
+ shift = value[SKEW];
+ if (shift > 1.)
+ shift = 1.;
+ shift *= value[DSTD];
- x = quantile( 0.9875, hist_cold ) + K_BASE;
- if ((value[KUPPER] + shift) > x)
- {
- value[KUPPER] = x;
- value[KLOWER] = x - shift; /** ??? COMPROBAR ***/
- }
- else
- {
- value[KUPPER] += shift;
- value[KLOWER] += shift;
- }
- }
+ x = quantile(0.9875, hist_cold) + K_BASE;
+ if ((value[KUPPER] + shift) > x) {
+ value[KUPPER] = x;
+ value[KLOWER] = x - shift; /* ??? COMPROBAR */
+ }
+ else {
+ value[KUPPER] += shift;
+ value[KLOWER] += shift;
+ }
+ }
- fprintf(stdout, " HISTOGRAM CLOUD SIGNATURE\n");
- fprintf(stdout, " Histogram classes: %d\n", hist_n);
- 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, " 97.50 quantile: %.2lf K\n", value[KUPPER]);
- fprintf(stdout, " 83.50 quantile: %.2lf K\n", value[KLOWER]);
+ G_message(" HISTOGRAM CLOUD SIGNATURE\n");
+ G_message(" Histogram classes: %d\n", hist_n);
+ G_message(" Mean temperature: %.2lf K\n", value[MEAN]);
+ G_message(" Standard deviation: %.2lf \n", value[DSTD]);
+ G_message(" Skewness: %.2lf \n", value[SKEW]);
+ G_message(" 97.50 quantile: %.2lf K\n", value[KUPPER]);
+ G_message(" 83.50 quantile: %.2lf K\n", value[KLOWER]);
+ }
+ else {
+ if (signa[KMEAN] < 295.) {
+ /* Retained warm and cold clouds */
+ G_message(" Scene with clouds\n");
+ warm_ambiguous = 0;
+ value[KUPPER] = 0.;
+ value[KLOWER] = 0.;
}
- else
- {
- if (signa[KMEAN] < 295.)
- {
- /* Retained warm and cold clouds */
- fprintf(stdout, " Scene with clouds\n");
- warm_ambiguous = 0;
- value[KUPPER] = 0.;
- value[KLOWER] = 0.;
- }
- else
- {
- /* Retained cold clouds */
- fprintf(stdout, " Scene cloud free\n");
- warm_ambiguous = 1;
- value[KUPPER] = 0.;
- value[KLOWER] = 0.;
- }
+ else {
+ /* Retained cold clouds */
+ G_message(" Scene cloud free\n");
+ warm_ambiguous = 1;
+ value[KUPPER] = 0.;
+ value[KLOWER] = 0.;
}
+ }
- /* SECOND FILTER ... */
- /* By-pass two processing but it retains warm and cold clouds */
- if (two_pass == 0)
- {
- warm_ambiguous = 0;
- value[KUPPER] = 0.;
- value[KLOWER] = 0.;
- }
- acca_second(verbose, out, band[BAND6],
- warm_ambiguous, value[KUPPER], value[KLOWER]);
- /* CATEGORIES: WARM_CLOUD, COLD_CLOUD, NULL (= NO_CLOUD) */
+ /* SECOND FILTER ... */
+ /* By-pass two processing but it retains warm and cold clouds */
+ if (two_pass == 0) {
+ warm_ambiguous = 0;
+ value[KUPPER] = 0.;
+ value[KLOWER] = 0.;
+ }
+ acca_second(verbose, out, band[BAND6],
+ warm_ambiguous, value[KUPPER], value[KLOWER]);
+ /* CATEGORIES: WARM_CLOUD, COLD_CLOUD, NULL (= NO_CLOUD) */
- return;
+ return;
}
void acca_first(int verbose, Gfile * out, Gfile band[],
- int with_shadow,
- int count[], int cold[], int warm[], double stats[])
+ int with_shadow,
+ int count[], int cold[], int warm[], double stats[])
{
- int i, row, col, nrows, ncols;
+ int i, row, col, nrows, ncols;
- char code;
- double pixel[5], nsdi, rat56, rat45;
+ char code;
+ double pixel[5], nsdi, rat56, rat45;
- /* Creation of output file */
- out->rast = G_allocate_raster_buf(CELL_TYPE);
- if ((out->fd = G_open_raster_new(out->name, CELL_TYPE)) < 0)
- G_fatal_error(_("Could not open <%s>"), out->name);
+ /* Creation of output file */
+ out->rast = G_allocate_raster_buf(CELL_TYPE);
+ if ((out->fd = G_open_raster_new(out->name, CELL_TYPE)) < 0)
+ G_fatal_error(_("Could not open <%s>"), out->name);
- /* ----- ----- */
- fprintf(stdout, "Pass one processing ... \n");
+ /* ----- ----- */
+ G_message("Pass one processing ... \n");
- stats[SUM_COLD] = 0.;
- stats[SUM_WARM] = 0.;
- stats[KMAX] = 0.;
- stats[KMIN] = 10000.;
+ stats[SUM_COLD] = 0.;
+ stats[SUM_WARM] = 0.;
+ stats[KMAX] = 0.;
+ stats[KMIN] = 10000.;
- nrows = G_window_rows();
- ncols = G_window_cols();
+ nrows = G_window_rows();
+ ncols = G_window_cols();
- for (row = 0; row < nrows; row++)
- {
- if (verbose)
- {
- G_percent(row, nrows, 2);
+ for (row = 0; row < nrows; row++) {
+ if (verbose) {
+ G_percent(row, nrows, 2);
+ }
+ for (i = BAND2; i <= BAND6; i++) {
+ if (G_get_d_raster_row(band[i].fd, band[i].rast, row) < 0)
+ G_fatal_error(_("Could not read row from <%s>"),
+ band[i].name);
+ }
+ for (col = 0; col < ncols; col++) {
+ code = NO_DEFINED;
+ /* Null when null pixel in any band */
+ for (i = BAND2; i <= BAND6; i++) {
+ if (G_is_d_null_value((void *)((DCELL *) band[i].rast + col))) {
+ code = NO_CLOUD;
+ break;
}
- for (i = BAND2; i <= BAND6; i++)
- {
- if (G_get_d_raster_row(band[i].fd, band[i].rast, row) < 0)
- G_fatal_error(_("Could not read row from <%s>"), band[i].name);
- }
- for (col = 0; col < ncols; col++)
- {
- code = NO_DEFINED;
- /* Null when null pixel in any band */
- for ( i = BAND2; i <= BAND6; i++ )
- {
- if (G_is_d_null_value((void *)((DCELL *) band[i].rast + col)))
- {
- code = NO_CLOUD;
- break;
- }
- pixel[i] = (double)((DCELL *) band[i].rast)[col];
- }
- /* Determina los pixeles de sombras */
- if (code == NO_DEFINED && with_shadow)
- {
- code = shadow_algorithm(pixel);
- }
- /* Analiza el valor de los pixeles no definidos */
- if (code == NO_DEFINED)
- {
- code = NO_CLOUD;
- count[TOTAL]++;
- nsdi = (pixel[BAND2] - pixel[BAND5]) /
- (pixel[BAND2] + pixel[BAND5]);
- /* ----------------------------------------------------- */
- /* Brightness Thershold: Eliminates dark images */
- if (pixel[BAND3] > th_1)
- {
- /* Normalized Snow Difference Index: Eliminates many types of snow */
- if (nsdi > th_2[0] && nsdi < th_2[1])
- {
- /* Temperature Thershold: Eliminates warm image features */
- if (pixel[BAND6] < th_3)
- {
- rat56 = (1 - pixel[BAND5]) * pixel[BAND6];
- /* 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);
- }
- 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 = NO_DEFINED; count[SOIL]++; }
- }
- else code = NO_DEFINED;
- }
- else code = (pixel[BAND5] < th_4_b) ? NO_CLOUD : NO_DEFINED;
- }
- else code = NO_CLOUD;
+ pixel[i] = (double)((DCELL *) band[i].rast)[col];
+ }
+ /* Determina los pixeles de sombras */
+ if (code == NO_DEFINED && with_shadow) {
+ code = shadow_algorithm(pixel);
+ }
+ /* Analiza el valor de los pixeles no definidos */
+ if (code == NO_DEFINED) {
+ code = NO_CLOUD;
+ count[TOTAL]++;
+ nsdi = (pixel[BAND2] - pixel[BAND5]) /
+ (pixel[BAND2] + pixel[BAND5]);
+ /* ----------------------------------------------------- */
+ /* Brightness Thershold: Eliminates dark images */
+ if (pixel[BAND3] > th_1) {
+ /* Normalized Snow Difference Index: Eliminates many types of snow */
+ if (nsdi > th_2[0] && nsdi < th_2[1]) {
+ /* Temperature Thershold: Eliminates warm image features */
+ if (pixel[BAND6] < th_3) {
+ rat56 = (1. - pixel[BAND5]) * pixel[BAND6];
+ /* 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);
}
- else
- {
- code = NO_CLOUD;
- if (nsdi > th_2_b) count[SNOW]++;
- }
+ 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 = NO_DEFINED;
+ count[SOIL]++;
+ }
}
- else
- {
- code = (pixel[BAND3] < th_1_b) ? NO_CLOUD : NO_DEFINED;
- }
- /* ----------------------------------------------------- */
+ else {
+ code = NO_DEFINED;
+ }
+ }
+ else {
+ code =
+ (pixel[BAND5] <
+ th_4_b) ? NO_CLOUD : NO_DEFINED;
+ }
}
- if (code == NO_CLOUD)
- {
- G_set_c_null_value((CELL *) out->rast + col, 1);
+ else {
+ code = NO_CLOUD;
}
- else
- {
- ((CELL *) out->rast)[col] = code;
- }
+ }
+ else {
+ code = NO_CLOUD;
+ if (nsdi > th_2_b)
+ count[SNOW]++;
+ }
}
- if (G_put_raster_row(out->fd, out->rast, CELL_TYPE) < 0)
- {
- G_fatal_error(_("Cannot write row to <%s>"), out->name);
+ else {
+ code = (pixel[BAND3] < th_1_b) ? NO_CLOUD : NO_DEFINED;
}
+ /* ----------------------------------------------------- */
+ }
+ if (code == NO_CLOUD) {
+ G_set_c_null_value((CELL *) out->rast + col, 1);
+ }
+ else {
+ ((CELL *) out->rast)[col] = code;
+ }
}
- /* ----- ----- */
+ if (G_put_raster_row(out->fd, out->rast, CELL_TYPE) < 0) {
+ G_fatal_error(_("Cannot write row to <%s>"), out->name);
+ }
+ }
- G_free(out->rast);
- G_close_cell(out->fd);
+ G_free(out->rast);
+ G_close_cell(out->fd);
- return;
+ return;
}
void acca_second(int verbose, Gfile * out, Gfile band,
- int warm_ambiguous, double upper, double lower)
+ int warm_ambiguous, double upper, double lower)
{
- int i, row, col, nrows, ncols;
- char *mapset;
+ int i, row, col, nrows, ncols;
+ char *mapset;
- int code;
- double temp;
- Gfile tmp;
+ int code;
+ double temp;
+ Gfile tmp;
- /* Open to read */
- mapset = G_find_cell2(out->name, "");
- if (mapset == NULL)
- G_fatal_error("cell file [%s] not found", out->name);
- out->rast = G_allocate_raster_buf(CELL_TYPE);
- if ((out->fd = G_open_cell_old(out->name, mapset)) < 0)
- G_fatal_error("Cannot open cell file [%s]", out->name);
+ /* Open to read */
+ mapset = G_find_cell2(out->name, "");
+ if (mapset == NULL)
+ G_fatal_error("cell file [%s] not found", out->name);
+ out->rast = G_allocate_raster_buf(CELL_TYPE);
+ if ((out->fd = G_open_cell_old(out->name, mapset)) < 0)
+ G_fatal_error("Cannot open cell file [%s]", out->name);
- /* Open to write */
- sprintf(tmp.name, "_%d.BBB", getpid()) ;
- tmp.rast = G_allocate_raster_buf(CELL_TYPE);
- if ((tmp.fd = G_open_raster_new(tmp.name, CELL_TYPE)) < 0)
- G_fatal_error(_("Could not open <%s>"), tmp.name);
+ /* Open to write */
+ sprintf(tmp.name, "_%d.BBB", getpid());
+ tmp.rast = G_allocate_raster_buf(CELL_TYPE);
+ if ((tmp.fd = G_open_raster_new(tmp.name, CELL_TYPE)) < 0)
+ G_fatal_error(_("Could not open <%s>"), tmp.name);
- if (upper == 0.)
- fprintf(stdout, "Removing ambiguous pixels ... \n");
- else
- fprintf(stdout, "Pass two processing ... \n");
+ if (upper == 0.)
+ G_message("Removing ambiguous pixels ... \n");
+ else
+ G_message("Pass two processing ... \n");
- nrows = G_window_rows();
- ncols = G_window_cols();
+ nrows = G_window_rows();
+ ncols = G_window_cols();
- for (row = 0; row < nrows; row++)
- {
- if (verbose)
- {
- G_percent(row, nrows, 2);
- }
- if (G_get_d_raster_row(band.fd, band.rast, row) < 0)
- G_fatal_error(_("Could not read from <%s>"), band.name);
- if (G_get_c_raster_row(out->fd, out->rast, row) < 0)
- G_fatal_error(_("Could not read from <%s>"), out->name);
+ for (row = 0; row < nrows; row++) {
+ if (verbose) {
+ G_percent(row, nrows, 2);
+ }
+ if (G_get_d_raster_row(band.fd, band.rast, row) < 0)
+ G_fatal_error(_("Could not read from <%s>"), band.name);
+ if (G_get_c_raster_row(out->fd, out->rast, row) < 0)
+ G_fatal_error(_("Could not read from <%s>"), out->name);
- for (col = 0; col < ncols; col++)
- {
- if (G_is_c_null_value((void *)((CELL *) out->rast + col)))
- {
- G_set_c_null_value((CELL *) tmp.rast + col, 1);
- }
- else
- {
- code = (int)((CELL *) out->rast)[col];
- /* Resolve ambiguous pixels */
- if (code == NO_DEFINED ||
- (code == WARM_CLOUD) && warm_ambiguous == 1)
- {
- temp = (double)((DCELL *) band.rast)[col];
- if (temp > upper)
- {
- G_set_c_null_value((CELL *) tmp.rast + col, 1);
- }
- else
- {
- ((CELL *) tmp.rast)[col] = (temp < lower) ? IS_COLD_CLOUD : IS_WARM_CLOUD;
- }
- }
- else
- /* Join warm (not ambiguous) and cold clouds */
- if (code == COLD_CLOUD ||
- (code == WARM_CLOUD) && warm_ambiguous == 0)
- {
- ((CELL *) tmp.rast)[col] = IS_COLD_CLOUD;
- }
- else
- ((CELL *) tmp.rast)[col] = IS_SHADOW;
- }
+ for (col = 0; col < ncols; col++) {
+ if (G_is_c_null_value((void *)((CELL *) out->rast + col))) {
+ G_set_c_null_value((CELL *) tmp.rast + col, 1);
+ }
+ else {
+ code = (int)((CELL *) out->rast)[col];
+ /* Resolve ambiguous pixels */
+ if (code == NO_DEFINED ||
+ (code == WARM_CLOUD) && warm_ambiguous == 1) {
+ temp = (double)((DCELL *) band.rast)[col];
+ if (temp > upper) {
+ G_set_c_null_value((CELL *) tmp.rast + col, 1);
+ }
+ else {
+ ((CELL *) tmp.rast)[col] =
+ (temp < lower) ? IS_COLD_CLOUD : IS_WARM_CLOUD;
+ }
}
- if (G_put_raster_row(tmp.fd, tmp.rast, CELL_TYPE) < 0)
- {
- G_fatal_error(_("Cannot write to <%s>"), tmp.name);
+ else
+ /* Join warm (not ambiguous) and cold clouds */
+ if (code == COLD_CLOUD ||
+ (code == WARM_CLOUD) && warm_ambiguous == 0) {
+ ((CELL *) tmp.rast)[col] = IS_COLD_CLOUD;
}
+ else
+ ((CELL *) tmp.rast)[col] = IS_SHADOW;
+ }
}
+ if (G_put_raster_row(tmp.fd, tmp.rast, CELL_TYPE) < 0) {
+ G_fatal_error(_("Cannot write to <%s>"), tmp.name);
+ }
+ }
- /* Finalización */
+ /* Finalización */
- G_free(tmp.rast);
- G_close_cell(tmp.fd);
+ G_free(tmp.rast);
+ G_close_cell(tmp.fd);
- G_free(out->rast);
- G_close_cell(out->fd);
+ G_free(out->rast);
+ G_close_cell(out->fd);
- G_remove("cats", out->name);
- G_remove("cell", out->name);
- G_remove("cellhd", out->name);
- G_remove("cell_misc", out->name);
- G_remove("hist", out->name);
+ G_remove("cats", out->name);
+ G_remove("cell", out->name);
+ G_remove("cellhd", out->name);
+ G_remove("cell_misc", out->name);
+ G_remove("hist", out->name);
- G_rename("cats", tmp.name, out->name);
- G_rename("cell", tmp.name, out->name);
- G_rename("cellhd", tmp.name, out->name);
- G_rename("cell_misc", tmp.name, out->name);
- G_rename("hist", tmp.name, out->name);
+ G_rename("cats", tmp.name, out->name);
+ G_rename("cell", tmp.name, out->name);
+ G_rename("cellhd", tmp.name, out->name);
+ G_rename("cell_misc", tmp.name, out->name);
+ G_rename("hist", tmp.name, out->name);
- return;
+ return;
}
/**********************************************************
@@ -467,16 +449,11 @@
int shadow_algorithm(double pixel[])
{
-// return NO_DEFINED;
+ if (pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. && pixel[BAND4] / pixel[BAND2] > 1. // Quita agua 1
+ // && (pixel[BAND3] - pixel[BAND5]) / (pixel[BAND3] + pixel[BAND5]) < 0.10
+ ) {
+ return IS_SHADOW;
+ }
- if( pixel[BAND3] < 0.07
- && (1 - pixel[BAND4]) * pixel[BAND6] > 240.
- && pixel[BAND4] / pixel[BAND2] > 1. // Quita agua 1
-// && (pixel[BAND3] - pixel[BAND5]) / (pixel[BAND3] + pixel[BAND5]) < 0.10
- )
- {
- return IS_SHADOW;
- }
-
- return NO_DEFINED;
+ return NO_DEFINED;
}
Modified: grass-addons/imagery/i.landsat.acca/description.html
===================================================================
--- grass-addons/imagery/i.landsat.acca/description.html 2010-09-28 01:15:57 UTC (rev 43722)
+++ grass-addons/imagery/i.landsat.acca/description.html 2010-09-28 12:12:02 UTC (rev 43723)
@@ -1,51 +1,30 @@
<H2>DESCRIPTION</H2>
-<EM>i.landsat.acca</EM> implements the Automated Cloud-Cover Assessment
-(ACCA) Algorithm from Irish (2000) with the constant values for pass filter
-one from Irish et al. (2006). To do this, it needs the band numbers 2, 3,
-4, 5, and 6 (or band 61 for Landsat-7 ETM+) which have already been
-processed from DN into reflectance and band-6 temperature (for example with
-<em>i.landsat.toar</em>).
+<EM>i.landsat.acca</EM> is to implement the Automated Cloud-Cover Assessment (ACCA) Algorithm from Irish (2000) with the constant values for pass filter one from Irish et al. (2006). To do, it need the band numbers 2, 3, 4, 5, and 6 (band 61 of Landsat-7 ETM+).</p>
-<p>
-The ACCA gives good results over most of the planet with the exception of
-ice sheets because ACCA operates on the premise that clouds are colder than
-the land surface they cover. The algorithm was designed to Landsat-7 ETM+
-but because of they use reflectances it usefull also for Landsat-4/5 TM.
+<p>The ACCA gives good results over most of the planet with the exception of ice sheets because ACCA operates on the premise that clouds are colder than the land surface they cover. The algorithm was designed to Landsat-7 ETM+ but because of they use reflectances it usefull also for Landsat-4/5 TM.</p>
+
<H2>NOTES</H2>
-<em>i.landsat.acca</em> works on the full raster, thus it does not take
-account the current region settings. <!-- this may have to change before the
-module can be added to the main source archive -->
+<em>i.landsat.acca</em> use full raster and it do not take account the current region.
-
<H2>EXAMPLES</H2>
-To run the standard ACCA algorithm with two pass processing (the <b>-2</b> flag)
-and filling cloud holes (the <b>-f</b> flag):
-<p>
-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>:
+<p>To make standard ACCA algorithm with two pass processing (flag 2) and filling cloud holes (flag f).</p>
+Reflectances in band rasters 226_62.toar.1, 226_62.toar.2 [...] and thermal band 226_62.toar.61 to output file 226_62.acca:</p>
<div class="code"><pre>
-i.landsat.acca -2 -f band_prefix=226_62.toar output=226_62.acca
+i.landsat.acca -f2 band=226_62.toar out=226_62.acca
</pre></div>
+
<H2>REFERENCES</H2>
<ol>
- <li>Irish R.R., Barker J.L., Goward S.N., and Arvidson T., 2006.
- Characterization of the Landsat-7 ETM+ Automated Cloud-Cover
- Assessment (ACCA) Algorithm. Photogrammetric Engineering and Remote
- Sensing vol. 72(10): 1179-1188.</li>
-
- <li>Irish, R.R., 2000. Landsat 7 Automatic Cloud Cover Assessment. In
- S.S. Shen and M.R. Descour (Eds.): Algorithms for Multispectral,
- Hyperspectral, and Ultraspectral Imagery VI. Proceedings of SPIE,
- 4049: 348-355.</li>
+ <li>Irish R.R., Barker J.L., Goward S.N., and Arvidson T., 2006. Characterization of the Landsat-7 ETM+ Automated Cloud-Cover Assessment (ACCA) Algorithm. Photogrammetric Engineering and Remote Sensing vol. 72(10): 1179-1188.</li>
+ <li>Irish, R.R., 2000. Landsat 7 Automatic Cloud Cover Assessment. In S.S. Shen and M.R. Descour (Eds.): Algorithms for Multispectral, Hyperspectral, and Ultraspectral Imagery VI. Proceedings of SPIE, 4049: 348-355.</li>
</ol>
@@ -60,8 +39,7 @@
<H2>AUTHOR</H2>
E. Jorge Tizado (ej.tizado unileon es)<br>
-Dept. Biodiversity and Environmental Management,
-University of León, Spain<BR>
+Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
<p>
<i>Last changed: $Date: 2007/07/07 00:00:00 $</i>
Modified: grass-addons/imagery/i.landsat.acca/local_proto.h
===================================================================
--- grass-addons/imagery/i.landsat.acca/local_proto.h 2010-09-28 01:15:57 UTC (rev 43722)
+++ grass-addons/imagery/i.landsat.acca/local_proto.h 2010-09-28 12:12:02 UTC (rev 43723)
@@ -21,18 +21,18 @@
typedef struct
{
- int fd;
- void * rast;
- char name[1024];
+ int fd;
+ void *rast;
+ char name[1024];
} Gfile;
-void acca_algorithm(int, Gfile *, Gfile [], int, int);
-void acca_first (int, Gfile *, Gfile [], int, int [], int [], int [], double []);
+void acca_algorithm(int, Gfile *, Gfile[], int, int);
+void acca_first(int, Gfile *, Gfile[], int, int[], int[], int[], double[]);
void acca_second(int, Gfile *, Gfile, int, double, double);
-int shadow_algorithm(double []);
+int shadow_algorithm(double[]);
void filter_holes(int, Gfile *);
Modified: grass-addons/imagery/i.landsat.acca/main.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/main.c 2010-09-28 01:15:57 UTC (rev 43722)
+++ grass-addons/imagery/i.landsat.acca/main.c 2010-09-28 12:12:02 UTC (rev 43723)
@@ -55,27 +55,27 @@
mapset = G_find_cell2(raster_name, "");
if (mapset == NULL) {
- G_warning(_("Raster map <%s> not found"), raster_name);
+ G_message("cell file [%s] not found", raster_name);
return -1;
}
if (G_legal_filename(raster_name) < 0) {
- G_warning(_("<%s> is an illegal file name"), raster_name);
+ G_message("[%s] is an illegal name", raster_name);
return -1;
}
if ((raster_fd = G_open_cell_old(raster_name, mapset)) < 0) {
- G_warning(_("Unable to open raster map <%s>"), raster_name);
+ G_message("Cannot open cell file [%s]", raster_name);
return -1;
}
if (G_get_cellhd(raster_name, mapset, &cellhd) < 0) {
- G_warning(_("Unable to read header of raster map <%s>"), raster_name);
+ G_message("Cannot read file header of [%s]", raster_name);
return -1;
}
if (G_set_window(&cellhd) < 0) {
- G_warning(_("Cannot reset current region"));
+ G_message("Unable to set region");
return -1;
}
if ((map_type = G_raster_map_type(raster_name, mapset)) != DCELL_TYPE) {
- G_warning(_("Map is not DCELL type (process DN to radiance first)"));
+ G_message("Map is not of DCELL_TYPE");
return -1;
}
return raster_fd;
@@ -95,8 +95,7 @@
struct Option *input, *output, *hist;
struct Flag *shadow, *sat5, *filter, *pass2;
char *in_name, *out_name;
- struct Categories cats;
- char title[RECORD_LEN];
+ double p;
Gfile band[5], out;
@@ -109,21 +108,25 @@
_("Landsat TM/ETM+ Automatic Cloud Cover Assessment (ACCA)");
input = G_define_option();
- input->key = "band_prefix";
+ input->key = _("band_prefix");
input->type = TYPE_STRING;
input->required = YES;
- input->gisprompt = "input,cell,raster";
+ input->gisprompt = _("input,cell,raster");
input->description =
_("Base name of the landsat band rasters ([band_prefix].[band_number])");
- output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output = G_define_option();
+ output->key = _("output");
+ output->type = TYPE_STRING;
output->required = YES;
+ output->gisprompt = _("output,cell,raster");
+ output->description = _("Output file name");
-
hist = G_define_option();
- hist->key = "histogram";
+ hist->key = _("histogram");
hist->type = TYPE_INTEGER;
hist->required = NO;
+ hist->gisprompt = _("input,integer");
hist->description =
_("Number of classes in the cloud temperature histogram");
hist->answer = "100";
@@ -131,23 +134,26 @@
sat5 = G_define_flag();
sat5->key = '5';
sat5->description = _("Landsat-5 TM");
+ sat5->answer = 0;
filter = G_define_flag();
filter->key = 'f';
- filter->description = _("Apply post-processing filter to remove small holes");
+ filter->description = _("Use final filter holes");
+ filter->answer = 0;
pass2 = G_define_flag();
pass2->key = '2';
- pass2->description = _("Use two-pass processing");
+ pass2->description = _("With pass two processing");
+ pass2->answer = 0;
shadow = G_define_flag();
shadow->key = 's';
- shadow->description = _("Include a category for cloud shadows");
+ shadow->description = _("Add class for cloud shadows");
+ shadow->answer = 0;
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
-
/* stores OPTIONS and FLAGS to variables */
hist_n = atoi(hist->answer);
@@ -160,7 +166,7 @@
snprintf(band[i].name, 127, "%s.%d%c", in_name, i + 2,
(i == BAND6 && !sat5->answer ? '1' : '\0'));
if ((band[i].fd = check_raster(band[i].name)) < 0) {
- G_fatal_error(_("Error in map name <%s>!"), band[i].name);
+ G_fatal_error(_("Error in filename [%s]!"), band[i].name);
}
band[i].rast = G_allocate_raster_buf(DCELL_TYPE);
}
@@ -169,7 +175,7 @@
snprintf(out.name, 127, "%s", out_name);
if (G_legal_filename(out_name) < 0)
- G_fatal_error(_("<%s> is an illegal file name"), out.name);
+ G_fatal_error(_("[%s] is an illegal name"), out.name);
/* --------------------------------------- */
@@ -189,22 +195,10 @@
G_close_cell(band[i].fd);
}
- /* write out map title and category labels */
- G_init_cats((CELL) 0, "", &cats);
- sprintf(title, "LANDSAT-%s Automatic Cloud Cover Assessment",
- sat5->answer ? "5 TM" : "7 ETM+");
- G_set_raster_cats_title(title, &cats);
+ // struct Categories cats;
+ // G_read_raster_cats(out.name, char *mapset, cats)
+ // G_write_raster_cats(out.name, &cats);
- G_set_cat(IS_SHADOW, "Shadow", &cats);
- G_set_cat(IS_COLD_CLOUD, "Cold cloud", &cats);
- G_set_cat(IS_WARM_CLOUD, "Warm cloud", &cats);
-
- if (G_write_cats(out.name, &cats) <= 0)
- G_warning(_("Cannot write category file for raster map <%s>"),
- out.name);
- G_free_cats(&cats);
-
- /* write out command line opts */
G_short_history(out.name, "raster", &history);
G_command_history(&history);
G_write_history(out.name, &history);
Modified: grass-addons/imagery/i.landsat.acca/tools.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/tools.c 2010-09-28 01:15:57 UTC (rev 43722)
+++ grass-addons/imagery/i.landsat.acca/tools.c 2010-09-28 12:12:02 UTC (rev 43723)
@@ -21,19 +21,21 @@
/* Global variable
allow use as parameter in the command line */
-int hist_n = 100; /* interval of real data 100/hist_n */
+int hist_n = 100; /* interval of real data 100/hist_n */
void hist_put(double t, int hist[])
{
- int i;
+ int i;
- /* scale factor */
- i = (int)(t * ((double)hist_n/100.));
+ /* scale factor */
+ i = (int)(t * ((double)hist_n / 100.));
- if (i < 1) i = 1;
- if (i > hist_n) i = hist_n;
+ if (i < 1)
+ i = 1;
+ if (i > hist_n)
+ i = hist_n;
- hist[i-1] += 1;
+ hist[i - 1] += 1;
}
/* histogram moment */
@@ -45,32 +47,30 @@
k = 0;
total = 0;
- hmean = 0.;
- for( i = 0; i < hist_n; i++ )
- {
- total += hist[i];
- hmean += (double)(i * hist[i]);
+ hmean = 0.;
+ for (i = 0; i < hist_n; i++) {
+ total += hist[i];
+ hmean += (double)(i * hist[i]);
}
- hmean /= (double)total; /* histogram mean */
+ hmean /= (double)total; /* histogram mean */
/*
- value = 0.;
- for( i = 0; i < hist_n; i++ )
- {
- cte = 1.;
- for( j = 0; j < n; j++ ) cte *= (i - hmean);
- value += cte * (double)hist[i]/(double)total;
- }
+ value = 0.;
+ for( i = 0; i < hist_n; i++ )
+ {
+ cte = 1.;
+ for( j = 0; j < n; j++ ) cte *= (i - hmean);
+ value += cte * (double)hist[i]/(double)total;
+ }
- /* remove scale factor *
- for( j = 0; j < n; j++ ) value /= ((double)hist_n/100.);
- */
+ /* remove scale factor *
+ for( j = 0; j < n; j++ ) value /= ((double)hist_n/100.);
+ */
value = 0.;
cte = 100. / ((double)hist_n * (double)(total - k));
- for( i = 0; i < hist_n; i++ )
- {
- value += (pow((i - hmean), n) * (double)hist[i] * cte);
+ for (i = 0; i < hist_n; i++) {
+ value += (pow((i - hmean), n) * (double)hist[i] * cte);
}
return value;
@@ -83,16 +83,15 @@
double value, mean;
total = 0;
- mean = 0.;
- for( i = 0; i < hist_n; i++ )
- {
- total += hist[i];
- mean += (double)(i * hist[i]);
+ mean = 0.;
+ for (i = 0; i < hist_n; i++) {
+ total += hist[i];
+ mean += (double)(i * hist[i]);
}
mean /= (double)total;
/* remove scale factor */
- return (mean / ((double)hist_n/100.));
+ return (mean / ((double)hist_n / 100.));
}
/* Real data quantile */
@@ -102,23 +101,20 @@
double value, qmax, qmin;
total = 0;
- for( i = 0; i < hist_n; i++ )
- {
- total += hist[i];
+ for (i = 0; i < hist_n; i++) {
+ total += hist[i];
}
qmax = 1.;
- for( i = hist_n - 1; i >= 0; i-- )
- {
- qmin = qmax - (double)hist[i]/(double)total;
- if( q >= qmin )
- {
- value = (q - qmin) / (qmax - qmin) + (i - 1);
- break;
- }
- qmax = qmin;
+ for (i = hist_n - 1; i >= 0; i--) {
+ qmin = qmax - (double)hist[i] / (double)total;
+ if (q >= qmin) {
+ value = (q - qmin) / (qmax - qmin) + (i - 1);
+ break;
+ }
+ qmax = qmin;
}
/* remove scale factor */
- return (value / ((double)hist_n/100.));
+ return (value / ((double)hist_n / 100.));
}
More information about the grass-commit
mailing list