forward port of r53209: removed sat zenit (given by NASA for all their sats in metadata); sun_elev -> solar_elev to use same terminology as MTL files; improvement reads all metadata

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

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

@@ -6,17 +6,20 @@ top-of-atmosphere reflectance and temperature (band 6 of the sensors
 TM and ETM+). Optionally, it can be used to calculate the at-surface
 radiance or reflectance with atmospheric correction (DOS method).
 
-<p>Usually, to do so the production date, the acquisition date, and the
+<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
 needed the gain (high or low) of the nine respective bands.
 
-<p>Optionally, the data can be read from header file (.met) for all
+<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 metfile are
-overwriten. This is necessary when the data in the metfile is
+product creation date are given the values of the metadata file are
+overwriten (only with .met files). This is necessary when the data in the .met file is
 incorrect or not accurate.
 
-<p><b>Attention</b>: Any null value or smaller than QCALmin in the input
+<p>
+<b>Attention</b>: Any null value or smaller than QCALmin in the input
 raster is set to null in the output raster and it is not included in
 the equations.
 
@@ -30,7 +33,8 @@ reflectance. The thermal band is first converted from QCAL to
 at-sensor radiance, and then to effective at-sensor temperature in
 Kelvin degrees.
 
-<p>Radiometric calibration converts QCAL to <b>at-sensor radiance</b>, a
+<p>
+Radiometric calibration converts QCAL to <b>at-sensor radiance</b>, a
 radiometric quantity measured in W/(m&sup2; * sr * &micro;m) using the
 equations:
 <ul>
@@ -43,7 +47,8 @@ where, <em>Lmax</em> and <em>Lmin</em> are the calibration constants,
 and <em>QCALmax</em> and <em>QCALmin</em> are the highest and the
 lowest points of the range of rescaled radiance in QCAL.
 
-<p>Then, to calculate <b>at-sensor reflectance</b> the equations are:
+<p>
+Then, to calculate <b>at-sensor reflectance</b> the equations are:
 
 <ul>
   <li> sun_radiance = [Esun * sin(e)] / (PI * d^2)</li>
@@ -68,7 +73,8 @@ bands):
   <li> reflectance = radiance / sun_radiance </li>
 </ul>
 
-<p><b>Note</b>: Other possibility to avoid negative values is set to zero
+<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>.
@@ -120,12 +126,8 @@ negative value.
 In verbose mode (flag <b>--verbose</b>), the program write basic
 satellite data and the parameters used in the transformations.
 
-<p>In L5_MTL mode (flag <b>-t</b>), the Landsat 5TM imagery that has a
-_MTL.txt metadata file can be processed. Landsat 7 ETM+ does not need
-a flag since .met and _MTL.txt are sufficient compatible for this
-sensor.
-
-<p>Production date is not an exact value but it is necessary to apply
+<p>
+Production date is not an exact value but it is necessary to apply
 correct calibration constants, which were changed in the dates:
 <ul>
   <li>Landsat-1 MSS: never </li>
@@ -147,23 +149,23 @@ files 293_39.61_toar and 293_39.62_toar:
 
 <div class="code"><pre>
 i.landsat.toar input_prefix=203_30. output_prefix=_toar \
-  sensor=tm7 met=p203r030_7x20010620.met
+  metfile=p203r030_7x20010620.met
 </pre></div>
 
 or
 
 <div class="code"><pre>
-i.landsat.toar input_prefix=203_30. output_prefix=_toar \
-  sensor=tm7 product_date=2004-06-07 date=2001-06-20 \
-  solar_elevation=64.3242970 gain="HHHLHLHHL"
+i.landsat.toar input_prefix=L5121060_06020060714. \
+  output_prefix=L5121060_06020060714_toar \
+  metfile=L5121060_06020060714_MTL.txt
 </pre></div>
 
 or
 
 <div class="code"><pre>
-i.landsat.toar input_prefix=L5121060_06020060714. \
-  output_prefix=L5121060_06020060714_toar sensor=tm5 \
-  metfile=L5121060_06020060714_MTL.txt -t
+i.landsat.toar input_prefix=203_30. output_prefix=_toar \
+  sensor=tm7 product_date=2004-06-07 date=2001-06-20 \
+  sun_elevation=64.3242970 gain="HHHLHLHHL"
 </pre></div>
 
 <h2>REFERENCES</h2>
@@ -171,27 +173,27 @@ i.landsat.toar input_prefix=L5121060_06020060714. \
 <ul>
   <li>Chander G., B.L. Markham and D.L. Helder, 2009: Remote Sensing of
     Environment, vol. 113</li>
-  
+
   <li>Chander G.H. and B. Markham, 2003.: IEEE Transactions On Geoscience And
     Remote Sensing, vol. 41, no. 11.</li>
-  
+
   <li>Chavez P.S., jr. 1996. Image-based atmospheric corrections -
     Revisited and Improved. Photogrammetric Engineering and Remote
     Sensing 62(9): 1025-1036.</li>
-  
+
   <li>Huang et al: At-Satellite Reflectance, 2002: A First Order Normalization
     Of Landsat 7 ETM+ Images.</li>
