v.what.strds: backported from grass71, https://trac.osgeo.org/grass/changeset/65686 and https://trac.osgeo.org/grass/changeset/65704

git-svn-id: https://svn.osgeo.org/grass/grass/branches/releasebranch_7_0@65705 15284696-431f-4ddb-bdfa-cd5b030d7da7

scripts/v.what.strds/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.what.strds
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

scripts/v.what.strds/testsuite/test_what_strds.py

@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+"""Test v.what.strds
+
+(C) 2014 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 Luca Delucchi
+"""
+
+from grass.gunittest.case import TestCase
+from grass.gunittest.gmodules import SimpleModule
+import grass.script as gscript
+
+
+class TestWhatStrds(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("r.mapcalc", expression="a_2 = 200",  overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_3 = 300",  overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_4 = 400",  overwrite=True)
+
+        cls.runModule("v.random", output="points", npoints=3, seed=1,
+                      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", 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_output(self):
+        self.assertModule("v.what.strds", input="points", strds="A",
+                          output="what_strds", overwrite=True)
+
+        maps = gscript.list_strings('vector')
+        self.assertIn('what_strds@{ma}'.format(ma=gscript.gisenv()['MAPSET']),
+                      maps)
+
+    def test_values(self):
+        self.assertModule("v.what.strds", input="points", strds="A",
+                          output="what_strds", overwrite=True)
+        db_sel = SimpleModule("v.db.select", map="what_strds")
+        self.assertModule(db_sel)
+        output = """cat|A_2001_01_01|A_2001_04_01|A_2001_07_01|A_2001_10_01
+1|100|200|300|400
+2|100|200|300|400
+3|100|200|300|400
+"""
+        self.assertLooksLike(output, db_sel.outputs.stdout)
+
+if __name__ == '__main__':
+    from grass.gunittest.main import test
+    test()

scripts/v.what.strds/v.what.strds.html

@@ -0,0 +1,28 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.what.strds</em> retrieves raster values from a given space-time raster datasets
+(STRDS) using a point vector map.
+
+<h2>NOTES</h2>
+
+TBD.
+
+<h2>EXAMPLES</h2>
+
+<div class="code"><pre>
+v.what.strds input=mypoints strds=mystrds output=newvector
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.what.rast.html">v.what.rast</a>,
+<a href="t.vect.observe.strds.html">t.vect.observe.strds</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Luca Delucchi
+
+<p>
+<i>Last changed: $Date$</i>

scripts/v.what.strds/v.what.strds.py

@@ -0,0 +1,233 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:       v.what.strds
+# AUTHOR(S):    Luca delucchi
+#
+# PURPOSE:      Uploads space time raster dataset values at positions of vector points to the table
+# COPYRIGHT:    (C) 2013 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: Uploads space time raster dataset values at positions of vector points to the table.
+#% keyword: vector
+#% keyword: temporal
+#% keyword: sampling
+#% keyword: position
+#% keyword: querying
+#% keyword: attribute table
+#%end
+
+#%option G_OPT_V_INPUT
+#%end
+
+#%option G_OPT_STRDS_INPUTS
+#% key: strds
+#%end
+
+#%option G_OPT_V_OUTPUT
+#%end
+
+#%option G_OPT_DB_WHERE
+#%end
+
+#%option G_OPT_T_WHERE
+#% key: t_where
+#%end
+
+import grass.script as grass
+import grass.temporal as tgis
+from grass.pygrass.utils import copy as gcopy
+from grass.pygrass.messages import Messenger
+from grass.pygrass.vector import Vector
+from grass.exceptions import CalledModuleError
+
+############################################################################
+
+
+class Sample(object):
+    def __init__(self, start=None, end=None, raster_names=None,
+                 strds_name=None):
+        self.start = start
+        self.end = end
+        if raster_names is not None:
+            self.raster_names = raster_names
+        else:
+            self.raster_names = []
+        self.strds_name = strds_name
+
+    def __str__(self):
+        return "Start: %s\nEnd: %s\nNames: %s\n" % (str(self.start),
+                                                    str(self.end),
+                                                    str(self.raster_names))
+
+    def printDay(self, date='start'):
+        if date == 'start':
+            return str(self.start).split(' ')[0].replace('-', '_')
+        elif date == 'end':
+            return str(self.end).split(' ')[0].replace('-', '_')
+        else:
+            grass.fatal("The values accepted by printDay in Sample are:"
+                        " 'start', 'end'")
+
+############################################################################
+
+
+def main():
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    strds = options["strds"]
+    where = options["where"]
+    tempwhere = options["t_where"]
+
+    if where == "" or where == " " or where == "\n":
+        where = None
+
+    overwrite = grass.overwrite()
+
+    quiet = True
+
+    if grass.verbosity() > 2:
+        quiet = False
+
+    # Check the number of sample strds and the number of columns
+    strds_names = strds.split(",")
+
+    # Make sure the temporal database exists
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+
+    samples = []
+
+    first_strds = tgis.open_old_stds(strds_names[0], "strds", dbif)
+    # Single space time raster dataset
+    if len(strds_names) == 1:
+        rows = first_strds.get_registered_maps("name,mapset,start_time,end_time",
+                                               tempwhere, "start_time",
+                                               dbif)
+
+        if not rows:
+            dbif.close()
+            grass.fatal(_("Space time raster dataset <%s> is empty") %
+                        first_strds.get_id())
+
+        for row in rows:
+            start = row["start_time"]
+            end = row["end_time"]
+            raster_maps = [row["name"] + "@" + row["mapset"], ]
+
+            s = Sample(start, end, raster_maps, first_strds.get_name())
+            samples.append(s)
+    else:
+        # Multiple space time raster datasets
+        for name in strds_names[1:]:
+            dataset = tgis.open_old_stds(name, "strds", dbif)
+            if dataset.get_temporal_type() != first_strds.get_temporal_type():
+                grass.fatal(_("Temporal type of space time raster "
+                              "datasets must be equal\n<%(a)s> of type "
+                              "%(type_a)s do not match <%(b)s> of type "
+                              "%(type_b)s" % {"a": first_strds.get_id(),
+                               "type_a": first_strds.get_temporal_type(),
+                               "b": dataset.get_id(),
+                               "type_b": dataset.get_temporal_type()}))
+
+        mapmatrizes = tgis.sample_stds_by_stds_topology("strds", "strds",
+                                                        strds_names,
+                                                        strds_names[0],
+                                                        False, None,
+                                                        "equal", False,
+                                                        False)
+
+        for i in xrange(len(mapmatrizes[0])):
+            isvalid = True
+            mapname_list = []
+            for mapmatrix in mapmatrizes:
+
+                entry = mapmatrix[i]
+
+                if entry["samples"]:
+                    sample = entry["samples"][0]
+                    name = sample.get_id()
+                    if name is None:
+                        isvalid = False
+                        break
+                    else:
+                        mapname_list.append(name)
+
+            if isvalid:
+                entry = mapmatrizes[0][i]
+                map = entry["granule"]
+
+                start, end = map.get_temporal_extent_as_tuple()
+                s = Sample(start, end, mapname_list, name)
+                samples.append(s)
+
+    # Get the layer and database connections of the input vector
+    gcopy(input, output, 'vect')
+
+    msgr = Messenger()
+    perc_curr = 0
+    perc_tot = len(samples)
+    pymap = Vector(output)
+    try:
+        pymap.open('r')
+    except:
+        dbif.close()
+        grass.fatal(_("It is not possible to open the new map %s" % output))
+
+    if len(pymap.dblinks) == 0:
+        try:
+            grass.run_command("v.db.addtable", map=output)
+        except CalledModuleError:
+            dbif.close()
+            grass.fatal(_("Impossible add table to vector %s" % output))
+    pymap.close()
+    for sample in samples:
+        raster_names = sample.raster_names
+        # Call v.what.rast for each raster map
+
+        for name in raster_names:
+            coltype = "DOUBLE PRECISION"
+            # Get raster map type
+            raster_map = tgis.RasterDataset(name)
+            raster_map.load()
+            if raster_map.metadata.get_datatype() == "CELL":
+                coltype = "INT"
+            day = sample.printDay()
+            column_name = "%s_%s" % (sample.strds_name, day)
+            column_string = "%s %s" % (column_name, coltype)
+            column_string.replace('.', '_')
+            try:
+                grass.run_command("v.db.addcolumn", map=output,
+                                  column=column_string,
+                                  overwrite=overwrite)
+            except CalledModuleError:
+                dbif.close()
+                grass.fatal(_("Unable to add column %s to vector map "
+                              "<%s> ") % (column_string, output))
+            try:
+                grass.run_command("v.what.rast", map=output, raster=name,
+                                  column=column_name, where=where,
+                                  quiet=quiet)
+            except CalledModuleError:
+                dbif.close()
+                grass.fatal(_("Unable to run v.what.rast for vector map"
+                              " <%s> and raster map <%s>") %
+                            (output, str(raster_names)))
+
+        msgr.percent(perc_curr, perc_tot, 1)
+        perc_curr += 1
+
+    dbif.close()
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()