[GRASS-SVN] r29969 - in grass-addons: . i.landsat.acca
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Feb 6 02:59:36 EST 2008
Author: neteler
Date: 2008-02-06 02:59:36 -0500 (Wed, 06 Feb 2008)
New Revision: 29969
Added:
grass-addons/i.landsat.acca/
grass-addons/i.landsat.acca/Makefile
grass-addons/i.landsat.acca/algorithm.c
grass-addons/i.landsat.acca/description.html
grass-addons/i.landsat.acca/local_proto.h
grass-addons/i.landsat.acca/main.c
grass-addons/i.landsat.acca/tools.c
Log:
E. Jorge Tizado: Automatic Cloud Cover Assessment (ACCA) (Irish 2000)
Added: grass-addons/i.landsat.acca/Makefile
===================================================================
--- grass-addons/i.landsat.acca/Makefile (rev 0)
+++ grass-addons/i.landsat.acca/Makefile 2008-02-06 07:59:36 UTC (rev 29969)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.landsat.acca
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/i.landsat.acca/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/i.landsat.acca/algorithm.c
===================================================================
--- grass-addons/i.landsat.acca/algorithm.c (rev 0)
+++ grass-addons/i.landsat.acca/algorithm.c 2008-02-06 07:59:36 UTC (rev 29969)
@@ -0,0 +1,479 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+#define SCALE 100.
+
+/* value and count */
+#define TOTAL 0
+#define WARM 1
+#define COLD 2
+#define SNOW 3
+#define SOIL 4
+
+/* signa */
+#define COVER 1
+#define SUM_COLD 0
+#define SUM_WARM 1
+#define KMEAN 2
+#define KMAX 3
+#define KMIN 4
+
+/* re-use value */
+#define KLOWER 0
+#define KUPPER 1
+#define MEAN 2
+#define SKEW 3
+#define VARI 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
+ *
+ **********************************************************/
+
+/*--------------------------------------------------------
+ CONSTANTS
+ Usar esta forma para que via extern puedan modificarse
+ como opciones desde el programa main.
+ ---------------------------------------------------------*/
+
+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_b = 0.8;
+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 */
+
+extern int hist_n;
+
+#define K_BASE 240.
+
+void acca_algorithm(int verbose, Gfile * out, Gfile band[],
+ 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;
+
+ /* Initialized varibles ... */
+ 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) */
+
+ 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. : value[0] /
+ (value[0] + (double)count[SOIL]));
+
+ /* 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];
+ }
+
+ 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]);
+
+ /* 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;
+ value[VARI] = moment( 2, hist_cold );
+ value[SKEW] = moment( 3, hist_cold );
+
+ if (value[SKEW] > 0.)
+ {
+ shift = value[SKEW];
+ if (shift > 1.) shift = 1.;
+ shift *= sqrt( value[VARI] );
+
+ 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", sqrt(value[VARI]));
+ fprintf(stdout, " Skewness: %.2lf\n", value[SKEW]);
+ fprintf(stdout, " 97.50 percentile: %.2lf K\n", value[KUPPER]);
+ fprintf(stdout, " 83.50 percentile: %.2lf K\n", value[KLOWER]);
+// fprintf(stdout, " 50.00 percentile: %.2lf K\n", quantile(0.5,hist_cold) + K_BASE);
+ }
+ 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.;
+ }
+ }
+
+ /* SECOND FILTER ... */
+ /* By-pass two processing by 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;
+}
+
+
+void acca_first(int verbose, Gfile * out, Gfile band[],
+ int with_shadow,
+ int count[], int cold[], int warm[], double stats[])
+{
+ int i, row, col, nrows, ncols;
+
+ 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);
+
+ /* ----- ----- */
+ fprintf(stdout, "Pass one processing ... \n");
+
+ stats[SUM_COLD] = 0.;
+ stats[SUM_WARM] = 0.;
+ stats[KMAX] = 0.;
+ stats[KMIN] = 10000.;
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ 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;
+ }
+ 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[BAND4]) * 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;
+ }
+ else { code = NO_CLOUD; if (nsdi > th_2_b) count[SNOW]++; }
+ }
+ 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);
+
+ return;
+}
+
+
+void acca_second(int verbose, Gfile * out, Gfile band,
+ int warm_ambiguous, double upper, double lower)
+{
+ int i, row, col, nrows, ncols;
+ char *mapset;
+
+ 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 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");
+
+ 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 (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
+ {
+ if (temp < lower)
+ ((CELL *) tmp.rast)[col] = IS_COLD_CLOUD;
+ else
+ ((CELL *) tmp.rast)[col] = 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;
+ }
+ }
+ if (G_put_raster_row(tmp.fd, tmp.rast, CELL_TYPE) < 0)
+ {
+ G_fatal_error(_("Cannot write to <%s>"), tmp.name);
+ }
+ }
+
+ /* Finalización */
+
+ G_free(tmp.rast);
+ G_close_cell(tmp.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_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;
+}
+
+/**********************************************************
+ *
+ * Cloud shadows
+ *
+ **********************************************************/
+
+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;
+ }
+
+ return NO_DEFINED;
+}
Property changes on: grass-addons/i.landsat.acca/algorithm.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/i.landsat.acca/description.html
===================================================================
--- grass-addons/i.landsat.acca/description.html (rev 0)
+++ grass-addons/i.landsat.acca/description.html 2008-02-06 07:59:36 UTC (rev 29969)
@@ -0,0 +1,45 @@
+<H2>DESCRIPTION</H2>
+
+<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>
+
+
+
+<H2>NOTES</H2>
+
+<em>i.landsat.acca</em> use full raster and it do not take account the current region.
+
+<H2>EXAMPLES</H2>
+
+<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 -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>
+</ol>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.landsat.toar.html">i.landsat.toar</A><br>
+</em>
+
+
+
+<H2>AUTHOR</H2>
+
+E. Jorge Tizado (ej.tizado unileon es)<br>
+Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
+
+<p>
+<i>Last changed: $Date: 2007/07/07 00:00:00 $</i>
Added: grass-addons/i.landsat.acca/local_proto.h
===================================================================
--- grass-addons/i.landsat.acca/local_proto.h (rev 0)
+++ grass-addons/i.landsat.acca/local_proto.h 2008-02-06 07:59:36 UTC (rev 29969)
@@ -0,0 +1,44 @@
+#ifndef _LOCAL_PROTO_H
+#define _LOCAL_PROTO_H
+
+
+#define BAND2 0
+#define BAND3 1
+#define BAND4 2
+#define BAND5 3
+#define BAND6 4
+
+#define NO_CLOUD 0
+#define IS_CLOUD 1
+#define COLD_CLOUD 30
+#define WARM_CLOUD 50
+
+#define NO_DEFINED 1
+#define IS_SHADOW 2
+#define IS_COLD_CLOUD 6
+#define IS_WARM_CLOUD 9
+
+
+typedef struct
+{
+ 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_second(int, Gfile *, Gfile, int, double, double);
+
+int shadow_algorithm(double []);
+
+void filter_holes(int, Gfile *);
+
+void hist_put(double t, int hist[]);
+double mean(int hist[]);
+double quantile(double q, int hist[]);
+double moment(int n, int hist[]);
+
+#endif
Property changes on: grass-addons/i.landsat.acca/local_proto.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/i.landsat.acca/main.c
===================================================================
--- grass-addons/i.landsat.acca/main.c (rev 0)
+++ grass-addons/i.landsat.acca/main.c 2008-02-06 07:59:36 UTC (rev 29969)
@@ -0,0 +1,207 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.landsat.acca
+ *
+ * AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
+ *
+ * PURPOSE: Landsat TM/ETM+ Automatic Cloud Cover Assessment
+ *
+ * COPYRIGHT: (C) 2008 by the GRASS Development Team
+ *
+ * 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>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+extern int hist_n;
+
+/*----------------------------------------------*
+ * Constant threshold of ACCA algorithm
+ *----------------------------------------------*/
+
+extern double th_1;
+extern double th_1_b;
+extern double th_2[];
+extern double th_2_b;
+extern double th_3;
+extern double th_4;
+extern double th_4_b;
+extern double th_5;
+extern double th_6;
+extern double th_7;
+extern double th_8;
+
+/*----------------------------------------------*
+ *
+ * Check a raster name y return fd of open file
+ *
+ *----------------------------------------------*/
+int check_raster( char *raster_name )
+{
+ struct Cell_head cellhd;
+ RASTER_MAP_TYPE map_type;
+ int raster_fd;
+ char *mapset;
+
+ mapset = G_find_cell2(raster_name, "");
+ if (mapset == NULL) {
+ G_message("cell file [%s] not found", raster_name);
+ return -1;
+ }
+ if (G_legal_filename(raster_name) < 0) {
+ G_message("[%s] is an illegal name", raster_name);
+ return -1;
+ }
+ if ((raster_fd = G_open_cell_old(raster_name, mapset)) < 0) {
+ G_message("Cannot open cell file [%s]", raster_name);
+ return -1;
+ }
+ if (G_get_cellhd(raster_name, mapset, &cellhd) < 0) {
+ G_message("Cannot read file header of [%s]", raster_name);
+ return -1;
+ }
+ if (G_set_window(&cellhd) < 0) {
+ G_message("Unable to set region");
+ return -1;
+ }
+ if ((map_type = G_raster_map_type(raster_name, mapset)) != DCELL_TYPE) {
+ G_message("Map is not of DCELL_TYPE");
+ return -1;
+ }
+ return raster_fd;
+}
+
+/*----------------------------------------------*
+ *
+ * MAIN FUNCTION
+ *
+ *----------------------------------------------*/
+int main(int argc, char *argv[])
+{
+ struct History history;
+ struct GModule *module;
+
+ int i, verbose = 1;
+ struct Option *input, *output, *hist;
+ struct Flag *shadow, *sat5, *filter, *pass2;
+ char *in_name, *out_name;
+ double p;
+
+ Gfile band[5], out;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]);
+
+ /* initialize module */
+ module = G_define_module();
+ module->description = _("Landsat TM/ETM+ Automatic Cloud Cover Assessment (ACCA)");
+
+ input = G_define_option();
+ input->key = _("band_prefix");
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->gisprompt = _("input,cell,raster");
+ input->description = _("Base name of the landsat band rasters ([band_prefix].[band_number])");
+
+ output = G_define_option();
+ output->key = _("output");
+ output->type = TYPE_STRING;
+ output->required = YES;
+ output->gisprompt = _("input,cell,raster");
+ output->description = _("Output file name");
+
+ hist = G_define_option();
+ 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";
+
+ sat5 = G_define_flag();
+ sat5->key = '5';
+ sat5->description = _("Landsat-5 TM");
+ sat5->answer = 0;
+
+ filter = G_define_flag();
+ filter->key = 'f';
+ filter->description = _("Use final filter holes");
+ filter->answer = 0;
+
+ pass2 = G_define_flag();
+ pass2->key = '2';
+ pass2->description = _("With pass two processing");
+ pass2->answer = 0;
+
+ shadow = G_define_flag();
+ shadow->key = 's';
+ 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 );
+ if( hist_n < 10) hist_n = 10;
+
+ in_name = input->answer;
+
+ for (i = BAND2; i <= BAND6; i++)
+ {
+ 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 filename [%s]!"), band[i].name);
+ }
+ band[i].rast = G_allocate_raster_buf(DCELL_TYPE);
+ }
+
+ out_name = output->answer;
+
+ snprintf(out.name, 127, "%s", out_name);
+ if (G_legal_filename(out_name) < 0)
+ G_fatal_error(_("[%s] is an illegal name"), out.name);
+
+ /* --------------------------------------- */
+
+// if( sat5 -> answer )
+// {
+// th_4 = 205.;
+// }
+// acca_test(verbose, &out, band);
+ acca_algorithm(verbose, &out, band,
+ pass2->answer, shadow->answer);
+
+ if (filter->answer)
+ filter_holes(verbose, &out);
+ /* --------------------------------------- */
+
+ for (i = BAND2; i <= BAND6; i++)
+ {
+ G_free(band[i].rast);
+ G_close_cell(band[i].fd);
+ }
+
+// struct Categories cats;
+// G_read_raster_cats(out.name, char *mapset, cats)
+// G_write_raster_cats(out.name, &cats);
+
+ G_short_history(out.name, "raster", &history);
+ G_command_history(&history);
+ G_write_history(out.name, &history);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/i.landsat.acca/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/i.landsat.acca/tools.c
===================================================================
--- grass-addons/i.landsat.acca/tools.c (rev 0)
+++ grass-addons/i.landsat.acca/tools.c 2008-02-06 07:59:36 UTC (rev 29969)
@@ -0,0 +1,356 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+
+/*--------------------------------------------------------
+ HISTOGRAM ANALYSIS
+ Define un factor de escala = hist_n/100 con objeto
+ de dividir el entero 1 por 100/hist_n partes y
+ aumentar la precision.
+
+ Afecta al almacenamiento en el histograma pero
+ modifica el calculo de quantiles y momentos.
+ --------------------------------------------------------*/
+
+/* Global variable
+ allow use as parameters of command line */
+int hist_n = 100; /* interval 100/hist_n */
+
+void hist_put(double t, int hist[])
+{
+ int i;
+
+ /* scale factor */
+ i = (int)(t * ((double)hist_n/100.));
+
+ if (i < 1) i = 1;
+ if (i > hist_n) i = hist_n;
+
+ hist[i - 1] += 1;
+}
+
+/* Media real no del histograma */
+double mean(int hist[])
+{
+ int i, total;
+ double value, mean;
+
+ total = 0;
+ 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.));
+}
+
+double moment(int n, int hist[])
+{
+ int i, j, total;
+ double value, hmean, temp;
+
+ total = 0;
+ hmean = 0.;
+ for( i = 0; i < hist_n; i++ )
+ {
+ total += hist[i];
+ hmean += (double)(i * hist[i]);
+ }
+ hmean /= (double)total;
+
+ value = 0.;
+ for( i = 0; i < hist_n; i++ )
+ {
+ temp = 1.;
+ for( j = 0; j < n; j++ ) temp *= (i - hmean);
+ value += temp * (double)hist[i]/(double)total;
+ }
+
+ /* remove scale factor */
+ for( j = 0; j < n; j++ ) value /= ((double)hist_n/100.);
+ return value;
+}
+
+double quantile(double q, int hist[])
+{
+ int i, total;
+ double value, qmax, qmin;
+
+ total = 0;
+ 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;
+ }
+
+ /* remove scale factor */
+ return (value / ((double)hist_n/100.));
+}
+
+
+
+/*--------------------------------------------------------
+ FILTER HOLES OF CLOUDS
+ This a >=50% filter of 3x3
+ if >= 50% vecinos cloud then pixel set to cloud
+ --------------------------------------------------------*/
+
+#define UNDEFINED -1
+
+int pval(void * rast, int i)
+{
+ void * ptr = (void *)((CELL *)rast + i);
+ if( G_is_c_null_value(ptr) )
+ return 0;
+ else
+ return (int)((CELL *) rast)[i];
+}
+
+void filter_holes(int verbose, Gfile * out)
+{
+ int row, col, nrows, ncols;
+ char *mapset;
+
+ void *arast, *brast, *crast;
+ int i, pixel[9], cold, warm, shadow, nulo, lim;
+
+ Gfile tmp;
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ if( nrows < 3 || ncols < 3 )
+ return;
+
+ /* Open to read */
+ mapset = G_find_cell2(out->name, "");
+ if (mapset == NULL)
+ G_fatal_error("cell file [%s] not found", out->name);
+ arast = G_allocate_raster_buf(CELL_TYPE);
+ brast = G_allocate_raster_buf(CELL_TYPE);
+ crast = 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);
+
+ fprintf(stdout, "Filling cloud holes ... \n");
+
+
+ for (row = 0; row < nrows; row++)
+ {
+ if (verbose)
+ {
+ G_percent(row, nrows, 2);
+ }
+ /* Read row values */
+ if (row != 0)
+ {
+ if (G_get_c_raster_row(out->fd, arast, row - 1) < 0)
+ G_fatal_error(_("Could not read from <%s>"), out->name);
+ }
+ if (G_get_c_raster_row(out->fd, brast, row) < 0)
+ {
+ G_fatal_error(_("Could not read from <%s>"), out->name);
+ }
+ if (row != (nrows - 1))
+ {
+ if (G_get_c_raster_row(out->fd, crast, row + 1) < 0)
+ G_fatal_error(_("Could not read from <%s>"), out->name);
+ }
+ /* Analysis of all pixels */
+ for (col = 0; col < ncols; col++)
+ {
+ pixel[0] = pval(brast, col);
+ if( pixel[0] == 0 )
+ {
+ if( row == 0 )
+ {
+ pixel[1] = -1;
+ pixel[2] = -1;
+ pixel[3] = -1;
+ if( col == 0 )
+ {
+ pixel[4] = -1;
+ pixel[5] = pval(brast, col + 1);
+ pixel[6] = -1;
+ pixel[7] = pval(crast, col);
+ pixel[8] = pval(crast, col + 1);
+ }
+ else if( col != (ncols - 1) )
+ {
+ pixel[4] = pval(brast, col - 1);
+ pixel[5] = pval(brast, col + 1);
+ pixel[6] = pval(crast, col - 1);
+ pixel[7] = pval(crast, col);
+ pixel[8] = pval(crast, col + 1);
+ }
+ else
+ {
+ pixel[4] = pval(brast, col - 1);
+ pixel[5] = -1;
+ pixel[6] = pval(crast, col - 1);
+ pixel[7] = pval(crast, col);
+ pixel[8] = -1;
+ }
+ }
+ else if( row != (nrows - 1) )
+ {
+ if( col == 0 )
+ {
+ pixel[1] = -1;
+ pixel[2] = pval(arast, col);
+ pixel[3] = pval(arast, col + 1);
+ pixel[4] = -1;
+ pixel[5] = pval(brast, col + 1);
+ pixel[6] = -1;
+ pixel[7] = pval(crast, col);
+ pixel[8] = pval(crast, col + 1);
+ }
+ else if( col != (ncols - 1) )
+ {
+ pixel[1] = pval(arast, col - 1);
+ pixel[2] = pval(arast, col);
+ pixel[3] = pval(arast, col + 1);
+ pixel[4] = pval(brast, col - 1);
+ pixel[5] = pval(brast, col + 1);
+ pixel[6] = pval(crast, col - 1);
+ pixel[7] = pval(crast, col);
+ pixel[8] = pval(crast, col + 1);
+ }
+ else
+ {
+ pixel[1] = pval(arast, col - 1);
+ pixel[2] = pval(arast, col);
+ pixel[3] = -1;
+ pixel[4] = pval(brast, col - 1);
+ pixel[5] = -1;
+ pixel[6] = pval(crast, col - 1);
+ pixel[7] = pval(crast, col);
+ pixel[8] = -1;
+ }
+ }
+ else
+ {
+ pixel[6] = -1;
+ pixel[7] = -1;
+ pixel[8] = -1;
+ if( col == 0 )
+ {
+ pixel[1] = -1;
+ pixel[2] = pval(arast, col);
+ pixel[3] = pval(arast, col + 1);
+ pixel[4] = -1;
+ pixel[5] = pval(brast, col + 1);
+ }
+ else if( col != (ncols - 1) )
+ {
+ pixel[1] = pval(arast, col - 1);
+ pixel[2] = pval(arast, col);
+ pixel[3] = pval(arast, col + 1);
+ pixel[4] = pval(brast, col - 1);
+ pixel[5] = pval(brast, col + 1);
+ }
+ else
+ {
+ pixel[1] = pval(arast, col - 1);
+ pixel[2] = pval(arast, col);
+ pixel[3] = -1;
+ pixel[4] = pval(brast, col - 1);
+ pixel[5] = -1;
+ }
+ }
+
+// pixel[1] = (row == 0 || col == 0) ? -1 : pval(arast, col - 1);
+// pixel[2] = (row == 0) ? -1 : pval(arast, col);
+// pixel[3] = (row == 0 || col == (ncols-1)) ? -1 : pval(arast, col + 1);
+// pixel[4] = (col == 0) ? -1 : pval(brast, col - 1);
+// pixel[5] = (col == (ncols-1)) ? -1 : pval(brast, col + 1);
+// pixel[6] = (row == (nrows-1) || col == 0) ? -1 : pval(crast, col - 1);
+// pixel[7] = (row == (nrows-1)) ? -1 : pval(crast, col);
+// pixel[8] = (row == (nrows-1) || col == (ncols-1)) ? -1 : pval(crast, col + 1);
+
+ cold = warm = shadow = nulo = 0;
+ for( i = 1; i < 9; i++ )
+ {
+ switch( pixel[i] )
+ {
+ case IS_COLD_CLOUD: cold++; break;
+ case IS_WARM_CLOUD: warm++; break;
+ case IS_SHADOW : shadow++; break;
+ default: nulo++; break;
+ }
+ }
+ lim = (int)(cold + warm + shadow + nulo)/2;
+
+ /* Entra pixel[0] = 0 */
+ if( nulo < lim )
+ {
+ if( shadow >= (cold + warm) ) pixel[0] = IS_SHADOW;
+ else pixel[0] = (warm > cold) ? IS_WARM_CLOUD : IS_COLD_CLOUD;
+ }
+ }
+ if( pixel[0] != 0 )
+ {
+ ((CELL *) tmp.rast)[col] = pixel[0];
+ }
+ else
+ {
+ G_set_c_null_value((CELL *) tmp.rast + col, 1);
+ }
+ }
+ if (G_put_raster_row(tmp.fd, tmp.rast, CELL_TYPE) < 0)
+ {
+ G_fatal_error(_("Cannot write to <%s>"), tmp.name);
+ }
+ }
+
+ G_free(arast);
+ G_free(brast);
+ G_free(crast);
+ G_close_cell(out->fd);
+
+ G_free(tmp.rast);
+ G_close_cell(tmp.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_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;
+}
+
Property changes on: grass-addons/i.landsat.acca/tools.c
___________________________________________________________________
Name: svn:eol-style
+ native
More information about the grass-commit
mailing list