[GRASS-SVN] r30494 - in grass-addons/imagery: . i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Mar 7 12:05:37 EST 2008
Author: ejtizado
Date: 2008-03-07 12:05:37 -0500 (Fri, 07 Mar 2008)
New Revision: 30494
Added:
grass-addons/imagery/i.landsat.toar/
grass-addons/imagery/i.landsat.toar/description.html
grass-addons/imagery/i.landsat.toar/landsat.c
grass-addons/imagery/i.landsat.toar/landsat.h
grass-addons/imagery/i.landsat.toar/landsat_set.c
grass-addons/imagery/i.landsat.toar/main.c
Removed:
grass-addons/imagery/i.landsat.toar/description.html
grass-addons/imagery/i.landsat.toar/landsat.c
grass-addons/imagery/i.landsat.toar/landsat.h
grass-addons/imagery/i.landsat.toar/landsat_set.c
grass-addons/imagery/i.landsat.toar/main.c
Log:
Move to imaginery/
Copied: grass-addons/imagery/i.landsat.toar (from rev 30401, grass-addons/i.landsat.toar)
Deleted: grass-addons/imagery/i.landsat.toar/description.html
===================================================================
--- grass-addons/i.landsat.toar/description.html 2008-02-29 15:19:34 UTC (rev 30401)
+++ grass-addons/imagery/i.landsat.toar/description.html 2008-03-07 17:05:37 UTC (rev 30494)
@@ -1,128 +0,0 @@
-<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.
-
-<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>
-
-<p>Optionally, the data can be read from header file (.met) for all Landsat MSS, TM and ETM+. However, if the solar elevation or the product creation date are given the values of metfile are overwriten. This is necessary when this data in metfile are incorrects or imprecisses.</p>
-
-<p><b>Attention</b>: Any null value or less than QCALmin in input raster is set to null in output raster and it is not included in the equations.</p>
-
-
-<H2>Uncorrected at-sensor values (method=uncorrected, default)</H2>
-
-<p>The standard geometric and radiometric correction results in a calibrated digital number (QCAL = DN) images. To further standardize the impact of illumination geometry, the QCAL images are converted first to at-sensor radiance and then to at-sensor reflectance. The thermal band is first converted from QCAL to at-sensor radiance, and then to effective at-sensor temperature in Kelvin degrees.</p>
-
-<p>Radiometric calibration convert QCAL to <b>at-sensor radiance</b>, a radiometric quantity measured in W/(m²·sr·µm) with the equations:
- <ul>
- <li> gain = (Lmax - Lmin) / (QCALmax - QCALmin)</li>
- <li> bias = Lmin - gain · QCALmin </li>
- <li> radiance = gain · QCAL + bias </li>
- </ul>
-where,
-<em>Lmax</em> and <em>Lmin</em> are the calibration constants, and
-<em>QCALmax</em> and <em>QCALmin</em> are the highest and the lowest points of the range of rescaled radiance in QCAL.
-</p>
-
-<p>Then, to calculate <b>at-sensor reflectance</b> the equations are:
- <ul>
- <li> sun_radiance = [Esun · sin(e)] / (PI · d^2)</li>
- <li> reflectance = radiance / sun_radiance </li>
- </ul>
-where,
-<em>d</em> is the earth-sun distance in astronomical units,
-<em>e</em> is the solar elevation angle, and
-<em>Esun</em> is the mean solar exoatmospheric irradiance in W/(m²·µm).
-</p>
-
-
-<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:
- <ul>
- <li> radiance = (uncorrected radiance - Lmin) </li>
- <li> reflectance = radiance / sun_radiance </li>
- </ul>
-</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:
- <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> reflectance = radiance / sun_radiance </li>
- </ul>
-where,
-<em>percent</em> is a value between 0.0 and 1.0 (usually 0.01),
-<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>
-
-<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>
-
-
-<H2>NOTES</H2>
-
-<p>In verbose mode (flag -v), the program write basic data of satellite and the parameters used in
-transformations.</p>
-
-<p>Production date is not a exact value and it is necessary to apply correct calibration constants, which were changed in the dates:
- <ul>
- <li>Landsat-1 MSS: never </li>
- <li>Landsat-2 MSS: July 16, 1975</li>
- <li>Landsat-3 MSS: June 1, 1978</li>
- <li>Landsat-4 MSS: August 26, 1982 and April 1, 1983</li>
- <li>Landsat-4 TM: August 1, 1983 and January 15, 1984</li>
- <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>
- </ul>
-</p>
-
-<H2>EXAMPLES</H2>
-
-<p>Transform digital numbers of Landsat-7 ETM+ in band rasters 203_30.1, 203_30.2
-[...] to uncorrected at-sensor reflectance in output files 203_30.toar.1,
-203_30.toar.2 [...] and at-sensor temperature in output files 293_39.toar.61 and 293_39.toar.62:</p>
-
-<div class="code"><pre>
-i.landsat.toar -7 band=203_30 met=p203r030_7x20010620.met
-</pre></div>
-
-<p>or</p>
-
-<div class="code"><pre>
-i.landsat.toar -7 band=203_30 product=2004-06-07 date=2001-06-20 solar=64.3242970 gain="HHHLHLHHL"
-</pre></div>
-
-<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>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>
-
-</ol>
-
-
-<H2>SEE ALSO</H2>
-
-<em>
-<A HREF="r.mapcalc.html">r.mapcalc</A><br>
-<A HREF="r.in.gdal.html">r.in.gdal</A><br>
-</em>
-
-
-<H2>AUTHOR</H2>
-
-E. Jorge Tizado (ej.tizado unileon es)<br>
-Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
-
-<p>
-<i>Last changed: $Date: 2007/07/07 00:00:00 $</i>
Copied: grass-addons/imagery/i.landsat.toar/description.html (from rev 30464, grass-addons/i.landsat.toar/description.html)
===================================================================
--- grass-addons/imagery/i.landsat.toar/description.html (rev 0)
+++ grass-addons/imagery/i.landsat.toar/description.html 2008-03-07 17:05:37 UTC (rev 30494)
@@ -0,0 +1,139 @@
+<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 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>
+
+<p>Optionally, the data can be read from header file (.met) for all Landsat MSS, TM and ETM+. However, if the solar elevation or the product creation date are given the values of metfile are overwriten. This is necessary when this data in metfile are incorrects or imprecisses.</p>
+
+<p><b>Attention</b>: Any null value or less than QCALmin in input raster is set to null in output raster and it is not included in the equations.</p>
+
+
+<H2>Uncorrected at-sensor values (method=uncorrected, default)</H2>
+
+<p>The standard geometric and radiometric correction results in a calibrated digital number (QCAL = DN) images. To further standardize the impact of illumination geometry, the QCAL images are converted first to at-sensor radiance and then to at-sensor reflectance. The thermal band is first converted from QCAL to at-sensor radiance, and then to effective at-sensor temperature in Kelvin degrees.</p>
+
+<p>Radiometric calibration convert QCAL to <b>at-sensor radiance</b>, a radiometric quantity measured in W/(m²·sr·µm) with the equations:
+ <ul>
+ <li> gain = (Lmax - Lmin) / (QCALmax - QCALmin)</li>
+ <li> bias = Lmin - gain · QCALmin </li>
+ <li> radiance = gain · QCAL + bias </li>
+ </ul>
+where,
+<em>Lmax</em> and <em>Lmin</em> are the calibration constants, and
+<em>QCALmax</em> and <em>QCALmin</em> are the highest and the lowest points of the range of rescaled radiance in QCAL.
+</p>
+
+<p>Then, to calculate <b>at-sensor reflectance</b> the equations are:
+ <ul>
+ <li> sun_radiance = [Esun · sin(e)] / (PI · d^2)</li>
+ <li> reflectance = radiance / sun_radiance </li>
+ </ul>
+where,
+<em>d</em> is the earth-sun distance in astronomical units,
+<em>e</em> is the solar elevation angle, and
+<em>Esun</em> is the mean solar exoatmospheric irradiance in W/(m²·µm).
+</p>
+
+
+<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 (not for thermal bands):
+ <ul>
+ <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=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 = 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,
+<em>percent</em> is a value between 0.0 and 1.0 (usually 0.01),
+<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.
+<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>
+
+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
+transformations.</p>
+
+<p>Production date is not a exact value and it is necessary to apply correct calibration constants, which were changed in the dates:
+ <ul>
+ <li>Landsat-1 MSS: never </li>
+ <li>Landsat-2 MSS: July 16, 1975</li>
+ <li>Landsat-3 MSS: June 1, 1978</li>
+ <li>Landsat-4 MSS: August 26, 1982 and April 1, 1983</li>
+ <li>Landsat-4 TM: August 1, 1983 and January 15, 1984</li>
+ <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>
+ </ul>
+</p>
+
+<H2>EXAMPLES</H2>
+
+<p>Transform digital numbers of Landsat-7 ETM+ in band rasters 203_30.1, 203_30.2
+[...] to uncorrected at-sensor reflectance in output files 203_30.toar.1,
+203_30.toar.2 [...] and at-sensor temperature in output files 293_39.toar.61 and 293_39.toar.62:</p>
+
+<div class="code"><pre>
+i.landsat.toar -7 band=203_30 met=p203r030_7x20010620.met
+</pre></div>
+
+<p>or</p>
+
+<div class="code"><pre>
+i.landsat.toar -7 band=203_30 product=2004-06-07 date=2001-06-20 solar=64.3242970 gain="HHHLHLHHL"
+</pre></div>
+
+<H2>REFERENCES</H2>
+<ol>
+ <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>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>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="r.mapcalc.html">r.mapcalc</A><br>
+<A HREF="r.in.gdal.html">r.in.gdal</A><br>
+</em>
+
+
+<H2>AUTHOR</H2>
+
+E. Jorge Tizado (ej.tizado unileon es)<br>
+Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
+
+<p>
+<i>Last changed: $Date: 2008/02/10 00:00:00 $</i>
Deleted: grass-addons/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass-addons/i.landsat.toar/landsat.c 2008-02-29 15:19:34 UTC (rev 30401)
+++ grass-addons/imagery/i.landsat.toar/landsat.c 2008-03-07 17:05:37 UTC (rev 30494)
@@ -1,96 +0,0 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<math.h>
-
-#include "landsat.h"
-
-#define PI 3.1415926535897932384626433832795
-#define R2D 57.295779513082320877
-#define D2R 0.017453292519943295769
-
-/****************************************************************************
- * PURPOSE: Calibrated Digital Number to at-satellite Radiance
- *****************************************************************************/
-double lsat_qcal2rad(int qcal, band_data * band)
-{
- return (double)(((double)qcal) * band->gain + band->bias);
-}
-
-/****************************************************************************
- * PURPOSE: Radiance of non-thermal band to at-satellite Reflectance
- *****************************************************************************/
-double lsat_rad2ref(double rad, band_data * band)
-{
- return (double)(rad / band->K2);
-}
-
-/****************************************************************************
- * PURPOSE: Radiance of thermal band to at-satellite Temperature
- *****************************************************************************/
-double lsat_rad2temp(double rad, band_data * band)
-{
- return (double)(band->K2 / log((band->K1 / rad) + 1.0));
-}
-
-/****************************************************************************
- * PURPOSE: Some band constants
- *
- * zenith = 90 - sun_elevation
- * sin( sun_elevation ) = cos( sun_zenith )
- *
- * 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
- *****************************************************************************/
-
-void lsat_bandctes(lsat_data * lsat, int i, char method, double aot, double percent)
-{
- double a, b, rad_sun;
- double TAUv, TAUz, Edown;
-
- a = (double)(PI * lsat->dist_es * lsat->dist_es);
- b = (double)(sin(D2R * lsat->sun_elev));
-
- /** 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);
-
- /** 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);
-
- 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;
- }
-
- /** 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;
- }
-}
-
Copied: grass-addons/imagery/i.landsat.toar/landsat.c (from rev 30464, grass-addons/i.landsat.toar/landsat.c)
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat.c (rev 0)
+++ grass-addons/imagery/i.landsat.toar/landsat.c 2008-03-07 17:05:37 UTC (rev 30494)
@@ -0,0 +1,164 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+
+#include "landsat.h"
+
+#define PI 3.1415926535897932384626433832795
+#define R2D 57.295779513082320877
+#define D2R 0.017453292519943295769
+
+/****************************************************************************
+ * PURPOSE: Calibrated Digital Number to at-satellite Radiance
+ *****************************************************************************/
+double lsat_qcal2rad(double qcal, band_data * band)
+{
+ return (double)(qcal * band->gain + band->bias);
+}
+
+/****************************************************************************
+ * PURPOSE: Radiance of non-thermal band to at-satellite Reflectance
+ *****************************************************************************/
+double lsat_rad2ref(double rad, band_data * band)
+{
+ return (double)(rad / band->K2);
+}
+
+/****************************************************************************
+ * PURPOSE: Radiance of thermal band to at-satellite Temperature
+ *****************************************************************************/
+double lsat_rad2temp(double rad, band_data * band)
+{
+ return (double)(band->K2 / log((band->K1 / rad) + 1.0));
+}
+
+/****************************************************************************
+ * PURPOSE: Some band constants
+ *
+ * zenith = 90 - sun_elevation
+ * sin( sun_elevation ) = cos( sun_zenith )
+ *
+ * lsat : satellite data
+ * i : band number
+ * method : level of atmospheric correction
+ * percent : percent of solar irradiance in path radiance
+ * dos : digital number of dark 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 sat_zenith, double rayleigh)
+{
+ 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;
+
+ 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));
+
+ /** 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);
+
+ lsat->band[i].K1 = 0.;
+ lsat->band[i].K2 = rad_sun;
+ }
+
+ /** 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));
+
+ 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;
+ }
+ }
+}
+
+
+
+
Deleted: grass-addons/imagery/i.landsat.toar/landsat.h
===================================================================
--- grass-addons/i.landsat.toar/landsat.h 2008-02-29 15:19:34 UTC (rev 30401)
+++ grass-addons/imagery/i.landsat.toar/landsat.h 2008-03-07 17:05:37 UTC (rev 30494)
@@ -1,62 +0,0 @@
-#ifndef _LANDSAT_H
-#define _LANDSAT_H
-
-#define MAX_BANDS 9
-
-#define UNCORRECTED 0
-#define CORRECTED 1
-#define SIMPLIFIED 2
-
-
-/*****************************************************
- * Landsat Structures
- *
- * Lmax and Lmin in W / (m² · sr · µm) -> Radiance
- * Esun in W / (m² · µm) -> Irradiance
- ****************************************************/
-
-typedef struct
-{
- int number; /* Band number */
- int code; /* Band code */
-
- double wavemax, wavemin; /* Wavelength in µm */
-
- 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 */
-
-} band_data;
-
-typedef struct
-{
- 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 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;
-
-
-/*****************************************************************************
- * Landsat Equations Prototypes
- *****************************************************************************/
-
-double lsat_qcal2rad(int, band_data *);
-double lsat_rad2ref(double, band_data *);
-double lsat_rad2temp(double, band_data *);
-
-void lsat_bandctes(lsat_data *, int, char, double, double);
-
-#endif
Copied: grass-addons/imagery/i.landsat.toar/landsat.h (from rev 30464, grass-addons/i.landsat.toar/landsat.h)
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat.h (rev 0)
+++ grass-addons/imagery/i.landsat.toar/landsat.h 2008-03-07 17:05:37 UTC (rev 30494)
@@ -0,0 +1,66 @@
+#ifndef _LANDSAT_H
+#define _LANDSAT_H
+
+#define UNCORRECTED 0
+#define CORRECTED 1
+#define DOS 10
+#define DOS1 12
+#define DOS2 14
+#define DOS2b 15
+#define DOS3 16
+#define DOS4 18
+
+
+/*****************************************************
+ * Landsat Structures
+ *
+ * Lmax and Lmin in W / (m² · sr · µm) -> Radiance
+ * Esun in W / (m² · µm) -> Irradiance
+ ****************************************************/
+
+#define MAX_BANDS 9
+
+typedef struct
+{
+ int number; /* Band number */
+ int code; /* Band code */
+
+ double wavemax, wavemin; /* Wavelength in µm */
+
+ 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 */
+
+} band_data;
+
+typedef struct
+{
+ 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 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;
+
+
+/*****************************************************************************
+ * Landsat Equations Prototypes
+ *****************************************************************************/
+
+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, int, double, double);
+
+#endif
Deleted: grass-addons/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass-addons/i.landsat.toar/landsat_set.c 2008-02-29 15:19:34 UTC (rev 30401)
+++ grass-addons/imagery/i.landsat.toar/landsat_set.c 2008-03-07 17:05:37 UTC (rev 30494)
@@ -1,437 +0,0 @@
-#include<stdio.h>
-#include<stdlib.h>
-#include<string.h>
-#include <math.h>
-
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-#include "local_proto.h"
-#include "earth_sun.h"
-
-void sensor_MSS(lsat_data * lsat)
-{
- 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 };
-
- 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;
-}
-
-void sensor_TM(lsat_data * lsat)
-{
- 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 };
-
- 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;
-}
-
-void sensor_ETM(lsat_data * lsat)
-{
- 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 };
-
- 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;
-}
-
-
-/** **********************************************
- ** Before access to this function ...
- ** store previously
- ** >>> adquisition date,
- ** >>> creation date, and
- ** >>> sun_elev
- ** **********************************************/
-
-/****************************************************************************
- * PURPOSE: Store values of Landsat-1 MSS
- * July 23, 1972 to January 6, 1978
- *****************************************************************************/
-void set_MSS1(lsat_data * lsat)
-{
- int i, j;
-
- /** Brian L. Markham and John L. Barker.
- EOSAT Landsat Technical Notes, No. 1, 1986 */
- /* Spectral radiances at detector */
- double lmax[] = { 248., 200., 176., 153. };
- double lmin[] = { 0., 0., 0., 0. };
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1852., 1584., 1276., 904. };
-
- lsat->number = 1;
- sensor_MSS( 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);
- }
- return;
-}
-
-/****************************************************************************
- * PURPOSE: Store values of Landsat-2 MSS
- * January 22, 1975 to February 25, 1982
- *****************************************************************************/
-void set_MSS2(lsat_data * lsat)
-{
- int i, j;
- double julian, *lmax, *lmin;
-
- /** Brian L. Markham and John L. Barker.
- EOSAT Landsat Technical Notes, No. 1, 1986 */
- /* Spectral radiances at detector */
- double Lmax[][4] = { { 210., 156., 140., 138. }, /* before July 16, 1975 */
- { 263., 176., 152., 130. } }; /* on or after July 16, 1975 */
- double Lmin[][4] = { { 10., 7., 7., 5. },
- { 8., 6., 6., 4. } };
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1856., 1559., 1269., 906. };
-
- julian = julian_char(lsat->creation);
- if (julian < julian_char("1975-07-16")) i = 0;
- else i = 1;
- lmax = Lmax[i];
- lmin = Lmin[i];
-
- lsat->number = 2;
- sensor_MSS( 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);
- }
- return;
-}
-
-/****************************************************************************
- * PURPOSE: Store values of Landsat-3 MSS
- * March 5, 1978 to March 31, 1983
- *
- * tiene una banda 8 thermal
- *****************************************************************************/
-void set_MSS3(lsat_data * lsat)
-{
- int i, j;
- double julian, *lmax, *lmin;
-
- /** Brian L. Markham and John L. Barker.
- EOSAT Landsat Technical Notes, No. 1, 1986 */
- /* Spectral radiances at detector */
- double Lmax[][4] = { { 220., 175., 145., 147. }, /* before June 1, 1978 */
- { 259., 179., 149., 128. } }; /* on or after June 1, 1978 */
- double Lmin[][4] = { { 4., 3., 3., 1. },
- { 4., 3., 3., 1. } };
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1860., 1571., 1289., 910. };
-
- julian = julian_char(lsat->creation);
- if (julian < julian_char("1978-06-01")) i = 0;
- else i = 1;
- lmax = Lmax[i];
- lmin = Lmin[i];
-
- lsat->number = 3;
- sensor_MSS( 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);
- }
- return;
-}
-
-/****************************************************************************
- * PURPOSE: Store values of Landsat-4 MSS/TM
- * July 16, 1982 to June 15, 2001
- *****************************************************************************/
-void set_MSS4(lsat_data * lsat)
-{
- int i, j;
- double julian, *lmax, *lmin;
-
- /** Brian L. Markham and John L. Barker.
- EOSAT Landsat Technical Notes, No. 1, 1986 */
- /* Spectral radiances at detector */
- double Lmax[][4] = { { 250., 180., 150., 133. }, /* before August 26, 1982 */
- { 230., 180., 130., 133. }, /* between */
- { 238., 164., 142., 116. } }; /* on or after April 1, 1983 */
- double Lmin[][4] = { { 2., 4., 4., 3. },
- { 2., 4., 4., 3. },
- { 4., 4., 5., 4. } };
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1851., 1593., 1260., 878. };
-
- julian = julian_char(lsat->creation);
- if (julian < julian_char("1982-08-26")) i = 0;
- else if (julian < julian_char("1983-03-31")) i = 1;
- else i = 2;
- lmax = Lmax[i];
- lmin = Lmin[i];
-
- lsat->number = 4;
- sensor_MSS( 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);
- }
- return;
-}
-
-void set_TM4(lsat_data * lsat)
-{
- int i, j;
- double julian, *lmax, *lmin;
-
- /** Brian L. Markham and John L. Barker.
- EOSAT Landsat Technical Notes, No. 1, 1986 */
- /* Spectral radiances at detector */
- double Lmax[][7] = { { 158.42, 308.17, 234.63, 224.32, 32.42, 15.64, 17.00 }, /* before August 1983 */
- { 142.86, 291.25, 225.00, 214.29, 30.00, 12.40, 15.93 }, /* before January 15, 1984 */
- { 152.10, 296.80, 204.30, 206.20, 27.19, 15.60, 14.38 } }; /* after Jaunary 15, 1984 */
- double Lmin[][7] = { { -1.52, -2.84, -1.17, -1.51, -0.37, 2.00, -0.15 },
- { 0.00, 0.00, 0.00, 0.00, 0.00, 4.84, 0.00 },
- { -1.50, -2.80, -1.20, -1.50, -0.37, 1.238, -0.15 } };
- /** Gyanesh Chander and Brian Markham.
- IEEE Transactions On Geoscience And Remote Sensing, Vol. 41, No. 11, November 2003 */
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1957., 1825., 1557., 1033., 214.9, 0., 80.72 };
- /* Thermal band calibration constants: K1 = 671.62 K2 = 1284.30 */
-
- julian = julian_char(lsat->creation);
- if (julian < julian_char("1983-08-01")) i = 0;
- else if (julian < julian_char("1984-01-15")) i = 1;
- else i = 2;
- lmax = Lmax[i];
- lmin = Lmin[i];
-
- lsat->number = 4;
- sensor_TM( 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 = 671.62;
- lsat->band[i].K2 = 1284.30;
- }
- }
- return;
-}
-
-
-/****************************************************************************
- * PURPOSE: Store values of Landsat-5 MSS/TM
- * March 1, 1984 to today
- *****************************************************************************/
-void set_MSS5(lsat_data * lsat)
-{
- int i, j;
- double julian, *lmax, *lmin;
-
- /** Brian L. Markham and John L. Barker.
- EOSAT Landsat Technical Notes, No. 1, 1986 */
- /* Spectral radiances at detector */
- double Lmax[][4] = { { 240., 170., 150., 127. }, /* before April 6, 1984 */
- { 268., 179., 159., 123. }, /* betweeen */
- { 268., 179., 148., 123. } }; /* after November 9, 1984 */
- double Lmin[][4] = { { 4., 3., 4., 2. },
- { 3., 3., 4., 3. },
- { 3., 3., 5., 3. } };
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1849., 1595., 1253., 870. };
-
- julian = julian_char(lsat->creation);
- if (julian < julian_char("1984-04-06")) i = 0;
- else if (julian < julian_char("1984-11-09")) i = 1;
- else i = 2;
- lmax = Lmax[i];
- lmin = Lmin[i];
-
- lsat->number = 5;
- sensor_MSS( 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);
- }
- return;
-}
-
-void set_TM5(lsat_data * lsat)
-{
- int i, j;
- double julian, *lmax, *lmin;
-
- /** 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 */
- { 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 },
- { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 },
- { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 } };
- /* Solar exoatmospheric spectral irradiances */
- double esun[] = { 1957., 1826., 1554., 1036., 215.0, 0., 80.67 };
- /* Thermal band calibration constants: K1 = 607.76 K2 = 1260.56 */
-
- julian = julian_char(lsat->creation);
- if (julian < julian_char("2003-05-04")) i = 0;
- else if (julian < julian_char("2007-04-02")) i = 1;
- else i = 2;
- lmax = Lmax[i];
- lmin = Lmin[i];
- if ( i == 2 ) {
- julian = julian_char(lsat->date); /* Yes, here acquisition date */
- if (julian >= julian_char("1992-01-01")) {
- lmax[0] = 193.0;
- lmax[1] = 365.0;
- }
- }
-
- lsat->number = 5;
- sensor_TM( 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 = 607.76;
- lsat->band[i].K2 = 1260.56;
- }
- }
- return;
-}
-
-
-/****************************************************************************
- * PURPOSE: Store values of Landsat-7 ETM+
- * April 15, 1999 to today
- *****************************************************************************/
-void set_ETM(lsat_data * lsat, char gain[])
-{
- int i, k, j;
- double julian, *lmax, *lmin;
-
- /** Richard Irish.
- Landsat 7. Science Data Users Handbook. Last update: February 17, 2007 */
- /* Spectral radiances at detector */
- /* - LOW GAIN - */
- double LmaxL[][8] = { { 297.5, 303.4, 235.5, 235.0, 47.70, 17.04, 16.600, 244.0 }, /* before July 1, 2000 */
- { 293.7, 300.9, 234.4, 241.1, 47.57, 17.04, 16.540, 243.1 } }; /* on or after July 1, 2000 */
- double LminL[][8] = { { -6.2, -6.0, -4.5, -4.5, -1.0, 0.0, -0.35, -5.0 },
- { -6.2, -6.4, -5.0, -5.1, -1.0, 0.0, -0.35, -4.7 } };
- /* - HIGH GAIN - */
- double LmaxH[][8] = { { 194.3, 202.4, 158.6, 157.5, 31.76, 12.65, 10.932, 158.4 },
- { 191.6, 196.5, 152.9, 157.4, 31.06, 12.65, 10.800, 158.3 } };
- double LminH[][8] = { { -6.2, -6.0, -4.5, -4.5, -1.0, 3.2, -0.35, -5.0 },
- { -6.2, -6.4, -5.0, -5.1, -1.0, 3.2, -0.35, -4.7 } };
- /* Solar exoatmospheric spectral irradiances */
- 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;
-
- lsat->number = 7;
- sensor_ETM( 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);
- 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;
-}
Copied: grass-addons/imagery/i.landsat.toar/landsat_set.c (from rev 30464, grass-addons/i.landsat.toar/landsat_set.c)
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat_set.c (rev 0)
+++ grass-addons/imagery/i.landsat.toar/landsat_set.c 2008-03-07 17:05:37 UTC (rev 30494)
@@ -0,0 +1,443 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+#include <math.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+#include "earth_sun.h"
+
+void sensor_MSS(lsat_data * lsat)
+{
+ int i;
+
+ /* 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");
+
+ 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;
+
+ /* 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");
+
+ 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;
+
+ /* 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+");
+
+ 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;
+}
+
+
+/** **********************************************
+ ** Before access to this function ...
+ ** store previously
+ ** >>> adquisition date,
+ ** >>> creation date, and
+ ** >>> sun_elev
+ ** **********************************************/
+
+/****************************************************************************
+ * PURPOSE: Store values of Landsat-1 MSS
+ * July 23, 1972 to January 6, 1978
+ *****************************************************************************/
+void set_MSS1(lsat_data * lsat)
+{
+ int i, j;
+
+ /** Brian L. Markham and John L. Barker.
+ EOSAT Landsat Technical Notes, No. 1, 1986 */
+ /* Spectral radiances at detector */
+ double lmax[] = { 248., 200., 176., 153. };
+ double lmin[] = { 0., 0., 0., 0. };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1852., 1584., 1276., 904. };
+
+ lsat->number = 1;
+ sensor_MSS( 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);
+ }
+ return;
+}
+
+/****************************************************************************
+ * PURPOSE: Store values of Landsat-2 MSS
+ * January 22, 1975 to February 25, 1982
+ *****************************************************************************/
+void set_MSS2(lsat_data * lsat)
+{
+ int i, j;
+ double julian, *lmax, *lmin;
+
+ /** Brian L. Markham and John L. Barker.
+ EOSAT Landsat Technical Notes, No. 1, 1986 */
+ /* Spectral radiances at detector */
+ double Lmax[][4] = { { 210., 156., 140., 138. }, /* before July 16, 1975 */
+ { 263., 176., 152., 130. } }; /* on or after July 16, 1975 */
+ double Lmin[][4] = { { 10., 7., 7., 5. },
+ { 8., 6., 6., 4. } };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1856., 1559., 1269., 906. };
+
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("1975-07-16")) i = 0;
+ else i = 1;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
+
+ lsat->number = 2;
+ sensor_MSS( 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);
+ }
+ return;
+}
+
+/****************************************************************************
+ * PURPOSE: Store values of Landsat-3 MSS
+ * March 5, 1978 to March 31, 1983
+ *
+ * tiene una banda 8 thermal
+ *****************************************************************************/
+void set_MSS3(lsat_data * lsat)
+{
+ int i, j;
+ double julian, *lmax, *lmin;
+
+ /** Brian L. Markham and John L. Barker.
+ EOSAT Landsat Technical Notes, No. 1, 1986 */
+ /* Spectral radiances at detector */
+ double Lmax[][4] = { { 220., 175., 145., 147. }, /* before June 1, 1978 */
+ { 259., 179., 149., 128. } }; /* on or after June 1, 1978 */
+ double Lmin[][4] = { { 4., 3., 3., 1. },
+ { 4., 3., 3., 1. } };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1860., 1571., 1289., 910. };
+
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("1978-06-01")) i = 0;
+ else i = 1;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
+
+ lsat->number = 3;
+ sensor_MSS( 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);
+ }
+ return;
+}
+
+/****************************************************************************
+ * PURPOSE: Store values of Landsat-4 MSS/TM
+ * July 16, 1982 to June 15, 2001
+ *****************************************************************************/
+void set_MSS4(lsat_data * lsat)
+{
+ int i, j;
+ double julian, *lmax, *lmin;
+
+ /** Brian L. Markham and John L. Barker.
+ EOSAT Landsat Technical Notes, No. 1, 1986 */
+ /* Spectral radiances at detector */
+ double Lmax[][4] = { { 250., 180., 150., 133. }, /* before August 26, 1982 */
+ { 230., 180., 130., 133. }, /* between */
+ { 238., 164., 142., 116. } }; /* on or after April 1, 1983 */
+ double Lmin[][4] = { { 2., 4., 4., 3. },
+ { 2., 4., 4., 3. },
+ { 4., 4., 5., 4. } };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1851., 1593., 1260., 878. };
+
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("1982-08-26")) i = 0;
+ else if (julian < julian_char("1983-03-31")) i = 1;
+ else i = 2;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
+
+ lsat->number = 4;
+ sensor_MSS( 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);
+ }
+ return;
+}
+
+void set_TM4(lsat_data * lsat)
+{
+ int i, j;
+ double julian, *lmax, *lmin;
+
+ /** Brian L. Markham and John L. Barker.
+ EOSAT Landsat Technical Notes, No. 1, 1986 */
+ /* Spectral radiances at detector */
+ double Lmax[][7] = { { 158.42, 308.17, 234.63, 224.32, 32.42, 15.64, 17.00 }, /* before August 1983 */
+ { 142.86, 291.25, 225.00, 214.29, 30.00, 12.40, 15.93 }, /* before January 15, 1984 */
+ { 152.10, 296.80, 204.30, 206.20, 27.19, 15.60, 14.38 } }; /* after Jaunary 15, 1984 */
+ double Lmin[][7] = { { -1.52, -2.84, -1.17, -1.51, -0.37, 2.00, -0.15 },
+ { 0.00, 0.00, 0.00, 0.00, 0.00, 4.84, 0.00 },
+ { -1.50, -2.80, -1.20, -1.50, -0.37, 1.238, -0.15 } };
+ /** Gyanesh Chander and Brian Markham.
+ IEEE Transactions On Geoscience And Remote Sensing, Vol. 41, No. 11, November 2003 */
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1957., 1825., 1557., 1033., 214.9, 0., 80.72 };
+ /* Thermal band calibration constants: K1 = 671.62 K2 = 1284.30 */
+
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("1983-08-01")) i = 0;
+ else if (julian < julian_char("1984-01-15")) i = 1;
+ else i = 2;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
+
+ lsat->number = 4;
+ sensor_TM( 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 = 671.62;
+ lsat->band[i].K2 = 1284.30;
+ }
+ }
+ return;
+}
+
+
+/****************************************************************************
+ * PURPOSE: Store values of Landsat-5 MSS/TM
+ * March 1, 1984 to today
+ *****************************************************************************/
+void set_MSS5(lsat_data * lsat)
+{
+ int i, j;
+ double julian, *lmax, *lmin;
+
+ /** Brian L. Markham and John L. Barker.
+ EOSAT Landsat Technical Notes, No. 1, 1986 */
+ /* Spectral radiances at detector */
+ double Lmax[][4] = { { 240., 170., 150., 127. }, /* before April 6, 1984 */
+ { 268., 179., 159., 123. }, /* betweeen */
+ { 268., 179., 148., 123. } }; /* after November 9, 1984 */
+ double Lmin[][4] = { { 4., 3., 4., 2. },
+ { 3., 3., 4., 3. },
+ { 3., 3., 5., 3. } };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1849., 1595., 1253., 870. };
+
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("1984-04-06")) i = 0;
+ else if (julian < julian_char("1984-11-09")) i = 1;
+ else i = 2;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
+
+ lsat->number = 5;
+ sensor_MSS( 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);
+ }
+ return;
+}
+
+void set_TM5(lsat_data * lsat)
+{
+ int i, j;
+ double julian, *lmax, *lmin;
+
+ /** 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 }, /* 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 },
+ { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 },
+ { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 } };
+ /* Solar exoatmospheric spectral irradiances */
+ double esun[] = { 1957., 1826., 1554., 1036., 215.0, 0., 80.67 };
+ /* Thermal band calibration constants: K1 = 607.76 K2 = 1260.56 */
+
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("2003-05-04")) i = 0;
+ else if (julian < julian_char("2007-04-02")) i = 1;
+ else i = 2;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
+ 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;
+ lmax[1] = 365.0;
+ }
+ }
+
+ lsat->number = 5;
+ sensor_TM( 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 = 607.76;
+ lsat->band[i].K2 = 1260.56;
+ }
+ }
+ return;
+}
+
+
+/****************************************************************************
+ * PURPOSE: Store values of Landsat-7 ETM+
+ * April 15, 1999 to today
+ *****************************************************************************/
+void set_ETM(lsat_data * lsat, char gain[])
+{
+ int i, k, j;
+ double julian, *lmax, *lmin;
+
+ /** Richard Irish.
+ Landsat 7. Science Data Users Handbook. Last update: February 17, 2007 */
+ /* Spectral radiances at detector */
+ /* - LOW GAIN - */
+ double LmaxL[][8] = { { 297.5, 303.4, 235.5, 235.0, 47.70, 17.04, 16.600, 244.0 }, /* before July 1, 2000 */
+ { 293.7, 300.9, 234.4, 241.1, 47.57, 17.04, 16.540, 243.1 } }; /* on or after July 1, 2000 */
+ double LminL[][8] = { { -6.2, -6.0, -4.5, -4.5, -1.0, 0.0, -0.35, -5.0 },
+ { -6.2, -6.4, -5.0, -5.1, -1.0, 0.0, -0.35, -4.7 } };
+ /* - HIGH GAIN - */
+ double LmaxH[][8] = { { 194.3, 202.4, 158.6, 157.5, 31.76, 12.65, 10.932, 158.4 },
+ { 191.6, 196.5, 152.9, 157.4, 31.06, 12.65, 10.800, 158.3 } };
+ double LminH[][8] = { { -6.2, -6.0, -4.5, -4.5, -1.0, 3.2, -0.35, -5.0 },
+ { -6.2, -6.4, -5.0, -5.1, -1.0, 3.2, -0.35, -4.7 } };
+ /* Solar exoatmospheric spectral irradiances */
+ 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;
+
+ lsat->number = 7;
+ sensor_ETM( 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);
+ 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;
+}
Deleted: grass-addons/imagery/i.landsat.toar/main.c
===================================================================
--- grass-addons/i.landsat.toar/main.c 2008-02-29 15:19:34 UTC (rev 30401)
+++ grass-addons/imagery/i.landsat.toar/main.c 2008-03-07 17:05:37 UTC (rev 30494)
@@ -1,436 +0,0 @@
-
-/****************************************************************************
- *
- * MODULE: i.landsat.toar
- *
- * AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
- *
- * PURPOSE: Calculate TOA Radiance or Reflectance and Kinetic Temperature
- * for Landsat 1/2/3/4/5 MS, 4/5 TM or 7 ETM+
- *
- * COPYRIGHT: (C) 2002, 2005 2008 by the GRASS Development Team
- *
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
- *
- *****************************************************************************/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-#include "local_proto.h"
-
-
-int main(int argc, char *argv[])
-{
- struct History history;
- struct GModule *module;
-
- struct Cell_head cellhd, window;
- char *mapset;
-
- void *inrast, *outrast;
- int infd, outfd;
- void *ptr;
- int nrows, ncols, row, col;
-
- 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;
-
- char band_in[127], band_out[127];
- int i, method, qcal;
- double rad, ref, percent;
- lsat_data lsat;
- char command[300];
-
- /* initialize GIS environment */
- G_gisinit(argv[0]);
-
- /* initialize module */
- module = G_define_module();
- module->description = _("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+");
- module->keywords = _("Top-Of-Atmosphere Radiance or Reflectance.\n Landsat-1 MSS, Landsat-2 MSS, Landsat-3 MSS, Landsat-4 MSS, Landsat-5 MSS,\n Landsat-4 TM, Landsat-5 TM,\n Landsat-7 ETM+");
-
- /* 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->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->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->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->description = _("Solar elevation in degrees");
-
- bgain = G_define_option();
- 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->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";
-
- 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";
-
- /* It defines the different flags */
- frad = G_define_flag();
- frad->key = 'r';
- frad->description = _("Output at-sensor radiance for all bands");
- frad->answer = 0;
-
- sat1 = G_define_flag();
- sat1->key = '1';
- sat1->description = _("Landsat-1 MSS");
- sat1->answer = 0;
-
- sat2 = G_define_flag();
- sat2->key = '2';
- sat2->description = _("Landsat-2 MSS");
- sat2->answer = 0;
-
- sat3 = G_define_flag();
- sat3->key = '3';
- sat3->description = _("Landsat-3 MSS");
- sat3->answer = 0;
-
- sat4 = G_define_flag();
- sat4->key = '4';
- sat4->description = _("Landsat-4 TM");
- sat4->answer = 0;
-
- sat5 = G_define_flag();
- sat5->key = '5';
- sat5->description = _("Landsat-5 TM");
- sat5->answer = 0;
-
- sat7 = G_define_flag();
- sat7->key = '7';
- sat7->description = _("Landsat-7 ETM+");
- sat7->answer = 0;
-
- msss = G_define_flag();
- msss->key = 's';
- msss->description = _("Modify sensor of Landsat-4/5 to MSS");
- msss->answer = 0;
-
- verbo = G_define_flag();
- verbo->key = 'v';
- verbo->description = _("Show parameters applied");
- verbo->answer = 0;
-
- /* options and afters parser */
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
-
- /*****************************************
- * ---------- START --------------------
- * Stores options and flag to variables
- *****************************************/
- met = metfn->answer;
- name = input->answer;
- atcor = metho->answer;
-
- if (adate->answer != NULL) {
- 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';
-
- lsat.sun_elev = elev->answer == NULL ? 0. : atof(elev->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"));
-
- /* 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"));
- }
- 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"));
- }
- }
- 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);
- }
-
-
- /*****************************************
- * ---------- 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);
- }
- }
- fprintf(stdout, "-------------------\n");
- fflush(stdout);
- }
-
- /*****************************************
- * ---------- CALCULUS -----------------
- *****************************************/
- G_get_window(&window);
-
- 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);
-
- mapset = G_find_cell2(band_in, "");
- if (mapset == NULL) {
- G_warning(_("raster file [%s] not found"), band_in);
- continue; }
-
- if (G_legal_filename(band_out) < 0)
- G_fatal_error(_("[%s] is an illegal name"), band_out);
-
- 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);
-
- /* 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 ((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);
-
- nrows = G_window_rows();
- ncols = G_window_cols();
- for (row = 0; row < nrows; row++)
- {
- if (verbose)
- G_percent(row, nrows, 2);
-
- if (G_get_raster_row(infd, inrast, row, in_data_type) < 0)
- G_fatal_error(_("Could not read from <%s>"), band_in);
-
- 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;
- }
- }
-
- 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);
-
- 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_set_window(&window);
- exit(EXIT_SUCCESS);
-}
Copied: grass-addons/imagery/i.landsat.toar/main.c (from rev 30464, grass-addons/i.landsat.toar/main.c)
===================================================================
--- grass-addons/imagery/i.landsat.toar/main.c (rev 0)
+++ grass-addons/imagery/i.landsat.toar/main.c 2008-03-07 17:05:37 UTC (rev 30494)
@@ -0,0 +1,577 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.landsat.toar
+ *
+ * AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
+ *
+ * PURPOSE: Calculate TOA Radiance or Reflectance and Kinetic Temperature
+ * for Landsat 1/2/3/4/5 MS, 4/5 TM or 7 ETM+
+ *
+ * COPYRIGHT: (C) 2002, 2005 2008 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+int main(int argc, char *argv[])
+{
+ struct History history;
+ struct GModule *module;
+
+ struct Cell_head cellhd, window;
+ char *mapset;
+
+ void *inrast, *outrast;
+ int infd, outfd;
+ void *ptr;
+ int nrows, ncols, row, col;
+
+ RASTER_MAP_TYPE in_data_type;
+
+ 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;
+
+ 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]);
+
+ /* initialize module */
+ module = G_define_module();
+ module->description = _("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+");
+ module->keywords = _("Top-Of-Atmosphere Radiance or Reflectance.\n Landsat-1 MSS, Landsat-2 MSS, Landsat-3 MSS, Landsat-4 MSS, Landsat-5 MSS,\n Landsat-4 TM, Landsat-5 TM,\n Landsat-7 ETM+");
+
+ /* 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->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->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->description = _("Image acquisition date (yyyy-mm-dd)");
+
+ elev = G_define_option();
+ 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->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->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,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 radiance in path radiance");
+ perc->description = _("Percent of solar radiance in path radiance");
+ perc->answer = "0.01";
+
+ 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->description = _("Output at-sensor radiance for all bands");
+ frad->answer = 0;
+
+ sat1 = G_define_flag();
+ sat1->key = '1';
+ sat1->description = _("Landsat-1 MSS");
+ sat1->answer = 0;
+
+ sat2 = G_define_flag();
+ sat2->key = '2';
+ sat2->description = _("Landsat-2 MSS");
+ sat2->answer = 0;
+
+ sat3 = G_define_flag();
+ sat3->key = '3';
+ sat3->description = _("Landsat-3 MSS");
+ sat3->answer = 0;
+
+ sat4 = G_define_flag();
+ sat4->key = '4';
+ sat4->description = _("Landsat-4 TM");
+ sat4->answer = 0;
+
+ sat5 = G_define_flag();
+ sat5->key = '5';
+ sat5->description = _("Landsat-5 TM");
+ sat5->answer = 0;
+
+ sat7 = G_define_flag();
+ sat7->key = '7';
+ sat7->description = _("Landsat-7 ETM+");
+ sat7->answer = 0;
+
+ msss = G_define_flag();
+ msss->key = 's';
+ msss->description = _("Set sensor of Landsat-4/5 to MSS");
+ msss->answer = 0;
+
+ verbo = G_define_flag();
+ verbo->key = 'v';
+ verbo->description = _("Show parameters applied");
+ verbo->answer = 0;
+
+ /* options and afters parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /*****************************************
+ * ---------- START --------------------
+ * Stores options and flag to variables
+ *****************************************/
+ met = metfn->answer;
+ name = input->answer;
+
+ if (adate->answer != NULL) {
+ G_strncpy(lsat.date, adate->answer, 11);
+ lsat.date[10] = '\0';
+ }
+ else
+ lsat.date[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);
+ percent = atof(perc->answer);
+ pixel = atoi(dark->answer);
+ sat_zenith = atof(satz->answer);
+ rayleigh = atof(atmo->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"));
+
+ /* 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);
+ 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 != 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"));
+ }
+ }
+ }
+
+ /*****************************************
+ * ------------ 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;
+
+ 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);
+ }
+
+ /*****************************************
+ * ------------ 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);
+ }
+
+ /*****************************************
+ * ------------ CALCULUS -----------------
+ *****************************************/
+
+ 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);
+
+ 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);
+
+ 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"));
+
+ /* controlling, if we can write the raster */
+ if (G_legal_filename(band_out) < 0)
+ G_fatal_error(_("[%s] is an illegal name"), band_out);
+
+ if ((outfd = G_open_raster_new(band_out, DCELL_TYPE)) < 0)
+ G_fatal_error(_("Could not open <%s>"), 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();
+ /* ================================================================= */
+ G_message("%s of %s to %s",
+ (frad->answer ? "Radiance"
+ : (lsat.band[i].thermal) ? "Temperature" : "Reflectance"),
+ band_in, band_out);
+
+ for (row = 0; row < nrows; row++)
+ {
+ if (verbose)
+ G_percent(row, nrows, 2);
+
+ 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]);
+ }
+
+ /* ================================================================= */
+
+ G_free(inrast);
+ G_close_cell(infd);
+ G_free(outrast);
+ G_close_cell(outfd);
+
+ char command[300];
+ sprintf(command, "r.colors map=%s color=grey", band_out);
+ system(command);
+
+ G_short_history(band_out, "raster", &history);
+
+ 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