[GRASS-SVN] r57936 - grass/trunk/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Oct 4 04:16:11 PDT 2013
Author: neteler
Date: 2013-10-04 04:16:10 -0700 (Fri, 04 Oct 2013)
New Revision: 57936
Modified:
grass/trunk/imagery/i.landsat.toar/earth_sun.c
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.h
grass/trunk/imagery/i.landsat.toar/landsat_met.c
grass/trunk/imagery/i.landsat.toar/landsat_set.c
grass/trunk/imagery/i.landsat.toar/local_proto.h
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 G6 version
Modified: grass/trunk/imagery/i.landsat.toar/earth_sun.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/earth_sun.c 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/earth_sun.c 2013-10-04 11:16:10 UTC (rev 57936)
@@ -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/trunk/imagery/i.landsat.toar/i.landsat.toar.html
===================================================================
--- grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html 2013-10-04 11:16:10 UTC (rev 57936)
@@ -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
@@ -204,10 +184,9 @@
<h2>SEE ALSO</h2>
<em>
-<a href="i.aster.toar.html">i.aster.toar</a>,
-<a href="i.atcorr.html">i.atcorr</a>,
-<a href="r.mapcalc.html">r.mapcalc</a>,
-<a href="r.in.gdal.html">r.in.gdal</a>
+ <a href="i.atcorr.html">i.atcorr</a>,
+ <a href="r.mapcalc.html">r.mapcalc</a>,
+ <a href="r.in.gdal.html">r.in.gdal</a>
</em>
<h2>AUTHOR</h2>
Modified: grass/trunk/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat.c 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/landsat.c 2013-10-04 11:16:10 UTC (rev 57936)
@@ -5,9 +5,9 @@
#include "landsat.h"
-#define PI M_PI
-#define R2D 180. / M_PI
-#define D2R M_PI / 180.
+#define PI M_PI
+#define R2D 180./M_PI
+#define D2R M_PI/180.
/****************************************************************************
* PURPOSE: Calibrated Digital Number to at-satellite Radiance
@@ -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/trunk/imagery/i.landsat.toar/landsat.h
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat.h 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/landsat.h 2013-10-04 11:16:10 UTC (rev 57936)
@@ -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/trunk/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat_met.c 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/landsat_met.c 2013-10-04 11:16:10 UTC (rev 57936)
@@ -9,10 +9,10 @@
#include "local_proto.h"
#include "earth_sun.h"
-#define PI M_PI
-#define D2R M_PI / 180.
+#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;
@@ -33,14 +42,16 @@
char *ptr;
if (ptrmet != NULL) {
- ptr = strstr(ptrmet, " VALUE ");
- if (ptr != NULL) {
- while (*ptr++ != '=') ;
- while (*ptr <= ' ' || *ptr == '\"')
+ ptr = strstr(ptrmet, " VALUE ");
+ if (ptr != NULL) {
+ while (*ptr++ != '=') ;
+ while (*ptr <= ' ')
ptr++;
- while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
+ if (*ptr == '\"')
+ ptr++;
+ while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
value[i++] = *ptr++;
- }
+ }
}
value[i] = '\0';
}
@@ -51,12 +62,12 @@
char *ptr = strstr(metadata, key);
if (ptr != NULL) {
- while (*ptr++ != '=') ;
-
- while (*ptr <= ' ' || *ptr == '\"')
+ while (*ptr++ != '=') ;
+ while (*ptr <= ' ')
ptr++;
-
- while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
+ if (*ptr == '\"')
+ ptr++;
+ while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
value[i++] = *ptr++;
}
value[i] = '\0';
@@ -70,83 +81,90 @@
void lsat_metadata(char *metafile, lsat_data * lsat)
{
FILE *f;
- char metadata[METADATA_SIZE];
+ char mtldata[METADATA_SIZE];
char key[MAX_STR], value[MAX_STR];
- void (*get_metadata) (const char [], char *, char []);
- int i, j, ver_meta;
+ 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);
- i = fread(metadata, METADATA_SIZE, 1, f);
+ /* 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;
- /* get version of the metadata file */
- get_metadata =
- (strstr(metadata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
- ver_meta = (strstr(metadata, "QCALMAX_BAND") != NULL) ? 0 : 1;
-
/* Fill with product metadata */
- get_metadata(metadata, "SPACECRAFT_ID", value);
+ get_mtldata(mtldata, "SPACECRAFT_ID", value);
if (value[0] == '\0') {
- get_metadata(metadata, "PLATFORMSHORTNAME", value);
+ get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
}
-
i = 0;
while (value[i] && (value[i] < '0' || value[i] > '9'))
i++;
lsat->number = atoi(value + i);
- get_metadata(metadata, "SENSOR_ID", value);
+ get_mtldata(mtldata, "SENSOR_ID", value);
if (value[0] == '\0') {
- get_metadata(metadata, "SENSORSHORTNAME", value);
+ get_mtldata(mtldata, "SENSORSHORTNAME", value);
}
chrncpy(lsat->sensor, value, 8);
- get_metadata(metadata, "DATE_ACQUIRED", value);
+ get_mtldata(mtldata, "DATE_ACQUIRED", value);
if (value[0] == '\0') {
- get_metadata(metadata, "ACQUISITION_DATE", value);
- if (value[0] == '\0') {
- get_metadata(metadata, "CALENDARDATE", value);
- }
+ get_mtldata(mtldata, "ACQUISITION_DATE", value);
+ if (value[0] == '\0') {
+ get_mtldata(mtldata, "CALENDARDATE", value);
+ }
}
- chrncpy(lsat->date, value, 10);
+ if (value[0] != '\0') {
+ date_replace_slash(value);
+ chrncpy(lsat->date, value, 10);
+ }
+ else
+ G_warning("Using adquisition date from the command line 'date'");
- get_metadata(metadata, "FILE_DATE", value);
+ get_mtldata(mtldata, "FILE_DATE", value);
if (value[0] == '\0') {
- get_metadata(metadata, "CREATION_TIME", value);
- if (value[0] == '\0') {
- get_metadata(metadata, "PRODUCTIONDATETIME", value);
- }
+ get_mtldata(mtldata, "CREATION_TIME", value);
+ if (value[0] == '\0') {
+ 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'");
- get_metadata(metadata, "SUN_AZIMUTH", value);
+ 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_metadata(metadata, "SUN_ELEVATION", value);
+ get_mtldata(mtldata, "SUN_ELEVATION", value);
if (value[0] == '\0') {
- get_metadata(metadata, "SolarElevation", value);
+ get_mtldata(mtldata, "SolarElevation", value);
}
lsat->sun_elev = atof(value);
if (lsat->sun_elev == 0.)
- G_warning("Sun elevation is %f", lsat->sun_elev);
+ G_warning("Sun elevation is %f", lsat->sun_elev);
- get_metadata(metadata, "SCENE_CENTER_TIME", value);
+ get_mtldata(mtldata, "SCENE_CENTER_TIME", value);
if (value[0] == '\0') {
- get_metadata(metadata, "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) {
@@ -175,13 +193,13 @@
{
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_meta], i, 1);
+ sprintf(key, fmt_gain[ver_mtl], i, 1);
if (i != 6)
- key[ver_meta == 0 ? 10 : 12] = '\0';
- get_metadata(metadata, key, value + j++);
+ key[ver_mtl == 0 ? 10 : 12] = '\0';
+ get_mtldata(mtldata, key, value + j++);
if (i == 6) {
- sprintf(key, fmt_gain[ver_meta], i, 2);
- get_metadata(metadata, key, value + j++);
+ sprintf(key, fmt_gain[ver_mtl], i, 2);
+ get_mtldata(mtldata, key, value + j++);
}
}
value[j] = '\0';
@@ -190,9 +208,8 @@
break;
}
case 8:
- set_LDCM(lsat);
+ set_OLI(lsat);
break;
-
default:
G_warning("Unable to recognize satellite platform [%d]",
lsat->number);
@@ -200,155 +217,168 @@
}
/* Update the information from metadata file, if infile */
- if (get_metadata == get_mtlformat) {
- if (ver_meta == 0) { /* now MTLold.txt */
+ if (get_mtldata == get_mtlformat) {
+ 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);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].lmax = atof(value);
sprintf(key, "LMIN_BAND%d", lsat->band[i].code);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].lmin = atof(value);
sprintf(key, "QCALMAX_BAND%d", lsat->band[i].code);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].qcalmax = atof(value);
sprintf(key, "QCALMIN_BAND%d", lsat->band[i].code);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].qcalmin = atof(value);
}
}
- else { /* now MTL.txt */
-
+ else {
G_verbose_message("Metada file is MTL file: new format");
- if (strstr(metadata, "RADIANCE_MAXIMUM_BAND") != NULL) {
+ /* 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_metadata(metadata, key, value);
+ 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_metadata(metadata, key, value);
+ 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_metadata(metadata, key, value);
+ 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_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].qcalmin = atof(value);
}
else {
sprintf(key, "RADIANCE_MAXIMUM_BAND_%d",
lsat->band[i].code);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].lmax = atof(value);
sprintf(key, "RADIANCE_MINIMUM_BAND_%d",
lsat->band[i].code);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].lmin = atof(value);
sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d",
lsat->band[i].code);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].qcalmax = atof(value);
sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d",
lsat->band[i].code);
- get_metadata(metadata, key, value);
+ 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(metadata, "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_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].gain = atof(value);
sprintf(key, fmt_radad[mode], lsat->band[i].code);
- get_metadata(metadata, key, value);
+ 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);
- get_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
lsat->band[i].K1 = atof(value);
sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
- get_metadata(metadata, key, value);
+ 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_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_metadata(metadata, key, value);
+ 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].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_metadata(metadata, "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++) {
sprintf(key, "Band%dGainSetting", lsat->band[i].code);
- get_metadata(metadata, key, value);
+ 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_metadata(metadata, key, value);
+ get_mtldata(mtldata, key, value);
if (value[0] == '\0') {
G_warning(key);
continue;
@@ -366,6 +396,5 @@
}
}
- fclose(f);
return;
}
Modified: grass/trunk/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat_set.c 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/landsat_set.c 2013-10-04 11:16:10 UTC (rev 57936)
@@ -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/trunk/imagery/i.landsat.toar/local_proto.h
===================================================================
--- grass/trunk/imagery/i.landsat.toar/local_proto.h 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/local_proto.h 2013-10-04 11:16:10 UTC (rev 57936)
@@ -17,6 +17,6 @@
void set_ETM(lsat_data *, char[]);
-void set_LDCM(lsat_data *);
+void set_OLI(lsat_data *);
#endif
Modified: grass/trunk/imagery/i.landsat.toar/main.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/main.c 2013-10-04 10:37:29 UTC (rev 57935)
+++ grass/trunk/imagery/i.landsat.toar/main.c 2013-10-04 11:16:10 UTC (rev 57936)
@@ -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]
* Adopted for GRASS 7 by Martin Landa <landa.martin gmail.com>
@@ -27,6 +27,8 @@
#include "local_proto.h"
+#define QCALMAX 65536 /* L1-7=256 but L8=65536 */
+
int main(int argc, char *argv[])
{
struct History history;
@@ -49,14 +51,14 @@
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;
int overwrite;
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 */
G_gisinit(argv[0]);
@@ -64,7 +66,7 @@
/* initialize module */
module = G_define_module();
module->description =
- _("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OT.");
+ _("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OLI");
G_add_keyword(_("imagery"));
G_add_keyword(_("radiometric conversion"));
G_add_keyword(_("radiance"));
@@ -100,29 +102,29 @@
sensor->key = "sensor";
sensor->type = TYPE_STRING;
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->description = _("Required only if 'metfile' not given (recommended for sanity)");
+ sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,oli8";
desc = NULL;
G_asprintf(&desc,
- "mss1;%s;mss2;%s;mss3;%s;mss4;%s;mss5;%s;tm4;%s;tm5;%s;tm7;%s;ot8;%s",
- _("Landsat-1 MSS"),
- _("Landsat-2 MSS"),
- _("Landsat-3 MSS"),
- _("Landsat-4 MSS"),
- _("Landsat-5 MSS"),
- _("Landsat-4 TM"),
- _("Landsat-5 TM"),
- _("Landsat-7 ETM+"), _("Landsat_8 OLI/TIRS"));
+ "mss1;%s;mss2;%s;mss3;%s;mss4;%s;mss5;%s;tm4;%s;tm5;%s;tm7;%s;oli8;%s",
+ _("Landsat-1 MSS"),
+ _("Landsat-2 MSS"),
+ _("Landsat-3 MSS"),
+ _("Landsat-4 MSS"),
+ _("Landsat-5 MSS"),
+ _("Landsat-4 TM"),
+ _("Landsat-5 TM"),
+ _("Landsat-7 ETM+"),
+ _("Landsat_8 OLI/TIRS"));
sensor->descriptions = desc;
- sensor->required = NO;
+ sensor->required = NO; /* perhaps YES for clarity */
sensor->guisection = _("Metadata");
metho = G_define_option();
metho->key = "method";
metho->type = TYPE_STRING;
metho->required = NO;
- metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
+ metho->options = "uncorrected,dos1,dos2,dos2b,dos3,dos4";
metho->label = _("Atmospheric correction method");
metho->description = _("Atmospheric correction method");
metho->answer = "uncorrected";
@@ -232,10 +234,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;
@@ -267,9 +269,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);
@@ -316,28 +319,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(_("Missing '%s' and/or '%s' for this satellite"),
+ 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)
@@ -354,13 +358,13 @@
set_MSS1(&lsat);
else
G_fatal_error(_("Unknown satellite type (defined by '%s')"),
- sensor->key);
+ sensorname);
}
- /*****************************************
- * ------------ PREPARATION --------------
- *****************************************/
- if (strcasecmp(metho->answer, "corrected") == 0)
+ /*****************************************
+ * ------------ PREPARATION --------------
+ *****************************************/
+ if (strcasecmp(metho->answer, "corrected") == 0) /* deleted 2013 */
method = CORRECTED;
else if (strcasecmp(metho->answer, "dos1") == 0)
method = DOS1;
@@ -376,21 +380,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;
+ 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);
if ((infd = Rast_open_old(band_in, "")) < 0)
@@ -399,6 +403,8 @@
G_set_window(&cellhd);
in_data_type = Rast_get_map_type(infd);
+ if (in_data_type < 0)
+ G_fatal_error(_("Unable to read data type of raster map <%s>"), band_in);
inrast = Rast_allocate_buf(in_data_type);
nrows = Rast_window_rows();
@@ -421,56 +427,62 @@
ptr = (void *)((DCELL *) inrast + col);
q = (int)*((DCELL *) ptr);
break;
+ default:
+ ptr = NULL;
+ q = -1.;
+ break;
}
if (!Rast_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], excluding DN > %d",
+ dn_dark[i], hist[dn_dark[i]], dn_mode[i], hist[dn_mode[i]],
+ dn_sat);
G_free(inrast);
Rast_close(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)));
+ 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",
@@ -478,16 +490,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) {
@@ -497,7 +508,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"),
@@ -515,8 +526,8 @@
}
/*****************************************
- * ------------ CALCULUS -----------------
- *****************************************/
+ * ------------ CALCULUS -----------------
+ *****************************************/
G_message(_("Calculating..."));
for (i = 0; i < lsat.bands; i++) {
@@ -539,6 +550,8 @@
}
in_data_type = Rast_get_map_type(infd);
+ if (in_data_type < 0)
+ G_fatal_error(_("Unable to read data type of raster map <%s>"), band_in);
Rast_get_cellhd(band_in, "", &cellhd);
/* set same size as original band raster */
@@ -551,17 +564,19 @@
if ((outfd = Rast_open_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 = Rast_allocate_buf(in_data_type);
outrast = Rast_allocate_buf(DCELL_TYPE);
nrows = Rast_window_rows();
ncols = Rast_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);
@@ -580,6 +595,10 @@
ptr = (void *)((DCELL *) inrast + col);
qcal = (double)((DCELL *) inrast)[col];
break;
+ default:
+ ptr = NULL;
+ qcal = -1.;
+ break;
}
if (Rast_is_null_value(ptr, in_data_type) ||
qcal < lsat.band[i].qcalmin) {
@@ -613,17 +632,20 @@
ref_mode = lsat_rad2ref(ref_mode, &lsat.band[i]);
}
- /* ================================================================= */
G_free(inrast);
+ G_free(outrast);
Rast_close(infd);
- G_free(outrast);
Rast_close(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 */
Rast_init_colors(&colors);
Rast_read_fp_range(band_out, G_mapset(), &range);
@@ -631,32 +653,36 @@
Rast_make_grey_scale_fp_colors(&colors, min, max);
Rast_write_colors(band_out, G_mapset(), &colors);
- /* Initialize the 'hist' structure with basic info */
+ /* Initialize the 'history' structure with basic info */
Rast_short_history(band_out, "raster", &history);
- /* Append a string to the 'history' structure */
Rast_append_format_history(&history,
" %s of Landsat-%d %s (method %s)",
- (frad->answer ? "Radiance" :
- (lsat.band[i].thermal ? "Temperature" :
- "Reflectance")),
+ (frad->
+ answer ? "Radiance" : (lsat.band[i].
+ thermal ?
+ "Temperature" :
+ "Reflectance")),
lsat.number, lsat.sensor, metho->answer);
Rast_append_history(&history,
- "----------------------------------------------------------------");
+ "-----------------------------------------------------------------");
Rast_append_format_history(&history,
- " Acquisition date ...................... %s",
- lsat.date);
+ " Acquisition date (and time) ........... %s (%.4lf h)",
+ lsat.date, lsat.time);
Rast_append_format_history(&history,
" Production date ....................... %s\n",
lsat.creation);
Rast_append_format_history(&history,
- " Earth-sun distance (d) ................ %.8lf",
+ " Earth-sun distance (d) ................ %.7lf",
lsat.dist_es);
Rast_append_format_history(&history,
+ " Sun elevation (and azimuth) ........... %.5lf (%.5lf)",
+ lsat.sun_elev, lsat.sun_az);
+ Rast_append_format_history(&history,
" Digital number (DN) range ............. %.0lf to %.0lf",
lsat.band[i].qcalmin,
lsat.band[i].qcalmax);
Rast_append_format_history(&history,
- " Calibration constants (Lmin to Lmax) .. %+.3lf to %+.3lf",
+ " Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
lsat.band[i].lmin, lsat.band[i].lmax);
Rast_append_format_history(&history,
" DN to Radiance (gain and bias) ........ %+.5lf and %+.5lf",
@@ -671,10 +697,10 @@
" Mean solar irradiance (ESUN) .......... %.3lf",
lsat.band[i].esun);
Rast_append_format_history(&history,
- " Reflectance = Radiance divided by ..... %.5lf",
+ " Radiance to Reflectance (divide by) ... %+.5lf",
lsat.band[i].K1);
if (method > DOS) {
- Rast_append_history(&history, " ");
+ Rast_append_format_history(&history, " ");
Rast_append_format_history(&history,
" Dark object (%4d pixels) DN = ........ %d",
pixel, dn_dark[i]);
@@ -684,7 +710,7 @@
}
}
Rast_append_history(&history,
- "-----------------------------------------------------------------");
+ "------------------------------------------------------------------");
Rast_command_history(&history);
Rast_write_history(band_out, &history);
@@ -692,7 +718,7 @@
if (lsat.band[i].thermal)
Rast_write_units(band_out, "Kelvin");
else if (frad->answer)
- Rast_write_units(band_out, "W/(m^2 sr µm)");
+ Rast_write_units(band_out, "W/(m^2 sr um)");
else
Rast_write_units(band_out, "unitless");
}
More information about the grass-commit
mailing list