[GRASS-SVN] r73051 - grass-addons/grass7/imagery/i.sentinel.mask

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Aug 6 07:13:33 PDT 2018


Author: Robifag
Date: 2018-08-06 07:13:33 -0700 (Mon, 06 Aug 2018)
New Revision: 73051

Added:
   grass-addons/grass7/imagery/i.sentinel.mask/i_sentinel_mask_ES.png
Modified:
   grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.html
   grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.py
Log:
Improved the manual page and cleaned the code following PEP8 and GRASS Python Scripting rules

Modified: grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.html
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.html	2018-08-06 07:25:33 UTC (rev 73050)
+++ grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.html	2018-08-06 14:13:33 UTC (rev 73051)
@@ -1,9 +1,13 @@
 <h2>DESCRIPTION</h2>
-<em>i.sentinel.mask</em> allows to automatically identify clouds and their shadows in Sentinel 2 images. 
+<em>i.sentinel.mask</em> allows to automatically identify clouds and their
+shadows in Sentinel 2 images. 
 <p>
-The implemented procedure consists essentially of values thresholds, comparisons and calculations between bands and they lead to two different 
-rough maps of clouds and shadows which require further improvements and elaborations (e.g. transformation from raster to vector, cleaning geometries,
- removing small areas, checking topology, etc.) carried out in the different steps of the procedure. 
+The implemented procedure consists essentially of an algorithm based on values
+thresholds, comparisons and calculations between bands which leads to two
+different rough maps of clouds and shadows. These require further improvements
+and elaborations (e.g. transformation from raster to vector, cleaning
+geometries, removing small areas, checking topology, etc.) carried out in the
+different steps of the procedure. 
 
 <table cellspacing="2" cellpadding="2" width="100%" border="0">
 	<tbody>
@@ -16,18 +20,23 @@
 </table>
 
 <p>
-The algorithm has been developed starting from rules found in literature (Parmes et. al 2017) and conveniently refined.<br>
-Regarding the detection of shadows, the algorithm has been developed to identify only the shadows of clouds on the ground but some misclassification can occur. Often shadows and water have in fact similar reflectance 
-values which can lead to erroneous classification of water bodies as shadows. Therefore, in order to increase the accuracy of 
-the final shadow mask, a control check is implemented. Clouds and shadows are spatially intersected in order to remove misclassified areas. 
-This means that all those shadow geometries which do not intersect a cloud geometry are removed.
+The algorithm has been developed starting from rules found in literature
+(Parmes et. al 2017) and conveniently refined.<br> Regarding the detection of
+shadows, the algorithm has been developed to identify only the shadows of
+clouds on the ground. Obviously, some misclassifications can occur. Often
+shadows and water have in fact, similar reflectance values which can lead to
+erroneous classification of water bodies as shadows. Therefore, in order to
+increase the accuracy of the final shadow mask, a control check is implemented.
+Clouds and shadows are spatially intersected in order to remove misclassified
+areas. This means that all those shadow geometries which do not intersect a
+cloud geometry are removed.
 
 <div align="center" style="margin: 10px">
 <a href="i_sentinel_mask_CS.png">
 <img src="i_sentinel_mask_CS.png" width="30%">
-</a><br>
-<i>Fig: Module General WorkFlow</i>
-</div>
+</a><br><br>
+<i>Fig: "Cleaning" pprocedure of the shadow masko</i>
+</div><br><br>
 <!--center>
 <img src="i_sentinel_mask_CS.png" width="30%">
 <br>
@@ -35,11 +44,16 @@
 </center-->
 
 <p>
-The algorithm works on reflectance values (Bottom of Atmosphere Reflectance) therefore the atmospheric correction has to be applied to all input bands (see <a href="i.atcorr.html">i.atcorr</a>)
+The algorithm works on reflectance values (Bottom of Atmosphere Reflectance)
+therefore, the atmospheric correction has to be applied to all input bands (see
+<a href="i.sentinel.preproc.html">i.sentinel.preproc</a> or <a
+href="i.atcorr.html">i.atcorr</a>)
 
 <p>
