[GRASS-SVN] r30464 - grass-addons/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Mar 4 09:50:46 EST 2008
Author: ejtizado
Date: 2008-03-04 09:50:46 -0500 (Tue, 04 Mar 2008)
New Revision: 30464
Modified:
grass-addons/i.landsat.toar/description.html
grass-addons/i.landsat.toar/landsat.c
grass-addons/i.landsat.toar/landsat.h
grass-addons/i.landsat.toar/landsat_set.c
grass-addons/i.landsat.toar/main.c
Log:
All DOS methods implemented
Modified: grass-addons/i.landsat.toar/description.html
===================================================================
--- grass-addons/i.landsat.toar/description.html 2008-03-04 14:29:12 UTC (rev 30463)
+++ grass-addons/i.landsat.toar/description.html 2008-03-04 14:50:46 UTC (rev 30464)
@@ -1,6 +1,6 @@
<H2>DESCRIPTION</H2>
-<EM>i.landsat.toar</EM> is to transform calibrated digital number of Landsat products to top-of-atmosphere radiance or top-of-atmosphere reflectance and temperature (band 6 of the sensors TM and ETM+). Optionally, to calculate a simplefied at-surface radiance or reflectance.
+<EM>i.landsat.toar</EM> is to transform calibrated digital number of Landsat products to top-of-atmosphere radiance or top-of-atmosphere reflectance and temperature (band 6 of the sensors TM and ETM+). Optionally, to calculate the at-surface radiance or reflectance with atmosferic correction (DOS method).
<p>Usually, to do this the production date, the acquisition date, and the solar elevation is needed. Moreover, also is needed for Landsat-7 ETM+ the gain (high or low) of the nine bands.</p>
@@ -38,21 +38,23 @@
<H2>Corrected at-sensor values (method=corrected)</H2>
-<p>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 also it can be negative. To avoid these possible negative values and set the minimum possible values at-sentor to zero, this method correct the radiance to output a corrected at-sensor values with the equations:
+<p>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 also it can be negative. To avoid these possible negative values and set the minimum possible values at-sentor to zero, this method correct the radiance to output a corrected at-sensor values with the equations (not for thermal bands):
<ul>
- <li> radiance = (uncorrected radiance - Lmin) </li>
+ <li> radiance = (uncorrected_radiance - Lmin) </li>
<li> reflectance = radiance / sun_radiance </li>
</ul>
</p>
+<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 r.mapcalc.</p>
-<H2>Simplified at-surface values (method=simplified)</H2>
-<p>Atmospheric correction and reflectance calibration remove the path radiance, i.e. the stray light from the atmosphere, and the spectral effect of solar illumination. To output these simple <b>at-surface radiance</b> and <b>at-surface reflectance</b>, the equations are:
+<H2>Simplified at-surface values (method=dos[1-4])</H2>
+
+<p>Atmospheric correction and reflectance calibration remove the path radiance, i.e. the stray light from the atmosphere, and the spectral effect of solar illumination. To output these simple <b>at-surface radiance</b> and <b>at-surface reflectance</b>, the equations are (not for thermal bands):
<ul>
- <li> sun_radiance = [Esun · sin(e) · TAUz + Esky] / (PI · d^2) </li>
- <li> rpath = Lmin - percent · sun_radiance </li>
- <li> radiance = (at-sensor radiance - rpath) / TAUv </li>
+ <li> sun_radiance = TAUv · [Esun · sin(e) · TAUz + Esky] / (PI · d^2) </li>
+ <li> radiance_path = radiance_dark - percent · sun_radiance </li>
+ <li> radiance = (at-sensor_radiance - radiance_path) </li>
<li> reflectance = radiance / sun_radiance </li>
</ul>
where,
@@ -60,11 +62,19 @@
<em>Esky</em> is the diffuse sky irradiance,
<em>TAUz</em> is the atmospheric transmittance along the path from the sun to the ground surface, and
<em>TAUv</em> is the atmospheric transmittance along the path from the ground surface to the sensor.
-</p>
+<em>radiance_dark</em> is the at-sensor radiance calculated from the darkest object, i.e. DN with a least 'dark_parameter' (usually 1000) pixels for the entire image.</p>
-<p>For this calculus, the values are TAUv = 1.0, Esky = 0. and TAUz = 1.0. When the method is simplified TAUz = sin(e) for all bands with maximum wave length less than 1., i.e. bands 4-6 MSS, 1-4 TM, and 1-4 ETM+.</p>
+The values are,
+ <ul>
+ <li>DOS1: TAUv = 1.0, TAUz = 1.0 and Esky = 0.0</li>
+ <li>DOS2: TAUv = 1.0, Esky = 0.0, and TAUz = sin(e) for all bands with maximum wave length less than 1. (i.e. bands 4-6 MSS, 1-4 TM, and 1-4 ETM+) other bands TAUz = 1.0</li>
+ <li>DOS3: TAUv = exp[-t/cos(sat_zenith)], TAUz = exp[-t/sin(e)], Esky = rayleigh</li>
+ <li>DOS4: TAUv = exp[-t/cos(sat_zenith)], TAUz = exp[-t/sin(e)], Esky = PI · radiance_dark </li>
+ </ul>
+<p><b>Attention</b>: Output radiance remain untouched (i.e. no set to 0. when it is negative) then they are possible negative values. However, output reflectance is set to 0. when is obtained a negative value.</p>
+
<H2>NOTES</H2>
<p>In verbose mode (flag -v), the program write basic data of satellite and the parameters used in
@@ -101,12 +111,13 @@
<H2>REFERENCES</H2>
<ol>
- <li>G.h Chander and B. Markham: IEEE Transactions On Geoscience And Remote Sensing, vol. 41, no. 11, November 2003.</li>
- <li>Chavez, P. S., jr. 1996. Image-based atmospheric corrections - Revisited and Improved. Photogrammetric Engineering and Remote Sensing 62 (9): 1025-1036.</li>
+ <li>Chander G.H. and B. Markham: IEEE Transactions On Geoscience And Remote Sensing, vol. 41, no. 11, November 2003.</li>
+ <li>Chavez P.S., jr. 1996. Image-based atmospheric corrections - Revisited and Improved. Photogrammetric Engineering and Remote Sensing 62 (9): 1025-1036.</li>
<li>Huang et al: At-Satellite Reflectance: A First Order Normalization Of Landsat 7 ETM+ Images. 2002.</li>
<li>R. Irish: <a href="http://ltpwww.gsfc.nasa.gov/IAS/handbook/handbook_toc.html">Landsat 7. Science Data Users Handbook. February 17, 2007.</a></li>
- <li>B. L. Markham and J.L. Barker: Landsat MSS and TM Post-Calibration Dynamic Ranges, Exoatmospheric Reflectances and At-Satellite Temperatures. EOSAT Landsat Technical Notes, No. 1, 1986</li>
- <li>M.S. Moran, R.D. Jackson, P.N. Slater and P.M. Teillet: Remote Sensing of Environment, vol. 41. 1992.</li>
+ <li>Markham B.L. and J.L. Barker: Landsat MSS and TM Post-Calibration Dynamic Ranges, Exoatmospheric Reflectances and At-Satellite Temperatures. EOSAT Landsat Technical Notes, No. 1, 1986</li>
+ <li>Moran M.S., R.D. Jackson, P.N. Slater and P.M. Teillet: Remote Sensing of Environment, vol. 41. 1992.</li>
+ <li> Song et al : Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects?. Remote Sensing of Environment, vol. 75. 2001. </li>
</ol>
@@ -125,4 +136,4 @@
Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
<p>
-<i>Last changed: $Date: 2007/07/07 00:00:00 $</i>
+<i>Last changed: $Date: 2008/02/10 00:00:00 $</i>
Modified: grass-addons/i.landsat.toar/landsat.c
===================================================================
--- grass-addons/i.landsat.toar/landsat.c 2008-03-04 14:29:12 UTC (rev 30463)
+++ grass-addons/i.landsat.toar/landsat.c 2008-03-04 14:50:46 UTC (rev 30464)
@@ -11,9 +11,9 @@
/****************************************************************************
* PURPOSE: Calibrated Digital Number to at-satellite Radiance
*****************************************************************************/
-double lsat_qcal2rad(int qcal, band_data * band)
+double lsat_qcal2rad(double qcal, band_data * band)
{
- return (double)(((double)qcal) * band->gain + band->bias);
+ return (double)(qcal * band->gain + band->bias);
}
/****************************************************************************
@@ -41,56 +41,124 @@
* lsat : satellite data
* i : band number
* method : level of atmospheric correction
- * aot : atmospheric optical thickness (usually 0.15)
* percent : percent of solar irradiance in path radiance
- *****************************************************************************/
+ * dos : digital number of dark object for DOS
+ *****************************************************************************/
-void lsat_bandctes(lsat_data * lsat, int i, char method, double aot, double percent)
+#define abs(x) (((x)>0)?(x):(-x))
+
+void lsat_bandctes(lsat_data * lsat, int i, char method,
+ double percent, int dos, double sat_zenith, double rayleigh)
{
- double a, b, rad_sun;
- double TAUv, TAUz, Edown;
+ double pi_d2, sin_e, cos_v, rad_sun;
+ /* TAUv = at. transmittance surface-sensor */
+ /* TAUz = at. transmittance sun-surface */
+ /* Edown = diffuse sky spectral irradiance */
+ double TAUv, TAUz, Edown;
- a = (double)(PI * lsat->dist_es * lsat->dist_es);
- b = (double)(sin(D2R * lsat->sun_elev));
+ pi_d2 = (double)(PI * lsat->dist_es * lsat->dist_es);
+ sin_e = (double)(sin(D2R * lsat->sun_elev));
+ cos_v = (double)(cos(D2R * sat_zenith));
- /** Radiance of the Sun */
- if (method == SIMPLIFIED) {
- TAUv = 1.; /* TAUv = at. transmittance surface-sensor */
- /* TAUz = at. transmittance sun-surface */
- TAUz = (lsat->band[i].wavemax < 1.) ? b : 1.;
- Edown = 0.; /* Edown = diffuse sky spectral irradiance */
- }
- else { /* UNCORRECTED or CORRECTED */
- TAUv = 1.;
- TAUz = 1.;
- Edown = 0.;
- }
- rad_sun = ((lsat->band[i].esun * b * TAUz + Edown) / a);
+ /** 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:
+ {
+ TAUv = 1.;
+ TAUz = (lsat->band[i].wavemax < 1.) ? sin_e : 1.;
+ Edown = 0.;
+ break;
+ }
+ case DOS2b:
+ {
+ TAUv = (lsat->band[i].wavemax < 1.) ? cos_v : 1.;
+ TAUz = (lsat->band[i].wavemax < 1.) ? sin_e : 1.;
+ Edown = 0.;
+ break;
+ }
+ case DOS3:
+ {
+ double t;
+ t = 2. / (lsat->band[i].wavemax + lsat->band[i].wavemin);
+ t = 0.008569 *t*t*t*t*(1+0.0113*t*t+0.000013*t*t*t*t);
+ TAUv = exp( -t / cos_v );
+ TAUz = exp( -t / sin_e );
+ Edown = rayleigh;
+ break;
+ }
+ 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].lmin;
+ double tv, Tv = 1.;
+ double Tz = 1.;
+ double Lp = 0.;
+ do
+ {
+ TAUz = Tz;
+ TAUv = Tv;
+ Lp = Ro - percent * TAUv * (lsat->band[i].esun * sin_e * TAUz + PI * Lp) / pi_d2;
+ Tz = 1 - (4 * pi_d2 * Lp) / (lsat->band[i].esun * sin_e);
+ 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.);
+ Edown = (Lp < 0. ? 0. : PI * Lp);
+ break;
+ }
+ default: /* DOS1 and Without atmosferic-correction */
+ TAUv = 1.;
+ TAUz = 1.;
+ Edown = 0.;
+ break;
+ }
+ rad_sun = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
+ G_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz, Edown);
- /** Radiance coefficients */
- lsat->band[i].gain = ((lsat->band[i].lmax - lsat->band[i].lmin) /
- (lsat->band[i].qcalmax - lsat->band[i].qcalmin));
- lsat->band[i].bias = (lsat->band[i].lmin -
- lsat->band[i].gain * lsat->band[i].qcalmin);
+ lsat->band[i].K1 = 0.;
+ lsat->band[i].K2 = rad_sun;
+ }
- if (method == CORRECTED)
- {
- lsat->band[i].bias -= lsat->band[i].lmin;
- }
- else if (method == SIMPLIFIED)
- {
- lsat->band[i].gain /= TAUv;
- lsat->band[i].bias -= (lsat->band[i].lmin - percent * rad_sun);
- lsat->band[i].bias /= TAUv;
- }
+ /** Digital number to radiance coefficients.
+ Whitout atmosferic calibration for thermal bands.
+ */
+ lsat->band[i].gain = ((lsat->band[i].lmax - lsat->band[i].lmin) /
+ (lsat->band[i].qcalmax - lsat->band[i].qcalmin));
- /** Reflectance coefficients, only NO thermal bands.
- Same variables are utilized to thermal constants
- */
- if (lsat->band[i].thermal == 0)
- {
- lsat->band[i].K1 = 0.;
- lsat->band[i].K2 = rad_sun;
- }
+ if (method == UNCORRECTED || lsat->band[i].thermal)
+ {
+ /* L = G * (DN - Qmin) + Lmin
+ -> bias = Lmin - G * Qmin */
+ lsat->band[i].bias = (lsat->band[i].lmin - lsat->band[i].gain * lsat->band[i].qcalmin);
+ }
+ 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;
+ }
+ }
}
+
+
+
Modified: grass-addons/i.landsat.toar/landsat.h
===================================================================
--- grass-addons/i.landsat.toar/landsat.h 2008-03-04 14:29:12 UTC (rev 30463)
+++ grass-addons/i.landsat.toar/landsat.h 2008-03-04 14:50:46 UTC (rev 30464)
@@ -1,11 +1,14 @@
#ifndef _LANDSAT_H
#define _LANDSAT_H
-#define MAX_BANDS 9
-
#define UNCORRECTED 0
#define CORRECTED 1
-#define SIMPLIFIED 2
+#define DOS 10
+#define DOS1 12
+#define DOS2 14
+#define DOS2b 15
+#define DOS3 16
+#define DOS4 18
/*****************************************************
@@ -15,37 +18,38 @@
* Esun in W / (m² · µm) -> Irradiance
****************************************************/
+#define MAX_BANDS 9
+
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 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 */
+ 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,
+ char thermal; /* Flag to thermal band */
+ double gain, bias; /* Gain and Bias of sensor */
+ double K1, K2; /* Thermal calibration constants,
or Rad2Ref constants */
} band_data;
typedef struct
{
- unsigned char number; /* Landsat number */
+ unsigned char number; /* Landsat number */
- char creation[11]; /* Image production date */
- char date[11]; /* Image acquisition date */
- double dist_es; /* Distance Earth-Sun */
- double sun_elev; /* Solar elevation */
+ char creation[11]; /* Image production date */
+ char date[11]; /* Image acquisition date */
+ double dist_es; /* Distance Earth-Sun */
+ double sun_elev; /* Solar elevation */
- char sensor[5]; /* Type of sensor: MSS, TM, ETM+ */
- int bands; /* Total number of bands */
- band_data band[MAX_BANDS]; /* Data for each band */
-
+ char sensor[5]; /* Type of sensor: MSS, TM, ETM+ */
+ int bands; /* Total number of bands */
+ band_data band[MAX_BANDS]; /* Data for each band */
} lsat_data;
@@ -53,10 +57,10 @@
* Landsat Equations Prototypes
*****************************************************************************/
-double lsat_qcal2rad(int, band_data *);
-double lsat_rad2ref(double, band_data *);
+double lsat_qcal2rad(double, band_data *);
+double lsat_rad2ref (double, band_data *);
double lsat_rad2temp(double, band_data *);
-void lsat_bandctes(lsat_data *, int, char, double, double);
+void lsat_bandctes(lsat_data *, int, char, double, int, double, double);
#endif
Modified: grass-addons/i.landsat.toar/landsat_set.c
===================================================================
--- grass-addons/i.landsat.toar/landsat_set.c 2008-03-04 14:29:12 UTC (rev 30463)
+++ grass-addons/i.landsat.toar/landsat_set.c 2008-03-04 14:50:46 UTC (rev 30464)
@@ -11,73 +11,79 @@
void sensor_MSS(lsat_data * lsat)
{
- int i;
+ int i;
- int band[] = { 1, 2, 3, 4 };
- int code[] = { 4, 5, 6, 7 };
- double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
- double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
+ /* green, red, near infrared, near infrared */
+ int band[] = { 1, 2, 3, 4 };
+ int code[] = { 4, 5, 6, 7 };
+ double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
+ double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
+ /* 79-82, 79-82, 79-82, 79-82 */
- strcpy(lsat->sensor, "MSS");
+ strcpy(lsat->sensor, "MSS");
- lsat->bands = 4;
- 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 = 0.;
- lsat->band[i].thermal = 0;
- }
- return;
+ lsat->bands = 4;
+ 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 = 0.;
+ lsat->band[i].thermal = 0;
+ }
+ return;
}
void sensor_TM(lsat_data * lsat)
{
- int i;
+ int i;
- int band[] = { 1, 2, 3, 4, 5, 6, 7 };
- double wmax[] = { 0.52, 0.60, 0.69, 0.90, 1.75, 12.50, 2.35 };
- double wmin[] = { 0.45, 0.52, 0.63, 0.76, 1.55, 10.40, 2.08 };
+ /* blue, green red, near infrared, shortwave IR, thermal IR, shortwave IR */
+ int band[] = { 1, 2, 3, 4, 5, 6, 7 };
+ double wmax[] = { 0.52, 0.60, 0.69, 0.90, 1.75, 12.50, 2.35 };
+ double wmin[] = { 0.45, 0.52, 0.63, 0.76, 1.55, 10.40, 2.08 };
+ /* 30, 30, 30, 30, 30, 120, 30 */
- strcpy(lsat->sensor, "TM");
+ strcpy(lsat->sensor, "TM");
- lsat->bands = 7;
- for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(band + i);
- lsat->band[i].wavemax = *(wmax + i);
- lsat->band[i].wavemin = *(wmin + i);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 0.;
- lsat->band[i].thermal = (lsat->band[i].number == 6);
- }
- return;
+ lsat->bands = 7;
+ for (i = 0; i < lsat->bands; i++) {
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = *(band + i);
+ lsat->band[i].wavemax = *(wmax + i);
+ lsat->band[i].wavemin = *(wmin + i);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 0.;
+ lsat->band[i].thermal = (lsat->band[i].number == 6 ? 1 : 0);
+ }
+ return;
}
void sensor_ETM(lsat_data * lsat)
{
- int i;
+ int i;
- int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
- int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
- 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, 2.09, 0.52 };
+ /* 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 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, 2.09, 0.52 };
+ /* 30, 30, 30, 30, 30, 60, 30, 15 */
- strcpy(lsat->sensor, "ETM+");
+ strcpy(lsat->sensor, "ETM+");
- lsat->bands = 9;
- 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 == 6);
- }
- return;
+ lsat->bands = 9;
+ 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 == 6 ? 1 : 0);
+ }
+ return;
}
@@ -337,7 +343,7 @@
/** Gyanesh Chander and Brian Markham.
IEEE Transactions On Geoscience And Remote Sensing, Vol. 41, No. 11, November 2003 */
/* Spectral radiances at detector */
- double Lmax[][7] = { { 152.10, 296.81, 204.30, 206.20, 27.19, 15.303, 14.38 }, /* on or before May 4, 2003 */
+ double Lmax[][7] = { { 152.10, 296.81, 204.30, 206.20, 27.19, 15.303, 14.38 }, /* before May 4, 2003 */
{ 193.00, 365.00, 264.00, 221.00, 30.20, 15.303, 16.50 }, /* after May 4, 2003 */
{ 169.00, 333.00, 264.00, 221.00, 30.20, 15.303, 16.50 } }; /* after April 2, 2007 */
double Lmin[][7] = { { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 },
@@ -353,7 +359,7 @@
else i = 2;
lmax = Lmax[i];
lmin = Lmin[i];
- if ( i == 2 ) {
+ if ( i == 2 ) { /* in Chander, Markham and Barsi 2007 */
julian = julian_char(lsat->date); /* Yes, here acquisition date */
if (julian >= julian_char("1992-01-01")) {
lmax[0] = 193.0;
@@ -406,32 +412,32 @@
double esun[] = { 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
/* 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;
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("2000-07-01")) k = 0;
+ else k = 1;
- lsat->number = 7;
- sensor_ETM( lsat );
+ lsat->number = 7;
+ sensor_ETM( lsat );
- lsat->dist_es = earth_sun(lsat->date);
+ 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);
- if (gain[i] == 'H' || gain[i] == 'h') {
- lmax = LmaxH[k];
- lmin = LminH[k];
- }
- else {
- lmax = LmaxL[k];
- lmin = LminL[k];
- }
- lsat->band[i].lmax = *(lmax + j);
- lsat->band[i].lmin = *(lmin + j);
- if (lsat->band[i].thermal ) {
- lsat->band[i].K1 = 666.09;
- lsat->band[i].K2 = 1282.71;
- }
+ 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];
+ }
+ else {
+ lmax = LmaxL[k];
+ lmin = LminL[k];
+ }
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
+ if (lsat->band[i].thermal) {
+ lsat->band[i].K1 = 666.09;
+ lsat->band[i].K2 = 1282.71;
+ }
}
return;
}
Modified: grass-addons/i.landsat.toar/main.c
===================================================================
--- grass-addons/i.landsat.toar/main.c 2008-03-04 14:29:12 UTC (rev 30463)
+++ grass-addons/i.landsat.toar/main.c 2008-03-04 14:50:46 UTC (rev 30464)
@@ -23,37 +23,36 @@
#include "local_proto.h"
-
int main(int argc, char *argv[])
{
- struct History history;
- struct GModule *module;
+ struct History history;
+ struct GModule *module;
- struct Cell_head cellhd, window;
- char *mapset;
+ struct Cell_head cellhd, window;
+ char *mapset;
- void *inrast, *outrast;
- int infd, outfd;
- void *ptr;
- int nrows, ncols, row, col;
+ void *inrast, *outrast;
+ int infd, outfd;
+ void *ptr;
+ int nrows, ncols, row, col;
- RASTER_MAP_TYPE in_data_type;
+ RASTER_MAP_TYPE in_data_type;
- int verbose = 1;
- struct Option *input, *metfn,
- *adate, *elev, *bgain, *pdate, *metho, * perc;
- char *name;
- char *met;
- char *atcor;
- struct Flag *sat1, *sat2, *sat3, *sat4, *sat5, *sat7, *msss,
- *verbo, *frad;
+ int verbose = 1;
+ struct Option *input, *metfn,
+ *adate, *pdate, *elev, *bgain, *metho, *perc, *dark,
+ *satz, *atmo;
+ char *name, *met;
+ struct Flag *sat1, *sat2, *sat3, *sat4, *sat5, *sat7, *msss,
+ *verbo, *frad;
- char band_in[127], band_out[127];
- int i, method, qcal;
- double rad, ref, percent;
- lsat_data lsat;
- char command[300];
+ lsat_data lsat;
+ char band_in[127], band_out[127];
+ int i, j, q, method, pixel, dn_dark[MAX_BANDS], dn_mode[MAX_BANDS];
+ double qcal, rad, ref, percent, ref_mode, sat_zenith, rayleigh;
+ unsigned long hist[256], h_max;
+
/* initialize GIS environment */
G_gisinit(argv[0]);
@@ -64,373 +63,515 @@
/* It defines the different parameters */
input = G_define_option();
- input->key = _("band_prefix");
- input->type = TYPE_STRING;
- input->required = YES;
- input->gisprompt = _("input,cell,raster");
+ input->key = _("band_prefix");
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->gisprompt = _("input,cell,raster");
input->description = _("Base name of the landsat band rasters (.#)");
metfn = G_define_option();
- metfn->key = _("metfile");
- metfn->type = TYPE_STRING;
- metfn->required = NO;
- metfn->gisprompt = _(".met file");
+ metfn->key = _("metfile");
+ metfn->type = TYPE_STRING;
+ metfn->required = NO;
+ metfn->gisprompt = _(".met file");
metfn->description = _("Landsat ETM+ or TM5 header file (.met)");
adate = G_define_option();
- adate->key = _("date");
- adate->type = TYPE_STRING;
- adate->required = NO;
- adate->gisprompt = _("image acquisition date");
+ adate->key = _("date");
+ adate->type = TYPE_STRING;
+ adate->required = NO;
+ adate->gisprompt = _("image acquisition date");
adate->description = _("Image acquisition date (yyyy-mm-dd)");
elev = G_define_option();
- elev->key = _("solar");
- elev->type = TYPE_DOUBLE;
- elev->required = NO;
- elev->gisprompt = _("solar elevation");
+ elev->key = _("solar_elevation");
+ elev->type = TYPE_DOUBLE;
+ elev->required = NO;
+ elev->gisprompt = _("solar elevation");
elev->description = _("Solar elevation in degrees");
bgain = G_define_option();
- bgain->key = _("gain");
- bgain->type = TYPE_STRING;
- bgain->required = NO;
- bgain->gisprompt = _("band gain");
+ bgain->key = _("gain");
+ bgain->type = TYPE_STRING;
+ bgain->required = NO;
+ bgain->gisprompt = _("band gain");
bgain->description = _("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
pdate = G_define_option();
- pdate->key = _("product_date");
- pdate->type = TYPE_STRING;
- pdate->required = NO;
- pdate->gisprompt = _("image production date");
+ pdate->key = _("product_date");
+ pdate->type = TYPE_STRING;
+ pdate->required = NO;
+ pdate->gisprompt = _("image production date");
pdate->description = _("Image creation date (yyyy-mm-dd)");
metho = G_define_option();
- metho->key = _("method");
- metho->type = TYPE_STRING;
- metho->required = NO;
- metho->options = "uncorrected,corrected,simplified";
- metho->gisprompt = _("atmosferic correction methods");
- metho->description = _("Atmosferic correction methods");
- metho->answer = "uncorrected";
+ metho->key = _("method");
+ metho->type = TYPE_STRING;
+ metho->required = NO;
+ metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
+ metho->gisprompt = _("atmosferic correction method");
+ metho->description = _("Atmosferic correction method");
+ metho->answer = "uncorrected";
- perc = G_define_option();
- perc->key = _("percent");
- perc->type = TYPE_DOUBLE;
- perc->required = NO;
- perc->gisprompt = _("percent of solar irradiance in path radiance");
- perc->description = _("Percent of solar irradiance in path radiance");
- perc->answer = "0.01";
+ perc = G_define_option();
+ perc->key = _("percent");
+ perc->type = TYPE_DOUBLE;
+ perc->required = NO;
+ perc->gisprompt = _("percent of solar radiance in path radiance");
+ perc->description = _("Percent of solar radiance in path radiance");
+ perc->answer = "0.01";
- /* It defines the different flags */
+ dark = G_define_option();
+ dark->key = _("pixel");
+ dark->type = TYPE_INTEGER;
+ dark->required = NO;
+ dark->gisprompt = _("pixels to dark object");
+ dark->description = _("Minimum pixels to consider digital number as dark object");
+ dark->answer = "1000";
+
+ satz = G_define_option();
+ satz->key = _("sat_zenith");
+ satz->type = TYPE_DOUBLE;
+ satz->required = NO;
+ satz->gisprompt = _("satellite zenith");
+ satz->description = _("Satellite zenith in degrees");
+ satz->answer = "8.2000";
+
+ atmo = G_define_option();
+ atmo->key = _("rayleigh");
+ atmo->type = TYPE_DOUBLE;
+ atmo->required = NO;
+ atmo->gisprompt = _("Rayleigh atmosphere");
+ atmo->description = _("Rayleigh atmosphere");
+ atmo->answer = "0.0";
+
+ /* It defines the different flags */
frad = G_define_flag();
- frad->key = 'r';
+ frad->key = 'r';
frad->description = _("Output at-sensor radiance for all bands");
- frad->answer = 0;
+ frad->answer = 0;
sat1 = G_define_flag();
- sat1->key = '1';
+ sat1->key = '1';
sat1->description = _("Landsat-1 MSS");
- sat1->answer = 0;
+ sat1->answer = 0;
sat2 = G_define_flag();
- sat2->key = '2';
+ sat2->key = '2';
sat2->description = _("Landsat-2 MSS");
- sat2->answer = 0;
+ sat2->answer = 0;
sat3 = G_define_flag();
- sat3->key = '3';
+ sat3->key = '3';
sat3->description = _("Landsat-3 MSS");
- sat3->answer = 0;
+ sat3->answer = 0;
sat4 = G_define_flag();
- sat4->key = '4';
+ sat4->key = '4';
sat4->description = _("Landsat-4 TM");
- sat4->answer = 0;
+ sat4->answer = 0;
sat5 = G_define_flag();
- sat5->key = '5';
+ sat5->key = '5';
sat5->description = _("Landsat-5 TM");
- sat5->answer = 0;
+ sat5->answer = 0;
sat7 = G_define_flag();
- sat7->key = '7';
+ sat7->key = '7';
sat7->description = _("Landsat-7 ETM+");
- sat7->answer = 0;
+ sat7->answer = 0;
msss = G_define_flag();
- msss->key = 's';
- msss->description = _("Modify sensor of Landsat-4/5 to MSS");
- msss->answer = 0;
+ msss->key = 's';
+ msss->description = _("Set sensor of Landsat-4/5 to MSS");
+ msss->answer = 0;
verbo = G_define_flag();
- verbo->key = 'v';
+ verbo->key = 'v';
verbo->description = _("Show parameters applied");
- verbo->answer = 0;
+ verbo->answer = 0;
/* options and afters parser */
if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ exit(EXIT_FAILURE);
/*****************************************
* ---------- START --------------------
* Stores options and flag to variables
*****************************************/
- met = metfn->answer;
- name = input->answer;
- atcor = metho->answer;
+ met = metfn->answer;
+ name = input->answer;
- if (adate->answer != NULL) {
- strncpy(lsat.date, adate->answer, 11);
- lsat.date[10] = '\0';
- }
- else
- lsat.date[0] = '\0';
+ if (adate->answer != NULL) {
+ G_strncpy(lsat.date, adate->answer, 11);
+ lsat.date[10] = '\0';
+ }
+ else
+ lsat.date[0] = '\0';
- if (pdate->answer != NULL) {
- strncpy(lsat.creation, pdate->answer, 11);
- lsat.creation[10] = '\0';
- }
- else
- lsat.creation[0] = '\0';
+ if (pdate->answer != NULL) {
+ G_strncpy(lsat.creation, pdate->answer, 11);
+ lsat.creation[10] = '\0';
+ }
+ else
+ lsat.creation[0] = '\0';
- lsat.sun_elev = elev->answer == NULL ? 0. : atof(elev->answer);
+ lsat.sun_elev = elev->answer == NULL ? 0. : atof(elev->answer);
+ percent = atof(perc->answer);
+ pixel = atoi(dark->answer);
+ sat_zenith = atof(satz->answer);
+ rayleigh = atof(atmo->answer);
- percent = atof(perc->answer);
+ if (met == NULL &&
+ (sat1->answer + sat2->answer + sat3->answer +
+ sat4->answer + sat5->answer + sat7->answer) != 1)
+ G_fatal_error(_("Select one and only one type of satellite"));
- if (met == NULL &&
- (sat1->answer + sat2->answer + sat3->answer +
- sat4->answer + sat5->answer + sat7->answer) != 1)
- G_fatal_error(_("Select one and only one type of satellite"));
-
/* Data from MET file: only Landsat-7 ETM+ and Landsat-5 TM */
- if (met != NULL) {
- if (sat7->answer) met_ETM(met, &lsat);
- else met_TM5(met, &lsat);
- fprintf(stdout, "Landsat-%d %s with data set in met file [%s]\n",
- lsat.number, lsat.sensor, met);
- if (elev->answer != NULL)
- lsat.sun_elev = atof(elev->answer); /* Overwrite sun elevation of met file */
- }
- /* Data from date and solar elevation */
- else if (adate->answer == NULL || elev->answer == NULL) {
- G_fatal_error(_("Lacking date or solar elevation for this satellite"));
- }
+ if (met != NULL) {
+ if (sat7->answer) met_ETM(met, &lsat); else met_TM5(met, &lsat);
+ G_message("Landsat-%d %s with data set in met file [%s]", lsat.number, lsat.sensor, met);
+ if (elev->answer != NULL)
+ lsat.sun_elev = atof(elev->answer); /* Overwrite solar elevation of met file */
+ }
+ /* Data from date and solar elevation */
+ else if (adate->answer == NULL || elev->answer == NULL) {
+ G_fatal_error(_("Lacking date or solar elevation for this satellite"));
+ }
else {
- if (sat7->answer) { /* Need gain */
- if (bgain->answer && strlen(bgain->answer) == 9) {
- set_ETM(&lsat, bgain->answer);
- fprintf(stdout, "Landsat 7 ETM+\n");
- }
- else {
- G_fatal_error(_("For Landsat-7 is necessary band gain with 9 (H/L) data"));
- }
+ if (sat7->answer) { /* Need gain */
+ if (bgain->answer != NULL && strlen(bgain->answer) == 9) {
+ set_ETM(&lsat, bgain->answer);
+ G_message("Landsat 7 ETM+");
+ }
+ else {
+ G_fatal_error(_("For Landsat-7 is necessary band gain with 9 (H/L) data"));
+ }
+ }
+ else { /* Not need gain */
+ if (sat5->answer) {
+ if (msss->answer) set_MSS5(&lsat); else set_TM5 (&lsat);
+ G_message("Landsat-5 %s", lsat.sensor);
+ }
+ else if (sat4->answer) {
+ if (msss->answer) set_MSS4(&lsat); else set_TM4 (&lsat);
+ G_message("Landsat-4 %s", lsat.sensor);
+ }
+ else if (sat3->answer) {
+ set_MSS3(&lsat);
+ G_message("Landsat-3 MSS");
+ }
+ else if (sat2->answer) {
+ set_MSS2(&lsat);
+ G_message("Landsat-2 MSS");
+ }
+ else if (sat1->answer) {
+ set_MSS1(&lsat);
+ G_message("Landsat-1 MSS");
+ }
+ else {
+ G_fatal_error(_("Lacking satellite type"));
+ }
+ }
}
- else { /* Not need gain */
- if (sat5->answer) {
- if (msss->answer) set_MSS5(&lsat);
- else set_TM5 (&lsat);
- fprintf(stdout, "Landsat-5 %s\n", lsat.sensor);
- }
- else if (sat4->answer) {
- if (msss->answer) set_MSS4(&lsat);
- else set_TM4 (&lsat);
- fprintf(stdout, "Landsat-4 %s\n", lsat.sensor);
- }
- else if (sat3->answer) {
- set_MSS3(&lsat);
- fprintf(stdout, "Landsat-3 MSS\n");
- }
- else if (sat2->answer) {
- set_MSS2(&lsat);
- fprintf(stdout, "Landsat-2 MSS\n");
- }
- else if (sat1->answer) {
- set_MSS1(&lsat);
- fprintf(stdout, "Landsat-1 MSS\n");
- }
- else
- {
- G_fatal_error(_("Lacking satellite type"));
- }
- }
- }
- /* Set method and calculate band constants */
- switch(atcor[0])
- {
- case 'c':
- method = CORRECTED;
- break;
- case 's':
- method = SIMPLIFIED;
- break;
- default:
- method = UNCORRECTED;
- break;
- }
- for (i = 0; i < lsat.bands; i++)
- {
- lsat_bandctes(&lsat, i, method, 0.0, percent);
- }
+ /*****************************************
+ * ------------ PREPARATION --------------
+ *****************************************/
+ if (G_strcasecmp(metho->answer, "corrected") == 0) method = CORRECTED;
+ else if (G_strcasecmp(metho->answer, "dos1") == 0) method = DOS1;
+ else if (G_strcasecmp(metho->answer, "dos2") == 0) method = DOS2;
+ else if (G_strcasecmp(metho->answer, "dos2b") == 0) method = DOS2b;
+ else if (G_strcasecmp(metho->answer, "dos3") == 0) method = DOS3;
+ else if (G_strcasecmp(metho->answer, "dos4") == 0) method = DOS4;
+ else 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;
- /*****************************************
- * ---------- VERBOSE ------------------
- *****************************************/
- if (verbo->answer) {
- fprintf(stdout, " ACQUISITION DATE %s [production date %s]\n", lsat.date, lsat.creation);
- fprintf(stdout, " earth-sun distance = %.8lf\n", lsat.dist_es);
- fprintf(stdout, " solar elevation angle = %.8lf\n", lsat.sun_elev);
- fprintf(stdout, " Method for at sensor values = %s\n",
- (atcor[0]=='c' ? "corrected" : (atcor[0]=='s' ? "simplified" : "uncorrected")));
- if (atcor[0] == 's')
- {
- fprintf(stdout, " percent of solar irradiance in path radiance = %.4lf\n", percent);
- }
- for (i = 0; i < lsat.bands; i++) {
- fprintf(stdout, "-------------------\n");
- fprintf(stdout, " BAND %d %s (code %d)\n",
- lsat.band[i].number, (lsat.band[i].thermal ? "thermal " : ""),
- lsat.band[i].code);
- fprintf(stdout, " calibrated digital number (DN): %.1lf to %.1lf\n",
- lsat.band[i].qcalmin, lsat.band[i].qcalmax);
- fprintf(stdout, " calibration constants (L): %.3lf to %.3lf\n",
- lsat.band[i].lmin, lsat.band[i].lmax);
- fprintf(stdout, " at-%s radiance = %.5lf · DN + %.5lf\n",
- (atcor[0] == 's' ? "surface" : "sensor"),
- lsat.band[i].gain, lsat.band[i].bias);
- if (lsat.band[i].thermal)
- {
- fprintf(stdout, " at-%s temperature = %.3lf / log[(%.3lf / radiance) + 1.0]\n",
- (atcor[0] == 's' ? "surface" : "sensor"),
- lsat.band[i].K2, lsat.band[i].K1);
- }
- else {
- fprintf(stdout, " mean solar exoatmospheric irradiance (ESUN): %.3lf\n",
- lsat.band[i].esun);
- fprintf(stdout, " at-%s reflectance = radiance / %.5lf\n",
- (atcor[0] == 's' ? "surface" : "sensor"),
- lsat.band[i].K2);
- }
+ for (i = 0; i < lsat.bands; i++)
+ {
+ dn_mode[i] = 0;
+ dn_dark[i] = (int)lsat.band[i].qcalmin;
+ /* Calculate dark pixel */
+ if (method > DOS && !lsat.band[i].thermal)
+ {
+ for (j = 0; j < 256; j++) hist[j] = 0L;
+
+ snprintf(band_in, 127, "%s.%d", name, lsat.band[i].code);
+ mapset = G_find_cell2(band_in, "");
+ if (mapset == NULL) {
+ G_warning(_("Raster file [%s] not found"), band_in);
+ continue;
+ }
+ if ((infd = G_open_cell_old(band_in, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), band_in);
+ if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), band_in);
+ if (G_set_window(&cellhd) < 0)
+ G_fatal_error(_("Unable to set region"));
+
+ in_data_type = G_raster_map_type(band_in, mapset);
+ inrast = G_allocate_raster_buf(in_data_type);
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ G_message("Calculating dark pixel of [%s] ... ", band_in);
+ for (row = 0; row < nrows; row++)
+ {
+ G_get_raster_row(infd, inrast, row, in_data_type);
+ for (col = 0; col < ncols; col++)
+ {
+ switch(in_data_type)
+ {
+ case CELL_TYPE:
+ ptr = (void *)((CELL *)inrast + col);
+ q = (int)*((CELL *)ptr);
+ break;
+ case FCELL_TYPE:
+ ptr = (void *)((FCELL *)inrast + col);
+ q = (int)*((FCELL *)ptr);
+ break;
+ case DCELL_TYPE:
+ ptr = (void *)((DCELL *)inrast + col);
+ q = (int)*((DCELL *)ptr);
+ break;
+ }
+ if (!G_is_null_value(ptr, in_data_type) &&
+ q >= lsat.band[i].qcalmin && q < 256)
+ hist[q]++;
+ }
+ }
+ /* DN of dark object */
+ for (j = lsat.band[i].qcalmin; j < 256; j++)
+ {
+ if (hist[j] >= 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 */
+ {
+// printf("%d-%ld\n", j, hist[j]);
+ if (hist[j] > h_max)
+ {
+ h_max = hist[j];
+ dn_mode[i] = j;
+ }
+ }
+ G_message("... DN = %.2d [%d] : mode %.2d [%d]",
+ dn_dark[i], hist[dn_dark[i]],
+ dn_mode[i], hist[dn_mode[i]],
+ hist[255] > hist[dn_mode[i]] ? ", excluding DN > 241" : "");
+
+ G_free(inrast);
+ G_close_cell(infd);
+ }
+ /* Calculate transformation constants */
+ lsat_bandctes(&lsat, i, method, percent, dn_dark[i], sat_zenith, rayleigh);
}
- fprintf(stdout, "-------------------\n");
- fflush(stdout);
- }
- /*****************************************
- * ---------- CALCULUS -----------------
- *****************************************/
- G_get_window(&window);
+ /*****************************************
+ * ------------ VERBOSE ------------------
+ *****************************************/
+ if (verbo->answer)
+ {
+ fprintf(stdout, " ACQUISITION DATE %s [production date %s]\n", lsat.date, lsat.creation);
+ fprintf(stdout, " earth-sun distance = %.8lf\n", lsat.dist_es);
+ fprintf(stdout, " solar elevation angle = %.8lf\n", lsat.sun_elev);
+ fprintf(stdout, " Method of calculus = %s\n",
+ (method==CORRECTED ? "CORRECTED"
+ : (method == UNCORRECTED ? "UNCORRECTED" : metho->answer)));
+ if (method > DOS) {
+ fprintf(stdout, " percent of solar irradiance in path radiance = %.4lf\n", percent);
+ }
+ for (i = 0; i < lsat.bands; i++)
+ {
+ fprintf(stdout, "-------------------\n");
+ fprintf(stdout, " BAND %d %s (code %d)\n",
+ lsat.band[i].number, (lsat.band[i].thermal ? "thermal " : ""),
+ lsat.band[i].code);
+ fprintf(stdout, " calibrated digital number (DN): %.1lf to %.1lf\n",
+ lsat.band[i].qcalmin, lsat.band[i].qcalmax);
+ fprintf(stdout, " calibration constants (L): %.3lf to %.3lf\n",
+ lsat.band[i].lmin, lsat.band[i].lmax);
+ fprintf(stdout, " at-%s radiance = %.5lf · DN + %.5lf\n",
+ (method > DOS ? "surface" : "sensor"),
+ lsat.band[i].gain, lsat.band[i].bias);
+ if (lsat.band[i].thermal) {
+ fprintf(stdout, " at-sensor temperature = %.3lf / log[(%.3lf / radiance) + 1.0]\n",
+ lsat.band[i].K2, lsat.band[i].K1);
+ }
+ else {
+ fprintf(stdout, " mean solar exoatmospheric irradiance (ESUN): %.3lf\n",
+ lsat.band[i].esun);
+ fprintf(stdout, " at-%s reflectance = radiance / %.5lf\n",
+ (method > DOS ? "surface" : "sensor"),
+ lsat.band[i].K2);
+ if (method > DOS) {
+ fprintf(stdout, " the darkness DN with a least %d pixels is %d\n", pixel, dn_dark[i] );
+ fprintf(stdout, " the mode of DN is %d\n", dn_mode[i] );
+ }
+ }
+ }
+ fprintf(stdout, "-------------------\n");
+ fflush(stdout);
+ }
- for (i = 0; i < lsat.bands; i++) {
- snprintf(band_in, 127, "%s.%d", name, lsat.band[i].code);
- snprintf(band_out, 127, "%s.toar.%d", name, lsat.band[i].code);
+ /*****************************************
+ * ------------ CALCULUS -----------------
+ *****************************************/
- mapset = G_find_cell2(band_in, "");
- if (mapset == NULL) {
- G_warning(_("raster file [%s] not found"), band_in);
- continue; }
+ for (i = 0; i < lsat.bands; i++)
+ {
+ snprintf(band_in, 127, "%s.%d", name, lsat.band[i].code);
+ snprintf(band_out, 127, "%s.toar.%d", name, lsat.band[i].code);
- if (G_legal_filename(band_out) < 0)
- G_fatal_error(_("[%s] is an illegal name"), band_out);
+ mapset = G_find_cell2(band_in, "");
+ if (mapset == NULL) {
+ G_warning(_("raster file [%s] not found"), band_in);
+ continue; }
- in_data_type = G_raster_map_type(band_in, mapset);
- if ((infd = G_open_cell_old(band_in, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), band_in);
+ in_data_type = G_raster_map_type(band_in, mapset);
+ if ((infd = G_open_cell_old(band_in, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), band_in);
- if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), band_in);
+ if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), band_in);
- /* set same size as original band raster */
- if (G_set_window(&cellhd) < 0)
- G_fatal_error(_("Unable to set region"));
+ /* set same size as original band raster */
+ if (G_set_window(&cellhd) < 0)
+ G_fatal_error(_("Unable to set region"));
- /* Allocate input and output buffer */
- inrast = G_allocate_raster_buf(in_data_type);
- outrast = G_allocate_raster_buf(DCELL_TYPE);
+ /* controlling, if we can write the raster */
+ if (G_legal_filename(band_out) < 0)
+ G_fatal_error(_("[%s] is an illegal name"), band_out);
- /* controlling, if we can write the raster */
- if ((outfd = G_open_raster_new(band_out, DCELL_TYPE)) < 0)
- G_fatal_error(_("Could not open <%s>"), band_out);
+ if ((outfd = G_open_raster_new(band_out, DCELL_TYPE)) < 0)
+ G_fatal_error(_("Could not open <%s>"), band_out);
- /* ================================================================= */
- G_message("%s of %s to %s",
- (frad->answer ? "Radiance" : "Reflectance"),
- band_in, band_out);
+ /* 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();
- for (row = 0; row < nrows; row++)
- {
- if (verbose)
- G_percent(row, nrows, 2);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ /* ================================================================= */
+ G_message("%s of %s to %s",
+ (frad->answer ? "Radiance"
+ : (lsat.band[i].thermal) ? "Temperature" : "Reflectance"),
+ band_in, band_out);
- if (G_get_raster_row(infd, inrast, row, in_data_type) < 0)
- G_fatal_error(_("Could not read from <%s>"), band_in);
+ for (row = 0; row < nrows; row++)
+ {
+ if (verbose)
+ G_percent(row, nrows, 2);
- for (col = 0; col < ncols; col++)
- {
- switch(in_data_type)
- {
- case CELL_TYPE:
- ptr = (void *)((CELL *)inrast + col);
- qcal = (int)((CELL *) inrast)[col];
- break;
- case FCELL_TYPE:
- ptr = (void *)((FCELL *)inrast + col);
- qcal = (int)((FCELL *) inrast)[col];
- break;
- case DCELL_TYPE:
- ptr = (void *)((DCELL *)inrast + col);
- qcal = (int)((DCELL *) inrast)[col];
- break;
- }
- if (G_is_null_value(ptr, in_data_type) ||
- qcal < lsat.band[i].qcalmin)
- {
- G_set_d_null_value((DCELL *)outrast + col, 1);
- }
- else
- {
- rad = lsat_qcal2rad(qcal, &lsat.band[i]);
- if (frad->answer)
- {
- ref = rad;
- }
- else
- {
- if (lsat.band[i].thermal)
- {
- ref = lsat_rad2temp(rad, &lsat.band[i]);
- }
- else
- {
- ref = lsat_rad2ref(rad, &lsat.band[i]);
- }
- }
- ((DCELL *) outrast)[col] = ref;
- }
- }
+ G_get_raster_row(infd, inrast, row, in_data_type);
+ for (col = 0; col < ncols; col++)
+ {
+ switch(in_data_type)
+ {
+ case CELL_TYPE:
+ ptr = (void *)((CELL *)inrast + col);
+ qcal = (double)((CELL *) inrast)[col];
+ break;
+ case FCELL_TYPE:
+ ptr = (void *)((FCELL *)inrast + col);
+ qcal = (double)((FCELL *) inrast)[col];
+ break;
+ case DCELL_TYPE:
+ ptr = (void *)((DCELL *)inrast + col);
+ qcal = (double)((DCELL *) inrast)[col];
+ break;
+ }
+ if (G_is_null_value(ptr, in_data_type) ||
+ qcal < lsat.band[i].qcalmin)
+ {
+ G_set_d_null_value((DCELL *)outrast + col, 1);
+ }
+ else
+ {
+ rad = lsat_qcal2rad(qcal, &lsat.band[i]);
+ if (frad->answer)
+ {
+ ref = rad;
+ }
+ else
+ {
+ if (lsat.band[i].thermal)
+ {
+ ref = lsat_rad2temp(rad, &lsat.band[i]);
+ }
+ else
+ {
+ ref = lsat_rad2ref(rad, &lsat.band[i]);
+ if (ref < 0. && method > DOS) ref = 0.;
+ }
+ }
+ ((DCELL *) outrast)[col] = ref;
+ }
+ }
+ if (G_put_raster_row(outfd, outrast, DCELL_TYPE) < 0)
+ G_fatal_error(_("Cannot write to <%s>"), band_out);
+ }
+ ref_mode = 0.;
+ if (method > DOS && !lsat.band[i].thermal)
+ {
+ ref_mode = lsat_qcal2rad(dn_mode[i], &lsat.band[i]);
+ ref_mode = lsat_rad2ref(ref_mode, &lsat.band[i]);
+ }
- if (G_put_raster_row(outfd, outrast, DCELL_TYPE) < 0)
- G_fatal_error(_("Cannot write to <%s>"), band_out);
- }
/* ================================================================= */
- G_free(inrast);
- G_close_cell(infd);
- G_free(outrast);
- G_close_cell(outfd);
+ G_free(inrast);
+ G_close_cell(infd);
+ G_free(outrast);
+ G_close_cell(outfd);
- sprintf(command, "r.colors map=%s color=grey", band_out);
- system(command);
+ char command[300];
+ sprintf(command, "r.colors map=%s color=grey", band_out);
+ system(command);
- G_short_history(band_out, "raster", &history);
- G_command_history(&history);
- G_write_history(band_out, &history);
- }
+ G_short_history(band_out, "raster", &history);
- G_set_window(&window);
- exit(EXIT_SUCCESS);
+ sprintf(history.edhist[0], " %s of Landsat-%d %s (method %s)", lsat.band[i].thermal ? "Temperature": "Reflectance",
+ lsat.number, lsat.sensor, metho->answer);
+ sprintf(history.edhist[1], "----------------------------------------------------------------");
+ sprintf(history.edhist[2], " Acquisition date ...................... %s", lsat.date);
+ sprintf(history.edhist[3], " Production date ....................... %s\n", lsat.creation);
+ sprintf(history.edhist[4], " Earth-sun distance (d) ................ %.8lf", lsat.dist_es);
+ sprintf(history.edhist[5], " 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", lsat.band[i].lmin, lsat.band[i].lmax);
+ sprintf(history.edhist[7], " 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], " Temperature (K1 and K2) ............... %.3lf and %.3lf", lsat.band[i].K1, lsat.band[i].K2);
+ history.edlinecnt = 9;
+ } else {
+ sprintf(history.edhist[8], " Mean solar irradiance (ESUN) .......... %.3lf", lsat.band[i].esun);
+ sprintf(history.edhist[9], " Reflectance = Radiance divided by ..... %.5lf", lsat.band[i].K2);
+ history.edlinecnt = 10;
+ if (method > DOS) {
+ sprintf(history.edhist[10], "");
+ sprintf(history.edhist[11], " Dark object (%4d pixels) DN = ........ %d", pixel, dn_dark[i]);
+ sprintf(history.edhist[12], " Mode in reflectance histogram ......... %.5lf", ref_mode);
+ history.edlinecnt = 13;
+ }
+ }
+ sprintf (history.edhist[history.edlinecnt], "-----------------------------------------------------------------");
+ history.edlinecnt++;
+
+ G_command_history(&history);
+ G_write_history(band_out, &history);
+ }
+
+ exit(EXIT_SUCCESS);
}
More information about the grass-commit
mailing list