Merge pull request #1 from mundialis/multi_module

i.sentinel_1 addon for Sentinel-1 SAR data processing

.flake8

@@ -0,0 +1,12 @@
+[flake8]
+
+# E501 line too long (83 > 79 characters)
+# F821 undefined name '_'
+
+exclude = .git
+max-line-length = 88
+per-file-ignores =
+    ./i.sentinel_1.change/i.sentinel_1.change.py: F821, E501
+    ./i.sentinel_1.download_asf/i.sentinel_1.download_asf.py: F821, E501
+    ./i.sentinel_1.import/i.sentinel_1.import.py: F821, E501
+    ./i.sentinel_1.mosaic/i.sentinel_1.mosaic.py: F821, E501

.github/workflows/flake8.yml

@@ -0,0 +1,25 @@
+name: Python Flake8 Code Quality
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  flake8:
+    runs-on: ubuntu-20.04
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Set up Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: 3.8
+
+      - name: Install
+        run: |
+          python -m pip install --upgrade pip
+          pip install flake8==3.8.4
+      - name: Run Flake8
+        run: |
+          flake8 --config=.flake8 --count --statistics --show-source --jobs=$(nproc) .

Makefile

@@ -0,0 +1,17 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.sentinel_1
+
+# note: to deactivate a module, just place a file "DEPRECATED" into the subdir
+ALL_SUBDIRS := ${sort ${dir ${wildcard */.}}}
+DEPRECATED_SUBDIRS := ${sort ${dir ${wildcard */DEPRECATED}}}
+RM_SUBDIRS := bin/ docs/ etc/ scripts/
+SUBDIRS_1 := $(filter-out $(DEPRECATED_SUBDIRS), $(ALL_SUBDIRS))
+SUBDIRS := $(filter-out $(RM_SUBDIRS), $(SUBDIRS_1))
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: parsubdirs htmldir
+
+install: installsubdirs
+	$(INSTALL_DATA) $(PGM).html $(INST_DIR)/docs/html/

i.sentinel_1.change/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.sentinel_1.change
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

i.sentinel_1.change/i.sentinel_1.change.html

@@ -0,0 +1,32 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.sentinel_1.change</em> is a GRASS GIS addon that extracts areas of changed backscatter values from
+two Sentinel-1 input scenes. The input backscatter values are expected in dB scale.
+The addon calculates two polarisation ratios (date2_vv/date1_vv and date2_vh/date1_vh),
+smoothes them with a median filter in  <a href="https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html">r.neighbors</a>
+and applies the user-defined thresholds 1+<b>change_threshold</b> and 1-<b>change_threshold</b> to identify
+areas of change. For example, a <b>change_threshold</b> of 0.5 would mean that areas with a ratio <= 0.5 are
+ identified as signal decrease, and areas with a ratio of >= 1.5 are identified as signal increase.
+
+<h2>EXAMPLE</h2>
+<div class="code"><pre>
+i.sentinel_1.change date1_vv=Date1_Sigma0_VV_log date1_vh=Date1_Sigma0_VH_log date2_vv=Date2_Sigma0_VV_log date2_vh=Date2_Sigma0_VH_log output=changed_areas change_threshold=0.75 min_size=1.0 --o
+</pre></div>
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.mapcalc.html">r.mapcalc</a>,
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html">r.neighbors</a>,
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.null.html">r.null</a>,
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.colors.html">r.colors</a>,
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.category.html">r.category</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Guido Riembauer, <a href="https://www.mundialis.de/">mundialis</a>
+
+<!--
+<p>
+<i>Last changed: $Date$</i>
+-->

i.sentinel_1.change/i.sentinel_1.change.py