-  
+
   <li>R. Irish: <a href="http://landsathandbook.gsfc.nasa.gov/orbit_coverage/">Landsat
       7. Science Data Users Handbook</a>. February 17, 2007; 15 May 2011.</li>
-  
+
   <li>Markham B.L. and J.L. Barker, 1986: Landsat MSS and TM Post-Calibration
     Dynamic Ranges, Exoatmospheric Reflectances and At-Satellite
     Temperatures. EOSAT Landsat Technical Notes, No. 1.</li>
-  
+
   <li>Moran M.S., R.D. Jackson, P.N. Slater and P.M. Teillet, 1992: Remote
     Sensing of Environment, vol. 41.</li>
-  
+
   <li>Song et al, 2001: Classification and Change Detection Using Landsat TM
     Data, When and How to Correct Atmospheric Effects? Remote Sensing
     of Environment, vol. 75.</li>
@@ -210,4 +212,5 @@ i.landsat.toar input_prefix=L5121060_06020060714. \
 E. Jorge Tizado  (ej.tizado unileon es), Dept. Biodiversity and Environmental Management,
 University of Le&oacute;n, Spain
 
-<p><i>Last changed: $Date$</i>
+<p>
+<i>Last changed: $Date$</i>

imagery/i.landsat.toar/landsat_met.c

@@ -1,5 +1,6 @@
 #include<stdio.h>
 #include<stdlib.h>
+#include<string.h>
 #include <math.h>
 
 #include <grass/gis.h>
@@ -8,278 +9,241 @@
 #include "local_proto.h"
 #include "earth_sun.h"
 
-#define ETM_MET_SIZE    5600	/* .met file size  5516 bytes */
-#define TM5_MET_SIZE    28700	/* .met file size 28686 bytes */
-#define MAX_STR      127
+#define MAX_STR         127
+#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)
+
+inline void chrncpy(char *dest, char src[], int n)
 {
     if (src == NULL)
-	n = 1;
+    {
+        dest[0] = '\0';
+    }
     else
-	strncpy(dest, src, n);
-    dest[n - 1] = '\0';
+    {
+        int i;
+        for (i = 0; i < n && src[i] != '\0' && src[i] != '\"'; i++) dest[i] = src[i];
+        dest[i] = '\0';
+    }
 }
 
+
 /****************************************************************************
- * PURPOSE:     Read values of Landsat-7 ETM+ from header (.met) file
+ * PURPOSE:     Read values of Landsat MSS/TM from header (.met) file
  *****************************************************************************/
-void get_value_met7(const char mettext[], char *text, char value[])
+void get_metdata(const char mettext[], char *text, char value[])
 {
     char *ptr;
 
-    value[0] = '\0';
+
 
     ptr = strstr(mettext, text);
     if (ptr == NULL)
-	return;
+    {
+        value[0] = 0;
+        return;
+    }
 
-    while (*ptr++ != '=') ;
-    sscanf(ptr, "%s", value);
+    ptr = strstr(ptr, " VALUE ");
+    if (ptr == NULL) return;
+
+    while (*ptr++ != '\"') ;
+    int i = 0;
+    while (*ptr != '\"' && i < MAX_STR) value[i++] = *ptr++;
+    value[i] = '\0';
 
     return;
 }
 
