[GRASS-SVN] r58097 - grass/trunk/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Oct 23 14:31:54 PDT 2013
Author: neteler
Date: 2013-10-23 14:31:53 -0700 (Wed, 23 Oct 2013)
New Revision: 58097
Modified:
grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html
grass/trunk/imagery/i.landsat.toar/landsat.c
grass/trunk/imagery/i.landsat.toar/landsat_met.c
grass/trunk/imagery/i.landsat.toar/main.c
Log:
i.landsat.toar: fixes for Landsat-8 metadata file support (author: E. Jorge Tizado); sync with r58093
Modified: grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html
===================================================================
--- grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html 2013-10-23 21:30:20 UTC (rev 58096)
+++ grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html 2013-10-23 21:31:53 UTC (rev 58097)
@@ -97,13 +97,66 @@
TAUz = exp[-t/sin(e)], Esky = PI * radiance_dark </li>
</ul>
-<b>Attention</b>: Output radiance remain untouched (i.e. no set to
-0. when it is negative) then they are possible negative
-values. However, output reflectance is set to 0. when is obtained a
-negative value.
+<b>Attention</b>: Output radiance remain untouched (i.e. no set to 0.0
+when it is negative) then they are possible negative values. However,
+output reflectance is set to 0.0 when is obtained a negative value.
<h2>NOTES</h2>
+<h3>On Landsat-8 metadata file </h3>
+
+NASA reports a structure of the L1G Metadata file
+(<a href="http://landsat.usgs.gov/documents/LDCM-DFCB-004.pdf">LDCM-DFCB-004.pdf</a>)
+for Landsat Data Continuity Mission (i.e. Landsat-8).
+
+<p>
+NASA retains in MIN_MAX_RADIANCE group the necessary information
+to transform Digital Numbers (DN) in radiance values. Then,
+<em>i.landsat.toar</em> replaces the possible standard values with the
+metadata values. The results match with the values reported by the
+metada file in RADIOMETRIC_RESCALING group.
+
+<p>
+Also, NASA reports the same values of reflectance for all bands
+in max-min values and in gain-bias values. This is strange that all
+bands have the same range of reflectance. Also, they wrote in the
+web page as to calculate reflectance directly from DN, first with
+RADIOMETRIC_RESCALING values and second divided by sin(sun_elevation).
+
+<p>
+This is a simple rescaling
+
+<ul>
+ <li> reflectance = radiance / sun_radiance = (DN * RADIANCE_MULT + RADIANCE_ADD) / sun_radiance</li>
+ <li> now reflectance = DN * REFLECTANCE_MULT + REFLECTANCE_ADD </li>
+ <li> then REFLECTANCE_MULT = RADIANCE_MULT / sun_radiance </li>
+ <li> and REFLECTANCE_ADD = RADIANCE_ADD / sun_radiance </li>
+</ul>
+
+<p>
+The problem arises when we need ESUN values (not provided) to
+compute sun_radiance and DOS. We assume that REFLECTANCE_MAXIMUM
+corresponds to the RADIANCE_MAXIMUM, then
+
+<ul>
+ <li> REFLECTANCE_MAXIMUM / sin(e) = RADIANCE_MAXIMUM / sun_radiance</li>
+ <li> Esun = (PI * d^2) * RADIANCE_MAXIMUM / REFLECTANCE_MAXIMUM </li>
+</ul>
+
+where <em>d</em> is the earth-sun distance provided by metadata file
+or computed inside the program.
+
+<p>
+The <em>i.landsat.toar</em> reverts back the NASA rescaling to continue
+using Lmax, Lmin, and Esun values to compute the constant to convert
+DN to radiance and radiance to reflectance with the "traditional"
+equations and simple atmospheric corrections.
+
+<b>Attention</b>: When MAXIMUM values are not provided,
+<em>i.landsat.toar</em> tries to calculate Lmax, Lmin, and Esun from
+RADIOMETRIC_RESCALING (in tests the results were the same).
+
+<h3>Calibration constants</h3>
In verbose mode (flag <b>--verbose</b>), the program write basic
satellite data and the parameters used in the transformations.
@@ -150,6 +203,13 @@
sun_elevation=64.3242970 gain="HHHLHLHHL"
</pre></div>
+or
+
+<div class="code"><pre>
+i.landsat.toar input_prefix=LC80160352013134LGN03_B output_prefix=toar \
+ metfile=LC80160352013134LGN03_MTL.txt sensor=oli8 date=2013-05-14
+</pre></div>
+
<h2>REFERENCES</h2>
<ul>
Modified: grass/trunk/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat.c 2013-10-23 21:30:20 UTC (rev 58096)
+++ grass/trunk/imagery/i.landsat.toar/landsat.c 2013-10-23 21:31:53 UTC (rev 58097)
@@ -87,8 +87,7 @@
double t;
t = 2. / (lsat->band[i].wavemax + lsat->band[i].wavemin);
- t = 0.008569 * t * t * t * t * (1 + 0.0113 * t * t +
- 0.000013 * t * t * t * t);
+ t = 0.008569 * t * t * t * t * (1 + 0.0113 * t * t + 0.000013 * t * t * t * t);
TAUv = exp(-t / cos_v);
TAUz = exp(-t / sin_e);
Edown = rayleigh;
@@ -110,13 +109,9 @@
do {
TAUz = Tz;
TAUv = Tv;
- Lp = Ro -
- percent * TAUv * (lsat->band[i].esun * sin_e * TAUz +
- PI * Lp) / pi_d2;
- Tz = 1. -
- (4. * pi_d2 * Lp) / (lsat->band[i].esun * sin_e);
+ Lp = Ro - percent * TAUv * (lsat->band[i].esun * sin_e * TAUz + PI * Lp) / pi_d2;
+ Tz = 1. - (4. * pi_d2 * Lp) / (lsat->band[i].esun * sin_e);
Tv = exp(sin_e * log(Tz) / cos_v);
- /* G_message("TAUv = %.5f (%.5f), TAUz = %.5f (%.5f) and Edown = %.5f\n", TAUv, Tv, TAUz, Tz, PI * Lp ); */
} while (TAUv != Tv && TAUz != Tz);
TAUz = (Tz < 1. ? Tz : 1.);
TAUv = (Tv < 1. ? Tv : 1.);
@@ -130,42 +125,27 @@
break;
}
lsat->band[i].K2 = 0.;
- lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2; /* rad_sun */
+ lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
if (method > DOS)
- G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
- TAUv, TAUz, Edown);
+ G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz, Edown);
}
/** 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);
+ lsat->band[i].gain = (lsat->band[i].lmax - lsat->band[i].lmin) / (lsat->band[i].qcalmax - lsat->band[i].qcalmin);
if (method == UNCORRECTED || lsat->band[i].thermal) {
/* L = G * (DN - Qmin) + Lmin
* -> bias = Lmin - G * Qmin
*/
- lsat->band[i].bias =
- (lsat->band[i].lmin - lsat->band[i].gain * lsat->band[i].qcalmin);
+ lsat->band[i].bias = (lsat->band[i].lmin - lsat->band[i].gain * lsat->band[i].qcalmin);
}
- /*
- else {
- if (method == CORRECTED) {
- // L = G * (DN - Qmin) + Lmin - Lmin
- -> bias = - G * Qmin *
- lsat->band[i].bias =
- -(lsat->band[i].gain * lsat->band[i].qcalmin);
- // Another possibility is cut when rad < 0 *
- }
- */
else if (method > DOS) {
/* L = Lsat - Lpath = G * DNsat + B - (G * + B - p * rad_sun)
* = G * DNsat - G * + p * rad_sun
* -> bias = p * rad_sun - G
*/
- lsat->band[i].bias =
- percent * lsat->band[i].K1 - lsat->band[i].gain * dark;
+ lsat->band[i].bias = percent * lsat->band[i].K1 - lsat->band[i].gain * dark;
}
}
Modified: grass/trunk/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat_met.c 2013-10-23 21:30:20 UTC (rev 58096)
+++ grass/trunk/imagery/i.landsat.toar/landsat_met.c 2013-10-23 21:31:53 UTC (rev 58097)
@@ -35,6 +35,7 @@
}
}
+/* OLD Metadata Files */
void get_metformat(const char metadata[], char *key, char value[])
{
int i = 0;
@@ -56,6 +57,17 @@
value[i] = '\0';
}
+double get_metdouble(const char metadata[], char *format, int code, char value[])
+{
+ char key[MAX_STR];
+
+ sprintf(key, format, code);
+ get_metformat(metadata, key, value);
+ return atof(value);
+}
+
+
+/* NEW Metadata Files */
void get_mtlformat(const char metadata[], char *key, char value[])
{
int i = 0;
@@ -73,7 +85,17 @@
value[i] = '\0';
}
+double get_mtldouble(const char metadata[], char *format, int code, char value[])
+{
+ char key[MAX_STR];
+
+ sprintf(key, format, code);
+ get_mtlformat(metadata, key, value);
+ return atof(value);
+}
+
+
/****************************************************************************
* PURPOSE: Read parameters from Landsat metadata files
*****************************************************************************/
@@ -84,7 +106,9 @@
char mtldata[METADATA_SIZE];
char key[MAX_STR], value[MAX_STR];
void (*get_mtldata) (const char[], char *, char[]);
+ void (*get_mtlreal) (const char[], char *, int, char[]);
int i, j, ver_mtl;
+ double X2;
/* store metadata in ram */
if ((f = fopen(metafile, "r")) == NULL)
@@ -93,8 +117,17 @@
(void)fclose(f);
/* set version of the metadata file */
- get_mtldata =
- (strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
+ /* get_mtldata = (strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat; */
+ if (strstr(mtldata, " VALUE ") != NULL)
+ {
+ get_mtldata = get_metformat;
+ get_mtlreal = get_metdouble;
+ }
+ else
+ {
+ get_mtldata = get_mtlformat;
+ get_mtlreal = get_mtldouble;
+ }
ver_mtl = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
/* Fill with product metadata */
@@ -238,120 +271,93 @@
else {
G_verbose_message("Metada file is MTL file: new format");
/* Other possible values in the metadata file */
- get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value); /* used after in calculus after */
+ get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value); /* Necessary after */
if (value[0] != '\0')
lsat->dist_es = atof(value);
- /* ----- */
- char *fmt_radmu[] =
- { "RADIANCE_MULTIPLICATIVE_FACTOR_BAND%d",
- "RADIANCE_MULT_BAND_%d" };
- char *fmt_radad[] =
- { "RADIANCE_ADDITIVE_FACTOR_BAND%d", "RADIANCE_ADD_BAND_%d" };
- char *fmt_k1cte[] =
- { "K1_CONSTANT_BAND%d", "K1_CONSTANT_BAND_%d" };
- char *fmt_k2cte[] =
- { "K2_CONSTANT_BAND%d", "K2_CONSTANT_BAND_%d" };
- char *fmt_refmu[] =
- { "REFLECTANCE_MULTIPLICATIVE_FACTOR_BAND%d",
- "REFLECTANCE_MULT_BAND_%d" };
- char *fmt_refad[] =
- { "REFLECTANCE_ADDITIVE_FACTOR_BAND%d",
- "REFLECTANCE_ADD_BAND_%d" };
- /* --NASA-- LDCM sample file: mode = 0; LDCM-DFCB-004.pdf file: mode = 1 */
- int mode =
- (strstr(mtldata, "RADIANCE_MULTIPLICATIVE_FACTOR_BAND") !=
- NULL) ? 0 : 1;
- /* ----- */
- if (strstr(mtldata, "RADIANCE_MAXIMUM_BAND") != NULL) {
- G_verbose_message
- ("RADIANCE & QUANTIZE from data of the metadata file");
+ if (strstr(mtldata, "RADIANCE_MAXIMUM_BAND") != NULL)
+ {
+ G_verbose_message("RADIANCE & QUANTIZE from MIN_MAX_(RADIANCE|PIXEL_VALUE)");
for (i = 0; i < lsat->bands; i++) {
if (lsat->number == 7 && lsat->band[i].thermal) {
- sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
- lsat->band[i].code - 60);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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);
+ sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].qcalmin = atof(value);
}
/* other possible values of each band */
if (lsat->band[i].thermal) {
- sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+ sprintf(key, "K1_CONSTANT_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].K1 = atof(value);
- sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+ sprintf(key, "K2_CONSTANT_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].K2 = atof(value);
}
+ else if (lsat->number == 8)
+ {
+ /* ESUN from REFLECTANCE and RADIANCE ADD_BAND */
+ sprintf(key, "REFLECTANCE_MAXIMUM_BAND_%d", lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ X2 = atof(value);
+ lsat->band[i].esun = (double)(PI * lsat->dist_es * lsat->dist_es * lsat->band[i].lmax) / X2;
+ }
}
+ if (lsat->number == 8)
+ G_warning("ESUN evaluated from REFLECTANCE_MAXIMUM_BAND");
}
else {
- G_verbose_message
- ("RADIANCE & QUANTIZE from radiometric rescaling group of the metadata file");
+ G_verbose_message("RADIANCE & QUANTIZE from RADIOMETRIC_RESCALING");
for (i = 0; i < lsat->bands; i++) {
- sprintf(key, fmt_radmu[mode], lsat->band[i].code);
+ sprintf(key, "RADIANCE_MULT_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].gain = atof(value);
- sprintf(key, fmt_radad[mode], lsat->band[i].code);
+ sprintf(key, "RADIANCE_ADD_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].bias = atof(value);
/* reversing to calculate the values of Lmin and Lmax ... */
- lsat->band[i].lmin =
- lsat->band[i].gain * lsat->band[i].qcalmin +
- lsat->band[i].bias;
- lsat->band[i].lmax =
- lsat->band[i].gain * lsat->band[i].qcalmax +
- lsat->band[i].bias;
+ lsat->band[i].lmin = lsat->band[i].gain * lsat->band[i].qcalmin + lsat->band[i].bias;
+ lsat->band[i].lmax = lsat->band[i].gain * lsat->band[i].qcalmax + lsat->band[i].bias;
/* ... qcalmax and qcalmin loaded in previous sensor_ function */
if (lsat->number == 8) {
if (lsat->band[i].thermal) {
- sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+ sprintf(key, "K1_CONSTANT_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].K1 = atof(value);
- sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+ sprintf(key, "K2_CONSTANT_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].K2 = atof(value);
}
else {
- sprintf(key, fmt_refmu[mode], lsat->band[i].code);
+ sprintf(key, "REFLECTANCE_MULT_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
lsat->band[i].K1 = atof(value);
- sprintf(key, fmt_refad[mode], lsat->band[i].code);
+ sprintf(key, "REFLECTANCE_ADD_BAND_%d", lsat->band[i].code);
get_mtldata(mtldata, key, value);
- lsat->band[i].K2 = atof(value);
- /* ESUN evaluates from metadafa file */
- lsat->band[i].esun =
- (double)(PI * lsat->dist_es * lsat->dist_es *
- lsat->band[i].bias) /
- lsat->band[i].K2;
+ lsat->band[i].K2 = atof(value);
+ /* ESUN from REFLECTANCE_ADD_BAND */
+ lsat->band[i].esun = (double)(PI * lsat->dist_es * lsat->dist_es * lsat->band[i].bias) / lsat->band[i].K2;
/*
double esun1 = (double) (PI * lsat->dist_es * lsat->dist_es * lsat->band[i].bias) / lsat->band[i].K2;
double esun2 = (double) (PI * lsat->dist_es * lsat->dist_es * lsat->band[i].gain) / lsat->band[i].K1;
@@ -360,8 +366,7 @@
}
}
}
- G_warning
- ("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
+ G_warning("ESUN evaluated from REFLECTANCE_ADDITIVE_FACTOR_BAND");
}
}
}
Modified: grass/trunk/imagery/i.landsat.toar/main.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/main.c 2013-10-23 21:30:20 UTC (rev 58096)
+++ grass/trunk/imagery/i.landsat.toar/main.c 2013-10-23 21:31:53 UTC (rev 58097)
@@ -224,7 +224,6 @@
named->description =
_("Input raster maps use as extension the number of the band instead the code");
- /* define the different flags */
print_meta = G_define_flag();
print_meta->key = 'p';
print_meta->description = _("Print output metadata info");
@@ -429,10 +428,9 @@
ptr = (void *)((DCELL *) inrast + col);
q = (int)*((DCELL *) ptr);
break;
- default:
- ptr = NULL;
- q = -1.;
- break;
+ default:
+ ptr = NULL;
+ q = -1.;
}
if (!Rast_is_null_value(ptr, in_data_type) &&
q >= lsat.band[i].qcalmin &&
@@ -599,7 +597,6 @@
default:
ptr = NULL;
qcal = -1.;
- break;
}
if (Rast_is_null_value(ptr, in_data_type) ||
qcal < lsat.band[i].qcalmin) {
@@ -635,8 +632,8 @@
G_free(inrast);
- G_free(outrast);
Rast_close(infd);
+ G_free(outrast);
Rast_close(outfd);
More information about the grass-commit
mailing list