[GRASS-SVN] r68397 - grass-addons/grass7/imagery/i.landsat8.swlst

svn_grass at osgeo.org svn_grass at osgeo.org
Sat May 7 14:32:41 PDT 2016


Author: nikosa
Date: 2016-05-07 14:32:41 -0700 (Sat, 07 May 2016)
New Revision: 68397

Modified:
   grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html
   grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py
Log:
Merged https://github.com/NikosAlexandris/i.landsat8.swlst/pull/1

Modified: grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html	2016-05-07 18:52:46 UTC (rev 68396)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.html	2016-05-07 21:32:41 UTC (rev 68397)
@@ -1,5 +1,5 @@
 <h2 id="description">DESCRIPTION</h2>
-<p><em>i.landsat8.swlst</em> is an implementation of a robust and practical Split-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>
+<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>
@@ -252,7 +252,7 @@
 </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>[0] Du, Chen; Ren, Huazhong; Qin, Qiming; Meng, Jinjie; Zhao, Shaohua. 2015. "A Practical Split-Window Algorithm for Estimating Land Surface Temperature from Landsat 8 Data." Remote Sens. 7, no. 1: 647-665. http://www.mdpi.com/2072-4292/7/1/647/htm</p></li>
 <li><p>[1] Huazhong Ren, Chen Du, Qiming Qin, Rongyuan Liu, Jinjie Meng, and Jing Li. "Atmospheric Water Vapor Retrieval from Landsat 8 and Its Validation." 3045--3048. IEEE, 2014.</p></li>
 <li><p>[2] Ren, H., Du, C., Liu, R., Qin, Q., Yan, G., Li, Z. L., & Meng, J. (2015). Atmospheric water vapor retrieval from Landsat 8 thermal infrared images. Journal of Geophysical Research: Atmospheres, 120(5), 1723-1738.</p></li>
 </ul>

