[GRASS-SVN] r44035 - grass-addons/imagery/i.landsat.toar
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Oct 24 13:22:46 EDT 2010
Author: martinl
Date: 2010-10-24 10:22:46 -0700 (Sun, 24 Oct 2010)
New Revision: 44035
Modified:
grass-addons/imagery/i.landsat.toar/landsat.c
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:
i.landsat.toar sync'ed with trunk
Modified: grass-addons/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat.c 2010-10-24 17:00:48 UTC (rev 44034)
+++ grass-addons/imagery/i.landsat.toar/landsat.c 2010-10-24 17:22:46 UTC (rev 44035)
@@ -99,10 +99,9 @@
{
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].qcalmin)
+ / (lsat->band[i].qcalmax - lsat->band[i].qcalmin) +
lsat->band[i].lmin;
double Tv = 1.;
double Tz = 1.;
@@ -131,8 +130,8 @@
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);
+ G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz,
+ Edown);
lsat->band[i].K1 = 0.;
lsat->band[i].K2 = rad_sun;
Modified: grass-addons/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat_met.c 2010-10-24 17:00:48 UTC (rev 44034)
+++ grass-addons/imagery/i.landsat.toar/landsat_met.c 2010-10-24 17:22:46 UTC (rev 44035)
@@ -28,7 +28,7 @@
{
char *ptr;
- value[0] = 0;
+ value[0] = '\0';
ptr = strstr(mettext, text);
if (ptr == NULL)
@@ -139,6 +139,7 @@
FILE *f;
char mettext[TM5_MET_SIZE];
char value[MAX_STR];
+
/* char metdate[MAX_STR]; */
if ((f = fopen(metfile, "r")) == NULL)
@@ -162,11 +163,12 @@
get_value_met(mettext, "SolarElevation", value);
- if (!value)
+ if (!value[0])
G_warning("Unable to read solar elevation from metadata file");
else
lsat->sun_elev = atof(value);
- G_debug(1, "met_TM5: value=[%s], SolarElevation = %.2f", value, lsat->sun_elev);
+ G_debug(1, "met_TM5: value=[%s], SolarElevation = %.2f", value,
+ lsat->sun_elev);
get_value_met(mettext, "PLATFORMSHORTNAME", value);
@@ -203,3 +205,81 @@
(void)fclose(f);
return;
}
+
+
+
+/****************************************************************************
+ * PURPOSE: Read values of Landsat TM5 from header (.mtl) file
+ *****************************************************************************/
+
+/****************************************************************************
+ * EXPLANATION: This module is a modification of the met_ETM() found before
+ * to allow TM5 from GLOVIS to use .MTL extension that responds
+ * near to perfectly to the .met parser. While L7 files using
+ * .MTL from GLOVIS can be processed as if having .met files
+ * seemlessly, TM5 using .MTL need to read basic info and
+ * additionally the LMIN, LMAX, QCALMIN, QCALMAX being explicitely
+ * provided in the .MTL as if in a .met file.
+ *****************************************************************************/
+void mtl_TM5(char *metfile, lsat_data * lsat)
+{
+ FILE *f;
+ char mettext[ETM_MET_SIZE];
+ char name[MAX_STR], value[MAX_STR];
+ int i;
+
+ static int code[] = { 1, 2, 3, 4, 5, 6, 7 };
+
+ if ((f = fopen(metfile, "r")) == NULL)
+ G_fatal_error(_("Metadata file <%s> not found"), metfile);
+
+ fread(mettext, 1, ETM_MET_SIZE, f);
+
+ /* --------------------------------------- */
+ get_value_met7(mettext, "SENSOR_ID", value);
+ if (value[1] == 'M')
+ chrncpy(lsat->sensor, value + 1, 4);
+ else
+ chrncpy(lsat->sensor, value + 1, 3);
+
+ if (lsat->creation[0] == 0) {
+ get_value_met7(mettext, "PRODUCT_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);
+
+ /* We still have to initialize most of the info */
+ /* So instead of rewriting a new function, we use set_TM5()... */
+ set_TM5(lsat);
+ /* ... and we rewrite the necessary 'a la Landsat 7' */
+
+ if (strcmp(lsat->sensor, "MSS") == 0)
+ lsat->bands = 4;
+ else
+ lsat->bands = 7;
+ for (i = 0; i < lsat->bands; i++) {
+ lsat->band[i].code = *(code + i);
+ 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;
+ }
+ (void)fclose(f);
+ return;
+}
Modified: grass-addons/imagery/i.landsat.toar/landsat_set.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/landsat_set.c 2010-10-24 17:00:48 UTC (rev 44034)
+++ grass-addons/imagery/i.landsat.toar/landsat_set.c 2010-10-24 17:22:46 UTC (rev 44035)
@@ -45,7 +45,8 @@
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");
+ if (!lsat->sensor)
+ strcpy(lsat->sensor, "TM");
lsat->bands = 7;
for (i = 0; i < lsat->bands; i++) {
@@ -141,8 +142,8 @@
/* 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 */
+ {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.},
@@ -189,8 +190,8 @@
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 */
+ {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.},
@@ -235,9 +236,9 @@
/* 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 */
+ {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.},
@@ -281,9 +282,9 @@
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 */
+ {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},
@@ -340,9 +341,9 @@
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 */
+ {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.},
@@ -387,9 +388,9 @@
/* 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 */
+ {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},
@@ -422,10 +423,11 @@
}
jbuf = julian_char("2004-04-04");
- if (julian >= jbuf)
+ if (julian >= jbuf) {
G_warning
("Using QCalMin=1.0 as a NLAPS product processed after 4/4/2004");
-
+ G_warning("Discard this message if using the L5_MTL (-t) flag");
+ }
lsat->number = 5;
sensor_TM(lsat);
@@ -462,15 +464,15 @@
/* 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 */
+ {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] = {
+ 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}
};
Modified: grass-addons/imagery/i.landsat.toar/local_proto.h
===================================================================
--- grass-addons/imagery/i.landsat.toar/local_proto.h 2010-10-24 17:00:48 UTC (rev 44034)
+++ grass-addons/imagery/i.landsat.toar/local_proto.h 2010-10-24 17:22:46 UTC (rev 44035)
@@ -6,6 +6,7 @@
void met_ETM(char *, lsat_data *);
void met_TM5(char *, lsat_data *);
+void mtl_TM5(char *, lsat_data *);
void set_MSS1(lsat_data *);
void set_MSS2(lsat_data *);
Modified: grass-addons/imagery/i.landsat.toar/main.c
===================================================================
--- grass-addons/imagery/i.landsat.toar/main.c 2010-10-24 17:00:48 UTC (rev 44034)
+++ grass-addons/imagery/i.landsat.toar/main.c 2010-10-24 17:22:46 UTC (rev 44035)
@@ -5,15 +5,16 @@
*
* AUTHOR(S): E. Jorge Tizado - ej.tizado at unileon.es
* Hamish Bowman (small grassification cleanups)
+ * Yann Chemin (v7 + L5TM _MTL.txt support)
*
* PURPOSE: Calculate TOA Radiance or Reflectance and Kinetic Temperature
* for Landsat 1/2/3/4/5 MS, 4/5 TM or 7 ETM+
*
- * COPYRIGHT: (C) 2002, 2005 2008 by the GRASS Development Team
+ * COPYRIGHT: (C) 2002, 2005, 2008, 2010 by the GRASS Development Team
*
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
+ * This program is free software under the GNU General
+ * Public License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
*
*****************************************************************************/
@@ -28,34 +29,32 @@
{
struct History history;
struct GModule *module;
-
- struct Cell_head cellhd, window;
+
+ struct Cell_head cellhd;
char *mapset;
void *inrast, *outrast;
int infd, outfd;
void *ptr;
int nrows, ncols, row, col;
-
+
RASTER_MAP_TYPE in_data_type;
-
- int verbose = TRUE;
- struct Option *input, *metfn, *sensor, *adate, *pdate, *elev,
- *bgain, *metho, *perc, *dark, *satz, *atmo;
- char *name, *met;
- struct Flag *msss, *verbo, *frad;
-
+
+ struct Option *band_prefix, *output_suffix, *metfn, *sensor, *adate, *pdate, *elev,
+ *bgain, *metho, *perc, *dark, *satz, *atmo;
+ char *basename, *met, *suffixname, *sensorname;
+ struct Flag *msss, *frad, *l5_mtl;
+
lsat_data lsat;
- char band_in[127], band_out[127];
+ char band_in[GNAME_MAX], band_out[GNAME_MAX];
int i, j, q, method, pixel, dn_dark[MAX_BANDS], dn_mode[MAX_BANDS];
- int sensor_id;
double qcal, rad, ref, percent, ref_mode, sat_zenith, rayleigh;
-
+
struct Colors colors;
struct FPRange range;
double min, max;
unsigned long hist[256], h_max;
-
+
/* initialize GIS environment */
G_gisinit(argv[0]);
@@ -63,79 +62,96 @@
module = G_define_module();
module->description =
_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+.");
- module->keywords =
- _("imagery, landsat, top-of-atmosphere radiance or reflectance");
+ module->keywords = _("imagery, landsat, top-of-atmosphere reflectance, dos-type simple atmospheric correction");
/* 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 (.#)");
+ band_prefix = G_define_option();
+ band_prefix->key = "band_prefix";
+ band_prefix->label = _("Base name of input raster bands");
+ band_prefix->description = _("Example: 'B.' for B.1, B.2, ...");
+ band_prefix->type = TYPE_STRING;
+ band_prefix->required = YES;
- metfn = G_define_option();
+ output_suffix = G_define_option();
+ output_suffix->key = "output_suffix";
+ output_suffix->label = _("Suffix for output raster maps");
+ output_suffix->description = _("Example: '_toar' generates B.1_toar, B.2_toar, ...");
+ output_suffix->type = TYPE_STRING;
+ output_suffix->required = YES;
+
+ metfn = G_define_standard_option(G_OPT_F_INPUT);
metfn->key = "metfile";
- metfn->type = TYPE_STRING;
metfn->required = NO;
- metfn->gisprompt = "old_file,file,file";
- metfn->description = _("Landsat ETM+ or TM5 header file (.met)");
+ metfn->description = _("Name of Landsat ETM+ or TM5 header file (.met/.mtl)");
+ metfn->guisection = _("Metadata");
+ sensor = G_define_option();
+ sensor->key = "sensor";
+ sensor->type = TYPE_STRING;
+ sensor->label = _("Spacecraft sensor");
+ sensor->description = _("Required only if 'metfile' not given");
+ sensor->options = "mss1,mss2,mss3,tm4,tm5,tm7";
+ sensor->descriptions =
+ _("mss1;Landsat-1 MSS;"
+ "mss2;Landsat-2 MSS;"
+ "mss3;Landsat-3 MSS;"
+ "tm4;Landsat-4 TM;"
+ "tm5;Landsat-5 TM;"
+ "tm7;Landsat-7 ETM+");
+ sensor->required = NO;
+ sensor->guisection = _("Metadata");
+
metho = G_define_option();
metho->key = "method";
metho->type = TYPE_STRING;
metho->required = NO;
metho->options = "uncorrected,corrected,dos1,dos2,dos2b,dos3,dos4";
- metho->description = _("Atmospheric correction method");
+ metho->label = _("Atmospheric correction method");
+ metho->description = _("Required only if 'metfile' not given");
metho->answer = "uncorrected";
-
- sensor = G_define_option();
- sensor->key = "sensor";
- sensor->type = TYPE_INTEGER;
- sensor->description = _("Spacecraft sensor");
- sensor->options = "1,2,3,4,5,7";
- sensor->descriptions =
- _("1;Landsat-1 MSS;"
- "2;Landsat-2 MSS;"
- "3;Landsat-3 MSS;"
- "4;Landsat-4 TM;"
- "5;Landsat-5 TM;"
- "7;Landsat-7 ETM+");
- sensor->required = NO;
-
+ metho->guisection = _("Metadata");
+
adate = G_define_option();
adate->key = "date";
adate->type = TYPE_STRING;
adate->required = NO;
adate->key_desc = "yyyy-mm-dd";
- adate->description = _("Image acquisition date (yyyy-mm-dd)");
-
+ adate->label = _("Image acquisition date (yyyy-mm-dd)");
+ adate->description = _("Required only if 'metfile' not given");
+ adate->guisection = _("Metadata");
+
elev = G_define_option();
elev->key = "solar_elevation";
elev->type = TYPE_DOUBLE;
elev->required = NO;
- elev->description = _("Solar elevation in degrees");
+ elev->label = _("Solar elevation in degrees");
+ elev->description = _("Required only if 'metfile' not given");
+ elev->guisection = _("Metadata");
+ pdate = G_define_option();
+ pdate->key = "product_date";
+ pdate->type = TYPE_STRING;
+ pdate->required = NO;
+ pdate->key_desc = "yyyy-mm-dd";
+ pdate->label = _("Image creation date (yyyy-mm-dd)");
+ pdate->description = _("Required only if 'metfile' not given");
+ pdate->guisection = _("Metadata");
+
bgain = G_define_option();
bgain->key = "gain";
bgain->type = TYPE_STRING;
bgain->required = NO;
- bgain->description =
+ bgain->label =
_("Gain (H/L) of all Landsat ETM+ bands (1-5,61,62,7,8)");
+ bgain->guisection = _("Settings");
- pdate = G_define_option();
- pdate->key = "product_date";
- pdate->type = TYPE_STRING;
- pdate->required = NO;
- pdate->key_desc = "yyyy-mm-dd";
- pdate->description = _("Image creation date (yyyy-mm-dd)");
-
perc = G_define_option();
perc->key = "percent";
perc->type = TYPE_DOUBLE;
perc->required = NO;
perc->description = _("Percent of solar radiance in path radiance");
perc->answer = "0.01";
+ perc->guisection = _("Settings");
dark = G_define_option();
dark->key = "pixel";
@@ -144,6 +160,7 @@
dark->description =
_("Minimum pixels to consider digital number as dark object");
dark->answer = "1000";
+ dark->guisection = _("Settings");
satz = G_define_option();
satz->key = "sat_zenith";
@@ -151,26 +168,30 @@
satz->required = NO;
satz->description = _("Satellite zenith in degrees");
satz->answer = "8.2000";
+ satz->guisection = _("Settings");
atmo = G_define_option();
atmo->key = "rayleigh";
atmo->type = TYPE_DOUBLE;
atmo->required = NO;
- atmo->description = _("Rayleigh atmosphere"); /* scattering coefficient? */
+ atmo->description = _("Rayleigh atmosphere"); /* scattering coefficient? */
atmo->answer = "0.0";
+ atmo->guisection = _("Settings");
/* define the different flags */
frad = G_define_flag();
frad->key = 'r';
frad->description = _("Output at-sensor radiance for all bands");
-
+
msss = G_define_flag();
msss->key = 's';
- msss->description = _("Set sensor of Landsat-4/5 to MSS");
+ msss->description = _("Set sensor of Landsat TM4/5 to MSS");
+ msss->guisection = _("Settings");
- verbo = G_define_flag();
- verbo->key = 'v';
- verbo->description = _("Show parameters applied");
+ l5_mtl = G_define_flag();
+ l5_mtl->key = 't';
+ l5_mtl->description = _("Landsat TM5 has a .mtl file instead of .met");
+ l5_mtl->guisection = _("Metadata");
/* options and afters parser */
if (G_parser(argc, argv))
@@ -182,13 +203,18 @@
* Stores options and flag to variables
*****************************************/
met = metfn->answer;
- name = input->answer;
-
+ basename = band_prefix->answer;
+ suffixname = output_suffix->answer;
+ sensorname = sensor -> answer ? sensor->answer: "";
+
+ G_zero(&lsat, sizeof(lsat));
+
if (adate->answer != NULL) {
strncpy(lsat.date, adate->answer, 11);
lsat.date[10] = '\0';
if (strlen(lsat.date) != 10)
- G_fatal_error(_("Illegal date format: [%s] (yyyy-mm-dd)"), lsat.date);
+ G_fatal_error(_("Illegal date format: [%s] (yyyy-mm-dd)"),
+ lsat.date);
}
else
lsat.date[0] = '\0';
@@ -197,7 +223,8 @@
strncpy(lsat.creation, pdate->answer, 11);
lsat.creation[10] = '\0';
if (strlen(lsat.creation) != 10)
- G_fatal_error(_("Illegal date format: [%s] (yyyy-mm-dd)"), lsat.creation);
+ G_fatal_error(_("Illegal date format: [%s] (yyyy-mm-dd)"),
+ lsat.creation);
}
else
lsat.creation[0] = '\0';
@@ -208,70 +235,69 @@
sat_zenith = atof(satz->answer);
rayleigh = atof(atmo->answer);
- if (sensor->answer)
- sensor_id = atoi(sensor->answer);
- else
- G_fatal_error(_("Must select type of satellite"));
-
/* Data from MET file: only Landsat-7 ETM+ and Landsat-5 TM */
if (met != NULL) {
- if (sensor_id == 7)
+ if (strcmp(sensorname, "tm7") == 0)
met_ETM(met, &lsat);
+ else if (l5_mtl->answer)
+ mtl_TM5(met, &lsat);
else
met_TM5(met, &lsat);
- G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number, lsat.sensor);
+ G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number,
+ lsat.sensor);
if (!lsat.sensor || lsat.number > 7 || lsat.number < 1)
G_fatal_error(_("Failed to identify satellite"));
- G_message(_("Landsat-%d %s with data set in met file [%s]"),
+ G_debug(1, "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"));
+ G_fatal_error(_("Lacking '%s' or '%s' for this satellite"),
+ adate->key, elev->key);
}
else {
- if (sensor_id == 7) { /* Need gain */
+ if (strcmp(sensorname, "tm7") == 0) { /* Need gain */
if (bgain->answer != NULL && strlen(bgain->answer) == 9) {
set_ETM(&lsat, bgain->answer);
- G_message("Landsat 7 ETM+");
+ G_debug(1, "Landsat 7 ETM+");
}
else {
G_fatal_error(_("Landsat-7 requires band gain with 9 (H/L) data"));
}
}
else { /* Not need gain */
- if (sensor_id == 5) {
+ if (strcmp(sensorname, "tm5") == 0) {
if (msss->answer)
set_MSS5(&lsat);
else
set_TM5(&lsat);
- G_message("Landsat-5 %s", lsat.sensor);
+ G_debug(1, "Landsat-5 %s", lsat.sensor);
}
- else if (sensor_id == 4) {
+ else if (strcmp(sensorname, "tm4") == 0) {
if (msss->answer)
set_MSS4(&lsat);
else
set_TM4(&lsat);
- G_message("Landsat-4 %s", lsat.sensor);
+ G_debug(1, "Landsat-4 %s", lsat.sensor);
}
- else if (sensor_id == 3) {
+ else if (strcmp(sensorname, "mss3") == 0) {
set_MSS3(&lsat);
- G_message("Landsat-3 MSS");
+ G_debug(1, "Landsat-3 MSS");
}
- else if (sensor_id == 2) {
+ else if (strcmp(sensorname, "mss2") == 0) {
set_MSS2(&lsat);
- G_message("Landsat-2 MSS");
+ G_debug(1, "Landsat-2 MSS");
}
- else if (sensor_id == 1) {
+ else if (strcmp(sensorname, "mss1") == 0) {
set_MSS1(&lsat);
- G_message("Landsat-1 MSS");
+ G_debug(1, "Landsat-1 MSS");
}
else {
- G_fatal_error(_("Unknown satellite type"));
+ G_fatal_error(_("Unknown satellite type (defined by '%s')"), sensor->key);
}
}
}
@@ -294,14 +320,14 @@
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;
@@ -311,18 +337,17 @@
for (j = 0; j < 256; j++)
hist[j] = 0L;
- snprintf(band_in, 127, "%s.%d", name, lsat.band[i].code);
+ sprintf(band_in, "%s%d", basename, lsat.band[i].code);
mapset = G_find_cell2(band_in, "");
if (mapset == NULL) {
G_warning(_("Raster map <%s> not found"), band_in);
continue;
}
- if ((infd = G_open_cell_old(band_in, mapset)) < 0)
+ if ((infd = G_open_cell_old(band_in, "")) < 0)
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
G_fatal_error(_("Unable to read header of raster map <%s>"), band_in);
- if (G_set_window(&cellhd) < 0)
- G_fatal_error(_("Cannot reset current region"));
+ G_set_window(&cellhd);
in_data_type = G_raster_map_type(band_in, mapset);
inrast = G_allocate_raster_buf(in_data_type);
@@ -330,7 +355,7 @@
nrows = G_window_rows();
ncols = G_window_cols();
- G_message("Calculating dark pixel of [%s] ... ", band_in);
+ 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++) {
@@ -355,7 +380,7 @@
}
/* DN of dark object */
for (j = lsat.band[i].qcalmin; j < 256; j++) {
- if (hist[j] >= pixel) {
+ if (hist[j] >= (unsigned int) pixel) {
dn_dark[i] = j;
break;
}
@@ -369,11 +394,11 @@
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_verbose_message("... DN = %.2d [%lu] : mode %.2d [%lu] %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);
@@ -384,69 +409,68 @@
}
if (strlen(lsat.creation) == 0)
- G_fatal_error(_("Unknown production date"));
+ G_fatal_error(_("Unknown production date (defined by '%s')"), pdate->key);
- /*****************************************
- * ------------ VERBOSE ------------------
- *****************************************/
- if (verbo->answer) {
- fprintf(stdout, " ACQUISITION DATE %s [production date %s]\n",
+ if (G_verbose() > G_verbose_std()) {
+ fprintf(stderr, " SENSOR: %s\n", lsat.sensor);
+ fprintf(stderr, " 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",
+ fprintf(stderr, " earth-sun distance = %.8lf\n", lsat.dist_es);
+ fprintf(stderr, " solar elevation angle = %.8lf\n", lsat.sun_elev);
+ fprintf(stderr, " Method of calculus = %s\n",
(method == CORRECTED ? "CORRECTED"
: (method == UNCORRECTED ? "UNCORRECTED" : metho->answer)));
if (method > DOS) {
- fprintf(stdout,
+ fprintf(stderr,
" 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",
+ fprintf(stderr, "-------------------\n");
+ fprintf(stderr, " BAND %d %s (code %d)\n",
lsat.band[i].number,
(lsat.band[i].thermal ? "thermal " : ""),
lsat.band[i].code);
- fprintf(stdout,
+ fprintf(stderr,
" 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",
+ fprintf(stderr, " calibration constants (L): %.3lf to %.3lf\n",
lsat.band[i].lmin, lsat.band[i].lmax);
- fprintf(stdout, " at-%s radiance = %.5lf * DN + %.5lf\n",
+ fprintf(stderr, " 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,
+ fprintf(stderr,
" at-sensor temperature = %.3lf / log[(%.3lf / radiance) + 1.0]\n",
lsat.band[i].K2, lsat.band[i].K1);
}
else {
- fprintf(stdout,
+ fprintf(stderr,
" mean solar exoatmospheric irradiance (ESUN): %.3lf\n",
lsat.band[i].esun);
- fprintf(stdout, " at-%s reflectance = radiance / %.5lf\n",
+ fprintf(stderr, " at-%s reflectance = radiance / %.5lf\n",
(method > DOS ? "surface" : "sensor"),
lsat.band[i].K2);
if (method > DOS) {
- fprintf(stdout,
+ fprintf(stderr,
" 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(stderr, " the mode of DN is %d\n", dn_mode[i]);
}
}
}
- fprintf(stdout, "-------------------\n");
- fflush(stdout);
+ fprintf(stderr, "-------------------\n");
+ fflush(stderr);
}
/*****************************************
* ------------ CALCULUS -----------------
*****************************************/
+ G_message(_("Calculating..."));
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);
+ sprintf(band_in, "%s%d", basename, lsat.band[i].code);
+ sprintf(band_out, "%s%d%s", basename, lsat.band[i].code, suffixname);
mapset = G_find_cell2(band_in, "");
if (mapset == NULL) {
@@ -454,16 +478,14 @@
continue;
}
- in_data_type = G_raster_map_type(band_in, mapset);
if ((infd = G_open_cell_old(band_in, mapset)) < 0)
G_fatal_error(_("Unable to open raster map <%s>"), band_in);
-
+ in_data_type = G_raster_map_type(band_in, mapset);
if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
G_fatal_error(_("Unable to read header of raster map <%s>"), band_in);
-
+
/* set same size as original band raster */
- if (G_set_window(&cellhd) < 0)
- G_fatal_error(_("Cannot reset current region"));
+ G_set_window(&cellhd);
/* controlling, if we can write the raster */
if (G_legal_filename(band_out) < 0)
@@ -479,14 +501,13 @@
nrows = G_window_rows();
ncols = G_window_cols();
/* ================================================================= */
- G_message(_("Writing %s of <%s> to <%s> ..."),
+ G_message(_("Writing %s of <%s> to <%s>..."),
(frad->answer ? _("radiance")
- : (lsat.band[i].thermal) ? _("temperature") : _("reflectance")),
+ : (lsat.band[i].
+ thermal) ? _("temperature") : _("reflectance")),
band_in, band_out);
-
for (row = 0; row < nrows; row++) {
- if (verbose)
- G_percent(row, nrows, 2);
+ G_percent(row, nrows, 2);
G_get_raster_row(infd, inrast, row, in_data_type);
for (col = 0; col < ncols; col++) {
@@ -526,9 +547,10 @@
((DCELL *) outrast)[col] = ref;
}
}
- if (G_put_raster_row(outfd, outrast, DCELL_TYPE) < 0)
- G_fatal_error(_("Failed writing raster map <%s> row %d"), band_out, row);
+ G_put_raster_row(outfd, outrast, DCELL_TYPE);
}
+ G_percent(1, 1, 1);
+
ref_mode = 0.;
if (method > DOS && !lsat.band[i].thermal) {
ref_mode = lsat_qcal2rad(dn_mode[i], &lsat.band[i]);
@@ -542,19 +564,19 @@
G_free(outrast);
G_close_cell(outfd);
-/* needed ?
- if (out_type != CELL_TYPE)
- G_quantize_fp_map_range(band_out, G_mapset(), 0., 360., 0, 360);
-*/
+ /* needed ?
+ if (out_type != CELL_TYPE)
+ G_quantize_fp_map_range(band_out, G_mapset(), 0., 360., 0, 360);
+ */
/* set grey255 colortable */
G_init_colors(&colors);
- G_read_fp_range(band_out, G_mapset(), &range);
- G_get_fp_range_min_max(&range, &min, &max);
- G_make_grey_scale_fp_colors(&colors, min, max);
+ G_read_fp_range(band_out, G_mapset(), &range);
+ G_get_fp_range_min_max(&range, &min, &max);
+ G_make_grey_scale_fp_colors(&colors, min, max);
G_write_colors(band_out, G_mapset(), &colors);
+ /* Initialize the 'hist' structure with basic info */
G_short_history(band_out, "raster", &history);
-
sprintf(history.edhist[0], " %s of Landsat-%d %s (method %s)",
lsat.band[i].thermal ? "Temperature" : "Reflectance",
lsat.number, lsat.sensor, metho->answer);
@@ -605,17 +627,14 @@
sprintf(history.edhist[history.edlinecnt],
"-----------------------------------------------------------------");
history.edlinecnt++;
-
+
G_command_history(&history);
G_write_history(band_out, &history);
- G_put_cell_title(band_out, history.edhist[0]);
-
- if(lsat.band[i].thermal)
+ if (lsat.band[i].thermal)
G_write_raster_units(band_out, "Kelvin");
/* else units = ...? */
/* set raster timestamp from acq date? (see r.timestamp module) */
-
}
exit(EXIT_SUCCESS);
More information about the grass-commit
mailing list