[GRASS-SVN] r60797 - grass/trunk/imagery/i.albedo
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 12 04:03:48 PDT 2014
Author: ychemin
Date: 2014-06-12 04:03:48 -0700 (Thu, 12 Jun 2014)
New Revision: 60797
Modified:
grass/trunk/imagery/i.albedo/bb_alb_landsat8.c
grass/trunk/imagery/i.albedo/main.c
Log:
better temporary fix
Modified: grass/trunk/imagery/i.albedo/bb_alb_landsat8.c
===================================================================
--- grass/trunk/imagery/i.albedo/bb_alb_landsat8.c 2014-06-11 21:49:13 UTC (rev 60796)
+++ grass/trunk/imagery/i.albedo/bb_alb_landsat8.c 2014-06-12 11:03:48 UTC (rev 60797)
@@ -3,14 +3,24 @@
* Temporary until a publication creates an algorithm
-* chan5 is OLI Band 6 (1.57-1.65)
* chan7 is OLI band 7 (2.11-2.29)
+ *
+ * Temporary better fix than weighted average
+ * ------------------------------------------
+ * r.regression.multi
+ * mapx=LC81270512014115LGN00.toar.1,LC81270512014115LGN00.toar.2,
+ * LC81270512014115LGN00.toar.3,LC81270512014115LGN00.toar.4,
+ * LC81270512014115LGN00.toar.5,LC81270512014115LGN00.toar.6,
+ * LC81270512014115LGN00.toar.7
+ * mapy=MCD43_2014113
+ *
*/
-double bb_alb_landsat8(double bluechan, double greenchan, double redchan,
- double nirchan, double chan5, double chan7)
+double bb_alb_landsat8(double shortbluechan, double bluechan, double greenchan, double redchan,
+ double nirchan, double chan5, double chan7)
{
double result;
- result =
- (0.06 * bluechan + 0.06 * greenchan + 0.03 * redchan +
- 0.03 * nirchan + 0.08 * chan5 + 0.18 * chan7)/0.44;
+ result = 0.058674+shortbluechan*2.153642+bluechan*(-2.242688)+
+ greenchan*(-0.520669)+redchan*0.622670+nirchan*0.129979+
+ chan5*(-0.047970)+chan7*0.152228;
return result;
}
Modified: grass/trunk/imagery/i.albedo/main.c
===================================================================
--- grass/trunk/imagery/i.albedo/main.c 2014-06-11 21:49:13 UTC (rev 60796)
+++ grass/trunk/imagery/i.albedo/main.c 2014-06-12 11:03:48 UTC (rev 60797)
@@ -9,7 +9,7 @@
* COPYRIGHT: (C) 2004-2008 by the GRASS Development Team
*
* This program is free software under the GNU Lesser General Public
- * License. Read the file COPYING that comes with GRASS for details.
+ * License. Read the file COPYING that comes with GRASS for details.
*
*****************************************************************************/
@@ -24,22 +24,23 @@
#define MAXFILES 8
double bb_alb_aster(double greenchan, double nirchan, double swirchan2,
- double swirchan3, double swirchan5, double swirchan6);
+ double swirchan3, double swirchan5, double swirchan6);
double bb_alb_landsat(double bluechan, double greenchan, double redchan,
- double nirchan, double chan5, double chan7);
+ double nirchan, double chan5, double chan7);
-double bb_alb_landsat8(double bluechan, double greenchan, double redchan,
- double nirchan, double chan5, double chan7);
+double bb_alb_landsat8(double shortbluechan, double bluechan,
+ double greenchan, double redchan,
+ double nirchan, double chan5, double chan7);
double bb_alb_noaa(double redchan, double nirchan);
double bb_alb_modis(double redchan, double nirchan, double chan3,
- double chan4, double chan5, double chan6, double chan7);
+ double chan4, double chan5, double chan6, double chan7);
int main(int argc, char *argv[])
{
- struct Cell_head cellhd; /*region+header info */
+ struct Cell_head cellhd; /*region+header info */
int nrows, ncols;
int row, col;
struct GModule *module;
@@ -47,12 +48,12 @@
struct Flag *flag1, *flag2, *flag3;
struct Flag *flag4, *flag5, *flag6;
struct Flag *flag7;
- struct History history; /*metadata */
- struct Colors colors; /*Color rules */
+ struct History history; /*metadata */
+ struct Colors colors; /*Color rules */
/************************************/
- char *name; /*input raster name */
- char *result; /*output raster name */
+ char *name; /*input raster name */
+ char *result; /*output raster name */
/*File Descriptors */
int nfiles;
int infd[MAXFILES];
@@ -65,7 +66,7 @@
void *inrast[MAXFILES];
unsigned char *outrast;
- RASTER_MAP_TYPE in_data_type[MAXFILES]; /* 0=numbers 1=text */
+ RASTER_MAP_TYPE in_data_type[MAXFILES]; /* 0=numbers 1=text */
RASTER_MAP_TYPE out_data_type = DCELL_TYPE;
CELL val1, val2;
@@ -119,7 +120,7 @@
flag4 = G_define_flag();
flag4->key = '8';
- flag4->description = _("Landsat 8 (6 input bands:2,3,4,5,6,7)");
+ flag4->description = _("Landsat 8 (7 input bands:1,2,3,4,5,6,7)");
flag5 = G_define_flag();
flag5->key = 'a';
@@ -129,21 +130,21 @@
flag6->key = 'c';
flag6->label = _("Agressive mode (Landsat)");
flag6->description =
- _("Albedo dry run to calculate some water to beach/sand/desert stretching, "
- "a kind of simple atmospheric correction");
+ _("Albedo dry run to calculate some water to beach/sand/desert stretching, "
+ "a kind of simple atmospheric correction");
flag7 = G_define_flag();
flag7->key = 'd';
flag7->label = _("Soft mode (Modis)");
flag7->description =
- _("Albedo dry run to calculate some water to beach/sand/desert stretching, "
- "a kind of simple atmospheric correction");
+ _("Albedo dry run to calculate some water to beach/sand/desert stretching, "
+ "a kind of simple atmospheric correction");
/* FMEO init nfiles */
nfiles = 1;
if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ exit(EXIT_FAILURE);
names = input->answers;
ptr = input->answers;
@@ -157,22 +158,22 @@
aster = (flag5->answer);
for (; *ptr != NULL; ptr++) {
- if (nfiles >= MAXFILES)
- G_fatal_error(_("Too many input maps. Only %d allowed."), MAXFILES);
- name = *ptr;
-
- /* Allocate input buffer */
- in_data_type[nfiles] = Rast_map_type(name, "");
- infd[nfiles] = Rast_open_old(name, "");
+ if (nfiles >= MAXFILES)
+ G_fatal_error(_("Too many input maps. Only %d allowed."), MAXFILES);
+ name = *ptr;
- Rast_get_cellhd(name, "", &cellhd);
+ /* Allocate input buffer */
+ in_data_type[nfiles] = Rast_map_type(name, "");
+ infd[nfiles] = Rast_open_old(name, "");
- inrast[nfiles] = Rast_allocate_buf(in_data_type[nfiles]);
- nfiles++;
+ Rast_get_cellhd(name, "", &cellhd);
+
+ inrast[nfiles] = Rast_allocate_buf(in_data_type[nfiles]);
+ nfiles++;
}
nfiles--;
if (nfiles <= 1)
- G_fatal_error(_("At least two raster maps are required"));
+ G_fatal_error(_("At least two raster maps are required"));
/* Allocate output buffer, use input map data_type */
nrows = Rast_window_rows();
@@ -186,217 +187,217 @@
/*This is correcting contrast for water/sand */
/*A poor man's atmospheric correction... */
for (i = 0; i < 100; i++)
- histogram[i] = 0;
+ histogram[i] = 0;
if (flag6->answer || flag7->answer) {
- DCELL de;
- DCELL d[MAXFILES];
+ DCELL de;
+ DCELL d[MAXFILES];
- /* Process pixels histogram */
- for (row = 0; row < nrows; row++) {
- G_percent(row, nrows, 2);
- /* read input map */
- for (i = 1; i <= nfiles; i++)
- Rast_get_row(infd[i], inrast[i], row, in_data_type[i]);
+ /* Process pixels histogram */
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ /* read input map */
+ for (i = 1; i <= nfiles; i++)
+ Rast_get_row(infd[i], inrast[i], row, in_data_type[i]);
- /*process the data */
- for (col = 0; col < ncols; col++) {
- for (i = 1; i <= nfiles; i++) {
- switch (in_data_type[i]) {
- case CELL_TYPE:
- d[i] = (double)((CELL *) inrast[i])[col];
- break;
- case FCELL_TYPE:
- d[i] = (double)((FCELL *) inrast[i])[col];
- break;
- case DCELL_TYPE:
- d[i] = (double)((DCELL *) inrast[i])[col];
- break;
- }
- }
- if (modis) {
- de = bb_alb_modis(d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
- }
- else if (avhrr) {
- de = bb_alb_noaa(d[1],d[2]);
- }
- else if (landsat) {
- de = bb_alb_landsat(d[1],d[2],d[3],d[4],d[5],d[6]);
- }
- else if (landsat8) {
- de = bb_alb_landsat8(d[1],d[2],d[3],d[4],d[5],d[6]);
- }
- else if (aster) {
- de = bb_alb_aster(d[1],d[2],d[3],d[4],d[5],d[6]);
- }
- if (Rast_is_d_null_value(&de)) {
- /*Do nothing */
- }
- else {
- int temp;
+ /*process the data */
+ for (col = 0; col < ncols; col++) {
+ for (i = 1; i <= nfiles; i++) {
+ switch (in_data_type[i]) {
+ case CELL_TYPE:
+ d[i] = (double)((CELL *) inrast[i])[col];
+ break;
+ case FCELL_TYPE:
+ d[i] = (double)((FCELL *) inrast[i])[col];
+ break;
+ case DCELL_TYPE:
+ d[i] = (double)((DCELL *) inrast[i])[col];
+ break;
+ }
+ }
+ if (modis) {
+ de = bb_alb_modis(d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
+ }
+ else if (avhrr) {
+ de = bb_alb_noaa(d[1],d[2]);
+ }
+ else if (landsat) {
+ de = bb_alb_landsat(d[1],d[2],d[3],d[4],d[5],d[6]);
+ }
+ else if (landsat8) {
+ de = bb_alb_landsat8(d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
+ }
+ else if (aster) {
+ de = bb_alb_aster(d[1],d[2],d[3],d[4],d[5],d[6]);
+ }
+ if (Rast_is_d_null_value(&de)) {
+ /*Do nothing */
+ }
+ else {
+ int temp;
- temp = (int)(de * 100);
- if (temp > 0) {
- histogram[temp] = histogram[temp] + 1.0;
- }
- }
- }
- }
+ temp = (int)(de * 100);
+ if (temp > 0) {
+ histogram[temp] = histogram[temp] + 1.0;
+ }
+ }
+ }
+ }
- G_message("Calculating histogram of albedo");
+ G_message("Calculating histogram of albedo");
- peak1 = 0;
- peak2 = 0;
- peak3 = 0;
- i_peak1 = 0;
- i_peak2 = 0;
- i_peak3 = 0;
- for (i = 0; i < 100; i++) {
- /*Search for peaks of datasets (1) */
- if (i <= 10) {
- if (histogram[i] > peak1) {
- peak1 = histogram[i];
- i_peak1 = i;
- }
- }
- /*Search for peaks of datasets (2) */
- if (i >= 10 && i <= 30) {
- if (histogram[i] > peak2) {
- peak2 = histogram[i];
- i_peak2 = i;
- }
- }
- /*Search for peaks of datasets (3) */
- if (i >= 30) {
- if (histogram[i] > peak3) {
- peak3 = histogram[i];
- i_peak3 = i;
- }
- }
- }
+ peak1 = 0;
+ peak2 = 0;
+ peak3 = 0;
+ i_peak1 = 0;
+ i_peak2 = 0;
+ i_peak3 = 0;
+ for (i = 0; i < 100; i++) {
+ /*Search for peaks of datasets (1) */
+ if (i <= 10) {
+ if (histogram[i] > peak1) {
+ peak1 = histogram[i];
+ i_peak1 = i;
+ }
+ }
+ /*Search for peaks of datasets (2) */
+ if (i >= 10 && i <= 30) {
+ if (histogram[i] > peak2) {
+ peak2 = histogram[i];
+ i_peak2 = i;
+ }
+ }
+ /*Search for peaks of datasets (3) */
+ if (i >= 30) {
+ if (histogram[i] > peak3) {
+ peak3 = histogram[i];
+ i_peak3 = i;
+ }
+ }
+ }
- bottom1a = 100000;
- bottom1b = 100000;
- bottom2a = 100000;
- bottom2b = 100000;
- bottom3a = 100000;
- bottom3b = 100000;
- i_bottom1a = 100;
- i_bottom1b = 100;
- i_bottom2a = 100;
- i_bottom2b = 100;
- i_bottom3a = 100;
- i_bottom3b = 100;
- /* Water histogram lower bound */
- for (i = 0; i < i_peak1; i++) {
- if (histogram[i] <= bottom1a) {
- bottom1a = histogram[i];
- i_bottom1a = i;
- }
- }
- /* Water histogram higher bound */
- for (i = i_peak2; i > i_peak1; i--) {
- if (histogram[i] <= bottom1b) {
- bottom1b = histogram[i];
- i_bottom1b = i;
- }
- }
- /* Land histogram lower bound */
- for (i = i_peak1; i < i_peak2; i++) {
- if (histogram[i] <= bottom2a) {
- bottom2a = histogram[i];
- i_bottom2a = i;
- }
- }
- /* Land histogram higher bound */
- for (i = i_peak3; i > i_peak2; i--) {
- if (histogram[i] < bottom2b) {
- bottom2b = histogram[i];
- i_bottom2b = i;
- }
- }
- /* Cloud/Snow histogram lower bound */
- for (i = i_peak2; i < i_peak3; i++) {
- if (histogram[i] < bottom3a) {
- bottom3a = histogram[i];
- i_bottom3a = i;
- }
- }
- /* Cloud/Snow histogram higher bound */
- for (i = 100; i > i_peak3; i--) {
- if (histogram[i] < bottom3b) {
- bottom3b = histogram[i];
- i_bottom3b = i;
- }
- }
- if (flag5->answer) {
- G_message("peak1 %d %d", peak1, i_peak1);
- G_message("bottom2b= %d %d", bottom2b, i_bottom2b);
- a = (0.36 - 0.05) / (i_bottom2b / 100.0 - i_peak1 / 100.0);
- b = 0.05 - a * (i_peak1 / 100.0);
- G_message("a= %f\tb= %f", a, b);
- }
- if (flag6->answer) {
- G_message("bottom1a %d %d", bottom1a, i_bottom1a);
- G_message("bottom2b= %d %d", bottom2b, i_bottom2b);
- a = (0.36 - 0.05) / (i_bottom2b / 100.0 - i_bottom1a / 100.0);
- b = 0.05 - a * (i_bottom1a / 100.0);
- G_message("a= %f\tb= %f", a, b);
- }
- } /*END OF FLAG1 */
+ bottom1a = 100000;
+ bottom1b = 100000;
+ bottom2a = 100000;
+ bottom2b = 100000;
+ bottom3a = 100000;
+ bottom3b = 100000;
+ i_bottom1a = 100;
+ i_bottom1b = 100;
+ i_bottom2a = 100;
+ i_bottom2b = 100;
+ i_bottom3a = 100;
+ i_bottom3b = 100;
+ /* Water histogram lower bound */
+ for (i = 0; i < i_peak1; i++) {
+ if (histogram[i] <= bottom1a) {
+ bottom1a = histogram[i];
+ i_bottom1a = i;
+ }
+ }
+ /* Water histogram higher bound */
+ for (i = i_peak2; i > i_peak1; i--) {
+ if (histogram[i] <= bottom1b) {
+ bottom1b = histogram[i];
+ i_bottom1b = i;
+ }
+ }
+ /* Land histogram lower bound */
+ for (i = i_peak1; i < i_peak2; i++) {
+ if (histogram[i] <= bottom2a) {
+ bottom2a = histogram[i];
+ i_bottom2a = i;
+ }
+ }
+ /* Land histogram higher bound */
+ for (i = i_peak3; i > i_peak2; i--) {
+ if (histogram[i] < bottom2b) {
+ bottom2b = histogram[i];
+ i_bottom2b = i;
+ }
+ }
+ /* Cloud/Snow histogram lower bound */
+ for (i = i_peak2; i < i_peak3; i++) {
+ if (histogram[i] < bottom3a) {
+ bottom3a = histogram[i];
+ i_bottom3a = i;
+ }
+ }
+ /* Cloud/Snow histogram higher bound */
+ for (i = 100; i > i_peak3; i--) {
+ if (histogram[i] < bottom3b) {
+ bottom3b = histogram[i];
+ i_bottom3b = i;
+ }
+ }
+ if (flag5->answer) {
+ G_message("peak1 %d %d", peak1, i_peak1);
+ G_message("bottom2b= %d %d", bottom2b, i_bottom2b);
+ a = (0.36 - 0.05) / (i_bottom2b / 100.0 - i_peak1 / 100.0);
+ b = 0.05 - a * (i_peak1 / 100.0);
+ G_message("a= %f\tb= %f", a, b);
+ }
+ if (flag6->answer) {
+ G_message("bottom1a %d %d", bottom1a, i_bottom1a);
+ G_message("bottom2b= %d %d", bottom2b, i_bottom2b);
+ a = (0.36 - 0.05) / (i_bottom2b / 100.0 - i_bottom1a / 100.0);
+ b = 0.05 - a * (i_bottom1a / 100.0);
+ G_message("a= %f\tb= %f", a, b);
+ }
+ } /*END OF FLAG1 */
/* End of processing histogram */
/* Process pixels */
for (row = 0; row < nrows; row++) {
- DCELL de;
- DCELL d[MAXFILES];
+ DCELL de;
+ DCELL d[MAXFILES];
- G_percent(row, nrows, 2);
- /* read input map */
- for (i = 1; i <= nfiles; i++)
- Rast_get_row(infd[i], inrast[i], row, in_data_type[i]);
+ G_percent(row, nrows, 2);
+ /* read input map */
+ for (i = 1; i <= nfiles; i++)
+ Rast_get_row(infd[i], inrast[i], row, in_data_type[i]);
- /*process the data */
- for (col = 0; col < ncols; col++) {
- for (i = 1; i <= nfiles; i++) {
- switch (in_data_type[i]) {
- case CELL_TYPE:
- d[i] = (double)((CELL *) inrast[i])[col];
- break;
- case FCELL_TYPE:
- d[i] = (double)((FCELL *) inrast[i])[col];
- break;
- case DCELL_TYPE:
- d[i] = (double)((DCELL *) inrast[i])[col];
- break;
- }
- }
- if (modis) {
- de = bb_alb_modis(d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
- }
- else if (avhrr) {
- de = bb_alb_noaa(d[1],d[2]);
- }
- else if (landsat) {
- de = bb_alb_landsat(d[1],d[2],d[3],d[4],d[5],d[6]);
- }
- else if (landsat8) {
- de = bb_alb_landsat8(d[1],d[2],d[3],d[4],d[5],d[6]);
- }
- else if (aster) {
- de = bb_alb_aster(d[1],d[2],d[3],d[4],d[5],d[6]);
- }
- if (flag6->answer || flag7->answer) {
- /* Post-Process Albedo */
- de = a * de + b;
- }
- ((DCELL *) outrast)[col] = de;
- }
- Rast_put_row(outfd, outrast, out_data_type);
+ /*process the data */
+ for (col = 0; col < ncols; col++) {
+ for (i = 1; i <= nfiles; i++) {
+ switch (in_data_type[i]) {
+ case CELL_TYPE:
+ d[i] = (double)((CELL *) inrast[i])[col];
+ break;
+ case FCELL_TYPE:
+ d[i] = (double)((FCELL *) inrast[i])[col];
+ break;
+ case DCELL_TYPE:
+ d[i] = (double)((DCELL *) inrast[i])[col];
+ break;
+ }
+ }
+ if (modis) {
+ de = bb_alb_modis(d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
+ }
+ else if (avhrr) {
+ de = bb_alb_noaa(d[1],d[2]);
+ }
+ else if (landsat) {
+ de = bb_alb_landsat(d[1],d[2],d[3],d[4],d[5],d[6]);
+ }
+ else if (landsat8) {
+ de = bb_alb_landsat8(d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
+ }
+ else if (aster) {
+ de = bb_alb_aster(d[1],d[2],d[3],d[4],d[5],d[6]);
+ }
+ if (flag6->answer || flag7->answer) {
+ /* Post-Process Albedo */
+ de = a * de + b;
+ }
+ ((DCELL *) outrast)[col] = de;
}
+ Rast_put_row(outfd, outrast, out_data_type);
+ }
for (i = 1; i <= nfiles; i++) {
- G_free(inrast[i]);
- Rast_close(infd[i]);
+ G_free(inrast[i]);
+ Rast_close(infd[i]);
}
G_free(outrast);
Rast_close(outfd);
More information about the grass-commit
mailing list