[GRASS-SVN] r57858 - grass/branches/releasebranch_6_4/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Sep 28 10:13:59 PDT 2013
Author: neteler
Date: 2013-09-28 10:13:58 -0700 (Sat, 28 Sep 2013)
New Revision: 57858
Modified:
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html
grass/branches/releasebranch_6_4/imagery/i.landsat.toar/earth_sun.c
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:
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-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/description.html 2013-09-28 17:13:58 UTC (rev 57858)
@@ -14,7 +14,7 @@
<p>
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
+if the solar elevation is given the value of the metadata file is
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.
@@ -60,26 +60,6 @@
units, <em>e</em> is the solar elevation angle, and <em>Esun</em> is
the mean solar exoatmospheric irradiance in W/(m² * µm).
-<h2>Corrected at-sensor values (method=corrected)</h2>
-
-At-sensor reflectance values range from zero to one, whereas at-sensor
-radiance must be greater or equal to zero. However, since Lmin can be
-a negative number then the at-sensor values can also be negative. To
-avoid these possible negative values and set the minimum possible
-values at-sensor to zero, this method corrects the radiance to output
-a corrected at-sensor values using the equations (not for thermal
-bands):
-<ul>
- <li> radiance = (uncorrected_radiance - Lmin) </li>
- <li> reflectance = radiance / sun_radiance </li>
-</ul>
-
-<p>
-<b>Note</b>: Other possibility to avoid negative values is set to zero
-this values (radiance and/or reflectance), but this option is ease
-with uncorrected method
-and <em><a href="r.mapcalc.html">r.mapcalc</a></em>.
-
<h2>Simplified at-surface values (method=dos[1-4])</h2>
Atmospheric correction and reflectance calibration remove the path
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/earth_sun.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/earth_sun.c 2013-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/earth_sun.c 2013-09-28 17:13:58 UTC (rev 57858)
@@ -1096,6 +1096,8 @@
R4 = ln_calc_series(earth_radius_r4, RADIUS_R4, t);
R5 = ln_calc_series(earth_radius_r5, RADIUS_R5, t);
- return (R0 + R1 * t + R2 * t * t + R3 * t * t * t + R4 * t * t * t * t +
- R5 * t * t * t * t * t);
+ return (R0 +
+ R1 * t +
+ R2 * t * t +
+ R3 * t * t * t + R4 * t * t * t * t + R5 * t * t * t * t * t);
}
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c 2013-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.c 2013-09-28 17:13:58 UTC (rev 57858)
@@ -30,7 +30,7 @@
*****************************************************************************/
double lsat_rad2temp(double rad, band_data * band)
{
- return (double)(band->K2 / log((band->K1 / rad) + 1.0));
+ return (double)(band->K2 / log((band->K1 / rad) + 1.));
}
/****************************************************************************
@@ -43,15 +43,15 @@
* i : band number
* method : level of atmospheric correction
* percent : percent of solar irradiance in path radiance
- * dos : digital number of dark object for DOS
+ * dark : digital number of object for DOS
*****************************************************************************/
#define abs(x) (((x)>0)?(x):(-x))
void lsat_bandctes(lsat_data * lsat, int i, char method,
- double percent, int dos, double rayleigh)
+ double percent, int dark, double rayleigh)
{
- double pi_d2, sin_e, cos_v, rad_sun;
+ double pi_d2, sin_e, cos_v;
/* TAUv = at. transmittance surface-sensor */
/* TAUz = at. transmittance sun-surface */
@@ -63,9 +63,9 @@
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
- */
+ * Radiance to reflectance coefficient, only NO thermal bands.
+ * K1 and K2 variables are also utilized as thermal constants
+ */
if (lsat->band[i].thermal == 0) {
switch (method) {
case DOS2:
@@ -97,10 +97,11 @@
case DOS4:
{
double Ro =
- (lsat->band[i].lmax - lsat->band[i].lmin) * (dos -
- lsat->band
- [i].qcalmin)
- / (lsat->band[i].qcalmax - lsat->band[i].qcalmin) +
+ (lsat->band[i].lmax - lsat->band[i].lmin) * (dark -
+ lsat->
+ band[i].
+ qcalmin) /
+ (lsat->band[i].qcalmax - lsat->band[i].qcalmin) +
lsat->band[i].lmin;
double Tv = 1.;
double Tz = 1.;
@@ -112,10 +113,10 @@
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);
+ 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( abs(TAUv - Tv) > 0.0000001 || abs(TAUz - Tz) > 0.0000001); */
} while (TAUv != Tv && TAUz != Tz);
TAUz = (Tz < 1. ? Tz : 1.);
TAUv = (Tv < 1. ? Tv : 1.);
@@ -128,41 +129,43 @@
Edown = 0.;
break;
}
- rad_sun = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
+ lsat->band[i].K2 = 0.;
+ lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2; /* rad_sun */
if (method > DOS)
G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
TAUv, TAUz, Edown);
-
- lsat->band[i].K1 = rad_sun;
- lsat->band[i].K2 = 0.;
}
/** 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));
+ * 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);
if (method == UNCORRECTED || lsat->band[i].thermal) {
/* L = G * (DN - Qmin) + Lmin
- -> bias = Lmin - G * Qmin */
+ * -> bias = Lmin - G * Qmin
+ */
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 * dark + B - p * rad_sun) =
- G * DNsat - G * dark + p * rad_sun
- -> bias = p * rad_sun - G * dark */
- lsat->band[i].bias = percent * rad_sun - lsat->band[i].gain * dos;
- }
+ /*
+ 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;
}
}
Modified: grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h 2013-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat.h 2013-09-28 17:13:58 UTC (rev 57858)
@@ -25,38 +25,38 @@
typedef struct
{
- int number; /* Band number */
- int code; /* Band code */
+ int number; /* Band number */
+ int code; /* Band code */
- double wavemax, wavemin; /* Wavelength in µm */
+ 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 */
+ double lmax, lmin; /* Spectral radiance */
+ double qcalmax, qcalmin; /* Quantized calibrated pixel */
- char thermal; /* Flag to thermal band */
- double gain, bias; /* Gain and Bias of sensor */
- double K1, K2; /* Thermal calibration or
- Rad2Ref constants */
-
+ char thermal; /* Flag to thermal band */
+ double gain, bias; /* Gain and Bias of sensor */
+ double K1, K2; /* Thermal calibration
+ * or Rad2Ref constants */
} band_data;
typedef struct
{
- int flag;
- unsigned char number; /* Landsat number */
+ int flag; /* Line-data or file-data */
+ unsigned char number; /* Landsat number */
- char creation[11]; /* Image production date */
- char date[11]; /* Image acquisition date */
+ char creation[11]; /* Image production date */
+ char date[11]; /* Image acquisition date */
+ double time; /* Image acquisition time */
- double dist_es; /* Distance Earth-Sun */
- double sun_elev; /* Sun elevation */
- double sun_az; /* Sun Azimuth */
- double time; /* Image Acquisition Time */
+ double dist_es; /* Distance Earth-Sun */
+ double sun_elev; /* Sun elevation */
+ double sun_az; /* Sun azimuth */
- char sensor[10]; /* Type of sensor: MSS, TM, ETM+, OLI/TIRS */
- int bands; /* Total number of bands */
- band_data band[MAX_BANDS]; /* Data for each band */
+ char sensor[10]; /* Sensor: MSS, TM, ETM+, OLI/TIRS */
+ int bands; /* Total number of bands */
+ band_data band[MAX_BANDS]; /* Data for each band */
+
} lsat_data;
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-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_met.c 2013-09-28 17:13:58 UTC (rev 57858)
@@ -12,7 +12,7 @@
#define PI M_PI
#define D2R M_PI/180.
-#define MAX_STR 127
+#define MAX_STR 256
#define METADATA_SIZE 65535 /* MTL.txt file size 65535 bytes */
#define TM5_MET_SIZE 28700 /* .met file size 28686 bytes */
@@ -26,6 +26,15 @@
dest[i] = '\0';
}
+inline void date_replace_slash(char *str)
+{
+ while (*str) {
+ if (*str == '/')
+ *str = '-';
+ *str++;
+ }
+}
+
void get_metformat(const char metadata[], char *key, char value[])
{
int i = 0;
@@ -77,12 +86,13 @@
void (*get_mtldata) (const char[], char *, char[]);
int i, j, ver_mtl;
+ /* store metadata in ram */
if ((f = fopen(metafile, "r")) == NULL)
G_fatal_error(_("Metadata file <%s> not found"), metafile);
-
i = fread(mtldata, METADATA_SIZE, 1, f);
+ (void)fclose(f);
- /* get version of the metadata file */
+ /* set version of the metadata file */
get_mtldata =
(strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
ver_mtl = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
@@ -92,7 +102,6 @@
if (value[0] == '\0') {
get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
}
-
i = 0;
while (value[i] && (value[i] < '0' || value[i] > '9'))
i++;
@@ -111,8 +120,10 @@
get_mtldata(mtldata, "CALENDARDATE", value);
}
}
- if (value[0] != '\0')
+ if (value[0] != '\0') {
+ date_replace_slash(value);
chrncpy(lsat->date, value, 10);
+ }
else
G_warning("Using adquisition date from the command line 'date'");
@@ -123,8 +134,10 @@
get_mtldata(mtldata, "PRODUCTIONDATETIME", value);
}
}
- if (value[0] != '\0')
+ if (value[0] != '\0') {
+ date_replace_slash(value);
chrncpy(lsat->creation, value, 10);
+ }
else
G_warning
("Using production date from the command line 'product_date'");
@@ -132,7 +145,7 @@
get_mtldata(mtldata, "SUN_AZIMUTH", value);
lsat->sun_az = atof(value);
if (lsat->sun_az == 0.)
- G_warning("Sun azimuth is %f", lsat->sun_az);
+ G_warning("Sun azimuth is %f", lsat->sun_az);
get_mtldata(mtldata, "SUN_ELEVATION", value);
if (value[0] == '\0') {
@@ -144,14 +157,14 @@
get_mtldata(mtldata, "SCENE_CENTER_TIME", value);
if (value[0] == '\0') {
- get_mtldata(mtldata, "SCENE_CENTER_SCAN_TIME", value);
+ get_mtldata(mtldata, "SCENE_CENTER_SCAN_TIME", value);
}
- /* Remove trailing 'z'*/
- value[strlen(value) - 1]='\0';
- /* Cast from hh:mm:ss into hh.hhh*/
- G_llres_scan(value, &lsat->time);
+ if (value[0] != '\0') {
+ value[strlen(value) - 1] = '\0'; /* Remove trailing 'z' */
+ G_llres_scan(value, &lsat->time); /* Cast from hh:mm:ss into hh.hhh */
+ }
if (lsat->time == 0.)
- G_warning("Time is %f", lsat->time);
+ G_warning("Scene time is %f", lsat->time);
/* Fill the data with the basic/default sensor parameters */
switch (lsat->number) {
@@ -195,9 +208,8 @@
break;
}
case 8:
- set_LDCM(lsat);
+ set_OLI(lsat);
break;
-
default:
G_warning("Unable to recognize satellite platform [%d]",
lsat->number);
@@ -206,7 +218,7 @@
/* Update the information from metadata file, if infile */
if (get_mtldata == get_mtlformat) {
- if (ver_mtl == 0) { /* now MTLold.txt */
+ if (ver_mtl == 0) {
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);
@@ -223,14 +235,38 @@
lsat->band[i].qcalmin = atof(value);
}
}
- else { /* now MTL.txt */
-
+ 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 */
+ 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");
for (i = 0; i < lsat->bands; i++) {
- if (lsat->band[i].thermal && lsat->number == 7) {
+ if (lsat->number == 7 && lsat->band[i].thermal) {
sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
lsat->band[i].code - 60);
get_mtldata(mtldata, key, value);
@@ -266,29 +302,20 @@
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);
+ 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 {
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);
@@ -296,17 +323,14 @@
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.;
+ /* 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;
- /* ----- */
+ /* ... 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);
@@ -317,34 +341,32 @@
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_refmu[mode], lsat->band[i].code);
+ get_mtldata(mtldata, key, value);
+ lsat->band[i].K1 = atof(value);
sprintf(key, fmt_refad[mode], 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) / (sin(D2R *
- lsat->
- sun_elev)
- *
- lsat->band[i].
- K2);
+ 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;
+ lsat->band[i].esun = (esun1 + esun2) / 2.;
+ */
}
}
}
- }
- /* Other possible values in file */
- get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value);
- if (value[0] != '\0') {
- lsat->dist_es = atof(value);
- G_verbose_message
+ G_warning
("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
}
}
}
else {
- G_verbose_message("Metada file is MET file");
+ G_verbose_message("Metadata file is MET file");
G_verbose_message
("RADIANCE & QUANTIZE from band setting of the metadata file");
for (i = 0; i < lsat->bands; i++) {
@@ -374,6 +396,5 @@
}
}
- (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-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/landsat_set.c 2013-09-28 17:13:58 UTC (rev 57858)
@@ -68,8 +68,10 @@
/* blue, green, red, near infrared, shortwave IR, thermal IR, shortwave IR, panchromatic */
int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
- double wmin[] = { 0.450, 0.525, 0.630, 0.75, 1.55, 10.40, 10.40, 2.09, 0.52 };
- double wmax[] = { 0.515, 0.605, 0.690, 0.90, 1.75, 12.50, 12.50, 2.35, 0.90 };
+ double wmin[] =
+ { 0.450, 0.525, 0.630, 0.75, 1.55, 10.40, 10.40, 2.09, 0.52 };
+ double wmax[] =
+ { 0.515, 0.605, 0.690, 0.90, 1.75, 12.50, 12.50, 2.35, 0.90 };
/* 30, 30, 30, 30, 30, 60 (after Feb. 25, 2010: 30), 30, 15 */
strcpy(lsat->sensor, "ETM+");
@@ -87,7 +89,7 @@
return;
}
-void sensor_LDCM(lsat_data * lsat)
+void sensor_OLI(lsat_data * lsat)
{
int i;
@@ -95,8 +97,12 @@
* 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");
@@ -107,7 +113,7 @@
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].qcalmax = 65535.;
lsat->band[i].qcalmin = 1.;
lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
}
@@ -444,7 +450,8 @@
lmax = Lmax[i];
lmin = Lmin[i];
- if (i == 2) { /* in Chander, Markham and Barsi 2007 */
+ /* in Chander, Markham and Barsi 2007 */
+ if (i == 2) {
julian = julian_char(lsat->date); /* Yes, here acquisition date */
if (julian >= julian_char("1992-01-01")) {
lmax[0] = 193.0;
@@ -516,10 +523,7 @@
/* Thermal band calibration constants: K1 = 666.09 K2 = 1282.71 */
julian = julian_char(lsat->creation);
- if (julian < julian_char("2000-07-01"))
- k = 0;
- else
- k = 1;
+ k = ((julian < julian_char("2000-07-01")) ? 0 : 1);
lsat->number = 7;
sensor_ETM(lsat);
@@ -552,29 +556,32 @@
* PURPOSE: Store values of Landsat-8 OLI/TIRS
* February 14, 2013
*****************************************************************************/
-void set_LDCM(lsat_data * lsat)
+void set_OLI(lsat_data * lsat)
{
int i, j;
double *lmax, *lmin;
/* Spectral radiances at detector */
-
- /* uncorrected values */
+ /* estimates */
double Lmax[][11] = {
- {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
+ {755.8, 770.7, 705.7, 597.7, 362.7, 91.4, 29.7, 673.3, 149.0, 22.0,
+ 22.0}
};
double Lmin[][11] = {
- {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
+ {-62.4, -63.6, -58.3, -49.4, -30.0, -7.5, -2.5, -55.6, -12.3, 0.1, 0.1}
};
+
/* Solar exoatmospheric spectral irradiances */
/* estimates */
- double esun[] = { 2062., 21931., 1990., 1688., 1037., 268.6, 94.6, 1892., 399.0, 0., 0. };
+ double esun[] =
+ { 2026.8, 2066.8, 1892.5, 1602.8, 972.6, 245.0, 79.7, 1805.5, 399.7,
+0., 0. };
lmax = Lmax[0];
lmin = Lmin[0];
lsat->number = 8;
- sensor_LDCM(lsat);
+ sensor_OLI(lsat);
lsat->dist_es = earth_sun(lsat->date);
@@ -584,8 +591,9 @@
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;
+ lsat->band[i].K1 = (lsat->band[i].number == 10 ? 774.89 : 480.89);
+ lsat->band[i].K2 =
+ (lsat->band[i].number == 10 ? 1321.08 : 1201.14);
}
}
G_debug(1, "Landsat-8 OLI/TIRS");
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-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/local_proto.h 2013-09-28 17:13:58 UTC (rev 57858)
@@ -17,6 +17,6 @@
void set_ETM(lsat_data *, char[]);
-void set_LDCM(lsat_data *);
+void set_OLI(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-09-27 14:28:32 UTC (rev 57857)
+++ grass/branches/releasebranch_6_4/imagery/i.landsat.toar/main.c 2013-09-28 17:13:58 UTC (rev 57858)
@@ -25,6 +25,8 @@
#include "local_proto.h"
+#define QCALMAX 65536 /* L1-7=256 but L8=65536 */
+
int main(int argc, char *argv[])
{
struct History history;
@@ -40,28 +42,29 @@
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;
lsat_data lsat;
char band_in[GNAME_MAX], band_out[GNAME_MAX];
- int i, j, q, method, pixel, dn_dark[MAX_BANDS], dn_mode[MAX_BANDS];
+ int i, j, q, method, pixel, dn_dark[MAX_BANDS], dn_mode[MAX_BANDS],
+ dn_sat;
double qcal, rad, ref, percent, ref_mode, rayleigh;
+ unsigned long hist[QCALMAX], h_max;
struct Colors colors;
struct FPRange range;
double min, max;
- unsigned long hist[256], h_max;
- /* initialize GIS environment */
+ /* It initializes GIS environment */
G_gisinit(argv[0]);
- /* initialize module */
+ /* It initializes module */
module = G_define_module();
module->description =
- _("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+.");
+ _("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");
@@ -93,7 +96,7 @@
sensor->label = _("Spacecraft sensor");
sensor->description =
_("Required only if 'metfile' not given (recommended for sanity)");
- sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,ot8";
+ sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,oli8";
sensor->descriptions =
_("mss1;Landsat_1 MSS;"
"mss2;Landsat_2 MSS;"
@@ -102,8 +105,8 @@
"mss5;Landsat_5 MSS;"
"tm4;Landsat_4 TM;"
"tm5;Landsat_5 TM;"
- "tm7;Landsat_7 ETM+;"
- "ot8;Landsat_8 OLI/TIRS");
+ "tm7;Landsat_7 ETM+;"
+ "oli8;Landsat_8 OLI/TIRS");
sensor->required = NO; /* perhaps YES for clarity */
sensor->guisection = _("Metadata");
@@ -111,7 +114,7 @@
metho->key = "method";
metho->type = TYPE_STRING;
metho->required = NO;
- metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
+ metho->options = "uncorrected,dos1,dos2,dos2b,dos3,dos4";
metho->label = _("Atmospheric correction method");
metho->description = _("Atmospheric correction method");
metho->answer = "uncorrected";
@@ -179,7 +182,7 @@
atmo->answer = "0.0";
atmo->guisection = _("Settings");
- /* define the different flags */
+ /* It defines the different flags */
frad = G_define_flag();
frad->key = 'r';
frad->description =
@@ -195,10 +198,10 @@
exit(EXIT_FAILURE);
- /*****************************************
- * ---------- START --------------------
- * Stores options and flag to variables
- *****************************************/
+ /*****************************************
+ * ---------- START --------------------
+ * Stores options and flags to variables
+ *****************************************/
met = metfn->answer;
inputname = input_prefix->answer;
outputname = output_prefix->answer;
@@ -228,9 +231,10 @@
pixel = atoi(dark->answer);
rayleigh = atof(atmo->answer);
- /* Data from metadata file */
- /* Unnecessary because G_zero filled, but for sanity */
- lsat.flag = NOMETADATAFILE;
+ /*
+ * Data from metadata file
+ */
+ lsat.flag = NOMETADATAFILE; /* Unnecessary because G_zero filled, but for sanity */
if (met != NULL) {
lsat.flag = METADATAFILE;
lsat_metadata(met, &lsat);
@@ -240,28 +244,29 @@
if (!lsat.sensor || lsat.number > 8 || lsat.number < 1)
G_fatal_error(_("Failed to identify satellite"));
- G_debug(1, "Landsat-%d %s with data set in met file [%s]",
+ G_debug(1, "Landsat-%d %s with data set in metadata file [%s]",
lsat.number, lsat.sensor, met);
- /* Overwrite solar elevation of metadata file */
- if (elev->answer != NULL)
+ if (elev->answer != NULL) {
lsat.sun_elev = atof(elev->answer);
+ G_warning("Overwriting solar elevation of metadata file");
+ }
}
- /* Data from date and solar elevation */
+ /*
+ * Data from command line
+ */
else if (adate->answer == NULL || elev->answer == NULL) {
G_fatal_error(_("Lacking '%s' and/or '%s' for this satellite"),
adate->key, elev->key);
}
else {
- /* 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"));
+ G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) characters"));
set_ETM(&lsat, bgain->answer);
}
- /* Not need gain */
- else if (strcmp(sensorname, "ot8") == 0)
- set_LDCM(&lsat);
+ else if (strcmp(sensorname, "oli8") == 0)
+ set_OLI(&lsat);
else if (strcmp(sensorname, "tm5") == 0)
set_TM5(&lsat);
else if (strcmp(sensorname, "tm4") == 0)
@@ -278,13 +283,13 @@
set_MSS1(&lsat);
else
G_fatal_error(_("Unknown satellite type (defined by '%s')"),
- sensor->key);
+ sensorname);
}
/*****************************************
* ------------ PREPARATION --------------
*****************************************/
- if (strcasecmp(metho->answer, "corrected") == 0)
+ if (strcasecmp(metho->answer, "corrected") == 0) /* deleted 2013 */
method = CORRECTED;
else if (strcasecmp(metho->answer, "dos1") == 0)
method = DOS1;
@@ -300,21 +305,21 @@
method = UNCORRECTED;
/*
- if (metho->answer[3] == 'r') method = CORRECTED;
- else if (metho->answer[3] == '1') method = DOS1;
- else if (metho->answer[3] == '2') method = (metho->answer[4] == '\0') ? DOS2 : DOS2b;
- else if (metho->answer[3] == '3') method = DOS3;
- else if (metho->answer[3] == '4') method = DOS4;
- else method = UNCORRECTED;
+ if (metho->answer[3] == '2') method = (metho->answer[4] == '\0') ? DOS2 : DOS2b;
+ else if (metho->answer[3] == '1') method = DOS1;
+ else if (metho->answer[3] == '3') method = DOS3;
+ else if (metho->answer[3] == '4') method = DOS4;
+ else method = UNCORRECTED;
*/
for (i = 0; i < lsat.bands; i++) {
dn_mode[i] = 0;
dn_dark[i] = (int)lsat.band[i].qcalmin;
- /* Calculate dark pixel */
+ dn_sat = (int)(0.90 * lsat.band[i].qcalmax);
+ /* Begin: calculate dark pixel */
if (method > DOS && !lsat.band[i].thermal) {
- for (j = 0; j < 256; j++)
- hist[j] = 0L;
+ for (q = 0; q <= lsat.band[i].qcalmax; q++)
+ hist[q] = 0L;
sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
mapset = G_find_cell2(band_in, "");
@@ -354,54 +359,56 @@
break;
}
if (!G_is_null_value(ptr, in_data_type) &&
- q >= lsat.band[i].qcalmin && q < 256)
+ q >= lsat.band[i].qcalmin &&
+ q <= lsat.band[i].qcalmax)
hist[q]++;
}
}
/* DN of dark object */
- for (j = lsat.band[i].qcalmin; j < 256; j++) {
+ for (j = lsat.band[i].qcalmin; j <= lsat.band[i].qcalmax; j++) {
if (hist[j] >= (unsigned int)pixel) {
dn_dark[i] = j;
break;
}
}
- /* Mode of DN */
+ /* Mode of DN (exclude potentially saturated) */
h_max = 0L;
- for (j = lsat.band[i].qcalmin; j < 241; j++) { /* Exclude potentially saturated < 240 */
+ for (j = lsat.band[i].qcalmin; j < dn_sat; j++) {
/* G_debug(5, "%d-%ld", j, hist[j]); */
if (hist[j] > h_max) {
h_max = hist[j];
dn_mode[i] = j;
}
}
- G_verbose_message("... DN = %.2d [%lu] : mode %.2d [%lu] %s",
- dn_dark[i], hist[dn_dark[i]],
- dn_mode[i], hist[dn_mode[i]],
- hist[255] >
- hist[dn_mode[i]] ? ", excluding DN > 241" : "");
+ G_verbose_message
+ ("... DN = %.2d [%lu] : mode %.2d [%lu], exclusing DN > %d",
+ dn_dark[i], hist[dn_dark[i]], dn_mode[i], hist[dn_mode[i]],
+ dn_sat);
G_free(inrast);
G_close_cell(infd);
}
+ /* End: calculate dark pixel */
+
/* Calculate transformation constants */
lsat_bandctes(&lsat, i, method, percent, dn_dark[i], rayleigh);
}
- /* 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);
+ /*
+ * 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,
+ fprintf(stderr, "\n 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, " Atmospheric correction: %s\n",
- (method == CORRECTED ? "CORRECTED"
- : (method == UNCORRECTED ? "UNCORRECTED" : metho->answer)));
+ (method == UNCORRECTED ? "UNCORRECTED" : metho->answer));
if (method > DOS) {
fprintf(stderr,
" percent of solar irradiance in path radiance = %.4lf\n",
@@ -409,16 +416,15 @@
}
for (i = 0; i < lsat.bands; i++) {
fprintf(stderr, "-------------------\n");
- fprintf(stderr, " BAND %d %s(code %d)\n",
- lsat.band[i].number,
+ fprintf(stderr, " BAND %d %s(code %d)\n", lsat.band[i].number,
(lsat.band[i].thermal ? "thermal " : ""),
lsat.band[i].code);
fprintf(stderr,
" calibrated digital number (DN): %.1lf to %.1lf\n",
lsat.band[i].qcalmin, lsat.band[i].qcalmax);
- fprintf(stderr, " calibration constants (L): %.3lf to %.3lf\n",
+ fprintf(stderr, " calibration constants (L): %.5lf to %.5lf\n",
lsat.band[i].lmin, lsat.band[i].lmax);
- fprintf(stderr, " at-%s radiance = %.8lf * DN + %.3lf\n",
+ fprintf(stderr, " at-%s radiance = %.8lf * DN + %.5lf\n",
(method > DOS ? "surface" : "sensor"), lsat.band[i].gain,
lsat.band[i].bias);
if (lsat.band[i].thermal) {
@@ -428,7 +434,7 @@
}
else {
fprintf(stderr,
- " mean solar exoatmospheric irradiance (ESUN): %.3lf\n",
+ " mean solar exoatmospheric irradiance (ESUN): %.5lf\n",
lsat.band[i].esun);
fprintf(stderr, " at-%s reflectance = radiance / %.5lf\n",
(method > DOS ? "surface" : "sensor"),
@@ -478,18 +484,19 @@
if ((outfd = G_open_raster_new(band_out, DCELL_TYPE)) < 0)
G_fatal_error(_("Unable to create raster map <%s>"), band_out);
- /* Allocate input and output buffer */
+ /* allocate input and output buffer */
inrast = G_allocate_raster_buf(in_data_type);
outrast = G_allocate_raster_buf(DCELL_TYPE);
nrows = G_window_rows();
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);
@@ -541,17 +548,20 @@
ref_mode = lsat_rad2ref(ref_mode, &lsat.band[i]);
}
- /* ================================================================= */
G_free(inrast);
G_close_cell(infd);
G_free(outrast);
G_close_cell(outfd);
- /* needed ?
- if (out_type != CELL_TYPE)
- G_quantize_fp_map_range(band_out, G_mapset(), 0., 360., 0, 360);
+
+ /*
+ * needed?
+ * if (out_type != CELL_TYPE)
+ * G_quantize_fp_map_range(band_out, G_mapset(), 0., 360., 0,
+ * 360);
*/
+
/* set grey255 colortable */
G_init_colors(&colors);
G_read_fp_range(band_out, G_mapset(), &range);
@@ -559,7 +569,7 @@
G_make_grey_scale_fp_colors(&colors, min, max);
G_write_colors(band_out, G_mapset(), &colors);
- /* Initialize the 'hist' structure with basic info */
+ /* 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->
@@ -568,51 +578,55 @@
"Reflectance")), lsat.number,
lsat.sensor, metho->answer);
sprintf(history.edhist[1],
- "----------------------------------------------------------------");
+ "-----------------------------------------------------------------");
sprintf(history.edhist[2],
- " Acquisition date ...................... %s", lsat.date);
+ " Acquisition date (and time) ........... %s (%.4lf h)",
+ lsat.date, lsat.time);
sprintf(history.edhist[3],
" Production date ....................... %s\n",
lsat.creation);
sprintf(history.edhist[4],
- " Earth-sun distance (d) ................ %.8lf",
+ " Earth-sun distance (d) ................ %.7lf",
lsat.dist_es);
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",
lsat.band[i].qcalmin, lsat.band[i].qcalmax);
- sprintf(history.edhist[6],
- " Calibration constants (Lmin to Lmax) .. %+.3lf to %+.3lf",
+ sprintf(history.edhist[7],
+ " Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
lsat.band[i].lmin, lsat.band[i].lmax);
- sprintf(history.edhist[7],
+ 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[8],
+ sprintf(history.edhist[9],
" Temperature (K1 and K2) ............... %.3lf and %.3lf",
lsat.band[i].K1, lsat.band[i].K2);
- history.edlinecnt = 9;
+ history.edlinecnt = 10;
}
else {
- sprintf(history.edhist[8],
+ sprintf(history.edhist[9],
" Mean solar irradiance (ESUN) .......... %.3lf",
lsat.band[i].esun);
- sprintf(history.edhist[9],
- " Reflectance = Radiance divided by ..... %.5lf",
+ sprintf(history.edhist[10],
+ " Radiance to Reflectance (divide by) ... %+.5lf",
lsat.band[i].K1);
- history.edlinecnt = 10;
+ history.edlinecnt = 11;
if (method > DOS) {
- sprintf(history.edhist[10], " ");
- sprintf(history.edhist[11],
+ sprintf(history.edhist[11], " ");
+ sprintf(history.edhist[12],
" Dark object (%4d pixels) DN = ........ %d", pixel,
dn_dark[i]);
- sprintf(history.edhist[12],
+ sprintf(history.edhist[13],
" Mode in reflectance histogram ......... %.5lf",
ref_mode);
- history.edlinecnt = 13;
+ history.edlinecnt = 14;
}
}
sprintf(history.edhist[history.edlinecnt],
- "-----------------------------------------------------------------");
+ "------------------------------------------------------------------");
history.edlinecnt++;
G_command_history(&history);
@@ -625,7 +639,7 @@
else
G_write_raster_units(band_out, "unitless");
- /* set raster timestamp from acq date? (see r.timestamp module) */
+ /* set raster timestamp from acq date? (see r.timestamp module) */
}
exit(EXIT_SUCCESS);
More information about the grass-commit
mailing list