i.landsat.toar: fixes for Landsat-8 metadata file support (author: E. Jorge Tizado); sync with G6 version

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

imagery/i.landsat.toar/earth_sun.c

@@ -1096,6 +1096,8 @@ double earth_sun(char *date)
     R4 = ln_calc_series(earth_radius_r4, RADIUS_R4, t);
     R5 = ln_calc_series(earth_radius_r5, RADIUS_R5, t);
 
-    return (R0 + R1 * t + R2 * t * t + R3 * t * t * t + R4 * t * t * t * t +
-	    R5 * t * t * t * t * t);
+    return (R0 +
+	    R1 * t +
+	    R2 * t * t +
+	    R3 * t * t * t + R4 * t * t * t * t + R5 * t * t * t * t * t);
 }

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

@@ -14,7 +14,7 @@ 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, ETM+ and OLI/TIRS. However,
-if the solar elevation is given the value of the metadata file are
+if the solar elevation is given the value of the metadata file is
 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.
@@ -60,26 +60,6 @@ where, <em>d</em> is the earth-sun distance in astronomical
 units, <em>e</em> is the solar elevation angle, and <em>Esun</em> is
 the mean solar exoatmospheric irradiance in W/(m&sup2; * &micro;m).
 
-<h2>Corrected at-sensor values (method=corrected)</h2>
-
-At-sensor reflectance values range from zero to one, whereas at-sensor
-radiance must be greater or equal to zero. However, since Lmin can be
-a negative number then the at-sensor values can also be negative. To
-avoid these possible negative values and set the minimum possible
-values at-sensor to zero, this method corrects the radiance to output
-a corrected at-sensor values using the equations (not for thermal
-bands):
-<ul>
-  <li> radiance = (uncorrected_radiance - Lmin) </li>
-  <li> reflectance = radiance / sun_radiance </li>
-</ul>
-
-<p>
-<b>Note</b>: Other possibility to avoid negative values is set to zero
-this values (radiance and/or reflectance), but this option is ease
-with uncorrected method
-and <em><a href="r.mapcalc.html">r.mapcalc</a></em>.
-
 <h2>Simplified at-surface values (method=dos[1-4])</h2>
 
 Atmospheric correction and reflectance calibration remove the path
@@ -204,10 +184,9 @@ i.landsat.toar input_prefix=203_30. output_prefix=_toar \
 <h2>SEE ALSO</h2>
 
 <em>
-<a href="i.aster.toar.html">i.aster.toar</a>,
-<a href="i.atcorr.html">i.atcorr</a>,
-<a href="r.mapcalc.html">r.mapcalc</a>,
-<a href="r.in.gdal.html">r.in.gdal</a>
+  <a href="i.atcorr.html">i.atcorr</a>,
+  <a href="r.mapcalc.html">r.mapcalc</a>,
+  <a href="r.in.gdal.html">r.in.gdal</a>
 </em>
 
 <h2>AUTHOR</h2>

imagery/i.landsat.toar/landsat.c

@@ -5,9 +5,9 @@
 
 #include "landsat.h"
 
-#define PI   M_PI
-#define R2D  180. / M_PI
-#define D2R  M_PI / 180.
+#define PI  M_PI
+#define R2D 180./M_PI
+#define D2R M_PI/180.
 
 /****************************************************************************
  * PURPOSE: Calibrated Digital Number to at-satellite Radiance
@@ -30,7 +30,7 @@ double lsat_rad2ref(double rad, band_data * band)
  *****************************************************************************/
 double lsat_rad2temp(double rad, band_data * band)
 {
-    return (double)(band->K2 / log((band->K1 / rad) + 1.0));
+    return (double)(band->K2 / log((band->K1 / rad) + 1.));
 }
 
 /****************************************************************************
@@ -43,15 +43,15 @@ double lsat_rad2temp(double rad, band_data * band)
  *         i : band number
  *    method : level of atmospheric correction
  *   percent : percent of solar irradiance in path radiance
- *       dos : digital number of dark object for DOS
+ *      dark : digital number of  object for DOS
   *****************************************************************************/
 
 #define abs(x)	(((x)>0)?(x):(-x))
 
 void lsat_bandctes(lsat_data * lsat, int i, char method,
-		   double percent, int dos, double rayleigh)
+		   double percent, int dark, double rayleigh)
 {
-    double pi_d2, sin_e, cos_v, rad_sun;
+    double pi_d2, sin_e, cos_v;
 
     /* TAUv  = at. transmittance surface-sensor */
     /* TAUz  = at. transmittance sun-surface    */
@@ -63,9 +63,9 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
     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
-    */
+	 * 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:
@@ -97,10 +97,11 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 	case DOS4:
 	    {
 		double Ro =
-		    (lsat->band[i].lmax - lsat->band[i].lmin) * (dos -
-								 lsat->band
-								 [i].qcalmin)
-		    / (lsat->band[i].qcalmax - lsat->band[i].qcalmin) +
+		    (lsat->band[i].lmax - lsat->band[i].lmin) * (dark -
+								 lsat->
+								 band[i].
+								 qcalmin) /
+		    (lsat->band[i].qcalmax - lsat->band[i].qcalmin) +
 		    lsat->band[i].lmin;
 		double Tv = 1.;
 		double Tz = 1.;
@@ -112,10 +113,10 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 		    Lp = Ro -
 			percent * TAUv * (lsat->band[i].esun * sin_e * TAUz +
 					  PI * Lp) / pi_d2;
-		    Tz = 1 - (4 * pi_d2 * Lp) / (lsat->band[i].esun * sin_e);
+		    Tz = 1. -
+			(4. * pi_d2 * Lp) / (lsat->band[i].esun * sin_e);
 		    Tv = exp(sin_e * log(Tz) / cos_v);
 		    /* G_message("TAUv = %.5f (%.5f), TAUz = %.5f (%.5f) and Edown = %.5f\n", TAUv, Tv, TAUz, Tz, PI * Lp ); */
-		    /* } while( abs(TAUv - Tv) > 0.0000001 || abs(TAUz - Tz) > 0.0000001); */
 		} while (TAUv != Tv && TAUz != Tz);
 		TAUz = (Tz < 1. ? Tz : 1.);
 		TAUv = (Tv < 1. ? Tv : 1.);
@@ -128,41 +129,43 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 	    Edown = 0.;
 	    break;
 	}