-void met_ETM(char *metfile, lsat_data * lsat)
+void lsat_metdata(char *metfile, lsat_data * lsat)
 {
     FILE *f;
-    char mettext[ETM_MET_SIZE];
+    char mettext[TM5_MET_SIZE];
     char name[MAX_STR], value[MAX_STR];
-    int i;
-
-    static int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
-    static int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
 
-    static double esun[] =
-	{ 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
+    /* char metdate[MAX_STR]; */
 
     if ((f = fopen(metfile, "r")) == NULL)
-	G_fatal_error(_("Metadata file <%s> not found"), metfile);
+    G_fatal_error(_("Metadata file <%s> not found"), metfile);
 
-    fread(mettext, 1, ETM_MET_SIZE, f);
+    fread(mettext, TM5_MET_SIZE, 1, f);
 
     /* --------------------------------------- */
-    lsat->number = 7;
+    get_metdata(mettext, "PLATFORMSHORTNAME", value);
+    chrncpy(name, value + 8, 1);
+    lsat->number = atoi(name);
 
-    get_value_met7(mettext, "SENSOR_ID", value);
-    chrncpy(lsat->sensor, value + 1, 5);
+    get_metdata(mettext, "SENSORSHORTNAME", value);
+    chrncpy(lsat->sensor, value + 1, 4);
 
-    if (lsat->creation[0] == 0) {
-	get_value_met7(mettext, "CREATION_TIME", value);
-	chrncpy(lsat->creation, value, 11);
-    }
+    get_metdata(mettext, "CALENDARDATE", value);
+    chrncpy(lsat->date, value, 10);
 
-    get_value_met7(mettext, "ACQUISITION_DATE", value);
-    chrncpy(lsat->date, value, 11);
-    lsat->dist_es = earth_sun(lsat->date);
+    if (lsat->creation[0] == 0)
+    {
+        get_metdata(mettext, "PRODUCTIONDATETIME", value);
+        if (!value[0])
+            G_fatal_error(_("Product creation date not in metadata file <%s>, input this data in the command line parameters"), metfile);
+        chrncpy(lsat->creation, value, 10);
+    }
 
-    get_value_met7(mettext, "SUN_ELEVATION", value);
-    lsat->sun_elev = atof(value);
+    if (lsat->sun_elev == 0)
+    {
+        get_metdata(mettext, "SolarElevation", value);
+        if (!value[0])
+           G_fatal_error(_("Unable to read solar elevation from metadata file in metadata file <%s>, input this data in the command line parameters"), metfile);
+        lsat->sun_elev = atof(value);
+    }
 
-    lsat->bands = 9;
-    for (i = 0; i < lsat->bands; i++) {
-	lsat->band[i].number = *(band + i);
-	lsat->band[i].code = *(code + i);
-	lsat->band[i].esun = *(esun + lsat->band[i].number - 1);
-	snprintf(name, MAX_STR, "LMAX_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].lmax = atof(value);
-	snprintf(name, MAX_STR, "LMIN_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].lmin = atof(value);
-	snprintf(name, MAX_STR, "QCALMAX_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].qcalmax = atof(value);
-	snprintf(name, MAX_STR, "QCALMIN_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].qcalmin = atof(value);
-	if (lsat->band[i].number == 6) {
-	    lsat->band[i].thermal = 1;
-	    lsat->band[i].K1 = 666.09;
-	    lsat->band[i].K2 = 1282.71;
-	}
-	else {
-	    lsat->band[i].thermal = 0;
-	}
+    /* Fill data with the sensor_XXX functions */
+    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;
+        default:
+            G_warning("Unable to recognize satellite platform [%d]", lsat->number);
+            break;
     }
 
+    /* --------------------------------------- */
     (void)fclose(f);
     return;
 }
 
+
 /****************************************************************************
- * PURPOSE:     Read values of Landsat MSS/TM from header (.met) file
+ * PURPOSE:     Read values of Landsat from MTL metadata (MTL.txt) file
  *****************************************************************************/
-void get_value_met(const char mettext[], char *text, char value[])
+
+void get_mtldata(const char mtltext[], char *text, char value[])
 {
     char *ptr;
-    int i;
 
-    value[0] = 0;
-
-    ptr = strstr(mettext, text);
+    ptr = strstr(mtltext, text);
     if (ptr == NULL)
-	return;
-
-    ptr = strstr(ptr, " VALUE ");
-    if (ptr == NULL)
-	return;
-
-    i = 0;
-    while (*ptr++ != '\"') ;
-    while (*ptr != '\"' && i < MAX_STR)
-	value[i++] = *ptr++;
-    value[i] = '\0';
-
-    return;
-}
-
-void met_TM5(char *metfile, lsat_data * lsat)
-{
-    FILE *f;
-    char mettext[TM5_MET_SIZE];
-    char value[MAX_STR];
-
-    /* char metdate[MAX_STR]; */
-
-    if ((f = fopen(metfile, "r")) == NULL)
-	G_fatal_error(_("Metadata file <%s> not found"), metfile);
-
-    fread(mettext, 1, TM5_MET_SIZE, f);
-
-    /* --------------------------------------- */
-    get_value_met(mettext, "CALENDARDATE", value);
-    chrncpy(lsat->date, value, 11);
-
-    if (lsat->creation[0] == 0) {
-	get_value_met(mettext, "PRODUCTIONDATETIME", value);
-	chrncpy(lsat->creation, value, 11);
+    {
+        value[0] = '\0';
+        return;
     }
 
-    if (lsat->creation[0] == 0)
-	G_fatal_error(_("Product creation date not in metadata file <%s>"),
-		      metfile);
-    G_debug(1, "met_TM5: Product creation date = [%s]", lsat->creation);
-
-
-    get_value_met(mettext, "SolarElevation", value);
-    if (!value[0])
-	G_warning("Unable to read solar elevation from metadata file");
-    else
-	lsat->sun_elev = atof(value);
-    G_debug(1, "met_TM5: value=[%s], SolarElevation = %.2f", value,
-	    lsat->sun_elev);
-
-
-    get_value_met(mettext, "PLATFORMSHORTNAME", value);
-    G_debug(1, "met_TM5: PLATFORMSHORTNAME=[%s]", value);
-    switch (value[8]) {
-    case '1':
-	set_MSS1(lsat);
-	break;
-    case '2':
-	set_MSS2(lsat);
-	break;
-    case '3':
-	set_MSS3(lsat);
-	break;
-    case '4':
-	get_value_met(mettext, "SENSORSHORTNAME", value);
-	if (value[0] == 'M')
-	    set_MSS4(lsat);
-	else
-	    set_TM4(lsat);
-	break;
-    case '5':
-	get_value_met(mettext, "SENSORSHORTNAME", value);
-	if (value[0] == 'M')
-	    set_MSS5(lsat);
-	else
-	    set_TM5(lsat);
-	break;
-    default:
-	G_warning("Unable to recognize satellite platform [%s]", value);
-	break;
-    }
+    while (*ptr++ != '=') ;
+    while (*ptr <= ' ' || *ptr == '\"') *ptr++;
+    int i = 0;
+    while (i < MAX_STR && *ptr != '\"' && *ptr > ' ') value[i++] = *ptr++;
+    value[i] = '\0';
 
-    (void)fclose(f);
     return;
 }
 
-
-
-/****************************************************************************
- * PURPOSE:     Read values of Landsat TM5 from header (.mtl) file
- *****************************************************************************/
-
-/****************************************************************************
- * EXPLANATION: This module is a modification of the met_ETM() found before
- *              to allow TM5 from GLOVIS to use .MTL extension that responds
- *              near to perfectly to the .met parser. While L7 files using
- *              .MTL from GLOVIS can be processed as if having .met files
- *              seemlessly, TM5 using .MTL need to read basic info and 
- *              additionally the LMIN, LMAX, QCALMIN, QCALMAX being explicitely
- *              provided in the .MTL as if in a .met file.
- *****************************************************************************/
-void mtl_TM5(char *metfile, lsat_data * lsat)
+void lsat_mtldata(char *mtlfile, lsat_data * lsat)
 {
     FILE *f;
-    char mettext[ETM_MET_SIZE];
+    char mtldata[METADATA_SIZE];
     char name[MAX_STR], value[MAX_STR];
     int i;
 
-    static int code[] = { 1, 2, 3, 4, 5, 6, 7 };
-
-    if ((f = fopen(metfile, "r")) == NULL)
-	G_fatal_error(_("Metadata file <%s> not found"), metfile);
+    if ((f = fopen(mtlfile, "r")) == NULL)
+       G_fatal_error(_("Metadata file <%s> not found"), mtlfile);
 
-    fread(mettext, 1, ETM_MET_SIZE, f);
+    fread(mtldata, METADATA_SIZE, 1, f);
 
     /* --------------------------------------- */
-    get_value_met7(mettext, "SENSOR_ID", value);
-    if (value[1] == 'M')
-	chrncpy(lsat->sensor, value + 1, 4);
-    else
-	chrncpy(lsat->sensor, value + 1, 3);
-    
-    if (lsat->creation[0] == 0) {
-	get_value_met7(mettext, "PRODUCT_CREATION_TIME", value);
-	chrncpy(lsat->creation, value, 11);
-    }
+    get_mtldata(mtldata, "SPACECRAFT_ID", value);
+    chrncpy(name, value + 7, 1);
+    lsat->number = atoi(name);
+
+    get_mtldata(mtldata, "SENSOR_ID", value);
+    chrncpy(lsat->sensor, value, 4);
 
-    get_value_met7(mettext, "ACQUISITION_DATE", value);
-    chrncpy(lsat->date, value, 11);
-    lsat->dist_es = earth_sun(lsat->date);
+    get_mtldata(mtldata, "ACQUISITION_DATE", value);
+    chrncpy(lsat->date, value, 10);
 
-    get_value_met7(mettext, "SUN_ELEVATION", value);
+    get_mtldata(mtldata, "CREATION_TIME", value);
+    chrncpy(lsat->creation, value, 10);
+
+    get_mtldata(mtldata, "SUN_ELEVATION", value);
     lsat->sun_elev = atof(value);
 
-    /* We still have to initialize most of the info */
-    /* So instead of rewriting a new function, we use set_TM5()... */
-    set_TM5(lsat);
-    /* ... and we rewrite the necessary 'a la Landsat 7' */
+    /* Fill data with the sensor_XXX functions */
+    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:
+            get_mtldata(mtldata, "BAND1_GAIN",  value);
+            get_mtldata(mtldata, "BAND2_GAIN",  value + 1);
+            get_mtldata(mtldata, "BAND3_GAIN",  value + 2);
+            get_mtldata(mtldata, "BAND4_GAIN",  value + 3);
+            get_mtldata(mtldata, "BAND5_GAIN",  value + 4);
+            get_mtldata(mtldata, "BAND6_GAIN1", value + 5);
+            get_mtldata(mtldata, "BAND6_GAIN2", value + 6);
+            get_mtldata(mtldata, "BAND7_GAIN",  value + 7);
+            get_mtldata(mtldata, "BAND8_GAIN",  value + 8);
+            value[9] = '\0';
+            set_ETM(lsat, value);
+            break;
+        default:
+            G_warning("Unable to recognize satellite platform [%d]", lsat->number);
+            break;
+    }
 
-    if (strcmp(lsat->sensor, "MSS") == 0)
-	lsat->bands = 4;
-    else
-	lsat->bands = 7;
+    /* Update the information from metadata file */
     for (i = 0; i < lsat->bands; i++) {
-	lsat->band[i].code = *(code + i);
-	snprintf(name, MAX_STR, "LMAX_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].lmax = atof(value);
-	snprintf(name, MAX_STR, "LMIN_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].lmin = atof(value);
-	snprintf(name, MAX_STR, "QCALMAX_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].qcalmax = atof(value);
-	snprintf(name, MAX_STR, "QCALMIN_BAND%d", lsat->band[i].code);
-	get_value_met7(mettext, name, value);
-	lsat->band[i].qcalmin = atof(value);
-	if (lsat->band[i].number == 6)
-	    lsat->band[i].thermal = 1;
+        snprintf(name, MAX_STR, "LMAX_BAND%d", lsat->band[i].code);
+        get_mtldata(mtldata, name, value);
+        lsat->band[i].lmax = atof(value);
+        snprintf(name, MAX_STR, "LMIN_BAND%d", lsat->band[i].code);
+        get_mtldata(mtldata, name, value);
+        lsat->band[i].lmin = atof(value);
+        snprintf(name, MAX_STR, "QCALMAX_BAND%d", lsat->band[i].code);
+        get_mtldata(mtldata, name, value);
+        lsat->band[i].qcalmax = atof(value);
+        snprintf(name, MAX_STR, "QCALMIN_BAND%d", lsat->band[i].code);
+        get_mtldata(mtldata, name, value);
+        lsat->band[i].qcalmin = atof(value);
     }
+    /* --------------------------------------- */
+
     (void)fclose(f);
     return;
 }
+

imagery/i.landsat.toar/landsat_set.c

@@ -18,7 +18,7 @@ void sensor_MSS(lsat_data * lsat)
     int code[] = { 4, 5, 6, 7 };
     double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
     double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
-    /* 79-82, 79-82, 79-82, 79-82 */
+    /* 68x83, 68x83, 68x83, 68x83 */
 
     strcpy(lsat->sensor, "MSS");
 
@@ -89,7 +89,7 @@ void sensor_ETM(lsat_data * lsat)
 
 
 /** **********************************************
- ** Before access to this function ...
+ ** Before access to these functions ...
  ** store previously
  ** >>> adquisition date,
  ** >>> creation date, and
@@ -98,21 +98,20 @@ void sensor_ETM(lsat_data * lsat)
 
 /****************************************************************************
  * PURPOSE:     Store values of Landsat-1 MSS
- *              July 23, 1972 to January 7, 1978
+ *              July 23, 1972 to January 6, 1978
  *****************************************************************************/
 void set_MSS1(lsat_data * lsat)
 {
     int i, j;
 
-    /** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
-        Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
+    /** USGS Calibration Parameter Files 2012 */
 
     /* Spectral radiances at detector */
     double lmax[] = { 248., 200., 176., 153. };
     double lmin[] = { 0., 0., 0., 0. };
 
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1823., 1559., 1276., 880.1 };
+    double esun[] = { 1824., 1570., 1249., 853.4 };
 
     lsat->number = 1;
     sensor_MSS(lsat);
@@ -125,33 +124,33 @@ void set_MSS1(lsat_data * lsat)
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-1 MSS");
     return;
 }
 
 /****************************************************************************
  * PURPOSE:     Store values of Landsat-2 MSS
- *              January 22, 1975 to February 25, 1982
+ *              January 22, 1975 to February 5, 1982
  *****************************************************************************/
 void set_MSS2(lsat_data * lsat)
 {
     int i, j;
     double julian, *lmax, *lmin;
 
-    /** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
-        Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
+    /** USGS Calibration Parameter Files 2012 */
 
     /* Spectral radiances at detector */
     double Lmax[][4] = {
 	{210., 156., 140., 138.},	/* before      July 16, 1975 */
-	{263., 176., 152., 130.333}	/* on or after July 16, 1975 */
+	{263., 176., 152., 130.}	/* on or after July 16, 1975 */
     };
     double Lmin[][4] = {
 	{10., 7., 7., 5.},
-	{8., 6., 6., 3.667}
+	{8., 6., 6., 4.}
     };
 
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1829., 1539., 1268., 886.6 };
+    double esun[] = { 1824., 1570., 1249., 853.4 };
 
     julian = julian_char(lsat->creation);
     if (julian < julian_char("1975-07-16"))
@@ -172,6 +171,7 @@ void set_MSS2(lsat_data * lsat)
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-2 MSS");
     return;
 }
 
@@ -186,8 +186,8 @@ void set_MSS3(lsat_data * lsat)
     int i, j;
     double julian, *lmax, *lmin;
 
-    /** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
-        Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
+    /** USGS Calibration Parameter Files 2012 */
+
     /* Spectral radiances at detector */
     double Lmax[][4] = {
 	{220., 175., 145., 147.},	/* before      June 1, 1978 */
@@ -198,7 +198,7 @@ void set_MSS3(lsat_data * lsat)
 	{4., 3., 3., 1.}
     };
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1839., 1555., 1291., 887.9 };
+    double esun[] = { 1824., 1570., 1249., 853.4 };
 
     julian = julian_char(lsat->creation);
     if (julian < julian_char("1978-06-01"))
@@ -219,6 +219,7 @@ void set_MSS3(lsat_data * lsat)
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-3 MSS");
     return;
 }
 
@@ -231,8 +232,7 @@ void set_MSS4(lsat_data * lsat)
     int i, j;
     double julian, *lmax, *lmin;
 
-    /** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
-        Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
+    /** USGS Calibration Parameter Files 2012 */
 
     /* Spectral radiances at detector */
     double Lmax[][4] = {
@@ -247,15 +247,15 @@ void set_MSS4(lsat_data * lsat)
     };
 
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1827., 1569., 1260., 866.4 };
+    double esun[] = { 1824., 1570., 1249., 853.4 };
 
     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];
 
@@ -270,6 +270,7 @@ void set_MSS4(lsat_data * lsat)
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-4 MSS");
     return;
 }
 
@@ -278,34 +279,32 @@ void set_TM4(lsat_data * lsat)
     int i, j;
     double julian, *lmax, *lmin;
 
-    /** Brian L. Markham and John L. Barker.
-        EOSAT Landsat Technical Notes, No. 1, 1986 */
+    /** USGS Calibration Parameter Files 2012 */
+
     /* Spectral radiances at detector */
     double Lmax[][7] = {
 	{158.42, 308.17, 234.63, 224.32, 32.42, 15.64, 17.00},	/* before August 1983      */
 	{142.86, 291.25, 225.00, 214.29, 30.00, 12.40, 15.93},	/* before January 15, 1984 */
-	{152.10, 296.81, 204.30, 206.20, 27.19, 15.3032, 14.38}	/* after  Jaunary 15, 1984 */
+	{171.00, 336.00, 254.00, 221.00, 31.40, 15.303, 16.60}	/* after  Jaunary 15, 1984 */
     };
     double Lmin[][7] = {
 	{-1.52, -2.84, -1.17, -1.51, -0.37, 2.00, -0.15},
 	{0.00, 0.00, 0.00, 0.00, 0.00, 4.84, 0.00},
-	{-1.50, -2.80, -1.20, -1.50, -0.37, 1.2378, -0.15}
+	{-1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15}
     };
 
-    /** Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009) */
-
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1983., 1795., 1539., 1028., 219.8, 0., 83.49 };
+    double esun[] = { 1957., 1825., 1557., 1033., 214.9, 0., 80.72 };
 
     /* Thermal band calibration constants: K1 = 671.62   K2 = 1284.30 */
 
     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];
 
