i.landsat.toar: fixes for Landsat-8 metadata file support (author: E. Jorge Tizado); sync with https://trac.osgeo.org/grass/changeset/58093

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@58097 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.landsat.toar/i.landsat.toar.html

@@ -97,13 +97,66 @@ The values are,
     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 @@ i.landsat.toar input_prefix=203_30. output_prefix=_toar \
   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>

imagery/i.landsat.toar/landsat.c

@@ -87,8 +87,7 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 		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 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 		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 @@ void lsat_bandctes(lsat_data * lsat, int i, char method,
 	    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;
     }
 }

imagery/i.landsat.toar/landsat_met.c

@@ -35,6 +35,7 @@ static void date_replace_slash(char *str)
     }
 }
 
+/* OLD Metadata Files */
 void get_metformat(const char metadata[], char *key, char value[])
 {
     int i = 0;
@@ -56,6 +57,17 @@ void get_metformat(const char metadata[], char *key, char value[])
     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,6 +85,16 @@ void get_mtlformat(const char metadata[], char *key, char value[])
     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 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
     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 lsat_metadata(char *metafile, lsat_data * lsat)
     (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 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
 	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 @@ void lsat_metadata(char *metafile, lsat_data * lsat)
 			}
 		    }
 		}
-		G_warning
-		    ("ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file");
+		G_warning("ESUN evaluated from REFLECTANCE_ADDITIVE_FACTOR_BAND");
 	    }
 	}
     }

imagery/i.landsat.toar/main.c

@@ -224,7 +224,6 @@ int main(int argc, char *argv[])
     named->description =
 	_("Input raster maps use as extension the number of the band instead the code");
 
-    /* define the different flags */
     print_meta = G_define_flag();
     print_meta->key = 'p';
     print_meta->description = _("Print output metadata info");
@@ -429,10 +428,9 @@ int main(int argc, char *argv[])
 			ptr = (void *)((DCELL *) inrast + col);
 			q = (int)*((DCELL *) ptr);
 			break;
-		   default:
-		        ptr = NULL;
-		        q = -1.;
-			break;
+		    default:
+			ptr = NULL;
+			q = -1.;
 		    }
 		    if (!Rast_is_null_value(ptr, in_data_type) &&
 			q >= lsat.band[i].qcalmin &&
@@ -599,7 +597,6 @@ int main(int argc, char *argv[])
 		default:
 		    ptr = NULL;
 		    qcal = -1.;
-		    break;
 		}
 		if (Rast_is_null_value(ptr, in_data_type) ||
 		    qcal < lsat.band[i].qcalmin) {
@@ -635,8 +632,8 @@ int main(int argc, char *argv[])
 
 
 	G_free(inrast);
-	G_free(outrast);
 	Rast_close(infd);
+	G_free(outrast);
 	Rast_close(outfd);