-	rad_sun = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
+	lsat->band[i].K2 = 0.;
+	lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;	/* rad_sun */
 	if (method > DOS)
 	    G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
 			      TAUv, TAUz, Edown);
-
-	lsat->band[i].K1 = rad_sun;
-	lsat->band[i].K2 = 0.;
     }
 
     /** 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));
+	 * 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);
 
     if (method == UNCORRECTED || lsat->band[i].thermal) {
 	/* L = G * (DN - Qmin) + Lmin
-	   -> bias = Lmin - G * Qmin    */
+	 *  -> bias = Lmin - G * Qmin    
+	 */
 	lsat->band[i].bias =
 	    (lsat->band[i].lmin - lsat->band[i].gain * lsat->band[i].qcalmin);
     }
-    else {
-	if (method == CORRECTED) {
-	    /* L = G * (DN - Qmin) + Lmin - Lmin
-	       -> bias = - G * Qmin */
-	    lsat->band[i].bias =
-		-(lsat->band[i].gain * lsat->band[i].qcalmin);
-	    /* Another possibility is cut when rad < 0 */
-	}
-	else if (method > DOS) {
-	    /* L = Lsat - Lpath =
-	       G * DNsat + B - (G * dark + B - p * rad_sun) =
-	       G * DNsat - G * dark + p * rad_sun
-	       -> bias = p * rad_sun - G * dark */
-	    lsat->band[i].bias = percent * rad_sun - lsat->band[i].gain * dos;
-	}
+    /*
+       else {
+       if (method == CORRECTED) {
+       // L = G * (DN - Qmin) + Lmin - Lmin
+       -> bias = - G * Qmin *
+       lsat->band[i].bias =
+       -(lsat->band[i].gain * lsat->band[i].qcalmin);
+       // Another possibility is cut when rad < 0 *
+       }
+     */
+    else if (method > DOS) {
+	/* L = Lsat - Lpath = G * DNsat + B - (G *  + B - p * rad_sun) 
+	 *   = G * DNsat - G *  + p * rad_sun
+	 *  -> bias = p * rad_sun - G 
+	 */
+	lsat->band[i].bias =
+	    percent * lsat->band[i].K1 - lsat->band[i].gain * dark;
     }
 }

imagery/i.landsat.toar/landsat.h

@@ -25,38 +25,38 @@
 
 typedef struct
 {
-    int number;			/* Band number                   */
-    int code;			/* Band code                     */
+    int number;			/* Band number */
+    int code;			/* Band code */
 
-    double wavemax, wavemin;	/* Wavelength in µm              */
+    double wavemax, wavemin;	/* Wavelength in µm */
 
-    double esun;		/* Mean solar irradiance         */
-    double lmax, lmin;		/* Spectral radiance             */
-    double qcalmax, qcalmin;	/* Quantized calibrated pixel    */
-
-    char thermal;		/* Flag to thermal band          */
-    double gain, bias;		/* Gain and Bias of sensor       */
-    double K1, K2;		/* Thermal calibration or
-				   Rad2Ref constants             */
+    double esun;		/* Mean solar irradiance */
+    double lmax, lmin;		/* Spectral radiance */
+    double qcalmax, qcalmin;	/* Quantized calibrated pixel */
 
+    char thermal;		/* Flag to thermal band */
+    double gain, bias;		/* Gain and Bias of sensor */
+    double K1, K2;		/* Thermal calibration 
+				 * or Rad2Ref constants */
 } band_data;
 
 typedef struct
 {
-    int flag;
-    unsigned char number;	/* Landsat number                */
+    int flag;			/* Line-data or file-data */
+    unsigned char number;	/* Landsat number */
+
+    char creation[11];		/* Image production date */
+    char date[11];		/* Image acquisition date */
+    double time;		/* Image acquisition time */
 
-    char creation[11];		/* Image production date         */
-    char date[11];		/* Image acquisition date        */
+    double dist_es;		/* Distance Earth-Sun */
+    double sun_elev;		/* Sun elevation */
+    double sun_az;		/* Sun azimuth */
 
-    double dist_es;		/* Distance Earth-Sun            */
-    double sun_elev;		/* Sun elevation                 */
-    double sun_az;		/* Sun Azimuth                   */
-    double time;		/* Image Acquisition Time        */
+    char sensor[10];		/* Sensor: MSS, TM, ETM+, OLI/TIRS */
+    int bands;			/* Total number of bands */
+    band_data band[MAX_BANDS];	/* Data for each band */
 
-    char sensor[10];		/* Type of sensor: MSS, TM, ETM+, OLI/TIRS */
-    int bands;			/* Total number of bands         */
-    band_data band[MAX_BANDS];	/* Data for each band            */
 } lsat_data;
 
 

imagery/i.landsat.toar/landsat_met.c

@@ -9,10 +9,10 @@
 #include "local_proto.h"
 #include "earth_sun.h"
 
-#define PI   M_PI
-#define D2R  M_PI / 180.
+#define PI  M_PI
+#define D2R M_PI/180.
 
-#define MAX_STR         127
+#define MAX_STR         256
 #define METADATA_SIZE   65535	/* MTL.txt file size  65535 bytes */
 #define TM5_MET_SIZE    28700	/* .met file size 28686 bytes */
 
@@ -26,6 +26,15 @@ static void chrncpy(char *dest, char src[], int n)
     dest[i] = '\0';
 }
 
+inline void date_replace_slash(char *str)
+{
+    while (*str) {
+	if (*str == '/')
+	    *str = '-';
+	str++;
+    }
+}
+
 void get_metformat(const char metadata[], char *key, char value[])
 {
     int i = 0;
@@ -33,14 +42,16 @@ void get_metformat(const char metadata[], char *key, char value[])
     char *ptr;
 
     if (ptrmet != NULL) {
-        ptr = strstr(ptrmet, " VALUE ");
-        if (ptr != NULL) {
-            while (*ptr++ != '=') ;
-            while (*ptr <= ' ' || *ptr == '\"')
+	ptr = strstr(ptrmet, " VALUE ");
+	if (ptr != NULL) {
+	    while (*ptr++ != '=') ;
+	    while (*ptr <= ' ')
+		ptr++;
+	    if (*ptr == '\"')
 		ptr++;
-            while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
+	    while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
 		value[i++] = *ptr++;
-        }
+	}
     }
     value[i] = '\0';
 }
