temporal modules: New module that creates a stvds with contours from a strds

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@64559 15284696-431f-4ddb-bdfa-cd5b030d7da7

temporal/Makefile

@@ -17,6 +17,7 @@ SUBDIRS = \
 	t.rast.aggregate \
 	t.rast.aggregate.ds \
 	t.rast.colors \
+	t.rast.contour \
 	t.rast.to.rast3 \
 	t.rast.univar \
 	t.rast.list \

temporal/t.rast.contour/Makefile

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

temporal/t.rast.contour/t.rast.contour.html

@@ -0,0 +1,43 @@
+<h2>DESCRIPTION</h2>
+
+<em>t.rast.contour</em> is designed to produce a space time vector dataset of 
+specified contours from a space time raster dataset. This module works as a front-end to
+<a href="r.contour.html">r.contour</a> and therefore supports all parameter
+of this module. Hence, all raster map layers in a space time raster dataset
+are successively passed to <a href="r.contour.html">r.contour</a> that computes 
+the contour lines. 
+Please refer to the <a href="r.contour.html">r.contour</a> documentation 
+for a detailed description.
+The new generated vector contour map
+layers will be registered in the output space time vector dataset, using 
+the same time stamps as their raster map layer origins.
+<p>
+This module supports the parallel processing of <a href="r.contour.html">r.contour</a>
+module instances. The number of parallel processes 
+can be set with the <em>nprocs</em> option. 
+However, this will only work in conjunction with the <em>-t</em>
+flag, that avoids the creation of attribute tables. 
+The parallel creation of attribute tables is not supported.
+<p>
+The <em>where</em> option allows to select subsets of the input space time raster
+dataset.
+<p>
+The flag <em>-n</em> can be used to force the registration of empty 
+vector map layers. Empty vector maps may occur in case that empty 
+raster map layers should be converted into vector map layers, or in case the
+chosen steps or contour levels are not present in the raster map layers.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.contour.html">r.contour</a>,
+<a href="t.vect.db.select.html">t.rast.db.select</a>,
+<a href="t.info.html">t.info</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+S&ouml;ren Gebbert, Geoinformatikb&uuml;ro Dassau
+
+<p><i>Last changed: $Date: 2014-12-03 20:39:26 +0100 (Mi, 03. Dez 2014) $</i>

temporal/t.rast.contour/t.rast.contour.py