-All necessary input bands (blue, green, red, nir, nir8a, swir11, swir12) must be imported in GRASS and specified one by one or using a text file.
-The text file has to be written following the syntax below: <em>variable=your_map</em>
+All necessary input bands (blue, green, red, nir, nir8a, swir11, swir12) must
+be imported in GRASS and specified one by one or using an input text file. The
+text file has to be written following the syntax below:
+<em>variable=your_map</em>
 
 <div class="code"><pre>
 blue=<em>your_blue_map</em>
@@ -51,39 +65,92 @@
 swir12=<em>your_swir12_map</em>
 </pre></div>
 
-Tha variables names (blue, green, red, nir, nir8a, swir11, swir12) have to be written precisely like in the example above (e.g. not Blue, nor BLUE but blue), 
+Tha variables names (blue, green, red, nir, nir8a, swir11, swir12) have to be
+written precisely like in the example above (e.g. not Blue, nor BLUE but blue),
 no spaces or special characters are permitted.
 
 <p>
-The final outputs are two different vector maps, one for clouds and one for shadows.
+The final outputs are two different vector maps, one for clouds and one for
+shadows.
 <p>
-The metadata file (MTD_TL.xml) is required only if both masks (cloud and shadow)
-are computed. The module retrieves from this file the sun azimuth and zenith necessary for the shadow mask cleaning phase 
+The metadata file (MTD_TL.xml or S2A_OPER_MTD_L1C_TL_MPS__*.xml) is required
+only if both masks (cloud and shadow) are computed. The module retrieves from
+this file the sun azimuth and zenith necessary for the shadow mask cleaning
+phase 
 <em>(see the schema above)</em>
 <p>
-If flag <b>-s</b> is given all selected bands are rescaled using the specified scale factor [<b>scale_fac</b>=<em>integer</em>]. By default the scale factor is set to 10000, 
-the QUANTIFICATION_VALUE from the metadata of Sentinel 2 images.
+If flag <b>-s</b> is given all selected bands are rescaled using the specified
+scale factor [<b>scale_fac</b>=<em>integer</em>]. By default the scale factor
+is set to 10000, the QUANTIFICATION_VALUE from the metadata of Sentinel 2
+images.
 <p>
-The module takes the current region settings into accout. To ignore the current region and set it from the whole image, the flag <b>-r</b> has to be given.
+The module takes the current region settings into accout. To ignore the current
+region and set it from the whole image, the flag <b>-r</b> has to be given.
 <p>
-The module allows to compute only the cloud mask or both cloud and shadow masks. If flag <b>-c</b> is given, only the cloud procedure will be performed. The computation 
-of cloud mask is mandatory for shadow mask creation. In fact cloud map is used during the cleaning phase of the shadow mask in order to remove misclassifications.
+The module allows to compute only the cloud mask or both cloud and shadow
+masks. If flag <b>-c</b> is given, only the cloud procedure will be performed.
+The computation of cloud mask is mandatory for shadow mask creation. In fact
+cloud map is used during the cleaning phase of the shadow mask in order to
+remove misclassifications.
 
 <h2>EXAMPLE</h2>