@@ -51,12 +62,12 @@ void get_mtlformat(const char metadata[], char *key, char value[])
     char *ptr = strstr(metadata, key);
 
     if (ptr != NULL) {
-        while (*ptr++ != '=') ;
-
-        while (*ptr <= ' ' || *ptr == '\"')
+	while (*ptr++ != '=') ;
+	while (*ptr <= ' ')
 	    ptr++;
-
-        while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
+	if (*ptr == '\"')
+	    ptr++;
+	while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
 	    value[i++] = *ptr++;
     }
     value[i] = '\0';
@@ -70,83 +81,90 @@ void get_mtlformat(const char metadata[], char *key, char value[])
 void lsat_metadata(char *metafile, lsat_data * lsat)
 {
     FILE *f;
-    char metadata[METADATA_SIZE];
+    char mtldata[METADATA_SIZE];
     char key[MAX_STR], value[MAX_STR];
-    void (*get_metadata) (const char [], char *, char []);
-    int i, j, ver_meta;
+    void (*get_mtldata) (const char[], char *, char[]);
+    int i, j, ver_mtl;
 
+    /* store metadata in ram */
     if ((f = fopen(metafile, "r")) == NULL)
 	G_fatal_error(_("Metadata file <%s> not found"), metafile);
+    i = fread(mtldata, METADATA_SIZE, 1, f);
+    (void)fclose(f);
 
-    i = fread(metadata, METADATA_SIZE, 1, f);
+    /* set version of the metadata file */
+    get_mtldata =
+	(strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
+    ver_mtl = (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_metadata(metadata, "SPACECRAFT_ID", value);
+    get_mtldata(mtldata, "SPACECRAFT_ID", value);
     if (value[0] == '\0') {
-        get_metadata(metadata, "PLATFORMSHORTNAME", value);
+	get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
     }
-    
     i = 0;
     while (value[i] && (value[i] < '0' || value[i] > '9'))
 	i++;
     lsat->number = atoi(value + i);
 
-    get_metadata(metadata, "SENSOR_ID", value);
+    get_mtldata(mtldata, "SENSOR_ID", value);
     if (value[0] == '\0') {
-        get_metadata(metadata, "SENSORSHORTNAME", value);
+	get_mtldata(mtldata, "SENSORSHORTNAME", value);
     }
     chrncpy(lsat->sensor, value, 8);
 
-    get_metadata(metadata, "DATE_ACQUIRED", value);
+    get_mtldata(mtldata, "DATE_ACQUIRED", value);
     if (value[0] == '\0') {
-        get_metadata(metadata, "ACQUISITION_DATE", value);
-        if (value[0] == '\0') {
-            get_metadata(metadata, "CALENDARDATE", value);
-        }
+	get_mtldata(mtldata, "ACQUISITION_DATE", value);
+	if (value[0] == '\0') {
+	    get_mtldata(mtldata, "CALENDARDATE", value);
+	}
+    }
+    if (value[0] != '\0') {
+	date_replace_slash(value);
+	chrncpy(lsat->date, value, 10);
     }
-    chrncpy(lsat->date, value, 10);
+    else
+	G_warning("Using adquisition date from the command line 'date'");
 
-    get_metadata(metadata, "FILE_DATE", value);
+    get_mtldata(mtldata, "FILE_DATE", value);
     if (value[0] == '\0') {
-        get_metadata(metadata, "CREATION_TIME", value);
-        if (value[0] == '\0') {
-            get_metadata(metadata, "PRODUCTIONDATETIME", value);
-        }
+	get_mtldata(mtldata, "CREATION_TIME", value);
+	if (value[0] == '\0') {
+	    get_mtldata(mtldata, "PRODUCTIONDATETIME", value);
+	}
     }
-    if (value[0] != '\0')
+    if (value[0] != '\0') {
+	date_replace_slash(value);
 	chrncpy(lsat->creation, value, 10);
+    }
     else
 	G_warning
 	    ("Using production date from the command line 'product_date'");
 
-    get_metadata(metadata, "SUN_AZIMUTH", value);
+    get_mtldata(mtldata, "SUN_AZIMUTH", value);
     lsat->sun_az = atof(value);
     if (lsat->sun_az == 0.)
-        G_warning("Sun azimuth is %f", lsat->sun_az);
+	G_warning("Sun azimuth is %f", lsat->sun_az);
 
-    get_metadata(metadata, "SUN_ELEVATION", value);
+    get_mtldata(mtldata, "SUN_ELEVATION", value);
     if (value[0] == '\0') {
-        get_metadata(metadata, "SolarElevation", value);
+	get_mtldata(mtldata, "SolarElevation", value);
     }
     lsat->sun_elev = atof(value);
     if (lsat->sun_elev == 0.)
-        G_warning("Sun elevation is %f", lsat->sun_elev);
+	G_warning("Sun elevation is %f", lsat->sun_elev);
 
-    get_metadata(metadata, "SCENE_CENTER_TIME", value);
+    get_mtldata(mtldata, "SCENE_CENTER_TIME", value);
     if (value[0] == '\0') {
-        get_metadata(metadata, "SCENE_CENTER_SCAN_TIME", value);
+	get_mtldata(mtldata, "SCENE_CENTER_SCAN_TIME", value);
+    }
+    if (value[0] != '\0') {
+	value[strlen(value) - 1] = '\0';	/* Remove trailing 'z' */
+	G_llres_scan(value, &lsat->time);	/* Cast from hh:mm:ss into hh.hhh */
     }
-    /* 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.)
-        G_warning("Time is %f", lsat->time);
+	G_warning("Scene time is %f", lsat->time);
 
     /* Fill the data with the basic/default sensor parameters */
     switch (lsat->number) {
@@ -175,13 +193,13 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
 	{
 	    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);
+		sprintf(key, fmt_gain[ver_mtl], i, 1);
 		if (i != 6)
-		    key[ver_meta == 0 ? 10 : 12] = '\0';
-		get_metadata(metadata, key, value + j++);
+		    key[ver_mtl == 0 ? 10 : 12] = '\0';
+		get_mtldata(mtldata, key, value + j++);
 		if (i == 6) {
-		    sprintf(key, fmt_gain[ver_meta], i, 2);
-		    get_metadata(metadata, key, value + j++);
+		    sprintf(key, fmt_gain[ver_mtl], i, 2);
+		    get_mtldata(mtldata, key, value + j++);
 		}
 	    }
 	    value[j] = '\0';
@@ -190,9 +208,8 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
 	    break;
 	}
     case 8:
-	set_LDCM(lsat);
+	set_OLI(lsat);
 	break;
-
     default:
 	G_warning("Unable to recognize satellite platform [%d]",
 		  lsat->number);
@@ -200,155 +217,168 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
     }
 
     /* Update the information from metadata file, if infile */
-    if (get_metadata == get_mtlformat) {
-	if (ver_meta == 0) {	/* now MTLold.txt */
+    if (get_mtldata == get_mtlformat) {
+	if (ver_mtl == 0) {
 	    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);
+		get_mtldata(mtldata, key, value);
 		lsat->band[i].lmax = atof(value);
 		sprintf(key, "LMIN_BAND%d", lsat->band[i].code);
-		get_metadata(metadata, key, value);
+		get_mtldata(mtldata, key, value);
 		lsat->band[i].lmin = atof(value);
 		sprintf(key, "QCALMAX_BAND%d", lsat->band[i].code);
-		get_metadata(metadata, key, value);
+		get_mtldata(mtldata, key, value);
 		lsat->band[i].qcalmax = atof(value);
 		sprintf(key, "QCALMIN_BAND%d", lsat->band[i].code);
-		get_metadata(metadata, key, value);
+		get_mtldata(mtldata, key, value);
 		lsat->band[i].qcalmin = atof(value);
 	    }
 	}
-	else {			/* now MTL.txt */
-
+	else {
 	    G_verbose_message("Metada file is MTL file: new format");
-	    if (strstr(metadata, "RADIANCE_MAXIMUM_BAND") != NULL) {
+	    /* Other possible values in the metadata file */
+	    get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value);	/* used after in calculus after */
+	    if (value[0] != '\0')
+		lsat->dist_es = atof(value);
+	    /* ----- */
+	    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_refmu[] =
+		{ "REFLECTANCE_MULTIPLICATIVE_FACTOR_BAND%d",
+    "REFLECTANCE_MULT_BAND_%d" };
+	    char *fmt_refad[] =
+		{ "REFLECTANCE_ADDITIVE_FACTOR_BAND%d",
+    "REFLECTANCE_ADD_BAND_%d" };
+	    /* --NASA-- LDCM sample file: mode = 0; LDCM-DFCB-004.pdf file: mode = 1 */
+	    int mode =
+		(strstr(mtldata, "RADIANCE_MULTIPLICATIVE_FACTOR_BAND") !=
+		 NULL) ? 0 : 1;
+	    /* ----- */
+	    if (strstr(mtldata, "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) {
+		    if (lsat->number == 7 && lsat->band[i].thermal) {
 			sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
 				lsat->band[i].code - 60);
-			get_metadata(metadata, key, value);
+			get_mtldata(mtldata, 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);
+			get_mtldata(mtldata, 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);
+			get_mtldata(mtldata, 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);
+			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmin = atof(value);
 		    }
 		    else {
 			sprintf(key, "RADIANCE_MAXIMUM_BAND_%d",
 				lsat->band[i].code);
-			get_metadata(metadata, key, value);
+			get_mtldata(mtldata, key, value);
 			lsat->band[i].lmax = atof(value);
 			sprintf(key, "RADIANCE_MINIMUM_BAND_%d",
 				lsat->band[i].code);
-			get_metadata(metadata, key, value);
+			get_mtldata(mtldata, key, value);
 			lsat->band[i].lmin = atof(value);
 			sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d",
 				lsat->band[i].code);
-			get_metadata(metadata, key, value);
+			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmax = atof(value);
 			sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d",
 				lsat->band[i].code);
-			get_metadata(metadata, key, value);
+			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmin = atof(value);
 		    }
+		    /* other possible values of each band */
+		    if (lsat->band[i].thermal) {
+			sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+			get_mtldata(mtldata, key, value);
+			lsat->band[i].K1 = atof(value);
+			sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+			get_mtldata(mtldata, key, value);
+			lsat->band[i].K2 = 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);
+		    get_mtldata(mtldata, key, value);
 		    lsat->band[i].gain = atof(value);
 		    sprintf(key, fmt_radad[mode], lsat->band[i].code);
-		    get_metadata(metadata, key, value);
+		    get_mtldata(mtldata, 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.;
+		    /* reversing to calculate the values of Lmin and Lmax ... */
 		    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;
-		    /* ----- */
+		    /* ... qcalmax and qcalmin loaded in previous sensor_ function */
 		    if (lsat->number == 8) {
 			if (lsat->band[i].thermal) {
 			    sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
-			    get_metadata(metadata, key, value);
+			    get_mtldata(mtldata, key, value);
 			    lsat->band[i].K1 = atof(value);
 			    sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
-			    get_metadata(metadata, key, value);
+			    get_mtldata(mtldata, 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_refmu[mode], lsat->band[i].code);
+			    get_mtldata(mtldata, key, value);
+			    lsat->band[i].K1 = atof(value);
 			    sprintf(key, fmt_refad[mode], lsat->band[i].code);
-			    get_metadata(metadata, key, value);
+			    get_mtldata(mtldata, key, value);
 			    lsat->band[i].K2 = atof(value);
+			    /* ESUN evaluates from metadafa file */
 			    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);
+				lsat->band[i].K2;
+			    /*
+			       double esun1 = (double) (PI * lsat->dist_es * lsat->dist_es * lsat->band[i].bias) / lsat->band[i].K2;
+			       double esun2 = (double) (PI * lsat->dist_es * lsat->dist_es * lsat->band[i].gain) / lsat->band[i].K1;
+			       lsat->band[i].esun = (esun1 + esun2) / 2.;
+			     */
 			}
 		    }
 		}
-	    }
-	    /* Other possible values in file */
-	    get_metadata(metadata, "EARTH_SUN_DISTANCE", value);
-	    if (value[0] != '\0') {
-		lsat->dist_es = atof(value);
-		G_verbose_message
+		G_warning
 		    ("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
 	    }
 	}
     }
     else {
-	G_verbose_message("Metada file is MET file");
+	G_verbose_message("Metadata 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);
+	    get_mtldata(mtldata, 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);
+	    get_mtldata(mtldata, key, value);
 	    if (value[0] == '\0') {
 		G_warning(key);
 		continue;
@@ -366,6 +396,5 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
 	}
     }
 
-    fclose(f);
     return;
 }

