Merge pull request #24 from griembauer/clouds_strds

t.sentinel.import: specify cloud STDS import

i.sentinel.import.worker/i.sentinel.import.worker.py

@@ -89,6 +89,16 @@
 # % required: no
 # %end
 
+# %option
+# % key: cloud_output
+# % type: string
+# % required: no
+# % multiple: no
+# % options: vector,raster
+# % answer: vector
+# % label: Create cloud mask as raster or vector map
+# %end
+
 # %flag
 # % key: r
 # % description: Reproject raster data using r.import if needed
@@ -248,6 +258,8 @@ def main():
         kwargs['metadata'] = options['metadata']
     if options["pattern_file"]:
         kwargs["pattern_file"] = options["pattern_file"]
+    if options["cloud_output"]:
+        kwargs["cloud_output"] = options["cloud_output"]
 
     kwargsstr = ""
     flagstr = ""

t.sentinel.import/t.sentinel.import.html

@@ -2,7 +2,15 @@ <h2>DESCRIPTION</h2>
 
 <em>t.sentinel.import</em> is a GRASS GIS addon Python script to
 download and import the Sentinel-2 scenes and create a STRDS.
-
+<p>
+Since Sentinel-2 <a href="https://forum.step.esa.int/t/info-introduction-of-additional-radiometric-offset-in-pb04-00-products/35431">Processing Baseline 4.0.0</a>,
+a systematic offset has been introduced to all reflectance values.<br>
+To account for this, the <b>offset</b> option allows the user to
+indicate how reflectance data should be corrected (e.g. -1000).
+<p>
+When using the <em>-c</em> flag to import cloud masks, the user may define the type and name of the <br>
+  cloud STDS to be imported by the <em>stvds_clouds</em>/<em>strds_clouds</em> option. If only the
+<em>-c</em> flag is given, a cloud STVDS is imported and named <em>[strds_output]_clouds</em>
 <h2>EXAMPLE</h2>
 
 <div class="code"><pre>

t.sentinel.import/t.sentinel.import.py

@@ -147,6 +147,22 @@
 # % gisprompt: new,stds,strds
 # %end
 
+# %option
+# % key: stvds_clouds
+# % type: string
+# % required: no
+# % multiple: no
+# % description: Name of the output cloudmask space time vector dataset. If not provided, it will be <strds_output>_clouds
+# %end
+
+# %option
+# % key: strds_clouds
+# % type: string
+# % required: no
+# % multiple: no
+# % description: Name of the output cloudmask space time raster dataset. If not provided, it will be <strds_output>_clouds
+# %end
+
 # %option
 # % key: directory
 # % type: string
@@ -186,13 +202,23 @@
 # % answer: 1
 # %end
 
+# %option
+# % key: offset
+# % type: integer
+# % required: no
+# % description: Offset to add to the Sentinel-2 bands to due to specific processing baseline (e.g. -1000)
+# %end
+
 # %rules
 # % collective: start, end, producttype
 # % excludes: s2names, start, end, producttype
 # % excludes: input_dir, s2names, start, end, producttype, settings, clouds
 # % required: input_dir, start, s2names
 # % requires: -a, sen2cor_path
 # % requires: -e, s2names
+# % requires: stvds_clouds, -c
+# % requires: strds_clouds, -c
+# % exclusive: stvds_clouds, strds_clouds
 # %end
 
 
@@ -201,6 +227,7 @@
 import multiprocessing as mp
 import os
 import re
+import shutil
 import sys
 
 import grass.script as grass
@@ -322,6 +349,11 @@ def main():
                     "\n" +
                     "g.extension i.zero2null")
 
+    if not grass.find_program('r.mapcalc.tiled', '--help'):
+        grass.fatal(_("The 'r.mapcalc.tiled' module was not found, install it first:") +
+                    "\n" +
+                    "g.extension r.mapcalc.tiled")
+
     # create temporary directory to download data
     if tmpdirectory:
         if not os.path.isdir(tmpdirectory):
@@ -457,6 +489,10 @@ def main():
             }
             if options["extent"] == "region":
                 import_kwargs["region"] = currentregion
