i.landsat.toar: sync to relbr6

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@56905 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.landsat.toar/i.landsat.toar.html

@@ -8,15 +8,16 @@ radiance or reflectance with atmospheric correction (DOS method).
 
 <p>
 Usually, to do so the production date, the acquisition date, and the
-solar elevation is needed. Moreover, for Landsat-7 ETM+ it is also
+solar elevation are needed. Moreover, for Landsat-7 ETM+ it is also
 needed the gain (high or low) of the nine respective bands.
 
 <p>
-Optionally (recommended), the data can be read from metadata file (.met or MTL.txt) for all
-Landsat MSS, TM and ETM+. However, if the solar elevation or the
-product creation date are given the values of the metadata file are
-overwriten. This is necessary when the data in the .met file is
-incorrect or not accurate.
+Optionally (recommended), the data can be read from metadata file
+(.met or MTL.txt) for all Landsat MSS, TM, ETM+ and OLI/TIRS. However,
+if the solar elevation is given the value of the metadata file are
+overwritten. This is necessary when the data in the .met file is
+incorrect or not accurate. Also, if acquisition or production dates are
+not found in the metadata file then the command line values are used.
 
 <p>
 <b>Attention</b>: Any null value or smaller than QCALmin in the input
@@ -138,6 +139,7 @@ correct calibration constants, which were changed in the dates:
   <li>Landsat-5 MSS: April 6, 1984 and November 9, 1984</li>
   <li>Landsat-5 TM:  May 4, 2003 and April, 2 2007</li>
   <li>Landsat-7 ETM+: July 1, 2000</li>
+  <li>Landsat-8 OLI/TIRS: launched in 2013</li>
 </ul>
 
 <h2>EXAMPLES</h2>

imagery/i.landsat.toar/landsat.c

@@ -5,9 +5,9 @@
 
 #include "landsat.h"
 
-#define PI   3.1415926535897932384626433832795
-#define R2D 57.295779513082320877
-#define D2R  0.017453292519943295769
+#define PI   M_PI
+#define R2D  180. / M_PI
+#define D2R  M_PI / 180.
 
 /****************************************************************************
  * PURPOSE: Calibrated Digital Number to at-satellite Radiance
@@ -22,7 +22,7 @@ double lsat_qcal2rad(double qcal, band_data * band)
  *****************************************************************************/
 double lsat_rad2ref(double rad, band_data * band)
 {
-    return (double)(rad / band->K2);
+    return (double)(rad / band->K1);
 }
 
 /****************************************************************************
@@ -62,10 +62,10 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
     sin_e = (double)(sin(D2R * lsat->sun_elev));
     cos_v = (double)(cos(D2R * (lsat->number < 4 ? 9.2 : 8.2)));
 
-	/** Global irradiance on the sensor.
-		Radiance to reflectance coefficient, only NO thermal bands.
-		K1 and K2 variables are also utilized as thermal constants
-	*/
+    /** Global irradiance on the sensor.
+	Radiance to reflectance coefficient, only NO thermal bands.
+	K1 and K2 variables are also utilized as thermal constants
+    */
     if (lsat->band[i].thermal == 0) {
 	switch (method) {
 	case DOS2:
@@ -129,16 +129,17 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 	    break;
 	}
 	rad_sun = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
-	G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz,
-			  Edown);
+	if (method > DOS)
+	    G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
+			      TAUv, TAUz, Edown);
 
-	lsat->band[i].K1 = 0.;
-	lsat->band[i].K2 = rad_sun;
+	lsat->band[i].K1 = rad_sun;
+	lsat->band[i].K2 = 0.;
     }
 
-	/** Digital number to radiance coefficients.
-		Whitout atmospheric calibration for thermal bands.
-	*/
+    /** Digital number to radiance coefficients.
+	Without atmospheric calibration for thermal bands.
+    */
     lsat->band[i].gain = ((lsat->band[i].lmax - lsat->band[i].lmin) /
 			  (lsat->band[i].qcalmax - lsat->band[i].qcalmin));
 

imagery/i.landsat.toar/landsat.h

@@ -30,14 +30,14 @@ typedef struct
 
     double wavemax, wavemin;	/* Wavelength in µm              */
 
+    double esun;		/* Mean solar irradiance         */
     double lmax, lmin;		/* Spectral radiance             */
     double qcalmax, qcalmin;	/* Quantized calibrated pixel    */
-    double esun;		/* Mean solar irradiance         */
 
     char thermal;		/* Flag to thermal band          */
     double gain, bias;		/* Gain and Bias of sensor       */
-    double K1, K2;		/* Thermal calibration constants,
-				   or Rad2Ref constants          */
+    double K1, K2;		/* Thermal calibration or
+				   Rad2Ref constants             */
 
 } band_data;
 
@@ -48,9 +48,10 @@ typedef struct
 
     char creation[11];		/* Image production date         */
     char date[11];		/* Image acquisition date        */
+
     double dist_es;		/* Distance Earth-Sun            */
-    double sun_elev;		/* Solar elevation               */
-    double sunaz;		/* Solar Azimuth                 */
+    double sun_elev;		/* Sun elevation                 */
+    double sun_az;		/* Sun Azimuth                   */
     double time;		/* Image Acquisition Time        */
 
     char sensor[10];		/* Type of sensor: MSS, TM, ETM+, OLI/TIRS */

imagery/i.landsat.toar/landsat_met.c

@@ -1,6 +1,6 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 #include <math.h>
 
 #include <grass/gis.h>
@@ -9,44 +9,55 @@
 #include "local_proto.h"
 #include "earth_sun.h"
 
+#define PI   M_PI
+#define D2R  M_PI / 180.
+
 #define MAX_STR         127
-#define METADATA_SIZE   65535  /* MTL.txt file size  65535 bytes */
-#define TM5_MET_SIZE    28700  /* .met file size 28686 bytes */
+#define METADATA_SIZE   65535	/* MTL.txt file size  65535 bytes */
+#define TM5_MET_SIZE    28700	/* .met file size 28686 bytes */
 
 
 inline void chrncpy(char *dest, char src[], int n)
 {
     int i;
-    for( i = 0 ; i < n && src[i] != '\0' && src[i] != '\"'; i++) dest[i] = src[i];
+
+    for (i = 0; i < n && src[i] != '\0' && src[i] != '\"'; i++)
+	dest[i] = src[i];
     dest[i] = '\0';
 }
 