-<div class="code">
-i.sentinel.mask -r -s input_file=/home/input_bands.txt cloud_mask=cloud_sen2 shadow_mask=shadow_sen2 mtd_file=/home/MTD_TL.xml scale_fac=1000
+<p>
+The example illustrates how to run i.sentinel.mask for a Sentinel 2A image
+(S2A_MSIL1C_20180713T155901_N0206_R097_T17SPV_20180713T211059.SAFE) in the
+North Carolina location.<br> Obviously, the image has been imported and
+atmospheric correction has been performed before running i.sentinel.mask .
+<div class="code"><pre>
+i.sentinel.mask -r input_file=path/input_cloud_mask.txt cloud_mask=cloud shadow_mask=shadow cloud_threshold=25000 shadow_threshold=5000 mtd_file=path/MTD_TL.xml
+</pre></div>
+<p>
+The input text file:
+<div class="code"><pre>
+blue=T17SPV_20180315T160021_B02_cor
+green=T17SPV_20180315T160021_B03_cor
+red=T17SPV_20180315T160021_B04_cor
+swir11=T17SPV_20180315T160021_B11_cor
+nir=T17SPV_20180315T160021_B08_cor
+swir12=T17SPV_20180315T160021_B12_cor
+nir8a=T17SPV_20180315T160021_B8A_cor
+</pre></div>
+<p>
+<b>-r</b> to set the computational region to the maximum image extente
+<br>
+<div align="center" style="margin: 10px">
+<a href="i_sentinel_mask_ES.png">
+<img src="i_sentinel_mask_ES.png" width="1422" height="565" border="0">
+</a><br><br>
+<i>Figure1 (left): Sentinel-2A Band 02 - Figure2 (right): Sentinel-2A Band 02
+with computed cloud and shadow masks</i>
 </div>
-<p>
--r to set the computational region to the maximum image extente and -s to rescale the input bands with the specified scale factor (in this case 1000)
 
+<h2>IMPORTANT NOTES</h2>
+i.sentinel.mask works for Sentinel 2 images whose names follow both the New
+Compact Naming Convention (e.g.
+S2A_MSIL1C_20170527T102031_N0205_R065_T32TMQ_20170527T102301.SAFE) and the Old
+format Naming Convention (e.g.
+S2A_OPER_PRD_MSIL1C_PDMC_20160930T155112_R079_V20160930T095022_20160930T095944.SAFE).
+Therefore, both the MTD_TL.xml and S2A_OPER_MTD_L1C_TL_MPS__*.xml file can be
+provided as input for the computation of shadow mask. Both files can be found
+in the <em>GRANULE</em> folder of the downloaded *.SAFE product.<br>
+For further information about the naming convention see <a
+href="https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/naming-convention">ESA
+Sentinel User Guide</a>.
+
 <h2>REFERENCE</h2>
 <ul>
 <li>Parmes et. al 2017</li>
 </ul>
 
+<h2>FOLLOW UP</h2>
+<ul>
+<li> Implement other existing algorithm of clouds and shadows detection
+</ul>
+
 <h2>SEE ALSO</h2>
 
 <em>
+<a href="i.sentinel.preproc.html">i.sentinel.preproc</a>,
 <a href="i.sentinel.download.html">i.sentinel.download</a>,
 <a href="i.sentinel.import.html">i.sentinel.import</a>,
 <a href="r.import.html">r.import</a>,

Modified: grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.py
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.py	2018-08-06 07:25:33 UTC (rev 73050)
+++ grass-addons/grass7/imagery/i.sentinel.mask/i.sentinel.mask.py	2018-08-06 14:13:33 UTC (rev 73051)
@@ -178,11 +178,7 @@
 def main ():
 
 
-    #import bands atmospherically corrected using arcsi (scale factor 1000 instead of 10000)
-    #############################################
-    # INPUT
-    #############################################
-    #temporary map names
+    # Temporary map names
     global tmp, t, mapset
     tmp = {}
     mapset = gscript.gisenv()['MAPSET']
@@ -201,18 +197,20 @@
     tmp["cl_shift"] = "cl_shift_" + processid
     tmp["overlay"] = "overlay_" + processid
     
-    #check temporary map names are not existing maps
+    # Check temporary map names are not existing maps
     for key, value in tmp.items():
         if gscript.find_file(value,
             element = 'vector',
             mapset = mapset)['file']:
-                gscript.fatal(("Temporary vector map <{}> already exists.").format(value))
+            gscript.fatal(('Temporary vector map <{}> \
+                already exists.').format(value))
         if gscript.find_file(value,
             element = 'cell',
             mapset = mapset)['file']:
