[GRASS-SVN] r44039 - in grass/trunk/imagery: . i.landsat.acca
i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Oct 25 17:34:01 EDT 2010
Author: martinl
Date: 2010-10-25 14:34:01 -0700 (Mon, 25 Oct 2010)
New Revision: 44039
Added:
grass/trunk/imagery/i.landsat.acca/
grass/trunk/imagery/i.landsat.acca/Makefile
grass/trunk/imagery/i.landsat.acca/algorithm.c
grass/trunk/imagery/i.landsat.acca/i.landsat.acca.html
grass/trunk/imagery/i.landsat.acca/local_proto.h
grass/trunk/imagery/i.landsat.acca/main.c
grass/trunk/imagery/i.landsat.acca/tools.c
Modified:
grass/trunk/imagery/i.landsat.toar/
grass/trunk/imagery/i.landsat.toar/main.c
Log:
i.landsat.acca from addons added to trunk
Property changes on: grass/trunk/imagery/i.landsat.acca
___________________________________________________________________
Added: svn:ignore
+ OBJ.*
Added: grass/trunk/imagery/i.landsat.acca/Makefile
===================================================================
--- grass/trunk/imagery/i.landsat.acca/Makefile (rev 0)
+++ grass/trunk/imagery/i.landsat.acca/Makefile 2010-10-25 21:34:01 UTC (rev 44039)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.landsat.acca
+
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass/trunk/imagery/i.landsat.acca/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass/trunk/imagery/i.landsat.acca/algorithm.c
===================================================================
--- grass/trunk/imagery/i.landsat.acca/algorithm.c (rev 0)
+++ grass/trunk/imagery/i.landsat.acca/algorithm.c 2010-10-25 21:34:01 UTC (rev 44039)
@@ -0,0 +1,472 @@
+/* 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>
+#include <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+#define SCALE 200.
+#define K_BASE 230.
+
+
+/* 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 DSTD 4
+
+
+/**********************************************************
+ *
+ * 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 Threshold */
+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 Threshold */
+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;
+
+void acca_algorithm(Gfile * out, Gfile band[],
+ int single_pass, int with_shadow, int cloud_signature)
+{
+ int i, count[5], hist_cold[hist_n], hist_warm[hist_n];
+ double max, value[5], signa[5], idesert, review_warm, 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;
+ }
+
+ /* FIRST FILTER ... */
+ acca_first(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] / ((double)count[SOIL]));
+
+ /* BAND-6 CLOUD SIGNATURE DEVELOPMENT */
+ if (idesert <= .5 || value[SNOW] > 0.01) {
+ /* Only the cold clouds are used
+ if snow or desert soil is present */
+ review_warm = 1;
+ }
+ else {
+ /* The cold and warm clouds are combined
+ and treated as a single population */
+ review_warm = 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] = SCALE * signa[SUM_COLD] / ((double)count[COLD]);
+ signa[COVER] = ((double)count[COLD]) / ((double)count[TOTAL]);
+
+ G_message(_("Preliminary scene analysis:"));
+ G_message(_("* Desert index: %.2lf"), idesert);
+ G_message(_("* Snow cover: %.2lf %%"), 100. * value[SNOW]);
+ G_message(_("* Cloud cover: %.2lf %%"), 100. * signa[COVER]);
+ G_message(_("* Temperature of clouds:"));
+ G_message(_("** Maximum: %.2lf "), signa[KMAX]);
+ G_message(_("** Mean (%s cloud) : %.2lf K"),
+ (review_warm ? "cold" : "all"), signa[KMEAN]);
+ G_message(_("** Minimum: %.2lf K"), 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.)) {
+ G_message(_("Histogram cloud signature:"));
+
+ 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);
+
+ G_message(_("* Mean temperature: %.2lf K"), value[MEAN]);
+ G_message(_("* Standard deviation: %.2lf"), value[DSTD]);
+ G_message(_("* Skewness: %.2lf"), value[SKEW]);
+ G_message(_("* Histogram classes: %d"), hist_n);
+
+ shift = value[SKEW];
+ if (shift > 1.)
+ shift = 1.;
+ else if (shift < 0.)
+ shift = 0.;
+
+ 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;
+
+ G_message(_("* 98.75 percentile: %.2lf K"), max);
+ G_message(_("* 97.50 percentile: %.2lf K"), value[KUPPER]);
+ G_message(_("* 83.50 percentile: %.2lf K"), value[KLOWER]);
+
+ /* step 17 & 18 */
+ if (shift > 0.) {
+ shift *= value[DSTD];
+
+ if ((value[KUPPER] + shift) > max) {
+ value[KLOWER] += (max - value[KUPPER]);
+ value[KUPPER] = max;
+ }
+ else {
+ value[KLOWER] += shift;
+ value[KUPPER] += shift;
+ }
+ }
+
+ G_message(_("Maximum temperature:"));
+ G_message(_("* Cold cloud: %.2lf K"), value[KUPPER]);
+ G_message(_("* Warn cloud: %.2lf K"), value[KLOWER]);
+ }
+ else {
+ if (signa[KMEAN] < 295.) {
+ /* Retained warm and cold clouds */
+ G_message(_("Result: Scene with clouds"));
+ review_warm = 0;
+ value[KUPPER] = 0.;
+ value[KLOWER] = 0.;
+ }
+ else {
+ /* Retained cold clouds */
+ G_message(_("Result: Scene cloud free"));
+ review_warm = 1;
+ value[KUPPER] = 0.;
+ value[KLOWER] = 0.;
+ }
+ }
+
+ /* SECOND FILTER ... */
+ /* By-pass two processing but it retains warm and cold clouds */
+ if (single_pass == TRUE) {
+ review_warm = -1.;
+ value[KUPPER] = 0.;
+ value[KLOWER] = 0.;
+ }
+ acca_second(out, band[BAND6],
+ review_warm, value[KUPPER], value[KLOWER]);
+ /* CATEGORIES: IS_WARM_CLOUD, IS_COLD_CLOUD, IS_SHADOW, NULL (= NO_CLOUD) */
+
+ return;
+}
+
+
+void acca_first(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;
+
+ /* Creation of output file */
+ out->rast = Rast_allocate_buf(CELL_TYPE);
+ if ((out->fd = Rast_open_new(out->name, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), out->name);
+
+ /* ----- ----- */
+ G_message(_("Processing first pass..."));
+
+ stats[SUM_COLD] = 0.;
+ stats[SUM_WARM] = 0.;
+ stats[KMAX] = 0.;
+ stats[KMIN] = 10000.;
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ for (i = BAND2; i <= BAND6; i++) {
+ Rast_get_d_row(band[i].fd, band[i].rast, row);
+ }
+ for (col = 0; col < ncols; col++) {
+ code = NO_DEFINED;
+ /* Null when null pixel in any band */
+ for (i = BAND2; i <= BAND6; i++) {
+ if (Rast_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]);
+ /* ----------------------------------------------------- */
+ /* step 1. Brightness Threshold: Eliminates dark images */
+ if (pixel[BAND3] > th_1) {
+ /* step 3. Normalized Snow Difference Index: Eliminates many types of snow */
+ if (nsdi > th_2[0] && nsdi < th_2[1]) {
+ /* step 5. Temperature Threshold: Eliminates warm image features */
+ if (pixel[BAND6] < th_3) {
+ rat56 = (1. - pixel[BAND5]) * pixel[BAND6];
+ /* step 6. Band 5/6 Composite: Eliminates numerous categories including ice */
+ if (rat56 < th_4) {
+ /* 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 = NO_DEFINED;
+ }
+ }
+ else {
+ code = NO_DEFINED;
+ count[SOIL]++;
+ }
+ }
+ else {
+ code = NO_DEFINED;
+ }
+ }
+ else {
+ /* step 7 */
+ code =
+ (pixel[BAND5] <
+ th_4_b) ? NO_CLOUD : NO_DEFINED;
+ }
+ }
+ else {
+ code = NO_CLOUD;
+ }
+ }
+ 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;
+ }
+ /* ----------------------------------------------------- */
+ }
+ if (code == NO_CLOUD) {
+ Rast_set_c_null_value((CELL *) out->rast + col, 1);
+ }
+ else {
+ ((CELL *) out->rast)[col] = code;
+ }
+ }
+ Rast_put_row(out->fd, out->rast, CELL_TYPE);
+ }
+ G_percent(1, 1, 1);
+
+ G_free(out->rast);
+ Rast_close(out->fd);
+
+ return;
+}
+
+
+void acca_second(Gfile * out, Gfile band,
+ int review_warm, double upper, double lower)
+{
+ int row, col, nrows, ncols;
+
+ int code;
+ double temp;
+ Gfile tmp;
+
+ /* Open to read */
+ if ((out->fd = Rast_open_old(out->name, "")) < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), out->name);
+
+ out->rast = Rast_allocate_buf(CELL_TYPE);
+
+ /* Open to write */
+ sprintf(tmp.name, "_%d.BBB", getpid());
+ tmp.rast = Rast_allocate_buf(CELL_TYPE);
+ if ((tmp.fd = Rast_open_new(tmp.name, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), tmp.name);
+
+ if (upper == 0.)
+ G_message(_("Removing ambiguous pixels..."));
+ else
+ G_message(_("Pass two processing..."));
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+
+ Rast_get_d_row(band.fd, band.rast, row);
+ Rast_get_c_row(out->fd, out->rast, row);
+
+ for (col = 0; col < ncols; col++) {
+ if (Rast_is_c_null_value((void *)((CELL *) out->rast + col))) {
+ Rast_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 && review_warm == 1)) {
+ temp = (double)((DCELL *) band.rast)[col];
+ if (temp > upper) {
+ Rast_set_c_null_value((CELL *) tmp.rast + col, 1);
+ }
+ else {
+ ((CELL *) tmp.rast)[col] =
+ (temp < lower) ? IS_WARM_CLOUD : IS_COLD_CLOUD;
+ }
+ }
+ else
+ /* Join warm (not ambiguous) and cold clouds */
+ if (code == COLD_CLOUD || code == WARM_CLOUD) {
+ ((CELL *) tmp.rast)[col] = (code == WARM_CLOUD &&
+ review_warm ==
+ 0) ? IS_WARM_CLOUD :
+ IS_COLD_CLOUD;
+ }
+ else
+ ((CELL *) tmp.rast)[col] = IS_SHADOW;
+ }
+ }
+ Rast_put_row(tmp.fd, tmp.rast, CELL_TYPE);
+ }
+ G_percent(1, 1, 1);
+
+ G_free(tmp.rast);
+ Rast_close(tmp.fd);
+
+ G_free(out->rast);
+ Rast_close(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[])
+{
+ /* 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)
+ /*
+ if (pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. &&
+ pixel[BAND4] / pixel[BAND2] > 1.)
+ */
+ {
+ return IS_SHADOW;
+ }
+
+ return NO_DEFINED;
+}
Property changes on: grass/trunk/imagery/i.landsat.acca/algorithm.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/imagery/i.landsat.acca/i.landsat.acca.html
===================================================================
--- grass/trunk/imagery/i.landsat.acca/i.landsat.acca.html (rev 0)
+++ grass/trunk/imagery/i.landsat.acca/i.landsat.acca.html 2010-10-25 21:34:01 UTC (rev 44039)
@@ -0,0 +1,61 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.landsat.acca</em> implements the <b>Automated Cloud-Cover
+Assessment</B> (ACCA) Algorithm from Irish (2000) with the constant
+values for pass filter one from Irish et al. (2006). To do this, it
+needs Landsat 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
+with <em><a href="i.landsat.toar.html">i.landsat.toar</a></em>).
+
+<p>
+The ACCA algorithm 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 for Landsat-7 ETM+ but because reflectance is used it is also
+useful for Landsat-4/5 TM.
+
+<h2>NOTES</h2>
+
+<em>i.landsat.acca</em> works in the current region settings.
+
+<h2>EXAMPLES</h2>
+
+Run the standard ACCA algorithm with filling of small cloud holes
+(the <b>-f</b> flag): 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>
+
+
+<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>
+</em>
+
+<h2>AUTHOR</h2>
+
+E. Jorge Tizado (ej.tizado unileon es), Dept. Biodiversity and Environmental Management, University of León, Spain
+
+<p>
+<i>Last changed: $Date$</i>
Property changes on: grass/trunk/imagery/i.landsat.acca/i.landsat.acca.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/imagery/i.landsat.acca/local_proto.h
===================================================================
--- grass/trunk/imagery/i.landsat.acca/local_proto.h (rev 0)
+++ grass/trunk/imagery/i.landsat.acca/local_proto.h 2010-10-25 21:34:01 UTC (rev 44039)
@@ -0,0 +1,43 @@
+#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[GNAME_MAX];
+
+} Gfile;
+
+
+void acca_algorithm(Gfile *, Gfile[], int, int, int);
+void acca_first(Gfile *, Gfile[], int, int[], int[], int[], double[]);
+void acca_second(Gfile *, Gfile, int, double, double);
+
+int shadow_algorithm(double[]);
+
+void filter_holes(Gfile *);
+
+void hist_put(double t, int hist[]);
+double quantile(double q, int hist[]);
+double moment(int n, int hist[], int k);
+
+#endif
Property changes on: grass/trunk/imagery/i.landsat.acca/local_proto.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass/trunk/imagery/i.landsat.acca/main.c
===================================================================
--- grass/trunk/imagery/i.landsat.acca/main.c (rev 0)
+++ grass/trunk/imagery/i.landsat.acca/main.c 2010-10-25 21:34:01 UTC (rev 44039)
@@ -0,0 +1,222 @@
+
+/****************************************************************************
+ *
+ * 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/raster.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)
+{
+ RASTER_MAP_TYPE map_type;
+ int raster_fd;
+
+ if ((raster_fd = Rast_open_old(raster_name, "")) < 0) {
+ G_fatal_error(_("Unable to open raster map <%s>"), raster_name);
+ }
+ /* Uncomment to work in full raster map
+ if (G_get_cellhd(raster_name, mapset, &cellhd) < 0) {
+ G_warning(_("Unable to read header of raster map <%s>"), raster_name);
+ return -1;
+ }
+ if (G_set_window(&cellhd) < 0) {
+ G_warning(_("Cannot reset current region"));
+ return -1;
+ }
+ */
+ if ((map_type = Rast_get_map_type(raster_fd)) != DCELL_TYPE) {
+ G_fatal_error(_("Input raster map <%s> is not floating point (process DN using i.landsat.toar to radiance first)"), raster_name);
+ }
+
+ return raster_fd;
+}
+
+/*----------------------------------------------*
+ *
+ * MAIN FUNCTION
+ *
+ *----------------------------------------------*/
+int main(int argc, char *argv[])
+{
+ struct History history;
+ struct GModule *module;
+
+ int i;
+ struct Option *band_prefix, *output, *hist, *b56c, *b45r;
+ struct Flag *shadow, *filter, *sat5, *pass2, *csig;
+ char *in_name, *out_name;
+ struct Categories cats;
+
+ CELL cell_shadow = IS_SHADOW, cell_cold_cloud = IS_COLD_CLOUD, cell_warm_cloud = IS_WARM_CLOUD;
+ Gfile band[5], out;
+
+ char title[1024];
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]);
+
+ /* initialize module */
+ module = G_define_module();
+ module->description =
+ _("Landsat TM/ETM+ Automatic Cloud Cover Assessment (ACCA).");
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("landsat"));
+ G_add_keyword(_("acca"));
+
+ band_prefix = G_define_option();
+ band_prefix->key = "band_prefix";
+ band_prefix->label = _("Base name of input raster bands");
+ band_prefix->description = _("Example: 'B.' for B.1, B.2, ...");
+ band_prefix->type = TYPE_STRING;
+ band_prefix->required = YES;
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+
+ b56c = G_define_option();
+ b56c->key = "b56composite";
+ b56c->type = TYPE_DOUBLE;
+ b56c->required = NO;
+ 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;
+ hist->required = NO;
+ hist->description =
+ _("Number of classes in the cloud temperature histogram");
+ hist->answer = "100";
+
+ sat5 = G_define_flag();
+ sat5->key = '5';
+ sat5->label = _("Data is Landsat-5 TM");
+ sat5->description = _("I.e. Thermal band is '.6' not '.61')");
+
+ filter = G_define_flag();
+ filter->key = 'f';
+ filter->description =
+ _("Apply post-processing filter to remove small holes");
+
+ csig = G_define_flag();
+ csig->key = 'x';
+ csig->description = _("Always use cloud signature (step 14)");
+
+ pass2 = G_define_flag();
+ pass2->key = '2';
+ pass2->description =
+ _("Bypass second-pass processing, and merge warm (not ambiguous) and cold clouds");
+
+ shadow = G_define_flag();
+ shadow->key = 's';
+ shadow->description = _("Include a category for cloud shadows");
+
+ 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 = band_prefix->answer;
+
+ for (i = BAND2; i <= BAND6; i++) {
+ sprintf(band[i].name, "%s%d%c", in_name, i + 2,
+ (i == BAND6 && !sat5->answer ? '1' : '\0'));
+ band[i].fd = check_raster(band[i].name);
+ band[i].rast = Rast_allocate_buf(DCELL_TYPE);
+ }
+
+ out_name = output->answer;
+
+ sprintf(out.name, "%s", out_name);
+ if (G_legal_filename(out_name) < 0)
+ G_fatal_error(_("<%s> is an illegal file name"), out.name);
+
+ /* --------------------------------------- */
+ th_4 = atof(b56c->answer);
+ th_7 = atof(b45r->answer);
+ acca_algorithm(&out, band, pass2->answer, shadow->answer,
+ csig->answer);
+
+ if (filter->answer)
+ filter_holes(&out);
+ /* --------------------------------------- */
+
+ for (i = BAND2; i <= BAND6; i++) {
+ G_free(band[i].rast);
+ Rast_close(band[i].fd);
+ }
+
+ /* write out map title and category labels */
+ Rast_init_cats("", &cats);
+ sprintf(title, "LANDSAT-%s Automatic Cloud Cover Assessment",
+ sat5->answer ? "5 TM" : "7 ETM+");
+ Rast_set_cats_title(title, &cats);
+
+ Rast_set_c_cat(&cell_shadow, &cell_shadow,
+ "Shadow", &cats);
+ Rast_set_c_cat(&cell_cold_cloud, &cell_cold_cloud,
+ "Cold cloud", &cats);
+ Rast_set_c_cat(&cell_warm_cloud, &cell_warm_cloud,
+ "Warm cloud", &cats);
+
+ Rast_write_cats(out.name, &cats);
+ Rast_free_cats(&cats);
+
+ /* write out command line opts */
+ Rast_short_history(out.name, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(out.name, &history);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/imagery/i.landsat.acca/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/imagery/i.landsat.acca/tools.c
===================================================================
--- grass/trunk/imagery/i.landsat.acca/tools.c (rev 0)
+++ grass/trunk/imagery/i.landsat.acca/tools.c 2010-10-25 21:34:01 UTC (rev 44039)
@@ -0,0 +1,303 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/raster.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 parameter in the command line */
+int hist_n = 100; /* interval of real data 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;
+}
+
+/* histogram moment */
+double moment(int n, int hist[], int k)
+{
+ int i, total;
+ double value, mean;
+
+ k = 0;
+
+ total = 0;
+ mean = 0.;
+ for (i = 0; i < hist_n; i++) {
+ total += hist[i];
+ mean += (double)(i * hist[i]);
+ }
+ mean /= ((double)total); /* histogram mean */
+
+ value = 0.;
+ for (i = 0; i < hist_n; i++) {
+ value += (pow((i - mean), n) * ((double)hist[i]));
+ }
+ value /= (double)(total - k);
+
+ return (value / pow((double)hist_n / 100., n) );
+}
+
+/* Real data quantile */
+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];
+ }
+
+ value = 0;
+ 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
+ --------------------------------------------------------*/
+
+int pval(void *rast, int i)
+{
+ void *ptr = (void *)((CELL *) rast + i);
+
+ if (Rast_is_c_null_value(ptr))
+ return 0;
+ else
+ return (int)((CELL *) rast)[i];
+}
+
+void filter_holes(Gfile * out)
+{
+ int row, col, nrows, ncols;
+
+ void *arast, *brast, *crast;
+ int i, pixel[9], cold, warm, shadow, nulo, lim;
+
+ Gfile tmp;
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ if (nrows < 3 || ncols < 3)
+ return;
+
+ /* Open to read */
+ if ((out->fd = Rast_open_old(out->name, "")) < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), out->name);
+
+ arast = Rast_allocate_buf(CELL_TYPE);
+ brast = Rast_allocate_buf(CELL_TYPE);
+ crast = Rast_allocate_buf(CELL_TYPE);
+
+ /* Open to write */
+ sprintf(tmp.name, "_%d.BBB", getpid());
+ tmp.rast = Rast_allocate_buf(CELL_TYPE);
+ if ((tmp.fd = Rast_open_new(tmp.name, CELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), tmp.name);
+
+ G_message(_("Filling small holes in clouds ..."));
+
+ /* Se puede acelerar creandolos nulos y luego arast = brast
+ brast = crast y cargando crast solamente
+ G_set_f_null_value(cell[2], ncols);
+ */
+
+ for (row = 0; row < nrows; row++) {
+ /* Read row values */
+ if (row != 0) {
+ Rast_get_c_row(out->fd, arast, row - 1);
+ Rast_get_c_row(out->fd, brast, row);
+ }
+ if (row != (nrows - 1)) {
+ Rast_get_c_row(out->fd, crast, row + 1);
+ }
+ /* 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;
+ }
+ }
+
+ 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 {
+ Rast_set_c_null_value((CELL *) tmp.rast + col, 1);
+ }
+ }
+ Rast_put_row(tmp.fd, tmp.rast, CELL_TYPE);
+ G_percent(row, nrows, 2);
+ }
+
+ G_free(arast);
+ G_free(brast);
+ G_free(crast);
+ Rast_close(out->fd);
+
+ G_free(tmp.rast);
+ Rast_close(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/trunk/imagery/i.landsat.acca/tools.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Property changes on: grass/trunk/imagery/i.landsat.toar
___________________________________________________________________
Added: svn:ignore
+ OBJ.*
Modified: grass/trunk/imagery/i.landsat.toar/main.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/main.c 2010-10-24 20:26:17 UTC (rev 44038)
+++ grass/trunk/imagery/i.landsat.toar/main.c 2010-10-25 21:34:01 UTC (rev 44039)
@@ -40,9 +40,9 @@
RASTER_MAP_TYPE in_data_type;
- struct Option *band_prefix, *output_suffix, *metfn, *sensor, *adate, *pdate, *elev,
+ struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate, *pdate, *elev,
*bgain, *metho, *perc, *dark, *satz, *atmo;
- char *basename, *met, *suffixname, *sensorname;
+ char *inputname, *met, *outputname, *sensorname;
struct Flag *msss, *frad, *l5_mtl;
lsat_data lsat;
@@ -68,19 +68,19 @@
G_add_keyword(_("dos-type simple atmospheric correction"));
/* It defines the different parameters */
- band_prefix = G_define_option();
- band_prefix->key = "band_prefix";
- band_prefix->label = _("Base name of input raster bands");
- band_prefix->description = _("Example: 'B.' for B.1, B.2, ...");
- band_prefix->type = TYPE_STRING;
- band_prefix->required = YES;
+ input_prefix = G_define_option();
+ input_prefix->key = "input_prefix";
+ input_prefix->label = _("Base name of input raster bands");
+ input_prefix->description = _("Example: 'B.' for B.1, B.2, ...");
+ input_prefix->type = TYPE_STRING;
+ input_prefix->required = YES;
- output_suffix = G_define_option();
- output_suffix->key = "output_suffix";
- output_suffix->label = _("Suffix for output raster maps");
- output_suffix->description = _("Example: '_toar' generates B.1_toar, B.2_toar, ...");
- output_suffix->type = TYPE_STRING;
- output_suffix->required = YES;
+ output_prefix = G_define_option();
+ output_prefix->key = "output_prefix";
+ output_prefix->label = _("Prefix for output raster maps");
+ output_prefix->description = _("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
+ output_prefix->type = TYPE_STRING;
+ output_prefix->required = YES;
metfn = G_define_standard_option(G_OPT_F_INPUT);
metfn->key = "metfile";
@@ -206,8 +206,8 @@
* Stores options and flag to variables
*****************************************/
met = metfn->answer;
- basename = band_prefix->answer;
- suffixname = output_suffix->answer;
+ inputname = input_prefix->answer;
+ outputname = output_prefix->answer;
sensorname = sensor -> answer ? sensor->answer: "";
G_zero(&lsat, sizeof(lsat));
@@ -340,7 +340,7 @@
for (j = 0; j < 256; j++)
hist[j] = 0L;
- sprintf(band_in, "%s%d", basename, lsat.band[i].code);
+ sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
if ((infd = Rast_open_old(band_in, "")) < 0)
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
Rast_get_cellhd(band_in, "", &cellhd);
@@ -466,8 +466,8 @@
G_message(_("Calculating..."));
for (i = 0; i < lsat.bands; i++) {
- sprintf(band_in, "%s%d", basename, lsat.band[i].code);
- sprintf(band_out, "%s%d%s", basename, lsat.band[i].code, suffixname);
+ sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
+ sprintf(band_out, "%s%d", outputname, lsat.band[i].code);
if ((infd = Rast_open_old(band_in, "")) < 0)
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
More information about the grass-commit
mailing list