Geografic-Information-System/imagery/i.aster.toar/main.c

/****************************************************************************
 *
 * MODULE:       i.aster.toar
 * AUTHOR(S):    Yann Chemin - yann.chemin@gmail.com
 * PURPOSE:      Calculate TOA Reflectance for Aster from DN.
 * 		 Input 9 bands (VNIR and SWIR).
 *
 * COPYRIGHT:    (C) 2002-2010 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/raster.h>
#include <grass/glocale.h>

#define MAXFILES 15

/* DN to radiance conversion factors */
double gain_aster(int band_number, int gain_code);

    /*Gain Code */
    /*0 - High (Not Applicable for band 10-14: TIR) */
    /*1 - Normal */
    /*2 - Low 1(Not Applicable for band 10-14: TIR) */
    /*3 - Low 2(Not Applicable for Band 1-3N/B & 10-14) */

/*sun exo-atmospheric irradiance */
#define KEXO1 1828.0
#define KEXO2 1559.0
#define KEXO3 1045.0
#define KEXO4 226.73
#define KEXO5 86.50
#define KEXO6 81.99
#define KEXO7 74.72
#define KEXO8 66.41
#define KEXO9 59.83

#define PI 3.1415926

double rad2ref_aster(double radiance, double doy, double sun_elevation,
		     double k_exo);

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 Option *input1, *input2;
    struct Flag *flag0, *flag1, *flag2;
    struct Flag *flag3, *flag4, *flag5;
    struct History history;	/*metadata */

    /************************************/
    char *name;			/*input raster name */
    char *result;		/*output raster name */

    /*Prepare new names for output files */
    char result0[GNAME_MAX], result1[GNAME_MAX];
    char result2[GNAME_MAX], result3[GNAME_MAX];
    char result4[GNAME_MAX], result5[GNAME_MAX];
    char result6[GNAME_MAX], result7[GNAME_MAX];
    char result8[GNAME_MAX], result9[GNAME_MAX];
    char result10[GNAME_MAX], result11[GNAME_MAX];
    char result12[GNAME_MAX], result13[GNAME_MAX];
    char result14[GNAME_MAX];

    /*File Descriptors */
    int infd[MAXFILES];
    int outfd[MAXFILES];
    char **names, **ptr;
    /* For some strange reason infd[0] cannot be used later */
    /* So nfiles is initialized with nfiles = 1 */
    int nfiles = 1;
    int i = 0, j = 0;
    int radiance = 0;
    void *inrast[MAXFILES];
    DCELL *outrast[MAXFILES];
    RASTER_MAP_TYPE in_data_type[MAXFILES];
    RASTER_MAP_TYPE out_data_type = DCELL_TYPE;	/* 0=numbers  1=text */
    double gain[MAXFILES], offset[MAXFILES];
    double kexo[MAXFILES];
    double doy, sun_elevation;

    /************************************/
    G_gisinit(argv[0]);
    module = G_define_module();
    G_add_keyword(_("imagery"));
    G_add_keyword(_("radiometric conversion"));
    G_add_keyword(_("Terra-ASTER"));
    G_add_keyword(_("radiance"));
    G_add_keyword(_("reflectance"));
    G_add_keyword(_("brightness temperature"));
    module->description =
	_("Calculates Top of Atmosphere Radiance/Reflectance/Brightness Temperature from ASTER DN.\n");

    /* Define the different options */
    input = G_define_standard_option(G_OPT_R_INPUTS);
    input->description = _("Names of ASTER DN layers (15 layers)");

    input1 = G_define_option();
    input1->key = "dayofyear";
    input1->type = TYPE_DOUBLE;
    input1->required = YES;
    input1->gisprompt = "value";
    input1->description = _("Day of Year of satellite overpass [0-366]");

    input2 = G_define_option();
    input2->key = "sun_elevation";
    input2->type = TYPE_DOUBLE;
    input2->required = YES;
    input2->gisprompt = "value";
    input2->description = _("Sun elevation angle (degrees, < 90.0)");

    output = G_define_standard_option(G_OPT_R_OUTPUT);
    output->description = _("Base name of the output layers (will add .x)");

    /* Define the different flags */
    flag0 = G_define_flag();
    flag0->key = 'r';
    flag0->description = _("Output is radiance (W/m2)");

    flag1 = G_define_flag();
    flag1->key = 'a';
    flag1->description = _("VNIR is High Gain");

    flag2 = G_define_flag();
    flag2->key = 'b';
    flag2->description = _("SWIR is High Gain");

    flag3 = G_define_flag();
    flag3->key = 'c';
    flag3->description = _("VNIR is Low Gain 1");

    flag4 = G_define_flag();
    flag4->key = 'd';
    flag4->description = _("SWIR is Low Gain 1");

    flag5 = G_define_flag();
    flag5->key = 'e';
    flag5->description = _("SWIR is Low Gain 2");

    /********************/
    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);

    names = input->answers;
    ptr = input->answers;
    doy = atof(input1->answer);
    sun_elevation = atof(input2->answer);
    result = output->answer;

    radiance = (flag0->answer);

    /********************/
    /*Prepare the output file names */

    /********************/
    sprintf(result0,"%s%s", result, ".1");
    sprintf(result1,"%s%s", result, ".2");
    sprintf(result2,"%s%s", result, ".3N");
    sprintf(result3,"%s%s", result, ".3B");
    sprintf(result4,"%s%s", result, ".4");
    sprintf(result5,"%s%s", result, ".5");
    sprintf(result6,"%s%s", result, ".6");
    sprintf(result7,"%s%s", result, ".7");
    sprintf(result8,"%s%s", result, ".8");
    sprintf(result9,"%s%s", result, ".9");
    sprintf(result10,"%s%s", result, ".10");
    sprintf(result11,"%s%s", result, ".11");
    sprintf(result12,"%s%s", result, ".12");
    sprintf(result13,"%s%s", result, ".13");
    sprintf(result14,"%s%s", result, ".14");
    /********************/
    /*Prepare radiance boundaries */

    /********************/
    int gain_code = 1;

    for (i = 0; i < MAXFILES; i++) {
	/*0 - High (Not Applicable for band 10-14: TIR) */
	/*1 - Normal */
	/*2 - Low 1(Not Applicable for band 10-14: TIR) */
	/*3 - Low 2(Not Applicable for Band 1-3N/B & 10-14) */
	if (flag1->answer && i <= 3)
	    gain_code = 0;
	if (flag2->answer && i >= 4 && i <= 9)
	    gain_code = 0;
	if (flag3->answer && i <= 3)
	    gain_code = 2;
	if (flag4->answer && i >= 4 && i <= 9)
	    gain_code = 2;
	if (flag5->answer && i >= 4 && i <= 9)
	    gain_code = 3;
	gain[i] = gain_aster(i, gain_code);
	/* Reset to NORMAL GAIN */
	gain_code = 1;
    }

    /********************/
    /*Prepare sun exo-atm irradiance */

    /********************/
    kexo[0] = KEXO1;
    kexo[1] = KEXO2;
    kexo[2] = KEXO3;
    kexo[3] = KEXO3;
    kexo[4] = KEXO4;
    kexo[5] = KEXO5;
    kexo[6] = KEXO6;
    kexo[7] = KEXO7;
    kexo[8] = KEXO8;
    kexo[9] = KEXO9;

    /********************/

    /********************/
    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-1] = Rast_map_type(name, "");
	/* For some strange reason infd[0] cannot be used later */
	/* So nfiles is initialized with nfiles = 1 */
	infd[nfiles] = Rast_open_old(name, "");

	Rast_get_cellhd(name, "", &cellhd);
	inrast[nfiles-1] = Rast_allocate_buf(in_data_type[nfiles-1]);
	nfiles++;
    }
    nfiles--;
    if (nfiles < MAXFILES) 
	G_fatal_error(_("The input band number should be 15"));

    /***************************************************/
    /* Allocate output buffer, use input map data_type */
    nrows = Rast_window_rows();
    ncols = Rast_window_cols();
    out_data_type = DCELL_TYPE;
    for (i = 0; i < MAXFILES; i++)
	outrast[i] = Rast_allocate_buf(out_data_type);

    outfd[1] = Rast_open_new(result0, 1);
    outfd[2] = Rast_open_new(result1, 1);
    outfd[3] = Rast_open_new(result2, 1);
    outfd[4] = Rast_open_new(result3, 1);
    outfd[5] = Rast_open_new(result4, 1);
    outfd[6] = Rast_open_new(result5, 1);
    outfd[7] = Rast_open_new(result6, 1);
    outfd[8] = Rast_open_new(result7, 1);
    outfd[9] = Rast_open_new(result8, 1);
    outfd[10] = Rast_open_new(result9, 1);
    outfd[11] = Rast_open_new(result10, 1);
    outfd[12] = Rast_open_new(result11, 1);
    outfd[13] = Rast_open_new(result12, 1);
    outfd[14] = Rast_open_new(result13, 1);
    outfd[15] = Rast_open_new(result14, 1);
    /* Process pixels */

    DCELL dout[MAXFILES];
    DCELL d[MAXFILES];

    for (row = 0; row < nrows; row++) {
	G_percent(row, nrows, 2);
	/* read input map */
	for (i = 1; i <= MAXFILES; i++)
	    Rast_get_row(infd[i], inrast[i-1], row, in_data_type[i-1]);

	/*process the data */
	for (col = 0; col < ncols; col++) {
	    for (i = 0; i < MAXFILES; 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 radiance mode or Thermal band */
		if (radiance || i >= 10) {
		    dout[i] = gain[i] * (d[i] - 1.0);
		}
		/* if reflectance default mode and Not Thermal Band */
		else {
		    dout[i] = gain[i] * (d[i] - 1.0);
		    dout[i] =
			rad2ref_aster(dout[i], doy, sun_elevation, kexo[i]);
		}
		outrast[i][col] = dout[i];
	    }
	}
	for (i = 1; i <= MAXFILES; i++)
	    Rast_put_row(outfd[i], outrast[i-1], out_data_type);
    }
    for (i = 1; i <= MAXFILES; i++) {
	G_free(inrast[i-1]);
	Rast_close(infd[i]);
	G_free(outrast[i-1]);
	Rast_close(outfd[i]);
    }
    exit(EXIT_SUCCESS);
}