[GRASS-SVN] r58094 - grass/branches/develbranch_6/imagery/i.landsat.toar

svn_grass at osgeo.org svn_grass at osgeo.org
Wed Oct 23 14:21:12 PDT 2013


Author: neteler
Date: 2013-10-23 14:21:12 -0700 (Wed, 23 Oct 2013)
New Revision: 58094

Modified:
   grass/branches/develbranch_6/imagery/i.landsat.toar/description.html
   grass/branches/develbranch_6/imagery/i.landsat.toar/landsat.c
   grass/branches/develbranch_6/imagery/i.landsat.toar/landsat_met.c
   grass/branches/develbranch_6/imagery/i.landsat.toar/main.c
Log:
i.landsat.toar: fixes for Landsat-8 metadata file support (author: E. Jorge Tizado)

Modified: grass/branches/develbranch_6/imagery/i.landsat.toar/description.html
===================================================================
--- grass/branches/develbranch_6/imagery/i.landsat.toar/description.html	2013-10-23 21:16:40 UTC (rev 58093)
+++ grass/branches/develbranch_6/imagery/i.landsat.toar/description.html	2013-10-23 21:21:12 UTC (rev 58094)
@@ -97,13 +97,66 @@
     TAUz = exp[-t/sin(e)], Esky = PI * radiance_dark </li>
 </ul>
 
-<b>Attention</b>: Output radiance remain untouched (i.e. no set to
-0. when it is negative) then they are possible negative
-values. However, output reflectance is set to 0. when is obtained a
-negative value.
+<b>Attention</b>: Output radiance remain untouched (i.e. no set to 0.0
+when it is negative) then they are possible negative values. However,
+output reflectance is set to 0.0 when is obtained a negative value.
 
 <h2>NOTES</h2>
 
+<h3>On Landsat-8 metadata file </h3>
+
+NASA reports a structure of the L1G Metadata file 
+(<a href="http://landsat.usgs.gov/documents/LDCM-DFCB-004.pdf‎">LDCM-DFCB-004.pdf</a>) 
+for Landsat Data Continuity Mission (i.e. Landsat-8).
+
+<p>
+NASA retains in MIN_MAX_RADIANCE group the necessary information 
+to transform Digital Numbers (DN) in radiance values. Then, 
+<em>i.landsat.toar</em> replaces the possible standard values with the 
+metadata values. The results match with the values reported by the 
+metada file in RADIOMETRIC_RESCALING group.
+
+<p>
+Also, NASA reports the same values of reflectance for all bands 
+in max-min values and in gain-bias values. This is strange that all 
+bands have the same range of reflectance. Also, they wrote in the 
+web page as to calculate reflectance directly from DN, first with 
+RADIOMETRIC_RESCALING values and second divided by sin(sun_elevation).
+
+<p>
+This is a simple rescaling 
+
+<ul>
+  <li> reflectance = radiance / sun_radiance = (DN * RADIANCE_MULT + RADIANCE_ADD) / sun_radiance</li>
+  <li> now reflectance = DN * REFLECTANCE_MULT + REFLECTANCE_ADD </li>
+  <li> then REFLECTANCE_MULT = RADIANCE_MULT / sun_radiance </li>
+  <li> and REFLECTANCE_ADD = RADIANCE_ADD / sun_radiance </li>
+</ul>
+
+<p>
+The problem arises when we need ESUN values (not provided) to 
+compute sun_radiance and DOS. We assume that REFLECTANCE_MAXIMUM 
+corresponds to the RADIANCE_MAXIMUM, then
+
+<ul>
+  <li> REFLECTANCE_MAXIMUM / sin(e) = RADIANCE_MAXIMUM / sun_radiance</li>
+  <li> Esun = (PI * d^2) * RADIANCE_MAXIMUM / REFLECTANCE_MAXIMUM </li>
+</ul>
+
+where <em>d</em> is the earth-sun distance provided by metadata file 
+or computed inside the program.
+
+<p>
+The <em>i.landsat.toar</em> reverts back the NASA rescaling to continue 
+using Lmax, Lmin, and Esun values to compute the constant to convert 
+DN to radiance and radiance to reflectance with the "traditional" 
+equations and simple atmospheric corrections.
+
+<b>Attention</b>: When MAXIMUM values are not provided,
+<em>i.landsat.toar</em> tries to calculate Lmax, Lmin, and Esun from 
+RADIOMETRIC_RESCALING (in tests the results were the same).
+
+<h3>Calibration constants</h3>
 In verbose mode (flag <b>--verbose</b>), the program write basic
 satellite data and the parameters used in the transformations.
 
