t.rast.patch (#51)

This new addon module patches raster maps that have gaps in time with subsequent maps (within a space time raster dataset).

grass7/temporal/t.rast.patch/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast.patch
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)

grass7/temporal/t.rast.patch/t.rast.patch.html

@@ -0,0 +1,45 @@
+<h2>DESCRIPTION</h2>
+This module patches raster maps that have gaps in time with subsequent maps (within a space time raster dataset)
+using <em>r.patch</em>. Hence it is a wrapper for <em>r.patch</em> in the temporal domain.
+The input of this module is a single space time raster dataset, the
+output is a single raster map layer. A subset of the input space time
+raster dataset can be selected using the <b>where</b> option. The
+sorting of the raster map layer can be set using the <b>sort</b>
+option. Be aware that the sorting of the maps significantly influences
+the result of the patch. By default the maps are
+sorted by <b>desc</b> by the <i>start_time</i> so that the newest raster map
+is the first input map in <b>r.patch</b>.
+<p>
+<em>t.rast.patch</em> is a simple wrapper for the raster module
+<b>r.patch</b>.
+
+<h2>EXAMPLE</h2>
+The example uses the North Carolina extra time series of MODIS Land Surface Temperature
+maps (<a href="https://grass.osgeo.org/download/sample-data/">download</a>).
+(The mapset has to be unzip in one of the North Carolina locations.)
+<p>
+Patching the MODIS Land Surface Temperature for 2016 (filling missing pixels by subsequent maps in the time series):
+<div class="code"><pre>
+t.rast.patch input=LST_Day_monthly@modis_lst output=LST_Day_patched_2016 \
+  where="start_time >= '2016-01' and start_time <= '2016-12'"
+r.info LST_Day_patched_2016
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.patch.html">r.patch</a>,
+<a href="t.rast.series.html">t.rast.series</a>,
+<a href="t.create.html">t.create</a>,
+<a href="t.info.html">t.info</a>
+</em>
+<p>
+<a href="http://grasswiki.osgeo.org/wiki/Temporal_data_processing">Temporal data processing Wiki</a>
+
+<h2>AUTHOR</h2>
+
+Anika Bettge, mundialis GmbH &amp; Co. KG
+
+<!--
+<p><i>Last changed: $Date$</i>
+-->

grass7/temporal/t.rast.patch/t.rast.patch.py

@@ -0,0 +1,159 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:    t.rast.patch
+# AUTHOR(S):    Anika Bettge
+#
+# PURPOSE:    Patches rasters that have gaps with subsequent maps in time within a space time raster dataset using r.patch
+# COPYRIGHT:    (C) 2019 by by mundialis 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: Patches multiple space time raster maps into a single raster map using r.patch.
+#% keyword: temporal
+#% keyword: aggregation
+#% keyword: series
+#% keyword: raster
+#% keyword: merge
+#% keyword: patching
+#% keyword: time
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#%end
+
+#%option G_OPT_T_WHERE
+#%end
+
+#%option G_OPT_R_OUTPUT
+#%end
+
+#%flag
+#% key: t
+#% description: Do not assign the space time raster dataset start and end time to the output map
+#%end
+
+#%flag
+#% key: z
+#% description: Use zero (0) for transparency instead of NULL
+#%end
+
+#%flag
+#% key: s
+#% description: Do not create color and category files
+#%end
+
+#%option
+#% key: sort
+#% description: Sort order (see sort parameter)
+#% options: asc,desc
+#% answer: desc
+#%end
+
+
+import grass.script as grass
+from grass.exceptions import CalledModuleError
+
+
+def main():
+    # lazy imports
+    import grass.temporal as tgis
+
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    where = options["where"]
+    sort = options["sort"]
+    add_time = flags["t"]
+    patch_s = flags["s"]
+    patch_z = flags["z"]
+
+    # Make sure the temporal database exists
+    tgis.init()
+
+    sp = tgis.open_old_stds(input, "strds")
+
+    rows = sp.get_registered_maps("id", where, "start_time", None)
+
+    if rows:
+
+        ordered_rasts = []
+        # newest images are first
+        if sort == 'desc':
+            rows_sorted = rows[::-1]
+        # older images are first
+        elif sort == 'asc':
+            rows_sorted = rows
+
+        for row in rows_sorted:
+            string = "%s" % (row["id"])
+            ordered_rasts.append(string)
+
+        patch_flags = ""
+        if patch_z:
+            patch_flags += "z"
+        if patch_s:
+            patch_flags += "s"
+
+        try:
+            grass.run_command("r.patch",
+                              overwrite=grass.overwrite(),
+                              input=(',').join(ordered_rasts),
+                              output=output,
+                              flags=patch_flags
+                              )
+        except CalledModuleError:
+            grass.fatal(_("%s failed. Check above error messages.") % 'r.patch')
+
+        if not add_time:
+
+            # We need to set the temporal extent from the subset of selected maps
+            maps = sp.get_registered_maps_as_objects(where=where, order="start_time", dbif=None)
+            first_map = maps[0]
+            last_map = maps[-1]
+            start_a, end_a = first_map.get_temporal_extent_as_tuple()
+            start_b, end_b = last_map.get_temporal_extent_as_tuple()
+
+            if end_b is None:
+                end_b = start_b
+
+            if first_map.is_time_absolute():
+                extent = tgis.AbsoluteTemporalExtent(start_time=start_a, end_time=end_b)
+            else:
+                extent = tgis.RelativeTemporalExtent(start_time=start_a, end_time=end_b,
+                                                     unit=first_map.get_relative_time_unit())
+
+            # Create the time range for the output map
+            if output.find("@") >= 0:
+                id = output
+            else:
+                mapset = grass.gisenv()["MAPSET"]
+                id = output + "@" + mapset
+
+            map = sp.get_new_map_instance(id)
+            map.load()
+
+            map.set_temporal_extent(extent=extent)
+
+            # Register the map in the temporal database
+            if map.is_in_db():
+                map.update_all()
+            else:
+                map.insert()
+
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()