[GRASS-SVN] r56905 - grass/trunk/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 24 08:34:19 PDT 2013
Author: mmetz
Date: 2013-06-24 08:34:19 -0700 (Mon, 24 Jun 2013)
New Revision: 56905
Modified:
grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html
grass/trunk/imagery/i.landsat.toar/landsat.c
grass/trunk/imagery/i.landsat.toar/landsat.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: sync to relbr6
Modified: grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html
===================================================================
--- grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/i.landsat.toar.html 2013-06-24 15:34:19 UTC (rev 56905)
@@ -8,15 +8,16 @@
<p>
Usually, to do so the production date, the acquisition date, and the
-solar elevation is needed. Moreover, for Landsat-7 ETM+ it is also
+solar elevation are needed. Moreover, for Landsat-7 ETM+ it is also
needed the gain (high or low) of the nine respective bands.
<p>
-Optionally (recommended), the data can be read from metadata file (.met or MTL.txt) for all
-Landsat MSS, TM and ETM+. However, if the solar elevation or the
-product creation date are given the values of the metadata file are
-overwriten. This is necessary when the data in the .met file is
-incorrect or not accurate.
+Optionally (recommended), the data can be read from metadata file
+(.met or MTL.txt) for all Landsat MSS, TM, ETM+ and OLI/TIRS. However,
+if the solar elevation is given the value of the metadata file are
+overwritten. This is necessary when the data in the .met file is
+incorrect or not accurate. Also, if acquisition or production dates are
+not found in the metadata file then the command line values are used.
<p>
<b>Attention</b>: Any null value or smaller than QCALmin in the input
@@ -138,6 +139,7 @@
<li>Landsat-5 MSS: April 6, 1984 and November 9, 1984</li>
<li>Landsat-5 TM: May 4, 2003 and April, 2 2007</li>
<li>Landsat-7 ETM+: July 1, 2000</li>
+ <li>Landsat-8 OLI/TIRS: launched in 2013</li>
</ul>
<h2>EXAMPLES</h2>
Modified: grass/trunk/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat.c 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/landsat.c 2013-06-24 15:34:19 UTC (rev 56905)
@@ -5,9 +5,9 @@
#include "landsat.h"
-#define PI 3.1415926535897932384626433832795
-#define R2D 57.295779513082320877
-#define D2R 0.017453292519943295769
+#define PI M_PI
+#define R2D 180. / M_PI
+#define D2R M_PI / 180.
/****************************************************************************
* PURPOSE: Calibrated Digital Number to at-satellite Radiance
@@ -22,7 +22,7 @@
*****************************************************************************/
double lsat_rad2ref(double rad, band_data * band)
{
- return (double)(rad / band->K2);
+ return (double)(rad / band->K1);
}
/****************************************************************************
@@ -62,10 +62,10 @@
sin_e = (double)(sin(D2R * lsat->sun_elev));
cos_v = (double)(cos(D2R * (lsat->number < 4 ? 9.2 : 8.2)));
- /** Global irradiance on the sensor.
- Radiance to reflectance coefficient, only NO thermal bands.
- K1 and K2 variables are also utilized as thermal constants
- */
+ /** Global irradiance on the sensor.
+ Radiance to reflectance coefficient, only NO thermal bands.
+ K1 and K2 variables are also utilized as thermal constants
+ */
if (lsat->band[i].thermal == 0) {
switch (method) {
case DOS2:
@@ -129,16 +129,17 @@
break;
}
rad_sun = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
- G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz,
- Edown);
+ if (method > DOS)
+ G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
+ TAUv, TAUz, Edown);
- lsat->band[i].K1 = 0.;
- lsat->band[i].K2 = rad_sun;
+ lsat->band[i].K1 = rad_sun;
+ lsat->band[i].K2 = 0.;
}
- /** Digital number to radiance coefficients.
- Whitout atmospheric calibration for thermal bands.
- */
+ /** Digital number to radiance coefficients.
+ Without atmospheric calibration for thermal bands.
+ */
lsat->band[i].gain = ((lsat->band[i].lmax - lsat->band[i].lmin) /
(lsat->band[i].qcalmax - lsat->band[i].qcalmin));
Modified: grass/trunk/imagery/i.landsat.toar/landsat.h
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat.h 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/landsat.h 2013-06-24 15:34:19 UTC (rev 56905)
@@ -30,14 +30,14 @@
double wavemax, wavemin; /* Wavelength in µm */
+ double esun; /* Mean solar irradiance */
double lmax, lmin; /* Spectral radiance */
double qcalmax, qcalmin; /* Quantized calibrated pixel */
- double esun; /* Mean solar irradiance */
char thermal; /* Flag to thermal band */
double gain, bias; /* Gain and Bias of sensor */
- double K1, K2; /* Thermal calibration constants,
- or Rad2Ref constants */
+ double K1, K2; /* Thermal calibration or
+ Rad2Ref constants */
} band_data;
@@ -48,9 +48,10 @@
char creation[11]; /* Image production date */
char date[11]; /* Image acquisition date */
+
double dist_es; /* Distance Earth-Sun */
- double sun_elev; /* Solar elevation */
- double sunaz; /* Solar Azimuth */
+ double sun_elev; /* Sun elevation */
+ double sun_az; /* Sun Azimuth */
double time; /* Image Acquisition Time */
char sensor[10]; /* Type of sensor: MSS, TM, ETM+, OLI/TIRS */
Modified: grass/trunk/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat_met.c 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/landsat_met.c 2013-06-24 15:34:19 UTC (rev 56905)
@@ -1,6 +1,6 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
#include <math.h>
#include <grass/gis.h>
@@ -9,44 +9,55 @@
#include "local_proto.h"
#include "earth_sun.h"
+#define PI M_PI
+#define D2R M_PI / 180.
+
#define MAX_STR 127
-#define METADATA_SIZE 65535 /* MTL.txt file size 65535 bytes */
-#define TM5_MET_SIZE 28700 /* .met file size 28686 bytes */
+#define METADATA_SIZE 65535 /* MTL.txt file size 65535 bytes */
+#define TM5_MET_SIZE 28700 /* .met file size 28686 bytes */
inline void chrncpy(char *dest, char src[], int n)
{
int i;
- for( i = 0 ; i < n && src[i] != '\0' && src[i] != '\"'; i++) dest[i] = src[i];
+
+ for (i = 0; i < n && src[i] != '\0' && src[i] != '\"'; i++)
+ dest[i] = src[i];
dest[i] = '\0';
}
-void get_oldformat( const char metadata[], char * key, char value[] )
+void get_metformat(const char metadata[], char *key, char value[])
{
int i = 0;
- char * ptr = strstr(metadata, key);
- if (ptr != NULL)
- {
- ptr = strstr(ptr, " VALUE ");
- if (ptr != NULL)
- {
+ char *ptrmet = strstr(metadata, key);
+ char *ptr;
+
+ if (ptrmet != NULL) {
+ ptr = strstr(ptrmet, " VALUE ");
+ if (ptr != NULL) {
while (*ptr++ != '=') ;
- while (*ptr <= ' ' || *ptr == '\"') *ptr++;
- while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ') value[i++] = *ptr++;
+ while (*ptr <= ' ' || *ptr == '\"')
+ ptr++;
+ while (i < MAX_STR && *ptr != '\"' && *ptr >= ' ')
+ value[i++] = *ptr++;
}
}
value[i] = '\0';
}
-void get_newformat( const char metadata[], char * key, char value[] )
+void get_mtlformat(const char metadata[], char *key, char value[])
{
int i = 0;
- char * ptr = strstr(metadata, key);
- if (ptr != NULL)
- {
+ char *ptr = strstr(metadata, key);
+
+ if (ptr != NULL) {
while (*ptr++ != '=') ;
- while (*ptr <= ' ' || *ptr == '\"') *ptr++;
- while (i < MAX_STR && *ptr != '\"' && *ptr > ' ') value[i++] = *ptr++;
+
+ while (*ptr <= ' ' || *ptr == '\"')
+ ptr++;
+
+ while (i < MAX_STR && *ptr != '\"' && *ptr > ' ')
+ value[i++] = *ptr++;
}
value[i] = '\0';
}
@@ -56,160 +67,305 @@
* PURPOSE: Read parameters from Landsat metadata files
*****************************************************************************/
-void lsat_metadata( char * metafile, lsat_data * lsat)
+void lsat_metadata(char *metafile, lsat_data * lsat)
{
FILE *f;
- char mtldata[METADATA_SIZE];
+ char metadata[METADATA_SIZE];
char key[MAX_STR], value[MAX_STR];
- void (*get_mtldata)( const char [], char *, char [] );
- int i, j, version;
+ void (*get_metadata) (const char [], char *, char []);
+ int i, j, ver_meta;
if ((f = fopen(metafile, "r")) == NULL)
- G_fatal_error(_("Metadata file <%s> not found"), metafile);
+ G_fatal_error(_("Metadata file <%s> not found"), metafile);
- i = fread(mtldata, METADATA_SIZE, 1, f);
- get_mtldata = (strstr(mtldata, " VALUE ") != NULL) ? get_oldformat : get_newformat;
- version = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
+ i = fread(metadata, METADATA_SIZE, 1, f);
+
+ /* 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_mtldata(mtldata, "SPACECRAFT_ID", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "PLATFORMSHORTNAME", value);
+ get_metadata(metadata, "SPACECRAFT_ID", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "PLATFORMSHORTNAME", value);
}
- i = 0; while( value[i] && (value[i] < '0' || value[i] > '9') ) i++;
+ i = 0;
+ while (value[i] && (value[i] < '0' || value[i] > '9'))
+ i++;
lsat->number = atoi(value + i);
- get_mtldata(mtldata, "SENSOR_ID", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "SENSORSHORTNAME", value);
+ get_metadata(metadata, "SENSOR_ID", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "SENSORSHORTNAME", value);
}
chrncpy(lsat->sensor, value, 8);
- get_mtldata(mtldata, "DATE_ACQUIRED", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "ACQUISITION_DATE", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "CALENDARDATE", value);
+ get_metadata(metadata, "DATE_ACQUIRED", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "ACQUISITION_DATE", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "CALENDARDATE", value);
}
}
chrncpy(lsat->date, value, 10);
- get_mtldata(mtldata, "FILE_DATE", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "CREATION_TIME", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "PRODUCTIONDATETIME", value);
+ get_metadata(metadata, "FILE_DATE", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "CREATION_TIME", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "PRODUCTIONDATETIME", value);
}
}
- chrncpy(lsat->creation, value, 10);
+ if (value[0] != '\0')
+ chrncpy(lsat->creation, value, 10);
+ else
+ G_warning
+ ("Using production date from the command line 'product_date'");
- get_mtldata(mtldata, "SUN_AZIMUTH", value);
- lsat->sunaz = atof(value);
- if( lsat->sunaz == 0. )
- G_warning("Sun azimuth is %f", lsat->sunaz);
+ get_metadata(metadata, "SUN_AZIMUTH", value);
+ lsat->sun_az = atof(value);
+ if (lsat->sun_az == 0.)
+ G_warning("Sun azimuth is %f", lsat->sun_az);
- get_mtldata(mtldata, "SUN_ELEVATION", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "SolarElevation", value);
+ get_metadata(metadata, "SUN_ELEVATION", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "SolarElevation", value);
}
lsat->sun_elev = atof(value);
- if( lsat->sun_elev == 0. )
+ if (lsat->sun_elev == 0.)
G_warning("Sun elevation is %f", lsat->sun_elev);
- get_mtldata(mtldata, "SCENE_CENTER_TIME", value);
- if( value[0] == '\0' )
- {
- get_mtldata(mtldata, "SCENE_CENTER_SCAN_TIME", value);
+ get_metadata(metadata, "SCENE_CENTER_TIME", value);
+ if (value[0] == '\0') {
+ get_metadata(metadata, "SCENE_CENTER_SCAN_TIME", value);
}
- /*Thanks Markus Metz !
- Remove trailing 'z'*/
+ /* Remove trailing 'z'*/
value[strlen(value) - 1]='\0';
/* Cast from hh:mm:ss into hh.hhh*/
G_llres_scan(value, &lsat->time);
- if( lsat->time == 0. )
+ if (lsat->time == 0.)
G_warning("Time is %f", lsat->time);
- /* Fill data with the basic sensor parameters */
- switch(lsat->number)
- {
- case 1:
- set_MSS1(lsat);
- break;
- case 2:
- set_MSS2(lsat);
- break;
- case 3:
- set_MSS3(lsat);
- break;
- case 4:
- if (lsat->sensor[0] == 'M')
- set_MSS4(lsat);
- else
- set_TM4(lsat);
- break;
- case 5:
- if (lsat->sensor[0] == 'M')
- set_MSS5(lsat);
- else
- set_TM5(lsat);
- break;
- case 7:
- {
- char * fmt_gain[] = { "BAND%d_GAIN%d ", " GAIN_BAND_%d_VCID_%d" };
- for( i = 1, j = 0; i < 9; i++ )
- {
- sprintf(key, fmt_gain[version], i, 1);
- if (i != 6) key[version == 0 ? 10 : 12] = '\0';
- get_mtldata(mtldata, key, value + j++);
- if (i == 6)
- {
- sprintf(key, fmt_gain[version], i, 2);
- get_mtldata(mtldata, key, value + j++);
- }
- }
- value[j] = '\0';
- G_debug(1, "ETM+ gain = [%s]", value);
- set_ETM(lsat, value);
- break;
- }
- default:
- G_warning("Unable to recognize satellite platform [%d]", lsat->number);
- break;
+ /* Fill the data with the basic/default sensor parameters */
+ switch (lsat->number) {
+ case 1:
+ set_MSS1(lsat);
+ break;
+ case 2:
+ set_MSS2(lsat);
+ break;
+ case 3:
+ set_MSS3(lsat);
+ break;
+ case 4:
+ if (lsat->sensor[0] == 'M')
+ set_MSS4(lsat);
+ else
+ set_TM4(lsat);
+ break;
+ case 5:
+ if (lsat->sensor[0] == 'M')
+ set_MSS5(lsat);
+ else
+ set_TM5(lsat);
+ break;
+ case 7:
+ {
+ char *fmt_gain[] = { "BAND%d_GAIN%d", " GAIN_BAND_%d_VCID_%d" };
+ for (i = 1, j = 0; i < 9; i++) {
+ sprintf(key, fmt_gain[ver_meta], i, 1);
+ if (i != 6)
+ key[ver_meta == 0 ? 10 : 12] = '\0';
+ get_metadata(metadata, key, value + j++);
+ if (i == 6) {
+ sprintf(key, fmt_gain[ver_meta], i, 2);
+ get_metadata(metadata, key, value + j++);
+ }
+ }
+ value[j] = '\0';
+ G_debug(1, "ETM+ gain = [%s]", value);
+ set_ETM(lsat, value);
+ break;
+ }
+ case 8:
+ set_LDCM(lsat);
+ break;
+
+ default:
+ G_warning("Unable to recognize satellite platform [%d]",
+ lsat->number);
+ break;
}
/* Update the information from metadata file, if infile */
- /* if( format == NEW_FORMAT ) */
- if( get_mtldata == get_newformat )
- {
- char * fmt_lmax[] = { "LMAX_BAND%d", "RADIANCE_MAXIMUM_BAND_%d" };
- char * fmt_lmin[] = { "LMIN_BAND%d", "RADIANCE_MINIMUM_BAND_%d" };
- char * fmt_qcalmax[] = { "QCALMAX_BAND%d", "QUANTIZE_CAL_MAX_BAND_%d" };
- char * fmt_qcalmin[] = { "QCALMIN_BAND%d", "QUANTIZE_CAL_MIN_BAND_%d" };
+ if (get_metadata == get_mtlformat) {
+ if (ver_meta == 0) { /* now MTLold.txt */
+ G_verbose_message("Metada file is MTL file: old format");
+ for (i = 0; i < lsat->bands; i++) {
+ sprintf(key, "LMAX_BAND%d", lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].lmax = atof(value);
+ sprintf(key, "LMIN_BAND%d", lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].lmin = atof(value);
+ sprintf(key, "QCALMAX_BAND%d", lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].qcalmax = atof(value);
+ sprintf(key, "QCALMIN_BAND%d", lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].qcalmin = atof(value);
+ }
+ }
+ else { /* now MTL.txt */
- for (i = 0; i < lsat->bands; i++) {
- sprintf(key, fmt_lmax[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].lmax = atof(value);
- sprintf(key, fmt_lmin[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].lmin = atof(value);
- sprintf(key, fmt_qcalmax[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].qcalmax = atof(value);
- sprintf(key, fmt_qcalmin[version], lsat->band[i].code);
- get_mtldata(mtldata, key, value);
- lsat->band[i].qcalmin = atof(value);
- }
+ G_verbose_message("Metada file is MTL file: new format");
+ if (strstr(metadata, "RADIANCE_MAXIMUM_BAND") != NULL) {
+ G_verbose_message
+ ("RADIANCE & QUANTIZE from data of the metadata file");
+ for (i = 0; i < lsat->bands; i++) {
+ if (lsat->band[i].thermal && lsat->number == 7) {
+ sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
+ lsat->band[i].code - 60);
+ get_metadata(metadata, 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);
+ 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);
+ 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);
+ lsat->band[i].qcalmin = atof(value);
+ }
+ else {
+ sprintf(key, "RADIANCE_MAXIMUM_BAND_%d",
+ lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].lmax = atof(value);
+ sprintf(key, "RADIANCE_MINIMUM_BAND_%d",
+ lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].lmin = atof(value);
+ sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d",
+ lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].qcalmax = atof(value);
+ sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d",
+ lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].qcalmin = atof(value);
+ }
+ }
+ }
+ else {
+ G_verbose_message
+ ("RADIANCE & QUANTIZE from radiometric rescaling group of the metadata file");
+ /* from LDCM sample file: mode = 0; from LDCM-DFCB-004.pdf file: mode = 1 */
+ int mode =
+ (strstr(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);
+ lsat->band[i].gain = atof(value);
+ sprintf(key, fmt_radad[mode], lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].bias = atof(value);
+ /* reverse to calculate the original values */
+ lsat->band[i].qcalmax =
+ (lsat->number < 8 ? 255. : 65535.);
+ lsat->band[i].qcalmin = 1.;
+ lsat->band[i].lmin =
+ lsat->band[i].gain * lsat->band[i].qcalmin +
+ lsat->band[i].bias;
+ lsat->band[i].lmax =
+ lsat->band[i].gain * lsat->band[i].qcalmax +
+ lsat->band[i].bias;
+ /* ----- */
+ if (lsat->number == 8) {
+ if (lsat->band[i].thermal) {
+ sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].K1 = atof(value);
+ sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].K2 = atof(value);
+ }
+ else {
+ lsat->band[i].K1 = 0.;
+ /* ESUN from metadafa file: bias/K2 seem better than gain/K1 */
+ sprintf(key, fmt_refad[mode], lsat->band[i].code);
+ get_metadata(metadata, key, value);
+ lsat->band[i].K2 = atof(value);
+ lsat->band[i].esun =
+ (double)(PI * lsat->dist_es * lsat->dist_es *
+ lsat->band[i].bias) /
+ (sin(D2R * lsat->sun_elev) *
+ lsat->band[i].K2);
+ }
+ }
+ }
+ }
+ /* Other possible values in file */
+ get_metadata(metadata, "EARTH_SUN_DISTANCE", value);
+ if (value[0] != '\0') {
+ lsat->dist_es = atof(value);
+ G_verbose_message
+ ("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
+ }
+ }
}
+ else {
+ G_verbose_message("Metada file is MET file");
+ G_verbose_message
+ ("RADIANCE & QUANTIZE from band setting of the metadata file");
+ for (i = 0; i < lsat->bands; i++) {
+ sprintf(key, "Band%dGainSetting", lsat->band[i].code);
+ get_metadata(metadata, 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);
+ if (value[0] == '\0') {
+ G_warning(key);
+ continue;
+ }
+ lsat->band[i].bias = atof(value);
- (void)fclose(f);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 1.;
+ lsat->band[i].lmin =
+ lsat->band[i].gain * lsat->band[i].qcalmin +
+ lsat->band[i].bias;
+ lsat->band[i].lmax =
+ lsat->band[i].gain * lsat->band[i].qcalmax +
+ lsat->band[i].bias;
+ }
+ }
+
+ fclose(f);
return;
}
Modified: grass/trunk/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/landsat_set.c 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/landsat_set.c 2013-06-24 15:34:19 UTC (rev 56905)
@@ -1,6 +1,6 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
#include <math.h>
#include <grass/gis.h>
@@ -15,10 +15,10 @@
/* green, red, near infrared, near infrared */
int band[] = { 1, 2, 3, 4 };
- int code[] = { 4, 5, 6, 7 };
+ int code[] = { 4, 5, 6, 7 }; /* corrected for MSS4 and MSS5 to 1,2,3,4 */
double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
- /* original: 79x57, now all 60 */
+ /* original: 79x57, now all 60 m */
strcpy(lsat->sensor, "MSS");
@@ -46,7 +46,7 @@
/* 30, 30, 30, 30, 30, 120 original, 60 resamples before Feb 25, 2010 and 30 after, 30 */
if (!lsat->sensor)
- strcpy(lsat->sensor, "TM");
+ strcpy(lsat->sensor, "TM");
lsat->bands = 7;
for (i = 0; i < lsat->bands; i++) {
@@ -65,11 +65,11 @@
{
int i;
- /* blue, green red, near infrared, shortwave IR, thermal IR, shortwave IR, panchromatic */
+ /* 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, 2.09, 0.52 };
- double wmax[] = { 0.515, 0.605, 0.690, 0.90, 1.75, 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 +87,7 @@
return;
}
-void sensor_OLI(lsat_data * lsat)
+void sensor_LDCM(lsat_data * lsat)
{
int i;
@@ -103,13 +103,13 @@
lsat->bands = 11;
for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(code + i);
- lsat->band[i].wavemax = *(wmax + i);
- lsat->band[i].wavemin = *(wmin + i);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 1.;
- lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = *(code + i);
+ lsat->band[i].wavemax = *(wmax + i);
+ lsat->band[i].wavemin = *(wmin + i);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 1.;
+ lsat->band[i].thermal = (lsat->band[i].number > 9 ? 1 : 0);
}
return;
}
@@ -180,9 +180,9 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("1975-07-16"))
- i = 0;
+ i = 0;
else
- i = 1;
+ i = 1;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -279,11 +279,11 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("1982-08-26"))
- i = 0;
+ i = 0;
else if (julian < julian_char("1983-03-31"))
- i = 1;
+ i = 1;
else
- i = 2;
+ i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -329,11 +329,11 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("1983-08-01"))
- i = 0;
+ i = 0;
else if (julian < julian_char("1984-01-15"))
- i = 1;
+ i = 1;
else
- i = 2;
+ i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -436,11 +436,11 @@
julian = julian_char(lsat->creation);
if (julian < julian_char("2003-05-04"))
- i = 0;
+ i = 0;
else if (julian < julian_char("2007-04-02"))
- i = 1;
+ i = 1;
else
- i = 2;
+ i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
@@ -453,9 +453,9 @@
}
jbuf = julian_char("2004-04-04");
- if (julian >= jbuf && !(lsat->flag & METADATAFILE) )
- {
- G_warning("Using QCalMin=1.0 as a NLAPS product processed after 04/04/2004");
+ if (julian >= jbuf && !(lsat->flag & METADATAFILE)) {
+ G_warning
+ ("Using QCalMin=1.0 as a NLAPS product processed after 04/04/2004");
}
lsat->number = 5;
sensor_TM(lsat);
@@ -463,9 +463,9 @@
lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
if (julian >= jbuf)
lsat->band[i].qcalmin = 1.;
+ j = lsat->band[i].number - 1;
lsat->band[i].esun = *(esun + j);
lsat->band[i].lmax = *(lmax + j);
lsat->band[i].lmin = *(lmin + j);
@@ -528,7 +528,6 @@
for (i = 0; i < lsat->bands; i++) {
j = lsat->band[i].number - 1;
- lsat->band[i].esun = *(esun + j);
if (gain[i] == 'H' || gain[i] == 'h') {
lmax = LmaxH[k];
lmin = LminH[k];
@@ -537,6 +536,7 @@
lmax = LmaxL[k];
lmin = LminL[k];
}
+ lsat->band[i].esun = *(esun + j);
lsat->band[i].lmax = *(lmax + j);
lsat->band[i].lmin = *(lmin + j);
if (lsat->band[i].thermal) {
@@ -552,42 +552,42 @@
* PURPOSE: Store values of Landsat-8 OLI/TIRS
* February 14, 2013
*****************************************************************************/
-void set_OLI(lsat_data * lsat)
+void set_LDCM(lsat_data * lsat)
{
int i, j;
- double julian, *lmax, *lmin;
+ double *lmax, *lmin;
- /* Spectral radiances at detector
- double Lmax[][4] = {
- {220., 175., 145., 147.},
- {259., 179., 149., 128.}
+ /* Spectral radiances at detector */
+
+ /* uncorrected values */
+ double Lmax[][11] = {
+ {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
};
- double Lmin[][4] = {
- {4., 3., 3., 1.},
- {4., 3., 3., 1.}
+ double Lmin[][11] = {
+ {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}
};
- * Solar exoatmospheric spectral irradiances
- double esun[] = { 1824., 1570., 1249., 853.4 };
+ /* Solar exoatmospheric spectral irradiances */
+ /* estimates */
+ double esun[] = { 2062., 21931., 1990., 1688., 1037., 268.6, 94.6, 1892., 399.0, 0., 0. };
- lmax = Lmax[i];
- lmin = Lmin[i];
+ lmax = Lmax[0];
+ lmin = Lmin[0];
- lsat->number = 3;
- sensor_MSS(lsat);
+ lsat->number = 8;
+ sensor_LDCM(lsat);
lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
- lsat->band[i].esun = *(esun + j);
- lsat->band[i].lmax = *(lmax + j);
- lsat->band[i].lmin = *(lmin + j);
- if (lsat->band[i].thermal) {
- lsat->band[i].K1 = 0;
- lsat->band[i].K2 = 0;
- }
+ j = lsat->band[i].number - 1;
+ lsat->band[i].esun = *(esun + j);
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
+ if (lsat->band[i].thermal) {
+ lsat->band[i].K1 = (lsat->band[i].number == 10 ? 666.09 : 697.56);
+ lsat->band[i].K2 = 1282.71;
+ }
}
- */
- G_debug(1, "Landsat-8 LDCM");
+ G_debug(1, "Landsat-8 OLI/TIRS");
return;
}
Modified: grass/trunk/imagery/i.landsat.toar/local_proto.h
===================================================================
--- grass/trunk/imagery/i.landsat.toar/local_proto.h 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/local_proto.h 2013-06-24 15:34:19 UTC (rev 56905)
@@ -17,4 +17,6 @@
void set_ETM(lsat_data *, char[]);
+void set_LDCM(lsat_data *);
+
#endif
Modified: grass/trunk/imagery/i.landsat.toar/main.c
===================================================================
--- grass/trunk/imagery/i.landsat.toar/main.c 2013-06-24 10:15:20 UTC (rev 56904)
+++ grass/trunk/imagery/i.landsat.toar/main.c 2013-06-24 15:34:19 UTC (rev 56905)
@@ -9,9 +9,9 @@
* Adopted for GRASS 7 by Martin Landa <landa.martin gmail.com>
*
* PURPOSE: Calculate TOA Radiance or Reflectance and Kinetic Temperature
- * for Landsat 1/2/3/4/5 MS, 4/5 TM or 7 ETM+
+ * for Landsat 1/2/3/4/5 MS, 4/5 TM, 7 ETM+, and 8 OLI/TIRS
*
- * COPYRIGHT: (C) 2006-2012 by the GRASS Development Team
+ * COPYRIGHT: (C) 2006-2013 by the GRASS Development Team
*
* This program is free software under the GNU General
* Public License (>=v2). Read the file COPYING that
@@ -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, *lsatmet;
+ struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate,
+ *pdate, *elev, *bgain, *metho, *perc, *dark, *atmo, *lsatmet;
char *inputname, *met, *outputname, *sensorname;
- struct Flag *frad, *print_meta;
+ struct Flag *frad, *print_meta, *named;
lsat_data lsat;
char band_in[GNAME_MAX], band_out[GNAME_MAX];
@@ -64,7 +64,7 @@
/* initialize 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+/OT.");
G_add_keyword(_("imagery"));
G_add_keyword(_("radiometric conversion"));
G_add_keyword(_("radiance"));
@@ -85,7 +85,8 @@
output_prefix = G_define_option();
output_prefix->key = "output_prefix";
output_prefix->label = _("Prefix for output raster maps");
- output_prefix->description = _("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
+ output_prefix->description =
+ _("Example: 'B.toar.' generates B.toar.1, B.toar.2, ...");
output_prefix->type = TYPE_STRING;
output_prefix->required = YES;
@@ -99,11 +100,11 @@
sensor->key = "sensor";
sensor->type = TYPE_STRING;
sensor->label = _("Spacecraft sensor");
- sensor->description = _("Required only if 'metfile' not given");
- sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7";
+ sensor->description = _("Required only if 'metfile' not given (recommended by sanity)");
+ sensor->options = "mss1,mss2,mss3,mss4,mss5,tm4,tm5,tm7,ot8";
desc = NULL;
G_asprintf(&desc,
- "mss1;%s;mss2;%s;mss3;%s;mss4;%s;mss5;%s;tm4;%s;tm5;%s;tm7;%s",
+ "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"),
@@ -111,7 +112,8 @@
_("Landsat-5 MSS"),
_("Landsat-4 TM"),
_("Landsat-5 TM"),
- _("Landsat-7 ETM+"));
+ _("Landsat-7 ETM+"),
+ _("Landsat_8 OLI/TIRS"));
sensor->descriptions = desc;
sensor->required = NO;
sensor->guisection = _("Metadata");
@@ -156,7 +158,7 @@
bgain->key = "gain";
bgain->type = TYPE_STRING;
bgain->required = NO;
- bgain->label = _("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
+ bgain->label = _("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
bgain->description = _("Required only if 'metfile' not given");
bgain->guisection = _("Settings");
@@ -164,7 +166,8 @@
perc->key = "percent";
perc->type = TYPE_DOUBLE;
perc->required = NO;
- perc->description = _("Percent of solar radiance in path radiance");
+ perc->label = _("Percent of solar radiance in path radiance");
+ perc->description = _("Required only if 'method' is any DOS");
perc->answer = "0.01";
perc->guisection = _("Settings");
@@ -172,8 +175,9 @@
dark->key = "pixel";
dark->type = TYPE_INTEGER;
dark->required = NO;
- dark->description =
+ dark->label =
_("Minimum pixels to consider digital number as dark object");
+ dark->description = _("Required only if 'method' is any DOS");
dark->answer = "1000";
dark->guisection = _("Settings");
@@ -181,7 +185,8 @@
atmo->key = "rayleigh";
atmo->type = TYPE_DOUBLE;
atmo->required = NO;
- atmo->description = _("Rayleigh atmosphere (diffuse sky irradiance)"); /* scattering coefficient? */
+ atmo->label = _("Rayleigh atmosphere (diffuse sky irradiance)"); /* scattering coefficient? */
+ atmo->description = _("Required only if 'method' is DOS3");
atmo->answer = "0.0";
atmo->guisection = _("Settings");
@@ -189,6 +194,7 @@
lsatmet->key = "lsatmet";
lsatmet->type = TYPE_STRING;
lsatmet->required = NO;
+ lsatmet->multiple = YES;
lsatmet->label = _("return value stored for a given metadata");
lsatmet->description = _("Required only if 'metfile' and -p given");
lsatmet->options = "number,creation,date,sun_elev,sensor,bands,sunaz,time";
@@ -209,8 +215,14 @@
/* define the different flags */
frad = G_define_flag();
frad->key = 'r';
- frad->description = _("Output at-sensor radiance for all bands");
+ frad->description =
+ _("Output at-sensor radiance instead of reflectance for all bands");
+ named = G_define_flag();
+ named->key = 'n';
+ named->description =
+ _("Input raster maps use as extension the number of the band instead the code");
+
/* define the different flags */
print_meta = G_define_flag();
print_meta->key = 'p';
@@ -228,7 +240,7 @@
met = metfn->answer;
inputname = input_prefix->answer;
outputname = output_prefix->answer;
- sensorname = sensor->answer ? sensor->answer: "";
+ sensorname = sensor->answer ? sensor->answer : "";
overwrite = G_check_overwrite(argc, argv);
@@ -251,53 +263,55 @@
}
lsat.sun_elev = elev->answer == NULL ? 0. : atof(elev->answer);
+
percent = atof(perc->answer);
pixel = atoi(dark->answer);
rayleigh = atof(atmo->answer);
/* Data from metadata file */
- /* Unnecessary because G_zero filled, but sanity */
+ /* Unnecessary because G_zero filled, but for sanity */
lsat.flag = NOMETADATAFILE;
- if (met != NULL)
- {
+ if (met != NULL) {
lsat.flag = METADATAFILE;
- lsat_metadata( met, &lsat );
- if(print_meta->answer) {
- if (lsatmet->answer == NULL) {
- G_fatal_error(_("Please use a metadata keyword with -p"));
- }
- if (strcmp(lsatmet->answer, "number") == 0) {
- fprintf(stdout,"%d\n",lsat.number);
- }
- if (strcmp(lsatmet->answer, "creation") == 0) {
- fprintf(stdout,"%s\n",lsat.creation);
- }
- if (strcmp(lsatmet->answer, "date") == 0) {
- fprintf(stdout,"%s\n",lsat.date);
- }
- if (strcmp(lsatmet->answer, "sun_elev") == 0) {
- fprintf(stdout,"%f\n",lsat.sun_elev);
- }
- if (strcmp(lsatmet->answer, "sensor") == 0) {
- fprintf(stdout,"%s\n",lsat.sensor);
- }
- if (strcmp(lsatmet->answer, "bands") == 0) {
- fprintf(stdout,"%d\n",lsat.bands);
- }
- if (strcmp(lsatmet->answer, "sunaz") == 0) {
- fprintf(stdout,"%f\n",lsat.sunaz);
- }
- if (strcmp(lsatmet->answer, "time") == 0) {
- fprintf(stdout,"%f\n",lsat.time);
- }
- exit(EXIT_SUCCESS);
+ lsat_metadata(met, &lsat);
+ if (print_meta->answer) {
+ if (lsatmet->answer == NULL) {
+ G_fatal_error(_("Please use a metadata keyword with -p"));
+ }
+ if (strcmp(lsatmet->answer, "number") == 0) {
+ fprintf(stdout,"%d\n",lsat.number);
+ }
+ if (strcmp(lsatmet->answer, "creation") == 0) {
+ fprintf(stdout,"%s\n",lsat.creation);
+ }
+ if (strcmp(lsatmet->answer, "date") == 0) {
+ fprintf(stdout,"%s\n",lsat.date);
+ }
+ if (strcmp(lsatmet->answer, "sun_elev") == 0) {
+ fprintf(stdout,"%f\n",lsat.sun_elev);
+ }
+ if (strcmp(lsatmet->answer, "sunaz") == 0) {
+ fprintf(stdout,"%f\n",lsat.sun_az);
+ }
+ if (strcmp(lsatmet->answer, "sensor") == 0) {
+ fprintf(stdout,"%s\n",lsat.sensor);
+ }
+ if (strcmp(lsatmet->answer, "bands") == 0) {
+ fprintf(stdout,"%d\n",lsat.bands);
+ }
+ if (strcmp(lsatmet->answer, "time") == 0) {
+ fprintf(stdout,"%f\n",lsat.time);
+ }
+ exit(EXIT_SUCCESS);
}
- G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number, lsat.sensor);
+ G_debug(1, "lsat.number = %d, lsat.sensor = [%s]",
+ lsat.number, lsat.sensor);
- if (!lsat.sensor || lsat.number > 7 || lsat.number < 1)
+ 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]", lsat.number, lsat.sensor, met);
+ G_debug(1, "Landsat-%d %s with data set in met file [%s]",
+ lsat.number, lsat.sensor, met);
/* Overwrite solar elevation of metadata file */
if (elev->answer != NULL)
@@ -305,7 +319,7 @@
}
/* Data from date and solar elevation */
else if (adate->answer == NULL || elev->answer == NULL) {
- G_fatal_error(_("Lacking '%s' or '%s' for this satellite"),
+ G_fatal_error(_("Missing '%s' and/or '%s' for this satellite"),
adate->key, elev->key);
}
else {
@@ -316,6 +330,8 @@
set_ETM(&lsat, bgain->answer);
}
/* Not need gain */
+ else if (strcmp(sensorname, "ot8") == 0)
+ set_LDCM(&lsat);
else if (strcmp(sensorname, "tm5") == 0)
set_TM5(&lsat);
else if (strcmp(sensorname, "tm4") == 0)
@@ -331,12 +347,13 @@
else if (strcmp(sensorname, "mss1") == 0)
set_MSS1(&lsat);
else
- G_fatal_error(_("Unknown satellite type (defined by '%s')"), sensor->key);
+ G_fatal_error(_("Unknown satellite type (defined by '%s')"),
+ sensor->key);
}
- /*****************************************
- * ------------ PREPARATION --------------
- *****************************************/
+ /*****************************************
+ * ------------ PREPARATION --------------
+ *****************************************/
if (strcasecmp(metho->answer, "corrected") == 0)
method = CORRECTED;
else if (strcasecmp(metho->answer, "dos1") == 0)
@@ -406,14 +423,14 @@
}
/* DN of dark object */
for (j = lsat.band[i].qcalmin; j < 256; j++) {
- if (hist[j] >= (unsigned int) pixel) {
+ if (hist[j] >= (unsigned int)pixel) {
dn_dark[i] = j;
break;
}
}
/* Mode of DN */
h_max = 0L;
- for (j = lsat.band[i].qcalmin; j < 241; j++) { /* Exclude ptentially saturated < 240 */
+ for (j = lsat.band[i].qcalmin; j < 241; j++) { /* Exclude potentially saturated < 240 */
/* G_debug(5, "%d-%ld", j, hist[j]); */
if (hist[j] > h_max) {
h_max = hist[j];
@@ -433,21 +450,24 @@
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);
+ */
if (G_verbose() > G_verbose_std()) {
- fprintf(stderr, " LANDSAT: %d SENSOR: %s\n", lsat.number, lsat.sensor);
+ fprintf(stderr, " LANDSAT: %d SENSOR: %s\n", lsat.number,
+ lsat.sensor);
fprintf(stderr, " ACQUISITION DATE %s [production date %s]\n",
lsat.date, lsat.creation);
- fprintf(stderr, " earth-sun distance = %.8lf\n", lsat.dist_es);
- fprintf(stderr, " solar elevation angle = %.8lf\n", lsat.sun_elev);
- fprintf(stderr, " Method of calculus = %s\n",
+ fprintf(stderr, " 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)));
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++) {
@@ -461,12 +481,12 @@
lsat.band[i].qcalmin, lsat.band[i].qcalmax);
fprintf(stderr, " calibration constants (L): %.3lf to %.3lf\n",
lsat.band[i].lmin, lsat.band[i].lmax);
- fprintf(stderr, " at-%s radiance = %.5lf * DN + %.5lf\n",
+ fprintf(stderr, " at-%s radiance = %.8lf * DN + %.3lf\n",
(method > DOS ? "surface" : "sensor"), lsat.band[i].gain,
lsat.band[i].bias);
if (lsat.band[i].thermal) {
fprintf(stderr,
- " at-sensor temperature = %.3lf / log[(%.3lf / radiance) + 1.0]\n",
+ " at-sensor temperature = %.5lf / log[(%.5lf / radiance) + 1.0]\n",
lsat.band[i].K2, lsat.band[i].K1);
}
else {
@@ -475,12 +495,12 @@
lsat.band[i].esun);
fprintf(stderr, " at-%s reflectance = radiance / %.5lf\n",
(method > DOS ? "surface" : "sensor"),
- lsat.band[i].K2);
+ lsat.band[i].K1);
if (method > DOS) {
fprintf(stderr,
" the darkness DN with a least %d pixels is %d\n",
pixel, dn_dark[i]);
- fprintf(stderr, " the mode of DN is %d\n", dn_mode[i]);
+ fprintf(stderr, " the DN mode is %d\n", dn_mode[i]);
}
}
}
@@ -488,13 +508,14 @@
fflush(stderr);
}
- /*****************************************
- * ------------ CALCULUS -----------------
- *****************************************/
+ /*****************************************
+ * ------------ CALCULUS -----------------
+ *****************************************/
G_message(_("Calculating..."));
for (i = 0; i < lsat.bands; i++) {
- sprintf(band_in, "%s%d", inputname, lsat.band[i].code);
+ sprintf(band_in, "%s%d", inputname,
+ (named->answer ? lsat.band[i].number : lsat.band[i].code));
sprintf(band_out, "%s%d", outputname, lsat.band[i].code);
if ((infd = Rast_open_old(band_in, "")) < 0)
@@ -531,10 +552,9 @@
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);
@@ -609,8 +629,9 @@
/* Append a string to the 'history' structure */
Rast_append_format_history(&history,
" %s of Landsat-%d %s (method %s)",
- 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,
"----------------------------------------------------------------");
@@ -644,7 +665,7 @@
lsat.band[i].esun);
Rast_append_format_history(&history,
" Reflectance = Radiance divided by ..... %.5lf",
- lsat.band[i].K2);
+ lsat.band[i].K1);
if (method > DOS) {
Rast_append_history(&history, " ");
Rast_append_format_history(&history,
@@ -663,8 +684,10 @@
if (lsat.band[i].thermal)
Rast_write_units(band_out, "Kelvin");
- /* else units = ...? */
- /* set raster timestamp from acq date? (see r.timestamp module) */
+ else if (frad->answer)
+ Rast_write_units(band_out, "W/(m^2 sr µm)");
+ else
+ Rast_write_units(band_out, "unitless");
}
exit(EXIT_SUCCESS);
More information about the grass-commit
mailing list