@@ -150,6 +203,13 @@
   sun_elevation=64.3242970 gain="HHHLHLHHL"
 </pre></div>
 
+or
+
+<div class="code"><pre>
+i.landsat.toar input_prefix=LC80160352013134LGN03_B output_prefix=toar \
+  metfile=LC80160352013134LGN03_MTL.txt sensor=oli8 date=2013-05-14
+</pre></div>
+
 <h2>REFERENCES</h2>
 
 <ul>

Modified: grass/branches/develbranch_6/imagery/i.landsat.toar/landsat.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.landsat.toar/landsat.c	2013-10-23 21:16:40 UTC (rev 58093)
+++ grass/branches/develbranch_6/imagery/i.landsat.toar/landsat.c	2013-10-23 21:21:12 UTC (rev 58094)
@@ -87,8 +87,7 @@
 		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);
+		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;
@@ -110,13 +109,9 @@
 		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);
+		    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 (TAUv != Tv && TAUz != Tz);
 		TAUz = (Tz < 1. ? Tz : 1.);
 		TAUv = (Tv < 1. ? Tv : 1.);
@@ -130,42 +125,27 @@
 	    break;
 	}
 	lsat->band[i].K2 = 0.;
-	lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;	/* rad_sun */
+	lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
 	if (method > DOS)
-	    G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n",
-			      TAUv, TAUz, Edown);
+	    G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz, Edown);
     }
 
     /** Digital number to radiance coefficients.
 	 * Without atmospheric 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);
+	lsat->band[i].bias = (lsat->band[i].lmin - lsat->band[i].gain * lsat->band[i].qcalmin);
     }
-    /*
-       else {
-       if (method == CORRECTED) {
-       // L = G * (DN - Qmin) + Lmin - Lmin
-       -> bias = - G * Qmin *
-       lsat->band[i].bias =
-       -(lsat->band[i].gain * lsat->band[i].qcalmin);
-       // Another possibility is cut when rad < 0 *
-       }
-     */
     else if (method > DOS) {
 	/* L = Lsat - Lpath = G * DNsat + B - (G *  + B - p * rad_sun) 
 	 *   = G * DNsat - G *  + p * rad_sun
 	 *  -> bias = p * rad_sun - G 
 	 */
-	lsat->band[i].bias =
-	    percent * lsat->band[i].K1 - lsat->band[i].gain * dark;
+	lsat->band[i].bias = percent * lsat->band[i].K1 - lsat->band[i].gain * dark;
     }
 }

Modified: grass/branches/develbranch_6/imagery/i.landsat.toar/landsat_met.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.landsat.toar/landsat_met.c	2013-10-23 21:16:40 UTC (rev 58093)
+++ grass/branches/develbranch_6/imagery/i.landsat.toar/landsat_met.c	2013-10-23 21:21:12 UTC (rev 58094)
@@ -35,6 +35,7 @@
     }
 }
 
+/* OLD Metadata Files */
 void get_metformat(const char metadata[], char *key, char value[])
 {
     int i = 0;
@@ -56,6 +57,17 @@
     value[i] = '\0';
 }
 