Modified: grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py
===================================================================
--- grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py	2016-05-07 18:52:46 UTC (rev 68396)
+++ grass-addons/grass7/imagery/i.landsat8.swlst/i.landsat8.swlst.py	2016-05-07 21:32:41 UTC (rev 68397)
@@ -350,6 +350,9 @@
     """
     grass.run_command('g.remove', flags='f', type="rast",
                       pattern='tmp.{pid}*'.format(pid=os.getpid()), quiet=True)
+    
+    if grass.find_file(name='MASK', element='cell')['file']:
+        r.mask(flags='r', verbose=True)
 
 
 def tmp_map_name(name):
@@ -552,8 +555,9 @@
     # 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)
+        run('g.rename', raster=(tmp_brightness_temperature, bt_output))
+        tmp_brightness_temperature = bt_output
+        del(bt_output)
 
     return tmp_brightness_temperature
 
@@ -575,21 +579,21 @@
            '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)
+    #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)
+    #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)
+    r.mask(raster=qa_band, maskcats=qa_pixel, flags='i', overwrite=True)
 
     # save for debuging
     #save_map(tmp_cloudmask)
 
-    del(qabits_expression)
-    del(cloud_masking_equation)
+    #del(qabits_expression)
+    #del(cloud_masking_equation)
 
 
 def replace_dummies(string, *args, **kwargs):
@@ -714,7 +718,7 @@
     
     # save land surface emissivity map?
     if emissivity_output:
-        run('g.copy', raster=(outname, emissivity_output))
+        run('g.rename', raster=(outname, emissivity_output))
 
 
 def determine_delta_emissivity(outname, landcover_map, delta_lse_expression):
@@ -746,7 +750,7 @@
 
     # save delta land surface emissivity map?
     if delta_emissivity_output:
-        run('g.copy', raster=(outname, delta_emissivity_output))
+        run('g.rename', raster=(outname, delta_emissivity_output))
 
 
 def get_cwv_window_means(outname, t1x, t1x_mean_expression):
@@ -835,7 +839,7 @@
 
     # save Column Water Vapor map?
     if cwv_output:
-        run('g.copy', raster=(outname, cwv_output))
+        run('g.rename', raster=(outname, cwv_output))
 
     # save for debuging
     #save_map(outname)
@@ -884,7 +888,7 @@
             source1=source1_cwv, source2=source2_cwv,
             history=history_cwv)
 
-        run('g.copy', raster=(outname, cwv_output))
+        run('g.rename', raster=(outname, cwv_output))
 
     del(cwv_expression)
     del(cwv_equation)
@@ -933,9 +937,14 @@
                                                   out_ti=t10,
                                                   in_tj=DUMMY_MAPCALC_STRING_T11,
                                                   out_tj=t11)
-
-    split_window_equation = equation.format(result=outname,
+    # Convert to Celsius?
+    if celsius:
+        split_window_expression = '({swe}) - 273.15'.format(swe=split_window_expression)
+        split_window_equation = equation.format(result=outname,
                                             expression=split_window_expression)
+    else:
+        split_window_equation = equation.format(result=outname,
+                                            expression=split_window_expression)
 
     grass.mapcalc(split_window_equation, overwrite=True)
 
@@ -945,26 +954,6 @@
     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
@@ -982,7 +971,7 @@
     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')
+    #tmp_lst = tmp_map_name('lst')
 
     # basic equation for mapcalc
     global equation, citation_lst
@@ -1023,7 +1012,10 @@
     
     # save Brightness Temperature maps?
     global brightness_temperature_prefix
-    brightness_temperature_prefix = options['prefix_bt']
+    if options['prefix_bt']:
+        brightness_temperature_prefix = options['prefix_bt']
+    else:
+        brightness_temperature_prefix = None
 
     global cwv_output
     cwv_window_size = int(options['window'])
@@ -1050,9 +1042,11 @@
     global info, null
     info = flags['i']
     keep_region = flags['k']
-    celsius = flags['c']
     timestamping = flags['t']
     null = flags['n']
+    
+    global celsius
+    celsius = flags['c']
 
     # ToDo:
     # shell = flags['g']
@@ -1070,11 +1064,11 @@
         # Improve below!
 
         if b10:
-            run('g.region', rast=b10)
+            run('g.region', rast=b10, align=b10)
             msg = msg.format(name=b10)
 
         elif t10:
-            run('g.region', rast=t10)
+            run('g.region', rast=t10, align=t10)
             msg = msg.format(name=t10)
 
         g.message(msg)
@@ -1155,12 +1149,17 @@
         if not average_emissivity_map:
             determine_average_emissivity(tmp_avg_lse, landcover_map,
                                          split_window_lst.average_lse_mapcalc)
+            if options['emissivity_out']:
+                tmp_avg_lse = options['emissivity_out']
+
         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)
+            if options['delta_emissivity_out']:
+                tmp_delta_lse = options['delta_emissivity_out']
 
     #
     # 4. Modified Split-Window Variance-Covariance Matrix > Column Water Vapor
@@ -1175,6 +1174,8 @@
     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())
+    if cwv_output:
+        tmp_cwv = cwv_output
 
     #
     # 5. Estimate Land Surface Temperature
@@ -1184,7 +1185,7 @@
         msg = '\n|* Will pick a random emissivity class!'
         grass.verbose(msg)
 
-    estimate_lst(tmp_lst, t10, t11,
+    estimate_lst(lst_output, t10, t11,
                  tmp_avg_lse, tmp_delta_lse, tmp_cwv,
                  split_window_lst.sw_lst_mapcalc)
 
@@ -1197,18 +1198,17 @@
 
     # time-stamping
     if timestamping:
-        add_timestamp(mtl_file, tmp_lst)
+        add_timestamp(mtl_file, lst_output)
 
         if cwv_output:
             add_timestamp(mtl_file, cwv_output)
 
-    # convert to celsius and apply color table?
+    # Apply color table
     if celsius:
-        kelvin_to_celsius(tmp_lst, tmp_lst)
-
+        run('r.colors', map=lst_output, color='celsius')
     else:
         # color table for kelvin
-        run('r.colors', map=tmp_lst, color='kelvin')
+        run('r.colors', map=lst_output, color='kelvin')
 
     # ToDo: helper function for r.support
     # strings for metadata
@@ -1231,13 +1231,13 @@
     source2_lst = landsat8_metadata.origin
 
     # history entry
-    run("r.support", map=tmp_lst, title=title_lst,
+    run("r.support", map=lst_output, 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))
+    #run("g.rename", rast=(tmp_lst, lst_output))
 
     # restore region
     if not keep_region:



More information about the grass-commit mailing list