-                gscript.fatal(("Temporary raster map <{}> already exists.").format(value))
+            gscript.fatal(('Temporary raster map <{}> \
+                already exists.').format(value))
 
-    #input file
+    # Input file
     mtd_file = options['mtd_file']
     bands = {} 
     if options['input_file']=='':
@@ -227,8 +225,15 @@
         input_file = options['input_file']
         for line in file(input_file):
             a = line.split('=')
-            if len(a) != 2 or a[0] not in ['blue', 'green', 'red', 'nir', 'nir8a', 'swir11', 'swir12' ]:
-                gscript.fatal("Syntax error in the txt file. See the manual for further information about the right syntax.")
+            if len(a) != 2 or a[0] not in ['blue',
+                'green',
+                'red',
+                'nir',
+                'nir8a',
+                'swir11',
+                'swir12' ]:
+                gscript.fatal('Syntax error in the txt file. See the manual \
+                    for further information about the right syntax.')
             a[1] = a[1].strip()
             bands[a[0]] = a[1]
     d = 'double'
@@ -240,20 +245,29 @@
     cloud_mask = options['cloud_mask']
     shadow_mask = options['shadow_mask']
 
-    if bands['blue'] == '' or bands['green'] == '' or bands['red'] == '' or bands['nir'] == '' or bands['nir8a'] == '' or bands['swir11'] == '' or bands['swir12'] == '':
-        gscript.fatal(("All input bands (blue, green, red, nir, nir8a, swir11, swir12) are required"))
+    if (bands['blue'] == '' or
+        bands['green'] == '' or
+        bands['red'] == '' or
+        bands['nir'] == '' or
+        bands['nir8a'] == ''or
+        bands['swir11'] == '' or
+        bands['swir12'] == ''):
+        gscript.fatal('All input bands (blue, green, red, nir, nir8a, swir11, \
+            swir12) are required')
 
     for key, value in bands.items():
         if not gscript.find_file(value,
             element = 'cell',
             mapset = mapset)['file']:
-                gscript.fatal(("Raster map <{}> not found.").format(value))
+            gscript.fatal(('Raster map <{}> not found.').format(value))
 
     if not flags["c"]:
         if options['mtd_file']== '':
-            gscript.fatal("Metadata file is required for shadow mask computation. Please specified it")
+            gscript.fatal('Metadata file is required for shadow mask \
+                computation. Please specified it')
         if options['shadow_mask']=='':
-            gscript.fatal("Output name is required for shadow mask. Please specified it")
+            gscript.fatal('Output name is required for shadow mask. \
+                Please specified it')
 
     if flags["r"]:
         gscript.use_temp_region()
@@ -260,9 +274,11 @@
         gscript.run_command('g.region',
             rast=bands.values(),
             flags='a')
-        gscript.message(_("--- The computational region has been temporarily set to image max extent ---"))
+        gscript.message(_('--- The computational region has been temporarily \
+            set to image max extent ---'))
     else:
-        gscript.warning(_("All subsequent operations will be limited to the current computational region"))
+        gscript.warning(_('All subsequent operations will be limited to the \
+            current computational region'))
 
     if flags["s"]:
         gscript.message(_('--- Start rescaling bands ---'))
@@ -279,8 +295,9 @@
     else:
         gscript.warning(_('Any rescale factor has been applied'))
         for key, b in bands.items():
-            if gscript.raster_info(b)['datatype'] != "DCELL" and gscript.raster_info(b)['datatype'] != "FCELL":
-                gscript.fatal("Raster maps must be DCELL o FCELL")
+            if (gscript.raster_info(b)['datatype'] != "DCELL" and
+                gscript.raster_info(b)['datatype'] != "FCELL"):
+                gscript.fatal('Raster maps must be DCELL o FCELL')
             else:
                 f_bands = bands
 
