-
Notifications
You must be signed in to change notification settings - Fork 13
/
i.landsat8.swlst.html
267 lines (264 loc) · 22.2 KB
/
i.landsat8.swlst.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
<h2 id="description">DESCRIPTION</h2>
<p><em>i.landsat8.swlst</em> is an implementation of a robust and practical Slit-Window (SW) algorithm estimating land surface temperature (LST), from the Thermal Infra-Red Sensor (TIRS) aboard Landsat 8 with an accuracy of better than 1.0 K. [1]</p>
<h3 id="overview">Overview</h3>
<p>The components of the algorithm estimating LST values are at-satellite brightness temperature (BT); land surface emissivity (LSE); and the coefficients of the main Split-Window equation (SWC).</p>
<p>LSEs are derived from an established look-up table linking the FROM-GLC classification scheme to average emissivities. The NDVI and the FVC are <em>not</em> computed each time an LST estimation is requested. Read [0] for details.</p>
<p>The SWC depend on each pixel's column water vapor (CWV). CWV values are retrieved based on a modified Split-Window Covariance-Variance Matrix Ratio method (MSWCVMR) [1, 2]. <strong>Note</strong>, the spatial discontinuity found in the images of the retrieved CWV, is attributed to the data gap in the images caused by stray light outside of the FOV of the TIRS instrument [2]. In addition, the size of the spatial window querying for CWV values in adjacent pixels, is a key parameter of the MSWCVMR method. It influences accuracy and performance.</p>
<p>At-satellite brightness temperatures are derived from the TIRS channels 10 and 11. Prior to any processing, these are filtered for clouds and their quantized digital numbers converted to at-satellite temperature values.</p>
<pre><code> +--------+ +--------------------------+
|Landsat8+--->Cloud screen & calibration|
+--------+ +---+--------+-------------+
| |
| |
+-v-+ +--v-+
|OLI| |TIRS|
+-+-+ +--+-+
| |
| |
+--v-+ +--v-------------------+ +-------------+
|NDVI| |Brightness temperature+-->MSWCVM method|
+----------+ +--+-+ +--+-------------------+ +----------+--+
|Land cover| | | |
+----------+ | | |
| +-v-+ +--v-------------------+ +------v--+
| |FVC| |Split Window Algorithm| |ColWatVap|
+---------------------v--+ +-+-+ +-------------------+--+ +------+--+
|Emissivity look|up table| | | |
+---------------------+--+ | | |
| +--v--------------------+ | +---------v--+
+------>Pixel emissivity ei, ej+--> | <--+Coefficients|
+-----------------------+ | +------------+
|
|
+---------------v--+
|LST and emissivity|
+------------------+</code></pre>
<pre><code> [ Figure 3 in [0]: Flowchart of retrieving LST from Landsat8 ]</code></pre>
<p>Hence, to produce an LST map, the algorithm requires at minimum:</p>
<ul>
<li>TIRS bands 10 and 11</li>
<li>the acquisition's metadata file (MTL)</li>
<li>a Finer Resolution Observation & Monitoring of Global Land Cover (FROM-GLC) product</li>
</ul>
<h3 id="details">Details</h3>
<p>A new refinement of the generalized split-window algorithm proposed by Wan (2014) [19] is added with a quadratic term of the difference amongst the brightness temperatures (Ti, Tj) of the adjacent thermal infrared channels, which can be expressed as (equation 2 in [0])</p>
<p><code>LST = b0 + [b1 + b2 * (1-e)/e + b3 * (De/e2)] * (Ti+T)/j2 + [b4 + b5 * (1-e)/e + b6 * (De/e2)] * (Ti-Tj)/2 + b7 * (Ti-Tj)^2</code></p>
<p>where:</p>
<ul>
<li><code>Ti</code> and <code>Tj</code> are Top of Atmosphere brightness temperatures measured in channels <code>i</code> (~11.0 microns) and <code>j</code> (~12.0 µm), respectively</li>
<li>from http://landsat.usgs.gov/band_designations_landsat_satellites.php:
<ul>
<li>Band 10, Thermal Infrared (TIRS) 1, 10.60-11.19, 100*(30)</li>
<li>Band 11, Thermal Infrared (TIRS) 2, 11.50-12.51, 100*(30)</li>
</ul></li>
<li>e is the average emissivity of the two channels (i.e., <code>e = 0.5 [ei + ej]</code>)</li>
<li>De is the channel emissivity difference (i.e., <code>De = ei - ej</code>)</li>
<li><code>bk</code> (k = 0, 1, ... 7) are the algorithm coefficients derived from a simulated dataset.</li>
</ul>
<blockquote>
<p>Note, however, <strong>the last quadratic term</strong> of the
Split-Window equation <strong>is applied only over barren land</strong>. [Reference
Required!]</p>
</blockquote>
<p>In the above equations,</p>
<ul>
<li><code>dk</code> (k = 0, 1...6) and <code>ek</code> (k = 1, 2, 3, 4) are the algorithm coefficients;</li>
<li><code>w</code> is the Column Water Vapor;</li>
<li><code>e</code> and <code>De</code> are the average emissivity and emissivity difference of two adjacent thermal channels, respectively, which are similar to Equation (2);</li>
<li>and <code>fk</code> (k = 0 and 1) is related to the influence of the atmospheric transmittance and emissivity, i.e., <code>fk = f(ei, ej, ti, tji)</code>.</li>
</ul>
<h3 id="comparing-to-other-split-window-algorithms">Comparing to other split-window algorithms</h3>
<p>From the paper:</p>
<blockquote>
<p>Note that the algorithm (Equation (6a)) proposed by Jimenez-Munoz et al. added column water vapor (CWV) directly to estimate LST. Rozenstein et al. used CWV to estimate the atmospheric transmittance (<code>ti</code>, <code>tj</code>) and optimize retrieval accuracy explicitly. Therefore, if the atmospheric CWV is unknown or cannot be obtained successfully, neither of the two algorithms in Equations (6a) and (6b) will work. By contrast, although the current algorithm also needs CWV to determine the coefficients, it still works for unknown CWVs because the coefficients are obtained regardless of the CWV, as shown in Table 1 [0].</p>
</blockquote>
<h2 id="notes">NOTES</h2>
<h3 id="cloud-masking">Cloud Masking</h3>
<p>The first important step of the algorithm is cloud screening. The module offers two ways to achieve this:</p>
<ol style="list-style-type: decimal">
<li>use of the Quality Assessment band and some user-defined QA pixel value</li>
<li>use an external cloud map as an inverted MASK</li>
</ol>
<h3 id="calibration-of-tirs-channels-10-11">Calibration of TIRS channels 10, 11</h3>
<h4 id="conversion-to-spectral-radiance">Conversion to Spectral Radiance</h4>
<p>Conversion of Digital Numbers to TOA Radiance. OLI and TIRS band data can be converted to TOA spectral radiance using the radiance rescaling factors provided in the metadata file:</p>
<p><code>Ll = ML * Qcal + AL</code></p>
<p>where:</p>
<ul>
<li><code>Ll</code> = TOA spectral radiance (Watts/( m2 * srad * microns))</li>
<li><code>ML</code> = Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number)</li>
<li><code>AL</code> = Band-specific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number)</li>
<li><code>Qcal</code> = Quantized and calibrated standard product pixel values (DN)</li>
</ul>
<h4 id="conversion-to-at-satellite-temperature">Conversion to at-Satellite Temperature</h4>
<p>Conversion to At-Satellite Brightness Temperature TIRS band data can be converted from spectral radiance to brightness temperature using the thermal constants provided in the metadata file:</p>
<p><code>T = K2 / ln((K1/Ll) + 1)</code></p>
<p>where:</p>
<ul>
<li><code>T</code> = At-satellite brightness temperature (K)</li>
<li><code>Ll</code> = TOA spectral radiance (Watts/(m^2 * srad * microns)), below 'DUMMY_RADIANCE'</li>
<li><code>K1</code> = Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the band number, 10 or 11)</li>
<li><code>K2</code> = Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the band number, 10 or 11)</li>
</ul>
<div class="figure">
<p><img src="images/LC81840332014146LGN00_B10.jpg"> <img src="images/LC81840332014146LGN00_B11.jpg"></p>
</div>
<h3 id="land-surface-emissivity">Land Surface Emissivity</h3>
<p>Determination of LSEs (overview of Section 3.2)</p>
<ol style="list-style-type: decimal">
<li><p>The FROM-GLC (30m) contains 10 types of land covers (cropland, forest, grassland, shrubland, wetland, waterbody, tundra, impervious, barren land and snow-ice).</p></li>
<li><p>Deriving emissivities for each land cover class by using different combinations of three BRDF kernel models (geometrical, volumetric and specular models)</p></li>
<li><p>Vegetation and ground emissivity spectra for the BRDF models selected from the MODIS University of California, Santa Barbara (UCSB) Emissivity Library</p></li>
<li><p>Estimating FVC (to obtain emissivity of land cover with temporal variation) from NDVI based on Carlson (1997) and Sobrino (2001)</p></li>
<li><p>Finally, establishing the average emissivity Look-Up table</p></li>
</ol>
<h3 id="column-water-vapor">Column Water Vapor</h3>
<p>Retrieving atmospheric CWV from Landsat8 TIRS data based on the modified split-window covariance and variance ratio (MSWCVR).</p>
<p>Algorithm Coefficients (overview of Section 3.1)</p>
<ol style="list-style-type: decimal">
<li><p>The CWV is divided into 5 sub-ranges with an overlap of 0.5 g/cm2 between 2 adjacent sub-ranges: [0.0, 2.5], [2.0, 3.5], [3.0, 4.5], [4.0, 5.5] and [5.0, 6.3] g/cm2.</p></li>
<li><p>The CWV is retrieved from a modified split-window covariance and variance ratio method.</p></li>
<li><p>However, given the somewhat unsuccessful CWV retrieval, a group of coefficients for the entire CWV range is calculated to ensure the spatial continuity of the LST product.</p></li>
</ol>
<h4 id="modified-split-window-covariance-variance-method">Modified Split-Window Covariance-Variance Method</h4>
<p>With a vital assumption that the atmosphere is unchanged over the neighboring pixels, the MSWCVR method relates the atmospheric CWV to the ratio of the upward transmittances in two thermal infrared bands, whereas the transmittance ratio can be calculated based on the TOA brightness temperatures of the two bands. Considering N adjacent pixels, the CWV in the MSWCVR method is estimated as:</p>
<ul>
<li><code>cwv = c0 + c1*(tj/ti) + c2*(tj/ti)^2</code> (3a)</li>
</ul>
<p>where:</p>
<ul>
<li><code>tj/ti</code> ~ <code>Rji = SUM [(Tik-Ti\_mean) \* (Tjk-Tj\_mean)] / SUM[(Tik-Tj\_mean)\^2]</code></li>
</ul>
<p>In Equation (3a):</p>
<ul>
<li><code>c0</code>, <code>c1</code> and <code>c2</code> are coefficients obtained from simulated data;</li>
<li><code>t</code> is the band effective atmospheric transmittance;</li>
<li><code>N</code> is the number of adjacent pixels (excluding water and cloud pixels) in a spatial window of size <code>n</code> (i.e., <code>N = n x n</code>);</li>
<li><code>Ti,k</code> and <code>Tj,k</code> are top of atmosphere brightness temperatures (K) of bands <code>i</code> and <code>j</code> for the <code>k</code>th pixel;</li>
<li><code>mean(Ti)</code> and <code>mean(Tj)</code> are the mean (or median -- not implemented yet) brightness temperatures of the <code>N</code> pixels for the two bands.</li>
</ul>
<p>TIRS channels are originally of 100m spatial resolution. However, bands 10 and 11 are resampled, via a cubic convolution filter, to 30m. Consequently, an appropriately sized spatial window is required for a meaningful CWV estimation attempt. The spatial window should be composed by a number of pixels stretching over an area that accounts for several adjacent <em>100m</em>-sized pixels. <strong>Note</strong>, while the CWV estimation accuracy increases with larger windows (up to a certain level), the performance (speed) of the module decreases greatly.</p>
<p>The regression coefficients:</p>
<ul>
<li><code>c0</code> = -9.674</li>
<li><code>c1</code> = 0.653</li>
<li><code>c2</code> = 9.087</li>
</ul>
<p>where obtained by:</p>
<ul>
<li>946 cloud-free TIGR atmospheric profiles,</li>
<li>the new high accurate atmospheric radiative transfer model MODTRAN 5.2</li>
<li>simulating the band effective atmospheric transmittance Model analysis indicated that this method will obtain a CWV RMSE of about 0.5 g/cm^2.</li>
</ul>
<p>The algorithm will not cause significant uncertainty to the final LST retrieval with known CWV, but it will lead some error to the LST result for the cases without input CWV. To reduce this effect, the authors are trying to find more representative profiles to optimize the current algorithm.</p>
<p>Details about the columnw water vapor retrieval can be found in:</p>
<p>Ren, H.; Du, C.; Qin, Q.; Liu, R.; Meng, J.; Li, J. Atmospheric water vapor retrieval from landsat 8 and its validation. In Proceedings of the IEEE International Geosciene and Remote Sensing Symposium (IGARSS), Quebec, QC, Canada, July 2014; pp. 3045--3048.</p>
<h3 id="split-window-algorithm">Split-Window Algorithm</h3>
<p>The algorithm removes the atmospheric effect through differential atmospheric absorption in the two adjacent thermal infrared channels centered at about 11 and 12 microns. The linear or non-linear combination of the brightness temperatures is finally applied for LST estimation based on the equation:</p>
<p><code>LST = b0 + + (b1 + b2 \* ((1-ae)/ae)) + + b3 \* (de/ae) \* ((t10 + t11)/2) + + (b4 + b5 \* ((1-ae)/ae) + b6 \* (de/ae\^2)) \* ((t10 - t11)/2) + + b7 \* (t10 - t11)\^2</code></p>
<p>To reduce the influence of the CWV error on the LST, for a CWV within the overlap of two adjacent CWV sub-ranges, we first use the coefficients from the two adjacent CWV sub-ranges to calculate the two initial temperatures and then use the average of the initial temperatures as the pixel LST.</p>
<p>For example, the LST pixel with a CWV of 2.1 g/cm2 is estimated by using the coefficients of [0.0, 2.5] and [2.0, 3.5]. This process initially reduces the <strong>delta-</strong>LSTinc and improves the spatial continuity of the LST product.</p>
<h2 id="examples">EXAMPLES</h2>
<p>At minimum, the module requires the following in order to derive a land surface temperature map:</p>
<ol style="list-style-type: decimal">
<li><p>The Landsat8 scene's acquisition metadata (MTL file)</p></li>
<li><p>Bands 10, 11 and QA</p></li>
<li><p>A FROM-GLC product for the same Path and Row as the Landsat scene to be processed</p></li>
</ol>
The shortest call for processing a complete Landsat8 scene normally is:
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC -n</code></pre>
</div>
<p>where:</p>
<ul>
<li><p><strong><code>mtl=</code></strong> the name of the MTL metadata file (normally with a <code>.txt</code> extension)</p></li>
<li><p><strong><code>prefix=</code></strong> the prefix of the band names imported in GRASS GIS' data base</p></li>
<li><p><strong><code>landcover=</code></strong> the name of the FROM-GLC map that covers the extent of the Landsat8 scene under processing</p></li>
<li><p>the <strong><code>n</code></strong> flag will set zero digital number values, which may represent NoData in the original bands, to NULL. This option is probably unnecessary for smaller regions in which there are no NoData pixels present.</p></li>
</ul>
<p>The pixel value 61440 is selected automatically to build a cloud mask. At the moment, only a single pixel value may be requested from the Quality Assessment band. For details, refer to [http://landsat.usgs.gov/L8QualityAssessmentBand.php USGS' webpage for Landsat8 Quality Assessment Band]</p>
<p><strong><code>window</code></strong> is an important option. It defines the size of the spatial window querying for column water vapor values. Small window sizes introduce a spatial discontinuation effect in the final LST image. Larger window sizes lead to more accurate results, at the cost of performance. However, too large window sizes should be avoided as they would include large variations of land and atmospheric conditions. In [2] it is stated:</p>
<blockquote>
<p>A small window size n (N = n * n, see equation (1a)) cannot ensure a high correlation between two bands' temperatures due to the instrument noise. In contrast, the size cannot be too large because the variations in the surface and atmospheric conditions become larger as the size increases.</p>
</blockquote>
<p>An example instructing a spatial window of size 9^2 is:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC window=9</code></pre>
</div>
<p>In order to restrict the processing in to the currently set computational region, the <strong><code>-k</code></strong> flag can be used:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC -k </code></pre>
</div>
<p>The Landsat8 scene's time and date of acquisition may be applied to the LST (and to the optionally requested CWV) map via the <strong><code>-t</code></strong> flag.</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC -k -t</code></pre>
</div>
<p>The output land surface temperature map maybe be delivered in Celsius degrees (units and appropriate color table) via the <strong><code>-c</code></strong> flag:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC -k -c</code></pre>
</div>
<div class="figure">
<div class="figure">
<img src="images/lst_window_5_detail_celsius_500px.jpg">
</div>
</div>
<p>A user defined map for clouds, instead of relying on the Quality Assessment band, can be used via the <strong><code>clouds</code></strong> option:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC clouds=Cloud_Map -k</code></pre>
</div>
<p>Using the <strong><code>prefix_bt</code></strong> option, the in-between at-satellite brightness temperature maps may be saved for re-use in sub-sequent trials via the <strong><code>t10</code></strong> and <strong><code>t11</code></strong> options. Using the <code>t10</code> and <code>t11</code> options, will skip the conversion from digital numbers for bands B10 and B11. Alternatively, any existing at-satellite brightness temperature maps maybe used via the <code>t10/11</code> options. For example using the <code>t11</code> option instead of <code>b11</code>:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL b10=B10 t11=AtSatellite_Temperature_11 landcover=FROM_GLC -k</code></pre>
</div>
<p>or using both <code>t10</code> and <code>t11</code>:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL t10=AtSatellite_Temperature_10 t11=AtSatellite_Temperature_11 landcover=FROM_GLC -k</code></pre>
</div>
<p>A faster run is achieved by using existing maps for all in-between processing steps: at-satellite temperatures, cloud and emissivity maps.</p>
<pre><code>* At-satellite temperature maps (optiones `t10`, `t11`) may be derived via
the i.landsat.toar module. Note that `i.landsat.toar` does not
process single bands selectively.
* The `clouds` option can be any user defined map. Essentialy, it applies
the given map as an inverted mask.
* The emissivity maps, derived by the module itself, can be saved once via
the `emissivity_out` and `delta_emissivity_out` options and used
afterwards via the **`emissivity`** and **`delta_emissivity`** options.
Expert users, however, may use emissivity maps from other sources
directly. An example command may be:</code></pre>
<div class="code">
<pre><code>i.landsat8.swlst t10=BT10 t11=BT11 clouds=Cloud_Map emissivity=Average_Emissivity_Map delta_emissivity=Delta_Emissivity_Map landcover=FROM_GLC -n</code></pre>
</div>
<p>Expert users may need to request for a <em>fixed</em> average surface emissivity, in order to perform the algorithm for a single land cover class (one from the classes defined in the FROM-GLC classification scheme) via the <code>emissivity_class</code> option. Consequently, <code>emissivity_class</code> cannot be used at the same time with the <code>landover</code> option.</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL b10=B10 t11=AtSatellite_Temperature_11 qab=BQA emissivity_class="Croplands" -c </code></pre>
</div>
<p>A <em>transparent</em> run-through of <em>what kind of</em> and <em>how</em> the module performs its computations, may be requested via the use of both the <strong><code>--v</code></strong> and <strong><code>-i</code></strong> flags:</p>
<div class="code">
<pre><code>i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC -i --v </code></pre>
</div>
<p>The above will print out a description for each individual processing step, as well as the actual mathematical epxressions applied via GRASS GIS' <code>r.mapcalc</code> module.</p>
<h3 id="example-figures">Example figures</h3>
<div class="figure">
<p><img src="images/lst_window_7.jpg"> <img src="images/lst_window_9.jpg"><img src="images/lst_window_11.jpg"></p>
</div>
<div class="figure">
<p><img src="images/lst_window_7_detail_500px.jpg"> <img src="images/lst_window_9_detail_500px.jpg"> <img src="images/lst_window_11_detail_500px.jpg"></p>
</div>
<h2 id="todo">TODO</h2>
<ul>
<li>Go through <a href="http://trac.osgeo.org/grass/wiki/Submitting/Python">Submitting Python</a></li>
<li>Proper command history tracking.</li>
<li>Deduplicate code where applicable</li>
<li>Test compiling in other systems</li>
<li>Improve documentation</li>
</ul>
<h2 id="references">REFERENCES</h2>
<ul>
<li><p>[0] Du, Chen; Ren, Huazhong; Qin, Qiming; Meng, Jinjie; Zhao, Shaohua. 2015. "A Practical Split-Window Algorithm for Estimating Land Surface Temperature from Landsat 8 Data." Remote Sens. 7, no. 1: 647-665. http://www.mdpi.com/2072-4292/7/1/647/htm</p></li>
<li><p>[1] Huazhong Ren, Chen Du, Qiming Qin, Rongyuan Liu, Jinjie Meng, and Jing Li. "Atmospheric Water Vapor Retrieval from Landsat 8 and Its Validation." 3045--3048. IEEE, 2014.</p></li>
<li><p>[2] Ren, H., Du, C., Liu, R., Qin, Q., Yan, G., Li, Z. L., & Meng, J. (2015). Atmospheric water vapor retrieval from Landsat 8 thermal infrared images. Journal of Geophysical Research: Atmospheres, 120(5), 1723-1738.</p></li>
</ul>
<h2 id="see-also">SEE ALSO</h2>
<p><em><a href="i.emissivity.html">i.emissivity</a></em></p>
<h2 id="authors">AUTHORS</h2>
<p>Nikos Alexandris</p>