[GRASS-SVN] r68288 - in grass-addons/grass7/imagery: . i.landsat8.swlst
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Apr 20 06:47:02 PDT 2016
Author: nikosa
Date: 2016-04-20 06:47:02 -0700 (Wed, 20 Apr 2016)
New Revision: 68288
Added:
grass-addons/grass7/imagery/i.landsat8.swlst/
grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B10.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B10_B11_absolute_differences.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B11.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_detail_hpfa_pansharpened.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_r4_g3_b2.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/Makefile
grass-addons/grass7/imagery/i.landsat8.swlst/README.md
grass-addons/grass7/imagery/i.landsat8.swlst/average_emissivity.csv
grass-addons/grass7/imagery/i.landsat8.swlst/column_water_vapor.py
grass-addons/grass7/imagery/i.landsat8.swlst/csv_to_dictionary.py
grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients.csv
grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients_complete_range.csv
grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html
grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py
grass-addons/grass7/imagery/i.landsat8.swlst/landsat8_mtl.py
grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_11_detail_500px.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_13_detail_500px.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_7_detail_500px.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_9_detail_500px.jpg
grass-addons/grass7/imagery/i.landsat8.swlst/map.png
grass-addons/grass7/imagery/i.landsat8.swlst/mtl.txt
grass-addons/grass7/imagery/i.landsat8.swlst/split_window_lst.py
grass-addons/grass7/imagery/i.landsat8.swlst/test_column_water_vapor.py
grass-addons/grass7/imagery/i.landsat8.swlst/test_landsat8_mtl.py
grass-addons/grass7/imagery/i.landsat8.swlst/test_split_window_lst.py
Log:
Practical split-window algorithm estimating Land Surface Temperature from Landsat 8 OLI/TIRS imagery
Added: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B10.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B10.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B10_B11_absolute_differences.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B10_B11_absolute_differences.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B11.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_B11.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_detail_hpfa_pansharpened.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_detail_hpfa_pansharpened.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_r4_g3_b2.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/LC81840332014146LGN00_r4_g3_b2.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/Makefile 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.landsat8.swlst
+
+ETCFILES = landsat8_mtl split_window_lst column_water_vapor csv_to_dictionary
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/Python.make
+
+default: script
Added: grass-addons/grass7/imagery/i.landsat8.swlst/README.md
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/README.md (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/README.md 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,212 @@
+*i.landsat8.swlst* is a GRASS GIS add-on, implementating a 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.
+
+
+Quick examples
+==============
+
+After installation (see section *Installation* below), from within a GRASS-GIS
+session, retrieve usage details via `i.landsat8.swlst --help`
+
+The shortest call for processing a complete Landsat8 scene normally is:
+
+```bash
+i.landsat8.swlst mtl=MTL prefix=B landcover=FROM_GLC
+```
+
+where:
+
+- `mtl=` the name of the MTL metadata file (normally with a `.txt` extension)
+- `prefix=` the prefix of the band names imported in GRASS GIS' data base
+- `landcover=` the name of a FROM-GLC map product that covers the extent of the
+ Landsat8 scene under processing. FORM-GLC products are available at
+ <http://data.ess.tsinghua.edu.cn/>.
+
+A faster call is to use existing maps for all in-between
+processing steps: at-satellite temperatures, cloud and emissivity maps.
+
+ * At-satellite temperature maps (options `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:
+
+```bash
+ i.landsat8.swlst t10=T10 t11=T11 clouds=Cloud_Map emissivity=Average_Emissivity_Map delta_emissivity=Delta_Emissivity_Map landcover=FROM_GLC -k -c
+```
+
+For details and more examples, read the manual.
+
+
+Description
+===========
+
+The components of the algorithm estimating LST values are at-satellite
+brightness temperature (BT); land surface emissivities (LSEs); and the
+coefficients of the main Split-Window equation (SWCs).
+
+**LSEs** are derived from an established look-up table linking the FROM-GLC
+classification scheme to average emissivities. The NDVI and the FVC are *not*
+computed each time an LST estimation is requested. Read [0] for details.
+
+The **SWCs** 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]. **Note**, 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. In [2]
+it is stated:
+
+> 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.
+
+At-satellite brightness temperatures are derived from the TIRS channels 10 and
+11. Prior to any processing, the raw digital numbers are filtered for clouds.
+
+To produce an LST map, the algorithm requires at minimum:
+
+- TIRS bands 10 and 11
+- the acquisition's metadata file (MTL)
+- a Finer Resolution Observation & Monitoring of Global Land Cover (FROM-GLC)
+ product
+
+Installation
+============
+
+## Requirements
+---------------
+
+see [GRASS Addons SVN repository, README file, Installation - Code
+Compilation](https://svn.osgeo.org/grass/grass-addons/README)
+
+## Steps
+
+Making the script `i.lansat8.swlst` available from within any GRASS-GIS ver.
+7.x session, may be done via the following steps:
+
+1. launch a GRASS-GIS’ ver. 7.x session
+
+2. navigate into the script’s source directory
+
+3. execute `make MODULE_TOPDIR=$GISBASE`
+
+
+
+Implementation notes
+====================
+
+- Created on Wed Mar 18 10:00:53 2015
+- First all-through execution: Tue May 12 21:50:42 EEST 2015
+
+
+## To Do
+
+[High Priority]
+
+- ~~Fix retrieval of adjacent subranges (exclude 6, get it only if there is no
+ match for subranges 1, 2, 3, 4, and 5)~~
+
+- Evaluate BIG mapcalc expressions -- are they correct? I guess so ;-)
+
+ - ~~Expression for Column Water Vapor~~
+ - ~~CWV output values range -- is it rational?~~ It was not. There is a
+ typo in paper [0]. The correct order of the coefficients is in papers [1,
+ 2].
+ - ~~Expression for Land Surface Temperature~~
+ - ~~LST output values range -- is it rational? At the moment, not!~~
+ Fixed. The main Split-Window equation was wrong.
+
+- ~~Why is the LST out of range when using a fixed land cover class?~~ Cloudy
+ pixels are, mainly, the reason. Better cloud masking is the solution.
+- ~~Why does the multi-step approach on deriving the CWV map differ from the
+ single big mapcalc expression?~~ **Fixed**
+- ~~Implement direct conversion of B10, B11 to brightness temperature values.~~
+ **Done**
+- ~~Get the FROM-GLC map,~~ **Found**
+- ~~implement mechanism to read land cover classes from it
+ and use'm to retrieve emissivities~~ **Done**
+- ~~How to use the FVC?~~ Don't. Just **use the Look-up table** (see [\*] for
+ details).
+- ~~Save average emissivity and delta emissivity maps for caching (re-use in
+ subsequent trials, huge time saver!)~~ **Implemented**
+
+[Mid]
+
+- ~~Redo the example screenshots for the manual after corrections for CWV, LST
+ equations.~~
+- Use existing i.emissivity? Not exactly compatible -- read paper for details.
+ Anyhow, options to input average and delta emissivity maps implemented.
+- Raster Row I/O -- Maybe *not* an option: see discussion with Pietro Zambelli
+- How to perform pixel value validity checks for in-between and end products?
+ `r.mapcalc` can't do this. Best to implement a test checking the values
+ on-the-fly while they are created. A C-function?
+
+[Low]
+
+- Deduplicate code in `split_window_lst` class, in functions
+`_build_average_emissivity_mapcalc()` and
+`_build_delta_emissivity_mapcalc()`
+- Implement a median window filter, as another option in addition to mean.
+- Profiling
+- Implement a complete cloud masking function using the BQA image. Support for
+ user requested confidence or types of clouds (?). Eg: options=
+ clouds,cirrus,high,low ?
+- Multi-Threading? Note, r.mapcalc is already.
+
+[\*] Details: the authors followed the CBEM method. Based on the FROM-GLC map,
+they derived the following look-up table (LUT):
+
+```
+Emissivity Class|TIRS10|TIRS11
+Cropland|0.971|0.968
+Forest|0.995|0.996
+Grasslands|0.97|0.971
+Shrublands|0.969|0.97
+Wetlands|0.992|0.998
+Waterbodies|0.992|0.998
+Tundra|0.98|0.984
+Impervious|0.973|0.981
+Barren Land|0.969|0.978
+Snow and ice|0.992|0.998
+```
+
+References
+==========
+
+- [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\#sthash.ba1pt9hj.dpuf
+
+- [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.
+
+- [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.
+
+- [3] FROM-GLC products, <http://data.ess.tsinghua.edu.cn/>
+
+Ευχαριστώ
+=========
+
+- Yann Chemin
+- Pietro Zambelli
+- StackExchange contributors
+- Sajid Pareeth
+- Georgios Alexandris, Anthoula Alexandri
+- Special thanks to author Huazhong Ren for commenting on questions (personal
+communication)
Added: grass-addons/grass7/imagery/i.landsat8.swlst/average_emissivity.csv
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/average_emissivity.csv (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/average_emissivity.csv 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,11 @@
+Emissivity Class|TIRS10|TIRS11|Class Code
+Cropland|0.971|0.968|(10, 11, 12, 13)
+Forest|0.995|0.996|(20, 21, 22, 23, 24)
+Grasslands|0.97|0.971|(30, 31, 32, 51, 72)
+Shrublands|0.969|0.97|(40, 71)
+Wetlands|0.992|0.998|
+Waterbodies|0.992|0.998|(50, 60, 61, 62, 63)
+Tundra|0.98|0.984|
+Impervious|0.973|0.981|(80, 81, 82)
+Barren Land|0.969|0.978|(90, 52, 91, 92, 93, 94, 95, 96)
+Snow and ice|0.992|0.998|(100, 101, 102)
Added: grass-addons/grass7/imagery/i.landsat8.swlst/column_water_vapor.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/column_water_vapor.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/column_water_vapor.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,439 @@
+#!/usr/bin/python\<nl>\
+# -*- coding: utf-8 -*-
+
+"""
+Determining atmospheric column water vapor based on
+Huazhong Ren, Chen Du, Qiming Qin, Rongyuan Liu, Jinjie Meng, Jing Li
+
+ at author nik | 2015-04-18 03:48:20
+"""
+
+# globals
+DUMMY_Ti_MEAN = 'Mean_Ti'
+DUMMY_Tj_MEAN = 'Mean_Tj'
+DUMMY_Rji = 'Ratio_ji'
+
+
+# helper functions
+def random_adjacent_pixel_values(pixel_modifiers):
+ """
+ """
+ import random
+ return [random.randint(250, 350) for dummy_idx in
+ range(len(pixel_modifiers))]
+
+
+class Column_Water_Vapor():
+ """
+ Retrieving atmospheric column water vapor from Landsat8 TIRS data based on
+ the modified split-window covariance and variance ratio (MSWCVR).
+
+ -------------------------------------------------------------------------
+ *Note,* this class produces valid expressions for GRASS GIS' mapcalc raster
+ processing module and does not directly compute column water vapor
+ estimations.
+ -------------------------------------------------------------------------
+
+ 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:
+
+ - cwv = c0 + c1 * (tj / ti) + c2 * (tj / ti)^2
+
+ - tj/ti ~ Rji = SUM [ ( Tik - Ti_mean ) * ( Tjk - Tj_mean ) ] /
+ SUM [ ( Tik - Tj_mean )^2 ]
+
+ In Equation (3a):
+
+ - c0, c1 and c2 are coefficients obtained from simulated data;
+
+ - τ is the band effective atmospheric transmittance;
+
+ - N is the number of adjacent pixels (excluding water and cloud pixels)
+ in a spatial window of size n (i.e., N = n × n);
+
+ - Ti,k and Tj,k are Top of Atmosphere brightness temperatures (K) of
+ bands i and j for the kth pixel;
+
+ - mean(Ti) and mean(Tj) are the mean or median brightness temperatures of
+ the N pixels for the two bands.
+
+
+ The regression coefficients:
+
+ ==================================================================
+
+ * NOTE, there is a typo in the paper
+
+ [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\#sthash.ba1pt9hj.dpuf
+
+ from which the equation's coefficients are (also) published.
+
+ The correct order of constants is as below, source from the
+ referenced paper below.
+
+ ==================================================================
+
+ - c2 = -9.674
+ - c1 = 0.653
+ - c0 = 9.087
+
+ where obtained by:
+
+ - 946 cloud-free TIGR atmospheric profiles,
+ - the new high accurate atmospheric radiative transfer model MODTRAN 5.2
+ - simulating the band effective atmospheric transmittance
+
+ Model analysis indicated that this method will obtain a CWV RMSE of about
+ 0.5 g/cm2. Details about the CWV retrieval can be found in:
+
+ 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.
+
+
+ Old reference:
+
+ 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.
+ """
+
+ def __init__(self, window_size, ti, tj):
+ """
+ """
+
+ # citation
+ self.citation = ('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.')
+
+ # model constants
+ self.c2 = -9.674
+ self.c1 = 0.653
+ self.c0 = 9.087
+
+ # window of N (= n by n) pixels, adjacent pixels
+ assert window_size % 2 != 0, "Window size should be an even number!"
+ assert window_size >= 7, "Window size should be larger than 5"
+ self.window_size = window_size
+
+ self.window_height = self.window_size
+ self.window_width = self.window_size
+ self.adjacent_pixels = self._derive_adjacent_pixels()
+
+ # maps for transmittance
+ self.ti = ti
+ self.tj = tj
+
+ # mapcalc modifiers to access neighborhood pixels
+ self.modifiers_ti = self._derive_modifiers(self.ti)
+ self.modifiers_tj = self._derive_modifiers(self.tj)
+ self.modifiers = zip(self.modifiers_ti, self.modifiers_tj)
+
+ # mapcalc expression for means
+ self.mean_ti_expression = self._mean_tirs_expression(self.modifiers_ti)
+ self.mean_tj_expression = self._mean_tirs_expression(self.modifiers_tj)
+
+ # mapcalc expression for medians -- ToDo
+ self.median_ti_expression = \
+ self._median_tirs_expression(self.modifiers_ti)
+ self.median_tj_expression = \
+ self._median_tirs_expression(self.modifiers_tj)
+
+ # mapcalc expression for ratio ji
+ self.ratio_ji_expression = self._ratio_ji_expression()
+
+ # mapcalc expression for column water vapor
+ self.column_water_vapor_expression = \
+ self._column_water_vapor_expression()
+
+ def __str__(self):
+ """
+ The object's self string
+ """
+ msg = 'Expression for r.mapcalc to determine column water vapor: '
+ return msg + str(self.column_water_vapor_expression)
+
+ def compute_column_water_vapor(self, tik, tjk):
+ """
+ Compute the column water vapor based on lists of input Ti and Tj
+ values.
+
+ This is a single value production function. It does not read or return
+ a map.
+ """
+ # feed with N pixels
+ ti_mean = sum(tik) / len(tik)
+ tj_mean = sum(tjk) / len(tjk)
+
+ # numerator: sum of all (Tik - Ti_mean) * (Tjk - Tj_mean)
+ numerator_ji_terms = []
+ for ti, tj in zip(tik, tjk):
+ numerator_ji_terms.append((ti - ti_mean) * (tj - tj_mean))
+
+ numerator_ji = sum(numerator_ji_terms) * 1.0
+
+ # denominator: sum of all (Tik - Tj_mean)^2
+ denominator_ji_terms = []
+ for ti in tik:
+ term = (ti - ti_mean)**2
+ denominator_ji_terms.append(term)
+ denominator_ji = sum(denominator_ji_terms) * 1.0
+
+ # ratio ji
+ ratio_ji = numerator_ji / denominator_ji
+
+ # column water vapor
+ cwv = self.c0 + self.c1 * (ratio_ji) + self.c2 * ((ratio_ji) ** 2)
+
+ # print '{c0} + {c1}*({rji}) + {c2}*({rji})^2 = '.format(c0=self.c0,
+ # c1=self.c1,
+ # rji=ratio_ji,
+ # c2=self.c2,
+ # cwv=cwv),
+ return cwv
+
+ def _derive_adjacent_pixels(self):
+ """
+ Derive a window/grid of "adjacent" pixels:
+
+ [-1, -1] [-1, 0] [-1, 1]
+ [ 0, -1] [ 0, 0] [ 0, 1]
+ [ 1, -1] [ 1, 0] [ 1, 1]
+ """
+ return [[col-1, row-1] for col in xrange(self.window_width)
+ for row in xrange(self.window_height)]
+
+ def _derive_modifiers(self, tx):
+ """
+ Return mapcalc map modifiers for adjacent pixels for the input map tx
+ """
+ return [tx + str(pixel) for pixel in self.adjacent_pixels]
+
+ def _mean_tirs_expression(self, modifiers):
+ """
+ Return mapcalc expression for window means based on the given mapcalc
+ pixel modifiers.
+ """
+ tx_mean_expression = '{sum_of_tx} / {length_of_tx}'
+ tx_sum = '(' + ' + '.join(modifiers) + ')'
+ tx_length = len(modifiers)
+
+ return tx_mean_expression.format(sum_of_tx=tx_sum,
+ length_of_tx=tx_length)
+
+ def _median_tirs_expression(self, modifiers):
+ """
+ Return mapcalc expression for window medians based on the given mapcalc
+ pixel modifiers.
+
+ r.mapcalc has a "median" function. Thus, just return the pixel
+ modifiers.
+ """
+ tx_median_expression = 'median({pixel_modifiers})'
+ # print tx_median_expression.format(pixel_modifiers=modifiers)
+
+ return tx_median_expression
+
+ def _numerator_for_ratio(self, mean_ti, mean_tj):
+ """
+ Numerator for Ratio ji. Use this function for the step-by-step approach
+ to estimate the column water vapor from within the main code (main
+ function) of the module i.landsat8.swlst
+ """
+ if not mean_ti:
+ mean_ti = 'Ti_mean'
+
+ if not mean_tj:
+ mean_tj = 'Tj_mean'
+
+ rji_numerator = '(' + '({Ti} - {Tim}) * ({Tj} - {Tjm})' + ')'
+
+ return ' + '.join([rji_numerator.format(Ti=mod_ti,
+ Tim=mean_ti,
+ Tj=mod_tj,
+ Tjm=mean_tj)
+ for mod_ti, mod_tj in self.modifiers])
+
+ def _numerator_for_ratio_big(self, **kwargs):
+ """
+ Numerator for Ratio ji. Requires two strings, one to represent the mean
+ of 'Ti's ('mean_ti') and one for the mean of 'Tj's ('mean_tj'). Use
+ this function for the big mapcalc expression.
+
+ Example:
+ _numerator_for_ratio_big(mean_ti='Some_String',
+ mean_tj='Another_String')
+ """
+ mean_ti = kwargs.get('mean_ti', 'ti_mean')
+ mean_tj = kwargs.get('mean_tj', 'tj_mean')
+
+ terms = '({Ti} - {Tim}) * ({Tj} - {Tjm})'
+ terms = ' + '.join([terms.format(Ti=mod_ti,
+ Tim=mean_ti,
+ Tj=mod_tj,
+ Tjm=mean_tj)
+ for mod_ti, mod_tj in self.modifiers])
+
+ return terms
+
+ def _denominator_for_ratio(self, mean_ti):
+ """
+ Denominator for Ratio ji. Use this function for the step-by-step
+ approach to estimate the column water vapor from within the main code
+ (main function) of the module i.landsat8.swlst
+ """
+ if not mean_ti:
+ mean_ti = 'Ti_mean'
+
+ rji_denominator = '({Ti} - {Tim})^2'
+
+ return ' + '.join([rji_denominator.format(Ti=mod,
+ Tim=mean_ti)
+ for mod in self.modifiers_ti])
+
+ def _denominator_for_ratio_big(self, **kwargs):
+ """
+ Denominator for Ratio ji. Use this function for the big mapcalc
+ expression.
+
+ Example:
+ _denominator_for_ratio_big(mean_ti='Some_String')
+ """
+ mean_ti = kwargs.get('mean_ti', 'ti_mean')
+ terms = '({Ti} - {Tim})^2'
+ terms = ' + '.join([terms.format(Ti=mod,
+ Tim=mean_ti)
+ for mod in self.modifiers_ti])
+
+ return terms
+
+ def _ratio_ji_expression(self):
+ """
+ Returns a mapcalc expression for the Ratio ji, part of the column water
+ vapor retrieval model.
+ """
+ rji_numerator = self._numerator_for_ratio(mean_ti=DUMMY_Ti_MEAN,
+ mean_tj=DUMMY_Tj_MEAN)
+
+ rji_denominator = self._denominator_for_ratio(mean_ti=DUMMY_Ti_MEAN)
+
+ rji = '( {numerator} ) / ( {denominator} )'
+ rji = rji.format(numerator=rji_numerator, denominator=rji_denominator)
+
+ return rji
+
+ def _column_water_vapor_expression(self):
+ """
+ Use this function for the step-by-step approach to estimate the column
+ water vapor from within the main code (main function) of the module
+ i.landsat8.swlst
+ """
+ cwv_expression = '({c0}) + ({c1}) * ({Rji}) + ({c2}) * ({Rji})^2'
+
+ return cwv_expression.format(c0=self.c0, c1=self.c1,
+ Rji=DUMMY_Rji, c2=self.c2)
+
+ def _big_cwv_expression(self):
+ """
+ Build and return a valid mapcalc expression for deriving a Column
+ Water Vapor map from Landsat8's brightness temperature channels
+ B10, B11 based on the MSWCVM method (see citation).
+ """
+ modifiers_ti = self._derive_modifiers(self.ti)
+ ti_sum = '(' + ' + '.join(modifiers_ti) + ')'
+ ti_length = len(modifiers_ti)
+ ti_mean = '{sum} / {length}'.format(sum=ti_sum, length=ti_length)
+
+ modifiers_tj = self._derive_modifiers(self.tj)
+ tj_sum = '(' + ' + '.join(modifiers_tj) + ')'
+ tj_length = len(modifiers_tj)
+ tj_mean = '{sum} / {length}'.format(sum=tj_sum, length=tj_length)
+
+ string_for_mean_ti = 'ti_mean'
+ string_for_mean_tj = 'tj_mean'
+
+ numerator = self._numerator_for_ratio_big(mean_ti=string_for_mean_ti,
+ mean_tj=string_for_mean_tj)
+
+ denominator = \
+ self._denominator_for_ratio_big(mean_ti=string_for_mean_ti)
+
+ cwv = ('eval('
+ '\ \n ti_mean = {tim},'
+ '\ \n'
+ '\ \n tj_mean = {tjm},'
+ '\ \n'
+ '\ \n numerator = {numerator},'
+ '\ \n'
+ '\ \n denominator = {denominator},'
+ '\ \n'
+ '\ \n rji = numerator / denominator,'
+ '\ \n'
+ '\ \n {c0} + {c1} * (rji) + {c2} * (rji)^2)')
+
+ cwv_expression = cwv.format(tim=ti_mean, tjm=tj_mean,
+ numerator=numerator,
+ denominator=denominator,
+ c0=self.c0, c1=self.c1, c2=self.c2)
+
+ return cwv_expression
+
+ def _big_cwv_expression_median():
+ """
+ Build and return a valid mapcalc expression for deriving a Column
+ Water Vapor map from Landsat8's brightness temperature channels
+ B10, B11 based on the MSWCVM method (see citation).
+ """
+ modifiers_ti = self._derive_modifiers(self.ti)
+ ti_median = 'median({modifiers}'.format(modifiers=modifiers_ti)
+
+ modifiers_tj = self._derive_modifiers(self.tj)
+ tj_median = 'median({modifiers}'.format(modifiers=modifiers_tj)
+
+ string_for_median_ti = 'ti_median'
+ string_for_median_tj = 'tj_median'
+
+ numerator = self._numerator_for_ratio_big(median_ti=string_for_median_ti,
+ median_tj=string_for_median_tj)
+
+ denominator = \
+ self._denominator_for_ratio_big(median_ti=string_for_median_ti)
+
+ cwv = ('eval('
+ '\ \n ti_median = {tim},'
+ '\ \n'
+ '\ \n tj_median = {tjm},'
+ '\ \n'
+ '\ \n numerator = {numerator},'
+ '\ \n'
+ '\ \n denominator = {denominator},'
+ '\ \n'
+ '\ \n rji = numerator / denominator,'
+ '\ \n'
+ '\ \n {c0} + {c1} * (rji) + {c2} * (rji)^2)')
+
+ cwv_expression = cwv.format(tim=ti_median, tjm=tj_median,
+ numerator=numerator,
+ denominator=denominator,
+ c0=self.c0, c1=self.c1, c2=self.c2)
+
+ return cwv_expression
+
+# reusable & stand-alone
+if __name__ == "__main__":
+ print ('Atmpspheric column water vapor retrieval '
+ 'from Landsat 8 TIRS data.'
+ ' (Running as stand-alone tool?)')
Added: grass-addons/grass7/imagery/i.landsat8.swlst/csv_to_dictionary.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/csv_to_dictionary.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/csv_to_dictionary.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,396 @@
+"""
+Convert csv data to a dictionary with namedtuples as values. However, a string
+csv-file look-a-like may be used directly. Currently, the latter option is
+used.
+
+ToDo:
+* Add usage examples!
+* Clean up which test to use for which csv file!
+* Clean up which transform function to use according to csv file!
+* Clean again...
+* Deduplication!
+"""
+
+# based on: <http://pastebin.com/tnyhmCJz>
+# see: <http://stackoverflow.com/q/29141609/1172302>
+
+# real data
+AE_STRING = '''Emissivity Class|TIRS10|TIRS11
+Cropland|0.971|0.968
+Forest|0.995|0.996
+Grasslands|0.97|0.971
+Shrublands|0.969|0.97
+Wetlands|0.992|0.998
+Waterbodies|0.992|0.998
+Tundra|0.98|0.984
+Impervious|0.973|0.981
+Barren Land|0.969|0.978
+Snow and ice|0.992|0.998'''
+
+CWV_STRING = '''Range|CWV|b0|b1|b2|b3|b4|b5|b6|b7|RMSE
+Range 1|(0.0, 2.5)|-2.78009|1.01408|0.15833|-0.34991|4.04487|3.55414|-8.88394|0.09152|0.34
+Range 2|(2.0, 3.5)|11.00824|0.95995|0.17243|-0.28852|7.11492|0.42684|-6.62025|-0.06381|0.60
+Range 3|(3.0, 4.5)|9.62610|0.96202|0.13834|-0.17262|7.87883|5.17910|-13.26611|-0.07603|0.71
+Range 4|(4.0, 5.5)|0.61258|0.99124|0.10051|-0.09664|7.85758|6.86626|-15.00742|-0.01185|0.86
+Range 5|(5.0, 6.3)|-0.34808|0.98123|0.05599|-0.03518|11.96444|9.06710|-14.74085|-0.20471|0.93
+Range 6|(0.0, 6.3)|-0.41165|1.00522|0.14543|-0.27297|4.06655|-6.92512|-18.27461|0.24468|0.87'''
+
+# required librairies
+import sys
+import csv
+from collections import namedtuple
+import random
+
+
+# helper functions
+def set_csvfile():
+ """
+ Set user defined csvfile, if any
+ """
+ if len(sys.argv) > 1:
+ return sys.argv[1]
+ else:
+ return False
+
+
+def is_number(value):
+ '''
+ Check if input is a number
+ '''
+ try:
+ float(value) # for int, long and float
+ except ValueError:
+ try:
+ complex(value) # for complex
+ except ValueError:
+ return False
+ return float(value)
+
+def to_tuple(string):
+ """
+ Convert string to tuple.
+ """
+ return tuple(map(float, string[1:-1].split(',')))
+
+def replace_dot_comma_space(string):
+ """
+ Source: <http://stackoverflow.com/a/9479972/1172302>
+ """
+ replacements = ('.', ''), (', ', '_'), (',', '_'), (' ', '_'), ('(', ''), (')', ''), ('/', '_')
+ return reduce(lambda alpha, omega: alpha.replace(*omega),
+ replacements, string)
+
+
+def csv_reader(csv_file):
+ '''
+ Transforms csv from a file into a multiline string. For example,
+ the following csv
+
+ --%<---
+ Emissivity Class|TIRS10|TIRS11
+ Cropland|0.971|0.968
+ Forest|0.995|0.996
+ Grasslands|0.97|0.971
+ Shrublands|0.969|0.97
+ Wetlands|0.992|0.998
+ Waterbodies|0.992|0.998
+ Tundra|0.98|0.984
+ Impervious|0.973|0.981
+ Barren Land|0.969|0.978
+ Snow and ice|0.992|0.998
+ --->%--
+
+ will be returned as:
+
+ """Emissivity Class|TIRS10|TIRS11
+ Cropland|0.971|0.968
+ Forest|0.995|0.996
+ Grasslands|0.97|0.971
+ Shrublands|0.969|0.97
+ Wetlands|0.992|0.998
+ Waterbodies|0.992|0.998
+ Tundra|0.98|0.984
+ Impervious|0.973|0.981
+ Barren Land|0.969|0.978
+ Snow and ice|0.992|0.998"""
+ '''
+ with open(csv_file, 'rb') as csvfile:
+ csvreader = csv.reader(csvfile, delimiter="|") # delimiter?
+ string = str()
+ for row in csvreader:
+ string += '\n' + str('|'.join(row))
+ string = string.strip('\n') # remove first newline!
+ return string
+
+
+def csv_to_dictionary(csv):
+ '''
+ Transform input from "special" csv into a python dictionary with namedtuples
+ as values. Note, "strings" of interest are hardcoded!
+
+ Also, fix the re-definition of the function transform(). See
+ <http://stackoverflow.com/q/30204197/1172302>
+ '''
+ # split input in rows
+ rows = csv.split('\n')
+ dictionary = {} # empty dictionary
+ fields = rows.pop(0).split('|')[1:] # header
+
+ strings = ('TIRS10', 'TIRS11')
+ if any(string in fields for string in strings):
+
+ def transform(row):
+ '''
+ Transform an input row in to a named tuple, then feed it in to a
+ dictionary.
+ '''
+ # split row in elements
+ elements = row.split('|')
+
+ # key: 1st column, replace
+ key = replace_dot_comma_space(elements[0])
+
+ # namedtuple
+ ect = namedtuple(key, [fields[0], fields[1]])
+
+ # feed namedtuples
+ ect.TIRS10 = is_number(elements[1])
+ ect.TIRS11 = is_number(elements[2])
+
+ # feed dictionary
+ dictionary[key] = dictionary.get(key, ect)
+
+ strings = ('b0', 'b1', 'b2', 'b3', 'b4', 'b5', 'b6', 'b7')
+ if any(string in fields for string in strings):
+
+ def transform(row):
+ '''
+ Transform an input row in to a named tuple, then feed it in to a
+ dictionary.
+ '''
+ # split row in elements
+ elements = row.split('|')
+
+ # key: 1st column, replace
+ key = replace_dot_comma_space(elements[0])
+
+ # *** small modification for the CWV field ***
+ fields[0] = 'cwv'
+
+ # named tuples
+ cwv = namedtuple(key,
+ [replace_dot_comma_space(fields[0]),
+ replace_dot_comma_space(fields[1]),
+ replace_dot_comma_space(fields[2]),
+ replace_dot_comma_space(fields[3]),
+ replace_dot_comma_space(fields[4]),
+ replace_dot_comma_space(fields[5]),
+ replace_dot_comma_space(fields[6]),
+ replace_dot_comma_space(fields[7]),
+ replace_dot_comma_space(fields[8]),
+ replace_dot_comma_space(fields[9])])
+
+ # feed named tuples
+ cwv.subrange = to_tuple(elements[1])
+ cwv.b0 = is_number(elements[2])
+ cwv.b1 = is_number(elements[3])
+ cwv.b2 = is_number(elements[4])
+ cwv.b3 = is_number(elements[5])
+ cwv.b4 = is_number(elements[6])
+ cwv.b5 = is_number(elements[7])
+ cwv.b6 = is_number(elements[8])
+ cwv.b7 = is_number(elements[9])
+ cwv.rmse = is_number(elements[10])
+ dictionary[key] = dictionary.get(key, cwv) # feed dictionary
+
+ map(transform, rows)
+ return dictionary
+
+
+def get_average_emissivities():
+ """
+ Read comma separated values for average emissivities and return a
+ dictionary wiht named tuples
+ """
+
+ try:
+ # read csv for average emissivities, convert to string
+ csvstring = csv_reader("average_emissivity.csv")
+
+ except:
+ csvstring = AE_STRING
+
+ # convert string to dictionary
+ average_emissivities = csv_to_dictionary(csvstring)
+
+ # return the dictionary with coefficients
+ return average_emissivities
+
+
+def get_column_water_vapor():
+ """
+ Read comma separated values for column water vapor coefficients and return
+ a dictionary wiht named tuples
+ """
+
+ try:
+ # read csv for average emissivities, convert to string
+ csvstring = csv_reader("cwv_coefficients.csv")
+ except:
+ csvstring = CWV_STRING
+
+ # convert string to dictionary
+ column_water_vapor_coefficients = csv_to_dictionary(csvstring)
+
+ # return the dictionary with coefficients
+ return column_water_vapor_coefficients
+
+
+# main
+def main():
+ """
+ Main function:
+ - reads a special csv file (or a multi-line string)
+ - converts and returns a dictionary which contains named tupples
+
+ - accepted csv are:
+ - average emissivity coefficients
+ - column water vapor
+ """
+
+ # user requested file?
+ global CSVFILE
+
+ if set_csvfile():
+ CSVFILE = set_csvfile()
+ print " * Reading comma separated values from:", CSVFILE
+
+ else:
+ raise IOError('Please define a file to read comma-separated-values from!')
+
+ # convert csv file to string
+ csvstring = csv_reader(CSVFILE)
+
+ # convert string to dictionary
+ coefficients_dictionary = csv_to_dictionary(csvstring) # csv < from string
+
+ # report on user requested file
+ if set_csvfile():
+ msg = ' > Dictionary with coefficients '
+ msg += str('(note, it contains named tuples):\n\n')
+ print msg, coefficients_dictionary
+
+ # return the dictionary with coefficients
+ return coefficients_dictionary
+
+
+# Test data
+def test_csvfile(infile):
+ '''
+ Test helper and main functions using as input a csv file.
+ '''
+ global CSVFILE
+ CSVFILE = infile
+ print "CSVFILE (global variable) = ", CSVFILE
+
+ print 'Test helper and main functions using as input a csv file.'
+ print
+
+ number = random.randint(1., 10.)
+ print " * Testing helper function 'is_number':", is_number(number)
+
+ if not infile:
+ csvfile = "average_emissivity.csv"
+ else:
+ csvfile = infile
+
+ print " * Testing 'csv_reader' on", csvfile, ":\n\n", csv_reader(csvfile)
+ print
+
+ csvstring = csv_reader(csvfile)
+ print " * Testing 'csv_to_dictionary':\n\n", csv_to_dictionary(csvstring)
+ print
+
+ d = csv_to_dictionary(csvstring)
+ somekey = random.choice(d.keys())
+ print "* Some random key:", somekey
+
+ fields = d[somekey]._fields
+ print "* Fields of namedtuple:", fields
+
+ random_field = random.choice(fields)
+ print "* Some random field:", random_field
+ # print "* Return values (namedtuple):", d[somekey].TIRS10, d[somekey].TIRS11
+ print "* Return values (namedtuple):", ('subrange', d[somekey].subrange,
+ 'b0', d[somekey].b0,
+ 'b1', d[somekey].b1,
+ 'b2', d[somekey].b2,
+ 'b3', d[somekey].b3,
+ 'b4', d[somekey].b4,
+ 'b5', d[somekey].b5,
+ 'b6', d[somekey].b6,
+ 'b7', d[somekey].b7,
+ 'rmse', d[somekey].rmse)
+
+#test_using_file(CSVFILE) # Ucomment to run test function!
+#CSVFILE = "cwv_coefficients.csv"
+#test_csvfile("cwv_coefficients.csv")
+#CSVFILE = ''
+
+
+def test(testdata):
+ '''
+ Test helper and main functions using as input a multi-line string.
+ '''
+ number = random.randint(1., 10.)
+ print " * Testing 'is_number':", is_number(number)
+ print
+
+ '''
+ Testing the process...
+ '''
+ d = csv_to_dictionary(testdata)
+ print "Dictionary is:\n", d
+ print
+
+ somekey = random.choice(d.keys())
+ print "Some random key:", somekey
+ print
+
+ fields = d[somekey]._fields
+ print "Fields of namedtuple:", fields
+ print
+
+ random_field = random.choice(fields)
+ print "Some random field:", random_field
+ print "Return values (namedtuple):", d[somekey].TIRS10, d[somekey].TIRS11
+
+testdata = '''LandCoverClass|TIRS10|TIRS11
+Cropland|0.971|0.968
+Forest|0.995|0.996
+Grasslands|0.970|0.971
+Shrublands|0.969|0.970
+Wetlands|0.992|0.998
+Waterbodies|0.992|0.998
+Tundra|0.980|0.984
+Impervious|0.973|0.981
+Barren_Land|0.969|0.978
+Snow_and_Ice|0.992|0.998'''
+
+#test(testdata) # Ucomment to run the test function!
+
+''' Output ------------------------------
+{'Wetlands': <class '__main__.Wetlands'>,
+ 'Snow_and_Ice': <class '__main__.Snow_and_Ice'>,
+ 'Impervious': <class '__main__.Impervious'>,
+ 'Grasslands': <class '__main__.Grasslands'>,
+ 'Shrublands': <class '__main__.Shrublands'>,
+ 'Cropland': <class '__main__.Cropland'>,
+ 'Tundra': <class '__main__.Tundra'>,
+ 'Barren_Land': <class '__main__.Barren_Land'>,
+ 'Forest': <class '__main__.Forest'>,
+ 'Waterbodies': <class '__main__.Waterbodies'>}
+------------------------------------ '''
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients.csv
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients.csv (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients.csv 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,7 @@
+"Range"|"CWV (g/cm2)"|"b0"|"b1"|"b2"|"b3"|"b4"|"b5"|"b6"|"b7"|"RMSE (K)"
+Range 1|(0.0, 2.5)|-2.78009|1.01408|0.15833|-0.34991|4.04487|3.55414|-8.88394|0.09152|0.34
+Range 2|(2.0, 3.5)|11.00824|0.95995|0.17243|-0.28852|7.11492|0.42684|-6.62025|-0.06381|0.60
+Range 3|(3.0, 4.5)|9.62610|0.96202|0.13834|-0.17262|7.87883|5.17910|-13.26611|-0.07603|0.71
+Range 4|(4.0, 5.5)|0.61258|0.99124|0.10051|-0.09664|7.85758|6.86626|-15.00742|-0.01185|0.86
+Range 5|(5.0, 6.3)|-0.34808|0.98123|0.05599|-0.03518|11.96444|9.06710|-14.74085|-0.20471|0.93
+Range 6|(0.0, 6.3)|-0.41165|1.00522|0.14543|-0.27297|4.06655|-6.92512|-18.27461|0.24468|0.87
Added: grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients_complete_range.csv
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients_complete_range.csv (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/cwv_coefficients_complete_range.csv 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,2 @@
+"Range"|"CWV (g/cm2)"|"b0"|"b1"|"b2"|"b3"|"b4"|"b5"|"b6"|"b7"|"RMSE (K)"
+Range 6|(0.0, 6.3)|-0.41165|1.00522|0.14543|-0.27297|4.06655|-6.92512|-18.27461|0.24468|0.87
Added: grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,262 @@
+<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>
+<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="LC81840332014146LGN00_B10.jpg"> <img src="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="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="lst_window_7.jpg"> <img src="lst_window_9.jpg"><img src="lst_window_11.jpg"></p>
+</div>
+<div class="figure">
+<p><img src="lst_window_7_detail_500px.jpg"> <img src="lst_window_9_detail_500px.jpg"> <img src="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#sthash.ba1pt9hj.dpuf</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>
Added: grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,1255 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+ MODULE: i.landsat8.swlst
+
+ AUTHOR(S): Nikos Alexandris <nik at nikosalexandris.net>
+ Created on Wed Mar 18 10:00:53 2015
+ First all-through execution: Tue May 12 21:50:42 EEST 2015
+
+ PURPOSE: A robust and practical Slit-Window (SW) algorithm estimating
+ land surface temperature, from the Thermal Infra-Red Sensor
+ (TIRS) aboard Landsat 8 with an accuracy of better than 1.0 K.
+
+ 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) linked to the Column Water Vapor.
+
+ The module's input parameters include:
+
+ - the brightness temperatures (Ti and Tj) of the two adjacent
+ TIRS channels,
+
+ - FROM-GLC land cover products and an emissivity look-up table,
+ which are a fraction of the FVC that can be estimated from the
+ red and near-infrared reflectance of the Operational Land
+ Imager (OLI).
+
+ The algorithm's flowchart (Figure 3 in the paper [0]) is:
+
+ +--------+ +--------------------------+
+ |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|
+ +------------------+
+
+ Sources:
+
+ [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#sthash.ba1pt9hj.dpuf>
+
+ [1] [Look below for the publised paper!] 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.
+
+ [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.
+
+ COPYRIGHT: (C) 2015 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.
+"""
+
+#%Module
+#% description: Practical split-window algorithm estimating Land Surface Temperature from Landsat 8 OLI/TIRS imagery (Du, Chen; Ren, Huazhong; Qin, Qiming; Meng, Jinjie; Zhao, Shaohua. 2015)
+#% keywords: imagery
+#% keywords: split window
+#% keywords: column water vapor
+#% keywords: land surface temperature
+#% keywords: lst
+#% keywords: landsat8
+#%End
+
+#%flag
+#% key: i
+#% description: Print out model equations, citation
+#%end
+
+#%flag
+#% key: k
+#% description: Keep current computational region settings
+#%end
+
+#%flag
+#% key: t
+#% description: Time-stamping the output LST (and optional CWV) map
+#%end
+
+#%flag
+#% key: c
+#% description: Convert LST output to celsius degrees, apply color table
+#%end
+
+#%flag
+#% key: n
+#% description: Set zero digital numbers in b10, b11 to NULL | ToDo: Perform in copy of input input maps!
+#%end
+
+#%option G_OPT_F_INPUT
+#% key: mtl
+#% key_desc: filename
+#% description: Landsat8 metadata file (MTL)
+#% required: no
+#%end
+
+#%option G_OPT_R_BASENAME_INPUT
+#% key: prefix
+#% key_desc: basename
+#% type: string
+#% label: OLI/TIRS band names prefix
+#% description: Prefix of Landsat8 OLI/TIRS band names
+#% required: no
+#%end
+
+##%rules
+##% collective: prefix, mtl
+##%end
+
+#%option G_OPT_R_INPUT
+#% key: b10
+#% key_desc: name
+#% description: TIRS 10 (10.60 - 11.19 microns)
+#% required : no
+#%end
+
+#%rules
+#% requires_all: b10, mtl
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: b11
+#% key_desc: name
+#% description: TIRS 11 (11.50 - 12.51 microns)
+#% required : no
+#%end
+
+#%rules
+#% requires_all: b11, mtl
+#%end
+
+#%option G_OPT_R_BASENAME_INPUT
+#% key: prefix_bt
+#% key_desc: basename
+#% type: string
+#% label: Prefix for output at-satellite brightness temperature maps (K)
+#% description: Prefix for brightness temperature maps (K)
+#% required: no
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: t10
+#% key_desc: name
+#% description: Brightness temperature (K) from band 10 | Overrides 'b10'
+#% required : no
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: t11
+#% key_desc: name
+#% description: Brightness temperature (K) from band 11 | Overrides 'b11'
+#% required : no
+#%end
+
+#%rules
+#% requires: b10, b11, t11
+#%end
+
+#%rules
+#% requires: b11, b10, t10
+#%end
+
+#%rules
+#% requires: t10, t11, b11
+#%end
+
+#%rules
+#% requires: t11, t10, b10
+#%end
+
+#%rules
+#% exclusive: b10, t10
+#%end
+
+#%rules
+#% exclusive: b11, t11
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: qab
+#% key_desc: name
+#% description: Landsat 8 Quality Assessment band
+#% required : no
+#%end
+
+#%option
+#% key: qapixel
+#% key_desc: pixelvalue
+#% description: Quality assessment pixel value for which to build a mask | Source: <http://landsat.usgs.gov/L8QualityAssessmentBand.php>.
+#% answer: 61440
+#% required: no
+#% multiple: yes
+#%end
+
+#%rules
+#% excludes: prefix, b10, b11, qab
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: clouds
+#% key_desc: name
+#% description: A raster map applied as an inverted MASK | Overrides 'qab'
+#% required : no
+#%end
+
+#%rules
+#% exclusive: qab, clouds
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: emissivity
+#% key_desc: name
+#% description: Land surface emissivity map | Expert use, overrides retrieving average emissivity from landcover
+#% required : no
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% key: emissivity_out
+#% key_desc: name
+#% description: Name for output emissivity map | For re-use as "emissivity=" input in subsequent trials with different spatial window sizes
+#% required: no
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: delta_emissivity
+#% key_desc: name
+#% description: Emissivity difference map for Landsat8 TIRS channels 10 and 11 | Expert use, overrides retrieving delta emissivity from landcover
+#% required : no
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% key: delta_emissivity_out
+#% key_desc: name
+#% description: Name for output delta emissivity map | For re-use as "delta_emissivity=" in subsequent trials with different spatial window sizes
+#% required: no
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: landcover
+#% key_desc: name
+#% description: FROM-GLC products covering the Landsat8 scene under processing. Source <http://data.ess.tsinghua.edu.cn/>.
+#% required : no
+#%end
+
+#%option
+#% key: emissivity_class
+#% key_desc: string
+#% description: Retrieve average emissivities only for a single land cover class (case sensitive) | Expert use
+#% options: Cropland, Forest, Grasslands, Shrublands, Wetlands, Waterbodies, Tundra, Impervious, Barren, Snow, Random
+#% required : no
+#%end
+
+#%rules
+#% required: landcover, emissivity_class
+#% exclusive: landcover, emissivity_class
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% key: lst
+#% key_desc: name
+#% description: Name for output Land Surface Temperature map
+#% required: yes
+#% answer: lst
+#%end
+
+#%option
+#% key: window
+#% key_desc: integer
+#% description: Odd number n sizing an n^2 spatial window for column water vapor retrieval | Increase to reduce spatial discontinuation in the final LST
+#% answer: 7
+#% required: yes
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% key: cwv
+#% key_desc: name
+#% description: Name for output Column Water Vapor map | Optional
+#% required: no
+#%end
+
+# required librairies
+import os
+import sys
+sys.path.insert(1, os.path.join(os.path.dirname(sys.path[0]),
+ 'etc', 'i.landsat8.swlst'))
+
+import atexit
+import grass.script as grass
+# from grass.exceptions import CalledModuleError
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.modules.shortcuts import raster as r
+# from grass.pygrass.raster.abstract import Info
+
+from split_window_lst import *
+from landsat8_mtl import Landsat8_MTL
+
+if "GISBASE" not in os.environ:
+ print "You must be in GRASS GIS to run this program."
+ sys.exit(1)
+
+# globals
+DUMMY_MAPCALC_STRING_RADIANCE = 'Radiance'
+DUMMY_MAPCALC_STRING_DN = 'DigitalNumber'
+DUMMY_MAPCALC_STRING_T10 = 'Input_T10'
+DUMMY_MAPCALC_STRING_T11 = 'Input_T11'
+DUMMY_MAPCALC_STRING_AVG_LSE = 'Input_AVG_LSE'
+DUMMY_MAPCALC_STRING_DELTA_LSE = 'Input_DELTA_LSE'
+DUMMY_MAPCALC_STRING_FROM_GLC = 'Input_FROMGLC'
+DUMMY_MAPCALC_STRING_CWV = 'Input_CWV'
+DUMMY_Ti_MEAN = 'Mean_Ti'
+DUMMY_Tj_MEAN = 'Mean_Tj'
+DUMMY_Rji = 'Ratio_ji'
+
+
+# helper functions
+def cleanup():
+ """
+ Clean up temporary maps
+ """
+ grass.run_command('g.remove', flags='f', type="rast",
+ pattern='tmp.{pid}*'.format(pid=os.getpid()), quiet=True)
+
+
+def tmp_map_name(name):
+ """
+ Return a temporary map name, for example:
+
+ tmp_avg_lse = tmp + '.avg_lse'
+ """
+ temporary_file = grass.tempfile()
+ tmp = "tmp." + grass.basename(temporary_file) # use its basename
+ return tmp + '.' + str(name)
+
+
+def run(cmd, **kwargs):
+ """
+ Pass required arguments to grass commands (?)
+ """
+ grass.run_command(cmd, quiet=True, **kwargs)
+
+
+def save_map(mapname):
+ """
+ Helper function to save some in-between maps, assisting in debugging
+ """
+ # run('r.info', map=mapname, flags='r')
+ run('g.copy', raster=(mapname, 'DebuggingMap'))
+
+
+def random_digital_numbers(count=2):
+ """
+ Return a user-requested amount of random Digital Number values for testing
+ purposes ranging in 12-bit
+ """
+ digital_numbers = []
+
+ for dn in range(0, count):
+ digital_numbers.append(random.randint(1, 2**12))
+
+ if count == 1:
+ return digital_numbers[0]
+
+ return digital_numbers
+
+
+def random_column_water_vapor_subrange():
+ """
+ Helper function, while coding and testing, returning a random column water
+ vapor key to assist in testing the module.
+ """
+ cwvkey = random.choice(COLUMN_WATER_VAPOUR.keys())
+ # COLUMN_WATER_VAPOUR[cwvkey].subrange
+ # COLUMN_WATER_VAPOUR[cwvkey].rmse
+ return cwvkey
+
+
+def random_column_water_vapor_value():
+ """
+ Helper function, while coding and testing, returning a random value for
+ column water vapor.
+ """
+ return random.uniform(0.0, 6.3)
+
+
+def extract_number_from_string(string):
+ """
+ Extract the (integer) number from a string. Meand to be used with band
+ names. For example:
+
+ print extract_number_from_string('B10')
+
+ will return
+
+ 10
+ """
+ import re
+ return str(re.findall(r"[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?",
+ string)[-1])
+
+
+def add_timestamp(mtl_filename, outname):
+ """
+ Retrieve metadata from MTL file.
+ """
+ import datetime
+ metadata = Landsat8_MTL(mtl_filename)
+
+ # required format is: day=integer month=string year=integer time=hh:mm:ss.dd
+ acquisition_date = str(metadata.date_acquired) ### FixMe ###
+ acquisition_date = datetime.datetime.strptime(acquisition_date, '%Y-%m-%d').strftime('%d %b %Y')
+ acquisition_time = str(metadata.scene_center_time)[0:8]
+ date_time_string = acquisition_date + ' ' + acquisition_time
+
+ #msg = "Date and time of acquisition: " + date_time_string
+ #grass.verbose(msg)
+
+ run('r.timestamp', map=outname, date=date_time_string)
+
+ del(date_time_string)
+
+
+def digital_numbers_to_radiance(outname, band, radiance_expression):
+ """
+ Convert Digital Number values to TOA Radiance. For details, see in Landsat8
+ class. Zero (0) DNs set to NULL here (not via the class' function).
+ """
+ if null:
+ msg = "\n|i Setting zero (0) Digital Numbers in {band} to NULL"
+ msg = msg.format(band=band)
+ g.message(msg)
+ run('r.null', map=band, setnull=0)
+
+ msg = "\n|i Rescaling {band} digital numbers to spectral radiance "
+ msg = msg.format(band=band)
+
+ if info:
+ msg += '| Expression: '
+ msg += radiance_expression
+ g.message(msg)
+ radiance_expression = replace_dummies(radiance_expression,
+ instring=DUMMY_MAPCALC_STRING_DN,
+ outstring=band)
+ radiance_equation = equation.format(result=outname,
+ expression=radiance_expression)
+
+ grass.mapcalc(radiance_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+ #run('r.univar', map=outname)
+
+ del(radiance_expression)
+ del(radiance_equation)
+
+
+def radiance_to_brightness_temperature(outname, radiance, temperature_expression):
+ """
+ Convert Spectral Radiance to At-Satellite Brightness Temperature. For
+ details see Landsat8 class.
+ """
+ temperature_expression = replace_dummies(temperature_expression,
+ instring=DUMMY_MAPCALC_STRING_RADIANCE,
+ outstring=radiance)
+
+ msg = "\n|i Converting spectral radiance to at-Satellite Temperature "
+ if info:
+ msg += "| Expression: " + str(temperature_expression)
+ g.message(msg)
+
+ temperature_equation = equation.format(result=outname,
+ expression=temperature_expression)
+
+ grass.mapcalc(temperature_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+ #run('r.univar', map=outname)
+
+ del(temperature_expression)
+ del(temperature_equation)
+
+
+def tirs_to_at_satellite_temperature(tirs_1x, mtl_file):
+ """
+ Helper function to convert TIRS bands 10 or 11 in to at-satellite
+ temperatures.
+
+ This function uses the pre-defined functions:
+
+ - extract_number_from_string()
+ - digital_numbers_to_radiance()
+ - radiance_to_brightness_temperature()
+
+ The inputs are:
+
+ - a name for the input tirs band (10 or 11)
+ - a Landsat8 MTL file
+
+ The output is a temporary at-Satellite Temperature map.
+ """
+ # which band number and MTL file
+ band_number = extract_number_from_string(tirs_1x)
+ tmp_radiance = tmp_map_name('radiance') + '.' + band_number
+ tmp_brightness_temperature = tmp_map_name('brightness_temperature') + '.' + \
+ band_number
+ landsat8 = Landsat8_MTL(mtl_file)
+
+ # rescale DNs to spectral radiance
+ radiance_expression = landsat8.toar_radiance(band_number)
+ digital_numbers_to_radiance(tmp_radiance, tirs_1x, radiance_expression)
+
+ # convert spectral radiance to at-satellite temperature
+ temperature_expression = landsat8.radiance_to_temperature(band_number)
+ radiance_to_brightness_temperature(tmp_brightness_temperature,
+ tmp_radiance,
+ temperature_expression)
+
+ del(radiance_expression)
+ del(temperature_expression)
+
+ # save Brightness Temperature map?
+ if brightness_temperature_prefix:
+ bt_output = brightness_temperature_prefix + band_number
+ run('g.copy', raster=(tmp_brightness_temperature, bt_output))
+ del(bt_output)
+
+ return tmp_brightness_temperature
+
+
+def mask_clouds(qa_band, qa_pixel):
+ """
+ ToDo:
+
+ - a better, independent mechanism for QA. --> see also Landsat8 class.
+ - support for multiple qa_pixel values (eg. input as a list of values)
+
+ Create and apply a cloud mask based on the Quality Assessment Band
+ (BQA.) Source: <http://landsat.usgs.gov/L8QualityAssessmentBand.php
+
+ See also:
+ http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/#Applying_the_Landsat_8_Quality_Assessment_%28QA%29_Band
+ """
+ msg = ('\n|i Masking for pixel values <{qap}> '
+ 'in the Quality Assessment band.'.format(qap=qa_pixel))
+ g.message(msg)
+
+ tmp_cloudmask = tmp_map_name('cloudmask')
+ qabits_expression = 'if({band} == {pixel}, 1, null())'.format(band=qa_band,
+ pixel=qa_pixel)
+
+ cloud_masking_equation = equation.format(result=tmp_cloudmask,
+ expression=qabits_expression)
+ grass.mapcalc(cloud_masking_equation)
+
+ r.mask(raster=tmp_cloudmask, flags='i', overwrite=True)
+
+ # save for debuging
+ #save_map(tmp_cloudmask)
+
+ del(qabits_expression)
+ del(cloud_masking_equation)
+
+
+def replace_dummies(string, *args, **kwargs):
+ """
+ Replace DUMMY_MAPCALC_STRINGS (see SplitWindowLST class for it)
+ with input maps ti, tj (here: t10, t11).
+
+ - in_ti and in_tj are the "input" strings, for example:
+ in_ti = 'Input_T10' and in_tj = 'Input_T11'
+
+ - out_ti and out_tj are the output strings which correspond to map
+ names, user-fed or in-between temporary maps, for example:
+ out_ti = t10 and out_tj = t11
+
+ or
+
+ out_ti = tmp_ti_mean and out_tj = tmp_ti_mean
+
+ (Idea sourced from: <http://stackoverflow.com/a/9479972/1172302>)
+ """
+ inout = set(['instring', 'outstring'])
+ # if inout.issubset(set(kwargs)): # alternative
+ if inout == set(kwargs):
+ instring = kwargs.get('instring', 'None')
+ outstring = kwargs.get('outstring', 'None')
+
+ # end comma important!
+ replacements = (str(instring), str(outstring)),
+
+ in_tij_out = set(['in_ti', 'out_ti', 'in_tj', 'out_tj'])
+
+ if in_tij_out == set(kwargs):
+ in_ti = kwargs.get('in_ti', 'None')
+ out_ti = kwargs.get('out_ti', 'None')
+ in_tj = kwargs.get('in_tj', 'None')
+ out_tj = kwargs.get('out_tj', 'None')
+
+ replacements = (in_ti, str(out_ti)), (in_tj, str(out_tj))
+
+ in_tijm_out = set(['in_ti', 'out_ti', 'in_tj', 'out_tj',
+ 'in_tim', 'out_tim', 'in_tjm', 'out_tjm'])
+
+ if in_tijm_out == set(kwargs):
+ in_ti = kwargs.get('in_ti', 'None')
+ out_ti = kwargs.get('out_ti', 'None')
+ in_tj = kwargs.get('in_tj', 'None')
+ out_tj = kwargs.get('out_tj', 'None')
+ in_tim = kwargs.get('in_tim', 'None')
+ out_tim = kwargs.get('out_tim', 'None')
+ in_tjm = kwargs.get('in_tjm', 'None')
+ out_tjm = kwargs.get('out_tjm', 'None')
+
+ replacements = (in_ti, str(out_ti)), (in_tj, str(out_tj)), \
+ (in_tim, str(out_tim)), (in_tjm, str(out_tjm))
+
+ in_cwv_out = set(['in_ti', 'out_ti', 'in_tj', 'out_tj', 'in_cwv',
+ 'out_cwv'])
+
+ if in_cwv_out == set(kwargs):
+ in_cwv = kwargs.get('in_cwv', 'None')
+ out_cwv = kwargs.get('out_cwv', 'None')
+ in_ti = kwargs.get('in_ti', 'None')
+ out_ti = kwargs.get('out_ti', 'None')
+ in_tj = kwargs.get('in_tj', 'None')
+ out_tj = kwargs.get('out_tj', 'None')
+
+ replacements = (in_ti, str(out_ti)), (in_tj, str(out_tj)), \
+ (in_cwv, str(out_cwv))
+
+ in_lst_out = set(['in_ti', 'out_ti', 'in_tj', 'out_tj', 'in_cwv',
+ 'out_cwv', 'in_avg_lse', 'out_avg_lse', 'in_delta_lse',
+ 'out_delta_lse'])
+
+ if in_lst_out == set(kwargs):
+ in_cwv = kwargs.get('in_cwv', 'None')
+ out_cwv = kwargs.get('out_cwv', 'None')
+ in_ti = kwargs.get('in_ti', 'None')
+ out_ti = kwargs.get('out_ti', 'None')
+ in_tj = kwargs.get('in_tj', 'None')
+ out_tj = kwargs.get('out_tj', 'None')
+ in_avg_lse = kwargs.get('in_avg_lse', 'None')
+ out_avg_lse = kwargs.get('out_avg_lse', 'None')
+ in_delta_lse = kwargs.get('in_delta_lse', 'None')
+ out_delta_lse = kwargs.get('out_delta_lse', 'None')
+
+ replacements = (in_ti, str(out_ti)), \
+ (in_tj, str(out_tj)), \
+ (in_cwv, str(out_cwv)), \
+ (in_avg_lse, str(out_avg_lse)), \
+ (in_delta_lse, str(out_delta_lse))
+
+ return reduce(lambda alpha, omega: alpha.replace(*omega),
+ replacements, string)
+
+
+def determine_average_emissivity(outname, landcover_map, avg_lse_expression):
+ """
+ Produce an average emissivity map based on FROM-GLC map covering the region
+ of interest.
+ """
+ msg = ('\n|i Determining average land surface emissivity based on a '
+ 'look-up table ')
+ if info:
+ msg += ('| Expression:\n\n {exp}')
+ msg = msg.format(exp=avg_lse_expression)
+ g.message(msg)
+
+ avg_lse_expression = replace_dummies(avg_lse_expression,
+ instring=DUMMY_MAPCALC_STRING_FROM_GLC,
+ outstring=landcover_map)
+
+ avg_lse_equation = equation.format(result=outname,
+ expression=avg_lse_expression)
+
+ grass.mapcalc(avg_lse_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ del(avg_lse_expression)
+ del(avg_lse_equation)
+
+ # save land surface emissivity map?
+ if emissivity_output:
+ run('g.copy', raster=(outname, emissivity_output))
+
+
+def determine_delta_emissivity(outname, landcover_map, delta_lse_expression):
+ """
+ Produce a delta emissivity map based on the FROM-GLC map covering the
+ region of interest.
+ """
+ msg = ('\n|i Determining delta land surface emissivity based on a '
+ 'look-up table ')
+ if info:
+ msg += ('| Expression:\n\n {exp}')
+ msg = msg.format(exp=delta_lse_expression)
+ g.message(msg)
+
+ delta_lse_expression = replace_dummies(delta_lse_expression,
+ instring=DUMMY_MAPCALC_STRING_FROM_GLC,
+ outstring=landcover_map)
+
+ delta_lse_equation = equation.format(result=outname,
+ expression=delta_lse_expression)
+
+ grass.mapcalc(delta_lse_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ del(delta_lse_expression)
+ del(delta_lse_equation)
+
+ # save delta land surface emissivity map?
+ if delta_emissivity_output:
+ run('g.copy', raster=(outname, delta_emissivity_output))
+
+
+def get_cwv_window_means(outname, t1x, t1x_mean_expression):
+ """
+
+ ***
+ This function is NOT used. It was part of an initial step-by-step approach,
+ while coding and testing. Kept for future plans!?
+ ***
+
+ Get window means for T1x
+ """
+ msg = ('\n |i Deriving window means from {Tx} ')
+ msg += ('using the expression:\n {exp}')
+ msg = msg.format(Tx=t1x, exp=t1x_mean_expression)
+ g.message(msg)
+
+ tx_mean_equation = equation.format(result=outname,
+ expression=t1x_mean_expression)
+ grass.mapcalc(tx_mean_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ del(t1x_mean_expression)
+ del(tx_mean_equation)
+
+ # save for debuging
+ #save_map(outname)
+
+
+def estimate_ratio_ji(outname, tmp_ti_mean, tmp_tj_mean, ratio_expression):
+ """
+
+ ***
+ This function is NOT used. It was part of an initial step-by-step approach,
+ while coding and testing. Kept for future plans!?
+ ***
+
+ Estimate Ratio ji for the Column Water Vapor retrieval equation.
+ """
+ msg = '\n |i Estimating ratio Rji...'
+ msg += '\n' + ratio_expression
+ g.message(msg)
+
+ ratio_expression = replace_dummies(ratio_expression,
+ in_ti=DUMMY_Ti_MEAN, out_ti=tmp_ti_mean,
+ in_tj=DUMMY_Tj_MEAN, out_tj=tmp_tj_mean)
+
+ ratio_equation = equation.format(result=outname,
+ expression=ratio_expression)
+
+ grass.mapcalc(ratio_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ # save for debuging
+ #save_map(outname)
+
+
+def estimate_column_water_vapor(outname, ratio, cwv_expression):
+ """
+
+ ***
+ This function is NOT used. It was part of an initial step-by-step approach,
+ while coding and testing. Kept for future plans!?
+ ***
+
+ """
+ msg = "\n|i Estimating atmospheric column water vapor "
+ msg += '| Mapcalc expression: '
+ msg += cwv_expression
+ g.message(msg)
+
+ cwv_expression = replace_dummies(cwv_expression,
+ instring=DUMMY_Rji,
+ outstring=ratio)
+
+ cwv_equation = equation.format(result=outname, expression=cwv_expression)
+
+ grass.mapcalc(cwv_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ # save Column Water Vapor map?
+ if cwv_output:
+ run('g.copy', raster=(outname, cwv_output))
+
+ # save for debuging
+ #save_map(outname)
+
+
+def estimate_cwv_big_expression(outname, t10, t11, cwv_expression):
+ """
+ Derive a column water vapor map using a single mapcalc expression based on
+ eval.
+
+ *** To Do: evaluate -- does it work correctly? *** !
+ """
+ msg = "\n|i Estimating atmospheric column water vapor "
+ if info:
+ msg += '| Expression:\n'
+ g.message(msg)
+
+ if info:
+ msg = replace_dummies(cwv_expression,
+ in_ti=t10, out_ti='T10',
+ in_tj=t11, out_tj='T11')
+ msg += '\n'
+ print msg
+
+ cwv_equation = equation.format(result=outname, expression=cwv_expression)
+ grass.mapcalc(cwv_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ # save Column Water Vapor map?
+ if cwv_output:
+
+ # strings for metadata
+ history_cwv = 'FixMe -- Column Water Vapor model: '
+ history_cwv += 'FixMe -- Add equation?'
+ title_cwv = 'Column Water Vapor'
+ description_cwv = 'Column Water Vapor'
+ units_cwv = 'g/cm^2'
+ source1_cwv = 'FixMe'
+ source2_cwv = 'FixMe'
+
+ # history entry
+ run("r.support", map=outname, title=title_cwv,
+ units=units_cwv, description=description_cwv,
+ source1=source1_cwv, source2=source2_cwv,
+ history=history_cwv)
+
+ run('g.copy', raster=(outname, cwv_output))
+
+ del(cwv_expression)
+ del(cwv_equation)
+
+
+def estimate_lst(outname, t10, t11, avg_lse_map, delta_lse_map, cwv_map, lst_expression):
+ """
+ Produce a Land Surface Temperature map based on a mapcalc expression
+ returned from a SplitWindowLST object.
+
+ Inputs are:
+
+ - brightness temperature maps t10, t11
+ - column water vapor map
+ - a temporary filename
+ - a valid mapcalc expression
+ """
+ msg = '\n|i Estimating land surface temperature '
+ if info:
+ msg += "| Expression:\n"
+ g.message(msg)
+
+ if info:
+ msg = lst_expression
+ msg += '\n'
+ print msg
+
+ # replace the "dummy" string...
+ if landcover_map:
+ split_window_expression = replace_dummies(lst_expression,
+ in_avg_lse=DUMMY_MAPCALC_STRING_AVG_LSE,
+ out_avg_lse=avg_lse_map,
+ in_delta_lse=DUMMY_MAPCALC_STRING_DELTA_LSE,
+ out_delta_lse=delta_lse_map,
+ in_cwv=DUMMY_MAPCALC_STRING_CWV,
+ out_cwv=cwv_map,
+ in_ti=DUMMY_MAPCALC_STRING_T10,
+ out_ti=t10,
+ in_tj=DUMMY_MAPCALC_STRING_T11,
+ out_tj=t11)
+ elif emissivity_class:
+ split_window_expression = replace_dummies(lst_expression,
+ in_cwv=DUMMY_MAPCALC_STRING_CWV,
+ out_cwv=cwv_map,
+ in_ti=DUMMY_MAPCALC_STRING_T10,
+ out_ti=t10,
+ in_tj=DUMMY_MAPCALC_STRING_T11,
+ out_tj=t11)
+
+ split_window_equation = equation.format(result=outname,
+ expression=split_window_expression)
+
+ grass.mapcalc(split_window_equation, overwrite=True)
+
+ if info:
+ run('r.info', map=outname, flags='r')
+
+ del(split_window_expression)
+ del(split_window_equation)
+
+
+def kelvin_to_celsius(outname, lst_map):
+ """
+ Convert Kelvin to Celsius
+ """
+ msg = '\n|i Converting LST output from Kelvin to Celsius degrees'
+ #g.message(msg)
+
+ kelvin_to_celsius_expression = '{lst} - 273.15'.format(lst=lst_map)
+ kelvin_to_celsius_equation = equation.format(result=lst_map,
+ expression=kelvin_to_celsius_expression)
+ grass.mapcalc(kelvin_to_celsius_equation, overwrite=True)
+
+ msg += '\n|i Applying the "Celsius" color table'
+ g.message(msg)
+ run('r.colors', map=lst_map, color='celsius')
+
+ del(kelvin_to_celsius_expression)
+ del(kelvin_to_celsius_equation)
+
+def main():
+ """
+ Main program
+ """
+
+ # Temporary filenames
+
+ # The following three are meant for a test step-by-step cwv estimation, see
+ # unused functions!
+
+ # tmp_ti_mean = tmp_map_name('ti_mean') # for cwv
+ # tmp_tj_mean = tmp_map_name('tj_mean') # for cwv
+ # tmp_ratio = tmp_map_name('ratio') # for cwv
+
+ tmp_avg_lse = tmp_map_name('avg_lse')
+ tmp_delta_lse = tmp_map_name('delta_lse')
+ tmp_cwv = tmp_map_name('cwv')
+ tmp_lst = tmp_map_name('lst')
+
+ # basic equation for mapcalc
+ global equation, citation_lst
+ equation = "{result} = {expression}"
+
+ # user input
+ mtl_file = options['mtl']
+
+ if not options['prefix']:
+ b10 = options['b10']
+ b11 = options['b11']
+ t10 = options['t10']
+ t11 = options['t11']
+
+ if not options['clouds']:
+ qab = options['qab']
+ cloud_map = False
+
+ else:
+ qab = False
+ cloud_map = options['clouds']
+
+ elif options['prefix']:
+ prefix = options['prefix']
+ b10 = prefix + '10'
+ b11 = prefix + '11'
+
+ if not options['clouds']:
+ qab = prefix + 'QA'
+ cloud_map = False
+
+ else:
+ cloud_map = options['clouds']
+ qab = False
+
+ qapixel = options['qapixel']
+ lst_output = options['lst']
+
+ # save Brightness Temperature maps?
+ global brightness_temperature_prefix
+ brightness_temperature_prefix = options['prefix_bt']
+
+ global cwv_output
+ cwv_window_size = int(options['window'])
+ assertion_for_cwv_window_size_msg = ('A spatial window of size 5^2 or less is not '
+ 'recommended. Please select a larger window. '
+ 'Refer to the manual\'s notes for details.')
+ assert cwv_window_size >= 7, assertion_for_cwv_window_size_msg
+ cwv_output = options['cwv']
+
+ # optional maps
+ average_emissivity_map = options['emissivity']
+ delta_emissivity_map = options['delta_emissivity']
+
+ # output for in-between maps?
+ global emissivity_output, delta_emissivity_output
+ emissivity_output = options['emissivity_out']
+ delta_emissivity_output = options['delta_emissivity_out']
+
+ global landcover_map, emissivity_class
+ landcover_map = options['landcover']
+ emissivity_class = options['emissivity_class']
+
+ # flags
+ global info, null
+ info = flags['i']
+ keep_region = flags['k']
+ celsius = flags['c']
+ timestamping = flags['t']
+ null = flags['n']
+
+ # ToDo:
+ # shell = flags['g']
+
+ #
+ # Pre-production actions
+ #
+
+ # Set Region
+ if not keep_region:
+ grass.use_temp_region() # safely modify the region
+ msg = "\n|! Matching region extent to map {name}"
+
+ # ToDo: check if extent-B10 == extent-B11? Unnecessary?
+ # Improve below!
+
+ if b10:
+ run('g.region', rast=b10)
+ msg = msg.format(name=b10)
+
+ elif t10:
+ run('g.region', rast=t10)
+ msg = msg.format(name=t10)
+
+ g.message(msg)
+
+ elif keep_region:
+ grass.warning(_('Operating on current region'))
+
+ #
+ # 1. Mask clouds
+ #
+
+ if cloud_map:
+ # user-fed cloud map?
+ msg = '\n|i Using {cmap} as a MASK'.format(cmap=cloud_map)
+ g.message(msg)
+ r.mask(raster=cloud_map, flags='i', overwrite=True)
+
+ else:
+ # using the quality assessment band and a "QA" pixel value
+ mask_clouds(qab, qapixel)
+
+ #
+ # 2. TIRS > Brightness Temperatures
+ #
+
+ if mtl_file:
+
+ # if MTL and b10 given, use it to compute at-satellite temperature t10
+ if b10:
+ # convert DNs to at-satellite temperatures
+ t10 = tirs_to_at_satellite_temperature(b10, mtl_file)
+
+ # likewise for b11 -> t11
+ if b11:
+ # convert DNs to at-satellite temperatures
+ t11 = tirs_to_at_satellite_temperature(b11, mtl_file)
+
+ #
+ # Initialise a SplitWindowLST object
+ #
+
+ split_window_lst = SplitWindowLST(emissivity_class)
+ citation_lst = split_window_lst.citation
+
+ #
+ # 3. Land Surface Emissivities
+ #
+
+ # use given fixed class?
+ if emissivity_class:
+
+ if split_window_lst.landcover_class is False:
+ # replace with meaningful error
+ g.warning('Unknown land cover class string! Note, this string '
+ 'input option is case sensitive.')
+
+ if emissivity_class == 'Random':
+ msg = "\n|! Random emissivity class selected > " + \
+ split_window_lst.landcover_class + ' '
+
+ else:
+ msg = '\n|! Retrieving average emissivities *only* for {eclass} '
+
+ if info:
+ msg += '| Average emissivities (channels 10, 11): '
+ msg += str(split_window_lst.emissivity_t10) + ', ' + \
+ str(split_window_lst.emissivity_t11)
+
+ msg = msg.format(eclass=split_window_lst.landcover_class)
+ g.message(msg)
+
+ # use the FROM-GLC map
+ elif landcover_map:
+
+ if average_emissivity_map:
+ tmp_avg_lse = average_emissivity_map
+
+ if not average_emissivity_map:
+ determine_average_emissivity(tmp_avg_lse, landcover_map,
+ split_window_lst.average_lse_mapcalc)
+ if delta_emissivity_map:
+ tmp_delta_lse = delta_emissivity_map
+
+ if not delta_emissivity_map:
+ determine_delta_emissivity(tmp_delta_lse, landcover_map,
+ split_window_lst.delta_lse_mapcalc)
+
+ #
+ # 4. Modified Split-Window Variance-Covariance Matrix > Column Water Vapor
+ #
+
+
+ if info:
+ msg = '\n|i Spatial window of size {n} for Column Water Vapor estimation: '
+ msg = msg.format(n=cwv_window_size)
+ g.message(msg)
+
+ cwv = Column_Water_Vapor(cwv_window_size, t10, t11)
+ citation_cwv = cwv.citation
+ estimate_cwv_big_expression(tmp_cwv, t10, t11, cwv._big_cwv_expression())
+
+ #
+ # 5. Estimate Land Surface Temperature
+ #
+
+ if info and emissivity_class == 'Random':
+ msg = '\n|* Will pick a random emissivity class!'
+ grass.verbose(msg)
+
+ estimate_lst(tmp_lst, t10, t11,
+ tmp_avg_lse, tmp_delta_lse, tmp_cwv,
+ split_window_lst.sw_lst_mapcalc)
+
+ #
+ # Post-production actions
+ #
+
+ # remove MASK
+ r.mask(flags='r', verbose=True)
+
+ # time-stamping
+ if timestamping:
+ add_timestamp(mtl_file, tmp_lst)
+
+ if cwv_output:
+ add_timestamp(mtl_file, cwv_output)
+
+ # convert to celsius and apply color table?
+ if celsius:
+ kelvin_to_celsius(tmp_lst, tmp_lst)
+
+ else:
+ # color table for kelvin
+ run('r.colors', map=tmp_lst, color='kelvin')
+
+ # ToDo: helper function for r.support
+ # strings for metadata
+ history_lst = '\n' + citation_lst
+ history_lst += '\n\n' + citation_cwv
+ history_lst += '\n\nSplit-Window model: '
+ history_lst += split_window_lst._equation # :wsw_lst_mapcalc
+ description_lst = ('Land Surface Temperature derived from a split-window algorithm. ')
+
+ if celsius:
+ title_lst = 'Land Surface Temperature (C)'
+ units_lst = 'Celsius'
+
+ else:
+ title_lst = 'Land Surface Temperature (K)'
+ units_lst = 'Kelvin'
+
+ landsat8_metadata = Landsat8_MTL(mtl_file)
+ source1_lst = landsat8_metadata.scene_id
+ source2_lst = landsat8_metadata.origin
+
+ # history entry
+ run("r.support", map=tmp_lst, title=title_lst,
+ units=units_lst, description=description_lst,
+ source1=source1_lst, source2=source2_lst,
+ history=history_lst)
+
+ # (re)name the LST product
+ run("g.rename", rast=(tmp_lst, lst_output))
+
+ # restore region
+ if not keep_region:
+ grass.del_temp_region() # restoring previous region settings
+ g.message("|! Original Region restored")
+
+ # print citation
+ if info:
+ print '\nSource: ' + citation_lst
+
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.landsat8.swlst/landsat8_mtl.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/landsat8_mtl.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/landsat8_mtl.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,273 @@
+#!/usr/bin/python\<nl>\
+# -*- coding: utf-8 -*-
+
+"""
+ at author nik |
+"""
+
+import sys
+from collections import namedtuple
+
+
+# globals
+MTLFILE = ''
+DUMMY_MAPCALC_STRING_RADIANCE = 'Radiance'
+DUMMY_MAPCALC_STRING_DN = 'DigitalNumber'
+
+
+# helper functions
+def set_mtlfile():
+ """
+ Set user defined MTL file, if any
+ """
+ if len(sys.argv) > 1:
+ return sys.argv[1]
+ else:
+ return False
+
+
+class Landsat8_MTL():
+ """
+ Retrieve metadata from a Landsat8 MTL file.
+ See <http://landsat.usgs.gov/Landsat8_Using_Product.php>.
+
+ ToDo:
+
+ - Implement toar_reflectance
+ - Implement mechanism to translate QA pixel values to QA bits, and vice
+ versa?
+ - Other Landsat8 related functions/algorithms?
+ """
+
+ def __init__(self, mtl_filename):
+ """
+ Initialise class object based on a Landsat8 MTL filename.
+ """
+ # read lines
+ with open(mtl_filename, 'r') as mtl_file:
+ mtl_lines = mtl_file.readlines()
+
+ # close and remove 'mtl_file'
+ mtl_file.close()
+ del(mtl_file)
+
+ # clean and convert MTL lines in to a named tuple
+ self.mtl = self._to_namedtuple(mtl_lines, 'metadata')
+ self._set_attributes()
+
+ # shorten LANDSAT_SCENE_ID, SENSOR_ID
+ self.scene_id = self.mtl.LANDSAT_SCENE_ID
+ self.sensor = self.mtl.SENSOR_ID
+
+ # bounding box related
+ self.corner_ul = (self.mtl.CORNER_UL_LAT_PRODUCT,
+ self.mtl.CORNER_UL_LON_PRODUCT)
+ self.corner_lr = (self.mtl.CORNER_LR_LAT_PRODUCT,
+ self.mtl.CORNER_LR_LON_PRODUCT)
+ self.corner_ul_projection = (self.mtl.CORNER_UL_PROJECTION_X_PRODUCT,
+ self.mtl.CORNER_UL_PROJECTION_Y_PRODUCT)
+ self.corner_lr_projection = (self.mtl.CORNER_LR_PROJECTION_X_PRODUCT,
+ self.mtl.CORNER_LR_PROJECTION_Y_PRODUCT)
+ self.cloud_cover = self.mtl.CLOUD_COVER
+
+ def _to_namedtuple(self, list_of_lines, name_for_tuple):
+ """
+ This function performs the following actions on the given
+ 'list_of_lines':
+ - excludes lines containing the strings 'GROUP' and 'END'
+ - removes whitespaces and doublequotes from strings
+ - converts list of lines in to a named tuple
+ """
+ import string
+
+ # exclude lines containing 'GROUP', 'END'
+ lines = [line.strip() for line in list_of_lines
+ if not any(x in line for x in ('GROUP', 'END'))]
+
+ # keep a copy, maybe useful?
+ self._mtl_lines = lines
+ del(list_of_lines)
+
+ # empty variables to hold values
+ field_names = []
+ field_values = []
+
+ # loop over lines, do some cleaning
+ for idx in range(len(lines)):
+
+ # split line in '='
+ line = lines[idx]
+ line_split = line.split('=')
+
+ # get field name & field value, clean whitespaces and "
+ field_name = line_split[0].strip()
+ field_names.append(field_name)
+ field_value = line_split[1].strip()
+ field_value = field_value.translate(string.maketrans("", "",), '"')
+ field_values.append(field_value)
+
+ # named tuple
+ named_tuple = namedtuple(name_for_tuple, field_names)
+
+ # return named tuple
+ return named_tuple(*field_values)
+
+ def _set_attributes(self):
+ """
+ Set all parsed field names and values, from the MTL file, fed to the
+ named tuple 'self.mtl', as attributes to the object.
+ """
+ for field in self.mtl._fields:
+ field_lowercase = field.lower()
+ field_value = getattr(self.mtl, field)
+ setattr(self, field_lowercase, field_value)
+
+ def __str__(self):
+ """
+ Return a string representation of the scene's id.
+ """
+ msg = 'Landsat8 scene ID:'
+ return msg + ' ' + self.scene_id
+
+ def _get_mtl_lines(self):
+ """
+ Return the "hidden" copy of the MTL lines before cleaning (lines
+ containing 'GROUP' or 'END' are though excluded).
+ """
+ return self._mtl_lines
+
+ def toar_radiance(self, bandnumber):
+ """
+ Note, this function returns a valid expression for GRASS GIS' r.mapcalc
+ raster processing module.
+
+ 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:
+
+ Lλ = ML * Qcal + AL
+
+ where:
+
+ - Lλ = TOA spectral radiance (Watts/( m2 * srad * μm))
+
+ - ML = Band-specific multiplicative rescaling factor from the metadata
+ (RADIANCE_MULT_BAND_x, where x is the band number)
+
+ - AL = Band-specific additive rescaling factor from the metadata
+ (RADIANCE_ADD_BAND_x, where x is the band number)
+
+ - Qcal = Quantized and calibrated standard product pixel values (DN)
+
+ Some code borrowed from
+ <https://github.com/micha-silver/grass-landsat8/blob/master/r.in.landsat8.py>
+ """
+ multiplicative_factor = getattr(self.mtl, ('RADIANCE_MULT_BAND_' +
+ str(bandnumber)))
+ # print "ML:", multiplicative_factor
+
+ additive_factor = getattr(self.mtl, 'RADIANCE_ADD_BAND_' +
+ str(bandnumber))
+ # print "AL:", additive_factor
+
+ formula = '{ML}*{DUMMY_DN} + {AL}'
+ mapcalc = formula.format(ML=multiplicative_factor,
+ DUMMY_DN=DUMMY_MAPCALC_STRING_DN,
+ AL=additive_factor)
+
+ return mapcalc
+
+ def toar_reflectance(self, bandnumber):
+ """
+ Note, this function returns a valid expression for GRASS GIS' r.mapcalc
+ raster processing module.
+
+ Conversion to TOA Reflectance OLI band data can also be converted to
+ TOA planetary reflectance using reflectance rescaling coefficients
+ provided in the product metadata file (MTL file). The following
+ equation is used to convert DN values to TOA reflectance for OLI data
+ as follows:
+
+ ρλ' = MρQcal + Aρ
+
+ where:
+
+ - ρλ' = TOA planetary reflectance, without correction for solar angle.
+ Note that ρλ' does not contain a correction for the sun angle.
+
+ - Mρ = Band-specific multiplicative rescaling factor from the metadata
+ (REFLECTANCE_MULT_BAND_x, where x is the band number)
+
+ - Aρ = Band-specific additive rescaling factor from the metadata
+ (REFLECTANCE_ADD_BAND_x, where x is the band number)
+
+ - Qcal = Quantized and calibrated standard product pixel values (DN)
+
+ TOA reflectance with a correction for the sun angle is then:
+
+ ρλ = ρλ' = ρλ' ### Fix This!
+ cos(θSZ) sin(θSE) ### Fix This!
+
+ where:
+
+ - ρλ = TOA planetary reflectance
+ - θSE = Local sun elevation angle. The scene center sun elevation angle
+ in degrees is provided in the metadata (SUN_ELEVATION).
+ - θSZ = Local solar zenith angle;
+ - θSZ = 90° - θSE
+
+ For more accurate reflectance calculations, per pixel solar angles
+ could be used instead of the scene center solar angle, but per pixel
+ solar zenith angles are not currently provided with the Landsat 8
+ products.
+ """
+ pass
+
+ def radiance_to_temperature(self, bandnumber):
+ """
+ Note, this function returns a valid expression for GRASS GIS' r.mapcalc
+ raster processing module.
+
+ 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:
+
+ T = K2 / ln( (K1/Lλ) + 1 )
+
+ where:
+
+ - T = At-satellite brightness temperature (K)
+
+ - Lλ = TOA spectral radiance (Watts/( m2 * srad * μm)), below
+ 'DUMMY_RADIANCE'
+
+ - K1 = Band-specific thermal conversion constant from the metadata
+ (K1_CONSTANT_BAND_x, where x is the band number, 10 or 11)
+
+ - K2 = Band-specific thermal conversion constant from the metadata
+ (K2_CONSTANT_BAND_x, where x is the band number, 10 or 11)
+ """
+ k2 = getattr(self.mtl, ('K2_CONSTANT_BAND_' + str(bandnumber)))
+ k1 = getattr(self.mtl, ('K1_CONSTANT_BAND_' + str(bandnumber)))
+
+ formula = '{K2} / ( log({K1} / {DUMMY_RADIANCE} + 1))'
+ mapcalc = formula.format(K2=k2,
+ K1=k1,
+ DUMMY_RADIANCE=DUMMY_MAPCALC_STRING_RADIANCE)
+
+ return mapcalc
+
+
+def main():
+ """
+ Main program.
+ """
+ if set_mtlfile():
+ MTLFILE = set_mtlfile()
+ print "| Reading metadata from:", MTLFILE
+ else:
+ MTLFILE = ''
+
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_11_detail_500px.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_11_detail_500px.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_13_detail_500px.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_13_detail_500px.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_7_detail_500px.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_7_detail_500px.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_9_detail_500px.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/lst_window_9_detail_500px.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/imagery/i.landsat8.swlst/map.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/imagery/i.landsat8.swlst/map.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/grass7/imagery/i.landsat8.swlst/mtl.txt
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/mtl.txt (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/mtl.txt 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,206 @@
+GROUP = L1_METADATA_FILE
+ GROUP = METADATA_FILE_INFO
+ ORIGIN = "Image courtesy of the U.S. Geological Survey"
+ REQUEST_ID = "0501405260188_00022"
+ LANDSAT_SCENE_ID = "LC81840332014146LGN00"
+ FILE_DATE = 2014-05-26T14:10:31Z
+ STATION_ID = "LGN"
+ PROCESSING_SOFTWARE_VERSION = "LPGS_2.3.0"
+ END_GROUP = METADATA_FILE_INFO
+ GROUP = PRODUCT_METADATA
+ DATA_TYPE = "L1T"
+ ELEVATION_SOURCE = "GLS2000"
+ OUTPUT_FORMAT = "GEOTIFF"
+ SPACECRAFT_ID = "LANDSAT_8"
+ SENSOR_ID = "OLI_TIRS"
+ WRS_PATH = 184
+ WRS_ROW = 33
+ NADIR_OFFNADIR = "NADIR"
+ TARGET_WRS_PATH = 184
+ TARGET_WRS_ROW = 33
+ DATE_ACQUIRED = 2014-05-26
+ SCENE_CENTER_TIME = 09:10:26.7368720Z
+ CORNER_UL_LAT_PRODUCT = 39.96125
+ CORNER_UL_LON_PRODUCT = 20.69909
+ CORNER_UR_LAT_PRODUCT = 39.93688
+ CORNER_UR_LON_PRODUCT = 23.39576
+ CORNER_LL_LAT_PRODUCT = 37.84718
+ CORNER_LL_LON_PRODUCT = 20.70789
+ CORNER_LR_LAT_PRODUCT = 37.82457
+ CORNER_LR_LON_PRODUCT = 23.32582
+ CORNER_UL_PROJECTION_X_PRODUCT = 474300.000
+ CORNER_UL_PROJECTION_Y_PRODUCT = 4423500.000
+ CORNER_UR_PROJECTION_X_PRODUCT = 704700.000
+ CORNER_UR_PROJECTION_Y_PRODUCT = 4423500.000
+ CORNER_LL_PROJECTION_X_PRODUCT = 474300.000
+ CORNER_LL_PROJECTION_Y_PRODUCT = 4188900.000
+ CORNER_LR_PROJECTION_X_PRODUCT = 704700.000
+ CORNER_LR_PROJECTION_Y_PRODUCT = 4188900.000
+ PANCHROMATIC_LINES = 15641
+ PANCHROMATIC_SAMPLES = 15361
+ REFLECTIVE_LINES = 7821
+ REFLECTIVE_SAMPLES = 7681
+ THERMAL_LINES = 7821
+ THERMAL_SAMPLES = 7681
+ FILE_NAME_BAND_1 = "LC81840332014146LGN00_B1.TIF"
+ FILE_NAME_BAND_2 = "LC81840332014146LGN00_B2.TIF"
+ FILE_NAME_BAND_3 = "LC81840332014146LGN00_B3.TIF"
+ FILE_NAME_BAND_4 = "LC81840332014146LGN00_B4.TIF"
+ FILE_NAME_BAND_5 = "LC81840332014146LGN00_B5.TIF"
+ FILE_NAME_BAND_6 = "LC81840332014146LGN00_B6.TIF"
+ FILE_NAME_BAND_7 = "LC81840332014146LGN00_B7.TIF"
+ FILE_NAME_BAND_8 = "LC81840332014146LGN00_B8.TIF"
+ FILE_NAME_BAND_9 = "LC81840332014146LGN00_B9.TIF"
+ FILE_NAME_BAND_10 = "LC81840332014146LGN00_B10.TIF"
+ FILE_NAME_BAND_11 = "LC81840332014146LGN00_B11.TIF"
+ FILE_NAME_BAND_QUALITY = "LC81840332014146LGN00_BQA.TIF"
+ METADATA_FILE_NAME = "LC81840332014146LGN00_MTL.txt"
+ BPF_NAME_OLI = "LO8BPF20140526084535_20140526092343.01"
+ BPF_NAME_TIRS = "LT8BPF20140526084141_20140526092436.01"
+ CPF_NAME = "L8CPF20140401_20140630.01"
+ RLUT_FILE_NAME = "L8RLUT20130211_20431231v09.h5"
+ END_GROUP = PRODUCT_METADATA
+ GROUP = IMAGE_ATTRIBUTES
+ CLOUD_COVER = 3.67
+ IMAGE_QUALITY_OLI = 9
+ IMAGE_QUALITY_TIRS = 9
+ ROLL_ANGLE = -0.001
+ SUN_AZIMUTH = 130.58922669
+ SUN_ELEVATION = 65.65724973
+ EARTH_SUN_DISTANCE = 1.0130385
+ GROUND_CONTROL_POINTS_MODEL = 521
+ GEOMETRIC_RMSE_MODEL = 5.833
+ GEOMETRIC_RMSE_MODEL_Y = 4.174
+ GEOMETRIC_RMSE_MODEL_X = 4.074
+ GROUND_CONTROL_POINTS_VERIFY = 187
+ GEOMETRIC_RMSE_VERIFY = 3.157
+ END_GROUP = IMAGE_ATTRIBUTES
+ GROUP = MIN_MAX_RADIANCE
+ RADIANCE_MAXIMUM_BAND_1 = 740.62347
+ RADIANCE_MINIMUM_BAND_1 = -61.16093
+ RADIANCE_MAXIMUM_BAND_2 = 758.40747
+ RADIANCE_MINIMUM_BAND_2 = -62.62954
+ RADIANCE_MAXIMUM_BAND_3 = 698.86597
+ RADIANCE_MINIMUM_BAND_3 = -57.71258
+ RADIANCE_MAXIMUM_BAND_4 = 589.32318
+ RADIANCE_MINIMUM_BAND_4 = -48.66650
+ RADIANCE_MAXIMUM_BAND_5 = 360.63638
+ RADIANCE_MINIMUM_BAND_5 = -29.78147
+ RADIANCE_MAXIMUM_BAND_6 = 89.68699
+ RADIANCE_MINIMUM_BAND_6 = -7.40638
+ RADIANCE_MAXIMUM_BAND_7 = 30.22932
+ RADIANCE_MINIMUM_BAND_7 = -2.49635
+ RADIANCE_MAXIMUM_BAND_8 = 666.95166
+ RADIANCE_MINIMUM_BAND_8 = -55.07708
+ RADIANCE_MAXIMUM_BAND_9 = 140.94489
+ RADIANCE_MINIMUM_BAND_9 = -11.63927
+ RADIANCE_MAXIMUM_BAND_10 = 22.00180
+ RADIANCE_MINIMUM_BAND_10 = 0.10033
+ RADIANCE_MAXIMUM_BAND_11 = 22.00180
+ RADIANCE_MINIMUM_BAND_11 = 0.10033
+ END_GROUP = MIN_MAX_RADIANCE
+ GROUP = MIN_MAX_REFLECTANCE
+ REFLECTANCE_MAXIMUM_BAND_1 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_1 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_2 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_2 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_3 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_3 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_4 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_4 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_5 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_5 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_6 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_6 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_7 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_7 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_8 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_8 = -0.099980
+ REFLECTANCE_MAXIMUM_BAND_9 = 1.210700
+ REFLECTANCE_MINIMUM_BAND_9 = -0.099980
+ END_GROUP = MIN_MAX_REFLECTANCE
+ GROUP = MIN_MAX_PIXEL_VALUE
+ QUANTIZE_CAL_MAX_BAND_1 = 65535
+ QUANTIZE_CAL_MIN_BAND_1 = 1
+ QUANTIZE_CAL_MAX_BAND_2 = 65535
+ QUANTIZE_CAL_MIN_BAND_2 = 1
+ QUANTIZE_CAL_MAX_BAND_3 = 65535
+ QUANTIZE_CAL_MIN_BAND_3 = 1
+ QUANTIZE_CAL_MAX_BAND_4 = 65535
+ QUANTIZE_CAL_MIN_BAND_4 = 1
+ QUANTIZE_CAL_MAX_BAND_5 = 65535
+ QUANTIZE_CAL_MIN_BAND_5 = 1
+ QUANTIZE_CAL_MAX_BAND_6 = 65535
+ QUANTIZE_CAL_MIN_BAND_6 = 1
+ QUANTIZE_CAL_MAX_BAND_7 = 65535
+ QUANTIZE_CAL_MIN_BAND_7 = 1
+ QUANTIZE_CAL_MAX_BAND_8 = 65535
+ QUANTIZE_CAL_MIN_BAND_8 = 1
+ QUANTIZE_CAL_MAX_BAND_9 = 65535
+ QUANTIZE_CAL_MIN_BAND_9 = 1
+ QUANTIZE_CAL_MAX_BAND_10 = 65535
+ QUANTIZE_CAL_MIN_BAND_10 = 1
+ QUANTIZE_CAL_MAX_BAND_11 = 65535
+ QUANTIZE_CAL_MIN_BAND_11 = 1
+ END_GROUP = MIN_MAX_PIXEL_VALUE
+ GROUP = RADIOMETRIC_RESCALING
+ RADIANCE_MULT_BAND_1 = 1.2235E-02
+ RADIANCE_MULT_BAND_2 = 1.2528E-02
+ RADIANCE_MULT_BAND_3 = 1.1545E-02
+ RADIANCE_MULT_BAND_4 = 9.7352E-03
+ RADIANCE_MULT_BAND_5 = 5.9575E-03
+ RADIANCE_MULT_BAND_6 = 1.4816E-03
+ RADIANCE_MULT_BAND_7 = 4.9937E-04
+ RADIANCE_MULT_BAND_8 = 1.1018E-02
+ RADIANCE_MULT_BAND_9 = 2.3283E-03
+ RADIANCE_MULT_BAND_10 = 3.3420E-04
+ RADIANCE_MULT_BAND_11 = 3.3420E-04
+ RADIANCE_ADD_BAND_1 = -61.17316
+ RADIANCE_ADD_BAND_2 = -62.64206
+ RADIANCE_ADD_BAND_3 = -57.72413
+ RADIANCE_ADD_BAND_4 = -48.67623
+ RADIANCE_ADD_BAND_5 = -29.78743
+ RADIANCE_ADD_BAND_6 = -7.40786
+ RADIANCE_ADD_BAND_7 = -2.49685
+ RADIANCE_ADD_BAND_8 = -55.08810
+ RADIANCE_ADD_BAND_9 = -11.64160
+ RADIANCE_ADD_BAND_10 = 0.10000
+ RADIANCE_ADD_BAND_11 = 0.10000
+ REFLECTANCE_MULT_BAND_1 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_2 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_3 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_4 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_5 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_6 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_7 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_8 = 2.0000E-05
+ REFLECTANCE_MULT_BAND_9 = 2.0000E-05
+ REFLECTANCE_ADD_BAND_1 = -0.100000
+ REFLECTANCE_ADD_BAND_2 = -0.100000
+ REFLECTANCE_ADD_BAND_3 = -0.100000
+ REFLECTANCE_ADD_BAND_4 = -0.100000
+ REFLECTANCE_ADD_BAND_5 = -0.100000
+ REFLECTANCE_ADD_BAND_6 = -0.100000
+ REFLECTANCE_ADD_BAND_7 = -0.100000
+ REFLECTANCE_ADD_BAND_8 = -0.100000
+ REFLECTANCE_ADD_BAND_9 = -0.100000
+ END_GROUP = RADIOMETRIC_RESCALING
+ GROUP = TIRS_THERMAL_CONSTANTS
+ K1_CONSTANT_BAND_10 = 774.89
+ K1_CONSTANT_BAND_11 = 480.89
+ K2_CONSTANT_BAND_10 = 1321.08
+ K2_CONSTANT_BAND_11 = 1201.14
+ END_GROUP = TIRS_THERMAL_CONSTANTS
+ GROUP = PROJECTION_PARAMETERS
+ MAP_PROJECTION = "UTM"
+ DATUM = "WGS84"
+ ELLIPSOID = "WGS84"
+ UTM_ZONE = 34
+ GRID_CELL_SIZE_PANCHROMATIC = 15.00
+ GRID_CELL_SIZE_REFLECTIVE = 30.00
+ GRID_CELL_SIZE_THERMAL = 30.00
+ ORIENTATION = "NORTH_UP"
+ RESAMPLING_OPTION = "CUBIC_CONVOLUTION"
+ END_GROUP = PROJECTION_PARAMETERS
+END_GROUP = L1_METADATA_FILE
+END
Added: grass-addons/grass7/imagery/i.landsat8.swlst/split_window_lst.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/split_window_lst.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/split_window_lst.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,745 @@
+# -*- coding: utf-8 -*-
+"""
+A class for the Split Window Algorithm for Land Surface Temperature estimation
+ at author: nik | Created on Wed Mar 18 11:28:45 2015
+"""
+
+# import average emissivities
+import random
+import csv_to_dictionary as coefficients
+from column_water_vapor import Column_Water_Vapor
+
+# globals
+EMISSIVITIES = coefficients.get_average_emissivities()
+COLUMN_WATER_VAPOR = coefficients.get_column_water_vapor()
+DUMMY_MAPCALC_STRING_T10 = 'Input_T10'
+DUMMY_MAPCALC_STRING_T11 = 'Input_T11'
+DUMMY_MAPCALC_STRING_AVG_LSE = 'Input_AVG_LSE'
+DUMMY_MAPCALC_STRING_DELTA_LSE = 'Input_DELTA_LSE'
+DUMMY_MAPCALC_STRING_FROM_GLC = 'Input_FROMGLC'
+DUMMY_MAPCALC_STRING_CWV = 'Input_CWV'
+
+# Remove from here, improve and use named tuples!
+FROM_GLC_CODES = [10, 11, 12, 13,
+ 20, 21, 22, 23, 24,
+ 30, 31, 32,
+ 51, 72,
+ 40, 71,
+ 60, 61, 62, 63,
+ 80, 81, 82,
+ 90, 52, 91, 92, 93, 94, 95, 96,
+ 100, 101, 102]
+
+FROM_GLC_LEGEND = {'Cropland': (10, 11, 12, 13),
+ 'Forest': (20, 21, 22, 23, 24),
+ 'Grasslands': (30, 31, 32, 51, 72),
+ 'Shrublands': (40, 71),
+ 'Waterbodies': (60, 61, 62, 63),
+ 'Tundra': (70,),
+ 'Impervious': (80, 81, 82),
+ 'Barren_Land': (90, 52, 91, 92, 93, 94, 95, 96),
+ 'Snow_and_ice': (100, 101, 102),
+ 'Cloud': (120,)}
+
+
+# helper functions
+def check_t1x_range(number):
+ """
+ Check if Brigthness Temperature (Kelvin degrees) values for T10, T11, lie
+ inside a reasonable range, eg [200, 330].
+
+ Note, the Digital Number values in bands B10 and B11, are 16-bit. The
+ actual data quantisation though is 12-bit.
+ """
+ if number < 200 or number > 330:
+ raise ValueError('The input value {t1x} for T1x is out of a '
+ 'reasonable range [200, 330]'.format(t1x=number))
+ else:
+ return True
+
+
+def check_cwv(cwv):
+ """
+ Check whether a column water value lies within a "valid" range. Which is?
+
+ - Questions to answer:
+
+ - What should happen when the value is out of range?
+ - Use subrange-6?
+ - If yes, how much tolerance for outliers < 0.0 and > 6.3 ? Testing
+ for +-.5
+
+ """
+ if cwv < 0.0 - .5 or cwv > 6.3 + .5:
+ raise ValueError('The column water vapor estimation is out of the '
+ 'expected range [0.0, 6.3]')
+ else:
+ return True
+
+
+class SplitWindowLST():
+ """
+ A class implementing the split-window algorithm for Landsat8 imagery
+
+ Inputs:
+
+ - The class itself requires only a string for 'landcover' which is:
+
+ 1) a fixed land cover class string (one from the classes defined in the
+ FROM-GLC legend)
+
+ 2) a land cover class code (integer) one from the classes defined in the
+ FROM-GLC classification scheme.
+
+ - Inputs for individual functions vary, look at their definitions.
+
+ Outputs:
+
+ - Valid expressions for GRASS GIS' r.mapcalc raster processing module
+ - Direct computation for... though not necessary, nor useful for GRASS GIS
+ modules directly?
+
+
+ Details
+
+ The algorithm removes the atmospheric effect through differential
+ atmospheric absorption in the two adjacent thermal infrared channels
+ centered at about 11 and 12 μm.
+
+ The linear or non-linear combination of the brightness temperatures is
+ finally applied for LST estimation based on the equation:
+
+ LST = b0 +
+ + ( b1 + b2 * ( 1 - ae ) / ae + b3 * de / ae^2 ) * ( t10 + t11 ) / 2 +
+ + ( b4 + b5 * ( 1 - ae ) / ae + b6 * de / ae^2 ) * ( t10 - t11 ) / 2 +
+ + b7 * ( t10 - t11 )^2
+
+ To reduce the influence of the CWV error on the LST, for a CWV within the
+ overlap of two adjacent CWV sub-ranges, the coefficients for the two
+ adjacent CWV sub-ranges are used to calculate the two initial temperatures.
+ Then, the average of these temperatures is assigned to as the pixel LST.
+
+ 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 δLSTinc and improves the spatial continuity of the LST product.
+ """
+
+ def __init__(self, landcover):
+ """
+ Create a class object for Split Window algorithm
+
+ Required inputs:
+ - B10
+ - B11 -- ToAR?
+ - land cover class?
+ - average emissivities for B10, B11
+ - subrange for column water vapor
+ """
+ # citation
+ self.citation = ('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.')
+
+ # basic equation (for __str__)
+ self._equation = ('[b0 + '
+ '(b1 + '
+ 'b2 * (1-ae) / ae + '
+ 'b3 * de / ae^2) * (t10 + t11) / 2 + '
+ '(b4 + '
+ 'b5 * (1-ae) / ae + '
+ 'b6 * de / ae^2) * (t10 - t11) / 2 + '
+ 'b7 * (t10 - t11)^2]')
+
+ # basic model (for... )
+ self._model = ('[{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]')
+
+ if landcover in EMISSIVITIES.keys() or landcover == 'Random':
+
+ # a fixed land cover class requested
+ self.landcover_class = landcover
+
+ # retrieve average emissivities for this class
+ emissivity_b10, emissivity_b11 = \
+ self._retrieve_average_emissivities(landcover)
+
+ # split avg emissivities of/for channels t10, t11
+ self.emissivity_t10 = float(emissivity_b10)
+ self.emissivity_t11 = float(emissivity_b11)
+
+ # average emissivity
+ self.average_emissivity = \
+ self._compute_average_emissivity(self.emissivity_t10,
+ self.emissivity_t11)
+
+ # delta emissivity
+ self.delta_emissivity = \
+ self._compute_delta_emissivity(self.emissivity_t10,
+ self.emissivity_t11)
+
+ else:
+
+ # if no fixed land cover class requested
+ self.landcover_class = False
+
+ # use mapcalc expressions instead, containing DUMMY strings for map
+ # names
+ self.average_lse_mapcalc = self._build_average_emissivity_mapcalc()
+ self.delta_lse_mapcalc = self._build_delta_emissivity_mapcalc()
+
+ # all-in-one split-window lst expression for mapcalc
+ self.sw_lst_mapcalc = self._build_swlst_mapcalc()
+
+ def __str__(self):
+ """
+ Return a string representation of the basic Split Window LST equation
+ """
+ equation = ' > The algorithm\'s basic equation: ' + self._equation
+
+ #if self.model:
+ # model = ' > The model: ' + self.model
+
+ #else:
+ # model = ' > The models:\n ' + ' a: ' + self._model_a + '\n ' + ' b: ' + self._model_b + '\n'
+
+ return equation #+ '\n' + model
+
+ def _landcover_string_validity(self, string):
+ """
+ Check whether the given string belongs to the list (keys) of known land
+ cover class names (to the FROM-GLC classification scheme) or is
+ identical to 'Random' and return, accordingly, True or False.
+ """
+ if string in FROM_GLC_LEGEND.keys():
+ return True
+ elif string == 'Random':
+ return True
+ else:
+ return False
+
+ def _retrieve_average_emissivities(self, emissivity_class):
+ """
+ Get land surface average emissivities from an emissivity look-up table.
+ This helper function returns a tuple.
+
+ Input is either one of the standard FROM-GLC land cover classes (...),
+ or one of its corresponding land cover class codes (...).
+
+ For testing purposes, the string "Random" is accepted to select a
+ random land surface emissivity class.
+ """
+
+ # land cover class code (integer)
+ if self.landcover_class and type(emissivity_class) == int:
+
+ assert self._landcover_string_validity(self.landcover_class), \
+ "Unknown land cover class name!"
+
+ if int(emissivity_class) in FROM_GLC_CODES:
+ landcover_code = int(emissivity_class)
+
+ # retrieving emissivity based on land cover class code
+ emissivity_class = [key for key in FROM_GLC_LEGEND
+ if landcover_code
+ in FROM_GLC_LEGEND[key]][0]
+
+ elif not int(emissivity_class) in FROM_GLC_LEGEND:
+ print ('The given land cover class code is not present in '
+ 'FROM-GLC map\'s legend!')
+
+ # random?
+ if type(emissivity_class) == str and emissivity_class == 'Random':
+ emissivity_class = random.choice(EMISSIVITIES.keys())
+ self.landcover_class = emissivity_class
+
+ # fields = EMISSIVITIES[emissivity_class]._fields
+ emissivity_b10 = EMISSIVITIES[emissivity_class].TIRS10
+ emissivity_b11 = EMISSIVITIES[emissivity_class].TIRS11
+
+ return (emissivity_b10, emissivity_b11)
+
+ def _compute_average_emissivity(self, emissivity_t10, emissivity_t11):
+ """
+ Return the average emissivity value for channels T10, T11.
+ """
+ average = float(0.5 * (emissivity_t10 + emissivity_t11))
+ return average
+
+ def _compute_delta_emissivity(self, emissivity_t10, emissivity_t11):
+ """
+ Return the difference of emissivity values for channels T10, T11.
+ """
+ delta = float(emissivity_t10 - emissivity_t11)
+ return delta
+
+ def _retrieve_adjacent_cwv_subranges(self, column_water_vapor):
+ """
+ Select and return adjacent subranges (string to be used as a dictionary
+ key) based on the atmospheric column water vapor estimation (float
+ ratio) ranging in (0.0, 6.3].
+
+ Input "cwv" is an estimation of the column water vapor (float ratio).
+ """
+ cwv = column_water_vapor
+ check_cwv(cwv) # check if float?
+
+ # a "subrange" generator
+ key_subrange_generator = ((key, COLUMN_WATER_VAPOR[key].subrange)
+ for key in COLUMN_WATER_VAPOR.keys())
+
+ # get all but the last -- using a list, after all!
+ subranges = list(key_subrange_generator)
+ # print " * Subranges:", subranges
+
+ # cwv in one or two subranges?
+ result = [range_x for range_x, (low, high) in subranges[:5]
+ if low < cwv < high]
+
+ # if one subrange, return a string
+ if len(result) == 1:
+ # self._cwv_subrange = result[0]
+ # self._cwv_subrange_a = self._cwv_subrange_b = False
+ return result[0]
+
+ # if two subranges, return a tuple
+ elif len(result) == 2:
+ # self._cwv_subrange = False
+ # self._cwv_subrange_a, self._cwv_subrange_b = tuple(result)
+ return result[0], result[1]
+
+ # what if it fails? -> subrange6
+ else:
+ # print " * Using the complete CWV range [0, 6.3]"
+ return subranges[5][0]
+
+ def _set_adjacent_cwv_subranges(self, column_water_vapor):
+ """
+ Set the retrieved cwv subranges as an attribute, though not a public
+ one.
+ """
+ result = self._retrieve_adjacent_cwv_subranges(column_water_vapor)
+ if len(result) == 1:
+ self._cwv_subrange = result[0]
+ self._cwv_subrange_a = self._cwv_subrange_b = False
+
+ elif len(result) == 2:
+ self._cwv_subrange = False
+ self._cwv_subrange_a, self._cwv_subrange_b = tuple(result)
+
+ # what to do for subrange6?
+
+ def _retrieve_cwv_coefficients(self, subrange):
+ """
+ Retrieve column water vapor coefficients for requested subrange
+ """
+ b0 = COLUMN_WATER_VAPOR[subrange].b0
+ b1 = COLUMN_WATER_VAPOR[subrange].b1
+ b2 = COLUMN_WATER_VAPOR[subrange].b2
+ b3 = COLUMN_WATER_VAPOR[subrange].b3
+ b4 = COLUMN_WATER_VAPOR[subrange].b4
+ b5 = COLUMN_WATER_VAPOR[subrange].b5
+ b6 = COLUMN_WATER_VAPOR[subrange].b6
+ b7 = COLUMN_WATER_VAPOR[subrange].b7
+
+ cwv_coefficients = (b0,
+ b1,
+ b2,
+ b3,
+ b4,
+ b5,
+ b6,
+ b7)
+
+ return cwv_coefficients
+
+ def _set_cwv_coefficients(self, subrange):
+ """
+ Set the coefficients as an attribute.
+ """
+ self.cwv_coefficients = self._retrieve_cwv_coefficients(subrange)
+
+ def get_cwv_coefficients(self):
+ """
+ Return the column water vapor coefficients from the object's attribute.
+ """
+ if self.cwv_coefficients:
+ return self.cwv_coefficients
+ else:
+ print "* CWV coefficients have not been set!"
+
+ def _retrieve_rmse(self, subrange):
+ """
+ Retrieve and set the associated RMSE for the column water vapor
+ coefficients for the subrange in question.
+ """
+ return COLUMN_WATER_VAPOR[subrange].rmse
+
+ def _set_rmse(self, subrange):
+ """
+ Set the RMSE as an attribute.
+ """
+ self.rmse = self._retrieve_rmse(subrange)
+
+ def report_rmse(self):
+ """
+ Report the associated R^2 value for the coefficients in question
+ """
+ msg = "Associated RMSE: "
+ return msg + str(self.rmse)
+
+ def compute_lst(self, t10, t11, coefficients):
+ """
+ Compute Land Surface Temperature based on the Split-Window algorithm.
+ Inputs are brightness temperatures measured in channels i(~11.0 μm)
+ and j (~12.0 μm).
+
+ *Note*, this is a single value computation function and does not read
+ or return a map.
+
+ LST = b0 +
+ + (b1 + b2 * ((1-ae)/ae) + b3 * (de/ae^2)) * ((t10 + t11)/2) +
+ + (b4 + b5 * ((1-ae)/ae) + b6 * (de/ae^2)) * ((t10 - t11)/2) +
+ + b7 * (t10 - t11)^2
+ """
+
+ # check validity of t10, t11
+ check_t1x_range(t10)
+ check_t1x_range(t11)
+
+ # check validity of subrange (string)?
+
+ # set cwv coefficients
+ b0, b1, b2, b3, b4, b5, b6, b7 = coefficients
+
+ # average and delta emissivity
+ avg = self.average_emissivity
+ delta = self.delta_emissivity
+
+ # addends
+ a = b0
+ b1 = b1 + b2 * (1-avg) / avg + b3 * delta / avg**2
+ b2 = (t10 + t11) / 2
+ b = b1 * b2
+ c1 = b4 + b5 * (1-avg) / avg + b6 * delta / avg**2
+ c2 = (t10 - t11) / 2
+ c = c1 * c2
+ d = b7 * (t10 - t11)**2
+
+ # land surface temperature
+ lst = a + b + c + d
+ return lst
+
+ def _set_lst(self):
+ """
+ Set the result of the single value computation function as an attribute
+ (?)
+ """
+ pass
+
+ def compute_average_lst(self, t10, t11, subrange_a, subrange_b):
+ """
+ Compute average LST, similar as the compute_lst function.
+ """
+
+ # retrieve coefficients for first subrange and compute lst for it
+ coefficients_a = self._retrieve_cwv_coefficients(subrange_a)
+ lst_a = self.compute_lst(t10, t11, coefficients_a)
+
+ # repeat for second subrange
+ coefficients_b = self._retrieve_cwv_coefficients(subrange_b)
+ lst_b = self.compute_lst(t10, t11, coefficients_b)
+
+ # average land surface temperature
+ return (lst_a + lst_b) / 2
+
+ def _set_average_lst(self):
+ """
+ Set the average LST pixel value as an attribute (?)
+ """
+ pass
+
+ def _build_average_emissivity_mapcalc(self):
+ """
+ ToDo: shorten the following!
+ """
+ # average emissivities
+ e10_t10, e10_t11 = self._retrieve_average_emissivities('Cropland')
+ avg_e10 = self._compute_average_emissivity(e10_t10, e10_t11)
+
+ e20_t10, e20_t11 = self._retrieve_average_emissivities('Forest')
+ avg_e20 = self._compute_average_emissivity(e20_t10, e20_t11)
+
+ e30_t10, e30_t11 = self._retrieve_average_emissivities('Grasslands')
+ avg_e30 = self._compute_average_emissivity(e30_t10, e30_t11)
+
+ e40_t10, e40_t11 = self._retrieve_average_emissivities('Shrublands')
+ avg_e40 = self._compute_average_emissivity(e40_t10, e40_t11)
+
+ e60_t10, e60_t11 = self._retrieve_average_emissivities('Waterbodies')
+ avg_e60 = self._compute_average_emissivity(e60_t10, e60_t11)
+
+ e80_t10, e80_t11 = self._retrieve_average_emissivities('Impervious')
+ avg_e80 = self._compute_average_emissivity(e80_t10, e80_t11)
+
+ e90_t10, e90_t11 = self._retrieve_average_emissivities('Barren_Land')
+ avg_e90 = self._compute_average_emissivity(e90_t10, e90_t11)
+
+ e100_t10, e100_t11 = self._retrieve_average_emissivities('Cropland')
+ avg_e100 = self._compute_average_emissivity(e100_t10, e100_t11)
+
+ expression = ('eval( class_10 = {landcover} >= 10 && {landcover} < 20,'
+ '\ \n class_20 = {landcover} >= 20 && {landcover} < 30,'
+ '\ \n class_30 = {landcover} == 72 || {landcover} >= 30 && {landcover} < 40,'
+ '\ \n class_40 = {landcover} >= 40 && {landcover} < 50,'
+ '\ \n class_50 = {landcover} >= 50 && {landcover} < 52,'
+ '\ \n class_60 = {landcover} >= 60 && {landcover} < 70,'
+ '\ \n class_70 = {landcover} >= 70 && {landcover} < 72,'
+ '\ \n class_80 = {landcover} >= 80 && {landcover} < 90,'
+ '\ \n class_90 = {landcover} == 52 || {landcover} >= 90 && {landcover} < 100,'
+ '\ \n class_100 = {landcover} >= 100 && {landcover} < 120,'
+ '\ \n if( class_10, {average_10},'
+ '\ \n if( class_20, {average_20},'
+ '\ \n if( class_30, {average_30},'
+ '\ \n if( class_40, {average_40},'
+ '\ \n if( class_50, {average_60},'
+ '\ \n if( class_60, {average_60},'
+ '\ \n if( class_70, {average_40},'
+ '\ \n if( class_80, {average_80},'
+ '\ \n if( class_90, {average_90},'
+ '\ \n if( class_100, {average_100},'
+ ' null() )))))))))))')
+
+ mapcalc = expression.format(landcover=DUMMY_MAPCALC_STRING_FROM_GLC,
+ average_10=avg_e10,
+ average_20=avg_e20,
+ average_30=avg_e30,
+ average_40=avg_e40,
+ average_60=avg_e60,
+ average_80=avg_e80,
+ average_90=avg_e90,
+ average_100=avg_e100)
+
+ return mapcalc
+
+ def _build_delta_emissivity_mapcalc(self):
+ """
+ ToDo: shorten the following!
+ """
+ # delta emissivities
+ e10_t10, e10_t11 = self._retrieve_average_emissivities('Cropland')
+ delta_e10 = self._compute_delta_emissivity(e10_t10, e10_t11)
+
+ e20_t10, e20_t11 = self._retrieve_average_emissivities('Forest')
+ delta_e20 = self._compute_delta_emissivity(e20_t10, e20_t11)
+
+ e30_t10, e30_t11 = self._retrieve_average_emissivities('Grasslands')
+ delta_e30 = self._compute_delta_emissivity(e30_t10, e30_t11)
+
+ e40_t10, e40_t11 = self._retrieve_average_emissivities('Shrublands')
+ delta_e40 = self._compute_delta_emissivity(e40_t10, e40_t11)
+
+ e60_t10, e60_t11 = self._retrieve_average_emissivities('Waterbodies')
+ delta_e60 = self._compute_delta_emissivity(e60_t10, e60_t11)
+
+ e80_t10, e80_t11 = self._retrieve_average_emissivities('Impervious')
+ delta_e80 = self._compute_delta_emissivity(e80_t10, e80_t11)
+
+ e90_t10, e90_t11 = self._retrieve_average_emissivities('Barren_Land')
+ delta_e90 = self._compute_delta_emissivity(e90_t10, e90_t11)
+
+ e100_t10, e100_t11 = self._retrieve_average_emissivities('Cropland')
+ delta_e100 = self._compute_delta_emissivity(e100_t10, e100_t11)
+
+ expression = ('eval( class_10 = {landcover} >= 10 && {landcover} < 20,'
+ '\ \n class_20 = {landcover} >= 20 && {landcover} < 30,'
+ '\ \n class_30 = {landcover} == 72 || {landcover} >= 30 && {landcover} < 40,'
+ '\ \n class_40 = {landcover} >= 40 && {landcover} < 50,'
+ '\ \n class_50 = {landcover} >= 50 && {landcover} < 52,'
+ '\ \n class_60 = {landcover} >= 60 && {landcover} < 70,'
+ '\ \n class_70 = {landcover} >= 70 && {landcover} < 72,'
+ '\ \n class_80 = {landcover} >= 80 && {landcover} < 90,'
+ '\ \n class_90 = {landcover} == 52 || {landcover} >= 90 && {landcover} < 100,'
+ '\ \n class_100 = {landcover} >= 100 && {landcover} < 120,'
+ '\ \n if( class_10, {delta_10},'
+ '\ \n if( class_20, {delta_20},'
+ '\ \n if( class_30, {delta_30},'
+ '\ \n if( class_40, {delta_40},'
+ '\ \n if( class_50, {delta_60},'
+ '\ \n if( class_60, {delta_60},'
+ '\ \n if( class_70, {delta_40},'
+ '\ \n if( class_80, {delta_80},'
+ '\ \n if( class_90, {delta_90},'
+ '\ \n if( class_100, {delta_100},'
+ ' null() )))))))))))')
+
+ mapcalc = expression.format(landcover=DUMMY_MAPCALC_STRING_FROM_GLC,
+ delta_10=delta_e10,
+ delta_20=delta_e20,
+ delta_30=delta_e30,
+ delta_40=delta_e40,
+ delta_60=delta_e60,
+ delta_80=delta_e80,
+ delta_90=delta_e90,
+ delta_100=delta_e100)
+
+ return mapcalc
+
+ def _build_model(self, coefficients):
+ """
+ Build model for __str__
+ """
+ b0, b1, b2, b3, b4, b5, b6, b7 = coefficients
+ model = self._model.format(b0=b0,
+ b1=b1,
+ b2=b2,
+ ae=self.average_emissivity,
+ de=self.delta_emissivity,
+ b3=b3,
+ b4=b4,
+ b5=b5,
+ b6=b6,
+ b7=b7,
+ t10=self.emissivity_t10,
+ t11=self.emissivity_t11)
+ return model
+
+ def _build_subrange_mapcalc(self, subrange):
+ """
+ Build formula for GRASS GIS' mapcalc for the given cwv subrange.
+
+ ToDo: Review and Improve the mechanism which selects emissivities from
+ either a fixed land cover class OR a land cover map.
+ """
+ # formula = '{c0} + {c1}*{dummy} + {c2}*{dummy}^2'
+ formula = ('{b0} + '
+ '({b1} + '
+ '({b2}) * ((1 - {ae}) / {ae}^2) + '
+ '({b3}) * ({de}/{ae}^2)) * (({DUMMY_T10}+{DUMMY_T11})/2) + '
+ '({b4} + '
+ '({b5}) * ((1 - {ae}) / {ae}) + '
+ '({b6}) * ({de}/{ae}^2)) * (({DUMMY_T10}-{DUMMY_T11})/2) + '
+ '({b7}) * ({DUMMY_T10} - {DUMMY_T11})^2')
+
+ # Implement mechanism to either select
+
+ try:
+ if self.landcover_class:
+ # print "Fixed land cover class"
+ emissivity_t10 = float(self.emissivity_t10)
+ emissivity_t11 = float(self.emissivity_t11)
+ avg_lse = self._compute_delta_emissivity(emissivity_t10,
+ emissivity_t11)
+ delta_lse = \
+ self._compute_delta_emissivity(emissivity_t10,
+ emissivity_t11)
+ except:
+ pass
+
+ if not self.landcover_class:
+
+ # This is required for when a fixed emissivity_class is used,
+ # instead of a FROM-GLC (landcover) map.
+
+ # print "Using the FROM-GLC map"
+ avg_lse = DUMMY_MAPCALC_STRING_AVG_LSE
+ delta_lse = DUMMY_MAPCALC_STRING_DELTA_LSE
+
+ coefficients = self._retrieve_cwv_coefficients(subrange)
+ b0, b1, b2, b3, b4, b5, b6, b7 = coefficients
+
+ mapcalc = formula.format(b0=b0,
+ b1=b1,
+ b2=b2,
+ ae=avg_lse,
+ de=delta_lse,
+ b3=b3,
+ b4=b4,
+ b5=b5,
+ b6=b6,
+ b7=b7,
+ DUMMY_T10=DUMMY_MAPCALC_STRING_T10,
+ DUMMY_T11=DUMMY_MAPCALC_STRING_T11)
+
+ return mapcalc
+
+ def _build_swlst_mapcalc(self):
+ """
+ Build and return a valid expression for GRASS GIS' r.mapcalc to
+ determine LST.
+ """
+ # subrange limits, low, high
+ low_1, high_1 = COLUMN_WATER_VAPOR['Range_1'].subrange
+ low_2, high_2 = COLUMN_WATER_VAPOR['Range_2'].subrange
+ low_3, high_3 = COLUMN_WATER_VAPOR['Range_3'].subrange
+ low_4, high_4 = COLUMN_WATER_VAPOR['Range_4'].subrange
+ low_5, high_5 = COLUMN_WATER_VAPOR['Range_5'].subrange
+ low_6, high_6 = COLUMN_WATER_VAPOR['Range_6'].subrange # unused
+
+ # build mapcalc expression for each subrange
+ expression_range_1 = self._build_subrange_mapcalc('Range_1')
+ expression_range_2 = self._build_subrange_mapcalc('Range_2')
+ expression_range_3 = self._build_subrange_mapcalc('Range_3')
+ expression_range_4 = self._build_subrange_mapcalc('Range_4')
+ expression_range_5 = self._build_subrange_mapcalc('Range_5')
+
+ # complete range
+ expression_range_6 = self._build_subrange_mapcalc('Range_6')
+
+ # build one big expression using mighty eval
+ expression = ('eval( sw_lst_1 = {exp_1},'
+ '\ \n sw_lst_2 = {exp_2},'
+ '\ \n sw_lst_12 = (sw_lst_1 + sw_lst_2) / 2,'
+ '\ \n sw_lst_3 = {exp_3},'
+ '\ \n sw_lst_23 = (sw_lst_2 + sw_lst_3) / 2,'
+ '\ \n sw_lst_4 = {exp_4},'
+ '\ \n sw_lst_34 = (sw_lst_3 + sw_lst_4) / 2,'
+ '\ \n sw_lst_5 = {exp_5},'
+ '\ \n sw_lst_45 = (sw_lst_4 + sw_lst_5) / 2,'
+ '\ \n sw_lst_6 = {exp_6},'
+ '\ \n in_range_1 = {low_1} < {DUMMY_CWV} && {DUMMY_CWV} < {high_1},'
+ '\ \n in_range_2 = {low_2} < {DUMMY_CWV} && {DUMMY_CWV} < {high_2},'
+ '\ \n in_range_3 = {low_3} < {DUMMY_CWV} && {DUMMY_CWV} < {high_3},'
+ '\ \n in_range_4 = {low_4} < {DUMMY_CWV} && {DUMMY_CWV} < {high_4},'
+ '\ \n in_range_5 = {low_5} < {DUMMY_CWV} && {DUMMY_CWV} < {high_5},'
+ '\ \n if( in_range_1 && in_range_2, sw_lst_12,'
+ '\ \n if( in_range_2 && in_range_3, sw_lst_23,'
+ '\ \n if( in_range_3 && in_range_4, sw_lst_34,'
+ '\ \n if( in_range_4 && in_range_5, sw_lst_45,'
+ '\ \n if( in_range_1, sw_lst_1,'
+ '\ \n if( in_range_2, sw_lst_2,'
+ '\ \n if( in_range_3, sw_lst_3,'
+ '\ \n if( in_range_4, sw_lst_4,'
+ '\ \n if( in_range_5, sw_lst_5,'
+ ' sw_lst_6 ))))))))))') # ' null() ))))))))))')
+
+ # replace keywords appropriately
+ swlst_expression = expression.format(exp_1=expression_range_1,
+ low_1=low_1,
+ DUMMY_CWV=DUMMY_MAPCALC_STRING_CWV,
+ high_1=high_1,
+ exp_2=expression_range_2,
+ low_2=low_2, high_2=high_2,
+ exp_3=expression_range_3,
+ low_3=low_3, high_3=high_3,
+ exp_4=expression_range_4,
+ low_4=low_4, high_4=high_4,
+ exp_5=expression_range_5,
+ low_5=low_5, high_5=high_5,
+ exp_6=expression_range_6)
+
+ return swlst_expression
+
+# reusable & stand-alone
+if __name__ == "__main__":
+ print ('Split-Window Algorithm for Estimating Land Surface Temperature '
+ 'from Landsat8 OLI/TIRS imagery.'
+ ' (Running as stand-alone tool?)\n')
Added: grass-addons/grass7/imagery/i.landsat8.swlst/test_column_water_vapor.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/test_column_water_vapor.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/test_column_water_vapor.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,116 @@
+#!/usr/bin/python\<nl>\
+# -*- coding: utf-8 -*-
+
+"""
+ at author nik |
+"""
+
+# required librairies
+import random
+from column_water_vapor import *
+
+
+# helper functions
+def random_adjacent_pixel_values(pixel_modifiers):
+ """
+ """
+ return [random.randint(250, 350) for dummy_idx in
+ range(len(pixel_modifiers))]
+
+
+def test_column_water_vapor():
+
+ print 'Equations for Column Water Vapor retrieval based on...'
+ print
+
+ print (' * Considering N adjacent pixels, the CWV in the MSWCVR method '
+ 'is estimated as:')
+ print
+ print (' cwv = c0 + c1 * (tj / ti) + c2 * (tj / ti)^2')
+ print
+ print (' where:\n\n - tj/ti ~ Rji = '
+ 'SUM [ ( Tik - Ti_mean ) * ( Tjk - Tj_mean ) ] / '
+ 'SUM [ ( Tik - Tj_mean )^2 ]')
+
+ print
+ print "Testing the Column_Water_Vapor class"
+ print
+
+ obj = Column_Water_Vapor(3, 'A', 'B')
+ print " | Testing the '__str__' method:\n\n ", obj
+ print
+
+ print " | Adjacent pixels:", obj.adjacent_pixels
+ print
+
+ print " | Map Ti:", obj.ti
+
+ print " | Map Tj:", obj.tj
+ print
+
+ print " | N pixels window modifiers for Ti:", obj.modifiers_ti
+
+ print " | N pixels window Modifiers for Tj:", obj.modifiers_tj
+
+ print " | Zipped modifiers_tij (used in a function for the Ratio ji):",
+ print obj.modifiers
+ print
+
+ print " ~ Random N pixel values for Ti:",
+ random_ti_values = random_adjacent_pixel_values(obj.modifiers_ti)
+ print random_ti_values
+
+ print " ~ Random N pixel values for Tj:",
+ random_tj_values = random_adjacent_pixel_values(obj.modifiers_ti)
+ print random_tj_values
+
+ print (' ~ Testing "compute_column_water_vapor" '
+ 'based on the above random values):'),
+ print obj.compute_column_water_vapor(random_ti_values, random_tj_values)
+ print
+
+ print " | Expression for Ti mean:", obj.mean_ti_expression
+
+ print " | Expression for Tj mean:", obj.mean_tj_expression
+ print
+
+ print " ~ Mean of random N pixel values for Ti:",
+ print sum(random_ti_values) / len(random_ti_values)
+
+ print " ~ Mean of random N pixel values for Tj:",
+ print sum(random_tj_values) / len(random_tj_values)
+ print
+
+ print " | Note, the following mapcalc expressions use dummy strings, meant to be replaced in the main program by the names of the maps in question"
+ print
+
+ print " | Expression for Numerator for Ratio (method):", obj._numerator_for_ratio('Ti_Mean', 'Tj_Mean')
+ print
+
+ print (' ~ Add example for Numerator based on '
+ 'Mean of random N pixel values for Ti, Tj:')
+ print
+
+
+ print " | Expression for Denominator for Ratio (method):", obj._denominator_for_ratio('Ti_Mean')
+ print
+
+ print (' ~ Add example for Denominator based on '
+ 'Mean of random N pixel values for Tj:')
+ print
+
+ print " | Ratio ji expression for mapcalc:", obj.ratio_ji_expression
+ print
+
+ print (' ~ Add example for Ratio ji expression based on '
+ 'numerator and denominator as defined above:')
+ print
+
+ print " | One big mapcalc expression:\n\n", obj._big_cwv_expression()
+ print
+
+# reusable & stand-alone
+if __name__ == "__main__":
+ print ('Testing the SplitWindowLST class')
+ print
+ test_column_water_vapor()
Added: grass-addons/grass7/imagery/i.landsat8.swlst/test_landsat8_mtl.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/test_landsat8_mtl.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/test_landsat8_mtl.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,47 @@
+#!/usr/bin/python\<nl>\
+# -*- coding: utf-8 -*-
+
+"""
+ at author nik |
+"""
+
+from landsat8_mtl import Landsat8_MTL
+
+
+MTLFILE = 'mtl.txt'
+
+
+def test(mtlfile):
+ """
+ Test the class.
+ """
+
+ print "! No file defined, testing with default MTl file!"
+ mtl = Landsat8_MTL(MTLFILE)
+ print
+
+ print "| Test the object's __str__ method:", mtl
+ print "| Test method _get_mtl_lines:\n ", mtl._get_mtl_lines()
+ print
+
+ print "| Basic metadata:"
+ print " > id (original field name is 'LANDSAT_SCENE_ID'):", mtl.scene_id
+ print " > WRS path:", mtl.wrs_path
+ print " > WRS row:", mtl.wrs_row
+ print " > Acquisition date:", mtl.date_acquired
+ print " > Scene's center time:", mtl.scene_center_time
+ print " > Upper left corner:", mtl.corner_ul
+ print " > Lower right corner:", mtl.corner_lr
+ print " > Upper left (projected):", mtl.corner_ul_projection
+ print " > Lower right (projected):", mtl.corner_lr_projection
+ print " > Cloud cover:", mtl.cloud_cover
+
+
+def main():
+ """
+ Main program.
+ """
+ test(MTLFILE)
+
+if __name__ == "__main__":
+ main()
Added: grass-addons/grass7/imagery/i.landsat8.swlst/test_split_window_lst.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/test_split_window_lst.py (rev 0)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/test_split_window_lst.py 2016-04-20 13:47:02 UTC (rev 68288)
@@ -0,0 +1,307 @@
+#!/usr/bin/python\<nl>\
+# -*- coding: utf-8 -*-
+
+"""
+ at author nik | 2015-04-18 03:48:20
+"""
+
+# required librairies
+import random
+import csv_to_dictionary as coefficients
+from split_window_lst import *
+
+# globals
+EMISSIVITIES = coefficients.get_average_emissivities()
+COLUMN_WATER_VAPOR = coefficients.get_column_water_vapor()
+
+# helper function(s)
+def random_digital_numbers(count=2):
+ """
+ Return a user-requested amount of random Digital Number values for testing
+ purposes ranging in 12-bit
+ """
+ digital_numbers = []
+
+ for dn in range(0, count):
+ digital_numbers.append(random.randint(1, 2**12))
+
+ if count == 1:
+ return digital_numbers[0]
+
+ return digital_numbers
+
+
+def random_brightness_temperature_values(count=2):
+ """
+ Return a user-requested amount of random Brightness Temperature values for testing
+ purposes ranging in [250, 330].
+ """
+ digital_numbers = []
+
+ for dn in range(0, count):
+ digital_numbers.append(random.randint(250, 330))
+
+ if count == 1:
+ return digital_numbers[0]
+
+ return digital_numbers
+
+def random_column_water_vapor():
+ """
+ Return a rational number ranging in [0.0, 6.3] to assisst in selecting
+ an atmospheric column water vapor subrange, as part of testing the
+ Split-WindowLST class.
+ """
+ return random.uniform(0.0 -1, 6.3 + 1)
+
+
+def test_split_window_lst():
+ """
+ Testing the SplitWindowLST class
+ """
+
+ #
+ # Helpers
+ #
+
+ print
+ print ">>> [Helper functions]"
+ print
+ t10 = random_brightness_temperature_values(1)
+ t11 = random.choice(((t10 + 50), (t10 - 50))) # check some failures as well
+ print " * Random brightness temperature values for T10, T11:", t10, t11
+ print (' * NOTE: Some out of a reasonable range T10, T11 values, which '
+ 'cause the current test to fail, are tolerated on-purpose in order '
+ 'to test the range checking function `check_t1x_range()`.')
+ print
+ print
+
+ #
+ # EMISSIVITIES
+ #
+
+ # get emissivities
+ print "[ EMISSIVITIES ]"
+ print
+ EMISSIVITIES = coefficients.get_average_emissivities()
+ print " * Dictionary for average emissivities:\n\n", EMISSIVITIES
+ print
+
+ somekey = random.choice(EMISSIVITIES.keys())
+ print " * Some random key from EMISSIVITIES:", somekey
+
+ fields = EMISSIVITIES[somekey]._fields
+ print " * Fields of namedtuple:", fields
+
+ random_field = random.choice(fields)
+ print " * Some random field:", random_field
+
+ command = 'EMISSIVITIES.[{key}].{field} ='
+ command = command.format(key=somekey, field=random_field)
+ print " * Example of retrieving values (named tuple): " + command,
+ print EMISSIVITIES[somekey].TIRS10, EMISSIVITIES[somekey].TIRS11
+
+ emissivity_b10 = EMISSIVITIES[somekey].TIRS10
+ print " * Average emissivity for B10:", emissivity_b10, "|Type:", type(emissivity_b10)
+ emissivity_b11 = EMISSIVITIES[somekey].TIRS11
+ print " * Average emissivity for B11:", emissivity_b11
+ print
+ print
+ print
+
+ #
+ # COLUMN_WATER_VAPOR
+ #
+
+ print "[ COLUMN_WATER_VAPOR ]"
+ print
+ print (' * NOTE: Some out of range values which cause the current test to '
+ 'fail, are tolerated on-purpose in order to check for the CWV '
+ 'range checking function. Check for the range of the '
+ 'random_column_water_vapor() function.')
+
+ COLUMN_WATER_VAPOR = coefficients.get_column_water_vapor()
+ print "\n * Dictionary for column water vapor coefficients:\n\n",
+ print COLUMN_WATER_VAPOR
+ print
+
+ cwvobj = Column_Water_Vapor(3, 'TiRS10', 'TiRS11')
+ print " * Retrieval of column water vapor via class, example: ", #, cwvobj
+ print " cwvobj = Column_Water_Vapour(3, MapName_for_T10, MapName_for_T11)"
+ print " * Mapcalc expression for it: ", cwvobj.column_water_vapor_expression
+ print
+
+ cwv = random_column_water_vapor()
+ print " * For the test, some random atmospheric column water vapor '(g/cm^2):", cwv
+
+ #
+ cwv_range_x = random.choice([key for key in COLUMN_WATER_VAPOR.keys()])
+
+ cwvfields = COLUMN_WATER_VAPOR[cwv_range_x]._fields
+ print " * Fields of namedtuple (same for all subranges):", cwvfields
+
+ random_cwvfield = random.choice(cwvfields)
+ print " * Some random field:", random_cwvfield
+
+ command = 'COLUMN_WATER_VAPOR.[{key}].{field} ='
+ command = command.format(key=cwv_range_x, field=random_cwvfield)
+ print " * Example of retrieving values (named tuple): " + command,
+ print COLUMN_WATER_VAPOR[cwv_range_x].subrange
+ print
+ print
+ print
+
+ #
+ # class
+ #
+
+ print "[ class SplitWindowLST ]"
+ print
+
+ # cwv_range_x = column_water_vapor_range(cwv)
+ # print " * Atmospheric column water vapor range:", cwv_range_x
+
+ swlst = SplitWindowLST(somekey) # somekey generated above
+ print "Create object and test '__str__' of SplitWindowLST() class:\n\n", swlst
+
+ print " > The 'citation' attribute:", swlst.citation
+ print
+
+ # get a column water vapor subrange
+ cwv_range_x = swlst._retrieve_adjacent_cwv_subranges(cwv)
+
+ # special case for Subrange 6
+ if cwv_range_x == 'Range_6':
+ msg = (' * The CWV value {cwv} falls outside of one of the subranges. '
+ 'Using the complete CWV range [0.0, 6.3] described as')
+ msg = msg.format(cwv=cwv)
+ print msg, cwv_range_x
+ else:
+ print " * The CWV value {cwv} falls inside:".format(cwv=cwv), cwv_range_x, "|Type:", type(cwv_range_x)
+
+### First, test all _retrieve functions ###
+
+ print
+ print "( Testing unpacking of values )"
+ print
+
+ if type(cwv_range_x) == str:
+
+ cwv_coefficients_x = swlst._retrieve_cwv_coefficients(cwv_range_x)
+ b0, b1, b2, b3, b4, b5, b6, b7 = cwv_coefficients_x
+
+ print " * Column Water Vapor coefficients (b0, b1, ..., b7) in <", cwv_range_x,
+
+ print "> :", b0, b1, b2, b3, b4, b5, b6, b7
+
+ print " * Model:", swlst._build_model(cwv_coefficients_x)
+
+ print " * Checking the '_retrieve_rmse' method:", swlst._retrieve_rmse(cwv_range_x)
+
+ print " * Testing the '_set_rmse' and 'report_rmse' methods:",
+
+ swlst._set_rmse(cwv_range_x)
+ print swlst.report_rmse()
+
+ print " * Testing the 'compute_lst' method:", swlst.compute_lst(t10, t11, cwv_coefficients_x)
+
+
+ elif type(cwv_range_x) == tuple and len(cwv_range_x) == 2:
+
+ print " * Two subranges returned:",
+
+ cwv_subrange_a, cwv_subrange_b = tuple(cwv_range_x)[0], tuple(cwv_range_x)[1]
+ print " Subrange a:", cwv_subrange_a, "| Subrange b:", cwv_subrange_b
+
+ #
+ # Subrange A
+ #
+
+ print
+ print " > Tests for subrange a"
+ print
+
+ coefficients_a = swlst._retrieve_cwv_coefficients(cwv_subrange_a)
+
+ b0, b1, b2, b3, b4, b5, b6, b7 = coefficients_a
+
+ print " * Column Water Vapor coefficients for", cwv_subrange_a,
+
+ print "> ", b0, b1, b2, b3, b4, b5, b6, b7
+
+ print " * Testing the '_set' and 'get' methods:",
+ swlst._set_cwv_coefficients(cwv_subrange_a) # does not return anything
+
+ print swlst.get_cwv_coefficients()
+
+ print " * Model:", swlst._build_model(coefficients_a)
+
+ print " * Checking the '_retrieve_rmse' method:", swlst._retrieve_rmse(cwv_subrange_a)
+ print " * Testing the '_set_rmse' and 'report_rmse' methods:",
+ swlst._set_rmse(cwv_subrange_a)
+
+ print swlst.report_rmse()
+
+
+
+ #
+ # Subrange B
+ #
+
+ print
+ print " > Tests for subrange b"
+ print
+
+ coefficients_b = swlst._retrieve_cwv_coefficients(cwv_subrange_b)
+
+ b0, b1, b2, b3, b4, b5, b6, b7 = coefficients_b
+
+ print " * Column Water Vapor coefficients for", cwv_subrange_b,
+
+ print "> ", b0, b1, b2, b3, b4, b5, b6, b7
+
+ print " * Testing the 'get' and '_set' methods:",
+ swlst._set_cwv_coefficients(cwv_subrange_b)
+
+ print swlst.get_cwv_coefficients()
+
+ print " * Model:", swlst._build_model(coefficients_b)
+
+ print " * Checking the '_retrieve_rmse' method:", swlst._retrieve_rmse(cwv_subrange_a)
+ print " * Testing the '_set_rmse' and 'report_rmse' methods:",
+ swlst._set_rmse(cwv_subrange_a)
+
+ print swlst.report_rmse()
+
+
+
+ #
+ # Average LST
+ #
+
+ print
+ print "( Computing average LST )"
+ print
+
+ print " * Compute the average LST: 'compute_average_lst()' >>>",
+
+ print swlst.compute_average_lst(t10, t11, cwv_subrange_a, cwv_subrange_b)
+
+ print
+
+ print
+ print "[ Subranges ]"
+ print
+
+ key_subrange_generator = ((key, COLUMN_WATER_VAPOR[key].subrange) for key in COLUMN_WATER_VAPOR.keys())
+
+
+ sw_lst_expression = swlst.sw_lst_mapcalc
+ print "Big expression:\n\n", sw_lst_expression
+
+
+# reusable & stand-alone
+if __name__ == "__main__":
+ print ('Testing the SplitWindowLST class')
+ print
+ test_split_window_lst()
More information about the grass-commit
mailing list