imagery/i.landsat.toar/landsat_set.c

@@ -68,8 +68,10 @@ void sensor_ETM(lsat_data * lsat)
     /* 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, 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 };
+    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 +89,7 @@ void sensor_ETM(lsat_data * lsat)
     return;
 }
 
-void sensor_LDCM(lsat_data * lsat)
+void sensor_OLI(lsat_data * lsat)
 {
     int i;
 
@@ -95,8 +97,12 @@ void sensor_LDCM(lsat_data * lsat)
      * cirrus, thermal infrared (TIR) 1, TIR 2 */
     int band[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
     int code[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
-    double wmin[] = { 0.433, 0.450, 0.525, 0.630, 0.845, 1.560, 2.100, 0.500, 1.360, 10.3, 11.5 };
-    double wmax[] = { 0.453, 0.515, 0.600, 0.680, 0.885, 1.660, 2.300, 0.680, 1.390, 11.3, 12.5 };
+    double wmin[] =
+	{ 0.433, 0.450, 0.525, 0.630, 0.845, 1.560, 2.100, 0.500, 1.360, 10.3,
+11.5 };
+    double wmax[] =
+	{ 0.453, 0.515, 0.600, 0.680, 0.885, 1.660, 2.300, 0.680, 1.390, 11.3,
+12.5 };
     /* 30, 30, 30, 30, 30, 30, 30, 15, 30, 100, 100 */
 
     strcpy(lsat->sensor, "OLI/TIRS");
@@ -107,7 +113,7 @@ void sensor_LDCM(lsat_data * lsat)
 	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].qcalmax = 65535.;
 	lsat->band[i].qcalmin = 1.;
 	lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
     }
