[GRASS-SVN] r43434 - grass-addons/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Sep 9 04:06:28 EDT 2010
Author: hamish
Date: 2010-09-09 08:06:28 +0000 (Thu, 09 Sep 2010)
New Revision: 43434
Modified:
grass-addons/imagery/i.landsat.toar/earth_sun.c
grass-addons/imagery/i.landsat.toar/landsat.c
grass-addons/imagery/i.landsat.toar/landsat.h
grass-addons/imagery/i.landsat.toar/landsat_met.c
grass-addons/imagery/i.landsat.toar/landsat_set.c
grass-addons/imagery/i.landsat.toar/local_proto.h
grass-addons/imagery/i.landsat.toar/main.c
Log:
run tools/grass_indent.sh
Modified: grass-addons/imagery/i.landsat.toar/earth_sun.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/earth_sun.c 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/earth_sun.c 2010-09-09 08:06:28 UTC (rev 43434)
@@ -1065,8 +1065,8 @@
b = 2 - a + (a / 4);
}
- return ((int)(365.25 * (year + 4716)) + (int)(30.6001 * (month + 1)) + day +
- b - 1524.5);
+ return ((int)(365.25 * (year + 4716)) + (int)(30.6001 * (month + 1)) +
+ day + b - 1524.5);
}
/* Get Julian day form Gregorian string yyyy-mm-dd */
Modified: grass-addons/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat.c 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/landsat.c 2010-09-09 08:06:28 UTC (rev 43434)
@@ -49,117 +49,121 @@
#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 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));
+ /* 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);
+ 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;
- lsat->band[i].K1 = 0.;
- lsat->band[i].K2 = rad_sun;
+ 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));
+ 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);
+ 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 == CORRECTED)
- {
- /* L = G * (DN - Qmin) + Lmin - Lmin
- -> bias = - G * Qmin */
- lsat->band[i].bias = - (lsat->band[i].gain * lsat->band[i].qcalmin);
- /* Another possibility is cut when rad < 0 */
- }
- else if (method > DOS )
- {
- /* L = Lsat - Lpath =
- G * DNsat + B - (G * dark + B - p * rad_sun) =
- G * DNsat - G * dark + p * rad_sun
- -> bias = p * rad_sun - G * dark */
- lsat->band[i].bias = percent * rad_sun - lsat->band[i].gain * dos;
- }
+ else if (method > DOS) {
+ /* L = Lsat - Lpath =
+ G * DNsat + B - (G * dark + B - p * rad_sun) =
+ G * DNsat - G * dark + p * rad_sun
+ -> bias = p * rad_sun - G * dark */
+ lsat->band[i].bias = percent * rad_sun - lsat->band[i].gain * dos;
}
+ }
}
-
-
-
-
Modified: grass-addons/imagery/i.landsat.toar/landsat.h
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat.h 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/landsat.h 2010-09-09 08:06:28 UTC (rev 43434)
@@ -22,34 +22,34 @@
typedef struct
{
- int number; /* Band number */
- int code; /* Band code */
+ int number; /* Band number */
+ int code; /* Band code */
- double wavemax, wavemin; /* Wavelength in µm */
+ double wavemax, wavemin; /* Wavelength in µm */
- double lmax, lmin; /* Spectral radiance */
- double qcalmax, qcalmin; /* Quantized calibrated pixel */
- double esun; /* Mean solar irradiance */
+ double lmax, lmin; /* Spectral radiance */
+ double qcalmax, qcalmin; /* Quantized calibrated pixel */
+ double esun; /* Mean solar irradiance */
- char thermal; /* Flag to thermal band */
- double gain, bias; /* Gain and Bias of sensor */
- double K1, K2; /* Thermal calibration constants,
- or Rad2Ref constants */
+ char thermal; /* Flag to thermal band */
+ double gain, bias; /* Gain and Bias of sensor */
+ double K1, K2; /* Thermal calibration constants,
+ or Rad2Ref constants */
} band_data;
typedef struct
{
- unsigned char number; /* Landsat number */
+ unsigned char number; /* Landsat number */
- char creation[11]; /* Image production date */
- char date[11]; /* Image acquisition date */
- double dist_es; /* Distance Earth-Sun */
- double sun_elev; /* Solar elevation */
+ char creation[11]; /* Image production date */
+ char date[11]; /* Image acquisition date */
+ double dist_es; /* Distance Earth-Sun */
+ double sun_elev; /* Solar elevation */
- char sensor[5]; /* Type of sensor: MSS, TM, ETM+ */
- int bands; /* Total number of bands */
- band_data band[MAX_BANDS]; /* Data for each band */
+ char sensor[5]; /* Type of sensor: MSS, TM, ETM+ */
+ int bands; /* Total number of bands */
+ band_data band[MAX_BANDS]; /* Data for each band */
} lsat_data;
@@ -58,7 +58,7 @@
*****************************************************************************/
double lsat_qcal2rad(double, band_data *);
-double lsat_rad2ref (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);
Modified: grass-addons/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat_met.c 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/landsat_met.c 2010-09-09 08:06:28 UTC (rev 43434)
@@ -8,15 +8,17 @@
#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 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)
+inline void chrncpy(char *dest, char *src, int n)
{
- if (src == NULL) n = 1;
- else strncpy(dest,src,n);
- dest[n-1] = '\0';
+ if (src == NULL)
+ n = 1;
+ else
+ strncpy(dest, src, n);
+ dest[n - 1] = '\0';
}
/****************************************************************************
@@ -25,10 +27,12 @@
void get_value_met7(const char mettext[], char *text, char value[])
{
char *ptr;
+
value[0] = 0;
ptr = strstr(mettext, text);
- if (ptr == NULL) return;
+ if (ptr == NULL)
+ return;
while (*ptr++ != '=') ;
sscanf(ptr, "%s", value);
@@ -43,10 +47,11 @@
char name[MAX_STR], value[MAX_STR];
int i, j;
- static int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
+ 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. };
+ 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);
@@ -57,11 +62,11 @@
lsat->number = 7;
get_value_met7(mettext, "SENSOR_ID", value);
- chrncpy(lsat->sensor, value+1, 5);
+ 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, "CREATION_TIME", value);
+ chrncpy(lsat->creation, value, 11);
}
get_value_met7(mettext, "ACQUISITION_DATE", value);
@@ -88,14 +93,14 @@
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;
- }
+ 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);
@@ -109,17 +114,21 @@
{
char *ptr;
int i;
+
value[0] = 0;
ptr = strstr(mettext, text);
- if (ptr == NULL) return;
+ if (ptr == NULL)
+ return;
ptr = strstr(ptr, " VALUE ");
- if (ptr == NULL) return;
+ if (ptr == NULL)
+ return;
i = 0;
while (*ptr++ != '\"') ;
- while( *ptr != '\"' && i < MAX_STR) value[i++]=*ptr++;
+ while (*ptr != '\"' && i < MAX_STR)
+ value[i++] = *ptr++;
value[i] = '\0';
return;
@@ -132,7 +141,7 @@
char metdate[MAX_STR], value[MAX_STR];
if ((f = fopen(metfile, "r")) == NULL) {
- G_fatal_error(_(".met file [%s] not found"), metfile);
+ G_fatal_error(_(".met file [%s] not found"), metfile);
}
fread(mettext, 1, TM5_MET_SIZE, f);
@@ -141,36 +150,42 @@
chrncpy(lsat->date, value, 11);
if (lsat->creation[0] == 0) {
- get_value_met(mettext, "PRODUCTIONDATETIME", value);
- chrncpy(lsat->creation, value, 11);
+ 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);
+ 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;
+ 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);
Modified: grass-addons/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat_set.c 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/landsat_set.c 2010-09-09 08:06:28 UTC (rev 43434)
@@ -11,79 +11,79 @@
void sensor_MSS(lsat_data * lsat)
{
- int i;
+ 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 */
+ /* green, red, near infrared, near infrared */
+ int band[] = { 1, 2, 3, 4 };
+ int code[] = { 4, 5, 6, 7 };
+ double wmax[] = { 0.6, 0.7, 0.8, 1.1 };
+ double wmin[] = { 0.5, 0.6, 0.7, 0.8 };
+ /* 79-82, 79-82, 79-82, 79-82 */
- strcpy(lsat->sensor, "MSS");
+ strcpy(lsat->sensor, "MSS");
- lsat->bands = 4;
- for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(code + i);
- lsat->band[i].wavemax = *(wmax + i);
- lsat->band[i].wavemin = *(wmin + i);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 0.;
- lsat->band[i].thermal = 0;
- }
- return;
+ lsat->bands = 4;
+ for (i = 0; i < lsat->bands; i++) {
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = *(code + i);
+ lsat->band[i].wavemax = *(wmax + i);
+ lsat->band[i].wavemin = *(wmin + i);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 0.;
+ lsat->band[i].thermal = 0;
+ }
+ return;
}
void sensor_TM(lsat_data * lsat)
{
- int i;
+ int i;
- /* 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 */
+ /* blue, green red, near infrared, shortwave IR, thermal IR, shortwave IR */
+ int band[] = { 1, 2, 3, 4, 5, 6, 7 };
+ double wmax[] = { 0.52, 0.60, 0.69, 0.90, 1.75, 12.50, 2.35 };
+ double wmin[] = { 0.45, 0.52, 0.63, 0.76, 1.55, 10.40, 2.08 };
+ /* 30, 30, 30, 30, 30, 120, 30 */
- strcpy(lsat->sensor, "TM");
+ strcpy(lsat->sensor, "TM");
- lsat->bands = 7;
- for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(band + i);
- lsat->band[i].wavemax = *(wmax + i);
- lsat->band[i].wavemin = *(wmin + i);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 0.; /* Modified in set_TM5 by date */
- lsat->band[i].thermal = (lsat->band[i].number == 6 ? 1 : 0);
- }
- return;
+ lsat->bands = 7;
+ for (i = 0; i < lsat->bands; i++) {
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = *(band + i);
+ lsat->band[i].wavemax = *(wmax + i);
+ lsat->band[i].wavemin = *(wmin + i);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 0.; /* Modified in set_TM5 by date */
+ lsat->band[i].thermal = (lsat->band[i].number == 6 ? 1 : 0);
+ }
+ return;
}
void sensor_ETM(lsat_data * lsat)
{
- int i;
+ 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 */
+ /* blue, green red, near infrared, shortwave IR, thermal IR, shortwave IR, panchromatic */
+ int band[] = { 1, 2, 3, 4, 5, 6, 6, 7, 8 };
+ int code[] = { 1, 2, 3, 4, 5, 61, 62, 7, 8 };
+ double wmax[] = { 0.515, 0.605, 0.690, 0.90, 1.75, 12.50, 2.35, 0.90 };
+ double wmin[] = { 0.450, 0.525, 0.630, 0.75, 1.55, 10.40, 2.09, 0.52 };
+ /* 30, 30, 30, 30, 30, 60, 30, 15 */
- strcpy(lsat->sensor, "ETM+");
+ strcpy(lsat->sensor, "ETM+");
- lsat->bands = 9;
- for (i = 0; i < lsat->bands; i++) {
- lsat->band[i].number = *(band + i);
- lsat->band[i].code = *(code + i);
- lsat->band[i].wavemax = *(wmax + i);
- lsat->band[i].wavemin = *(wmin + i);
- lsat->band[i].qcalmax = 255.;
- lsat->band[i].qcalmin = 1.;
- lsat->band[i].thermal = (lsat->band[i].number == 6 ? 1 : 0);
- }
- return;
+ lsat->bands = 9;
+ for (i = 0; i < lsat->bands; i++) {
+ lsat->band[i].number = *(band + i);
+ lsat->band[i].code = *(code + i);
+ lsat->band[i].wavemax = *(wmax + i);
+ lsat->band[i].wavemin = *(wmin + i);
+ lsat->band[i].qcalmax = 255.;
+ lsat->band[i].qcalmin = 1.;
+ lsat->band[i].thermal = (lsat->band[i].number == 6 ? 1 : 0);
+ }
+ return;
}
@@ -107,20 +107,20 @@
Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
/* Spectral radiances at detector */
double lmax[] = { 248., 200., 176., 153. };
- double lmin[] = { 0., 0., 0., 0. };
+ double lmin[] = { 0., 0., 0., 0. };
/* Solar exoatmospheric spectral irradiances */
double esun[] = { 1823., 1559., 1276., 880.1 };
lsat->number = 1;
- sensor_MSS( lsat );
+ 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);
+ 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;
}
@@ -137,29 +137,33 @@
/** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
/* Spectral radiances at detector */
- double Lmax[][4] = { { 210., 156., 140., 138. }, /* before July 16, 1975 */
- { 263., 176., 152., 130.333 } }; /* on or after July 16, 1975 */
- double Lmin[][4] = { { 10., 7., 7., 5. },
- { 8., 6., 6., 3.667 } };
+ double Lmax[][4] = { {210., 156., 140., 138.}, /* before July 16, 1975 */
+ {263., 176., 152., 130.333}
+ }; /* on or after July 16, 1975 */
+ double Lmin[][4] = { {10., 7., 7., 5.},
+ {8., 6., 6., 3.667}
+ };
/* Solar exoatmospheric spectral irradiances */
double esun[] = { 1829., 1539., 1268., 886.6 };
julian = julian_char(lsat->creation);
- if (julian < julian_char("1975-07-16")) i = 0;
- else i = 1;
+ if (julian < julian_char("1975-07-16"))
+ i = 0;
+ else
+ i = 1;
lmax = Lmax[i];
lmin = Lmin[i];
lsat->number = 2;
- sensor_MSS( lsat );
+ 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);
+ 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;
}
@@ -178,29 +182,33 @@
/** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
/* 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. } };
+ 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[] = { 1839., 1555., 1291., 887.9 };
julian = julian_char(lsat->creation);
- if (julian < julian_char("1978-06-01")) i = 0;
- else i = 1;
+ if (julian < julian_char("1978-06-01"))
+ i = 0;
+ else
+ i = 1;
lmax = Lmax[i];
lmin = Lmin[i];
lsat->number = 3;
- sensor_MSS( lsat );
+ 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);
+ 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;
}
@@ -217,32 +225,37 @@
/** Markham and Barker. EOSAT Landsat Technical Notes, No. 1, 1986;
Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009)*/
/* 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. } };
+ 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[] = { 1827., 1569., 1260., 866.4 };
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;
+ 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 );
+ 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);
+ 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;
}
@@ -255,38 +268,44 @@
/** 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.81, 204.30, 206.20, 27.19, 15.3032, 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.2378, -0.15 } };
+ 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.81, 204.30, 206.20, 27.19, 15.3032, 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.2378, -0.15}
+ };
+
/** Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009) */
/* Solar exoatmospheric spectral irradiances */
double esun[] = { 1983., 1795., 1539., 1028., 219.8, 0., 83.49 };
/* 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;
+ 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 );
+ sensor_TM(lsat);
lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
+ j = lsat->band[i].number - 1;
lsat->band[i].esun = *(esun + j);
lsat->band[i].lmax = *(lmax + j);
lsat->band[i].lmin = *(lmin + j);
- if (lsat->band[i].thermal ) {
- lsat->band[i].K1 = 671.62;
- lsat->band[i].K2 = 1284.30;
- }
+ if (lsat->band[i].thermal) {
+ lsat->band[i].K1 = 671.62;
+ lsat->band[i].K2 = 1284.30;
+ }
}
return;
}
@@ -304,32 +323,37 @@
/** 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. } };
+ 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[] = { 1824., 1570., 1249., 853.4 };
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;
+ 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 );
+ 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);
+ 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;
}
@@ -342,49 +366,58 @@
/** 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 } };
+ 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}
+ };
+
/** Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009) */
/* Solar exoatmospheric spectral irradiances */
double esun[] = { 1983., 1796., 1536., 1031., 220.0, 0., 83.44 };
/* 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;
+ 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;
- }
+ 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;
+ }
}
- jbuf = julian_char("2004-04-04");
- if (julian >= jbuf) G_warning("Using QCalMin=1.0 as a NLAPS product processed after 4/4/2004");
+ jbuf = julian_char("2004-04-04");
+ if (julian >= jbuf)
+ G_warning
+ ("Using QCalMin=1.0 as a NLAPS product processed after 4/4/2004");
lsat->number = 5;
- sensor_TM( lsat );
+ sensor_TM(lsat);
lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
- if (julian >= jbuf) lsat->band[i].qcalmin = 1.;
- 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;
- }
+ j = lsat->band[i].number - 1;
+ if (julian >= jbuf)
+ lsat->band[i].qcalmin = 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;
}
@@ -403,46 +436,54 @@
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.60, 244.0 }, /* before July 1, 2000 */
- { 293.7, 300.9, 234.4, 241.1, 47.57, 17.04, 16.54, 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 } };
+ double LmaxL[][8] = { {297.5, 303.4, 235.5, 235.0, 47.70, 17.04, 16.60, 244.0}, /* before July 1, 2000 */
+ {293.7, 300.9, 234.4, 241.1, 47.57, 17.04, 16.54, 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.80, 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 } };
+ 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.80, 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}
+ };
+
/** Chander, Markham and Helder. Remote Sensing of Environment, 113 (2009) */
/* Solar exoatmospheric spectral irradiances */
double esun[] = { 1997., 1812., 1533., 1039., 230.8, 0., 84.90, 1362. };
/* Thermal band calibration constants: K1 = 666.09 K2 = 1282.71 */
- julian = julian_char(lsat->creation);
- if (julian < julian_char("2000-07-01")) k = 0;
- else k = 1;
+ julian = julian_char(lsat->creation);
+ if (julian < julian_char("2000-07-01"))
+ k = 0;
+ else
+ k = 1;
- lsat->number = 7;
- sensor_ETM( lsat );
+ lsat->number = 7;
+ sensor_ETM(lsat);
- lsat->dist_es = earth_sun(lsat->date);
+ lsat->dist_es = earth_sun(lsat->date);
for (i = 0; i < lsat->bands; i++) {
- j = lsat->band[i].number - 1;
- lsat->band[i].esun = *(esun + j);
- if (gain[i] == 'H' || gain[i] == 'h') {
- lmax = LmaxH[k];
- lmin = LminH[k];
- }
- else {
- lmax = LmaxL[k];
- lmin = LminL[k];
- }
- lsat->band[i].lmax = *(lmax + j);
- lsat->band[i].lmin = *(lmin + j);
- if (lsat->band[i].thermal) {
- lsat->band[i].K1 = 666.09;
- lsat->band[i].K2 = 1282.71;
- }
+ j = lsat->band[i].number - 1;
+ lsat->band[i].esun = *(esun + j);
+ if (gain[i] == 'H' || gain[i] == 'h') {
+ lmax = LmaxH[k];
+ lmin = LminH[k];
+ }
+ else {
+ lmax = LmaxL[k];
+ lmin = LminL[k];
+ }
+ lsat->band[i].lmax = *(lmax + j);
+ lsat->band[i].lmin = *(lmin + j);
+ if (lsat->band[i].thermal) {
+ lsat->band[i].K1 = 666.09;
+ lsat->band[i].K2 = 1282.71;
+ }
}
return;
}
Modified: grass-addons/imagery/i.landsat.toar/local_proto.h
===================================================================
--- grass-addons/imagery/i.landsat.toar/local_proto.h 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/local_proto.h 2010-09-09 08:06:28 UTC (rev 43434)
@@ -4,8 +4,8 @@
#include <string.h>
#include "landsat.h"
-void met_ETM (char *, lsat_data *);
-void met_TM5 (char *, lsat_data *);
+void met_ETM(char *, lsat_data *);
+void met_TM5(char *, lsat_data *);
void set_MSS1(lsat_data *);
void set_MSS2(lsat_data *);
@@ -13,8 +13,8 @@
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 []);
+void set_TM4(lsat_data *);
+void set_TM5(lsat_data *);
+void set_ETM(lsat_data *, char[]);
#endif
Modified: grass-addons/imagery/i.landsat.toar/main.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/main.c 2010-09-09 07:51:43 UTC (rev 43433)
+++ grass-addons/imagery/i.landsat.toar/main.c 2010-09-09 08:06:28 UTC (rev 43434)
@@ -25,553 +25,591 @@
int main(int argc, char *argv[])
{
- struct History history;
- struct GModule *module;
+ struct History history;
+ struct GModule *module;
- struct Cell_head cellhd, window;
- char *mapset;
+ struct Cell_head cellhd, window;
+ char *mapset;
- void *inrast, *outrast;
- int infd, outfd;
- void *ptr;
- int nrows, ncols, row, col;
+ void *inrast, *outrast;
+ int infd, outfd;
+ void *ptr;
+ int nrows, ncols, row, col;
- RASTER_MAP_TYPE in_data_type;
+ RASTER_MAP_TYPE in_data_type;
- int verbose = 1;
- struct Option *input, *metfn,
- *adate, *pdate, *elev, *bgain, *metho, *perc, *dark,
- *satz, *atmo;
- char *name, *met;
- struct Flag *sat1, *sat2, *sat3, *sat4, *sat5, *sat7, *msss,
- *verbo, *frad;
+ int verbose = 1;
+ struct Option *input, *metfn,
+ *adate, *pdate, *elev, *bgain, *metho, *perc, *dark, *satz, *atmo;
+ char *name, *met;
+ struct Flag *sat1, *sat2, *sat3, *sat4, *sat5, *sat7, *msss,
+ *verbo, *frad;
- 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;
+ 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;
+ 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+");
+ 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->key = _("band_prefix");
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->gisprompt = _("input,cell,raster");
input->description = _("Base name of the landsat band rasters (.#)");
metfn = G_define_option();
- metfn->key = _("metfile");
- metfn->type = TYPE_STRING;
- metfn->required = NO;
- metfn->gisprompt = _(".met file");
+ metfn->key = _("metfile");
+ metfn->type = TYPE_STRING;
+ metfn->required = NO;
+ metfn->gisprompt = _(".met file");
metfn->description = _("Landsat ETM+ or TM5 header file (.met)");
adate = G_define_option();
- adate->key = _("date");
- adate->type = TYPE_STRING;
- adate->required = NO;
- adate->gisprompt = _("image acquisition date");
+ adate->key = _("date");
+ adate->type = TYPE_STRING;
+ adate->required = NO;
+ adate->gisprompt = _("image acquisition date");
adate->description = _("Image acquisition date (yyyy-mm-dd)");
elev = G_define_option();
- elev->key = _("solar_elevation");
- elev->type = TYPE_DOUBLE;
- elev->required = NO;
- elev->gisprompt = _("solar elevation");
+ elev->key = _("solar_elevation");
+ elev->type = TYPE_DOUBLE;
+ elev->required = NO;
+ elev->gisprompt = _("solar elevation");
elev->description = _("Solar elevation in degrees");
bgain = G_define_option();
- bgain->key = _("gain");
- bgain->type = TYPE_STRING;
- bgain->required = NO;
- bgain->gisprompt = _("band gain");
- bgain->description = _("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
+ bgain->key = _("gain");
+ bgain->type = TYPE_STRING;
+ bgain->required = NO;
+ bgain->gisprompt = _("band gain");
+ bgain->description =
+ _("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
pdate = G_define_option();
- pdate->key = _("product_date");
- pdate->type = TYPE_STRING;
- pdate->required = NO;
- pdate->gisprompt = _("image production date");
+ pdate->key = _("product_date");
+ pdate->type = TYPE_STRING;
+ pdate->required = NO;
+ pdate->gisprompt = _("image production date");
pdate->description = _("Image creation date (yyyy-mm-dd)");
metho = G_define_option();
- metho->key = _("method");
- metho->type = TYPE_STRING;
- metho->required = NO;
- metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
- metho->gisprompt = _("atmosferic correction method");
+ 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";
+ 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";
+ 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";
+ 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";
+ 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";
+ 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 */
+ /* It defines the different flags */
frad = G_define_flag();
- frad->key = 'r';
+ frad->key = 'r';
frad->description = _("Output at-sensor radiance for all bands");
- frad->answer = 0;
+ frad->answer = 0;
sat1 = G_define_flag();
- sat1->key = '1';
+ sat1->key = '1';
sat1->description = _("Landsat-1 MSS");
- sat1->answer = 0;
+ sat1->answer = 0;
sat2 = G_define_flag();
- sat2->key = '2';
+ sat2->key = '2';
sat2->description = _("Landsat-2 MSS");
- sat2->answer = 0;
+ sat2->answer = 0;
sat3 = G_define_flag();
- sat3->key = '3';
+ sat3->key = '3';
sat3->description = _("Landsat-3 MSS");
- sat3->answer = 0;
+ sat3->answer = 0;
sat4 = G_define_flag();
- sat4->key = '4';
+ sat4->key = '4';
sat4->description = _("Landsat-4 TM");
- sat4->answer = 0;
+ sat4->answer = 0;
sat5 = G_define_flag();
- sat5->key = '5';
+ sat5->key = '5';
sat5->description = _("Landsat-5 TM");
- sat5->answer = 0;
+ sat5->answer = 0;
sat7 = G_define_flag();
- sat7->key = '7';
+ sat7->key = '7';
sat7->description = _("Landsat-7 ETM+");
- sat7->answer = 0;
+ sat7->answer = 0;
msss = G_define_flag();
- msss->key = 's';
+ msss->key = 's';
msss->description = _("Set sensor of Landsat-4/5 to MSS");
- msss->answer = 0;
+ msss->answer = 0;
verbo = G_define_flag();
- verbo->key = 'v';
+ verbo->key = 'v';
verbo->description = _("Show parameters applied");
- verbo->answer = 0;
+ verbo->answer = 0;
/* options and afters parser */
if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ exit(EXIT_FAILURE);
/*****************************************
* ---------- START --------------------
* Stores options and flag to variables
*****************************************/
- met = metfn->answer;
- name = input->answer;
+ met = metfn->answer;
+ name = input->answer;
- if (adate->answer != NULL) {
- strncpy(lsat.date, adate->answer, 11);
- lsat.date[10] = '\0';
- }
- else
- lsat.date[0] = '\0';
+ if (adate->answer != NULL) {
+ strncpy(lsat.date, adate->answer, 11);
+ lsat.date[10] = '\0';
+ }
+ else
+ lsat.date[0] = '\0';
- if (pdate->answer != NULL) {
- strncpy(lsat.creation, pdate->answer, 11);
- lsat.creation[10] = '\0';
- }
- else
- lsat.creation[0] = '\0';
+ if (pdate->answer != NULL) {
+ 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);
+ 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"));
+ 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 */
+ 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"));
+ }
}
- /* 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 { /* Not need gain */
+ if (sat5->answer) {
+ if (msss->answer)
+ set_MSS5(&lsat);
+ else
+ set_TM5(&lsat);
+ G_message("Landsat-5 %s", lsat.sensor);
+ }
+ else if (sat4->answer) {
+ if (msss->answer)
+ set_MSS4(&lsat);
+ else
+ set_TM4(&lsat);
+ G_message("Landsat-4 %s", lsat.sensor);
+ }
+ else if (sat3->answer) {
+ set_MSS3(&lsat);
+ G_message("Landsat-3 MSS");
+ }
+ else if (sat2->answer) {
+ set_MSS2(&lsat);
+ G_message("Landsat-2 MSS");
+ }
+ else if (sat1->answer) {
+ set_MSS1(&lsat);
+ G_message("Landsat-1 MSS");
+ }
+ else {
+ G_fatal_error(_("Lacking satellite type"));
+ }
}
- else {
- 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 (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;
+ // 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;
+ 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"));
+ 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);
+ 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();
+ 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] %s",
- dn_dark[i], hist[dn_dark[i]],
- dn_mode[i], hist[dn_mode[i]],
- hist[255] > hist[dn_mode[i]] ? ", excluding DN > 241" : "");
+ G_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] %s",
+ dn_dark[i], hist[dn_dark[i]],
+ dn_mode[i], hist[dn_mode[i]],
+ hist[255] >
+ hist[dn_mode[i]] ? ", excluding DN > 241" : "");
- G_free(inrast);
- G_close_cell(infd);
- }
- /* Calculate transformation constants */
- lsat_bandctes(&lsat, i, method, percent, dn_dark[i], sat_zenith, rayleigh);
+ 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 (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, " percent of solar irradiance in path radiance = %.4lf\n", percent);
+ 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]);
}
- 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);
+ }
}
+ 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);
+ 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; }
+ mapset = G_find_cell2(band_in, "");
+ if (mapset == NULL) {
+ G_warning(_("raster file [%s] not found"), band_in);
+ continue;
+ }
- in_data_type = G_raster_map_type(band_in, mapset);
- if ((infd = G_open_cell_old(band_in, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), band_in);
+ in_data_type = G_raster_map_type(band_in, mapset);
+ if ((infd = G_open_cell_old(band_in, mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), band_in);
- if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), band_in);
+ if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), band_in);
- /* set same size as original band raster */
- if (G_set_window(&cellhd) < 0)
- G_fatal_error(_("Unable to set region"));
+ /* set same size as original band raster */
+ if (G_set_window(&cellhd) < 0)
+ G_fatal_error(_("Unable to set region"));
- /* controlling, if we can write the raster */
- if (G_legal_filename(band_out) < 0)
- G_fatal_error(_("[%s] is an illegal name"), band_out);
+ /* controlling, if we can write the raster */
+ if (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);
+ 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);
+ /* 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);
+ 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);
+ 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;
- }
+ 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]);
}
- if (G_put_raster_row(outfd, outrast, DCELL_TYPE) < 0)
- G_fatal_error(_("Cannot write to <%s>"), band_out);
+ else {
+ ref = lsat_rad2ref(rad, &lsat.band[i]);
+ if (ref < 0. && method > DOS)
+ ref = 0.;
+ }
+ }
+ ((DCELL *) outrast)[col] = ref;
}
- ref_mode = 0.;
- if (method > DOS && !lsat.band[i].thermal)
- {
- ref_mode = lsat_qcal2rad(dn_mode[i], &lsat.band[i]);
- ref_mode = lsat_rad2ref(ref_mode, &lsat.band[i]);
- }
+ }
+ if (G_put_raster_row(outfd, outrast, DCELL_TYPE) < 0)
+ G_fatal_error(_("Cannot write to <%s>"), band_out);
+ }
+ 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);
+ 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);
+ char command[300];
- G_short_history(band_out, "raster", &history);
+ sprintf(command, "r.colors map=%s color=grey", band_out);
+ system(command);
- 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_short_history(band_out, "raster", &history);
- G_command_history(&history);
- G_write_history(band_out, &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++;
- exit(EXIT_SUCCESS);
+ G_command_history(&history);
+ G_write_history(band_out, &history);
+ }
+
+ exit(EXIT_SUCCESS);
}
More information about the grass-commit
mailing list