-void get_oldformat( const char metadata[], char * key, char value[] )
+void get_metformat(const char metadata[], char *key, char value[])
 {
     int i = 0;
-    char * ptr = strstr(metadata, key);
-    if (ptr != NULL)
-    {
-        ptr = strstr(ptr, " VALUE ");
-        if (ptr != NULL)
-        {
+    char *ptrmet = strstr(metadata, key);
+    char *ptr;
+
+    if (ptrmet != NULL) {
+        ptr = strstr(ptrmet, " VALUE ");
+        if (ptr != NULL) {
             while (*ptr++ != '=') ;
-            while (*ptr <= ' ' || *ptr == '\"') *ptr++;
-            while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ') value[i++] = *ptr++;
+            while (*ptr <= ' ' || *ptr == '\"')
+		ptr++;
+            while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
+		value[i++] = *ptr++;
         }
     }
     value[i] = '\0';
 }
 
-void get_newformat( const char metadata[], char * key, char value[] )
+void get_mtlformat(const char metadata[], char *key, char value[])
 {
     int i = 0;
-    char * ptr = strstr(metadata, key);
-    if (ptr != NULL)
-    {
+    char *ptr = strstr(metadata, key);
+
+    if (ptr != NULL) {
         while (*ptr++ != '=') ;
-        while (*ptr <= ' ' || *ptr == '\"') *ptr++;
-        while (i < MAX_STR && *ptr != '\"' && *ptr > ' ') value[i++] = *ptr++;
+
+        while (*ptr <= ' ' || *ptr == '\"')
+	    ptr++;
+
+        while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
+	    value[i++] = *ptr++;
     }
     value[i] = '\0';
 }
@@ -56,160 +67,305 @@ void get_newformat( const char metadata[], char * key, char value[] )
  * PURPOSE:     Read parameters from Landsat metadata files
  *****************************************************************************/
 
-void lsat_metadata( char * metafile, lsat_data * lsat)
+void lsat_metadata(char *metafile, lsat_data * lsat)
 {
     FILE *f;
-    char mtldata[METADATA_SIZE];
+    char metadata[METADATA_SIZE];
     char key[MAX_STR], value[MAX_STR];
-    void (*get_mtldata)( const char [], char *, char [] );
-    int i, j, version;
+    void (*get_metadata) (const char [], char *, char []);
+    int i, j, ver_meta;
 
     if ((f = fopen(metafile, "r")) == NULL)
-       G_fatal_error(_("Metadata file <%s> not found"), metafile);
+	G_fatal_error(_("Metadata file <%s> not found"), metafile);
+
+    i = fread(metadata, METADATA_SIZE, 1, f);
 
-    i = fread(mtldata, METADATA_SIZE, 1, f);
-    get_mtldata = (strstr(mtldata, " VALUE ") != NULL) ? get_oldformat : get_newformat;
-    version = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
+    /* get version of the metadata file */
+    get_metadata = 
+	(strstr(metadata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
+    ver_meta = (strstr(metadata, "QCALMAX_BAND") != NULL) ? 0 : 1;
     
     /* Fill with product metadata */
-    get_mtldata(mtldata, "SPACECRAFT_ID", value);
-    if( value[0] == '\0' )
-    {
-        get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
+    get_metadata(metadata, "SPACECRAFT_ID", value);
+    if (value[0] == '\0') {
+        get_metadata(metadata, "PLATFORMSHORTNAME", value);
     }
     
-    i = 0; while( value[i] && (value[i] < '0' || value[i] > '9') ) i++;
+    i = 0;
+    while (value[i] && (value[i] < '0' || value[i] > '9'))
+	i++;
     lsat->number = atoi(value + i);
 
-    get_mtldata(mtldata, "SENSOR_ID", value);
-    if( value[0] == '\0' )
-    {
-        get_mtldata(mtldata, "SENSORSHORTNAME", value);
+    get_metadata(metadata, "SENSOR_ID", value);
+    if (value[0] == '\0') {
+        get_metadata(metadata, "SENSORSHORTNAME", value);
     }
     chrncpy(lsat->sensor, value, 8);
 
-    get_mtldata(mtldata, "DATE_ACQUIRED", value);
-    if( value[0] == '\0' )
-    {
-        get_mtldata(mtldata, "ACQUISITION_DATE", value);
-        if( value[0] == '\0' )
-        {
-            get_mtldata(mtldata, "CALENDARDATE", value);
+    get_metadata(metadata, "DATE_ACQUIRED", value);
+    if (value[0] == '\0') {
+        get_metadata(metadata, "ACQUISITION_DATE", value);
+        if (value[0] == '\0') {
+            get_metadata(metadata, "CALENDARDATE", value);
         }
     }
     chrncpy(lsat->date, value, 10);
 
-    get_mtldata(mtldata, "FILE_DATE", value);
-    if( value[0] == '\0' )
-    {
-        get_mtldata(mtldata, "CREATION_TIME", value);
-        if( value[0] == '\0' )
-        {
-            get_mtldata(mtldata, "PRODUCTIONDATETIME", value);
+    get_metadata(metadata, "FILE_DATE", value);
+    if (value[0] == '\0') {
+        get_metadata(metadata, "CREATION_TIME", value);
+        if (value[0] == '\0') {
+            get_metadata(metadata, "PRODUCTIONDATETIME", value);
         }
     }
-    chrncpy(lsat->creation, value, 10);
+    if (value[0] != '\0')
+	chrncpy(lsat->creation, value, 10);
+    else
+	G_warning
+	    ("Using production date from the command line 'product_date'");
 
-    get_mtldata(mtldata, "SUN_AZIMUTH", value);
-    lsat->sunaz = atof(value);
-    if( lsat->sunaz == 0. )
-        G_warning("Sun azimuth is %f", lsat->sunaz);
+    get_metadata(metadata, "SUN_AZIMUTH", value);
+    lsat->sun_az = atof(value);
+    if (lsat->sun_az == 0.)
+        G_warning("Sun azimuth is %f", lsat->sun_az);
 
-    get_mtldata(mtldata, "SUN_ELEVATION", value);
-    if( value[0] == '\0' )
-    {
-        get_mtldata(mtldata, "SolarElevation", value);
+    get_metadata(metadata, "SUN_ELEVATION", value);
+    if (value[0] == '\0') {
+        get_metadata(metadata, "SolarElevation", value);
     }
     lsat->sun_elev = atof(value);
-    if( lsat->sun_elev == 0. )
+    if (lsat->sun_elev == 0.)
         G_warning("Sun elevation is %f", lsat->sun_elev);
 
-    get_mtldata(mtldata, "SCENE_CENTER_TIME", value);
-    if( value[0] == '\0' )
-    {
-        get_mtldata(mtldata, "SCENE_CENTER_SCAN_TIME", value);
+    get_metadata(metadata, "SCENE_CENTER_TIME", value);
+    if (value[0] == '\0') {
+        get_metadata(metadata, "SCENE_CENTER_SCAN_TIME", value);
     }
-    /*Thanks Markus Metz !
-    Remove trailing 'z'*/
+    /* Remove trailing 'z'*/
     value[strlen(value) - 1]='\0';
     /* Cast from hh:mm:ss into hh.hhh*/
     G_llres_scan(value, &lsat->time);
-    if( lsat->time == 0. )
+    if (lsat->time == 0.)
         G_warning("Time is %f", lsat->time);
 
-    /* Fill data with the basic sensor parameters */
-    switch(lsat->number)
-    {
-        case 1:
-            set_MSS1(lsat);
-            break;
-        case 2:
-            set_MSS2(lsat);
-            break;
-        case 3:
-            set_MSS3(lsat);
-            break;
-        case 4:
-            if (lsat->sensor[0] == 'M')
-                set_MSS4(lsat);
-            else
-                set_TM4(lsat);
-            break;
-        case 5:
-            if (lsat->sensor[0] == 'M')
-                set_MSS5(lsat);
-            else
-                set_TM5(lsat);
-            break;
-        case 7:
-        {
-            char * fmt_gain[] = { "BAND%d_GAIN%d ", " GAIN_BAND_%d_VCID_%d" };
-            for( i = 1, j = 0; i < 9; i++ )
-            {
-                sprintf(key, fmt_gain[version], i, 1);
-                if (i != 6) key[version == 0 ? 10 : 12] = '\0';
-                get_mtldata(mtldata, key,  value + j++);
-                if (i == 6)
-                {
-                    sprintf(key, fmt_gain[version], i, 2);
-                    get_mtldata(mtldata, key,  value + j++);
-                }
-            }
-            value[j] = '\0';
-            G_debug(1, "ETM+ gain = [%s]", value);
-            set_ETM(lsat, value);
-            break;
-        }
-        default:
-            G_warning("Unable to recognize satellite platform [%d]", lsat->number);
-            break;
+    /* Fill the data with the basic/default sensor parameters */
+    switch (lsat->number) {
+    case 1:
+	set_MSS1(lsat);
+	break;
+    case 2:
+	set_MSS2(lsat);
+	break;
+    case 3:
+	set_MSS3(lsat);
+	break;
+    case 4:
+	if (lsat->sensor[0] == 'M')
+	    set_MSS4(lsat);
+	else
+	    set_TM4(lsat);
+	break;
+    case 5:
+	if (lsat->sensor[0] == 'M')
+	    set_MSS5(lsat);
+	else
+	    set_TM5(lsat);
+	break;
+    case 7:
+	{
+	    char *fmt_gain[] = { "BAND%d_GAIN%d", " GAIN_BAND_%d_VCID_%d" };
+	    for (i = 1, j = 0; i < 9; i++) {
+		sprintf(key, fmt_gain[ver_meta], i, 1);
+		if (i != 6)
+		    key[ver_meta == 0 ? 10 : 12] = '\0';
+		get_metadata(metadata, key, value + j++);
+		if (i == 6) {
+		    sprintf(key, fmt_gain[ver_meta], i, 2);
+		    get_metadata(metadata, key, value + j++);
+		}
+	    }
+	    value[j] = '\0';
+	    G_debug(1, "ETM+ gain = [%s]", value);
+	    set_ETM(lsat, value);
+	    break;
+	}
+    case 8:
+	set_LDCM(lsat);
+	break;
+
+    default:
+	G_warning("Unable to recognize satellite platform [%d]",
+		  lsat->number);
+	break;
     }
 
     /* Update the information from metadata file, if infile */
-    /* if( format == NEW_FORMAT ) */
-    if( get_mtldata == get_newformat )
-    {
-        char * fmt_lmax[]     = { "LMAX_BAND%d",    "RADIANCE_MAXIMUM_BAND_%d" };
-        char * fmt_lmin[]     = { "LMIN_BAND%d",    "RADIANCE_MINIMUM_BAND_%d" };
-        char * fmt_qcalmax[]  = { "QCALMAX_BAND%d", "QUANTIZE_CAL_MAX_BAND_%d" };
-        char * fmt_qcalmin[]  = { "QCALMIN_BAND%d", "QUANTIZE_CAL_MIN_BAND_%d" };
-
-        for (i = 0; i < lsat->bands; i++) {
-            sprintf(key, fmt_lmax[version], lsat->band[i].code);
-            get_mtldata(mtldata, key, value);
-            lsat->band[i].lmax = atof(value);
-            sprintf(key, fmt_lmin[version], lsat->band[i].code);
-            get_mtldata(mtldata, key, value);
-            lsat->band[i].lmin = atof(value);
-            sprintf(key, fmt_qcalmax[version], lsat->band[i].code);
-            get_mtldata(mtldata, key, value);
-            lsat->band[i].qcalmax = atof(value);
-            sprintf(key, fmt_qcalmin[version], lsat->band[i].code);
-            get_mtldata(mtldata, key, value);
-            lsat->band[i].qcalmin = atof(value);
-        }
+    if (get_metadata == get_mtlformat) {
+	if (ver_meta == 0) {	/* now MTLold.txt */
+	    G_verbose_message("Metada file is MTL file: old format");
+	    for (i = 0; i < lsat->bands; i++) {
+		sprintf(key, "LMAX_BAND%d", lsat->band[i].code);
+		get_metadata(metadata, key, value);
+		lsat->band[i].lmax = atof(value);
+		sprintf(key, "LMIN_BAND%d", lsat->band[i].code);
+		get_metadata(metadata, key, value);
+		lsat->band[i].lmin = atof(value);
+		sprintf(key, "QCALMAX_BAND%d", lsat->band[i].code);
+		get_metadata(metadata, key, value);
+		lsat->band[i].qcalmax = atof(value);
+		sprintf(key, "QCALMIN_BAND%d", lsat->band[i].code);
+		get_metadata(metadata, key, value);
+		lsat->band[i].qcalmin = atof(value);
+	    }
+	}
+	else {			/* now MTL.txt */
+
+	    G_verbose_message("Metada file is MTL file: new format");
+	    if (strstr(metadata, "RADIANCE_MAXIMUM_BAND") != NULL) {
+		G_verbose_message
+		    ("RADIANCE & QUANTIZE from data of the metadata file");
+		for (i = 0; i < lsat->bands; i++) {
+		    if (lsat->band[i].thermal && lsat->number == 7) {
+			sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
+				lsat->band[i].code - 60);
+			get_metadata(metadata, key, value);
+			lsat->band[i].lmax = atof(value);
+			sprintf(key, "RADIANCE_MINIMUM_BAND_6_VCID_%d",
+				lsat->band[i].code - 60);
+			get_metadata(metadata, key, value);
+			lsat->band[i].lmin = atof(value);
+			sprintf(key, "QUANTIZE_CAL_MAX_BAND_6_VCID_%d",
+				lsat->band[i].code - 60);
+			get_metadata(metadata, key, value);
+			lsat->band[i].qcalmax = atof(value);
+			sprintf(key, "QUANTIZE_CAL_MIN_BAND_6_VCID_%d",
+				lsat->band[i].code - 60);
+			get_metadata(metadata, key, value);
+			lsat->band[i].qcalmin = atof(value);
+		    }
+		    else {
+			sprintf(key, "RADIANCE_MAXIMUM_BAND_%d",
+				lsat->band[i].code);
+			get_metadata(metadata, key, value);
+			lsat->band[i].lmax = atof(value);
+			sprintf(key, "RADIANCE_MINIMUM_BAND_%d",
+				lsat->band[i].code);
+			get_metadata(metadata, key, value);
+			lsat->band[i].lmin = atof(value);
+			sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d",
+				lsat->band[i].code);
+			get_metadata(metadata, key, value);
+			lsat->band[i].qcalmax = atof(value);
+			sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d",
+				lsat->band[i].code);
+			get_metadata(metadata, key, value);
+			lsat->band[i].qcalmin = atof(value);
+		    }
+		}
+	    }
+	    else {
+		G_verbose_message
+		    ("RADIANCE & QUANTIZE from radiometric rescaling group of the metadata file");
+		/* from LDCM sample file: mode = 0; from LDCM-DFCB-004.pdf file: mode = 1 */
+		int mode =
+		    (strstr(metadata, "RADIANCE_MULTIPLICATIVE_FACTOR_BAND") !=
+		     NULL) ? 0 : 1;
+		char *fmt_radmu[] =
+		    { "RADIANCE_MULTIPLICATIVE_FACTOR_BAND%d",
+	"RADIANCE_MULT_BAND_%d" };
+		char *fmt_radad[] =
+		    { "RADIANCE_ADDITIVE_FACTOR_BAND%d",
+	"RADIANCE_ADD_BAND_%d" };
+		char *fmt_k1cte[] =
+		    { "K1_CONSTANT_BAND%d", "K1_CONSTANT_BAND_%d" };
+		char *fmt_k2cte[] =
+		    { "K2_CONSTANT_BAND%d", "K2_CONSTANT_BAND_%d" };
+		char *fmt_refad[] =
+		    { "REFLECTANCE_ADDITIVE_FACTOR_BAND%d",
+	"REFLECTANCE_ADD_BAND_%d" };
+
+		for (i = 0; i < lsat->bands; i++) {
+		    sprintf(key, fmt_radmu[mode], lsat->band[i].code);
+		    get_metadata(metadata, key, value);
+		    lsat->band[i].gain = atof(value);
+		    sprintf(key, fmt_radad[mode], lsat->band[i].code);
+		    get_metadata(metadata, key, value);
+		    lsat->band[i].bias = atof(value);
+		    /* reverse to calculate the original values */
+		    lsat->band[i].qcalmax =
+			(lsat->number < 8 ? 255. : 65535.);
+		    lsat->band[i].qcalmin = 1.;
+		    lsat->band[i].lmin =
+			lsat->band[i].gain * lsat->band[i].qcalmin +
+			lsat->band[i].bias;
+		    lsat->band[i].lmax =
+			lsat->band[i].gain * lsat->band[i].qcalmax +
+			lsat->band[i].bias;
+		    /* ----- */
+		    if (lsat->number == 8) {
+			if (lsat->band[i].thermal) {
+			    sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+			    get_metadata(metadata, key, value);
+			    lsat->band[i].K1 = atof(value);
+			    sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+			    get_metadata(metadata, key, value);
+			    lsat->band[i].K2 = atof(value);
+			}
+			else {
+			    lsat->band[i].K1 = 0.;
+			    /* ESUN from metadafa file: bias/K2 seem better than gain/K1 */
+			    sprintf(key, fmt_refad[mode], lsat->band[i].code);
+			    get_metadata(metadata, key, value);
+			    lsat->band[i].K2 = atof(value);
+			    lsat->band[i].esun =
+				(double)(PI * lsat->dist_es * lsat->dist_es *
+					 lsat->band[i].bias) /
+					 (sin(D2R * lsat->sun_elev) *
+					 lsat->band[i].K2);
+			}
+		    }
+		}
+	    }
+	    /* Other possible values in file */
+	    get_metadata(metadata, "EARTH_SUN_DISTANCE", value);
+	    if (value[0] != '\0') {
+		lsat->dist_es = atof(value);
+		G_verbose_message
+		    ("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
+	    }
+	}
+    }
+    else {
+	G_verbose_message("Metada file is MET file");
+	G_verbose_message
+	    ("RADIANCE & QUANTIZE from band setting of the metadata file");
+	for (i = 0; i < lsat->bands; i++) {
+	    sprintf(key, "Band%dGainSetting", lsat->band[i].code);
+	    get_metadata(metadata, key, value);
+	    if (value[0] == '\0') {
+		G_warning(key);
+		continue;
+	    }
+	    lsat->band[i].gain = atof(value);
+	    sprintf(key, "Band%dBiasSetting", lsat->band[i].code);
+	    get_metadata(metadata, key, value);
+	    if (value[0] == '\0') {
+		G_warning(key);
+		continue;
+	    }
+	    lsat->band[i].bias = atof(value);
+
+	    lsat->band[i].qcalmax = 255.;
+	    lsat->band[i].qcalmin = 1.;
+	    lsat->band[i].lmin =
+		lsat->band[i].gain * lsat->band[i].qcalmin +
+		lsat->band[i].bias;
+	    lsat->band[i].lmax =
+		lsat->band[i].gain * lsat->band[i].qcalmax +
+		lsat->band[i].bias;
+	}
     }
 
-    (void)fclose(f);
+    fclose(f);
     return;
 }

imagery/i.landsat.toar/landsat_set.c

@@ -1,6 +1,6 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 #include <math.h>
 
 #include <grass/gis.h>
@@ -15,10 +15,10 @@ void sensor_MSS(lsat_data * lsat)
 
     /* green, red, near infrared, near infrared */
     int band[] = { 1, 2, 3, 4 };
-    int code[] = { 4, 5, 6, 7 };
+    int code[] = { 4, 5, 6, 7 };	/* corrected for MSS4 and MSS5 to 1,2,3,4 */
     double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
     double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
-    /* original: 79x57, now all 60 */
+    /* original: 79x57, now all 60 m */
 
     strcpy(lsat->sensor, "MSS");
 
@@ -46,7 +46,7 @@ void sensor_TM(lsat_data * lsat)
     /* 30, 30, 30, 30, 30, 120 original, 60 resamples before Feb 25, 2010 and 30 after, 30 */
 
     if (!lsat->sensor)
-      strcpy(lsat->sensor, "TM");
+	strcpy(lsat->sensor, "TM");
 
     lsat->bands = 7;
     for (i = 0; i < lsat->bands; i++) {
@@ -65,11 +65,11 @@ void sensor_ETM(lsat_data * lsat)
 {
     int i;
 
-    /* blue, green red, near infrared, shortwave IR, thermal IR, shortwave IR, panchromatic */
+    /* blue, green, red, near infrared, shortwave IR, thermal IR, shortwave IR, panchromatic */
     int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
     int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
-    double wmin[] = { 0.450, 0.525, 0.630, 0.75, 1.55, 10.40, 2.09, 0.52 };
-    double wmax[] = { 0.515, 0.605, 0.690, 0.90, 1.75, 12.50, 2.35, 0.90 };
+    double wmin[] = { 0.450, 0.525, 0.630, 0.75, 1.55, 10.40, 10.40, 2.09, 0.52 };
+    double wmax[] = { 0.515, 0.605, 0.690, 0.90, 1.75, 12.50, 12.50, 2.35, 0.90 };
     /* 30, 30, 30, 30, 30, 60 (after Feb. 25, 2010: 30), 30, 15 */
 
     strcpy(lsat->sensor, "ETM+");
@@ -87,7 +87,7 @@ void sensor_ETM(lsat_data * lsat)
     return;
 }
 
-void sensor_OLI(lsat_data * lsat)
+void sensor_LDCM(lsat_data * lsat)
 {
     int i;
 
@@ -103,13 +103,13 @@ void sensor_OLI(lsat_data * lsat)
 
     lsat->bands = 11;
     for (i = 0; i < lsat->bands; i++) {
-        lsat->band[i].number = *(band + i);
-        lsat->band[i].code = *(code + i);
-        lsat->band[i].wavemax = *(wmax + i);
-        lsat->band[i].wavemin = *(wmin + i);
-        lsat->band[i].qcalmax = 255.;
-        lsat->band[i].qcalmin = 1.;
-        lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
+	lsat->band[i].number = *(band + i);
+	lsat->band[i].code = *(code + i);
+	lsat->band[i].wavemax = *(wmax + i);
+	lsat->band[i].wavemin = *(wmin + i);
+	lsat->band[i].qcalmax = 255.;
+	lsat->band[i].qcalmin = 1.;
+	lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
     }
     return;
 }
@@ -180,9 +180,9 @@ void set_MSS2(lsat_data * lsat)
 
     julian = julian_char(lsat->creation);
     if (julian < julian_char("1975-07-16"))
-        i = 0;
+	i = 0;
     else
-        i = 1;
+	i = 1;
 
     lmax = Lmax[i];
     lmin = Lmin[i];
@@ -279,11 +279,11 @@ void set_MSS4(lsat_data * lsat)
 
     julian = julian_char(lsat->creation);
     if (julian < julian_char("1982-08-26"))
-	    i = 0;
+	i = 0;
     else if (julian < julian_char("1983-03-31"))
-	    i = 1;
+	i = 1;
     else
-	    i = 2;
+	i = 2;
 
     lmax = Lmax[i];
     lmin = Lmin[i];
@@ -329,11 +329,11 @@ void set_TM4(lsat_data * lsat)
 
     julian = julian_char(lsat->creation);
     if (julian < julian_char("1983-08-01"))
-	    i = 0;
+	i = 0;
     else if (julian < julian_char("1984-01-15"))
-	    i = 1;
+	i = 1;
     else
-	    i = 2;
+	i = 2;
 
     lmax = Lmax[i];
     lmin = Lmin[i];
@@ -436,11 +436,11 @@ void set_TM5(lsat_data * lsat)
 
     julian = julian_char(lsat->creation);
     if (julian < julian_char("2003-05-04"))
-	    i = 0;
+	i = 0;
     else if (julian < julian_char("2007-04-02"))
-	    i = 1;
+	i = 1;
     else
-	    i = 2;
+	i = 2;
 
     lmax = Lmax[i];
     lmin = Lmin[i];
@@ -453,9 +453,9 @@ void set_TM5(lsat_data * lsat)
     }
 
     jbuf = julian_char("2004-04-04");
-    if (julian >= jbuf && !(lsat->flag & METADATAFILE) )
-    {
-        G_warning("Using QCalMin=1.0 as a NLAPS product processed after 04/04/2004");
+    if (julian >= jbuf && !(lsat->flag & METADATAFILE)) {
+	G_warning
+	    ("Using QCalMin=1.0 as a NLAPS product processed after 04/04/2004");
     }
     lsat->number = 5;
     sensor_TM(lsat);
@@ -463,9 +463,9 @@ void set_TM5(lsat_data * lsat)
     lsat->dist_es = earth_sun(lsat->date);
 
     for (i = 0; i < lsat->bands; i++) {
-	j = lsat->band[i].number - 1;
 	if (julian >= jbuf)
 	    lsat->band[i].qcalmin = 1.;
+	j = lsat->band[i].number - 1;
 	lsat->band[i].esun = *(esun + j);
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
@@ -528,7 +528,6 @@ void set_ETM(lsat_data * lsat, char gain[])
 
     for (i = 0; i < lsat->bands; i++) {
 	j = lsat->band[i].number - 1;
-	lsat->band[i].esun = *(esun + j);
 	if (gain[i] == 'H' || gain[i] == 'h') {
 	    lmax = LmaxH[k];
 	    lmin = LminH[k];
@@ -537,6 +536,7 @@ void set_ETM(lsat_data * lsat, char gain[])
 	    lmax = LmaxL[k];
 	    lmin = LminL[k];
 	}
+	lsat->band[i].esun = *(esun + j);
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
 	if (lsat->band[i].thermal) {
@@ -552,42 +552,42 @@ void set_ETM(lsat_data * lsat, char gain[])
  * PURPOSE:     Store values of Landsat-8 OLI/TIRS
  *              February 14, 2013
  *****************************************************************************/
-void set_OLI(lsat_data * lsat)
+void set_LDCM(lsat_data * lsat)
 {
     int i, j;
-    double julian, *lmax, *lmin;
+    double *lmax, *lmin;
 
-    /* Spectral radiances at detector
-    double Lmax[][4] = {
-        {220., 175., 145., 147.},
-        {259., 179., 149., 128.}
+    /* Spectral radiances at detector */
+    
+    /* uncorrected values */ 
+    double Lmax[][11] = {
+	{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
     };
-    double Lmin[][4] = {
-        {4., 3., 3., 1.},
-        {4., 3., 3., 1.}
+    double Lmin[][11] = {
+	{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
     };
-    * Solar exoatmospheric spectral irradiances
-    double esun[] = { 1824., 1570., 1249., 853.4 };
+    /* Solar exoatmospheric spectral irradiances */
+    /* estimates */
+    double esun[] = { 2062., 21931., 1990., 1688., 1037., 268.6, 94.6, 1892., 399.0, 0., 0. };
 
-    lmax = Lmax[i];
-    lmin = Lmin[i];
+    lmax = Lmax[0];
+    lmin = Lmin[0];
 
-    lsat->number = 3;
-    sensor_MSS(lsat);
+    lsat->number = 8;
+    sensor_LDCM(lsat);
 
     lsat->dist_es = earth_sun(lsat->date);
 
     for (i = 0; i < lsat->bands; i++) {
-        j = lsat->band[i].number - 1;
-        lsat->band[i].esun = *(esun + j);
-        lsat->band[i].lmax = *(lmax + j);
-        lsat->band[i].lmin = *(lmin + j);
-        if (lsat->band[i].thermal) {
-            lsat->band[i].K1 = 0;
-            lsat->band[i].K2 = 0;
-        }
+	j = lsat->band[i].number - 1;
+	lsat->band[i].esun = *(esun + j);
+	lsat->band[i].lmax = *(lmax + j);
+	lsat->band[i].lmin = *(lmin + j);
+	if (lsat->band[i].thermal) {
+	    lsat->band[i].K1 = (lsat->band[i].number == 10 ? 666.09 : 697.56);
+	    lsat->band[i].K2 = 1282.71;
+	}
     }
-    */
-    G_debug(1, "Landsat-8 LDCM");
+    G_debug(1, "Landsat-8 OLI/TIRS");
     return;
 }

imagery/i.landsat.toar/local_proto.h

@@ -17,4 +17,6 @@ void set_TM5(lsat_data *);
 
 void set_ETM(lsat_data *, char[]);
 
+void set_LDCM(lsat_data *);
+
 #endif

imagery/i.landsat.toar/main.c

@@ -9,9 +9,9 @@
  *               Adopted for GRASS 7 by Martin Landa <landa.martin gmail.com>
  *
  * PURPOSE:      Calculate TOA Radiance or Reflectance and Kinetic Temperature
- *               for Landsat 1/2/3/4/5 MS, 4/5 TM or 7 ETM+
+ *               for Landsat 1/2/3/4/5 MS, 4/5 TM, 7 ETM+, and 8 OLI/TIRS
  *
- * COPYRIGHT:    (C) 2006-2012 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2006-2013 by the GRASS Development Team
  *
  *               This program is free software under the GNU General
  *               Public License (>=v2). Read the file COPYING that
@@ -42,10 +42,10 @@ int main(int argc, char *argv[])
 
     RASTER_MAP_TYPE in_data_type;
 
-    struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate, *pdate, *elev,
-	*bgain, *metho, *perc, *dark, *atmo, *lsatmet;
+    struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate,
+	*pdate, *elev, *bgain, *metho, *perc, *dark, *atmo, *lsatmet;
     char *inputname, *met, *outputname, *sensorname;
-    struct Flag *frad, *print_meta;
+    struct Flag *frad, *print_meta, *named;
 
     lsat_data lsat;
     char band_in[GNAME_MAX], band_out[GNAME_MAX];
@@ -64,7 +64,7 @@ int main(int argc, char *argv[])
     /* initialize module */
     module = G_define_module();
     module->description =
-	_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+.");
+	_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OT.");
     G_add_keyword(_("imagery"));
     G_add_keyword(_("radiometric conversion"));
     G_add_keyword(_("radiance"));
@@ -85,7 +85,8 @@ int main(int argc, char *argv[])
     output_prefix = G_define_option();
     output_prefix->key = "output_prefix";
     output_prefix->label = _("Prefix for output raster maps");
-    output_prefix->description = _("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
+    output_prefix->description =
+	_("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
     output_prefix->type = TYPE_STRING;
     output_prefix->required = YES;
 
@@ -99,11 +100,11 @@ int main(int argc, char *argv[])
     sensor->key = "sensor";
     sensor->type = TYPE_STRING;
     sensor->label = _("Spacecraft sensor");
-    sensor->description = _("Required only if 'metfile' not given");
-    sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7";
+    sensor->description = _("Required only if 'metfile' not given (recommended by sanity)");
+    sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,ot8";
     desc = NULL;
     G_asprintf(&desc,
-	        "mss1;%s;mss2;%s;mss3;%s;mss4;%s;mss5;%s;tm4;%s;tm5;%s;tm7;%s",
+	        "mss1;%s;mss2;%s;mss3;%s;mss4;%s;mss5;%s;tm4;%s;tm5;%s;tm7;%s;ot8;%s",
 	        _("Landsat-1 MSS"),
 	        _("Landsat-2 MSS"),
 	        _("Landsat-3 MSS"),
@@ -111,7 +112,8 @@ int main(int argc, char *argv[])
 	        _("Landsat-5 MSS"),
 	        _("Landsat-4 TM"),
 	        _("Landsat-5 TM"),
-	        _("Landsat-7 ETM+"));
+	        _("Landsat-7 ETM+"),
+		_("Landsat_8 OLI/TIRS"));
     sensor->descriptions = desc;
     sensor->required = NO;
     sensor->guisection = _("Metadata");
@@ -156,7 +158,7 @@ int main(int argc, char *argv[])
     bgain->key = "gain";
     bgain->type = TYPE_STRING;
     bgain->required = NO;
-    bgain->label = 	_("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
+    bgain->label = _("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
     bgain->description = _("Required only if 'metfile' not given");
     bgain->guisection = _("Settings");
 
@@ -164,7 +166,8 @@ int main(int argc, char *argv[])
     perc->key = "percent";
     perc->type = TYPE_DOUBLE;
     perc->required = NO;
-    perc->description = _("Percent of solar radiance in path radiance");
+    perc->label = _("Percent of solar radiance in path radiance");
+    perc->description = _("Required only if 'method' is any DOS");
     perc->answer = "0.01";
     perc->guisection = _("Settings");
 
@@ -172,8 +175,9 @@ int main(int argc, char *argv[])
     dark->key = "pixel";
     dark->type = TYPE_INTEGER;
     dark->required = NO;
-    dark->description =
+    dark->label =
 	_("Minimum pixels to consider digital number as dark object");
+    dark->description = _("Required only if 'method' is any DOS");
     dark->answer = "1000";
     dark->guisection = _("Settings");
 
@@ -181,7 +185,8 @@ int main(int argc, char *argv[])
     atmo->key = "rayleigh";
     atmo->type = TYPE_DOUBLE;
     atmo->required = NO;
-    atmo->description = _("Rayleigh atmosphere (diffuse sky irradiance)");	/* scattering coefficient? */
+    atmo->label = _("Rayleigh atmosphere (diffuse sky irradiance)");	/* scattering coefficient? */
+    atmo->description = _("Required only if 'method' is DOS3");
     atmo->answer = "0.0";
     atmo->guisection = _("Settings");
 
@@ -189,6 +194,7 @@ int main(int argc, char *argv[])
     lsatmet->key = "lsatmet";
     lsatmet->type = TYPE_STRING;
     lsatmet->required = NO;
+    lsatmet->multiple = YES;
     lsatmet->label = 	_("return value stored for a given metadata");
     lsatmet->description = _("Required only if 'metfile' and -p given");
     lsatmet->options = "number,creation,date,sun_elev,sensor,bands,sunaz,time";
@@ -209,7 +215,13 @@ int main(int argc, char *argv[])
     /* define the different flags */
     frad = G_define_flag();
     frad->key = 'r';
-    frad->description = _("Output at-sensor radiance for all bands");
+    frad->description =
+	_("Output at-sensor radiance instead of reflectance for all bands");
+
+    named = G_define_flag();
+    named->key = 'n';
+    named->description =
+	_("Input raster maps use as extension the number of the band instead the code");
 
     /* define the different flags */
     print_meta = G_define_flag();
@@ -228,7 +240,7 @@ int main(int argc, char *argv[])
     met = metfn->answer;
     inputname = input_prefix->answer;
     outputname = output_prefix->answer;
-    sensorname = sensor->answer ? sensor->answer: "";
+    sensorname = sensor->answer ? sensor->answer : "";
 
     overwrite = G_check_overwrite(argc, argv);
 
@@ -251,53 +263,55 @@ int main(int argc, char *argv[])
     }
 
     lsat.sun_elev = elev->answer == NULL ? 0. : atof(elev->answer);
+
     percent = atof(perc->answer);
     pixel = atoi(dark->answer);
     rayleigh = atof(atmo->answer);
 
     /* Data from metadata file */
-    /* Unnecessary because G_zero filled, but sanity */
+    /* Unnecessary because G_zero filled, but for sanity */
     lsat.flag = NOMETADATAFILE;
-    if (met != NULL)
-    {
+    if (met != NULL) {
         lsat.flag = METADATAFILE;
-        lsat_metadata( met, &lsat );
-	if(print_meta->answer) {
-		if (lsatmet->answer == NULL) {
-			G_fatal_error(_("Please use a metadata keyword with -p"));
-		}
-		if (strcmp(lsatmet->answer, "number") == 0) {
-			fprintf(stdout,"%d\n",lsat.number);
-		}
-		if (strcmp(lsatmet->answer, "creation") == 0) {
-			fprintf(stdout,"%s\n",lsat.creation);
-		}
-		if (strcmp(lsatmet->answer, "date") == 0) {
-			fprintf(stdout,"%s\n",lsat.date);
-		}
-		if (strcmp(lsatmet->answer, "sun_elev") == 0) {
-			fprintf(stdout,"%f\n",lsat.sun_elev);
-		}
-		if (strcmp(lsatmet->answer, "sensor") == 0) {
-			fprintf(stdout,"%s\n",lsat.sensor);
-		}
-		if (strcmp(lsatmet->answer, "bands") == 0) {
-			fprintf(stdout,"%d\n",lsat.bands);
-		}
-		if (strcmp(lsatmet->answer, "sunaz") == 0) {
-			fprintf(stdout,"%f\n",lsat.sunaz);
-		}
-		if (strcmp(lsatmet->answer, "time") == 0) {
-			fprintf(stdout,"%f\n",lsat.time);
-		}
-    		exit(EXIT_SUCCESS);
+        lsat_metadata(met, &lsat);
+	if (print_meta->answer) {
+	    if (lsatmet->answer == NULL) {
+		G_fatal_error(_("Please use a metadata keyword with -p"));
+	    }
+	    if (strcmp(lsatmet->answer, "number") == 0) {
+		fprintf(stdout,"%d\n",lsat.number);
+	    }
+	    if (strcmp(lsatmet->answer, "creation") == 0) {
+		fprintf(stdout,"%s\n",lsat.creation);
+	    }
+	    if (strcmp(lsatmet->answer, "date") == 0) {
+		fprintf(stdout,"%s\n",lsat.date);
+	    }
+	    if (strcmp(lsatmet->answer, "sun_elev") == 0) {
+		fprintf(stdout,"%f\n",lsat.sun_elev);
+	    }
+	    if (strcmp(lsatmet->answer, "sunaz") == 0) {
+		fprintf(stdout,"%f\n",lsat.sun_az);
+	    }
+	    if (strcmp(lsatmet->answer, "sensor") == 0) {
+		fprintf(stdout,"%s\n",lsat.sensor);
+	    }
+	    if (strcmp(lsatmet->answer, "bands") == 0) {
+		fprintf(stdout,"%d\n",lsat.bands);
+	    }
+	    if (strcmp(lsatmet->answer, "time") == 0) {
+		fprintf(stdout,"%f\n",lsat.time);
+	    }
+	    exit(EXIT_SUCCESS);
 	}
-        G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number, lsat.sensor);
+        G_debug(1, "lsat.number = %d, lsat.sensor = [%s]",
+	        lsat.number, lsat.sensor);
         
-        if (!lsat.sensor || lsat.number > 7 || lsat.number < 1)
+        if (!lsat.sensor || lsat.number > 8 || lsat.number < 1)
             G_fatal_error(_("Failed to identify satellite"));
 
-        G_debug(1, "Landsat-%d %s with data set in met file [%s]", lsat.number, lsat.sensor, met);
+        G_debug(1, "Landsat-%d %s with data set in met file [%s]",
+	        lsat.number, lsat.sensor, met);
 
         /* Overwrite solar elevation of metadata file */
         if (elev->answer != NULL)
@@ -305,7 +319,7 @@ int main(int argc, char *argv[])
     }
     /* Data from date and solar elevation */
     else if (adate->answer == NULL || elev->answer == NULL) {
-	G_fatal_error(_("Lacking '%s' or '%s' for this satellite"),
+	G_fatal_error(_("Missing '%s' and/or '%s' for this satellite"),
 		      adate->key, elev->key);
     }
     else {
@@ -316,6 +330,8 @@ int main(int argc, char *argv[])
             set_ETM(&lsat, bgain->answer);
         }
         /* Not need gain */
+	else if (strcmp(sensorname, "ot8") == 0)
+	    set_LDCM(&lsat);
         else if (strcmp(sensorname, "tm5") == 0)
             set_TM5(&lsat);
         else if (strcmp(sensorname, "tm4") == 0)
@@ -331,12 +347,13 @@ int main(int argc, char *argv[])
         else if (strcmp(sensorname, "mss1") == 0)
             set_MSS1(&lsat);
         else
-            G_fatal_error(_("Unknown satellite type (defined by '%s')"), sensor->key);
+            G_fatal_error(_("Unknown satellite type (defined by '%s')"),
+	                  sensor->key);
     }
 
-	/*****************************************
-	* ------------ PREPARATION --------------
-	*****************************************/
+    /*****************************************
+    * ------------ PREPARATION --------------
+    *****************************************/
     if (strcasecmp(metho->answer, "corrected") == 0)
 	method = CORRECTED;
     else if (strcasecmp(metho->answer, "dos1") == 0)
@@ -406,14 +423,14 @@ int main(int argc, char *argv[])
 	    }
 	    /* DN of dark object */
 	    for (j = lsat.band[i].qcalmin; j < 256; j++) {
-	      if (hist[j] >= (unsigned int) pixel) {
+		if (hist[j] >= (unsigned int)pixel) {
 		    dn_dark[i] = j;
 		    break;
 		}
 	    }
 	    /* Mode of DN */
 	    h_max = 0L;
-	    for (j = lsat.band[i].qcalmin; j < 241; j++) {	/* Exclude ptentially saturated < 240 */
+	    for (j = lsat.band[i].qcalmin; j < 241; j++) {	/* Exclude potentially saturated < 240 */
 		/* G_debug(5, "%d-%ld", j, hist[j]); */
 		if (hist[j] > h_max) {
 		    h_max = hist[j];
@@ -433,21 +450,24 @@ int main(int argc, char *argv[])
 	lsat_bandctes(&lsat, i, method, percent, dn_dark[i], rayleigh);
     }
 
+    /* unnecessary or necessary with more checking as acquisition date,... 
     if (strlen(lsat.creation) == 0)
 	G_fatal_error(_("Unknown production date (defined by '%s')"), pdate->key);
+    */
 
     if (G_verbose() > G_verbose_std()) {
-	fprintf(stderr, " LANDSAT: %d SENSOR: %s\n", lsat.number, lsat.sensor);
+	fprintf(stderr, " LANDSAT: %d SENSOR: %s\n", lsat.number,
+	        lsat.sensor);
 	fprintf(stderr, " ACQUISITION DATE %s [production date %s]\n",
 		lsat.date, lsat.creation);
-	fprintf(stderr, "   earth-sun distance    = %.8lf\n", lsat.dist_es);
-	fprintf(stderr, "   solar elevation angle = %.8lf\n", lsat.sun_elev);
-	fprintf(stderr, "   Method of calculus = %s\n",
+	fprintf(stderr, "   Earth-sun distance    = %.8lf\n", lsat.dist_es);
+	fprintf(stderr, "   Solar elevation angle = %.8lf\n", lsat.sun_elev);
+	fprintf(stderr, "   Atmospheric correction = %s\n",
 		(method == CORRECTED ? "CORRECTED"
 		 : (method == UNCORRECTED ? "UNCORRECTED" : metho->answer)));
 	if (method > DOS) {
 	    fprintf(stderr,
-		    "   percent of solar irradiance in path radiance = %.4lf\n",
+		    "   Percent of solar irradiance in path radiance = %.4lf\n",
 		    percent);
 	}
 	for (i = 0; i < lsat.bands; i++) {
@@ -461,12 +481,12 @@ int main(int argc, char *argv[])
 		    lsat.band[i].qcalmin, lsat.band[i].qcalmax);
 	    fprintf(stderr, "   calibration constants (L): %.3lf to %.3lf\n",
 		    lsat.band[i].lmin, lsat.band[i].lmax);
-	    fprintf(stderr, "   at-%s radiance = %.5lf * DN + %.5lf\n",
+	    fprintf(stderr, "   at-%s radiance = %.8lf * DN + %.3lf\n",
 		    (method > DOS ? "surface" : "sensor"), lsat.band[i].gain,
 		    lsat.band[i].bias);
 	    if (lsat.band[i].thermal) {
 		fprintf(stderr,
-			"   at-sensor temperature = %.3lf / log[(%.3lf / radiance) + 1.0]\n",
+			"   at-sensor temperature = %.5lf / log[(%.5lf / radiance) + 1.0]\n",
 			lsat.band[i].K2, lsat.band[i].K1);
 	    }
 	    else {
@@ -475,12 +495,12 @@ int main(int argc, char *argv[])
 			lsat.band[i].esun);
 		fprintf(stderr, "   at-%s reflectance = radiance / %.5lf\n",
 			(method > DOS ? "surface" : "sensor"),
-			lsat.band[i].K2);
+			lsat.band[i].K1);
 		if (method > DOS) {
 		    fprintf(stderr,
 			    "   the darkness DN with a least %d pixels is %d\n",
 			    pixel, dn_dark[i]);
-		    fprintf(stderr, "   the mode of DN is %d\n", dn_mode[i]);
+		    fprintf(stderr, "   the DN mode is %d\n", dn_mode[i]);
 		}
 	    }
 	}
@@ -488,13 +508,14 @@ int main(int argc, char *argv[])
 	fflush(stderr);
     }
 
-	/*****************************************
-	 * ------------ CALCULUS -----------------
-	 *****************************************/
+    /*****************************************
+     * ------------ CALCULUS -----------------
+     *****************************************/
 
     G_message(_("Calculating..."));
     for (i = 0; i < lsat.bands; i++) {
-	sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
+	sprintf(band_in, "%s%d", inputname,
+	        (named->answer ? lsat.band[i].number : lsat.band[i].code));
 	sprintf(band_out, "%s%d", outputname, lsat.band[i].code);
 	
 	if ((infd = Rast_open_old(band_in, "")) < 0)
@@ -531,10 +552,9 @@ int main(int argc, char *argv[])
 	ncols = Rast_window_cols();
 	/* ================================================================= */
 	G_important_message(_("Writing %s of <%s> to <%s>..."),
-		  (frad->answer ? _("radiance")
-		   : (lsat.band[i].
-		      thermal) ? _("temperature") : _("reflectance")),
-		  band_in, band_out);
+		            (frad->answer ? _("radiance")
+		             : (lsat.band[i].thermal) ? _("temperature")
+			        : _("reflectance")), band_in, band_out);
 	for (row = 0; row < nrows; row++) {
 	    G_percent(row, nrows, 2);
 
@@ -609,8 +629,9 @@ int main(int argc, char *argv[])
 	/* Append a string to the 'history' structure */
 	Rast_append_format_history(&history,
 				   " %s of Landsat-%d %s (method %s)",
-				   lsat.band[i].
-				   thermal ? "Temperature" : "Reflectance",
+				   (frad->answer ? "Radiance" : 
+				    (lsat.band[i].thermal ? "Temperature" : 
+				     "Reflectance")),
 				   lsat.number, lsat.sensor, metho->answer);
 	Rast_append_history(&history,
 			    "----------------------------------------------------------------");
@@ -644,7 +665,7 @@ int main(int argc, char *argv[])
 				       lsat.band[i].esun);
 	    Rast_append_format_history(&history,
 				       " Reflectance = Radiance divided by ..... %.5lf",
-				       lsat.band[i].K2);
+				       lsat.band[i].K1);
 	    if (method > DOS) {
 		Rast_append_history(&history, " ");
 		Rast_append_format_history(&history,
@@ -663,8 +684,10 @@ int main(int argc, char *argv[])
 
 	if (lsat.band[i].thermal)
 	    Rast_write_units(band_out, "Kelvin");
-	/* else units = ...? */
-	/* set raster timestamp from acq date? (see r.timestamp module) */
+	else if (frad->answer)
+	    Rast_write_units(band_out, "W/(m^2 sr µm)");
+	else
+	    Rast_write_units(band_out, "unitless");
     }
 
     exit(EXIT_SUCCESS);