@@ -444,7 +450,8 @@ void set_TM5(lsat_data * lsat)
 
     lmax = Lmax[i];
     lmin = Lmin[i];
-    if (i == 2) {		/* in Chander, Markham and Barsi 2007 */
+    /* in Chander, Markham and Barsi 2007 */
+    if (i == 2) {
 	julian = julian_char(lsat->date);	/* Yes, here acquisition date */
 	if (julian >= julian_char("1992-01-01")) {
 	    lmax[0] = 193.0;
@@ -516,10 +523,7 @@ void set_ETM(lsat_data * lsat, char gain[])
     /*  Thermal band calibration constants: K1 = 666.09   K2 = 1282.71 */
 
     julian = julian_char(lsat->creation);
-    if (julian < julian_char("2000-07-01"))
-	k = 0;
-    else
-	k = 1;
+    k = ((julian < julian_char("2000-07-01")) ? 0 : 1);
 
     lsat->number = 7;
     sensor_ETM(lsat);
@@ -552,29 +556,32 @@ void set_ETM(lsat_data * lsat, char gain[])
  * PURPOSE:     Store values of Landsat-8 OLI/TIRS
  *              February 14, 2013
  *****************************************************************************/
-void set_LDCM(lsat_data * lsat)
+void set_OLI(lsat_data * lsat)
 {
     int i, j;
     double *lmax, *lmin;
 
     /* Spectral radiances at detector */
-    
-    /* uncorrected values */ 
+    /* estimates */
     double Lmax[][11] = {
-	{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
+	{755.8, 770.7, 705.7, 597.7, 362.7, 91.4, 29.7, 673.3, 149.0, 22.0,
+	 22.0}
     };
     double Lmin[][11] = {
-	{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
+	{-62.4, -63.6, -58.3, -49.4, -30.0, -7.5, -2.5, -55.6, -12.3, 0.1, 0.1}
     };
+
     /* Solar exoatmospheric spectral irradiances */
     /* estimates */
-    double esun[] = { 2062., 21931., 1990., 1688., 1037., 268.6, 94.6, 1892., 399.0, 0., 0. };
+    double esun[] =
+	{ 2026.8, 2066.8, 1892.5, 1602.8, 972.6, 245.0, 79.7, 1805.5, 399.7,
+0., 0. };
 
     lmax = Lmax[0];
     lmin = Lmin[0];
 
     lsat->number = 8;
-    sensor_LDCM(lsat);
+    sensor_OLI(lsat);
 
     lsat->dist_es = earth_sun(lsat->date);
 
@@ -584,8 +591,9 @@ void set_LDCM(lsat_data * lsat)
 	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;
+	    lsat->band[i].K1 = (lsat->band[i].number == 10 ? 774.89 : 480.89);
+	    lsat->band[i].K2 =
+		(lsat->band[i].number == 10 ? 1321.08 : 1201.14);
 	}
     }
     G_debug(1, "Landsat-8 OLI/TIRS");

imagery/i.landsat.toar/local_proto.h

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

imagery/i.landsat.toar/main.c

@@ -3,7 +3,7 @@
  *
  * MODULE:       i.landsat.toar
  *
- * AUTHOR(S):    E. Jorge Tizado - ej.tizado@unileon.es
+ * AUTHOR(S):    E. Jorge Tizado - ej.tizado at unileon.es
  *               Hamish Bowman (small grassification cleanups)
  *               Yann Chemin (v7 + L5TM _MTL.txt support) [removed after update]
  *               Adopted for GRASS 7 by Martin Landa <landa.martin gmail.com>
@@ -27,6 +27,8 @@
 
 #include "local_proto.h"
 
+#define QCALMAX   65536		/* L1-7=256 but L8=65536 */
+
 int main(int argc, char *argv[])
 {
     struct History history;
@@ -49,14 +51,14 @@ int main(int argc, char *argv[])
 
     lsat_data lsat;
     char band_in[GNAME_MAX], band_out[GNAME_MAX];
-    int i, j, q, method, pixel, dn_dark[MAX_BANDS], dn_mode[MAX_BANDS];
+    int i, j, q, method, pixel, dn_dark[MAX_BANDS], dn_mode[MAX_BANDS], dn_sat;
     int overwrite;
     double qcal, rad, ref, percent, ref_mode, rayleigh;
+    unsigned long hist[QCALMAX], h_max;
 
     struct Colors colors;
     struct FPRange range;
     double min, max;
-    unsigned long hist[256], h_max;
 
     /* initialize GIS environment */
     G_gisinit(argv[0]);
@@ -64,7 +66,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+/OT.");
+	_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OLI");
     G_add_keyword(_("imagery"));
     G_add_keyword(_("radiometric conversion"));
     G_add_keyword(_("radiance"));
@@ -100,29 +102,29 @@ 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 (recommended for sanity)");
-    sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,ot8";
+    sensor->description = _("Required only if 'metfile' not given (recommended for sanity)");
+    sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,oli8";
     desc = NULL;
     G_asprintf(&desc,
-	       "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"),
-	       _("Landsat-4 MSS"),
-	       _("Landsat-5 MSS"),
-	       _("Landsat-4 TM"),
-	       _("Landsat-5 TM"),
-	       _("Landsat-7 ETM+"), _("Landsat_8 OLI/TIRS"));
+	        "mss1;%s;mss2;%s;mss3;%s;mss4;%s;mss5;%s;tm4;%s;tm5;%s;tm7;%s;oli8;%s",
+	        _("Landsat-1 MSS"),
+	        _("Landsat-2 MSS"),
+	        _("Landsat-3 MSS"),
+	        _("Landsat-4 MSS"),
+	        _("Landsat-5 MSS"),
+	        _("Landsat-4 TM"),
+	        _("Landsat-5 TM"),
+	        _("Landsat-7 ETM+"),
+	        _("Landsat_8 OLI/TIRS"));
     sensor->descriptions = desc;
-    sensor->required = NO;
+    sensor->required = NO;	/* perhaps YES for clarity */
     sensor->guisection = _("Metadata");
 
     metho = G_define_option();
     metho->key = "method";
     metho->type = TYPE_STRING;
     metho->required = NO;
-    metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
+    metho->options = "uncorrected,dos1,dos2,dos2b,dos3,dos4";
     metho->label = _("Atmospheric correction method");
     metho->description = _("Atmospheric correction method");
     metho->answer = "uncorrected";
@@ -232,10 +234,10 @@ int main(int argc, char *argv[])
 	exit(EXIT_FAILURE);
 
 
-    /*****************************************
-     * ---------- START --------------------
-     * Stores options and flag to variables
-     *****************************************/
+  /*****************************************
+  * ---------- START --------------------
+  * Stores options and flags to variables
+  *****************************************/
     met = metfn->answer;
     inputname = input_prefix->answer;
     outputname = output_prefix->answer;
@@ -267,9 +269,10 @@ int main(int argc, char *argv[])
     pixel = atoi(dark->answer);
     rayleigh = atof(atmo->answer);
 
-    /* Data from metadata file */
-    /* Unnecessary because G_zero filled, but for sanity */
-    lsat.flag = NOMETADATAFILE;
+    /*
+     * Data from metadata file
+     */
+    lsat.flag = NOMETADATAFILE;	/* Unnecessary because G_zero filled, but for sanity */
     if (met != NULL) {
 	lsat.flag = METADATAFILE;
 	lsat_metadata(met, &lsat);
@@ -316,28 +319,29 @@ int main(int argc, char *argv[])
 	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]",
+	G_debug(1, "Landsat-%d %s with data set in metadata file [%s]",
 		lsat.number, lsat.sensor, met);
 
-	/* Overwrite solar elevation of metadata file */
-	if (elev->answer != NULL)
+	if (elev->answer != NULL) {
 	    lsat.sun_elev = atof(elev->answer);
+	    G_warning("Overwriting solar elevation of metadata file");
+	}
     }
-    /* Data from date and solar elevation */
+    /*
+     * Data from command line
+     */
     else if (adate->answer == NULL || elev->answer == NULL) {
-	G_fatal_error(_("Missing '%s' and/or '%s' for this satellite"),
+	G_fatal_error(_("Lacking '%s' and/or '%s' for this satellite"),
 		      adate->key, elev->key);
     }
     else {
-	/* Need gain */
 	if (strcmp(sensorname, "tm7") == 0) {
 	    if (bgain->answer == NULL || strlen(bgain->answer) != 9)
-		G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) data"));
+		G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) characters"));
 	    set_ETM(&lsat, bgain->answer);
 	}
-	/* Not need gain */
-	else if (strcmp(sensorname, "ot8") == 0)
-	    set_LDCM(&lsat);
+	else if (strcmp(sensorname, "oli8") == 0)
+	    set_OLI(&lsat);
 	else if (strcmp(sensorname, "tm5") == 0)
 	    set_TM5(&lsat);
 	else if (strcmp(sensorname, "tm4") == 0)
@@ -354,13 +358,13 @@ int main(int argc, char *argv[])
 	    set_MSS1(&lsat);
 	else
 	    G_fatal_error(_("Unknown satellite type (defined by '%s')"),
-			  sensor->key);
+			  sensorname);
     }
 
