[GRASS-SVN] r55376 - grass/branches/releasebranch_6_4/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Mar 14 05:01:19 PDT 2013
Author: neteler
Date: 2013-03-14 05:01:18 -0700 (Thu, 14 Mar 2013)
New Revision: 55376
Modified:
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_set.c
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/local_proto.h
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c
Log:
Jorge Tizado: bugfix for thermal band and metadata file
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html 2013-03-14 12:01:18 UTC (rev 55376)
@@ -8,15 +8,16 @@
<p>
Usually, to do so the production date, the acquisition date, and the
-solar elevation is needed. Moreover, for Landsat-7 ETM+ it is also
+solar elevation are needed. Moreover, for Landsat-7 ETM+ it is also
needed the gain (high or low) of the nine respective bands.
<p>
-Optionally (recommended), the data can be read from metadata file (.met or MTL.txt) for all
-Landsat MSS, TM and ETM+. However, if the solar elevation or the
-product creation date are given the values of the metadata file are
-overwriten. This is necessary when the data in the .met file is
-incorrect or not accurate.
+Optionally (recommended), the data can be read from metadata file
+(.met or MTL.txt) for all Landsat MSS, TM, ETM+ and OLI/TIRS. However,
+if the solar elevation is given the value of the metadata file are
+overwritten. This is necessary when the data in the .met file is
+incorrect or not accurate. Also, if acquisition or production dates are
+not found in the metadata file then the command line values are used.
<p>
<b>Attention</b>: Any null value or smaller than QCALmin in the input
@@ -138,6 +139,7 @@
<li>Landsat-5 MSS: April 6, 1984 and November 9, 1984</li>
<li>Landsat-5 TM: May 4, 2003 and April, 2 2007</li>
<li>Landsat-7 ETM+: July 1, 2000</li>
+ <li>Landsat-8 OLI/TIRS: launched in 2013</li>
</ul>
<h2>EXAMPLES</h2>
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c 2013-03-14 12:01:18 UTC (rev 55376)
@@ -5,9 +5,9 @@
#include "landsat.h"
-#define PI 3.1415926535897932384626433832795
-#define R2D 57.295779513082320877
-#define D2R 0.017453292519943295769
+#define PI M_PI
+#define R2D 180./M_PI
+#define D2R M_PI/180.
/****************************************************************************
* PURPOSE: Calibrated Digital Number to at-satellite Radiance
@@ -22,7 +22,7 @@
*****************************************************************************/
double lsat_rad2ref(double rad, band_data * band)
{
- return (double)(rad / band->K2);
+ return (double)(rad / band->K1);
}
/****************************************************************************
@@ -62,10 +62,9 @@
sin_e = (double)(sin(D2R * lsat->sun_elev));
cos_v = (double)(cos(D2R * (lsat->number < 4 ? 9.2 : 8.2)));
- /** Global irradiance on the sensor.
- Radiance to reflectance coefficient, only NO thermal bands.
- K1 and K2 variables are also utilized as thermal constants
- */
+ /** Global irradiance on the sensor.
+ Radiance to reflectance coefficient, only NO thermal bands.
+ */
if (lsat->band[i].thermal == 0) {
switch (method) {
case DOS2:
@@ -129,16 +128,17 @@
break;
}
rad_sun = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
- G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz,
- Edown);
+ if (method > DOS)
+ G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
+ TAUv, TAUz, Edown);
- lsat->band[i].K1 = 0.;
- lsat->band[i].K2 = rad_sun;
+ lsat->band[i].K1 = rad_sun;
+ lsat->band[i].K2 = 0.;
}
- /** Digital number to radiance coefficients.
- Whitout atmospheric calibration for thermal bands.
- */
+ /** Digital number to radiance coefficients.
+ Without atmospheric calibration for thermal bands.
+ */
lsat->band[i].gain = ((lsat->band[i].lmax - lsat->band[i].lmin) /
(lsat->band[i].qcalmax - lsat->band[i].qcalmin));
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h 2013-03-14 12:01:18 UTC (rev 55376)
@@ -30,14 +30,14 @@
double wavemax, wavemin; /* Wavelength in µm */
+ double esun; /* Mean solar irradiance */
double lmax, lmin; /* Spectral radiance */
double qcalmax, qcalmin; /* Quantized calibrated pixel */
- double esun; /* Mean solar irradiance */
char thermal; /* Flag to thermal band */
double gain, bias; /* Gain and Bias of sensor */
- double K1, K2; /* Thermal calibration constants,
- or Rad2Ref constants */
+ double K1, K2; /* Thermal calibration or
+ Rad2Ref constants */
} band_data;
@@ -48,8 +48,9 @@
char creation[11]; /* Image production date */
char date[11]; /* Image acquisition date */
+
double dist_es; /* Distance Earth-Sun */
- double sun_elev; /* Solar elevation */
+ double sun_elev; /* Sun elevation */
char sensor[10]; /* Type of sensor: MSS, TM, ETM+, OLI/TIRS */
int bands; /* Total number of bands */
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c 2013-03-14 12:01:18 UTC (rev 55376)
@@ -9,44 +9,56 @@
#include "local_proto.h"
#include "earth_sun.h"
+#define PI M_PI
+#define D2R M_PI/180.
+
#define MAX_STR 127
-#define METADATA_SIZE 65535 /* MTL.txt file size 65535 bytes */
-#define TM5_MET_SIZE 28700 /* .met file size 28686 bytes */
+#define METADATA_SIZE 65535 /* MTL.txt file size 65535 bytes */
+#define TM5_MET_SIZE 28700 /* .met file size 28686 bytes */
inline void chrncpy(char *dest, char src[], int n)
{
int i;
- for( i = 0 ; i < n && src[i] != '\0' && src[i] != '\"'; i++) dest[i] = src[i];
+
+ for (i = 0; i < n && src[i] != '\0' && src[i] != '\"'; i++)
+ dest[i] = src[i];
dest[i] = '\0';
}
-void get_oldformat( const char metadata[], char * key, char value[] )
+void get_metformat(const char metadata[], char *key, char value[])
{
int i = 0;
- char * ptr = strstr(metadata, key);
- if (ptr != NULL)
- {
- ptr = strstr(ptr, " VALUE ");
- if (ptr != NULL)
- {
- while (*ptr++ != '=') ;
- while (*ptr <= ' ' || *ptr == '\"') *ptr++;
- while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ') value[i++] = *ptr++;
- }
+ char *ptr = strstr(metadata, key);
+
+ if (ptr != NULL) {
+ ptr = strstr(ptr, " VALUE ");
+ if (ptr != NULL) {
+ while (*ptr++ != '=') ;
+ while (*ptr <= ' ')
+ *ptr++;
+ if (*ptr == '\"')
+ *ptr++;
+ while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
+ value[i++] = *ptr++;
+ }
}
value[i] = '\0';
}
-void get_newformat( const char metadata[], char * key, char value[] )
+void get_mtlformat(const char metadata[], char *key, char value[])
{
int i = 0;
- char * ptr = strstr(metadata, key);
- if (ptr != NULL)
- {
- while (*ptr++ != '=') ;
- while (*ptr <= ' ' || *ptr == '\"') *ptr++;
- while (i < MAX_STR && *ptr != '\"' && *ptr > ' ') value[i++] = *ptr++;
+ char *ptr = strstr(metadata, key);
+
+ if (ptr != NULL) {
+ while (*ptr++ != '=') ;
+ while (*ptr <= ' ')
+ *ptr++;
+ if (*ptr == '\"')
+ *ptr++;
+ while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
+ value[i++] = *ptr++;
}
value[i] = '\0';
}
@@ -56,142 +68,295 @@
* PURPOSE: Read parameters from Landsat metadata files
*****************************************************************************/
-void lsat_metadata( char * metafile, lsat_data * lsat)
+void lsat_metadata(char *metafile, lsat_data * lsat)
{
FILE *f;
char mtldata[METADATA_SIZE];
char key[MAX_STR], value[MAX_STR];
- void (*get_mtldata)( const char [], char *, char [] );
- int i, j, version;
+ void (*get_mtldata) (const char[], char *, char[]);
+ int i, j, ver_mtl;
if ((f = fopen(metafile, "r")) == NULL)
- G_fatal_error(_("Metadata file <%s> not found"), metafile);
+ G_fatal_error(_("Metadata file <%s> not found"), metafile);
i = fread(mtldata, METADATA_SIZE, 1, f);
- get_mtldata = (strstr(mtldata, " VALUE ") != NULL) ? get_oldformat : get_newformat;
- version = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
-
+
+ /* get version of the metadata file */
+ get_mtldata =
+ (strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
+ ver_mtl = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
+
/* Fill with product metadata */
get_mtldata(mtldata, "SPACECRAFT_ID", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
}
-
- i = 0; while( value[i] && (value[i] < '0' || value[i] > '9') ) i++;
+
+ i = 0;
+ while (value[i] && (value[i] < '0' || value[i] > '9'))
+ i++;
lsat->number = atoi(value + i);
get_mtldata(mtldata, "SENSOR_ID", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "SENSORSHORTNAME", value);
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "SENSORSHORTNAME", value);
}
chrncpy(lsat->sensor, value, 8);
get_mtldata(mtldata, "DATE_ACQUIRED", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "ACQUISITION_DATE", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "CALENDARDATE", value);
- }
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "ACQUISITION_DATE", value);
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "CALENDARDATE", value);
+ }
}
- chrncpy(lsat->date, value, 10);
+ if (value[0] != '\0')
+ chrncpy(lsat->date, value, 10);
+ else
+ G_warning("Using adquisition date from the command line 'date'");
get_mtldata(mtldata, "FILE_DATE", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "CREATION_TIME", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "PRODUCTIONDATETIME", value);
- }
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "CREATION_TIME", value);
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "PRODUCTIONDATETIME", value);
+ }
}
- chrncpy(lsat->creation, value, 10);
+ if (value[0] != '\0')
+ chrncpy(lsat->creation, value, 10);
+ else
+ G_warning
+ ("Using production date from the command line 'product_date'");
get_mtldata(mtldata, "SUN_ELEVATION", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "SolarElevation", value);
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "SolarElevation", value);
}
lsat->sun_elev = atof(value);
- if( lsat->sun_elev == 0. )
- G_warning("Sun elevation is %f", lsat->sun_elev);
+ if (lsat->sun_elev == 0.)
+ G_warning("Sun elevation is %f", lsat->sun_elev);
- /* Fill data with the basic sensor parameters */
- switch(lsat->number)
- {
- case 1:
- set_MSS1(lsat);
- break;
- case 2:
- set_MSS2(lsat);
- break;
- case 3:
- set_MSS3(lsat);
- break;
- case 4:
- if (lsat->sensor[0] == 'M')
- set_MSS4(lsat);
- else
- set_TM4(lsat);
- break;
- case 5:
- if (lsat->sensor[0] == 'M')
- set_MSS5(lsat);
- else
- set_TM5(lsat);
- break;
- case 7:
- {
- char * fmt_gain[] = { "BAND%d_GAIN%d ", " GAIN_BAND_%d_VCID_%d" };
- for( i = 1, j = 0; i < 9; i++ )
- {
- sprintf(key, fmt_gain[version], i, 1);
- if (i != 6) key[version == 0 ? 10 : 12] = '\0';
- get_mtldata(mtldata, key, value + j++);
- if (i == 6)
- {
- sprintf(key, fmt_gain[version], i, 2);
- get_mtldata(mtldata, key, value + j++);
- }
- }
- value[j] = '\0';
- G_debug(1, "ETM+ gain = [%s]", value);
- set_ETM(lsat, value);
- break;
- }
- default:
- G_warning("Unable to recognize satellite platform [%d]", lsat->number);
- break;
+ /* Fill the data with the basic/default sensor parameters */
+ switch (lsat->number) {
+ case 1:
+ set_MSS1(lsat);
+ break;
+ case 2:
+ set_MSS2(lsat);
+ break;
+ case 3:
+ set_MSS3(lsat);
+ break;
+ case 4:
+ if (lsat->sensor[0] == 'M')
+ set_MSS4(lsat);
+ else
+ set_TM4(lsat);
+ break;
+ case 5:
+ if (lsat->sensor[0] == 'M')
+ set_MSS5(lsat);
+ else
+ set_TM5(lsat);
+ break;
+ case 7:
+ {
+ char *fmt_gain[] = { "BAND%d_GAIN%d", " GAIN_BAND_%d_VCID_%d" };
+ for (i = 1, j = 0; i < 9; i++) {
+ sprintf(key, fmt_gain[ver_mtl], i, 1);
+ if (i != 6)
+ key[ver_mtl == 0 ? 10 : 12] = '\0';
+ get_mtldata(mtldata, key, value + j++);
+ if (i == 6) {
+ sprintf(key, fmt_gain[ver_mtl], i, 2);
+ get_mtldata(mtldata, key, value + j++);
+ }
+ }
+ value[j] = '\0';
+ G_debug(1, "ETM+ gain = [%s]", value);
+ set_ETM(lsat, value);
+ break;
+ }
+ case 8:
+ set_LDCM(lsat);
+ break;
+
+ default:
+ G_warning("Unable to recognize satellite platform [%d]",
+ lsat->number);
+ break;
}
/* Update the information from metadata file, if infile */
-// if( format == NEW_FORMAT )
- if( get_mtldata == get_newformat )
- {
- char * fmt_lmax[] = { "LMAX_BAND%d", "RADIANCE_MAXIMUM_BAND_%d" };
- char * fmt_lmin[] = { "LMIN_BAND%d", "RADIANCE_MINIMUM_BAND_%d" };
- char * fmt_qcalmax[] = { "QCALMAX_BAND%d", "QUANTIZE_CAL_MAX_BAND_%d" };
- char * fmt_qcalmin[] = { "QCALMIN_BAND%d", "QUANTIZE_CAL_MIN_BAND_%d" };
+ if (get_mtldata == get_mtlformat) {
+ if (ver_mtl == 0) { /* now MTLold.txt */
+ G_verbose_message("Metada file is MTL file: old format");
+ for (i = 0; i < lsat->bands; i++) {
+ sprintf(key, "LMAX_BAND%d", lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].lmax = atof(value);
+ sprintf(key, "LMIN_BAND%d", lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].lmin = atof(value);
+ sprintf(key, "QCALMAX_BAND%d", lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].qcalmax = atof(value);
+ sprintf(key, "QCALMIN_BAND%d", lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].qcalmin = atof(value);
+ }
+ }
+ else { /* now MTL.txt */
- for (i = 0; i < lsat->bands; i++) {
- sprintf(key, fmt_lmax[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].lmax = atof(value);
- sprintf(key, fmt_lmin[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].lmin = atof(value);
- sprintf(key, fmt_qcalmax[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].qcalmax = atof(value);
- sprintf(key, fmt_qcalmin[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].qcalmin = atof(value);
- }
+ G_verbose_message("Metada file is MTL file: new format");
+ 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) {
+ sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
+ lsat->band[i].code - 60);
+ 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_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_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_mtldata(mtldata, key, value);
+ lsat->band[i].qcalmin = atof(value);
+ }
+ else {
+ sprintf(key, "RADIANCE_MAXIMUM_BAND_%d",
+ lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].lmax = atof(value);
+ sprintf(key, "RADIANCE_MINIMUM_BAND_%d",
+ lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].lmin = atof(value);
+ sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d",
+ lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].qcalmax = atof(value);
+ sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d",
+ lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].qcalmin = atof(value);
+ }
+ }
+ }
+ else {
+ G_verbose_message
+ ("RADIANCE & QUANTIZE from radiometric rescaling group of the metadata file");
+ /* from LDCM sample file: mode = 0; from LDCM-DFCB-004.pdf file: mode = 1 */
+ int mode =
+ (strstr(mtldata, "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_mtldata(mtldata, key, value);
+ lsat->band[i].gain = atof(value);
+ sprintf(key, fmt_radad[mode], lsat->band[i].code);
+ 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.;
+ lsat->band[i].lmin =
+ lsat->band[i].gain * lsat->band[i].qcalmin +
+ lsat->band[i].bias;
+ lsat->band[i].lmax =
+ lsat->band[i].gain * lsat->band[i].qcalmax +
+ lsat->band[i].bias;
+ /* ----- */
+ if (lsat->number == 8) {
+ if (lsat->band[i].thermal) {
+ sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+ get_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 {
+ lsat->band[i].K1 = 0.;
+ /* ESUN from metadafa file: bias/K2 seem better than gain/K1 */
+ sprintf(key, fmt_refad[mode], lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].K2 = atof(value);
+ lsat->band[i].esun =
+ (double)(PI * lsat->dist_es * lsat->dist_es *
+ lsat->band[i].bias) / (sin(D2R *
+ lsat->
+ sun_elev)
+ *
+ lsat->band[i].
+ K2);
+ }
+ }
+ }
+ }
+ /* Other possible values in file */
+ get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value);
+ if (value[0] != '\0') {
+ lsat->dist_es = atof(value);
+ G_verbose_message(1,
+ "ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
+ }
+ }
}
+ else {
+ G_verbose_message("Metada file is MET file");
+ G_verbose_message
+ ("RADIANCE & QUANTIZE from band setting of the metadata file");
+ for (i = 0; i < lsat->bands; i++) {
+ sprintf(key, "Band%dGainSetting", lsat->band[i].code);
+ get_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_mtldata(mtldata, key, value);
+ if (value[0] == '\0') {
+ G_warning(key);
+ continue;
+ }
+ lsat->band[i].bias = atof(value);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 1.;
+ lsat->band[i].lmin =
+ lsat->band[i].gain * lsat->band[i].qcalmin +
+ lsat->band[i].bias;
+ lsat->band[i].lmax =
+ lsat->band[i].gain * lsat->band[i].qcalmax +
+ lsat->band[i].bias;
+ }
+ }
+
(void)fclose(f);
return;
}
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_set.c 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_set.c 2013-03-14 12:01:18 UTC (rev 55376)
@@ -15,10 +15,10 @@
/* green, red, near infrared, near infrared */
int band[] = { 1, 2, 3, 4 };
- int code[] = { 4, 5, 6, 7 };
+ int code[] = { 4, 5, 6, 7 }; /* corrected for MSS4 y MSS5 to 1,2,3,4 */
double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
- /* original: 79x57, now all 60 */
+ /* original: 79x57, now all 60 m */
strcpy(lsat->sensor, "MSS");
@@ -46,7 +46,7 @@
/* 30, 30, 30, 30, 30, 120 original, 60 resamples before Feb 25, 2010 and 30 after, 30 */
if (!lsat->sensor)
- strcpy(lsat->sensor, "TM");
+ strcpy(lsat->sensor, "TM");
lsat->bands = 7;
for (i = 0; i < lsat->bands; i++) {
@@ -87,7 +87,7 @@
return;
}
-void sensor_OLI(lsat_data * lsat)
+void sensor_LDCM(lsat_data * lsat)
{
int i;
@@ -95,21 +95,25 @@
* 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");
lsat->bands = 11;
for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(code + i);
- lsat->band[i].wavemax = *(wmax + i);
- lsat->band[i].wavemin = *(wmin + i);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 1.;
- lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = *(code + i);
+ lsat->band[i].wavemax = *(wmax + i);
+ lsat->band[i].wavemin = *(wmin + i);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 1.;
+ lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
}
return;
}
@@ -180,9 +184,9 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("1975-07-16"))
- i = 0;
+ i = 0;
else
- i = 1;
+ i = 1;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -279,11 +283,11 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("1982-08-26"))
- i = 0;
+ i = 0;
else if (julian < julian_char("1983-03-31"))
- i = 1;
+ i = 1;
else
- i = 2;
+ i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -329,11 +333,11 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("1983-08-01"))
- i = 0;
+ i = 0;
else if (julian < julian_char("1984-01-15"))
- i = 1;
+ i = 1;
else
- i = 2;
+ i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -436,11 +440,11 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("2003-05-04"))
- i = 0;
+ i = 0;
else if (julian < julian_char("2007-04-02"))
- i = 1;
+ i = 1;
else
- i = 2;
+ i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -453,9 +457,9 @@
}
jbuf = julian_char("2004-04-04");
- if (julian >= jbuf && !(lsat->flag & METADATAFILE) )
- {
- G_warning("Using QCalMin=1.0 as a NLAPS product processed after 04/04/2004");
+ if (julian >= jbuf && !(lsat->flag & METADATAFILE)) {
+ G_warning
+ ("Using QCalMin=1.0 as a NLAPS product processed after 04/04/2004");
}
lsat->number = 5;
sensor_TM(lsat);
@@ -463,9 +467,9 @@
lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
if (julian >= jbuf)
lsat->band[i].qcalmin = 1.;
+ j = lsat->band[i].number - 1;
lsat->band[i].esun = *(esun + j);
lsat->band[i].lmax = *(lmax + j);
lsat->band[i].lmin = *(lmin + j);
@@ -528,7 +532,6 @@
for (i = 0; i < lsat->bands; i++) {
j = lsat->band[i].number - 1;
- lsat->band[i].esun = *(esun + j);
if (gain[i] == 'H' || gain[i] == 'h') {
lmax = LmaxH[k];
lmin = LminH[k];
@@ -537,6 +540,7 @@
lmax = LmaxL[k];
lmin = LminL[k];
}
+ lsat->band[i].esun = *(esun + j);
lsat->band[i].lmax = *(lmax + j);
lsat->band[i].lmin = *(lmin + j);
if (lsat->band[i].thermal) {
@@ -552,42 +556,39 @@
* PURPOSE: Store values of Landsat-8 OLI/TIRS
* February 14, 2013
*****************************************************************************/
-void set_OLI(lsat_data * lsat)
+void set_LDCM(lsat_data * lsat)
{
int i, j;
double julian, *lmax, *lmin;
- /* Spectral radiances at detector
- double Lmax[][4] = {
- {220., 175., 145., 147.},
- {259., 179., 149., 128.}
+ /* Spectral radiances at detector */
+ double Lmax[][11] = {
+ {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.} // <<<<<<<<<<<<<< valores incorrectos
};
- double Lmin[][4] = {
- {4., 3., 3., 1.},
- {4., 3., 3., 1.}
+ double Lmin[][11] = {
+ {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.} // <<<<<<<<<<<<<< valores incorrectos
};
- * Solar exoatmospheric spectral irradiances
- double esun[] = { 1824., 1570., 1249., 853.4 };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 2062., 21931., 1990., 1688., 1037., 268.6, 94.6, 1892., 399.0, 0., 0. }; // <<<<<<<<<<<<<< estimados
- lmax = Lmax[i];
- lmin = Lmin[i];
+ lmax = Lmax[0];
+ lmin = Lmin[0];
- lsat->number = 3;
- sensor_MSS(lsat);
+ lsat->number = 8;
+ sensor_LDCM(lsat);
lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
- lsat->band[i].esun = *(esun + j);
- lsat->band[i].lmax = *(lmax + j);
- lsat->band[i].lmin = *(lmin + j);
- if (lsat->band[i].thermal) {
- lsat->band[i].K1 = 0;
- lsat->band[i].K2 = 0;
- }
+ j = lsat->band[i].number - 1;
+ lsat->band[i].esun = *(esun + j);
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
+ if (lsat->band[i].thermal) {
+ lsat->band[i].K1 = (lsat->band[i].number == 10 ? 666.09 : 697.56);
+ lsat->band[i].K2 = 1282.71;
+ }
}
- */
- G_debug(1, "Landsat-8 LDCM");
+ G_debug(1, "Landsat-8 OLI/TIRS");
return;
}
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/local_proto.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/local_proto.h 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/local_proto.h 2013-03-14 12:01:18 UTC (rev 55376)
@@ -17,4 +17,7 @@
void set_ETM(lsat_data *, char[]);
+void set_LDCM(lsat_data *);
+
+
#endif
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c 2013-03-14 10:41:43 UTC (rev 55375)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c 2013-03-14 12:01:18 UTC (rev 55376)
@@ -8,9 +8,9 @@
* 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+
+ * for Landsat 1/2/3/4/5 MS, 4/5 TM, 7 ETM+, and 8 OLI/TIRS
*
- * COPYRIGHT: (C) 2006-2012 by the GRASS Development Team
+ * COPYRIGHT: (C) 2006-2013 by the GRASS Development Team
*
* This program is free software under the GNU General
* Public License (>=v2). Read the file COPYING that
@@ -40,10 +40,10 @@
RASTER_MAP_TYPE in_data_type;
- struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate, *pdate, *elev,
- *bgain, *metho, *perc, *dark, *atmo;
+ struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate,
+ *pdate, *elev, *bgain, *metho, *perc, *dark, *atmo;
char *inputname, *met, *outputname, *sensorname;
- struct Flag *frad;
+ struct Flag *frad, *named;
lsat_data lsat;
char band_in[GNAME_MAX], band_out[GNAME_MAX];
@@ -62,7 +62,8 @@
module = G_define_module();
module->description =
_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+.");
- module->keywords = _("imagery, landsat, top-of-atmosphere reflectance, dos-type simple atmospheric correction");
+ module->keywords =
+ _("imagery, landsat, top-of-atmosphere reflectance, dos-type simple atmospheric correction");
/* It defines the different parameters */
input_prefix = G_define_option();
@@ -75,7 +76,8 @@
output_prefix = G_define_option();
output_prefix->key = "output_prefix";
output_prefix->label = _("Prefix for output raster maps");
- output_prefix->description = _("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
+ output_prefix->description =
+ _("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
output_prefix->type = TYPE_STRING;
output_prefix->required = YES;
@@ -89,18 +91,20 @@
sensor->key = "sensor";
sensor->type = TYPE_STRING;
sensor->label = _("Spacecraft sensor");
- sensor->description = _("Required only if 'metfile' not given");
- sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7";
+ sensor->description =
+ _("Required only if 'metfile' not given (recommended by sanity)");
+ sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,ot8";
sensor->descriptions =
- _("mss1;Landsat-1 MSS;"
- "mss2;Landsat-2 MSS;"
- "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;
+ _("mss1;Landsat_1 MSS;"
+ "mss2;Landsat_2 MSS;"
+ "mss3;Landsat_3 MSS;"
+ "mss4;Landsat_4 MSS;"
+ "mss5;Landsat_5 MSS;"
+ "tm4;Landsat_4 TM;"
+ "tm5;Landsat_5 TM;"
+ "tm7;Landsat_7 ETM+;"
+ "ot8;Landsat_8 OLI/TIRS");
+ sensor->required = NO; /* perhaps YES for clarity */
sensor->guisection = _("Metadata");
metho = G_define_option();
@@ -123,10 +127,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");
@@ -143,7 +147,7 @@
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");
@@ -151,7 +155,8 @@
perc->key = "percent";
perc->type = TYPE_DOUBLE;
perc->required = NO;
- perc->description = _("Percent of solar radiance in path radiance");
+ perc->label = _("Percent of solar radiance in path radiance");
+ perc->description = _("Required only if 'method' is any DOS");
perc->answer = "0.01";
perc->guisection = _("Settings");
@@ -159,8 +164,9 @@
dark->key = "pixel";
dark->type = TYPE_INTEGER;
dark->required = NO;
- dark->description =
+ dark->label =
_("Minimum pixels to consider digital number as dark object");
+ dark->description = _("Required only if 'method' is any DOS");
dark->answer = "1000";
dark->guisection = _("Settings");
@@ -169,14 +175,21 @@
atmo->type = TYPE_DOUBLE;
atmo->required = NO;
atmo->description = _("Rayleigh atmosphere (diffuse sky irradiance)"); /* scattering coefficient? */
+ atmo->description = _("Required only if 'method' is DOS3");
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");
+ frad->description =
+ _("Output at-sensor radiance instead reflectance for all bands");
+ named = G_define_flag();
+ named->key = 'n';
+ named->description =
+ _("Input raster maps use as extension the number of the band instead the code");
+
/* options and afters parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -189,7 +202,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));
@@ -210,57 +223,62 @@
}
lsat.sun_elev = elev->answer == NULL ? 0. : atof(elev->answer);
+
percent = atof(perc->answer);
pixel = atoi(dark->answer);
rayleigh = atof(atmo->answer);
/* Data from metadata file */
- /* Unnecessary because G_zero filled, but sanity */
+ /* Unnecessary because G_zero filled, but by sanity */
lsat.flag = NOMETADATAFILE;
- if (met != NULL)
- {
- lsat.flag = METADATAFILE;
- lsat_metadata( 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"));
+ if (met != NULL) {
+ lsat.flag = METADATAFILE;
+ lsat_metadata(met, &lsat);
+ G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number,
+ lsat.sensor);
- G_debug(1, "Landsat-%d %s with data set in met file [%s]", lsat.number, lsat.sensor, met);
+ if (!lsat.sensor || lsat.number > 8 || lsat.number < 1)
+ G_fatal_error(_("Failed to identify satellite"));
- /* Overwrite solar elevation of metadata file */
- if (elev->answer != NULL)
- lsat.sun_elev = atof(elev->answer);
+ 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) {
- G_fatal_error(_("Lacking '%s' or '%s' for this satellite"),
+ G_fatal_error(_("Lacking '%s' and/or '%s' for this satellite"),
adate->key, elev->key);
}
else {
- /* Need gain */
+ /* 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"));
- set_ETM(&lsat, bgain->answer);
- }
- /* Not need gain */
- else 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);
+ if (bgain->answer == NULL || strlen(bgain->answer) != 9)
+ G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) data"));
+ set_ETM(&lsat, bgain->answer);
+ }
+ /* Not need gain */
+ else if (strcmp(sensorname, "ot8") == 0)
+ set_LDCM(&lsat);
+ else if (strcmp(sensorname, "tm5") == 0)
+ set_TM5(&lsat);
+ else if (strcmp(sensorname, "tm4") == 0)
+ 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);
}
/*****************************************
@@ -307,7 +325,8 @@
if ((infd = G_open_cell_old(band_in, "")) < 0)
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
- G_fatal_error(_("Unable to read header of raster map <%s>"), band_in);
+ G_fatal_error(_("Unable to read header of raster map <%s>"),
+ band_in);
G_set_window(&cellhd);
in_data_type = G_raster_map_type(band_in, mapset);
@@ -341,7 +360,7 @@
}
/* DN of dark object */
for (j = lsat.band[i].qcalmin; j < 256; j++) {
- if (hist[j] >= (unsigned int) pixel) {
+ if (hist[j] >= (unsigned int)pixel) {
dn_dark[i] = j;
break;
}
@@ -368,16 +387,19 @@
lsat_bandctes(&lsat, i, method, percent, dn_dark[i], rayleigh);
}
- 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, lsat.sensor);
+ fprintf(stderr, " LANDSAT: %d SENSOR: %s\n", lsat.number,
+ lsat.sensor);
fprintf(stderr, " ACQUISITION DATE %s [production date %s]\n",
lsat.date, lsat.creation);
fprintf(stderr, " earth-sun distance = %.8lf\n", lsat.dist_es);
fprintf(stderr, " solar elevation angle = %.8lf\n", lsat.sun_elev);
- fprintf(stderr, " Method of calculus = %s\n",
+ fprintf(stderr, " Atmospheric correction: %s\n",
(method == CORRECTED ? "CORRECTED"
: (method == UNCORRECTED ? "UNCORRECTED" : metho->answer)));
if (method > DOS) {
@@ -387,7 +409,7 @@
}
for (i = 0; i < lsat.bands; i++) {
fprintf(stderr, "-------------------\n");
- fprintf(stderr, " BAND %d %s (code %d)\n",
+ fprintf(stderr, " BAND %d %s(code %d)\n",
lsat.band[i].number,
(lsat.band[i].thermal ? "thermal " : ""),
lsat.band[i].code);
@@ -396,12 +418,12 @@
lsat.band[i].qcalmin, lsat.band[i].qcalmax);
fprintf(stderr, " calibration constants (L): %.3lf to %.3lf\n",
lsat.band[i].lmin, lsat.band[i].lmax);
- fprintf(stderr, " at-%s radiance = %.5lf * DN + %.5lf\n",
+ fprintf(stderr, " at-%s radiance = %.8lf * DN + %.3lf\n",
(method > DOS ? "surface" : "sensor"), lsat.band[i].gain,
lsat.band[i].bias);
if (lsat.band[i].thermal) {
fprintf(stderr,
- " at-sensor temperature = %.3lf / log[(%.3lf / radiance) + 1.0]\n",
+ " at-sensor temperature = %.5lf / log[(%.5lf / radiance) + 1.0]\n",
lsat.band[i].K2, lsat.band[i].K1);
}
else {
@@ -410,12 +432,12 @@
lsat.band[i].esun);
fprintf(stderr, " at-%s reflectance = radiance / %.5lf\n",
(method > DOS ? "surface" : "sensor"),
- lsat.band[i].K2);
+ lsat.band[i].K1);
if (method > DOS) {
fprintf(stderr,
" the darkness DN with a least %d pixels is %d\n",
pixel, dn_dark[i]);
- fprintf(stderr, " the mode of DN is %d\n", dn_mode[i]);
+ fprintf(stderr, " the DN mode is %d\n", dn_mode[i]);
}
}
}
@@ -423,13 +445,14 @@
fflush(stderr);
}
- /*****************************************
- * ------------ CALCULUS -----------------
- *****************************************/
+ /*****************************************
+ * ------------ CALCULUS -----------------
+ *****************************************/
G_message(_("Calculating..."));
for (i = 0; i < lsat.bands; i++) {
- sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
+ sprintf(band_in, "%s%d", inputname,
+ (named->answer ? lsat.band[i].number : lsat.band[i].code));
sprintf(band_out, "%s%d", outputname, lsat.band[i].code);
mapset = G_find_cell2(band_in, "");
@@ -442,7 +465,8 @@
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
in_data_type = G_raster_map_type(band_in, mapset);
if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
- G_fatal_error(_("Unable to read header of raster map <%s>"), band_in);
+ G_fatal_error(_("Unable to read header of raster map <%s>"),
+ band_in);
/* set same size as original band raster */
G_set_window(&cellhd);
@@ -462,10 +486,10 @@
ncols = G_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);
@@ -538,8 +562,11 @@
/* Initialize the 'hist' structure with basic info */
G_short_history(band_out, "raster", &history);
sprintf(history.edhist[0], " %s of Landsat-%d %s (method %s)",
- lsat.band[i].thermal ? "Temperature" : "Reflectance",
- lsat.number, lsat.sensor, metho->answer);
+ (frad->
+ answer ? "Radiance" : (lsat.band[i].
+ thermal ? "Temperature" :
+ "Reflectance")), lsat.number,
+ lsat.sensor, metho->answer);
sprintf(history.edhist[1],
"----------------------------------------------------------------");
sprintf(history.edhist[2],
@@ -571,7 +598,7 @@
lsat.band[i].esun);
sprintf(history.edhist[9],
" Reflectance = Radiance divided by ..... %.5lf",
- lsat.band[i].K2);
+ lsat.band[i].K1);
history.edlinecnt = 10;
if (method > DOS) {
sprintf(history.edhist[10], " ");
@@ -593,7 +620,11 @@
if (lsat.band[i].thermal)
G_write_raster_units(band_out, "Kelvin");
- /* else units = ...? */
+ else if (frad->answer)
+ G_write_raster_units(band_out, "W/(m^2 sr um)");
+ else
+ G_write_raster_units(band_out, "unitless");
+
/* set raster timestamp from acq date? (see r.timestamp module) */
}
More information about the grass-commit
mailing list