@@ -324,6 +323,7 @@ void set_TM4(lsat_data * lsat)
 	    lsat->band[i].K2 = 1284.30;
 	}
     }
+    G_debug(1, "Landsat-4 TM");
     return;
 }
 
@@ -337,8 +337,8 @@ void set_MSS5(lsat_data * lsat)
     int i, j;
     double julian, *lmax, *lmin;
 
-    /** Brian L. Markham and John L. Barker.
-        EOSAT Landsat Technical Notes, No. 1, 1986 */
+    /** USGS Calibration Parameter Files 2012 */
+
     /* Spectral radiances at detector */
     double Lmax[][4] = {
 	{240., 170., 150., 127.},	/* before   April 6, 1984    */
@@ -357,7 +357,7 @@ void set_MSS5(lsat_data * lsat)
     julian = julian_char(lsat->creation);
     if (julian < julian_char("1984-04-06"))
 	i = 0;
-    else if (julian < julian_char("1984-11-09"))
+    else if (julian < julian_char("1984-11-08"))
 	i = 1;
     else
 	i = 2;
@@ -375,6 +375,7 @@ void set_MSS5(lsat_data * lsat)
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-5 MSS");
     return;
 }
 
@@ -383,8 +384,7 @@ void set_TM5(lsat_data * lsat)
     int i, j;
     double julian, *lmax, *lmin, jbuf;
 