@@ -0,0 +1,199 @@
+#!/usr/bin/env python3
+#
+############################################################################
+#
+# MODULE:      i.sentinel_1.change
+# AUTHOR(S):   Guido Riembauer, <riembauer at mundialis.de>
+#
+# PURPOSE:     extracts changes from Sentinel-1 Scenes from two dates
+# COPYRIGHT:   (C) 2021-2022 mundialis GmbH & Co. KG, and the GRASS Development Team
+#
+#  This program is free software; you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation; either version 2 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#############################################################################
+# %Module
+# % description: Extracts changes from Sentinel-1 Scenes from two dates.
+# % keyword: imagery
+# % keyword: Sentinel
+# % keyword: satellite
+# % keyword: change
+# %End
+
+# %option G_OPT_R_INPUT
+# % key: date1_vv
+# % label: Input raster map of first date in VV polarization (dB)
+# %end
+
+# %option G_OPT_R_INPUT
+# % key: date1_vh
+# % label: Input raster map of first date in VH polarization (dB)
+# %end
+
+# %option G_OPT_R_INPUT
+# % key: date2_vv
+# % label: Input raster map of second date in VV polarization (dB)
+# %end
+
+# %option G_OPT_R_INPUT
+# % key: date2_vh
+# % label: Input raster map of second date in VH polarization (dB)
+# %end
+
+# %option
+# % key: change_threshold
+# % type: double
+# % label: Minimum change in ratio (date2/date1) map to be considered (all areas with values from (1-<change_threshold>) to (1+<change_threshold>) will be considered unchanged)
+# % options: 0.0-1.0
+# % answer: 0.75
+# %end
+
+# %option
+# % key: min_size
+# % type: double
+# % label: Minimum size of identified changed areas in ha
+# % answer: 1.0
+# %end
+
+# %option G_OPT_R_OUTPUT
+# % key: output
+# % label: Output raster map containing areas of change
+# %end
+
+
+import grass.script as grass
+import os
+import atexit
+
+rm_rasters = []
+
+
+def cleanup():
+    nuldev = open(os.devnull, "w")
+    kwargs = {"flags": "f", "quiet": True, "stderr": nuldev}
+    for rmrast in rm_rasters:
+        if grass.find_file(name=rmrast, element="raster")["file"]:
+            grass.run_command("g.remove", type="raster", name=rmrast, **kwargs)
+
+
+def main():
+
+    global rm_rasters
+    date1_vv = options["date1_vv"]
+    date1_vh = options["date1_vh"]
+    date2_vv = options["date2_vv"]
+    date2_vh = options["date2_vh"]
+    minsize_ha = options["min_size"]
+    change_thresh = float(options["change_threshold"])
+    output = options["output"]
+
+    pid_str = str(os.getpid())
+    tmp_ratio_vv = "ratio_vv_{}".format(pid_str)
+    tmp_ratio_vh = "ratio_vh_{}".format(pid_str)
+    rm_rasters.extend([tmp_ratio_vv, tmp_ratio_vh])
+    # calculate ratios of natural values (10^(<db_val/10>)) to get positive
+    # numbers only
+    exp_ratio_vv = "{} = float(10^({}/10))/float(10^({}/10))".format(
+        tmp_ratio_vv, date2_vv, date1_vv
+    )
+    exp_ratio_vh = "{} = float(10^({}/10))/float(10^({}/10))".format(
+        tmp_ratio_vh, date2_vh, date1_vh
+    )
+    for exp in [exp_ratio_vv, exp_ratio_vh]:
+        grass.message(_("Calculating ratio raster {}...").format(exp.split("=")[0]))
+        grass.run_command("r.mapcalc", expression=exp, quiet=True)
+
+    result_rasts = []
+    for idx, tuple in enumerate([(tmp_ratio_vv, "vv"), (tmp_ratio_vh, "vh")]):
+        ratio = tuple[0]
+        grass.message(_("Smoothing ratio raster {}...").format(ratio))
+
+        ratio_smoothed = "{}_tmp_{}_smoothed".format(ratio, idx)
+        rm_rasters.append(ratio_smoothed)
+        grass.run_command(
+            "r.neighbors",
+            input=ratio,
+            method="median",
+            size=9,
+            output=ratio_smoothed,
+            quiet=True,
+        )
+
+        # apply thresholds
+        upper_thresh_orig = 1.0 + change_thresh
+        lower_thresh_orig = 1.0 - change_thresh
+        change_raster = f"{ratio_smoothed}_changes"
+        rm_rasters.append(change_raster)
+        result_rasts.append(change_raster)
+        ch_exp = (
+            f"{change_raster} = if({ratio}>={upper_thresh_orig},1,"
+            f"if({ratio}<={lower_thresh_orig},2,null()))"
+        )
+        grass.run_command("r.mapcalc", expression=ch_exp, quiet=True)
+
+    patched = f"changes_patched_{pid_str}"
+    rm_rasters.append(patched)
+    grass.run_command("r.patch", input=result_rasts, output=patched, quiet=True)
+
+    no_contradictions = f"no_contradictions_{pid_str}"
+    rm_rasters.append(no_contradictions)
+    contr_exp = (
+        f"{no_contradictions} = if(({result_rasts[0]}+"
+        f"{result_rasts[1]})==3,null(),{patched})"
+    )
+    grass.run_command("r.mapcalc", expression=contr_exp, quiet=True)
+    try:
+        grass.run_command(
+            "r.reclass.area",
+            input=no_contradictions,
+            output=output,
+            mode="greater",
+            value=minsize_ha,
+            quiet=True,
+        )
+    except Exception as e:
+        # reclass.area fails if there are no areas larger than minsize_ha
+        grass.warning(
+            _(
+                f"An Exception occured in r.reclass.area: {e}\n"
+                f"No areas larger than {minsize_ha} ha found. "
+                "Producing an output raster without changes..."
+            )
+        )
+        grass.run_command("r.mapcalc", expression=f"{output} = null()", quiet=True)
+
+    # adapt the colors
+    grass.run_command("r.null", null=0, map=output, quiet=True)
+    colors_new = ["0 255:255:255", "1 0:200:0", "2 200:0:0"]
+    color_str = "\n".join(colors_new)
+    col_proc = grass.feed_command("r.colors", map=output, rules="-", quiet=True)
+    col_proc.stdin.write(color_str.encode())
+    col_proc.stdin.close()
+    col_proc.wait()
+
+    # adapt the category labels
+    labels = ["0:No signficant Change", "1:Signal Increase", "2:Signal Decrease"]
+    category_text = "\n".join(labels)
+
+    # assign labels
+    cat_proc = grass.feed_command("r.category", map=output, rules="-", separator=":")
+    cat_proc.stdin.write(category_text.encode())
+    cat_proc.stdin.close()
+    # feed_command does not wait until finished
+    cat_proc.wait()
+
+    grass.message(_(f"Successfully created output raster map <{output}>"))
+    return 0
+
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    atexit.register(cleanup)
+    main()