[GRASS-SVN] r53209 - grass-addons/grass6/imagery/i.landsat.toar

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Sep 18 09:29:18 PDT 2012


Author: ejtizado
Date: 2012-09-18 09:29:17 -0700 (Tue, 18 Sep 2012)
New Revision: 53209

Modified:
   grass-addons/grass6/imagery/i.landsat.toar/description.html
   grass-addons/grass6/imagery/i.landsat.toar/landsat_met.c
   grass-addons/grass6/imagery/i.landsat.toar/landsat_set.c
   grass-addons/grass6/imagery/i.landsat.toar/local_proto.h
   grass-addons/grass6/imagery/i.landsat.toar/main.c
Log:
Update to read metadata from .met and MTL.txt for all landsats

Modified: grass-addons/grass6/imagery/i.landsat.toar/description.html
===================================================================
--- grass-addons/grass6/imagery/i.landsat.toar/description.html	2012-09-18 15:09:24 UTC (rev 53208)
+++ grass-addons/grass6/imagery/i.landsat.toar/description.html	2012-09-18 16:29:17 UTC (rev 53209)
@@ -12,10 +12,10 @@
 needed the gain (high or low) of the nine respective bands.
 
 <p>
-Optionally, the data can be read from header file (.met) for all
+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>
@@ -127,12 +127,6 @@
 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
 correct calibration constants, which were changed in the dates:
 <ul>
@@ -155,23 +149,23 @@
 
 <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>

Modified: grass-addons/grass6/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass-addons/grass6/imagery/i.landsat.toar/landsat_met.c	2012-09-18 15:09:24 UTC (rev 53208)
+++ grass-addons/grass6/imagery/i.landsat.toar/landsat_met.c	2012-09-18 16:29:17 UTC (rev 53209)
@@ -1,5 +1,6 @@
 #include<stdio.h>
 #include<stdlib.h>
+#include<string.h>
 #include <math.h>
 
 #include <grass/gis.h>
@@ -8,199 +9,239 @@
 #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 };
+    /* char metdate[MAX_STR]; */
 
-    static double esun[] =
-	{ 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
-
     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);
+
+    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, "ACQUISITION_DATE", value);
-    chrncpy(lsat->date, value, 11);
-    lsat->dist_es = earth_sun(lsat->date);
+    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);
+    }
 
-    get_value_met7(mettext, "SUN_ELEVATION", value);
-    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;
+    {
+        value[0] = '\0';
+        return;
+    }
 
-    ptr = strstr(ptr, " VALUE ");
-    if (ptr == NULL)
-	return;
-
-    i = 0;
-    while (*ptr++ != '\"') ;
-    while (*ptr != '\"' && i < MAX_STR)
-	value[i++] = *ptr++;
+    while (*ptr++ != '=') ;
+    while (*ptr <= ' ' || *ptr == '\"') *ptr++;
+    int i = 0;
+    while (i < MAX_STR && *ptr != '\"' && *ptr > ' ') value[i++] = *ptr++;
     value[i] = '\0';
 
     return;
 }
 
-void met_TM5(char *metfile, lsat_data * lsat)
+void lsat_mtldata(char *mtlfile, lsat_data * lsat)
 {
     FILE *f;
-    char mettext[TM5_MET_SIZE];
-    char value[MAX_STR];
+    char mtldata[METADATA_SIZE];
+    char name[MAX_STR], value[MAX_STR];
+    int i;
 
-    /* char metdate[MAX_STR]; */
+    if ((f = fopen(mtlfile, "r")) == NULL)
+       G_fatal_error(_("Metadata file <%s> not found"), mtlfile);
 
-    if ((f = fopen(metfile, "r")) == NULL)
-	G_fatal_error(_("Metadata file <%s> not found"), metfile);
+    fread(mtldata, METADATA_SIZE, 1, f);
 
-    fread(mettext, 1, TM5_MET_SIZE, f);
-
     /* --------------------------------------- */
-    get_value_met(mettext, "CALENDARDATE", value);
-    chrncpy(lsat->date, value, 11);
+    get_mtldata(mtldata, "SPACECRAFT_ID", value);
+    chrncpy(name, value + 7, 1);
+    lsat->number = atoi(name);
 
-    if (lsat->creation[0] == 0) {
-	get_value_met(mettext, "PRODUCTIONDATETIME", value);
-	chrncpy(lsat->creation, value, 11);
-    }
+    get_mtldata(mtldata, "SENSOR_ID", value);
+    chrncpy(lsat->sensor, value, 4);
 
-    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_mtldata(mtldata, "ACQUISITION_DATE", value);
+    chrncpy(lsat->date, value, 10);
 
+    get_mtldata(mtldata, "CREATION_TIME", value);
+    chrncpy(lsat->creation, value, 10);
 
-    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_mtldata(mtldata, "SUN_ELEVATION", value);
+    lsat->sun_elev = atof(value);
 
+    // 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;
+    }
 