+double get_metdouble(const char metadata[], char *format, int code, char value[])
+{
+    char key[MAX_STR];
+    
+    sprintf(key, format, code);
+    get_metformat(metadata, key, value);
+    return atof(value);			    
+}
+
+
+/* NEW Metadata Files */
 void get_mtlformat(const char metadata[], char *key, char value[])
 {
     int i = 0;
@@ -73,7 +85,17 @@
     value[i] = '\0';
 }
 
+double get_mtldouble(const char metadata[], char *format, int code, char value[])
+{
+    char key[MAX_STR];
+    
+    sprintf(key, format, code);
+    get_mtlformat(metadata, key, value);
+    return atof(value);			    
+}
 
+
+
 /****************************************************************************
  * PURPOSE:     Read parameters from Landsat metadata files
  *****************************************************************************/
@@ -84,7 +106,9 @@
     char mtldata[METADATA_SIZE];
     char key[MAX_STR], value[MAX_STR];
     void (*get_mtldata) (const char[], char *, char[]);
+    void (*get_mtlreal) (const char[], char *, int, char[]);
     int i, j, ver_mtl;
+    double X2;
 
     /* store metadata in ram */
     if ((f = fopen(metafile, "r")) == NULL)
@@ -93,8 +117,17 @@
     (void)fclose(f);
 
     /* set version of the metadata file */