-    /** Gyanesh Chander and Brian Markham.
-        IEEE Transactions On Geoscience And Remote Sensing, Vol. 41, No. 11, November 2003 */
+    /** USGS Calibration Parameter Files 2012 */
 
     /* Spectral radiances at detector */
     double Lmax[][7] = {
@@ -398,20 +398,18 @@ void set_TM5(lsat_data * lsat)
 	{-1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15}
     };
 
-	/** Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009) */
-
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1983., 1796., 1536., 1031., 220.0, 0., 83.44 };
+    double esun[] = { 1957., 1826., 1554., 1036., 215.0, 0., 80.67 };
 
     /* Thermal band calibration constants: K1 = 607.76   K2 = 1260.56 */
 
     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];
     if (i == 2) {		/* in Chander, Markham and Barsi 2007 */
@@ -445,21 +443,21 @@ void set_TM5(lsat_data * lsat)
 	    lsat->band[i].K2 = 1260.56;
 	}
     }
+    G_debug(1, "Landsat-5 TM");
     return;
 }
 
 
 /****************************************************************************
  * PURPOSE:     Store values of Landsat-7 ETM+
- *              April 15, 1999 to today
+ *              April 15, 1999 to May 31, 2003 (SLC failure)
  *****************************************************************************/
 void set_ETM(lsat_data * lsat, char gain[])
 {
     int i, k, j;
     double julian, *lmax, *lmin;
 
-    /** Richard Irish.
-        Landsat 7. Science Data Users Handbook. Last update: February 17, 2007 */
+    /** USGS Calibration Parameter Files 2012 */
 
     /* Spectral radiances at detector */
     /* - LOW GAIN - */
@@ -481,10 +479,8 @@ void set_ETM(lsat_data * lsat, char gain[])
 	{-6.2, -6.4, -5.0, -5.1, -1.0, 3.2, -0.35, -4.7}
     };
 