-    /*****************************************
-    * ------------ PREPARATION --------------
-    *****************************************/
-    if (strcasecmp(metho->answer, "corrected") == 0)
+	/*****************************************
+	* ------------ PREPARATION --------------
+	*****************************************/
+    if (strcasecmp(metho->answer, "corrected") == 0)	/* deleted 2013 */
 	method = CORRECTED;
     else if (strcasecmp(metho->answer, "dos1") == 0)
 	method = DOS1;
@@ -376,21 +380,21 @@ int main(int argc, char *argv[])
 	method = UNCORRECTED;
 
     /*
-       if (metho->answer[3] == 'r')            method = CORRECTED;
-       else if (metho->answer[3] == '1')       method = DOS1;
-       else if (metho->answer[3] == '2')       method = (metho->answer[4] == '\0') ? DOS2 : DOS2b;
-       else if (metho->answer[3] == '3')       method = DOS3;
-       else if (metho->answer[3] == '4')       method = DOS4;
+       if (metho->answer[3] == '2') method = (metho->answer[4] == '\0') ? DOS2 : DOS2b;
+       else if (metho->answer[3] == '1') method = DOS1;
+       else if (metho->answer[3] == '3') method = DOS3;
+       else if (metho->answer[3] == '4') method = DOS4;
        else method = UNCORRECTED;
      */
 
     for (i = 0; i < lsat.bands; i++) {
 	dn_mode[i] = 0;
 	dn_dark[i] = (int)lsat.band[i].qcalmin;
-	/* Calculate dark pixel */
+	dn_sat = (int)(0.90 * lsat.band[i].qcalmax);
+	/* Begin: calculate dark pixel */
 	if (method > DOS && !lsat.band[i].thermal) {
-	    for (j = 0; j < 256; j++)
-		hist[j] = 0L;
+	    for (q = 0; q <= lsat.band[i].qcalmax; q++)
+		hist[q] = 0L;
 
 	    sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
 	    if ((infd = Rast_open_old(band_in, "")) < 0)
@@ -399,6 +403,8 @@ int main(int argc, char *argv[])
 	    G_set_window(&cellhd);
 
 	    in_data_type = Rast_get_map_type(infd);
+	    if (in_data_type < 0)
+		G_fatal_error(_("Unable to read data type of raster map <%s>"), band_in);
 	    inrast = Rast_allocate_buf(in_data_type);
 
 	    nrows = Rast_window_rows();
@@ -421,56 +427,62 @@ int main(int argc, char *argv[])
 			ptr = (void *)((DCELL *) inrast + col);
 			q = (int)*((DCELL *) ptr);
 			break;
+		   default:
+		        ptr = NULL;
+		        q = -1.;
+			break;
 		    }
 		    if (!Rast_is_null_value(ptr, in_data_type) &&
-			q >= lsat.band[i].qcalmin && q < 256)
+			q >= lsat.band[i].qcalmin &&
+			q <= lsat.band[i].qcalmax)
 			hist[q]++;
 		}
 	    }
 	    /* DN of dark object */