@@ -0,0 +1,265 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:       t.rast.contour
+# AUTHOR(S):    Soeren Gebbert
+#
+# PURPOSE:      Produces a space time vector dataset of specified contours from a space time raster dataset.
+#
+# COPYRIGHT:    (C) 2015 by the GRASS Development Team
+#
+#               This program is free software under the GNU General Public
+#               License (version 2). Read the file COPYING that comes with GRASS
+#               for details.
+#
+#############################################################################
+
+#%module
+#% description: Produces a space time vector dataset of specified contours from a space time raster dataset.
+#% keyword: temporal
+#% keyword: contour
+#% keyword: raster
+#% keyword: vector
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#%end
+
+#%option G_OPT_STVDS_OUTPUT
+#%end
+
+#%option G_OPT_T_WHERE
+#%end
+
+#%option
+#% key: basename
+#% type: string
+#% label: Basename of the new generated output maps
+#% description: A numerical suffix separated by an underscore will be attached to create a unique identifier
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: step
+#% type: double
+#% description: Increment between contour levels
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: levels
+#% type: double
+#% description: List of contour levels
+#% required: no
+#% multiple: yes
+#%end
+
+#%option
+#% key: minlevel
+#% type:  double
+#% description: Minimum contour level
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: maxlevel
+#% type: double
+#% description: Maximum contour level
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: cut
+#% type: integer
+#% description: Minimum number of points for a contour line (0 -> no limit)
+#% required: no
+#% multiple: no
+#% answer: 0
+#%end
+
+#%option
+#% key: nprocs
+#% type: integer
+#% description: Number of r.contour processes to run in parallel, more than 1 process works only in conjunction with flag -t
+#% required: no
+#% multiple: no
+#% answer: 1
+#%end
+
+#%flag
+#% key: n
+#% description: Register empty vector maps
+#%end
+
+#%flag
+#% key: t
+#% description: Do not create attribute tables
+#%end
+
+
+import sys
+import copy
+import grass.script as gscript
+import grass.temporal as tgis
+import grass.pygrass.modules as pymod
+from grass.exceptions import FatalError
+
+
+############################################################################
+
+def main(options, flags):
+
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    where = options["where"]
+    base = options["basename"]
+    nprocs = int(options["nprocs"])
+    step = options["step"]
+    levels = options["levels"]
+    minlevel = options["minlevel"]
+    maxlevel = options["maxlevel"]
+    cut = options["cut"]
+
+    register_null = flags["n"]
+    t_flag = flags["t"]
+    
+
+    # Make sure the temporal database exists
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+
+    overwrite = gscript.overwrite()
+
+    sp = tgis.open_old_stds(input, "strds", dbif)
+    maps = sp.get_registered_maps_as_objects(where=where, dbif=dbif)
+
+    if not maps:
+        dbif.close()
+        gscript.warning(_("Space time raster dataset <%s> is empty") % sp.get_id())
+        return
+
+    # Check the new stvds
+    new_sp = tgis.check_new_stds(output, "stvds", dbif=dbif,
+                                 overwrite=overwrite)
+                                               
+    # Setup the flags
+    flags = ""
+    if t_flag is True:
+        flags += "t"
+    
+    # Configure the r.to.vect module
+    contour_module = pymod.Module("r.contour", input="dummy",
+                                   output="dummy", run_=False,
+                                   finish_=False, flags=flags,
+                                   overwrite=overwrite,
+                                   quiet=True)
+
+    if step:
+        contour_module.inputs.step = float(step)
+    if minlevel:
+        contour_module.inputs.minlevel = float(minlevel)
+    if maxlevel:
+        contour_module.inputs.maxlevel = float(maxlevel)
+    if levels:
+        contour_module.inputs.levels = levels.split(",")
+    if cut:
+        contour_module.inputs.cut = int(cut)
+
+    # The module queue for parallel execution, except if attribute tables should
+    # be created. Then force single process use
+    if t_flag is False:
+        if nprocs > 1:
+            nprocs = 1
+            gscript.warning(_("The number of parellel r.contour processes was "\
+                              "reduced to 1 because of the table attribute "\
+                              "creation"))
+    process_queue = pymod.ParallelModuleQueue(int(nprocs))
+
+    count = 0
+    num_maps = len(maps)
+    new_maps = []
+
+    # run r.to.vect all selected maps
+    for map in maps:
+        count += 1
+        map_name = "%s_%i" % (base, count)
+        new_map = tgis.open_new_map_dataset(map_name, None, type="vector",
+                                            temporal_extent=map.get_temporal_extent(),
+                                            overwrite=overwrite, dbif=dbif)
+        new_maps.append(new_map)
+
+        mod = copy.deepcopy(contour_module)
+        mod(input=map.get_id(), output=new_map.get_id())
+        sys.stderr.write(mod.get_bash() + "\n")
+        process_queue.put(mod)
+
+        if count%10 == 0:
+            gscript.percent(count, num_maps, 1)
+
+    # Wait for unfinished processes
+    process_queue.wait()
+
+    # Open the new space time vector dataset
+    ttype, stype, title, descr = sp.get_initial_values()
+    new_sp = tgis.open_new_stds(output, "stvds", ttype, title,
+                                descr, stype, dbif, overwrite)
+    # collect empty maps to remove them
+    num_maps = len(new_maps)
+    empty_maps = []
+
+    # Register the maps in the database
+    count = 0
+    for map in new_maps:
+        count += 1
+
+        if count%10 == 0:
+            gscript.percent(count, num_maps, 1)
+
+        # Do not register empty maps
+        try:
+            if map.load() is not True:
+                continue
+        except FatalError:
+            continue
+        if map.metadata.get_number_of_primitives() == 0:
+            if not register_null:
+                empty_maps.append(map)
+                continue
+
+        # Insert map in temporal database
+        map.insert(dbif)
+        new_sp.register_map(map, dbif)
+
+    # Update the spatio-temporal extent and the metadata table entries
+    new_sp.update_from_registered_maps(dbif)
+    gscript.percent(1, 1, 1)
+
+    # Remove empty maps
+    if len(empty_maps) > 0:
+        names = ""
+        count = 0
+        for map in empty_maps:
+            if count == 0:
+                count += 1
+                names += "%s" % (map.get_name())
+            else:
+                names += ",%s" % (map.get_name())
+
+        gscript.run_command("g.remove", flags='f', type='vector', name=names, 
+                            quiet=True)
+
+    dbif.close()
+
+############################################################################
+
+if __name__ == "__main__":
+    options, flags = gscript.parser()
+    main(options, flags)