-	/** Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009) */
-
     /* Solar exoatmospheric spectral irradiances */
-    double esun[] = { 1997., 1812., 1533., 1039., 230.8, 0., 84.90, 1362. };
+    double esun[] = { 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
 
     /*  Thermal band calibration constants: K1 = 666.09   K2 = 1282.71 */
 
@@ -517,5 +513,6 @@ void set_ETM(lsat_data * lsat, char gain[])
 	    lsat->band[i].K2 = 1282.71;
 	}
     }
+    G_debug(1, "Landsat-7 ETM+");
     return;
 }

imagery/i.landsat.toar/local_proto.h

@@ -4,9 +4,8 @@
 #include <string.h>
 #include "landsat.h"
 
-void met_ETM(char *, lsat_data *);
-void met_TM5(char *, lsat_data *);
-void mtl_TM5(char *, lsat_data *);
+void lsat_mtldata(char *, lsat_data *);
+void lsat_metdata(char *, lsat_data *);
 
 void set_MSS1(lsat_data *);
 void set_MSS2(lsat_data *);
@@ -16,6 +15,7 @@ void set_MSS5(lsat_data *);
 
 void set_TM4(lsat_data *);
 void set_TM5(lsat_data *);
+
 void set_ETM(lsat_data *, char[]);
 
 #endif

