[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