-	    for (j = lsat.band[i].qcalmin; j < 256; j++) {
+	    for (j = lsat.band[i].qcalmin; j <= lsat.band[i].qcalmax; j++) {
 		if (hist[j] >= (unsigned int)pixel) {
 		    dn_dark[i] = j;
 		    break;
 		}
 	    }
-	    /* Mode of DN */
+	    /* Mode of DN (exclude potentially saturated) */
 	    h_max = 0L;
-	    for (j = lsat.band[i].qcalmin; j < 241; j++) {	/* Exclude potentially saturated < 240 */
+	    for (j = lsat.band[i].qcalmin; j < dn_sat; j++) {
 		/* G_debug(5, "%d-%ld", j, hist[j]); */
 		if (hist[j] > h_max) {
 		    h_max = hist[j];
 		    dn_mode[i] = j;
 		}
 	    }
-	    G_verbose_message("... DN = %.2d [%lu] : mode %.2d [%lu] %s",
-			      dn_dark[i], hist[dn_dark[i]],
-			      dn_mode[i], hist[dn_mode[i]],
-			      hist[255] >
-			      hist[dn_mode[i]] ? ", excluding DN > 241" : "");
+	    G_verbose_message
+		("... DN = %.2d [%lu] : mode %.2d [%lu], excluding DN > %d",
+		 dn_dark[i], hist[dn_dark[i]], dn_mode[i], hist[dn_mode[i]],
+		 dn_sat);
 
 	    G_free(inrast);
 	    Rast_close(infd);
 	}
+	/* End: calculate dark pixel */
+
 	/* Calculate transformation constants */
 	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);
+    /*
+     * 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,
+	fprintf(stderr, "\n 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, "   Atmospheric correction = %s\n",
-		(method == CORRECTED ? "CORRECTED"
-		 : (method == UNCORRECTED ? "UNCORRECTED" : metho->answer)));
+	fprintf(stderr, "   Atmospheric correction: %s\n",
+		(method == UNCORRECTED ? "UNCORRECTED" : metho->answer));
 	if (method > DOS) {
 	    fprintf(stderr,
 		    "   Percent of solar irradiance in path radiance = %.4lf\n",
@@ -478,16 +490,15 @@ int main(int argc, char *argv[])
 	}
 	for (i = 0; i < lsat.bands; i++) {
 	    fprintf(stderr, "-------------------\n");
-	    fprintf(stderr, " BAND %d %s (code %d)\n",
-		    lsat.band[i].number,
+	    fprintf(stderr, " BAND %d %s(code %d)\n", lsat.band[i].number,
 		    (lsat.band[i].thermal ? "thermal " : ""),
 		    lsat.band[i].code);
 	    fprintf(stderr,
 		    "   calibrated digital number (DN): %.1lf to %.1lf\n",
 		    lsat.band[i].qcalmin, lsat.band[i].qcalmax);
-	    fprintf(stderr, "   calibration constants (L): %.3lf to %.3lf\n",
+	    fprintf(stderr, "   calibration constants (L): %.5lf to %.5lf\n",
 		    lsat.band[i].lmin, lsat.band[i].lmax);
-	    fprintf(stderr, "   at-%s radiance = %.8lf * DN + %.3lf\n",
+	    fprintf(stderr, "   at-%s radiance = %.8lf * DN + %.5lf\n",
 		    (method > DOS ? "surface" : "sensor"), lsat.band[i].gain,
 		    lsat.band[i].bias);
 	    if (lsat.band[i].thermal) {
@@ -497,7 +508,7 @@ int main(int argc, char *argv[])
 	    }
 	    else {
 		fprintf(stderr,
-			"   mean solar exoatmospheric irradiance (ESUN): %.3lf\n",
+			"   mean solar exoatmospheric irradiance (ESUN): %.5lf\n",
 			lsat.band[i].esun);
 		fprintf(stderr, "   at-%s reflectance = radiance / %.5lf\n",
 			(method > DOS ? "surface" : "sensor"),
@@ -515,8 +526,8 @@ int main(int argc, char *argv[])
     }
 
     /*****************************************
-     * ------------ CALCULUS -----------------
-     *****************************************/
+    * ------------ CALCULUS -----------------
+    *****************************************/
 
     G_message(_("Calculating..."));
     for (i = 0; i < lsat.bands; i++) {
@@ -539,6 +550,8 @@ int main(int argc, char *argv[])
 	}
 
 	in_data_type = Rast_get_map_type(infd);