-    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;
+    // Update the information from metadata file
+    for (i = 0; i < lsat->bands; i++) {
+        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;
@@ -208,78 +249,3 @@
 
 
 
-/****************************************************************************
- * 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)
-{
-    FILE *f;
-    char mettext[ETM_MET_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);
-
-    fread(mettext, 1, ETM_MET_SIZE, 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_value_met7(mettext, "ACQUISITION_DATE", value);
-    chrncpy(lsat->date, value, 11);
-    lsat->dist_es = earth_sun(lsat->date);
-
-    get_value_met7(mettext, "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' */
-
-    if (strcmp(lsat->sensor, "MSS") == 0)
-	lsat->bands = 4;
-    else
-	lsat->bands = 7;
-    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;
-    }
-    (void)fclose(f);
-    return;
-}

Modified: grass-addons/grass6/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass-addons/grass6/imagery/i.landsat.toar/landsat_set.c	2012-09-18 15:09:24 UTC (rev 53208)
+++ grass-addons/grass6/imagery/i.landsat.toar/landsat_set.c	2012-09-18 16:29:17 UTC (rev 53209)
@@ -18,7 +18,7 @@
     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 @@
 
 
 /** **********************************************
- ** Before access to this function ...
+ ** Before access to these functions ...
  ** store previously
  ** >>> adquisition date,
  ** >>> creation date, and
@@ -98,21 +98,20 @@
 
 /****************************************************************************
  * 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 @@
 	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 @@
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-2 MSS");
     return;
 }
 
@@ -186,8 +186,8 @@
     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 @@
 	{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 @@
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-3 MSS");
     return;
 }
 
@@ -231,8 +232,7 @@
     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 @@
     };
 
     /* 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 @@
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-4 MSS");
     return;
 }
 
@@ -278,34 +279,32 @@
     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 @@
 	    lsat->band[i].K2 = 1284.30;
 	}
     }
+    G_debug(1, "Landsat-4 TM");
     return;
 }
 
@@ -337,8 +337,8 @@
     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 @@
     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 @@
 	lsat->band[i].lmax = *(lmax + j);
 	lsat->band[i].lmin = *(lmin + j);
     }
+    G_debug(1, "Landsat-5 MSS");
     return;
 }
 
@@ -383,8 +384,7 @@
     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 @@
 	{-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 @@
 	    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 @@
 	{-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 @@
 	    lsat->band[i].K2 = 1282.71;
 	}
     }
+    G_debug(1, "Landsat-7 ETM+");
     return;
 }

Modified: grass-addons/grass6/imagery/i.landsat.toar/local_proto.h
===================================================================
--- grass-addons/grass6/imagery/i.landsat.toar/local_proto.h	2012-09-18 15:09:24 UTC (rev 53208)
+++ grass-addons/grass6/imagery/i.landsat.toar/local_proto.h	2012-09-18 16:29:17 UTC (rev 53209)
@@ -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_TM4(lsat_data *);
 void set_TM5(lsat_data *);
+
 void set_ETM(lsat_data *, char[]);
 
 #endif

Modified: grass-addons/grass6/imagery/i.landsat.toar/main.c
===================================================================
--- grass-addons/grass6/imagery/i.landsat.toar/main.c	2012-09-18 15:09:24 UTC (rev 53208)
+++ grass-addons/grass6/imagery/i.landsat.toar/main.c	2012-09-18 16:29:17 UTC (rev 53209)
@@ -4,13 +4,13 @@
  * MODULE:       i.landsat.toar
  *
  * AUTHOR(S):    E. Jorge Tizado - ej.tizado at 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]
  *
  * 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) 2012 by the GRASS Development Team
  *
  *               This program is free software under the GNU General
  *               Public License (>=v2). Read the file COPYING that
@@ -43,7 +43,7 @@
     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];
@@ -82,35 +82,37 @@
     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");
 
+    metho = G_define_option();
+    metho->key = "method";
+    metho->type = TYPE_STRING;
+    metho->required = NO;
+    metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
+    metho->label = _("Atmospheric correction method");
+    metho->description = _("Atmospheric correction method");
+    metho->answer = "uncorrected";
+    metho->guisection = _("Metadata");
+
     sensor = G_define_option();
     sensor->key = "sensor";
     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";
     sensor->descriptions =
 	_("mss1;Landsat-1 MSS;"
 	  "mss2;Landsat-2 MSS;"
-	  "mss3;Landsat-3 MSS;"
-	  "tm4;Landsat-4 TM;"
-	  "tm5;Landsat-5 TM;"
+      "mss3;Landsat-3 MSS;"
+      "mss4;Landsat-4 MSS;"
+      "mss5;Landsat-5 MSS;"
+      "tm4;Landsat-4 TM;"
+      "tm5;Landsat-5 TM;"
 	  "tm7;Landsat-7 ETM+");
     sensor->required = NO;
     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->label = _("Atmospheric correction method");
-    metho->description = _("Required only if 'metfile' not given");
-    metho->answer = "uncorrected";
-    metho->guisection = _("Metadata");
-
     adate = G_define_option();
     adate->key = "date";
     adate->type = TYPE_STRING;
@@ -121,10 +123,10 @@
     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");
 
@@ -141,8 +143,8 @@
     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();
@@ -166,25 +168,16 @@
     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");
 
+
     /* define the different flags */
     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))
 	exit(EXIT_FAILURE);