@@ -289,10 +306,11 @@
         gscript.message(_(fb))
         stats = gscript.parse_command('r.univar', flags='g', map=fb)
         raster_max[key] = (float(stats['max']))
-    gscript.message(_('--- Computed maximum value: {} ---'.format(raster_max.values())))
+    gscript.message(_('--- Computed maximum value: {} ---'.format(
+        raster_max.values())))
     gscript.message(_('--- Statistics have been computed! ---'))
 
-    #### start of Clouds detection  (some rules from litterature)###
+    # Start of Clouds detection  (some rules from litterature)
     gscript.message(_('--- Start clouds detection procedure ---'))
     gscript.message(_('--- Computing cloud mask... ---'))
     first_rule = '(({} > (0.08*{})) && ({} > (0.08*{})) && ({} > (0.08*{})))'.format(
@@ -342,10 +360,10 @@
         tool='rmarea',
         threshold=cloud_threshold)
     gscript.message(_('--- Finish cloud detection procedure ---'))
-    ### end of Clouds detection ####
+    # End of Clouds detection
 
     if not flags["c"]:
-        ### start of shadows detection ###
+        # Start of shadows detection
         gscript.message(_('--- Start shadows detection procedure ---'))
         gscript.message(_('--- Computing shadow mask... ---'))
         sixth_rule = '((({} > {}) && ({} < {}) && ({} < 0.1) && ({} < 0.1)) \
@@ -387,14 +405,13 @@
             tool='rmarea',
             threshold=shadow_threshold)
         gscript.message(_('--- Finish Shadows detection procedure ---'))
-        ### end of shadows detection ###
+        # End of shadows detection
 
-        #####################################################################
         # START shadows cleaning Procedure (remove shadows misclassification)
-        #####################################################################
-        ### start shadow mask preparation ###
+        # Start shadow mask preparation
 
-        gscript.message(_('--- Start removing misclassification from the shadow mask ---'))
+        gscript.message(_('--- Start removing misclassification from \
+            the shadow mask ---'))
         gscript.message(_('--- Data preparation... ---'))
         gscript.run_command('v.centroids',
             input=tmp["shadow_temp_mask"],
@@ -436,8 +453,8 @@
             map=tmp["addcat"],
             columns='value')
 
-        ### end shadow mask preparation ###
-        ### start cloud mask preparation ###
+        # End shadow mask preparation
+        # Start cloud mask preparation
 
         gscript.run_command('v.db.droptable',
             map=cloud_mask,
@@ -446,10 +463,12 @@
             map=cloud_mask,
             columns='value')
 
-        ### end cloud mask preparation ###
-        ### shift cloud mask using dE e dN ###
-        ### start reading mean sun zenith and azimuth from xml file to compute dE and dN automatically ###
-        gscript.message(_('--- Reading mean sun zenith and azimuth from MTD_TL.xml file to compute clouds shift ---'))
+        # End cloud mask preparation
+        # Shift cloud mask using dE e dN
+        # Start reading mean sun zenith and azimuth from xml file to compute 
+        #dE and dN automatically
+        gscript.message(_('--- Reading mean sun zenith and azimuth from \
+            MTD_TL.xml file to compute clouds shift ---'))
         xml_tree = et.parse(mtd_file)
         root = xml_tree.getroot()
         ZA = []
@@ -459,12 +478,17 @@
                 ZA.append (subelem.text)
         z = float(ZA[0])
         a = float(ZA[1])
-        gscript.message(_('--- the mean sun Zenith is: {:.3f} deg ---'.format(z)))
-        gscript.message(_('--- the mean sun Azimuth is: {:.3f} deg ---'.format(a)))
+        gscript.message(_('--- the mean sun Zenith is: \
+            {:.3f} deg ---'.format(z)))
+        gscript.message(_('--- the mean sun Azimuth is: \
+            {:.3f} deg ---'.format(a)))
 