-    get_mtldata =
-	(strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat;
+    /* get_mtldata = (strstr(mtldata, " VALUE ") != NULL) ? get_metformat : get_mtlformat; */
+    if (strstr(mtldata, " VALUE ") != NULL)
+    {
+	get_mtldata = get_metformat;
+	get_mtlreal = get_metdouble;
+    }
+    else
+    {
+	get_mtldata = get_mtlformat;
+	get_mtlreal = get_mtldouble;
+    }
     ver_mtl = (strstr(mtldata, "QCALMAX_BAND") != NULL) ? 0 : 1;
 
     /* Fill with product metadata */
@@ -238,120 +271,93 @@
 	else {
 	    G_verbose_message("Metada file is MTL file: new format");
 	    /* Other possible values in the metadata file */
-	    get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value);	/* used after in calculus after */
+	    get_mtldata(mtldata, "EARTH_SUN_DISTANCE", value);	/* Necessary after */
 	    if (value[0] != '\0')
 		lsat->dist_es = atof(value);
-	    /* ----- */
-	    char *fmt_radmu[] =
-		{ "RADIANCE_MULTIPLICATIVE_FACTOR_BAND%d",
-    "RADIANCE_MULT_BAND_%d" };
-	    char *fmt_radad[] =
-		{ "RADIANCE_ADDITIVE_FACTOR_BAND%d", "RADIANCE_ADD_BAND_%d" };
-	    char *fmt_k1cte[] =
-		{ "K1_CONSTANT_BAND%d", "K1_CONSTANT_BAND_%d" };
-	    char *fmt_k2cte[] =
-		{ "K2_CONSTANT_BAND%d", "K2_CONSTANT_BAND_%d" };
-	    char *fmt_refmu[] =
-		{ "REFLECTANCE_MULTIPLICATIVE_FACTOR_BAND%d",
-    "REFLECTANCE_MULT_BAND_%d" };
-	    char *fmt_refad[] =
-		{ "REFLECTANCE_ADDITIVE_FACTOR_BAND%d",
-    "REFLECTANCE_ADD_BAND_%d" };
-	    /* --NASA-- LDCM sample file: mode = 0; LDCM-DFCB-004.pdf file: mode = 1 */
-	    int mode =
-		(strstr(mtldata, "RADIANCE_MULTIPLICATIVE_FACTOR_BAND") !=
-		 NULL) ? 0 : 1;
-	    /* ----- */
-	    if (strstr(mtldata, "RADIANCE_MAXIMUM_BAND") != NULL) {
-		G_verbose_message
-		    ("RADIANCE & QUANTIZE from data of the metadata file");
+	    if (strstr(mtldata, "RADIANCE_MAXIMUM_BAND") != NULL) 
+	    {
+		G_verbose_message("RADIANCE & QUANTIZE from  MIN_MAX_(RADIANCE|PIXEL_VALUE)");
 		for (i = 0; i < lsat->bands; i++) {
 		    if (lsat->number == 7 && lsat->band[i].thermal) {
-			sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d",
-				lsat->band[i].code - 60);
+			sprintf(key, "RADIANCE_MAXIMUM_BAND_6_VCID_%d", lsat->band[i].code - 60);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].lmax = atof(value);
-			sprintf(key, "RADIANCE_MINIMUM_BAND_6_VCID_%d",
-				lsat->band[i].code - 60);
+			sprintf(key, "RADIANCE_MINIMUM_BAND_6_VCID_%d", lsat->band[i].code - 60);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].lmin = atof(value);
-			sprintf(key, "QUANTIZE_CAL_MAX_BAND_6_VCID_%d",
-				lsat->band[i].code - 60);
+			sprintf(key, "QUANTIZE_CAL_MAX_BAND_6_VCID_%d", lsat->band[i].code - 60);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmax = atof(value);
-			sprintf(key, "QUANTIZE_CAL_MIN_BAND_6_VCID_%d",
-				lsat->band[i].code - 60);
+			sprintf(key, "QUANTIZE_CAL_MIN_BAND_6_VCID_%d", lsat->band[i].code - 60);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmin = atof(value);
 		    }
 		    else {
-			sprintf(key, "RADIANCE_MAXIMUM_BAND_%d",
-				lsat->band[i].code);
+			sprintf(key, "RADIANCE_MAXIMUM_BAND_%d", lsat->band[i].code);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].lmax = atof(value);
-			sprintf(key, "RADIANCE_MINIMUM_BAND_%d",
-				lsat->band[i].code);
+			sprintf(key, "RADIANCE_MINIMUM_BAND_%d", lsat->band[i].code);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].lmin = atof(value);
-			sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d",
-				lsat->band[i].code);
+			sprintf(key, "QUANTIZE_CAL_MAX_BAND_%d", lsat->band[i].code);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmax = atof(value);
-			sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d",
-				lsat->band[i].code);
+			sprintf(key, "QUANTIZE_CAL_MIN_BAND_%d", lsat->band[i].code);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].qcalmin = atof(value);
 		    }
 		    /* other possible values of each band */
 		    if (lsat->band[i].thermal) {
-			sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+			sprintf(key, "K1_CONSTANT_BAND_%d", lsat->band[i].code);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].K1 = atof(value);
-			sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+			sprintf(key, "K2_CONSTANT_BAND_%d", lsat->band[i].code);
 			get_mtldata(mtldata, key, value);
 			lsat->band[i].K2 = atof(value);
 		    }
+		    else if (lsat->number == 8)
+		    {
+			/* ESUN from  REFLECTANCE and RADIANCE ADD_BAND */
+			sprintf(key, "REFLECTANCE_MAXIMUM_BAND_%d", lsat->band[i].code);
+			get_mtldata(mtldata, key, value);
+			X2 = atof(value);
+			lsat->band[i].esun = (double)(PI * lsat->dist_es * lsat->dist_es * lsat->band[i].lmax) / X2;
+		    }
 		}
