[GRASS-SVN] r58093 - grass/branches/releasebranch_6_4/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Oct 23 14:16:40 PDT 2013
Author: neteler
Date: 2013-10-23 14:16:40 -0700 (Wed, 23 Oct 2013)
New Revision: 58093
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_met.c
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c
Log:
i.landsat.toar: fixes for Landsat-8 metadata file support (author: E. Jorge Tizado)
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html 2013-10-23 16:17:37 UTC (rev 58092)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html 2013-10-23 21:16:40 UTC (rev 58093)
@@ -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/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c 2013-10-23 16:17:37 UTC (rev 58092)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c 2013-10-23 21:16:40 UTC (rev 58093)
@@ -22,6 +22,7 @@
*****************************************************************************/
double lsat_rad2ref(double rad, band_data * band)
{
+
return (double)(rad / band->K1);
}
@@ -87,8 +88,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 +110,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 +126,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/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c 2013-10-23 16:17:37 UTC (rev 58092)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c 2013-10-23 21:16:40 UTC (rev 58093)
@@ -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/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c 2013-10-23 16:17:37 UTC (rev 58092)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c 2013-10-23 21:16:40 UTC (rev 58093)
@@ -3,7 +3,7 @@
*
* MODULE: i.landsat.toar
*
- * AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
+ * AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
* Hamish Bowman (small grassification cleanups)
* Yann Chemin (v7 + L5TM _MTL.txt support) [removed after update]
*
@@ -32,7 +32,7 @@
struct History history;
struct GModule *module;
- struct Cell_head cellhd;
+ struct Cell_head cellhd, orig_cellhd;
char *mapset;
void *inrast, *outrast;
@@ -42,10 +42,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, *named;
+ struct Flag *frad, *print_meta, *named;
lsat_data lsat;
char band_in[GNAME_MAX], band_out[GNAME_MAX];
@@ -66,7 +66,7 @@
module->description =
_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OLI");
module->keywords =
- _("imagery, landsat, top-of-atmosphere reflectance, dos-type simple atmospheric correction");
+ _("imagery, Landsat, radiance, reflectance, brightness temperature, atmospheric correction");
/* It defines the different parameters */
input_prefix = G_define_option();
@@ -182,7 +182,7 @@
atmo->answer = "0.0";
atmo->guisection = _("Settings");
- /* It defines the different flags */
+ /* define the different flags */
frad = G_define_flag();
frad->key = 'r';
frad->description =
@@ -193,6 +193,10 @@
named->description =
_("Input raster maps use as extension the number of the band instead the code");
+ print_meta = G_define_flag();
+ print_meta->key = 'p';
+ print_meta->description = _("Print output metadata info");
+
/* options and afters parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -207,6 +211,7 @@
outputname = output_prefix->answer;
sensorname = sensor->answer ? sensor->answer : "";
+ G_get_window(&orig_cellhd);
G_zero(&lsat, sizeof(lsat));
if (adate->answer != NULL) {
@@ -238,8 +243,18 @@
if (met != NULL) {
lsat.flag = METADATAFILE;
lsat_metadata(met, &lsat);
- G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number,
- lsat.sensor);
+ if (print_meta->answer) {
+ G_message("Landsat-%d %s\n", lsat.number, lsat.sensor);
+ G_message("Number of bands = %d\n", lsat.bands);
+ G_message("Acquisition date = %s\n", lsat.date);
+ G_message("Acquisition time (h) = %f\n", lsat.time);
+ G_message("Sun elevation (degrees) = %f\n", lsat.sun_elev);
+ G_message("Sun azimuth (degrees) = %f\n", lsat.sun_az);
+ G_message("Product creation date = %s\n", lsat.creation);
+ exit(EXIT_SUCCESS);
+ }
+ G_debug(1, "lsat.number = %d, lsat.sensor = [%s]",
+ lsat.number, lsat.sensor);
if (!lsat.sensor || lsat.number > 8 || lsat.number < 1)
G_fatal_error(_("Failed to identify satellite"));
@@ -327,16 +342,17 @@
G_warning(_("Raster map <%s> not found"), band_in);
continue;
}
- 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_set_window(&cellhd);
+ if ((infd = G_open_cell_old(band_in, "")) < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), band_in);
in_data_type = G_raster_map_type(band_in, mapset);
if (in_data_type < 0)
- G_fatal_error(_("Unable to read data type of raster map <%s>"), band_in);
+ G_fatal_error(_("Unable to read data type of raster map <%s>"),
+ band_in);
inrast = G_allocate_raster_buf(in_data_type);
nrows = G_window_rows();
@@ -359,10 +375,9 @@
ptr = (void *)((DCELL *) inrast + col);
q = (int)*((DCELL *) ptr);
break;
- default:
- ptr = NULL;
- q = -1.;
- break;
+ default:
+ ptr = NULL;
+ q = -1.;
}
if (!G_is_null_value(ptr, in_data_type) &&
q >= lsat.band[i].qcalmin &&
@@ -411,13 +426,13 @@
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, " Earth-sun distance = %.8lf\n", lsat.dist_es);
+ fprintf(stderr, " Solar elevation angle = %.8lf\n", lsat.sun_elev);
fprintf(stderr, " Atmospheric correction: %s\n",
(method == UNCORRECTED ? "UNCORRECTED" : metho->answer));
if (method > DOS) {
fprintf(stderr,
- " percent of solar irradiance in path radiance = %.4lf\n",
+ " Percent of solar irradiance in path radiance = %.4lf\n",
percent);
}
for (i = 0; i < lsat.bands; i++) {
@@ -473,17 +488,18 @@
continue;
}
+ /* set same size as original band raster */
+ if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
+ G_fatal_error(_("Unable to read header of raster map <%s>"),
+ band_in);
+ G_set_window(&cellhd);
if ((infd = G_open_cell_old(band_in, mapset)) < 0)
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
in_data_type = G_raster_map_type(band_in, mapset);
if (in_data_type < 0)
- G_fatal_error(_("Unable to read data type of 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>"),
+ G_fatal_error(_("Unable to read data type of raster map <%s>"),
band_in);
- /* set same size as original band raster */
- G_set_window(&cellhd);
/* controlling, if we can write the raster */
if (G_legal_filename(band_out) < 0)
@@ -500,11 +516,9 @@
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);
@@ -526,7 +540,6 @@
default:
ptr = NULL;
qcal = -1.;
- break;
}
if (G_is_null_value(ptr, in_data_type) ||
qcal < lsat.band[i].qcalmin) {
@@ -584,56 +597,40 @@
/* Initialize the 'history' structure with basic info */
G_short_history(band_out, "raster", &history);
sprintf(history.edhist[0], " %s of Landsat-%d %s (method %s)",
- (frad->
- answer ? "Radiance" : (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],
- " Acquisition date (and time) ........... %s (%.4lf h)",
+ sprintf(history.edhist[2], " Acquisition date (and time) ........... %s (%.4lf h)",
lsat.date, lsat.time);
- sprintf(history.edhist[3],
- " Production date ....................... %s\n",
+ sprintf(history.edhist[3], " Production date ....................... %s\n",
lsat.creation);
- sprintf(history.edhist[4],
- " Earth-sun distance (d) ................ %.7lf",
+ sprintf(history.edhist[4], " Earth-sun distance (d) ................ %.7lf",
lsat.dist_es);
- sprintf(history.edhist[5],
- " Sun elevation (and azimuth) ........... %.5lf (%.5lf)",
+ sprintf(history.edhist[5], " Sun elevation (and azimuth) ........... %.5lf (%.5lf)",
lsat.sun_elev, lsat.sun_az);
- sprintf(history.edhist[6],
- " Digital number (DN) range ............. %.0lf to %.0lf",
+ sprintf(history.edhist[6], " Digital number (DN) range ............. %.0lf to %.0lf",
lsat.band[i].qcalmin, lsat.band[i].qcalmax);
- sprintf(history.edhist[7],
- " Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
+ sprintf(history.edhist[7], " Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
lsat.band[i].lmin, lsat.band[i].lmax);
- sprintf(history.edhist[8],
- " DN to Radiance (gain and bias) ........ %+.5lf and %+.5lf",
+ sprintf(history.edhist[8], " DN to Radiance (gain and bias) ........ %+.5lf and %+.5lf",
lsat.band[i].gain, lsat.band[i].bias);
if (lsat.band[i].thermal) {
- sprintf(history.edhist[9],
- " Temperature (K1 and K2) ............... %.3lf and %.3lf",
+ sprintf(history.edhist[9], " Temperature (K1 and K2) ............... %.3lf and %.3lf",
lsat.band[i].K1, lsat.band[i].K2);
history.edlinecnt = 10;
}
else {
sprintf(history.edhist[9],
- " Mean solar irradiance (ESUN) .......... %.3lf",
- lsat.band[i].esun);
+ " Mean solar irradiance (ESUN) .......... %.3lf", lsat.band[i].esun);
sprintf(history.edhist[10],
- " Radiance to Reflectance (divide by) ... %+.5lf",
- lsat.band[i].K1);
+ " Radiance to Reflectance (divide by) ... %+.5lf", lsat.band[i].K1);
history.edlinecnt = 11;
if (method > DOS) {
sprintf(history.edhist[11], " ");
- sprintf(history.edhist[12],
- " Dark object (%4d pixels) DN = ........ %d", pixel,
- dn_dark[i]);
- sprintf(history.edhist[13],
- " Mode in reflectance histogram ......... %.5lf",
- ref_mode);
+ sprintf(history.edhist[12], " Dark object (%4d pixels) DN = ........ %d", pixel, dn_dark[i]);
+ sprintf(history.edhist[13], " Mode in reflectance histogram ......... %.5lf", ref_mode);
history.edlinecnt = 14;
}
}
@@ -653,6 +650,7 @@
/* set raster timestamp from acq date? (see r.timestamp module) */
}
+ G_set_window(&orig_cellhd);
exit(EXIT_SUCCESS);
}
More information about the grass-commit
mailing list