-        ### stop reading mean sun zenith and azimuth from xml file to compute dE and dN automatically ###
-        ### start computing the east and north shift for clouds and the overlapping area between clouds and shadows at steps of 100m ###
-        gscript.message(_('--- Start computing the east and north clouds shift at steps of 100m of clouds height---'))
+        # Stop reading mean sun zenith and azimuth from xml file to compute dE 
+        #and dN automatically
+        # Start computing the east and north shift for clouds and the 
+        #overlapping area between clouds and shadows at steps of 100m
+        gscript.message(_('--- Start computing the east and north clouds \
+            shift at steps of 100m of clouds height---'))
         H = 1000
         dH = 100
         HH = []
@@ -511,10 +535,11 @@
             area2 = gscript.parse_key_val(area, sep='|')
             AA.append(float(area2['total area']))
 
-        # find the maximum overlapping area between clouds and shadows#
+        # Find the maximum overlapping area between clouds and shadows
         index_maxAA = numpy.argmax(AA)
 
-        # clouds are shifted using the clouds height corresponding to the maximum overlapping area then are intersect with shadows#
+        # Clouds are shifted using the clouds height corresponding to the
+        #maximum overlapping area then are intersect with shadows
         gscript.run_command('v.transform',
             input=cloud_mask,
             output=tmp["cl_shift"],
@@ -531,9 +556,12 @@
             operator='intersects',
             quiet=True)
 
-        gscript.message(_('--- the estimated clouds height is: {} m ---'.format(HH[index_maxAA])))
-        gscript.message(_('--- the estimated east shift is: {:.2f} m ---'.format(dE[index_maxAA])))
-        gscript.message(_('--- the estimated north shift is: {:.2f} m ---'.format(dN[index_maxAA])))
+        gscript.message(_('--- the estimated clouds height is: \
+            {} m ---'.format(HH[index_maxAA])))
+        gscript.message(_('--- the estimated east shift is: \
+            {:.2f} m ---'.format(dE[index_maxAA])))
+        gscript.message(_('--- the estimated north shift is: \
+            {:.2f} m ---'.format(dN[index_maxAA])))
     else:
         gscript.warning(_('No shadow mask will be computed'))
 
@@ -540,18 +568,23 @@
 def cleanup():
     if flags["r"]:
         gscript.del_temp_region()
-    gscript.message(_('--- The computational region has been reset to the previous one---'))
+    gscript.message(_('--- The computational region has been reset to \
+        the previous one---'))
     if flags["t"]:
         gscript.message(_('--- No temporary files have been deleted ---'))
     else:
         for key, value in tmp.items():
-            if gscript.find_file(value, element = 'vector', mapset = mapset)['file']:
+            if gscript.find_file(value,
+                element = 'vector',
+                mapset = mapset)['file']:
                 gscript.run_command("g.remove",
                     flags="f",
                     type='vector',
                     name=",".join([tmp[m] for m in tmp.keys()]),
                     quiet=True)
-            if gscript.find_file(value, element = 'cell', mapset = mapset)['file']:
+            if gscript.find_file(value,
+                element = 'cell',
+                mapset = mapset)['file']:
                 gscript.run_command("g.remove",
                     flags="f",
                     type='raster',

Added: grass-addons/grass7/imagery/i.sentinel.mask/i_sentinel_mask_ES.png
===================================================================
(Binary files differ)

Index: grass-addons/grass7/imagery/i.sentinel.mask/i_sentinel_mask_ES.png
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.mask/i_sentinel_mask_ES.png	2018-08-06 07:25:33 UTC (rev 73050)
+++ grass-addons/grass7/imagery/i.sentinel.mask/i_sentinel_mask_ES.png	2018-08-06 14:13:33 UTC (rev 73051)

Property changes on: grass-addons/grass7/imagery/i.sentinel.mask/i_sentinel_mask_ES.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property


More information about the grass-commit mailing list