[GRASS-SVN] r33564 - in grass/trunk/imagery: . i.albedo i.qc.modis
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Sep 27 11:11:06 EDT 2008
Author: ychemin
Date: 2008-09-27 11:11:05 -0400 (Sat, 27 Sep 2008)
New Revision: 33564
Added:
grass/trunk/imagery/i.albedo/
grass/trunk/imagery/i.albedo/Makefile
grass/trunk/imagery/i.albedo/bb_alb_aster.c
grass/trunk/imagery/i.albedo/bb_alb_landsat.c
grass/trunk/imagery/i.albedo/bb_alb_modis.c
grass/trunk/imagery/i.albedo/bb_alb_noaa.c
grass/trunk/imagery/i.albedo/description.html
grass/trunk/imagery/i.albedo/functions.h
grass/trunk/imagery/i.albedo/main.c
grass/trunk/imagery/i.qc.modis/
grass/trunk/imagery/i.qc.modis/Makefile
grass/trunk/imagery/i.qc.modis/description.html
grass/trunk/imagery/i.qc.modis/main.c
grass/trunk/imagery/i.qc.modis/qc250a.c
grass/trunk/imagery/i.qc.modis/qc250b.c
grass/trunk/imagery/i.qc.modis/qc250c.c
grass/trunk/imagery/i.qc.modis/qc250d.c
grass/trunk/imagery/i.qc.modis/qc250e.c
grass/trunk/imagery/i.qc.modis/qc250f.c
grass/trunk/imagery/i.qc.modis/qc500a.c
grass/trunk/imagery/i.qc.modis/qc500c.c
grass/trunk/imagery/i.qc.modis/qc500d.c
grass/trunk/imagery/i.qc.modis/qc500e.c
Log:
Added i.albedo and i.qc.modis
Added: grass/trunk/imagery/i.albedo/Makefile
===================================================================
--- grass/trunk/imagery/i.albedo/Makefile (rev 0)
+++ grass/trunk/imagery/i.albedo/Makefile 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.albedo
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+ EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd
Added: grass/trunk/imagery/i.albedo/bb_alb_aster.c
===================================================================
--- grass/trunk/imagery/i.albedo/bb_alb_aster.c (rev 0)
+++ grass/trunk/imagery/i.albedo/bb_alb_aster.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,16 @@
+/* Broadband albedo Aster
+ * After Liang, S.L., 2001.
+ * Narrowband to broadband conversion of land surface albedo 1 Algorithms.
+ * Remote Sensing of Environment. 2001, 76, 213-238.
+ * Input: Ref1, ref3, Ref5, Ref6, Ref8, Ref9
+ */
+double bb_alb_aster(double greenchan, double nirchan, double swirchan2,
+ double swirchan3, double swirchan5, double swirchan6)
+{
+ double result;
+
+ result =
+ (0.484 * greenchan + 0.335 * nirchan - 0.324 * swirchan2 +
+ 0.551 * swirchan3 + 0.305 * swirchan5 - 0.367 * swirchan6 - 0.0015);
+ return result;
+}
Property changes on: grass/trunk/imagery/i.albedo/bb_alb_aster.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/imagery/i.albedo/bb_alb_landsat.c
===================================================================
--- grass/trunk/imagery/i.albedo/bb_alb_landsat.c (rev 0)
+++ grass/trunk/imagery/i.albedo/bb_alb_landsat.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,13 @@
+/* Broadband albedo Landsat 5TM and 7ETM+
+ * (maybe others too but not sure)
+ */
+double bb_alb_landsat(double bluechan, double greenchan, double redchan,
+ double nirchan, double chan5, double chan7)
+{
+ double result;
+
+ result =
+ (0.293 * bluechan + 0.274 * greenchan + 0.233 * redchan +
+ 0.156 * nirchan + 0.033 * chan5 + 0.011 * chan7);
+ return result;
+}
Property changes on: grass/trunk/imagery/i.albedo/bb_alb_landsat.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/imagery/i.albedo/bb_alb_modis.c
===================================================================
--- grass/trunk/imagery/i.albedo/bb_alb_modis.c (rev 0)
+++ grass/trunk/imagery/i.albedo/bb_alb_modis.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,14 @@
+/*
+ * Broadband albedo MODIS
+ */
+double bb_alb_modis(double redchan, double nirchan, double chan3,
+ double chan4, double chan5, double chan6, double chan7)
+{
+ double result;
+
+ result =
+ (0.22831 * redchan + 0.15982 * nirchan +
+ 0.09132 * (chan3 + chan4 + chan5) + 0.10959 * chan6 +
+ 0.22831 * chan7);
+ return result;
+}
Property changes on: grass/trunk/imagery/i.albedo/bb_alb_modis.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/imagery/i.albedo/bb_alb_noaa.c
===================================================================
--- grass/trunk/imagery/i.albedo/bb_alb_noaa.c (rev 0)
+++ grass/trunk/imagery/i.albedo/bb_alb_noaa.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,11 @@
+/*
+ * Broadband albedo NOAA AVHRR 14
+ * (maybe others too but not sure)
+ */
+double bb_alb_noaa(double redchan, double nirchan)
+{
+ double result;
+
+ result = (0.035 + 0.545 * nirchan - 0.32 * redchan);
+ return result;
+}
Property changes on: grass/trunk/imagery/i.albedo/bb_alb_noaa.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass/trunk/imagery/i.albedo/description.html
===================================================================
--- grass/trunk/imagery/i.albedo/description.html (rev 0)
+++ grass/trunk/imagery/i.albedo/description.html 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,24 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.albedo</EM> calculates the Albedo, that is the Shortwave surface reflectance in the range of 0.3-3 micro-meters.
+It takes input of individual bands of surface reflectance from Modis, AVHRR, Landsat or Aster and calculates the Albedo for those.
+This is an precursor to r.sun and any Energy-Balance processing.
+<H2>NOTES</H2>
+It assumes MODIS product surface reflectance in [0;10000]
+<H2>TODO</H2>
+Maybe change input requirement of MODIS to [0.0-1.0]?
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="r.sun.html">r.sun</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, International Rice Research Institute, The Philippines<BR>
+
+
+<p>
+<i>Last changed: $Date: 2006/10/08 11:41:43 $</i>
Added: grass/trunk/imagery/i.albedo/functions.h
===================================================================
--- grass/trunk/imagery/i.albedo/functions.h (rev 0)
+++ grass/trunk/imagery/i.albedo/functions.h 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,13 @@
+/* These are the headers for the RS functions */
+/* 2004 */
+
+/* BB_albedo functions */
+double bb_alb_aster(double greenchan, double redchan, double nirchan,
+ double swirchan1, double swirchan2, double swirchan3,
+ double swirchan4, double swirchan5, double swirchan6);
+double bb_alb_landsat(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);
Added: grass/trunk/imagery/i.albedo/main.c
===================================================================
--- grass/trunk/imagery/i.albedo/main.c (rev 0)
+++ grass/trunk/imagery/i.albedo/main.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,461 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.albedo
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Calculate Broadband Albedo (0.3-3 Micrometers)
+ * from Surface Reflectance (Modis, AVHRR, Landsat, Aster).
+ *
+ * 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.
+ *
+ *****************************************************************************/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#define MAXFILES 8
+
+double bb_alb_aster(double greenchan, double nirchan, double swirchan2,
+ double swirchan3, double swirchan5, double swirchan6);
+double bb_alb_landsat(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);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /*region+header info */
+
+ char *mapset; /*mapset name */
+
+ int nrows, ncols;
+
+ int row, col;
+
+ struct GModule *module;
+
+ struct Option *input, *output;
+
+ struct Flag *flag1, *flag2, *flag3;
+
+ struct Flag *flag4, *flag5, *flag6;
+
+ struct History history; /*metadata */
+
+ struct Colors colors; /*Color rules */
+
+ /************************************/
+ /* FMEO Declarations**************** */
+ char *name; /*input raster name */
+
+ char *result; /*output raster name */
+
+ /*File Descriptors */
+ int nfiles;
+
+ int infd[MAXFILES];
+
+ int outfd;
+
+ char **names;
+
+ char **ptr;
+
+ int ok;
+
+ int i = 0, j = 0;
+
+ int modis = 0, aster = 0, avhrr = 0, landsat = 0;
+
+ void *inrast[MAXFILES];
+
+ unsigned char *outrast;
+
+ int data_format; /* 0=double 1=float 2=32bit signed int 5=8bit unsigned int (ie text) */
+
+ RASTER_MAP_TYPE in_data_type[MAXFILES]; /* 0=numbers 1=text */
+
+ RASTER_MAP_TYPE out_data_type = DCELL_TYPE;
+
+ /************************************/
+
+ /************************************/
+ int histogram[100];
+
+ /* Albedo correction coefficients*** */
+ double a, b;
+
+ /************************************/
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("Albedo,surface reflectance,r.sun");
+ module->description =
+ _("Broad Band Albedo from Surface Reflectance.\n NOAA AVHRR(n), Modis(m), Landsat(l), Aster(a)\n");
+
+ /* Define the different options */
+
+ input = G_define_standard_option(G_OPT_R_INPUT);
+ input->multiple = YES;
+ input->description = _("Names of surface reflectance layers");
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->description = _("Name of the BB_Albedo layer");
+
+ /* Define the different flags */
+
+ flag1 = G_define_flag();
+ flag1->key = _('m');
+ flag1->description = _("Modis (7 input bands:1,2,3,4,5,6,7)");
+
+ flag2 = G_define_flag();
+ flag2->key = _('n');
+ flag2->description = _("NOAA AVHRR (2 input bands:1,2)");
+
+ flag3 = G_define_flag();
+ flag3->key = _('l');
+ flag3->description = _("Landsat (6 input bands:1,2,3,4,5,7)");
+
+ flag4 = G_define_flag();
+ flag4->key = _('a');
+ flag4->description = _("Aster (6 input bands:1,3,5,6,8,9)");
+
+ flag5 = G_define_flag();
+ flag5->key = _('c');
+ flag5->description =
+ _("Albedo dry run to calculate some water to beach/sand/desert stretching, a kind of simple atmospheric correction. Agressive mode (Landsat).");
+
+ flag6 = G_define_flag();
+ flag6->key = _('d');
+ flag6->description =
+ _("Albedo dry run to calculate some water to beach/sand/desert stretching, a kind of simple atmospheric correction. Soft mode (Modis).");
+
+ /* FMEO init nfiles */
+ nfiles = 1;
+
+ /********************/
+ if (G_parser(argc, argv))
+ exit(-1);
+
+ ok = 1;
+ names = input->answers;
+ ptr = input->answers;
+
+ result = output->answer;
+
+ modis = (flag1->answer);
+ avhrr = (flag2->answer);
+ landsat = (flag3->answer);
+ aster = (flag4->answer);
+
+ /***************************************************/
+ if (G_legal_filename(result) < 0) {
+ G_fatal_error(_("[%s] is an illegal name"), result);
+ }
+
+ /***************************************************/
+
+ /***************************************************/
+ for (; *ptr != NULL; ptr++) {
+ if (nfiles >= MAXFILES)
+ G_fatal_error(_("%s - too many ETa files. Only %d allowed"),
+ G_program_name(), MAXFILES);
+ name = *ptr;
+ /* find map in mapset */
+ mapset = G_find_cell2(name, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), name);
+ ok = 0;
+ }
+ if (!ok) {
+ continue;
+ }
+ infd[nfiles] = G_open_cell_old(name, mapset);
+ if (infd[nfiles] < 0) {
+ ok = 0;
+ continue;
+ }
+ /* Allocate input buffer */
+ in_data_type[nfiles] = G_raster_map_type(name, mapset);
+ if ((infd[nfiles] = G_open_cell_old(name, mapset)) < 0) {
+ G_fatal_error(_("Cannot open cell file [%s]"), name);
+ }
+ if ((G_get_cellhd(name, mapset, &cellhd)) < 0) {
+ G_fatal_error(_("Cannot read file header of [%s]"), name);
+ }
+ inrast[nfiles] = G_allocate_raster_buf(in_data_type[nfiles]);
+ nfiles++;
+ }
+ nfiles--;
+ if (nfiles <= 1) {
+ G_fatal_error(_("The min specified input map is two (that is NOAA AVHRR)"));
+ }
+
+ /***************************************************/
+ /* Allocate output buffer, use input map data_type */
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ outrast = G_allocate_raster_buf(out_data_type);
+
+ /* Create New raster files */
+ if ((outfd = G_open_raster_new(result, 1)) < 0) {
+ G_fatal_error(_("Could not open <%s>"), result);
+ }
+ /*START ALBEDO HISTOGRAM STRETCH */
+ /*This is correcting contrast for water/sand */
+ /*A poor man's atmospheric correction... */
+ for (i = 0; i < 100; i++) {
+ histogram[i] = 0;
+ }
+ if (flag5->answer || flag6->answer) {
+ 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++) {
+ if ((G_get_raster_row
+ (infd[i], inrast[i], row, in_data_type[i])) < 0) {
+ G_fatal_error(_("Could not read from <%s>"), name);
+ }
+ }
+ /*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 (aster) {
+ de = bb_alb_aster(d[1], d[2], d[3], d[4], d[5], d[6]);
+ }
+ if (G_is_d_null_value(&de)) {
+ /*Do nothing */
+ }
+ else {
+ int temp;
+
+ temp = (int)(de * 100);
+ if (temp > 0) {
+ histogram[temp] = histogram[temp] + 1.0;
+ }
+ }
+ }
+ }
+ G_message("Histogram of Albedo\n");
+ int peak1, peak2, peak3;
+
+ int i_peak1, i_peak2, i_peak3;
+
+ 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;
+ }
+ }
+ }
+ int bottom1a, bottom1b;
+
+ int bottom2a, bottom2b;
+
+ int bottom3a, bottom3b;
+
+ int i_bottom1a, i_bottom1b;
+
+ int i_bottom2a, i_bottom2b;
+
+ int i_bottom3a, i_bottom3b;
+
+ 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\n", peak1, i_peak1);
+ G_message("bottom2b= %d %d\n", 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\n", a, b);
+ }
+ if (flag6->answer) {
+ G_message("bottom1a %d %d\n", bottom1a, i_bottom1a);
+ G_message("bottom2b= %d %d\n", 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\n", a, b);
+ }
+ } /*END OF FLAG1 */
+ /* End of processing histogram */
+
+ /*******************/
+ /* Process pixels */
+ for (row = 0; row < nrows; row++) {
+ DCELL de;
+
+ DCELL d[MAXFILES];
+
+ G_percent(row, nrows, 2);
+ /* read input map */
+ for (i = 1; i <= nfiles; i++) {
+ if ((G_get_raster_row(infd[i], inrast[i], row, in_data_type[i])) <
+ 0) {
+ G_fatal_error(_("Could not read from <%s>"), name);
+ }
+ }
+ /*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 (aster) {
+ de = bb_alb_aster(d[1], d[2], d[3], d[4], d[5], d[6]);
+ }
+ if (flag5->answer || flag6->answer) {
+ /* Post-Process Albedo */
+ de = a * de + b;
+ }
+ ((DCELL *) outrast)[col] = de;
+ }
+ if (G_put_raster_row(outfd, outrast, out_data_type) < 0)
+ G_fatal_error(_("Cannot write to <%s>"), result);
+ }
+ for (i = 1; i <= nfiles; i++) {
+ G_free(inrast[i]);
+ G_close_cell(infd[i]);
+ }
+ G_free(outrast);
+ G_close_cell(outfd);
+
+ /* Color table from 0.0 to 1.0 */
+ G_init_colors(&colors);
+ G_add_color_rule(0.0, 0, 0, 0, 1.0, 255, 255, 255, &colors);
+ /* Metadata */
+ G_short_history(result, "raster", &history);
+ G_command_history(&history);
+ G_write_history(result, &history);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass/trunk/imagery/i.qc.modis/Makefile
===================================================================
--- grass/trunk/imagery/i.qc.modis/Makefile (rev 0)
+++ grass/trunk/imagery/i.qc.modis/Makefile 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.qc.modis
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+ EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd
Added: grass/trunk/imagery/i.qc.modis/description.html
===================================================================
--- grass/trunk/imagery/i.qc.modis/description.html (rev 0)
+++ grass/trunk/imagery/i.qc.modis/description.html 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,63 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.qc.modis</EM> Extracts Requested Quality Assessment flags from Modis 09 A and Q products.
+
+<EM>MODLAND QA Bits 250m/500m bits[0-1]</EM>
+[00]= class 1: Corrected product produced at ideal quality -- all bands
+[01]= class 2: Corrected product produced at less than idel quality -- some or all bands
+[10]= class 3: Corrected product NOT produced due to cloud effect -- all bands
+[11]= class 4: Corrected product NOT produced due to other reasons -- some or all bands maybe be fill value (Note that a value of [11] overrides a value of [01])
+
+<EM> Cloud State 250m Unsigned Int bits[2-3] </EM>
+[00]= class 1: Clear -- No clouds
+[01]= class 2: Cloudy
+[10]= class 3: Mixed
+[11]= class 4: Not Set ; Assumed Clear
+
+<EM> Band-wise Data Quality 250m Unsigned Int bits[4-7][8-11]</EM>
+<EM> Band-wise Data Quality 500m long Int bits[2-5][6-9][10-13][14-17][18-21][22-25][26-29]</EM>
+[0000]= class 1: highest quality
+[0111]= class 2: noisy detector
+[1000]= class 3: dead detector; data interpolated in L1B
+[1001]= class 4: solar zenith ge 86 degrees
+[1010]= class 5: solar zenith ge 85 and lt 86 degrees
+[1011]= class 6: missing input
+[1100] class 7: internal constant used in place of climatological data for at least one atmospheric constant
+[1101]= class 8: correction out of bounds, pixel constrained to extreme allowable value
+[1110]= class 9: L1B data faulty
+[1111]= class 10: not processed due to deep ocean or cloud
+Class 11: Combination of bits unused
+
+<EM>Atmospheric correction 250m/500m bit[12]/[30]</EM>
+[00]= class 1: Not Corrected product
+[01]= class 2: Corrected product
+
+<EM>Adjacency correction 250m/500m bit[13][31]</EM>
+[00]= class 1: Not Corrected product
+[01]= class 2: Corrected product
+
+<EM>Different orbit from 500m product, 250m Unsigned Int bit[14]</EM>
+[00]= class 1: same orbit as 500m
+[01]= class 2: different orbit from 500m
+
+
+
+<H2>NOTES</H2>
+
+
+<H2>TODO</H2>
+Testing is required.
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.vi.html">i.vi</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+Yann Chemin, International Rice Research Institute, The Philippines.<BR>
+
+
+<p>
+<i>Last changed: $Date: 2008/07/28 15:33:42 $</i>
Added: grass/trunk/imagery/i.qc.modis/main.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/main.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/main.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,326 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.qc.modis
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Converts Quality Control indicators into human readable classes
+ * for Modis surface reflectance products 250m/500m
+ * (MOD09Q/MOD09A)
+ *
+ * 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>
+
+ /* 250m Products (MOD09Q) */
+int qc250a(unsigned int pixel);
+
+int qc250b(unsigned int pixel);
+
+int qc250c(unsigned int pixel, int bandno);
+
+int qc250d(unsigned int pixel);
+
+int qc250e(unsigned int pixel);
+
+int qc250f(unsigned int pixel);
+
+
+ /* 500m Products (MOD09A) */
+int qc500a(long int pixel);
+
+int qc500c(long int pixel, int bandno);
+
+int qc500d(long int pixel);
+
+int qc500e(long int pixel);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /*region+header info */
+
+ char *mapset; /*mapset name */
+
+ int nrows, ncols;
+
+ int row, col;
+
+ char *qcflag; /*Switch for particular index */
+
+ struct GModule *module;
+
+ struct Option *input1, *input2, *input_band, *output;
+
+ struct Flag *flag1, *flag2;
+
+ struct History history; /*metadata */
+
+ struct Colors colors; /*Color rules */
+
+
+
+ /************************************/
+ /* FMEO Declarations**************** */
+ char *name; /*input raster name */
+
+ char *result; /*output raster name */
+
+
+ /*File Descriptors */
+ int infd;
+
+ int outfd;
+
+ char *qcchan;
+
+ int bandno;
+
+ int i = 0, j = 0;
+
+ void *inrast;
+
+ CELL * outrast;
+ RASTER_MAP_TYPE data_type_output = CELL_TYPE;
+ RASTER_MAP_TYPE data_type_qcchan;
+
+
+ /************************************/
+ G_gisinit(argv[0]);
+ module = G_define_module();
+ module->keywords = _("QC, Quality Control, surface reflectance, Modis");
+ module->description =
+ _("Extract quality control parameters from Modis QC layers");
+
+ /* Define the different options */
+ input1 = G_define_option();
+ input1->key = _("qcname");
+ input1->type = TYPE_STRING;
+ input1->required = YES;
+ input1->gisprompt = _("Name of QC type to extract");
+ input1->description =
+ _("Name of QC: modland_qa_bits, cloud, data_quality, atcorr, adjcorr, diff_orbit_from_500m");
+ input1->answer = _("modland_qa_bits");
+ input2 = G_define_standard_option(G_OPT_R_INPUT);
+ input2->description =
+ _("Name of the surface reflectance QC layer [bit array]");
+ input_band = G_define_option();
+ input_band->key = "band";
+ input_band->type = TYPE_INTEGER;
+ input_band->required = NO;
+ input_band->gisprompt = "old,value";
+ input_band->description =
+ _("Band number of Modis product 250m=[1,2],500m=[1-7]");
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->key = _("output");
+ output->description =
+ _("Name of the output QC type classification layer");
+ output->answer = _("qc");
+ flag1 = G_define_flag();
+ flag1->key = 'A';
+ flag1->description =
+ _("QC for MOD09A product @ 500m instead of default MOD09Q at 250m");
+
+
+ /********************/
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+ qcflag = input1->answer;
+ qcchan = input2->answer;
+ if (input_band->answer) {
+ bandno = atoi(input_band->answer);
+ }
+ result = output->answer;
+
+
+ /***************************************************/
+ if ((!strcoll(qcflag, "cloud") && flag1->answer) ||
+ (!strcoll(qcflag, "diff_orbit_from_500m") && flag1->answer)) {
+ G_fatal_error("Those flags cannot work with MOD09A @ 500m products");
+ }
+ if (!strcoll(qcflag, "data_quality")) {
+ if (bandno < 1 || bandno > 7)
+ G_fatal_error("band number out of allowed range [1-7]");
+ if (!flag1->answer && bandno > 2)
+ G_fatal_error("250m band number is out of allowed range [1,2]");
+ }
+
+
+ /***************************************************/
+ mapset = G_find_cell2(qcchan, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), qcchan);
+ }
+ data_type_qcchan = G_raster_map_type(qcchan, mapset);
+ if ((infd = G_open_cell_old(qcchan, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), qcchan);
+ if (G_get_cellhd(qcchan, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), qcchan);
+ inrast = G_allocate_raster_buf(data_type_qcchan);
+
+
+ /***************************************************/
+ G_debug(3, "number of rows %d", cellhd.rows);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ outrast = G_allocate_raster_buf(data_type_output);
+
+ /* Create New raster files */
+ if ((outfd = G_open_raster_new(result, data_type_output)) < 0)
+ G_fatal_error(_("Could not open <%s>"), result);
+
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ CELL c;
+ unsigned int qc250chan;
+
+ CELL qc500chan;
+ G_percent(row, nrows, 2);
+ if (G_get_raster_row(infd, inrast, row, data_type_qcchan) < 0)
+ G_fatal_error(_("Could not read from <%s>"), qcchan);
+
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ switch (data_type_qcchan) {
+ case CELL_TYPE:
+ c = (int)((CELL *) inrast)[col];
+ break;
+ case FCELL_TYPE:
+ c = (int)((FCELL *) inrast)[col];
+ break;
+ case DCELL_TYPE:
+ c = (int)((DCELL *) inrast)[col];
+ break;
+ }
+ if (flag1->answer) {
+
+ /* This is a MOD09A at 500m product, QC layer is 32-bit */
+ qc500chan = c;
+ }
+ else {
+
+ /* This is a MOD09A at 250m product, QC layer is 16-bit */
+ qc250chan = (unsigned int)((CELL *) inrast)[col];
+ } if (G_is_c_null_value(&c)) {
+ G_set_c_null_value(&outrast[col], 1);
+ }
+ else {
+
+ /*calculate modland QA bits extraction */
+ if (!strcoll(qcflag, "modland_qa_bits")) {
+ if (flag1->answer) {
+
+ /* 500m product */
+ c = qc500a(qc500chan);
+ }
+ else {
+
+ /* 250m product */
+ c = qc250a(qc250chan);
+ }
+ outrast[col] = c;
+ }
+ else
+
+ /*calculate cloud state */
+ if (!strcoll(qcflag, "cloud")) {
+
+ /* ONLY 250m product! */
+ c = qc250b(qc250chan);
+ outrast[col] = c;
+ }
+ else
+
+ /*calculate modland QA bits extraction */
+ if (!strcoll(qcflag, "data_quality")) {
+ if (flag1->answer) {
+
+ /* 500m product */
+ c = qc500c(qc500chan, bandno);
+ }
+ else {
+
+ /* 250m product */
+ c = qc250c(qc250chan, bandno);
+ }
+ outrast[col] = c;
+ }
+ else
+
+ /*calculate atmospheric correction flag */
+ if (!strcoll(qcflag, "atcorr")) {
+ if (flag1->answer) {
+
+ /* 500m product */
+ c = qc500d(qc500chan);
+ }
+ else {
+
+ /* 250m product */
+ c = qc250d(qc250chan);
+ }
+ outrast[col] = c;
+ }
+ else
+
+ /*calculate adjacency correction flag */
+ if (!strcoll(qcflag, "adjcorr")) {
+ if (flag1->answer) {
+
+ /* 500m product */
+ c = qc500e(qc500chan);
+ }
+ else {
+
+ /* 250m product */
+ c = qc250e(qc250chan);
+ }
+ outrast[col] = c;
+ }
+ else
+
+ /*calculate different orbit from 500m flag */
+ if (!strcoll(qcflag, "diff_orbit_from_500m")) {
+
+ /* ONLY 250m product! */
+ c = qc250f(qc500chan);
+ outrast[col] = c;
+ }
+ else {
+
+ /* Signal user that the flag name is badly written */
+ /* therefore not understood by the application */
+ G_fatal_error
+ ("Unknown flag name, please check spelling");
+ }
+ }
+ }
+ if (G_put_raster_row(outfd, outrast, data_type_output) < 0)
+ G_fatal_error(_("Cannot write to output raster file"));
+ }
+ G_free(inrast);
+ G_close_cell(infd);
+ G_free(outrast);
+ G_close_cell(outfd);
+
+ /* Color from -1.0 to +1.0 in grey */
+ G_init_colors(&colors);
+ G_add_color_rule(0, 0, 0, 0, 10, 255, 255, 255, &colors);
+ G_short_history(result, "raster", &history);
+ G_command_history(&history);
+ G_write_history(result, &history);
+ exit(EXIT_SUCCESS);
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc250a.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc250a.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc250a.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,39 @@
+/* MODLAND QA Bits 250m Unsigned Int bits[0-1]
+ * 00 -> class 1: Corrected product produced at ideal quality -- all bands
+ * 01 -> class 2: Corrected product produced at less than idel quality -- some or all bands
+ * 10 -> class 3: Corrected product NOT produced due to cloud effect -- all bands
+ * 11 -> class 4: Corrected product NOT produced due to other reasons -- some or all bands mayb be fill value (Note that a value of [11] overrides a value of [01])
+ */
+int qc250a(unsigned int pixel)
+{
+ unsigned int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x02) {
+
+ /* Non-Corrected product */
+ if (qctemp & 0x01) {
+ class = 4; /*other reasons */
+ }
+ else {
+ class = 3; /*because of clouds */
+ }
+ }
+ else {
+
+ /* Corrected product */
+ if (qctemp & 0x01) {
+ class = 2; /*less than ideal quality */
+ }
+ else {
+ class = 1; /*ideal quality */
+ }
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc250b.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc250b.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc250b.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,36 @@
+/* Cloud State 250m Unsigned Int bits[2-3]
+ * 00 -> class 1: Clear -- No clouds
+ * 01 -> class 2: Cloudy
+ * 10 -> class 3: Mixed
+ * 11 -> class 4: Not Set ; Assumed Clear
+ */
+int qc250b(unsigned int pixel)
+{
+ unsigned int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 2; /*bits [2-3] become [0-1] */
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x02) { /* 1 ? */
+ if (qctemp & 0x01) {
+ class = 4; /*Not Set ; Assumed Clear */
+ }
+ else {
+ class = 3; /*Mixed clouds */
+ }
+ }
+ else { /* 0 ? */
+ if (qctemp & 0x01) {
+ class = 2; /*Cloudy */
+ }
+ else {
+ class = 1; /*Clear -- No Clouds */
+ }
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc250c.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc250c.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc250c.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,103 @@
+/* Band-wise Data Quality 250m Unsigned Int bits[0-1]
+ * 0000 -> class 1: highest quality
+ * 0111 -> class 2: noisy detector
+ * 1000 -> class 3: dead detector; data interpolated in L1B
+ * 1001 -> class 4: solar zenith >= 86 degrees
+ * 1010 -> class 5: solar zenith >= 85 and < 86 degrees
+ * 1011 -> class 6: missing input
+ * 1100 -> class 7: internal constant used in place of climatological data for at least one atmospheric constant
+ * 1101 -> class 8: correction out of bounds, pixel constrained to extreme allowable value
+ * 1110 -> class 9: L1B data faulty
+ * 1111 -> class 10: not processed due to deep ocean or cloud
+ * Class 11: Combination of bits unused
+ */
+int qc250c(unsigned int pixel, int bandno)
+{
+ unsigned int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 4 * bandno; /* bitshift [4-7] or [8-11] to [0-3] */
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x08) { /* 1??? */
+ if (qctemp & 0x04) { /* 11?? */
+ if (qctemp & 0x02) { /* 111? */
+ if (qctemp & 0x01) { /* 1111 */
+ class = 10; /*Not proc.ocean/cloud */
+ }
+ else { /* 1110 */
+ class = 9; /*L1B faulty data */
+ }
+ }
+ else { /* 110? */
+ if (qctemp & 0x01) { /* 1101 */
+ class = 8; /*corr. out of bounds */
+ }
+ else { /* 1100 */
+ class = 7; /*internal constant used */
+ }
+ }
+ }
+ else {
+ if (qctemp & 0x02) { /* 101? */
+ if (qctemp & 0x01) { /* 1011 */
+ class = 6; /*missing input */
+ }
+ else { /* 1010 */
+ class = 5; /*solarzen>=85&&<86 */
+ }
+ }
+ else { /* 100? */
+ if (qctemp & 0x01) { /* 1001 */
+ class = 4; /*solar zenith angle>=86 */
+ }
+ else { /* 1000 */
+ class = 3; /*Dead detector */
+ }
+ }
+ }
+ }
+ else { /* 0??? */
+ if (qctemp & 0x04) { /* 01?? */
+ if (qctemp & 0x02) { /* 011? */
+ if (qctemp & 0x01) { /* 0111 */
+ class = 2; /*Noisy detector */
+ }
+ else { /* 0110 */
+ class = 11; /*Unused */
+ }
+ }
+ else { /* 010? */
+ if (qctemp & 0x01) { /* 0101 */
+ class = 11; /*Unused */
+ }
+ else { /* 0100 */
+ class = 11; /*Unused */
+ }
+ }
+ }
+ else { /* 00?? */
+ if (qctemp & 0x02) { /* 001? */
+ if (qctemp & 0x01) { /* 0011 */
+ class = 11; /*Unused */
+ }
+ else { /* 0010 */
+ class = 11; /*Unused */
+ }
+ }
+ else { /* 000? */
+ if (qctemp & 0x01) { /* 0001 */
+ class = 11; /*Unused */
+ }
+ else { /* 0000 */
+ class = 1; /*Highest quality */
+ }
+ }
+ }
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc250d.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc250d.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc250d.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,24 @@
+/* Atmospheric correction 250m Unsigned Int bit[12]
+ * 00 -> class 1: Not Corrected product
+ * 01 -> class 2: Corrected product
+ */
+int qc250d(unsigned int pixel)
+{
+ unsigned int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 12; /* bit no 12 becomes 0 */
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x01) {
+ class = 2; /*Corrected */
+ }
+ else {
+ class = 1; /*Not corrected */
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc250e.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc250e.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc250e.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,24 @@
+/* Adjacency correction 250m Unsigned Int bit[13]
+ * 00 -> class 1: Not Corrected product
+ * 01 -> class 2: Corrected product
+ */
+int qc250e(unsigned int pixel)
+{
+ unsigned int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 13; /* bit no 13 becomes 0 */
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x01) {
+ class = 2; /*Corrected */
+ }
+ else {
+ class = 1; /*Not corrected */
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc250f.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc250f.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc250f.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,24 @@
+/* Different orbit from 500m product, 250m Unsigned Int bit[14]
+ * 00 -> class 1: same orbit as 500m
+ * 01 -> class 2: different orbit from 500m
+ */
+int qc250f(unsigned int pixel)
+{
+ unsigned int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 14; /* bit no 14 becomes 0 */
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x01) {
+ class = 2; /*different orbit */
+ }
+ else {
+ class = 1; /*same orbit */
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc500a.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc500a.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc500a.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,39 @@
+/* MODLAND QA Bits 500m long int bits[0-1]
+ * 00 -> class 1: Corrected product produced at ideal quality -- all bands
+ * 01 -> class 2: Corrected product produced at less than idel quality -- some or all bands
+ * 10 -> class 3: Corrected product NOT produced due to cloud effect -- all bands
+ * 11 -> class 4: Corrected product NOT produced due to other reasons -- some or all bands mayb be fill value (Note that a value of [11] overrides a value of [01])
+ */
+int qc500a(long int pixel)
+{
+ long int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x02) {
+
+ /* Non-Corrected product */
+ if (qctemp & 0x01) {
+ class = 4; /*other reasons */
+ }
+ else {
+ class = 3; /*because of clouds */
+ }
+ }
+ else {
+
+ /* Corrected product */
+ if (qctemp & 0x01) {
+ class = 2; /*less than ideal quality */
+ }
+ else {
+ class = 1; /*ideal quality */
+ }
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc500c.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc500c.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc500c.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,104 @@
+/* Band-wise Data Quality 500m long Int
+ * bits[2-5][6-9][10-13][14-17][18-21][22-25][26-29]
+ * 0000 -> class 1: highest quality
+ * 0111 -> class 2: noisy detector
+ * 1000 -> class 3: dead detector; data interpolated in L1B
+ * 1001 -> class 4: solar zenith >= 86 degrees
+ * 1010 -> class 5: solar zenith >= 85 and < 86 degrees
+ * 1011 -> class 6: missing input
+ * 1100 -> class 7: internal constant used in place of climatological data for at least one atmospheric constant
+ * 1101 -> class 8: correction out of bounds, pixel constrained to extreme allowable value
+ * 1110 -> class 9: L1B data faulty
+ * 1111 -> class 10: not processed due to deep ocean or cloud
+ * Class 11: Combination of bits unused
+ */
+int qc500c(long int pixel, int bandno)
+{
+ long int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 2 + (4 * bandno - 1); /* bitshift [] to [0-3] etc. */
+ swab(&swabfrom, &swabto, 4);
+ qctemp = swabto;
+ if (qctemp & 0x08) { /* 1??? */
+ if (qctemp & 0x04) { /* 11?? */
+ if (qctemp & 0x02) { /* 111? */
+ if (qctemp & 0x01) { /* 1111 */
+ class = 10; /*Not proc.ocean/cloud */
+ }
+ else { /* 1110 */
+ class = 9; /*L1B faulty data */
+ }
+ }
+ else { /* 110? */
+ if (qctemp & 0x01) { /* 1101 */
+ class = 8; /*corr. out of bounds */
+ }
+ else { /* 1100 */
+ class = 7; /*internal constant used */
+ }
+ }
+ }
+ else {
+ if (qctemp & 0x02) { /* 101? */
+ if (qctemp & 0x01) { /* 1011 */
+ class = 6; /*missing input */
+ }
+ else { /* 1010 */
+ class = 5; /*solarzen>=85&&<86 */
+ }
+ }
+ else { /* 100? */
+ if (qctemp & 0x01) { /* 1001 */
+ class = 4; /*solar zenith angle>=86 */
+ }
+ else { /* 1000 */
+ class = 3; /*Dead detector */
+ }
+ }
+ }
+ }
+ else { /* 0??? */
+ if (qctemp & 0x04) { /* 01?? */
+ if (qctemp & 0x02) { /* 011? */
+ if (qctemp & 0x01) { /* 0111 */
+ class = 2; /*Noisy detector */
+ }
+ else { /* 0110 */
+ class = 11; /*Unused */
+ }
+ }
+ else { /* 010? */
+ if (qctemp & 0x01) { /* 0101 */
+ class = 11; /*Unused */
+ }
+ else { /* 0100 */
+ class = 11; /*Unused */
+ }
+ }
+ }
+ else { /* 00?? */
+ if (qctemp & 0x02) { /* 001? */
+ if (qctemp & 0x01) { /* 0011 */
+ class = 11; /*Unused */
+ }
+ else { /* 0010 */
+ class = 11; /*Unused */
+ }
+ }
+ else { /* 000? */
+ if (qctemp & 0x01) { /* 0001 */
+ class = 11; /*Unused */
+ }
+ else { /* 0000 */
+ class = 1; /*Highest quality */
+ }
+ }
+ }
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc500d.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc500d.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc500d.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,24 @@
+/* Atmospheric correction 500m long Int bit[30]
+ * 00 -> class 1: Not Corrected product
+ * 01 -> class 2: Corrected product
+ */
+int qc500d(long int pixel)
+{
+ long int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 30; /* bit no 30 becomes 0 */
+ swab(&swabfrom, &swabto, 1);
+ qctemp = swabto;
+ if (qctemp & 0x01) {
+ class = 2; /*Corrected */
+ }
+ else {
+ class = 1; /*Not corrected */
+ }
+ return class;
+}
+
+
Added: grass/trunk/imagery/i.qc.modis/qc500e.c
===================================================================
--- grass/trunk/imagery/i.qc.modis/qc500e.c (rev 0)
+++ grass/trunk/imagery/i.qc.modis/qc500e.c 2008-09-27 15:11:05 UTC (rev 33564)
@@ -0,0 +1,24 @@
+/* Adjacency correction 500m long Int bit[31]
+ * 00 -> class 1: Not Corrected product
+ * 01 -> class 2: Corrected product
+ */
+int qc500e(long int pixel)
+{
+ long int swabfrom, swabto, qctemp;
+
+ int class;
+
+ swabfrom = pixel;
+ swabfrom >> 31; /* bit no 31 becomes 0 */
+ swab(&swabfrom, &swabto, 1);
+ qctemp = swabto;
+ if (qctemp & 0x01) {
+ class = 2; /*Corrected */
+ }
+ else {
+ class = 1; /*Not corrected */
+ }
+ return class;
+}
+
+
More information about the grass-commit
mailing list