[GRASS-SVN] r43728 - grass-addons/imagery/i.landsat.acca
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Sep 29 05:34:05 EDT 2010
Author: ejtizado
Date: 2010-09-29 09:34:05 +0000 (Wed, 29 Sep 2010)
New Revision: 43728
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:
more standard ACCA processing and new/change flags
Modified: grass-addons/imagery/i.landsat.acca/algorithm.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/algorithm.c 2010-09-29 00:58:17 UTC (rev 43727)
+++ grass-addons/imagery/i.landsat.acca/algorithm.c 2010-09-29 09:34:05 UTC (rev 43728)
@@ -74,10 +74,10 @@
#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 cloud_signature)
{
int i, count[5], hist_cold[hist_n], hist_warm[hist_n];
- double x, value[5], signa[5], idesert, warm_ambiguous, shift;
+ double max, value[5], signa[5], idesert, review_warm, shift;
/* Reset variables ... */
for (i = 0; i < 5; i++) {
@@ -100,20 +100,18 @@
value[0] = (double)(count[WARM] + count[COLD]);
idesert =
- (value[0] == 0. ? 0. : 1. / (1. + (double)count[SOIL]) / value[0]);
+ (value[0] == 0. ? 0. : 1. / (1. + ((double)count[SOIL]) / value[0]));
/* BAND-6 CLOUD SIGNATURE DEVELOPMENT */
- if (value[SNOW] > 0.01)
- // if (idesert <= .5 || 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;
+ review_warm = 1;
}
else {
/* The cold and warm clouds are combined
and treated as a single population */
- warm_ambiguous = 0;
+ review_warm = 0;
count[COLD] += count[WARM];
value[COLD] += value[WARM];
signa[SUM_COLD] += signa[SUM_WARM];
@@ -124,62 +122,68 @@
signa[KMEAN] = (signa[SUM_COLD] / (double)count[COLD]) * SCALE;
signa[COVER] = (double)count[COLD] / (double)count[TOTAL];
- 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]);
+ G_message("> PRELIMINARY SCENE ANALYSIS\n");
+ G_message("> Desert index: %.3lf", idesert);
+ G_message("> Snow cover: %.3lf %%", 100. * value[SNOW]);
+ G_message("> Cloud cover: %.3lf %%", 100. * signa[COVER]);
+ G_message("> Temperature of clouds");
+ G_message("> Maximum: %.2lf K\n", signa[KMAX]);
+ G_message("> Mean (%s cloud) : %.2lf K\n",
+ (review_warm ? "cold" : "all"), signa[KMEAN]);
+ G_message("> Minimum: %.2lf K\n", signa[KMIN]);
- /* WARNING: variable 'value' reutilizada con nuevos valores */
+ /* WARNING: re-use of the variable 'value' with new meaning */
- 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 ): */
+ /* step 14 */
+ if (cloud_signature ||
+ (idesert <= .5 && signa[COVER] > 0.004 && signa[KMEAN] < 295.)) {
+ value[MEAN] = quantile(0.5, hist_cold) + K_BASE; /* mean(hist_cold) + K_BASE; */
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];
+ shift = value[SKEW];
+ if (shift > 1.)
+ shift = 1.;
+ else if (shift < 0.)
+ shift = 0.;
- x = quantile(0.9875, hist_cold) + K_BASE;
- if ((value[KUPPER] + shift) > x) {
- value[KUPPER] = x;
- value[KLOWER] = x - shift; /* ??? COMPROBAR */
+ shift *= value[DSTD];
+
+ max = quantile(0.9875, hist_cold) + K_BASE;
+ value[KUPPER] = quantile(0.975, hist_cold) + K_BASE;
+ value[KLOWER] = quantile(0.835, hist_cold) + K_BASE;
+ /* step 17 & 18 */
+ if (shift > 0.) {
+ if ((value[KUPPER] * shift) > max) {
+ value[KLOWER] += (value[KUPPER] * (shift - 1.));
+ value[KUPPER] = max;
}
else {
- value[KUPPER] += shift;
- value[KLOWER] += shift;
+ value[KLOWER] *= shift;
+ value[KUPPER] *= shift;
}
}
- 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]);
+ 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("> Cold cloud: maximun %.2lf K\n", value[KUPPER]);
+ G_message("> Warn cloud: maximun %.2lf K\n", value[KLOWER]);
}
else {
if (signa[KMEAN] < 295.) {
/* Retained warm and cold clouds */
G_message(" Scene with clouds\n");
- warm_ambiguous = 0;
+ review_warm = 0;
value[KUPPER] = 0.;
value[KLOWER] = 0.;
}
else {
/* Retained cold clouds */
G_message(" Scene cloud free\n");
- warm_ambiguous = 1;
+ review_warm = 1;
value[KUPPER] = 0.;
value[KLOWER] = 0.;
}
@@ -187,14 +191,14 @@
/* SECOND FILTER ... */
/* By-pass two processing but it retains warm and cold clouds */
- if (two_pass == 0) {
- warm_ambiguous = 0;
+ if (two_pass == 1) {
+ review_warm = -1.;
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) */
+ review_warm, value[KUPPER], value[KLOWER]);
+ /* CATEGORIES: IS_WARM_CLOUD, IS_COLD_CLOUD, IS_SHADOW, NULL (= NO_CLOUD) */
return;
}
@@ -212,10 +216,10 @@
/* 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);
+ G_fatal_error(_("Unable to create raster map <%s>"), out->name);
/* ----- ----- */
- G_message("Pass one processing ... \n");
+ G_message("Pass one processing ...");
stats[SUM_COLD] = 0.;
stats[SUM_WARM] = 0.;
@@ -226,13 +230,10 @@
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);
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ band[i].name, row);
}
for (col = 0; col < ncols; col++) {
code = NO_DEFINED;
@@ -333,7 +334,10 @@
}
}
if (G_put_raster_row(out->fd, out->rast, CELL_TYPE) < 0) {
- G_fatal_error(_("Cannot write row to <%s>"), out->name);
+ G_fatal_error(_("Failed writing raster map <%s> row %d"),
+ out->name, row);
+
+ G_percent(row, nrows, 2);
}
}
@@ -345,7 +349,7 @@
void acca_second(int verbose, Gfile * out, Gfile band,
- int warm_ambiguous, double upper, double lower)
+ int review_warm, double upper, double lower)
{
int i, row, col, nrows, ncols;
char *mapset;
@@ -357,21 +361,21 @@
/* Open to read */
mapset = G_find_cell2(out->name, "");
if (mapset == NULL)
- G_fatal_error("cell file [%s] not found", out->name);
+ G_fatal_error(_("Raster map <%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);
+ G_fatal_error(_("Unable to open raster map <%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);
+ G_fatal_error(_("Unable to create raster map <%s>"), tmp.name);
if (upper == 0.)
- G_message("Removing ambiguous pixels ... \n");
+ G_message(_("Removing ambiguous pixels ..."));
else
- G_message("Pass two processing ... \n");
+ G_message(_("Pass two processing ..."));
nrows = G_window_rows();
ncols = G_window_cols();
@@ -381,9 +385,11 @@
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);
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ band.name, row);
if (G_get_c_raster_row(out->fd, out->rast, row) < 0)
- G_fatal_error(_("Could not read from <%s>"), out->name);
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ out->name, row);
for (col = 0; col < ncols; col++) {
if (G_is_c_null_value((void *)((CELL *) out->rast + col))) {
@@ -393,28 +399,30 @@
code = (int)((CELL *) out->rast)[col];
/* Resolve ambiguous pixels */
if (code == NO_DEFINED ||
- (code == WARM_CLOUD) && warm_ambiguous == 1) {
+ (code == WARM_CLOUD && review_warm == 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;
+ (temp < lower) ? IS_WARM_CLOUD : IS_COLD_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;
+ /* Joint 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;
}
}
if (G_put_raster_row(tmp.fd, tmp.rast, CELL_TYPE) < 0) {
- G_fatal_error(_("Cannot write to <%s>"), tmp.name);
+ G_fatal_error(_("Cannot write to raster map <%s>"), tmp.name);
}
}
@@ -449,9 +457,15 @@
int shadow_algorithm(double pixel[])
{
- 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
- ) {
+ /*
+ if (pixel[BAND3] < 0.07 && (1 - pixel[BAND4]) * pixel[BAND6] > 240. &&
+ pixel[BAND4] / pixel[BAND2] > 1.) {
+ */
+ /* 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) {
return IS_SHADOW;
}
Modified: grass-addons/imagery/i.landsat.acca/description.html
===================================================================
--- grass-addons/imagery/i.landsat.acca/description.html 2010-09-29 00:58:17 UTC (rev 43727)
+++ grass-addons/imagery/i.landsat.acca/description.html 2010-09-29 09:34:05 UTC (rev 43728)
@@ -1,30 +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>
+<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>).
-<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>
+<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.
-
-
<H2>NOTES</H2>
-<em>i.landsat.acca</em> use full raster and it do not take account the current region.
+<em>i.landsat.acca</em> works in the current region settings.
<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>
+To run the standard ACCA algorithm with 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>:
<div class="code"><pre>
-i.landsat.acca -f2 band=226_62.toar out=226_62.acca
+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>
+ <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>
@@ -39,7 +54,8 @@
<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>
+<i>Last changed: $Date: 2010/09/29 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-29 00:58:17 UTC (rev 43727)
+++ grass-addons/imagery/i.landsat.acca/local_proto.h 2010-09-29 09:34:05 UTC (rev 43728)
@@ -28,7 +28,7 @@
} Gfile;
-void acca_algorithm(int, Gfile *, Gfile[], int, int);
+void acca_algorithm(int, Gfile *, Gfile[], int, int, int);
void acca_first(int, Gfile *, Gfile[], int, int[], int[], int[], double[]);
void acca_second(int, Gfile *, Gfile, int, double, double);
Modified: grass-addons/imagery/i.landsat.acca/main.c
===================================================================
--- grass-addons/imagery/i.landsat.acca/main.c 2010-09-29 00:58:17 UTC (rev 43727)
+++ grass-addons/imagery/i.landsat.acca/main.c 2010-09-29 09:34:05 UTC (rev 43728)
@@ -55,27 +55,29 @@
mapset = G_find_cell2(raster_name, "");
if (mapset == NULL) {
- G_message("cell file [%s] not found", raster_name);
+ G_warning(_("Raster map <%s> not found"), raster_name);
return -1;
}
if (G_legal_filename(raster_name) < 0) {
- G_message("[%s] is an illegal name", raster_name);
+ G_warning(_("<%s> is an illegal file 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);
+ G_warning(_("Unable to open raster map <%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;
- }
+ /* 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 = G_raster_map_type(raster_name, mapset)) != DCELL_TYPE) {
- G_message("Map is not of DCELL_TYPE");
+ G_warning(_("Map is not DCELL type (process DN to radiance first)"));
return -1;
}
return raster_fd;
@@ -92,10 +94,11 @@
struct GModule *module;
int i, verbose = 1;
- struct Option *input, *output, *hist;
- struct Flag *shadow, *sat5, *filter, *pass2;
+ struct Option *input, *output, *hist, *b56c;
+ struct Flag *shadow, *filter, *sat5, *pass2, *csig;
char *in_name, *out_name;
- double p;
+ struct Categories cats;
+ char title[RECORD_LEN];
Gfile band[5], out;
@@ -108,22 +111,29 @@
_("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_option();
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
output->key = _("output");
output->type = TYPE_STRING;
output->required = YES;
output->gisprompt = _("output,cell,raster");
output->description = _("Output file name");
+ b56c = G_define_option();
+ b56c->key = "b56composite";
+ b56c->type = TYPE_DOUBLE;
+ b56c->required = NO;
+ b56c->description = _("Value for step 6: B56composite");
+ b56c->answer = "225.";
+
hist = G_define_option();
- hist->key = _("histogram");
+ hist->key = "histogram";
hist->type = TYPE_INTEGER;
hist->required = NO;
hist->gisprompt = _("input,integer");
@@ -133,22 +143,30 @@
sat5 = G_define_flag();
sat5->key = '5';
- sat5->description = _("Landsat-5 TM");
+ sat5->description =
+ _("Landsat-5 TM (i.e. thermal band is '.6' not '.61')");
sat5->answer = 0;
filter = G_define_flag();
filter->key = 'f';
- filter->description = _("Use final filter holes");
+ filter->description =
+ _("Apply post-processing filter to remove small holes");
filter->answer = 0;
+ csig = G_define_flag();
+ csig->key = 'x';
+ csig->description = _("Always use cloud signature (step 14)");
+ csig->answer = 0;
+
pass2 = G_define_flag();
pass2->key = '2';
- pass2->description = _("With pass two processing");
+ pass2->description =
+ _("By-pass second processing, and join warm (not ambiguous) and cold clouds");
pass2->answer = 0;
shadow = G_define_flag();
shadow->key = 's';
- shadow->description = _("Add class for cloud shadows");
+ shadow->description = _("Include a category for cloud shadows");
shadow->answer = 0;
if (G_parser(argc, argv))
@@ -166,7 +184,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 filename [%s]!"), band[i].name);
+ G_fatal_error(_("Error in map name <%s>!"), band[i].name);
}
band[i].rast = G_allocate_raster_buf(DCELL_TYPE);
}
@@ -175,17 +193,13 @@
snprintf(out.name, 127, "%s", out_name);
if (G_legal_filename(out_name) < 0)
- G_fatal_error(_("[%s] is an illegal name"), out.name);
+ G_fatal_error(_("<%s> is an illegal file name"), out.name);
/* --------------------------------------- */
+ th_4 = atof(b56c->answer);
+ acca_algorithm(verbose, &out, band, pass2->answer, shadow->answer,
+ csig->answer);
- // 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);
/* --------------------------------------- */
@@ -195,10 +209,22 @@
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);
+ /* 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);
+ 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-29 00:58:17 UTC (rev 43727)
+++ grass-addons/imagery/i.landsat.acca/tools.c 2010-09-29 09:34:05 UTC (rev 43728)
@@ -41,36 +41,23 @@
/* histogram moment */
double moment(int n, int hist[], int k)
{
- int i, j, total;
- double value, hmean, cte;
+ int i, total;
+ double value, mean, cte;
k = 0;
total = 0;
- hmean = 0.;
+ mean = 0.;
for (i = 0; i < hist_n; i++) {
total += hist[i];
- hmean += (double)(i * hist[i]);
+ mean += (double)(i * hist[i]);
}
- hmean /= (double)total; /* histogram mean */
+ mean /= (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;
- }
-
- /* 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);
+ value += (pow((i - mean), n) * (double)hist[i] * cte);
}
return value;
@@ -118,3 +105,229 @@
/* 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 (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(_("Raster map <%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("Unable to open raster map <%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(_("Unable to create raster map <%s>"), tmp.name);
+
+ G_message(_("Filling small holes in cloud ..."));
+
+ /* 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) {
+ if (G_get_c_raster_row(out->fd, arast, row - 1) < 0)
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ out->name, row - 1);
+ }
+ if (G_get_c_raster_row(out->fd, brast, row) < 0) {
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ out->name, row);
+ }
+ if (row != (nrows - 1)) {
+ if (G_get_c_raster_row(out->fd, crast, row + 1) < 0)
+ G_fatal_error(_("Unable to read raster map <%s> row %d"),
+ out->name, 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 {
+ 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(_("Failed writing raster map <%s> row %d"),
+ tmp.name, row);
+
+ G_percent(row, nrows, 2);
+ }
+
+ 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;
+}
+
More information about the grass-commit
mailing list