+		if (lsat->number == 8)
+		    G_warning("ESUN evaluated from REFLECTANCE_MAXIMUM_BAND");
 	    }
 	    else {
-		G_verbose_message
-		    ("RADIANCE & QUANTIZE from radiometric rescaling group of the metadata file");
+		G_verbose_message("RADIANCE & QUANTIZE from RADIOMETRIC_RESCALING");
 		for (i = 0; i < lsat->bands; i++) {
-		    sprintf(key, fmt_radmu[mode], lsat->band[i].code);
+		    sprintf(key, "RADIANCE_MULT_BAND_%d", lsat->band[i].code);
 		    get_mtldata(mtldata, key, value);
 		    lsat->band[i].gain = atof(value);
-		    sprintf(key, fmt_radad[mode], lsat->band[i].code);
+		    sprintf(key, "RADIANCE_ADD_BAND_%d", lsat->band[i].code);
 		    get_mtldata(mtldata, key, value);
 		    lsat->band[i].bias = atof(value);
 		    /* reversing to calculate the values of Lmin and Lmax ... */
-		    lsat->band[i].lmin =
-			lsat->band[i].gain * lsat->band[i].qcalmin +
-			lsat->band[i].bias;
-		    lsat->band[i].lmax =
-			lsat->band[i].gain * lsat->band[i].qcalmax +
-			lsat->band[i].bias;
+		    lsat->band[i].lmin = lsat->band[i].gain * lsat->band[i].qcalmin + lsat->band[i].bias;
+		    lsat->band[i].lmax = lsat->band[i].gain * lsat->band[i].qcalmax + lsat->band[i].bias;
 		    /* ... qcalmax and qcalmin loaded in previous sensor_ function */
 		    if (lsat->number == 8) {
 			if (lsat->band[i].thermal) {
-			    sprintf(key, fmt_k1cte[mode], lsat->band[i].code);
+			    sprintf(key, "K1_CONSTANT_BAND_%d", lsat->band[i].code);
 			    get_mtldata(mtldata, key, value);
 			    lsat->band[i].K1 = atof(value);
-			    sprintf(key, fmt_k2cte[mode], lsat->band[i].code);
+			    sprintf(key, "K2_CONSTANT_BAND_%d", lsat->band[i].code);
 			    get_mtldata(mtldata, key, value);
 			    lsat->band[i].K2 = atof(value);
 			}
 			else {
-			    sprintf(key, fmt_refmu[mode], lsat->band[i].code);
+			    sprintf(key, "REFLECTANCE_MULT_BAND_%d", lsat->band[i].code);
 			    get_mtldata(mtldata, key, value);
 			    lsat->band[i].K1 = atof(value);
-			    sprintf(key, fmt_refad[mode], lsat->band[i].code);
+			    sprintf(key, "REFLECTANCE_ADD_BAND_%d", lsat->band[i].code);
 			    get_mtldata(mtldata, key, value);
-			    lsat->band[i].K2 = atof(value);
-			    /* ESUN evaluates from metadafa file */
-			    lsat->band[i].esun =
-				(double)(PI * lsat->dist_es * lsat->dist_es *
-					 lsat->band[i].bias) /
-				lsat->band[i].K2;
+			    lsat->band[i].K2 = atof(value);			    
+			    /* ESUN from  REFLECTANCE_ADD_BAND */
+			    lsat->band[i].esun = (double)(PI * lsat->dist_es * lsat->dist_es * lsat->band[i].bias) / lsat->band[i].K2;
 			    /*
 			       double esun1 = (double) (PI * lsat->dist_es * lsat->dist_es * lsat->band[i].bias) / lsat->band[i].K2;
 			       double esun2 = (double) (PI * lsat->dist_es * lsat->dist_es * lsat->band[i].gain) / lsat->band[i].K1;
@@ -360,8 +366,7 @@
 			}
 		    }
 		}
-		G_warning
-		    ("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
+		G_warning("ESUN evaluated from REFLECTANCE_ADDITIVE_FACTOR_BAND");
 	    }
 	}
     }

Modified: grass/branches/develbranch_6/imagery/i.landsat.toar/main.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.landsat.toar/main.c	2013-10-23 21:16:40 UTC (rev 58093)
+++ grass/branches/develbranch_6/imagery/i.landsat.toar/main.c	2013-10-23 21:21:12 UTC (rev 58094)
@@ -3,7 +3,7 @@
  *
  * MODULE:       i.landsat.toar
  *
- * AUTHOR(S):    E. Jorge Tizado - ej.tizado at unileon.es
+ * AUTHOR(S):    E. Jorge Tizado - ej.tizado at unileon.es
  *               Hamish Bowman (small grassification cleanups)
  *               Yann Chemin (v7 + L5TM _MTL.txt support) [removed after update]
  *
@@ -32,7 +32,7 @@
     struct History history;
     struct GModule *module;
 
-    struct Cell_head cellhd;
+    struct Cell_head cellhd, orig_cellhd;
     char *mapset;
 
     void *inrast, *outrast;
@@ -42,10 +42,10 @@
 
     RASTER_MAP_TYPE in_data_type;
 
-    struct Option *input_prefix, *output_prefix, *metfn,
-	*sensor, *adate, *pdate, *elev, *bgain, *metho, *perc, *dark, *atmo;
+    struct Option *input_prefix, *output_prefix, *metfn, *sensor, *adate,
+	*pdate, *elev, *bgain, *metho, *perc, *dark, *atmo;
     char *inputname, *met, *outputname, *sensorname;
-    struct Flag *frad, *named;
+    struct Flag *frad, *print_meta, *named;
 
     lsat_data lsat;
     char band_in[GNAME_MAX], band_out[GNAME_MAX];
@@ -66,7 +66,7 @@
     module->description =
 	_("Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OLI");
     module->keywords =
-	_("imagery, landsat, top-of-atmosphere reflectance, dos-type simple atmospheric correction");
+	_("imagery, Landsat, radiance, reflectance, brightness temperature, atmospheric correction");
 
     /* It defines the different parameters */
     input_prefix = G_define_option();
@@ -182,7 +182,7 @@
     atmo->answer = "0.0";
     atmo->guisection = _("Settings");
 
-    /* It defines the different flags */
+    /* define the different flags */
     frad = G_define_flag();
     frad->key = 'r';
     frad->description =
@@ -193,6 +193,10 @@
     named->description =
 	_("Input raster maps use as extension the number of the band instead the code");
 
+    print_meta = G_define_flag();
+    print_meta->key = 'p';
+    print_meta->description = _("Print output metadata info");
+
     /* options and afters parser */
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
@@ -207,6 +211,7 @@
     outputname = output_prefix->answer;
     sensorname = sensor->answer ? sensor->answer : "";
 
+    G_get_window(&orig_cellhd);
     G_zero(&lsat, sizeof(lsat));
 
     if (adate->answer != NULL) {
@@ -238,8 +243,18 @@
     if (met != NULL) {
 	lsat.flag = METADATAFILE;
 	lsat_metadata(met, &lsat);
-	G_debug(1, "lsat.number = %d, lsat.sensor = [%s]", lsat.number,
-		lsat.sensor);
+	if (print_meta->answer) {
+	    G_message("Landsat-%d %s\n", lsat.number, lsat.sensor);
+	    G_message("Number of bands = %d\n", lsat.bands);
+	    G_message("Acquisition date = %s\n", lsat.date);
+	    G_message("Acquisition time (h) = %f\n", lsat.time);
+	    G_message("Sun elevation (degrees) = %f\n", lsat.sun_elev);
+	    G_message("Sun azimuth (degrees) = %f\n", lsat.sun_az);
+	    G_message("Product creation date = %s\n", lsat.creation);
+	    exit(EXIT_SUCCESS);
+	}
+	G_debug(1, "lsat.number = %d, lsat.sensor = [%s]",
+		lsat.number, lsat.sensor);
 
 	if (!lsat.sensor || lsat.number > 8 || lsat.number < 1)
 	    G_fatal_error(_("Failed to identify satellite"));
@@ -327,16 +342,17 @@
 		G_warning(_("Raster map <%s> not found"), band_in);
 		continue;
 	    }
-	    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);
 	    G_set_window(&cellhd);
+	    if ((infd = G_open_cell_old(band_in, "")) < 0)
+		G_fatal_error(_("Unable to open raster map <%s>"), band_in);
 
 	    in_data_type = G_raster_map_type(band_in, mapset);
 	    if (in_data_type < 0)
-	        G_fatal_error(_("Unable to read data type of raster map <%s>"), band_in);
+		G_fatal_error(_("Unable to read data type of raster map <%s>"),
+			      band_in);
 	    inrast = G_allocate_raster_buf(in_data_type);
 
 	    nrows = G_window_rows();
@@ -359,10 +375,9 @@
 			ptr = (void *)((DCELL *) inrast + col);
 			q = (int)*((DCELL *) ptr);
 			break;
-		   default:
-		        ptr = NULL;
-		        q = -1.;
-			break;
+		    default:
+			ptr = NULL;
+			q = -1.;
 		    }
 		    if (!G_is_null_value(ptr, in_data_type) &&
 			q >= lsat.band[i].qcalmin &&
@@ -411,13 +426,13 @@
 		lsat.sensor);
 	fprintf(stderr, " ACQUISITION DATE %s [production date %s]\n",
 		lsat.date, lsat.creation);
-	fprintf(stderr, "   earth-sun distance    = %.8lf\n", lsat.dist_es);
-	fprintf(stderr, "   solar elevation angle = %.8lf\n", lsat.sun_elev);
+	fprintf(stderr, "   Earth-sun distance    = %.8lf\n", lsat.dist_es);
+	fprintf(stderr, "   Solar elevation angle = %.8lf\n", lsat.sun_elev);
 	fprintf(stderr, "   Atmospheric correction: %s\n",
 		(method == UNCORRECTED ? "UNCORRECTED" : metho->answer));
 	if (method > DOS) {
 	    fprintf(stderr,
-		    "   percent of solar irradiance in path radiance = %.4lf\n",
+		    "   Percent of solar irradiance in path radiance = %.4lf\n",
 		    percent);
 	}
 	for (i = 0; i < lsat.bands; i++) {
@@ -473,17 +488,18 @@
 	    continue;
 	}
 
+	/* set same size as original band raster */
+	if (G_get_cellhd(band_in, mapset, &cellhd) < 0)
+	    G_fatal_error(_("Unable to read header of raster map <%s>"),
+			  band_in);
+	G_set_window(&cellhd);
 	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 (in_data_type < 0)
-	    G_fatal_error(_("Unable to read data type of 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>"),
+	    G_fatal_error(_("Unable to read data type of raster map <%s>"),
 			  band_in);
 
-	/* set same size as original band raster */
-	G_set_window(&cellhd);
 
 	/* controlling, if we can write the raster */
 	if (G_legal_filename(band_out) < 0)
@@ -500,11 +516,9 @@
 	ncols = G_window_cols();
 
 	G_important_message(_("Writing %s of <%s> to <%s>..."),
-			    (frad->
-			     answer ? _("radiance") : (lsat.band[i].
-						       thermal) ?
-			     _("temperature") : _("reflectance")), band_in,
-			    band_out);
+			    (frad->answer ? _("radiance")
+			     : (lsat.band[i].thermal) ? _("temperature") :
+			     _("reflectance")), band_in, band_out);
 	for (row = 0; row < nrows; row++) {
 	    G_percent(row, nrows, 2);
 
@@ -526,7 +540,6 @@
 		default:
 		    ptr = NULL;
 		    qcal = -1.;
-		    break;
 		}
 		if (G_is_null_value(ptr, in_data_type) ||
 		    qcal < lsat.band[i].qcalmin) {
@@ -584,56 +597,40 @@
 	/* Initialize the 'history' structure with basic info */
 	G_short_history(band_out, "raster", &history);
 	sprintf(history.edhist[0], " %s of Landsat-%d %s (method %s)",
-		(frad->
-		 answer ? "Radiance" : (lsat.band[i].
-					thermal ? "Temperature" :
-					"Reflectance")), lsat.number,
-		lsat.sensor, metho->answer);
+		(frad->answer ? "Radiance"
+		 : (lsat.band[i].thermal ? "Temperature" : "Reflectance")),
+		lsat.number, lsat.sensor, metho->answer);
 	sprintf(history.edhist[1],
 		"-----------------------------------------------------------------");
-	sprintf(history.edhist[2],
-		" Acquisition date (and time) ........... %s (%.4lf h)",
+	sprintf(history.edhist[2], " Acquisition date (and time) ........... %s (%.4lf h)",
 		lsat.date, lsat.time);
-	sprintf(history.edhist[3],
-		" Production date ....................... %s\n",
+	sprintf(history.edhist[3], " Production date ....................... %s\n",
 		lsat.creation);
-	sprintf(history.edhist[4],
-		" Earth-sun distance (d) ................ %.7lf",
+	sprintf(history.edhist[4], " Earth-sun distance (d) ................ %.7lf",
 		lsat.dist_es);
-	sprintf(history.edhist[5],
-		" Sun elevation (and azimuth) ........... %.5lf (%.5lf)",
+	sprintf(history.edhist[5], " Sun elevation (and azimuth) ........... %.5lf (%.5lf)",
 		lsat.sun_elev, lsat.sun_az);
-	sprintf(history.edhist[6],
-		" Digital number (DN) range ............. %.0lf to %.0lf",
+	sprintf(history.edhist[6], " Digital number (DN) range ............. %.0lf to %.0lf",
 		lsat.band[i].qcalmin, lsat.band[i].qcalmax);
-	sprintf(history.edhist[7],
-		" Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
+	sprintf(history.edhist[7], " Calibration constants (Lmin to Lmax) .. %+.5lf to %+.5lf",
 		lsat.band[i].lmin, lsat.band[i].lmax);
-	sprintf(history.edhist[8],
-		" DN to Radiance (gain and bias) ........ %+.5lf and %+.5lf",
+	sprintf(history.edhist[8], " 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[9],
-		    " Temperature (K1 and K2) ............... %.3lf and %.3lf",
+	    sprintf(history.edhist[9], " Temperature (K1 and K2) ............... %.3lf and %.3lf",
 		    lsat.band[i].K1, lsat.band[i].K2);
 	    history.edlinecnt = 10;
 	}
 	else {
 	    sprintf(history.edhist[9],
-		    " Mean solar irradiance (ESUN) .......... %.3lf",
-		    lsat.band[i].esun);
+		    " Mean solar irradiance (ESUN) .......... %.3lf", lsat.band[i].esun);
 	    sprintf(history.edhist[10],
-		    " Radiance to Reflectance (divide by) ... %+.5lf",
-		    lsat.band[i].K1);
+		    " Radiance to Reflectance (divide by) ... %+.5lf", lsat.band[i].K1);
 	    history.edlinecnt = 11;
 	    if (method > DOS) {
 		sprintf(history.edhist[11], " ");
-		sprintf(history.edhist[12],
-			" Dark object (%4d pixels) DN = ........ %d", pixel,
-			dn_dark[i]);
-		sprintf(history.edhist[13],
-			" Mode in reflectance histogram ......... %.5lf",
-			ref_mode);
+		sprintf(history.edhist[12], " Dark object (%4d pixels) DN = ........ %d", pixel, dn_dark[i]);
+		sprintf(history.edhist[13], " Mode in reflectance histogram ......... %.5lf", ref_mode);
 		history.edlinecnt = 14;
 	    }
 	}
@@ -653,6 +650,7 @@
 
 	/*  set raster timestamp from acq date? (see r.timestamp module)  */
     }
+    G_set_window(&orig_cellhd);
 
     exit(EXIT_SUCCESS);
 }



More information about the grass-commit mailing list