@@ -197,7 +190,7 @@
     met = metfn->answer;
     inputname = input_prefix->answer;
     outputname = output_prefix->answer;
-    sensorname = sensor -> answer ? sensor->answer: "";
+    sensorname = sensor->answer ? sensor->answer: "";
 
     G_zero(&lsat, sizeof(lsat));
 
@@ -226,24 +219,30 @@
     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);
+    /* 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, "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);
 
-	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 */
+        /* 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) {
@@ -252,61 +251,48 @@
     }
     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);
+        }
+    }
 
 	/*****************************************
 	* ------------ PREPARATION --------------
 	*****************************************/
-    if (G_strcasecmp(metho->answer, "corrected") == 0)
+    if (strcasecmp(metho->answer, "corrected") == 0)
 	method = CORRECTED;
-    else if (G_strcasecmp(metho->answer, "dos1") == 0)
+    else if (strcasecmp(metho->answer, "dos1") == 0)
 	method = DOS1;
-    else if (G_strcasecmp(metho->answer, "dos2") == 0)
+    else if (strcasecmp(metho->answer, "dos2") == 0)
 	method = DOS2;
-    else if (G_strcasecmp(metho->answer, "dos2b") == 0)
+    else if (strcasecmp(metho->answer, "dos2b") == 0)
 	method = DOS2b;
-    else if (G_strcasecmp(metho->answer, "dos3") == 0)
+    else if (strcasecmp(metho->answer, "dos3") == 0)
 	method = DOS3;
-    else if (G_strcasecmp(metho->answer, "dos4") == 0)
+    else if (strcasecmp(metho->answer, "dos4") == 0)
 	method = DOS4;
     else
 	method = UNCORRECTED;
@@ -402,7 +388,7 @@
 	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, " LANSAT: %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);



More information about the grass-commit mailing list