+	if (in_data_type < 0)
+	    G_fatal_error(_("Unable to read data type of raster map <%s>"), band_in);
 	Rast_get_cellhd(band_in, "", &cellhd);
 
 	/* set same size as original band raster */
@@ -551,17 +564,19 @@ int main(int argc, char *argv[])
 	if ((outfd = Rast_open_new(band_out, DCELL_TYPE)) < 0)
 	    G_fatal_error(_("Unable to create raster map <%s>"), band_out);
 
-	/* Allocate input and output buffer */
+	/* allocate input and output buffer */
 	inrast = Rast_allocate_buf(in_data_type);
 	outrast = Rast_allocate_buf(DCELL_TYPE);
 
 	nrows = Rast_window_rows();
 	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);
 
@@ -580,6 +595,10 @@ int main(int argc, char *argv[])
 		    ptr = (void *)((DCELL *) inrast + col);
 		    qcal = (double)((DCELL *) inrast)[col];
 		    break;
+		default:
+		    ptr = NULL;
+		    qcal = -1.;
+		    break;
 		}
 		if (Rast_is_null_value(ptr, in_data_type) ||
 		    qcal < lsat.band[i].qcalmin) {
@@ -613,17 +632,20 @@ int main(int argc, char *argv[])
 	    ref_mode = lsat_rad2ref(ref_mode, &lsat.band[i]);
 	}
 
-	/* ================================================================= */
 
 	G_free(inrast);
-	Rast_close(infd);
 	G_free(outrast);
+	Rast_close(infd);
 	Rast_close(outfd);
 
-	/* needed ?
-	   if (out_type != CELL_TYPE)
-	   G_quantize_fp_map_range(band_out, G_mapset(), 0., 360., 0, 360);
+
+	/*
+	 * needed?
+	 * if (out_type != CELL_TYPE)
+	 * G_quantize_fp_map_range(band_out, G_mapset(), 0., 360., 0,
+	 * 360);
 	 */
+
 	/* set grey255 colortable */
 	Rast_init_colors(&colors);
 	Rast_read_fp_range(band_out, G_mapset(), &range);
@@ -631,32 +653,36 @@ int main(int argc, char *argv[])
 	Rast_make_grey_scale_fp_colors(&colors, min, max);
 	Rast_write_colors(band_out, G_mapset(), &colors);
 
-	/* Initialize the 'hist' structure with basic info */
+	/* Initialize the 'history' structure with basic info */
 	Rast_short_history(band_out, "raster", &history);
-	/* Append a string to the 'history' structure */
 	Rast_append_format_history(&history,
 				   " %s of Landsat-%d %s (method %s)",
-				   (frad->answer ? "Radiance" :
-				    (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,
-			    "----------------------------------------------------------------");
+			    "-----------------------------------------------------------------");
 	Rast_append_format_history(&history,
-				   " Acquisition date ...................... %s",
-				   lsat.date);
+				   " Acquisition date (and time) ........... %s (%.4lf h)",
+				   lsat.date, lsat.time);
 	Rast_append_format_history(&history,
 				   " Production date ....................... %s\n",
 				   lsat.creation);
 	Rast_append_format_history(&history,
-				   " Earth-sun distance (d) ................ %.8lf",
+				   " Earth-sun distance (d) ................ %.7lf",
 				   lsat.dist_es);
 	Rast_append_format_history(&history,
+				   " Sun elevation (and azimuth) ........... %.5lf (%.5lf)",
+				   lsat.sun_elev, lsat.sun_az);
+	Rast_append_format_history(&history,
 				   " Digital number (DN) range ............. %.0lf to %.0lf",
 				   lsat.band[i].qcalmin,
 				   lsat.band[i].qcalmax);
 	Rast_append_format_history(&history,
-				   " Calibration constants (Lmin to Lmax) .. %+.3lf to %+.3lf",
+				   " Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
 				   lsat.band[i].lmin, lsat.band[i].lmax);
 	Rast_append_format_history(&history,
 				   " DN to Radiance (gain and bias) ........ %+.5lf and %+.5lf",
@@ -671,10 +697,10 @@ int main(int argc, char *argv[])
 				       " Mean solar irradiance (ESUN) .......... %.3lf",
 				       lsat.band[i].esun);
 	    Rast_append_format_history(&history,
-				       " Reflectance = Radiance divided by ..... %.5lf",
+				       " Radiance to Reflectance (divide by) ... %+.5lf",
 				       lsat.band[i].K1);
 	    if (method > DOS) {
-		Rast_append_history(&history, " ");
+		Rast_append_format_history(&history, " ");
 		Rast_append_format_history(&history,
 					   " Dark object (%4d pixels) DN = ........ %d",
 					   pixel, dn_dark[i]);
@@ -684,7 +710,7 @@ int main(int argc, char *argv[])
 	    }
 	}
 	Rast_append_history(&history,
-			    "-----------------------------------------------------------------");
+			    "------------------------------------------------------------------");
 
 	Rast_command_history(&history);
 	Rast_write_history(band_out, &history);
@@ -692,7 +718,7 @@ int main(int argc, char *argv[])
 	if (lsat.band[i].thermal)
 	    Rast_write_units(band_out, "Kelvin");
 	else if (frad->answer)
-	    Rast_write_units(band_out, "W/(m^2 sr µm)");
+	    Rast_write_units(band_out, "W/(m^2 sr um)");
 	else
 	    Rast_write_units(band_out, "unitless");
     }