+            if flags["c"]:
+                import_kwargs["cloud_output"] = "vector"
+                if options["strds_clouds"]:
+                    import_kwargs["cloud_output"] = "raster"
             if single_folders is True:
                 directory = os.path.join(download_dir, subfolder)
             else:
@@ -496,8 +532,42 @@ def main():
             cloudlist.append(vect)
             grass.run_command('g.copy', vector=vect + '@' + new_mapset + ',' + vect)
         for rast in grass.parse_command('g.list', type='raster', mapset=new_mapset):
-            maplist.append(rast)
-            grass.run_command('g.copy', raster=rast + '@' + new_mapset + ',' + rast)
+            grass.run_command('g.copy',
+                              raster=rast + '@' + new_mapset + ',' + rast)
+            if "CLOUDS" in rast:
+                cloudlist.append(rast)
+            else:
+                maplist.append(rast)
+
+                if options["offset"]:
+                    # save the description.json
+                    tmp_desc_dir = os.path.join(tmpdirectory, "descriptions_json")
+                    if not os.path.isdir(tmp_desc_dir):
+                        try:
+                            os.makedirs(tmp_desc_dir)
+                        except:
+                            grass.fatal(_(f"Unable to create directory {tmp_desc_dir}"))
+
+                    desc_file_save = os.path.join(tmp_desc_dir, f"{rast}_description.json")
+                    desc_file_in = os.path.join(json_standard_folder, rast, "description.json")
+                    shutil.copy(desc_file_in, desc_file_save)
+
+                    # calculate offset (metadata in cell_misc will be lost)
+                    tmp_rast = f"rast_tmp_{os.getpid()}"
+                    # clipping to 0 to keep the value within the valid 0-10000 range
+                    mapc_exp = (f"{tmp_rast} = if({rast} + {options['offset']} < 0, "
+                                f"0, {rast} + {options['offset']} )")
+                    grass.run_command(f"r.mapcalc.tiled",
+                                      expression=mapc_exp,
+                                      nprocs=nprocs_final,
+                                      quiet=True)
+                    grass.run_command("g.copy",
+                                      raster=f"{tmp_rast},{rast}",
+                                      overwrite=True,
+                                      quiet=True)
+                    # copy the description.json back
+                    shutil.copy(desc_file_save, desc_file_in)
+
         grass.utils.try_rmdir(os.path.join(gisdbase, location, new_mapset))
     # space time dataset
     grass.message(_("Creating STRDS of Sentinel scenes ..."))
@@ -534,11 +604,20 @@ def main():
         # remove registerfile
         grass.try_remove(registerfile)
 
-        if flags['c']:
-            stvdsclouds = strds + '_clouds'
+        if flags["c"]:
+            dtype="vector"
+            stds_type = "stvds"
+            if options["strds_clouds"]:
+                stdsclouds = options["strds_clouds"]
+                dtype="raster"
+                stds_type = "strds"
+            elif options["stvds_clouds"]:
+                stdsclouds = options["stvds_clouds"]
+            else:
+                stdsclouds = strds + '_clouds'
             grass.run_command(
-                't.create', output=stvdsclouds, title="Sentinel-2 clouds",
-                desc="Sentinel-2 clouds", quiet=True, type='stvds')
+                't.create', output=stdsclouds, title="Sentinel-2_sen2cor_clouds",
+                desc="Sentinel-2_sen2cor_clouds", quiet=True, type=stds_type)
             registerfileclouds = grass.tempfile()
             fileclouds = open(registerfileclouds, 'w')
             for imp_clouds in cloudlist:
@@ -549,8 +628,8 @@ def main():
                 fileclouds.write("%s|%s %s\n" % (imp_clouds, date_str2, clock_str2))
             fileclouds.close()
             grass.run_command(
-                't.register', type='vector', input=stvdsclouds, file=registerfileclouds, quiet=True)
-            grass.message("<%s> is created" % (stvdsclouds))
+                't.register', type=dtype, input=stdsclouds, file=registerfileclouds, quiet=True)
+            grass.message("<%s> is created" % (stdsclouds))
             # remove registerfile
             grass.try_remove(registerfileclouds)