imagery/i.landsat.toar/main.c

@@ -4,14 +4,14 @@
  * MODULE:       i.landsat.toar
  *
  * AUTHOR(S):    E. Jorge Tizado - ej.tizado@unileon.es
- *		 Hamish Bowman (small grassification cleanups)
- *               Yann Chemin (v7 + L5TM _MTL.txt support)
+ *               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>
  *
  * PURPOSE:      Calculate TOA Radiance or Reflectance and Kinetic Temperature
  *               for Landsat 1/2/3/4/5 MS, 4/5 TM or 7 ETM+
  *
- * COPYRIGHT:    (C) 2002, 2005, 2008, 2010 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2002-2012 by the GRASS Development Team
  *
  *               This program is free software under the GNU General
  *               Public License (>=v2). Read the file COPYING that
@@ -45,8 +45,8 @@ int main(int argc, char *argv[])
     struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate, *pdate, *elev,
 	*bgain, *metho, *perc, *dark, *atmo;
     char *inputname, *met, *outputname, *sensorname;
-    struct Flag *msss, *frad, *l5_mtl;
-    
+    struct Flag *frad;
+
     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];
@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
     metfn = G_define_standard_option(G_OPT_F_INPUT);
     metfn->key = "metfile";
     metfn->required = NO;
-    metfn->description = _("Name of Landsat ETM+ or TM5 header file (.met/MTL.txt)");
+    metfn->description = _("Name of Landsat metadata file (.met or MTL.txt)");
     metfn->guisection = _("Metadata");
 
     sensor = G_define_option();
@@ -98,13 +98,15 @@ int main(int argc, char *argv[])
     sensor->type = TYPE_STRING;
     sensor->label = _("Spacecraft sensor");
     sensor->description = _("Required only if 'metfile' not given");
-    sensor->options = "mss1,mss2,mss3,tm4,tm5,tm7";
+    sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7";
     desc = NULL;
     G_asprintf(&desc,
-	        "mss1;%s;mss2;%s;mss3;%s;tm4;%s;tm5;%s;tm7;%s",
+	        "mss1;%s;mss2;%s;mss3;%s;mss4;%smss5;%stm4;%s;tm5;%s;tm7;%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+"));
@@ -118,10 +120,10 @@ int main(int argc, char *argv[])
     metho->required = NO;
     metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
     metho->label = _("Atmospheric correction method");
-    metho->description = _("Required only if 'metfile' not given");
+    metho->description = _("Atmospheric correction method");
     metho->answer = "uncorrected";
     metho->guisection = _("Metadata");
-    
+
     adate = G_define_option();
     adate->key = "date";
     adate->type = TYPE_STRING;
@@ -130,12 +132,12 @@ int main(int argc, char *argv[])
     adate->label = _("Image acquisition date (yyyy-mm-dd)");
     adate->description = _("Required only if 'metfile' not given");
     adate->guisection = _("Metadata");
-    
+
     elev = G_define_option();
-    elev->key = "solar_elevation";
+    elev->key = "sun_elevation";
     elev->type = TYPE_DOUBLE;
     elev->required = NO;
-    elev->label = _("Solar elevation in degrees");
+    elev->label = _("Sun elevation in degrees");
     elev->description = _("Required only if 'metfile' not given");
     elev->guisection = _("Metadata");
 
@@ -152,8 +154,8 @@ 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");
 
     perc = G_define_option();
@@ -177,7 +179,7 @@ int main(int argc, char *argv[])
     atmo->key = "rayleigh";
     atmo->type = TYPE_DOUBLE;
     atmo->required = NO;
-    atmo->description = _("Rayleigh atmosphere");	/* scattering coefficient? */
+    atmo->description = _("Rayleigh atmosphere (diffuse sky irradiance)");	/* scattering coefficient? */
     atmo->answer = "0.0";
     atmo->guisection = _("Settings");
 
@@ -185,16 +187,6 @@ int main(int argc, char *argv[])
     frad = G_define_flag();
     frad->key = 'r';
     frad->description = _("Output at-sensor radiance for all bands");
-    
-    msss = G_define_flag();
-    msss->key = 's';
-    msss->description = _("Set sensor of Landsat TM4/5 to MSS");
-    msss->guisection = _("Settings");
-
-    l5_mtl = G_define_flag();
-    l5_mtl->key = 't';
-    l5_mtl->description = _("Landsat ETM+/TM5 has a MTL.txt file instead of .met");
-    l5_mtl->guisection = _("Metadata");
 
     /* options and afters parser */
     if (G_parser(argc, argv))