temporal/t.rast.contour/testsuite/test_convert.py

@@ -0,0 +1,221 @@
+"""Test t.rast.contour
+
+(C) 2015 by the GRASS Development Team
+This program is free software under the GNU General Public
+License (>=v2). Read the file COPYING that comes with GRASS
+for details.
+
+@author Soeren Gebbert
+"""
+
+import subprocess
+from grass.gunittest.case import TestCase
+from grass.gunittest.gmodules import SimpleModule
+
+class TestRasterContour(TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        """Initiate the temporal GIS and set the region
+        """
+        cls.use_temp_region()
+        cls.runModule("g.region", s=0, n=80, w=0, e=120, b=0, t=50, res=10, res3=10)
+
+        cls.runModule("r.mapcalc", expression="a_1 = row()", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_2 = col()", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_3 = row()", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_4 = col()", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_5 = null()", overwrite=True)
+
+        cls.runModule("t.create", type="strds", temporaltype="absolute",  
+                                 output="A", title="A test", description="A test",  
+                                 overwrite=True)
+        cls.runModule("t.register",  flags="i", type="raster", input="A",  
+                                     maps="a_1,a_2,a_3,a_4,a_5", start="2001-01-01", 
+                                     increment="3 months", overwrite=True)
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region
+        """
+        cls.runModule("t.remove", flags="rf", type="strds",  
+                                   inputs="A")
+        cls.del_temp_region()
+        
+    def tearDown(self):
+        """Remove generated data"""
+        self.runModule("t.remove", flags="rf", type="stvds",  
+                                   inputs="result")
+
+    def test_register_empty_maps(self):
+        self.assertModule("t.rast.contour", input="A", output="result", 
+                          levels=(1,2,3,4,5,6,7,8), flags="n",
+                          basename="test",
+                          nprocs=1, overwrite=True, verbose=True)
+
+        #self.assertModule("t.info",  type="stvds", flags="g",  input="result")
+
+        tinfo_string="""start_time=2001-01-01 00:00:00
+                        end_time=2002-04-01 00:00:00
+                        granularity=3 months
+                        map_time=interval
+                        number_of_maps=5
+                        lines=28
+                        primitives=28
+                        nodes=56"""
+
+        info = SimpleModule("t.info", flags="g", type="stvds", input="result")
+        self.assertModuleKeyValue(module=info, reference=tinfo_string, precision=2, sep="=")
+
+    def test_simple(self):
+        """Do not register empty maps"""
+        self.assertModule("t.rast.contour", input="A", output="result", 
+                          levels=(1,2,3,4,5,6,7,8),
+                          basename="test",
+                          nprocs=1, overwrite=True, verbose=True)
+
+        #self.assertModule("t.info",  type="stvds", flags="g",  input="result")
+
+        tinfo_string="""start_time=2001-01-01 00:00:00
+                        end_time=2002-01-01 00:00:00
+                        granularity=3 months
+                        map_time=interval
+                        number_of_maps=4
+                        lines=28
+                        primitives=28
+                        nodes=56"""
+
+        info = SimpleModule("t.info", flags="g", type="stvds", input="result")
+        self.assertModuleKeyValue(module=info, reference=tinfo_string, precision=2, sep="=")
+
+    def test_where(self):
+        """Use where statement and do not register empty maps"""
+        self.assertModule("t.rast.contour", input="A", output="result", 
+                          levels=(1,2,3,4,5,6,7,8),
+                          basename="test",
+                          where="start_time > '2001-02-01'",
+                          nprocs=1, overwrite=True, verbose=True)
+
+        #self.assertModule("t.info",  type="stvds", flags="g",  input="result")
+
+        tinfo_string="""start_time=2001-04-01 00:00:00
+                        end_time=2002-01-01 00:00:00
+                        granularity=3 months
+                        map_time=interval
+                        number_of_maps=3
+                        lines=21
+                        primitives=21
+                        nodes=42"""
+
+        info = SimpleModule("t.info", flags="g", type="stvds", input="result")
+        self.assertModuleKeyValue(module=info, reference=tinfo_string, precision=2, sep="=")
+
+    def test_parallel(self):
+        """Run 4 contour processes do not create attribute tables"""
+        self.assertModule("t.rast.contour", input="A", output="result", 
+                          levels=(1,2,3,4,5,6,7,8),
+                          basename="test", flags="t",
+                          nprocs=4, overwrite=True, verbose=True)
+
+        #self.assertModule("t.info",  type="stvds", flags="g",  input="result")
+
+        tinfo_string="""start_time=2001-01-01 00:00:00
+                        end_time=2002-01-01 00:00:00
+                        granularity=3 months
+                        map_time=interval
+                        number_of_maps=4
+                        lines=28
+                        primitives=28
+                        nodes=56"""
+
+        info = SimpleModule("t.info", flags="g", type="stvds", input="result")
+        self.assertModuleKeyValue(module=info, reference=tinfo_string, precision=2, sep="=")
+
+    def test_parallel_cut(self):
+        """Do not register empty maps"""
+        self.assertModule("t.rast.contour", input="A", output="result", 
+                          levels=(1,2,3,4,5,6,7,8), cut=2,
+                          basename="test", flags="t",
+                          nprocs=4, overwrite=True, verbose=True)
+
+        #self.assertModule("t.info",  type="stvds", flags="g",  input="result")
+
+        tinfo_string="""start_time=2001-01-01 00:00:00
+                        end_time=2002-01-01 00:00:00
+                        granularity=3 months
+                        map_time=interval
+                        number_of_maps=4
+                        lines=28
+                        primitives=28
+                        nodes=56"""
+
+        info = SimpleModule("t.info", flags="g", type="stvds", input="result")
+        self.assertModuleKeyValue(module=info, reference=tinfo_string, precision=2, sep="=")
+
+    def test_where_step(self):
+        """Use where statement and do not register empty maps"""
+        self.assertModule("t.rast.contour", input="A", output="result", 
+                          step=1, minlevel=1, maxlevel=8,
+                          basename="test",
+                          where="start_time > '2001-02-01'",
+                          nprocs=1, overwrite=True, verbose=True)
+
+        #self.assertModule("t.info",  type="stvds", flags="g",  input="result")
+
+        tinfo_string="""start_time=2001-04-01 00:00:00
+                        end_time=2002-01-01 00:00:00
+                        granularity=3 months
+                        map_time=interval
+                        number_of_maps=3
+                        lines=21
+                        primitives=21
+                        nodes=42"""
+
+        info = SimpleModule("t.info", flags="g", type="stvds", input="result")
+        self.assertModuleKeyValue(module=info, reference=tinfo_string, precision=2, sep="=")
+
+
+class TestRasterContourFails(TestCase):
+    @classmethod
+    def setUpClass(cls):
+        """Initiate the temporal GIS and set the region
+        """
+        cls.use_temp_region()
+        cls.runModule("g.region", s=0, n=80, w=0, e=120, b=0, t=50, res=10, res3=10)
+
+        cls.runModule("r.mapcalc", expression="a_1 = 100", overwrite=True)
+
+        cls.runModule("t.create",  type="strds", temporaltype="absolute",  
+                                 output="A", title="A test", description="A test",  
+                                 overwrite=True)
+        cls.runModule("t.register",  flags="i", type="raster", input="A",  
+                                     maps="a_1", start="2001-01-01", 
+                                     increment="3 months", overwrite=True)
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region
+        """
+        cls.runModule("t.remove", flags="rf", type="strds",
+                      inputs="A")
+        cls.del_temp_region()
+
+    def test_error_handling(self):
+        # No basename
+        self.assertModuleFail("t.rast.contour", input="A", output="result", 
+                              step=1, minlevel=1, maxlevel=8,
+                              where="start_time > '2001-02-01'",
+                              nprocs=1, overwrite=True, verbose=True)
+
+
+    def test_empty_strds(self):
+        """Test for empty strds"""
+        self.assertModule("t.rast.contour", input="A", output="result",
+                          basename="test",
+                          step=1, minlevel=1, maxlevel=8,
+                          where="start_time > '2010-01-01'",
+                          nprocs=1, overwrite=True, verbose=True)
+
+if __name__ == '__main__':
+    from grass.gunittest.main import test
+    test()