[GRASS-SVN] r29915 - grass-addons/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jan 31 08:23:36 EST 2008
Author: neteler
Date: 2008-01-31 08:23:35 -0500 (Thu, 31 Jan 2008)
New Revision: 29915
Added:
grass-addons/i.landsat.toar/landsat_met.c
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/local_proto.h
grass-addons/i.landsat.toar/main.c
Log:
E. Jorge Tizado: support for all Landsats and simple atmosferic correction
Modified: grass-addons/i.landsat.toar/description.html
===================================================================
--- grass-addons/i.landsat.toar/description.html 2008-01-31 11:15:48 UTC (rev 29914)
+++ grass-addons/i.landsat.toar/description.html 2008-01-31 13:23:35 UTC (rev 29915)
@@ -1,71 +1,112 @@
<H2>DESCRIPTION</H2>
-<EM>i.landsat.toar</EM> is to transform digital number (DN) of
-Landsat-4/5 TM or Landsat 7 ETM+ product to reflectance at top of
-atmosphere (band 1, 2, 3, 4, 5, 7 or 8) and temperature (band 6).
+<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>For Landsat-4/5 TM need date and solar.<br>For Landsat-7 ETM+
-either metfile or date, solar, and gain</p>
+<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>For flag -f, make distinction between acquisition and processing dates.</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>
-<H2>NOTES</H2>
+<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>
-<p>The standard geometric and radiometric correction results in
-digital number (DN) images. To further standardize the impact of
-illumination geometry, the DN images are converted first to
-at-satellite radiance and then to at-satellite reflectance.<br>
-The thermal band is first converted from DN to at-satellite radiance,
-and then to effective at-satellite temperature in Kelvin degrees.</p>
-<p>In verbose mode, the program write the parameters used in
-transformation to at-satellite radiance, to at-satellite reflectance
-and constant to remove the atmospheric effect using DOS
-method.<br><br>
+<H2>Uncorrected at-sensor values (method=uncorrected, default)</H2>
-To remove atmospheric effect, use
+<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>
-<div class="code"><pre>
-r.mapcalc 'refc = if(ref < (a*dp+b), b, ref-a*dp)'
-</pre></div>
-being <em>refc</em> the reflectance corrected without atmospheric
-effect, <em>ref</em> the at-satellite reflectance, and <em>dp</em> the
-dark pixel in DN scale.</p>
-
-
-<p>
-Other equations are:<br>
+<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 * DN + bias </li>
- <li> cte_ref = (PI * d_earth_sun * d_earth_sun) / (ESUN * cos(90 - sun_elevation) </li>
- <li> reflectance = cte_ref * radiance </li>
- <li> a = cte_ref * gain</li>
- <li> b = cte_ref * bias</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>More details on how to convert DN to at-satellite reflectance are
-provided by Markham and Barker (1986), Irish (2000), and Huang et
-al. (2002).</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 from DOS2 method (Chavez, 1996): TAUz = sin(e) [for the bands with maximum wave length less than 1., i.e. bands 4/5/6 MSS, 1-4 TM, and 1-4 ETM+] otherwise TAUz = 1.07, TAUv = 1.0, and Esky = 0.</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>Transfor digital numbers from band raster 203_30.1, 203_30.2,
-... to reflectance at top-of-atmosphere in output files 203_30.toar.1,
-203_30.toar.2, ...</p>
+<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><a href="http://landcover.usgs.gov/pdf/image_preprocessing.pdf">http://landcover.usgs.gov/pdf/image_preprocessing.pdf</a></li>
- <li>Gyanesh Chander (SAIC/EDC/USGS) and Brian Markham
- (LPSO/GSFC/NASA): Revised L5 TM Radiometric Calibration
- Procedures and Postcalibration Dynamic Ranges </li>
+ <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>
@@ -80,7 +121,8 @@
<H2>AUTHOR</H2>
-E. Jorge Tizado, University of León, Spain (ejtizado ono com)<BR>
+E. Jorge Tizado (ej.tizado unileon es)<br>
+Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
<p>
-<i>Last changed: $Date: 2006/11/17 00:00:00 $</i>
+<i>Last changed: $Date: 2007/07/07 00:00:00 $</i>
Modified: grass-addons/i.landsat.toar/landsat.c
===================================================================
--- grass-addons/i.landsat.toar/landsat.c 2008-01-31 11:15:48 UTC (rev 29914)
+++ grass-addons/i.landsat.toar/landsat.c 2008-01-31 13:23:35 UTC (rev 29915)
@@ -2,71 +2,95 @@
#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 Radiance
+ * PURPOSE: Calibrated Digital Number to at-satellite Radiance
*****************************************************************************/
-double lsat_qcal2rad(int qcal, double lmax, double lmin, double qcalmax, double qcalmin)
+double lsat_qcal2rad(int qcal, band_data * band)
{
- double gain, bias;
-
- gain = ((lmax - lmin) / (qcalmax - qcalmin));
- bias = lmin - gain * qcalmin;
-
- return (double)(((double)qcal) * gain + bias);
+ return (double)(((double)qcal) * band->gain + band->bias);
}
/****************************************************************************
- * PURPOSE: Radiance of non-thermal band to TOA Reflectance
+ * PURPOSE: Radiance of non-thermal band to at-satellite Reflectance
*****************************************************************************/
-double lsat_rad2ref(double rad, double dist_es, double sun_elev, double esun)
+double lsat_rad2ref(double rad, band_data * band)
{
- double a, b;
-
- a = PI * dist_es * dist_es;
- b = esun * cos((90 - sun_elev) * D2R);
-
- return (double)((rad * a) / b);
+ return (double)(rad / band->K2);
}
-double lsat_refrad_ratio(double dist_es, double sun_elev, double esun)
+/****************************************************************************
+ * PURPOSE: Radiance of thermal band to at-satellite Temperature
+ *****************************************************************************/
+double lsat_rad2temp(double rad, band_data * band)
{
- double a, b;
-
- a = PI * dist_es * dist_es;
- b = esun * cos((90 - sun_elev) * D2R);
-
- return (double)(a / b);
+ return (double)(band->K2 / log((band->K1 / rad) + 1.0));
}
/****************************************************************************
- * PURPOSE: Radiance of thermal band to Kinetic Temperature
+ * 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
*****************************************************************************/
-double lsat_rad2temp(double rad, double K1, double K2)
-{
- return (double)(K2 / log((K1 / rad) + 1.0));
-}
-/****************************************************************************
- * PUPOSE: Calibrated digital number of Landsat 5 to TOA Reflectance
- *****************************************************************************
-double lsat_qcal2ref_tm5(int qcal, int band)
+void lsat_bandctes(lsat_data * lsat, int i, char method, double aot, double percent)
{
- static const double slope[] =
- { 0.9398, 1.7731, 1.5348, 1.4239, 0.9828, 0., 1.3017 };
- static const double inter[] =
- { 4.2934, 4.7289, 3.9796, 7.0320, 7.0185, 0., 7.6568 };
+ double a, b, rad_sun;
+ double TAUv, TAUz, Edown;
- static const double gain[] =
- { 0.7756863, 0.7956862, 0.6192157, 0.6372549, 0.1257255, 0., 0.0437255 };
- static const double bias[] =
- { -6.1999969, -6.3999939, -5.0000000, -5.1000061, -0.9999981, 0., -0.3500004 };
+ a = (double)(PI * lsat->dist_es * lsat->dist_es);
+ b = (double)(sin(D2R * lsat->sun_elev));
- --band;
- return (gain[band] * ((double)(qcal) * slope[band] + inter[band]) +
- bias[band]);
+ /** 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;
+ }
}
-*/
+
Modified: grass-addons/i.landsat.toar/landsat.h
===================================================================
--- grass-addons/i.landsat.toar/landsat.h 2008-01-31 11:15:48 UTC (rev 29914)
+++ grass-addons/i.landsat.toar/landsat.h 2008-01-31 13:23:35 UTC (rev 29915)
@@ -3,24 +3,46 @@
#define MAX_BANDS 9
-/* Band data */
+#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 */
+ int code; /* Band code */
+
+ double wavemax, wavemin; /* Wavelength in µm */
+
double lmax, lmin; /* Spectral radiance */
double qcalmax, qcalmin; /* Quantized calibrated pixel */
- double esun; /* Solar espectral irradiance */
+ 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;
-/* Landsat data */
typedef struct
{
- char date[11]; /* Satelite image date */
+ 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 */
- double K1, K2; /* Thermal calibration constant */
+
+ char sensor[5]; /* Type of sensor: MSS, TM, ETM+ */
int bands; /* Total number of bands */
band_data band[MAX_BANDS]; /* Data for each band */
@@ -31,9 +53,10 @@
* Landsat Equations Prototypes
*****************************************************************************/
-double lsat_qcal2rad(int qcal, double lmax, double lmin, double qcalmax, double qcalmin);
-double lsat_rad2ref(double rad, double dist_es, double sun_elev, double esun);
-double lsat_refrad_ratio(double dist_es, double sun_elev, double esun);
-double lsat_rad2temp(double rad, double K1, double K2);
+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
Added: grass-addons/i.landsat.toar/landsat_met.c
===================================================================
--- grass-addons/i.landsat.toar/landsat_met.c (rev 0)
+++ grass-addons/i.landsat.toar/landsat_met.c 2008-01-31 13:23:35 UTC (rev 29915)
@@ -0,0 +1,177 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include <math.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+#include "earth_sun.h"
+
+#define ETM_MET_SIZE 5600 /* .met file size 5516 bytes */
+#define TM5_MET_SIZE 28700 /* .met file size 28686 bytes */
+#define MAX_STR 127
+
+inline void chrncpy(char * dest, char * src, int n)
+{
+ if (src == NULL) n = 1;
+ else strncpy(dest,src,n);
+ dest[n-1] = '\0';
+}
+
+/****************************************************************************
+ * PURPOSE: Read values of Landsat-7 ETM+ from header (.met) file
+ *****************************************************************************/
+void get_value_met7(const char mettext[], char *text, char value[])
+{
+ char *ptr;
+ value[0] = 0;
+
+ ptr = strstr(mettext, text);
+ if (ptr == NULL) return;
+
+ while (*ptr++ != '=') ;
+ sscanf(ptr, "%s", value);
+
+ return;
+}
+
+void met_ETM(char *metfile, lsat_data * lsat)
+{
+ FILE *f;
+ char mettext[ETM_MET_SIZE];
+ char name[MAX_STR], value[MAX_STR];
+ int i, j;
+
+ static int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
+ static int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
+
+ static double esun[] = { 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
+
+ if ((f = fopen(metfile, "r")) == NULL) {
+ G_fatal_error(_(".met file [%s] not found"), metfile);
+ }
+ fread(mettext, 1, ETM_MET_SIZE, f);
+
+ /* --------------------------------------- */
+ lsat->number = 7;
+
+ get_value_met7(mettext, "SENSOR_ID", value);
+ chrncpy(lsat->sensor, value+1, 5);
+
+ if (lsat->creation[0] == 0) {
+ get_value_met7(mettext, "CREATION_TIME", value);
+ chrncpy(lsat->creation, value, 11);
+ }
+
+ get_value_met7(mettext, "ACQUISITION_DATE", value);
+ chrncpy(lsat->date, value, 11);
+ lsat->dist_es = earth_sun(lsat->date);
+
+ get_value_met7(mettext, "SUN_ELEVATION", value);
+ lsat->sun_elev = atof(value);
+
+ 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].esun = *(esun + lsat->band[i].number - 1);
+ snprintf(name, MAX_STR, "LMAX_BAND%d", lsat->band[i].code);
+ get_value_met7(mettext, name, value);
+ lsat->band[i].lmax = atof(value);
+ snprintf(name, MAX_STR, "LMIN_BAND%d", lsat->band[i].code);
+ get_value_met7(mettext, name, value);
+ lsat->band[i].lmin = atof(value);
+ snprintf(name, MAX_STR, "QCALMAX_BAND%d", lsat->band[i].code);
+ get_value_met7(mettext, name, value);
+ lsat->band[i].qcalmax = atof(value);
+ snprintf(name, MAX_STR, "QCALMIN_BAND%d", lsat->band[i].code);
+ get_value_met7(mettext, name, value);
+ lsat->band[i].qcalmin = atof(value);
+ if (lsat->band[i].number == 6) {
+ lsat->band[i].thermal = 1;
+ lsat->band[i].K1 = 666.09;
+ lsat->band[i].K2 = 1282.71;
+ }
+ else {
+ lsat->band[i].thermal = 0;
+ }
+ }
+
+ (void)fclose(f);
+ return;
+}
+
+/****************************************************************************
+ * PURPOSE: Read values of Landsat MSS/TM from header (.met) file
+ *****************************************************************************/
+void get_value_met(const char mettext[], char *text, char value[])
+{
+ char *ptr;
+ value[0] = 0;
+
+ ptr = strstr(mettext, text);
+ if (ptr == NULL) return;
+
+ ptr = strstr(ptr, " VALUE ");
+ if (ptr == NULL) return;
+
+ while (*ptr++ != '\"') ;
+ sscanf(ptr, "%s", value);
+ ptr = value;
+ do { if (*ptr == '\"') *ptr = '\0'; } while (*ptr++);
+
+ return;
+}
+
+void met_TM5(char *metfile, lsat_data * lsat)
+{
+ FILE *f;
+ char mettext[TM5_MET_SIZE];
+ char metdate[MAX_STR], value[MAX_STR];
+
+ if ((f = fopen(metfile, "r")) == NULL) {
+ G_fatal_error(_(".met file [%s] not found"), metfile);
+ }
+ fread(mettext, 1, TM5_MET_SIZE, f);
+
+ /* --------------------------------------- */
+ get_value_met(mettext, "CALENDARDATE", value);
+ chrncpy(lsat->date, value, 11);
+
+ if (lsat->creation[0] == 0) {
+ get_value_met(mettext, "PRODUCTIONDATETIME", value);
+ chrncpy(lsat->creation, value, 11);
+ }
+
+ if (lsat->creation[0] == 0 )
+ G_fatal_error(_("Product creation date not in .met file [%s]"), metfile);
+
+ get_value_met(mettext, "SolarElevation", value);
+ lsat->sun_elev = atof(value);
+
+ get_value_met(mettext, "PLATFORMSHORTNAME", value);
+ switch(value[8])
+ {
+ case '1':
+ set_MSS1(lsat);
+ break;
+ case '2':
+ set_MSS2(lsat);
+ break;
+ case '3':
+ set_MSS3(lsat);
+ break;
+ case '4':
+ get_value_met(mettext, "SENSORSHORTNAME", value);
+ if (value[0] = 'M') set_MSS4(lsat); else set_TM4(lsat);
+ break;
+ case '5':
+ get_value_met(mettext, "SENSORSHORTNAME", value);
+ if (value[0] = 'M') set_MSS5(lsat); else set_TM5(lsat);
+ break;
+ }
+
+ (void)fclose(f);
+ return;
+}
Property changes on: grass-addons/i.landsat.toar/landsat_met.c
___________________________________________________________________
Name: svn:eol-style
+ native
Modified: grass-addons/i.landsat.toar/landsat_set.c
===================================================================
--- grass-addons/i.landsat.toar/landsat_set.c 2008-01-31 11:15:48 UTC (rev 29914)
+++ grass-addons/i.landsat.toar/landsat_set.c 2008-01-31 13:23:35 UTC (rev 29915)
@@ -6,170 +6,375 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-#include "landsat.h"
+#include "local_proto.h"
#include "earth_sun.h"
-#define MET_SIZE 5600 /* .met file size 5516 bytes */
-#define MAX_STR 127
+void sensor_MSS(lsat_data * lsat)
+{
+ int i;
-/* utility for read met file */
-void get_value_met(char *mettext, char *text, char *value)
+ int band[] = { 4, 5, 6, 7 };
+ int code[] = { 1, 2, 3, 4 };
+ 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)
{
- char *ptr;
+ int i;
- ptr = strstr(mettext, text);
- if (ptr == NULL)
- return;
+ 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 };
- while (*ptr++ != '=') ;
- sscanf(ptr, "%s", value);
+ strcpy(lsat->sensor, "TM");
+
+ lsat->bands = 7;
+ for (i = 0; i < lsat->bands; i++) {
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = lsat->band[i].number;
+ 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: Read values of Landsat-7 ETM+ from header (.met) file
+ * PURPOSE: Store values of Landsat-1 MSS
+ * July 23, 1972 to January 6, 1978
*****************************************************************************/
-void met_ETM(char *metfile, lsat_data * lsat)
+void set_MSS1(lsat_data * lsat)
{
- FILE *f;
- char mettext[MET_SIZE];
- char name[MAX_STR], value[MAX_STR];
int i, j;
- static int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
- static int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
+ /** 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. };
- static double esun[] =
- { 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
+ lsat->number = 1;
+ sensor_MSS( lsat );
+ lsat->dist_es = earth_sun(lsat->date);
- if ((f = fopen(metfile, "r")) == NULL) {
- G_fatal_error(_(".met file '%s' not found"), metfile);
+ 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);
}
- fread(mettext, 1, MET_SIZE, f);
+ return;
+}
- get_value_met(mettext, "ACQUISITION_DATE", value);
- strncpy(lsat->date, value, 11);
- lsat->dist_es = earth_sun(lsat->date);
+/****************************************************************************
+ * 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;
- get_value_met(mettext, "SUN_ELEVATION", value);
- lsat->sun_elev = atof(value);
+ /** 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. };
- lsat->K1 = 666.09;
- lsat->K2 = 1282.71;
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("1975-07-16")) i = 0;
+ else i = 1;
+ lmax = Lmax[i];
+ lmin = Lmin[i];
- lsat->bands = 9;
+ lsat->number = 2;
+ sensor_MSS( lsat );
+
+ lsat->dist_es = earth_sun(lsat->date);
+
for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(code + i);
- lsat->band[i].esun = *(esun + lsat->band[i].number - 1);
-
- snprintf(name, MAX_STR, "LMAX_BAND%d", lsat->band[i].code);
- get_value_met(mettext, name, value);
- lsat->band[i].lmax = atof(value);
- snprintf(name, MAX_STR, "LMIN_BAND%d", lsat->band[i].code);
- get_value_met(mettext, name, value);
- lsat->band[i].lmin = atof(value);
- snprintf(name, MAX_STR, "QCALMAX_BAND%d", lsat->band[i].code);
- get_value_met(mettext, name, value);
- lsat->band[i].qcalmax = atof(value);
- snprintf(name, MAX_STR, "QCALMIN_BAND%d", lsat->band[i].code);
- get_value_met(mettext, name, value);
- lsat->band[i].qcalmin = atof(value);
+ j = lsat->band[i].number - 1;
+ lsat->band[i].esun = *(esun + j);
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
}
-
- (void)fclose(f);
return;
}
-
/****************************************************************************
- * PURPOSE: Store values of Landsat-4 TM
+ * PURPOSE: Store values of Landsat-3 MSS
+ * March 5, 1978 to March 31, 1983
*
- * date: adquisition date of image
- * elevation: solar elevation angle
+ * tiene una banda 8 thermal
*****************************************************************************/
-void set_TM4(lsat_data * lsat, char date[], double elevation)
+void set_MSS3(lsat_data * lsat)
{
- int i;
+ int i, j;
+ double julian, *lmax, *lmin;
- static int band[] = { 1, 2, 3, 4, 5, 6, 7 };
- static double esun[] = { 1958., 1828., 1559., 1045., 219.1, 0., 74.57 };
- static double lmax[] = { 185., 342., 245., 270., 36., 0., 19. };
- static double lmin[] = { -1.5, -3.1, -2.7, -2.5, -0.45, 0., -0.3 };
+ /** 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. };
- strncpy(lsat->date, date, 11);
- lsat->sun_elev = elevation;
+ 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);
- lsat->K1 = 607.76;
- lsat->K2 = 1260.56;
+ 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;
+}
- lsat->bands = 7;
+/****************************************************************************
+ * 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++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = lsat->band[i].number;
- lsat->band[i].esun = *(esun + lsat->band[i].number - 1);
- lsat->band[i].lmax = *(lmax + lsat->band[i].number - 1);
- lsat->band[i].lmin = *(lmin + lsat->band[i].number - 1);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 0.;
+ j = lsat->band[i].number - 1;
+ lsat->band[i].esun = *(esun + j);
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
}
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 TM
- *
- * date: adquisition date of image
- * elevation: solar elevation angle
- * after: flag for processing date
+ * PURPOSE: Store values of Landsat-5 MSS/TM
+ * March 1, 1984 to today
*****************************************************************************/
-void set_TM5(lsat_data * lsat, char date[], double elevation, char after)
+void set_MSS5(lsat_data * lsat)
{
- int i;
- double *lmax, *lmin;
+ int i, j;
+ double julian, *lmax, *lmin;
- static int band[] = { 1, 2, 3, 4, 5, 6, 7 };
- static double esun[] = { 1957., 1829., 1557., 1047., 219.3, 0., 74.52 };
+ /** 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. };
- /* Spectral radiances at detecter before May 4, 2003 */
- static double lmaxb[] =
- { 152.10, 296.81, 204.30, 206.20, 27.19, 15.303, 14.38 };
- static double lminb[] =
- { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 };
+ 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];
- /* Spectral radiances at detecter after May 5, 2003 */
- static double lmaxa[] =
- { 193.00, 365.00, 264.00, 221.00, 30.20, 15.303, 16.50 };
- static double lmina[] =
- { -1.52, -2.84, -1.17, -1.51, -0.37, 1.2378, -0.15 };
+ lsat->number = 5;
+ sensor_MSS( lsat );
- strncpy(lsat->date, date, 11);
- lsat->sun_elev = elevation;
lsat->dist_es = earth_sun(lsat->date);
- lsat->K1 = 607.76;
- lsat->K2 = 1260.56;
+ 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;
+}
- lsat->bands = 7;
- if (after != 0) {
- lmax = lmaxa;
- lmin = lmina;
+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 }, /* on or after May 4, 2003 */
+ { 169.00, 333.00, 264.00, 221.00, 30.20, 15.303, 16.50 } }; /* on or 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;
+ }
}
- else {
- lmax = lmaxb;
- lmin = lminb;
- }
+
+ lsat->number = 5;
+ sensor_TM( lsat );
+
+ lsat->dist_es = earth_sun(lsat->date);
+
for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = lsat->band[i].number;
- lsat->band[i].esun = *(esun + lsat->band[i].number - 1);
- lsat->band[i].lmax = *(lmax + lsat->band[i].number - 1);
- lsat->band[i].lmin = *(lmin + lsat->band[i].number - 1);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 0.;
+ j = lsat->band[i].number - 1;
+ lsat->band[i].esun = *(esun + j);
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
+ if (lsat->band[i].thermal ) {
+ lsat->band[i].K1 = 607.76;
+ lsat->band[i].K2 = 1260.56;
+ }
}
return;
}
@@ -177,81 +382,56 @@
/****************************************************************************
* PURPOSE: Store values of Landsat-7 ETM+
- *
- * date: adquisition date of image
- * elevation: solar elevation angle
- * gain: nine H/L chars for band gain
- * before: flag for processing date
+ * April 15, 1999 to today
*****************************************************************************/
-void set_ETM(lsat_data * lsat, char date[], double elevation, char gain[],
- char before)
+void set_ETM(lsat_data * lsat, char gain[])
{
- int i;
- double *lmax, *lmin;
+ int i, k, j;
+ double julian, *lmax, *lmin;
- static int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
- static int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
+ /** 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 */
- static double esun[] =
- { 1969., 1840., 1551., 1044., 225.7, 0., 82.07, 1368. };
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("2000-07-01")) k = 0;
+ else k = 1;
- /* Spectral radiances at detector before July 1, 2000 */
- /* Low gain */
- static double lmaxLb[] =
- { 297.5, 303.4, 235.5, 235.0, 47.70, 17.04, 16.600, 244.0 };
- static double lminLb[] = { -6.2, -6.0, -4.5, -4.5, -1.0, 0.0, -0.35, -5.0 };
- /* High gain */
- static double lmaxHb[] =
- { 194.3, 202.4, 158.6, 157.5, 31.76, 12.65, 10.932, 158.4 };
- static double lminHb[] = { -6.2, -6.0, -4.5, -4.5, -1.0, 3.2, -0.35, -5.0 };
+ lsat->number = 7;
+ sensor_ETM( lsat );
- /* Spectral radiances at detector after July 1, 2000 */
- /* Low gain */
- static double lmaxLa[] =
- { 293.7, 300.9, 234.4, 241.1, 47.57, 17.04, 16.540, 243.1 };
- static double lminLa[] = { -6.2, -6.4, -5.0, -5.1, -1.0, 0.0, -0.35, -4.7 };
- /* High gain */
- static double lmaxHa[] =
- { 191.6, 196.5, 152.9, 157.4, 31.06, 12.65, 10.800, 158.3 };
- static double lminHa[] = { -6.2, -6.4, -5.0, -5.1, -1.0, 3.2, -0.35, -4.7 };
-
-
- strncpy(lsat->date, date, 11);
- lsat->sun_elev = elevation;
lsat->dist_es = earth_sun(lsat->date);
- lsat->K1 = 666.09;
- lsat->K2 = 1282.71;
-
- 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].esun = *(esun + lsat->band[i].number - 1);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 1.;
- if (before != 0) {
- if (gain[i] == 'H' || gain[i] == 'h') {
- lmax = lmaxHb;
- lmin = lminHb;
- }
- else {
- lmax = lmaxLb;
- lmin = lminLb;
- }
- }
- else {
- if (gain[i] == 'H' || gain[i] == 'h') {
- lmax = lmaxHa;
- lmin = lminHa;
- }
- else {
- lmax = lmaxLa;
- lmin = lminLa;
- }
- }
- lsat->band[i].lmax = *(lmax + lsat->band[i].number - 1);
- lsat->band[i].lmin = *(lmin + lsat->band[i].number - 1);
+ 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/local_proto.h
===================================================================
--- grass-addons/i.landsat.toar/local_proto.h 2008-01-31 11:15:48 UTC (rev 29914)
+++ grass-addons/i.landsat.toar/local_proto.h 2008-01-31 13:23:35 UTC (rev 29915)
@@ -1,7 +1,20 @@
+#ifndef _LOCAL_PROTO_H
+#define _LOCAL_PROTO_H
+
+#include <string.h>
#include "landsat.h"
void met_ETM (char *, lsat_data *);
-void set_TM4 (lsat_data *, char, double);
-void set_TM5 (lsat_data *, char, double, char);
-void set_ETM (lsat_data *, char, double, char, char);
+void met_TM5 (char *, lsat_data *);
+void set_MSS1(lsat_data *);
+void set_MSS2(lsat_data *);
+void set_MSS3(lsat_data *);
+void set_MSS4(lsat_data *);
+void set_MSS5(lsat_data *);
+
+void set_TM4 (lsat_data *);
+void set_TM5 (lsat_data *);
+void set_ETM (lsat_data *, char []);
+
+#endif
Modified: grass-addons/i.landsat.toar/main.c
===================================================================
--- grass-addons/i.landsat.toar/main.c 2008-01-31 11:15:48 UTC (rev 29914)
+++ grass-addons/i.landsat.toar/main.c 2008-01-31 13:23:35 UTC (rev 29915)
@@ -1,14 +1,14 @@
/****************************************************************************
*
- * MODULE: r.landsat.toar
+ * MODULE: i.landsat.toar
*
- * AUTHOR(S): E. Jorge Tizado - ejtizado at unileon.es
+ * AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
*
- * PURPOSE: Calculate TOA Reflectance and Kinectic Temperature
- * for Landsat 4/5 TM or 7 ETM+
+ * 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 by the GRASS Development Team
+ * 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
@@ -18,41 +18,39 @@
#include <stdio.h>
#include <stdlib.h>
-#include <string.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include "local_proto.h"
-#include "landsat.h"
+
int main(int argc, char *argv[])
{
+ struct History history;
+ struct GModule *module;
+
struct Cell_head cellhd, window;
- char *mapset; /* mapset name */
+ char *mapset;
- void *inrast; /* input buffer */
- unsigned char *outrast; /* output buffer */
- int nrows, ncols;
- int row, col;
- int infd, outfd; /* file descriptor */
+ void *inrast, *outrast;
+ int infd, outfd;
+ void *ptr;
+ int nrows, ncols, row, col;
RASTER_MAP_TYPE in_data_type;
- RASTER_MAP_TYPE out_data_type = DCELL_TYPE;
- struct History history; /* holds meta-data (title, comments,..) */
- struct GModule *module; /* GRASS module for parsing arguments */
-
int verbose = 1;
- struct Option *input, *output, *metfn, *date, *elev, *bgain;
- char *name; /* band raster name */
- char *met; /* met filename */
- double elevation; /* solar elevation */
- struct Flag *sat4, *sat5, *sat7, *flag, *param;
+ 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, qcal;
- double gain, bias, cref, rad, ref;
-
+ int i, method, qcal;
+ double rad, ref, percent;
lsat_data lsat;
char command[300];
@@ -61,37 +59,31 @@
/* initialize module */
module = G_define_module();
- module->description = _("Calculates reflectance at top of atmosphere or temperature for Landsat");
+ 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 (.?)");
+ 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 = _("Header File [.met] (only ETM+)");
+ metfn->description = _("Landsat ETM+ or TM5 header file (.met)");
- bgain = G_define_option();
- bgain->key = _("gain");
- bgain->type = TYPE_STRING;
- bgain->required = NO;
- bgain->gisprompt = _("band gain");
- bgain->description = _("Gain of bands 1,2,3,4,5,61,62,7,8 (only ETM+)");
- bgain->answer = "HHHLHLHHL";
+ adate = G_define_option();
+ adate->key = _("date");
+ adate->type = TYPE_STRING;
+ adate->required = NO;
+ adate->gisprompt = _("image acquisition adate");
+ adate->description = _("Image acquisition adate (yyyy-mm-dd)");
- date = G_define_option();
- date->key = _("date");
- date->type = TYPE_STRING;
- date->required = NO;
- date->gisprompt = _("image acquisition date");
- date->description = _("Image acquisition date (yyyy-mm-dd)");
-
elev = G_define_option();
elev->key = _("solar");
elev->type = TYPE_DOUBLE;
@@ -99,107 +91,233 @@
elev->gisprompt = _("solar elevation");
elev->description = _("Solar elevation in degrees");
- /* Define the different flags */
+ 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 adate");
+ pdate->description = _("Image creation adate (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 (specify date and solar)");
+ sat4->description = _("Landsat-4 TM");
sat4->answer = 0;
sat5 = G_define_flag();
sat5->key = '5';
- sat5->description = _("Landsat-5 TM (specify date and solar)");
+ sat5->description = _("Landsat-5 TM");
sat5->answer = 0;
sat7 = G_define_flag();
sat7->key = '7';
- sat7->description = _("Landsat-7 ETM+ (specify either met or date, solar, and gain)");
+ sat7->description = _("Landsat-7 ETM+");
sat7->answer = 0;
- flag = G_define_flag();
- flag->key = 'f';
- flag->description = _("LANDSAT-5: product creation after instead of before May 5, 2003.\n\n\tLANDSAT-7: product creation before instead of after July 1, 2000");
- flag->answer = 0;
+ msss = G_define_flag();
+ msss->key = 's';
+ msss->description = _("Modify sensor of Landsat-4/5 to MSS");
+ msss->answer = 0;
- param = G_define_flag();
- param->key = 'v';
- param->description = _("Show parameters applied");
+ verbo = G_define_flag();
+ verbo->key = 'v';
+ verbo->description = _("Show parameters applied");
+ verbo->answer = 0;
- /* options and flags parser */
+ /* options and afters parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- /* ----- START -------------------- */
- /* stores options and flags to variables */
+ /*****************************************
+ * ---------- START --------------------
+ * Stores options and flag to variables
+ *****************************************/
met = metfn->answer;
name = input->answer;
- elevation = 0.;
+ atcor = metho->answer;
- if( !sat4->answer && !sat5->answer && !sat7->answer)
- G_fatal_error(_("Need satellite type"));
+ if (adate->answer != NULL) {
+ strncpy(lsat.date, adate->answer, 11);
+ lsat.date[10] = '\0';
+ }
+ else
+ lsat.date[0] = '\0';
- if (sat7->answer && met != NULL) {
- met_ETM(met, &lsat);
- fprintf(stdout, "Landsat-7 ETM+ with data in met file [%s]\n", met);
+ if (pdate->answer != NULL) {
+ strncpy(lsat.creation, pdate->answer, 11);
+ lsat.creation[10] = '\0';
}
- else if (date->answer == NULL || elev->answer == NULL) {
- G_fatal_error(_("Need date and solar elevation"));
+ 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 adate and solar elevation */
+ else if (adate->answer == NULL || elev->answer == NULL) {
+ G_fatal_error(_("Lacking adate and solar elevation for this satellite"));
+ }
else {
- elevation = atof(elev->answer);
- if (sat7->answer) {
- if (bgain->answer && strlen(bgain->answer) > 8) {
- set_ETM(&lsat, date->answer, elevation, bgain->answer, flag->answer);
- fprintf(stdout, "Landsat 7 ETM+ %s July 1, 2000\n", ((flag->answer) ? "before" : "after"));
+ 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(_("Need band gain with 9 (H/L) data for Landsat-7"));
+ G_fatal_error(_("For Landsat-7 is necessary band gain with 9 (H/L) data"));
}
}
- else {
- if (sat4->answer) {
- set_TM4(&lsat, date->answer, elevation);
- fprintf(stdout, "Landsat-4 TM\n");
- }
- else {
- set_TM5(&lsat, date->answer, elevation, flag->answer);
- fprintf(stdout, "Landsat-5 TM %s May 5, 2003\n", ((flag->answer) ? "after" : "before"));
- }
+ 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"));
+ }
}
}
- if (param->answer) {
- fprintf(stdout, " ACQUISITION_DATE %s\n", lsat.date);
+ /* Now calculate band constants */
+ for (i = 0; i < lsat.bands; i++) {
+ switch(atcor[0])
+ {
+ case 'c':
+ method = CORRECTED;
+ break;
+ case 's':
+ method = SIMPLIFIED;
+ break;
+ default:
+ method = UNCORRECTED;
+ break;
+ }
+ 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, " solar elevation angle = %.8lf\n", lsat.sun_elev);
+ 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 (code %d)\n", lsat.band[i].number, lsat.band[i].code);
- if (lsat.band[i].number == 6) {
- fprintf(stdout, " const K1 = %lf\n", lsat.K1);
- fprintf(stdout, " const K2 = %lf\n", lsat.K2);
+ 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, " calibrated digital number (QCAL): %.1lf to %.1lf\n", lsat.band[i].qcalmin,
- lsat.band[i].qcalmax);
- fprintf(stdout, " radiance at-detector (L): %.3lf to %.3lf ", lsat.band[i].lmin, lsat.band[i].lmax);
-
- gain = (lsat.band[i].lmax - lsat.band[i].lmin) / (lsat.band[i].qcalmax - lsat.band[i].qcalmin);
- bias = lsat.band[i].lmin - gain * lsat.band[i].qcalmin;
- fprintf(stdout, "... radiance = %.8lf DN + %.8lf\n", gain, bias);
-
- cref = lsat_refrad_ratio(lsat.dist_es, lsat.sun_elev, lsat.band[i].esun);
- fprintf(stdout, " exoatmospheric irradiance (ESUN): %.5lf ", lsat.band[i].esun);
- fprintf(stdout, "... reflectance = %.8lf · radiance\n", cref);
-
- fprintf(stdout, " atmospheric effect: a = %.8lf, b = %.8lf\n", cref * gain, cref * bias);
- fprintf(stdout, " r.mapcalc '%s.%drc=if(%s.%dr<(%.8lf*dp+%.8lf),%.8lf,%s.%dr-%.8lf*dp)'\n",
- name, lsat.band[i].number, name, lsat.band[i].number, cref*gain, cref*bias, cref*bias, name, lsat.band[i].number, cref*gain);
+ 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++) {
@@ -208,19 +326,16 @@
mapset = G_find_cell2(band_in, "");
if (mapset == NULL) {
- G_warning(_("cell file [%s] not found"), band_in);
- continue;
- }
+ 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);
- /* determine the inputmap type (CELL/FCELL/DCELL) */
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);
- /* controlling, if we can open input raster */
if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
G_fatal_error(_("Cannot read file header of [%s]"), band_in);
@@ -230,65 +345,84 @@
/* Allocate input and output buffer */
inrast = G_allocate_raster_buf(in_data_type);
- nrows = G_window_rows();
- ncols = G_window_cols();
- outrast = G_allocate_raster_buf(out_data_type);
+ outrast = G_allocate_raster_buf(DCELL_TYPE);
/* controlling, if we can write the raster */
- if ((outfd = G_open_raster_new(band_out, out_data_type)) < 0)
+ if ((outfd = G_open_raster_new(band_out, DCELL_TYPE)) < 0)
G_fatal_error(_("Could not open <%s>"), band_out);
/* ================================================================= */
- /* ----- CORE ----- */
- G_message("Reflectance of %s to %s", band_in, band_out);
- for (row = 0; row < nrows; row++) {
+ 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)
+ 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++) {
-
- qcal = (int)((CELL *) inrast)[col];
-
- if (qcal < lsat.band[i].qcalmin) {
- ref = -1.;
- }
- else {
- rad = lsat_qcal2rad(qcal, lsat.band[i].lmax, lsat.band[i].lmin, lsat.band[i].qcalmax, lsat.band[i].qcalmin);
- if (lsat.band[i].number == 6) {
- ref = lsat_rad2temp(rad, lsat.K1, lsat.K2);
- }
- else {
- ref = lsat_rad2ref(rad, lsat.dist_es, lsat.sun_elev, lsat.band[i].esun);
- }
- }
- ((DCELL *) outrast)[col] = ref;
+ 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, out_data_type) < 0)
+ if (G_put_raster_row(outfd, outrast, DCELL_TYPE) < 0)
G_fatal_error(_("Cannot write to <%s>"), band_out);
}
- /* ----- END CORE ----- */
/* ================================================================= */
G_free(inrast);
+ G_close_cell(infd);
G_free(outrast);
- G_close_cell(infd);
G_close_cell(outfd);
- /* set map color to grey */
sprintf(command, "r.colors map=%s color=grey", band_out);
system(command);
- /* set -1. to null */
- G_message("REMEMBER: -1 is a NULL value");
-/* sprintf(command, "r.null map=%s setnull=-1", band_out); */
-/* system(command); */
-
- /* add command line incantation to history file */
G_short_history(band_out, "raster", &history);
G_command_history(&history);
G_write_history(band_out, &history);
More information about the grass-commit
mailing list