@@ -208,12 +200,12 @@ 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);
-    
+
     G_zero(&lsat, sizeof(lsat));
-    
+
     if (adate->answer != NULL) {
 	strncpy(lsat.date, adate->answer, 11);
 	lsat.date[10] = '\0';
@@ -239,24 +231,30 @@ int main(int argc, char *argv[])
     pixel = atoi(dark->answer);
     rayleigh = atof(atmo->answer);
 
-    /* Data from MET file: only Landsat-7 ETM+ and Landsat-5 TM  */
-    if (met != NULL) {
-	if (strcmp(sensorname, "tm7") == 0)
-	    met_ETM(met, &lsat);
-	else if (l5_mtl->answer)
-	    mtl_TM5(met, &lsat);
-	else
-	    met_TM5(met, &lsat);
-
-	G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number,
-		lsat.sensor);
-	if (!lsat.sensor || lsat.number > 7 || 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);
-	if (elev->answer != NULL)
-	    lsat.sun_elev = atof(elev->answer);	/* Overwrite solar elevation of met file */
+    /* Data from metadata file */
+    if (met != NULL)
+    {
+        i = strlen(met);
+        if (strcmp(met + i - 7, "MTL.txt") == 0)
+        {
+            lsat_mtldata(met, &lsat);
+        }
+        else if (strcmp(met + i - 4, ".met") == 0)
+        {
+            if (strcmp(sensorname, "tm7") == 0)
+                lsat_mtldata(met, &lsat);  /* .met of Landsat-7 = new MTL file */
+            else
+                lsat_metdata(met, &lsat);
+        }
+        G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number, lsat.sensor);
+        if (!lsat.sensor || lsat.number > 7 || 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);
+
+        /* Overwrite solar elevation of metadata file */
+        if (elev->answer != NULL)
+            lsat.sun_elev = atof(elev->answer);
     }
     /* Data from date and solar elevation */
     else if (adate->answer == NULL || elev->answer == NULL) {
@@ -265,45 +263,32 @@ int main(int argc, char *argv[])
     }
     else {
 	if (strcmp(sensorname, "tm7") == 0) {	/* Need gain */
-	    if (bgain->answer != NULL && strlen(bgain->answer) == 9) {
-		set_ETM(&lsat, bgain->answer);
-		G_debug(1, "Landsat 7 ETM+");
-	    }
-	    else {
-		G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) data"));
-	    }
-	}
-	else {			/* Not need gain */
-	    if (strcmp(sensorname, "tm5") == 0) {
-		if (msss->answer)
-		    set_MSS5(&lsat);
-		else
-		    set_TM5(&lsat);
-		G_debug(1, "Landsat-5 %s", lsat.sensor);
-	    }
-	    else if (strcmp(sensorname, "tm4") == 0) {
-		if (msss->answer)
-		    set_MSS4(&lsat);
-		else
-		    set_TM4(&lsat);
-		G_debug(1, "Landsat-4 %s", lsat.sensor);
-	    }
-	    else if (strcmp(sensorname, "mss3") == 0) {
-		set_MSS3(&lsat);
-		G_debug(1, "Landsat-3 MSS");
-	    }
-	    else if (strcmp(sensorname, "mss2") == 0) {
-		set_MSS2(&lsat);
-		G_debug(1, "Landsat-2 MSS");
-	    }
-	    else if (strcmp(sensorname, "mss1") == 0) {
-		set_MSS1(&lsat);
-		G_debug(1, "Landsat-1 MSS");
-	    }
-	    else {
-		G_fatal_error(_("Unknown satellite type (defined by '%s')"), sensor->key);
-	    }
-	}
+        if (bgain->answer != NULL && strlen(bgain->answer) == 9) {
+            set_ETM(&lsat, bgain->answer);
+            G_debug(1, "Landsat 7 ETM+");
+        }
+        else {
+            G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) data"));
+        }
+    }
+	else {  /* Not need gain */
+        if (strcmp(sensorname, "tm5") == 0)
+            set_TM5(&lsat);
+        else if (strcmp(sensorname, "tm4") == 0)
+            set_TM4(&lsat);
+        else if (strcmp(sensorname, "mss5") == 0)
+            set_MSS5(&lsat);
+        else if (strcmp(sensorname, "mss4") == 0)
+            set_MSS4(&lsat);
+        else if (strcmp(sensorname, "mss3") == 0)
+            set_MSS3(&lsat);
+        else if (strcmp(sensorname, "mss2") == 0)
+            set_MSS2(&lsat);
+        else if (strcmp(sensorname, "mss1") == 0)
+            set_MSS1(&lsat);
+        else
+            G_fatal_error(_("Unknown satellite type (defined by '%s')"), sensor->key);
+        }
     }
 
 	/*****************************************
@@ -409,7 +394,7 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Unknown production date (defined by '%s')"), pdate->key);
 
     if (G_verbose() > G_verbose_std()) {
-	fprintf(stderr, " SENSOR: %s\n", 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);