Implemented a simple r.series wrapper for space time raster datasets.

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

lib/python/temporal/abstract_datasets.py

@@ -464,6 +464,34 @@ class abstract_space_time_dataset(abstract_dataset):
         self.metadata.set_title(title)
         self.metadata.set_description(description)
 
+    def get_registered_maps(self, dbif=None, where = None):
+        """Return sqlite rows with ids of all registered maps. In case nothing found None is returned"""
+
+        connect = False
+
+        if dbif == None:
+            dbif = sql_database_interface()
+            dbif.connect()
+            connect = True
+
+        rows = None
+
+        if self.get_map_register():
+            sql = "SELECT id FROM " + self.get_map_register()
+            if where:
+                sql += " WHERE %s" % (where)
+            try:
+                dbif.cursor.execute(sql)
+                rows = dbif.cursor.fetchall()
+            except:
+                core.error("Unable to get map ids from register table <" + self.get_map_register() + ">")
+                raise
+
+        if connect == True:
+            dbif.close()
+
+        return rows
+
     def delete(self, dbif=None):
         """Delete a space time dataset from the database"""
         # First we need to check if maps are registered in this dataset and
@@ -479,22 +507,19 @@ class abstract_space_time_dataset(abstract_dataset):
             connect = True
 
         if self.get_map_register():
-            sql = "SELECT id FROM " + self.get_map_register()
+            rows = self.get_registered_maps(dbif)
+            # Unregister each registered map in the table
+            if rows:
+                for row in rows:
+                    # Unregister map
+                    map = self.get_new_map_instance(row["id"])
+                    self.unregister_map(map, dbif)
             try:
-                dbif.cursor.execute(sql)
-                rows = dbif.cursor.fetchall()
-                # Unregister each registered map in the table
-                if rows:
-                    for row in rows:
-                        # Unregister map
-                        map = self.get_new_map_instance(row["id"])
-                        self.unregister_map(map, dbif)
-
                 # Drop remove the map register table
                 sql = "DROP TABLE " + self.get_map_register()
                 dbif.cursor.execute(sql)
             except:
-                core.error("Unable to unregister maps from register table <" + self.get_map_register() + ">")
+                core.error("Unable to drop table <" + self.get_map_register() + ">")
                 raise
 
         # Remove the primary key, the foreign keys will be removed by trigger
@@ -513,7 +538,6 @@ class abstract_space_time_dataset(abstract_dataset):
             In case the map is already registered this function will break with a warning
             and return False
         """
-
         connect = False
 
         if dbif == None:
@@ -563,6 +587,8 @@ class abstract_space_time_dataset(abstract_dataset):
             # Create a unique id
             uuid_rand = "map_" + str(uuid.uuid4()).replace("-", "")
 
+            map_register_table = uuid_rand + "_" + self.get_type() + "_register"
+            
             # Read the SQL template
             sql = open(os.path.join(sql_path, "map_stds_register_table_template.sql"), 'r').read()
             # Create the raster, raster3d and vector tables
@@ -574,11 +600,20 @@ class abstract_space_time_dataset(abstract_dataset):
             try:
                 dbif.cursor.executescript(sql)
             except:
-                core.error("Unable to create the space time " + map.get_type() +\
-                " dataset register table for " + map.get_type() + " map <" + map.get_id())
-                raise
+                try:
+                    # Drop stds register table
+                    sql = "DROP TABLE " + map_register_table
+                    dbif.cursor.execute(sql)
+                except:
+                    core.error(_("Unable to drop table <%s>" % (map_register_table)))
+                    raise
+                try:
+                    dbif.cursor.executescript(sql_script)
+                except:
+                    core.error("Unable to create the space time " + map.get_type() +\
+                    " dataset register table for " + map.get_type() + " map <" + map.get_id())
+                    raise
 
-            map_register_table = uuid_rand + "_" + self.get_type() + "_register"
             # Set the stds register table name and put it into the DB
             map.set_stds_register(map_register_table)
             map.metadata.update(dbif)
@@ -587,6 +622,8 @@ class abstract_space_time_dataset(abstract_dataset):
 
         # We need to create the table and register it
         if stds_register_table == None:
+            # Create table name
+            stds_register_table = stds_name + "_" + stds_mapset + "_" + map.get_type() + "_register"
             # Read the SQL template
             sql = open(os.path.join(sql_path, "stds_map_register_table_template.sql"), 'r').read()
             # Create the raster, raster3d and vector tables
@@ -603,9 +640,19 @@ class abstract_space_time_dataset(abstract_dataset):
             try:
                 dbif.cursor.executescript(sql_script)
             except:
-                core.error("Unable to create the " + map.get_type() +\
-                " map register table for space time " + map.get_type() + " dataset <" + map.get_id())
-                raise
+                try:
+                    # Drop map register table
+                    sql = "DROP TABLE " + stds_register_table
+                    dbif.cursor.execute(sql)
+                except:
+                    core.error(_("Unable to drop table <%s>" % (stds_register_table)))
+                    raise
+                try:
+                    dbif.cursor.executescript(sql_script)
+                except:
+                    core.error("Unable to create the " + map.get_type() +\
+                    " map register table for space time " + map.get_type() + " dataset <" + map.get_id())
+                    raise
 
             # Trigger have been disabled due to peformance issues while registration
             ## We need raster specific trigger
@@ -623,8 +670,6 @@ class abstract_space_time_dataset(abstract_dataset):
 
             #dbif.cursor.executescript(sql_script)
 
-            stds_register_table = stds_name + "_" + stds_mapset + "_" + map.get_type() + "_register"
-
             # Set the map register table name and put it into the DB
             self.set_map_register(stds_register_table)
             self.metadata.update(dbif)

temporal/Makefile

@@ -8,6 +8,8 @@ SUBDIRS = \
 	t.time.abs \
 	t.time.rel \
 	tr.register \
+	tr.series \
+	tr.extract \
 	tr3.register \
 	tv.register \
 	tr.unregister \

temporal/tr.extract/Makefile

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

temporal/tr.extract/test.tr.extract.sh

@@ -0,0 +1,61 @@
+# This is a test to register and unregister raster maps in
+# space time raster dataset.
+# The raster maps will be registered in different space time raster
+# datasets
+
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster with r.mapcalc and create two space time raster datasets
+# with relative and absolute time
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+r.mapcalc --o expr="prec_1 = rand(0, 550)"
+r.mapcalc --o expr="prec_2 = rand(0, 450)"
+r.mapcalc --o expr="prec_3 = rand(0, 320)"
+r.mapcalc --o expr="prec_4 = rand(0, 510)"
+r.mapcalc --o expr="prec_5 = rand(0, 300)"
+r.mapcalc --o expr="prec_6 = rand(0, 650)"
+
+# The first @test
+# We create the space time raster datasets and register the raster maps with absolute time interval
+
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs1 gran="1 senconds" title="A test" descr="A test"
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs2 gran="1 minutes" title="A test" descr="A test"
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs3 gran="1 hours" title="A test" descr="A test"
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs4 gran="1 days" title="A test" descr="A test"
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs5 gran="1 weeks" title="A test" descr="A test"
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs6 gran="1 months" title="A test" descr="A test"
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs7 gran="1 years" title="A test" descr="A test"
+
+tr.register --v dataset=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="1 seconds"
+t.info type=strds dataset=precip_abs1
+tr.unregister --v dataset=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.info type=strds dataset=precip_abs1
+
+tr.register --v dataset=precip_abs2 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="20 seconds, 5 minutes"
+t.info type=strds dataset=precip_abs2
+
+tr.register --v dataset=precip_abs3 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="8 hours"
+t.info type=strds dataset=precip_abs3
+tr.unregister --v maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.info type=strds dataset=precip_abs3
+
+tr.register dataset=precip_abs4 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="3 days"
+t.info type=strds dataset=precip_abs4
+
+tr.register dataset=precip_abs5 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="4 weeks"
+t.info type=strds dataset=precip_abs5
+
+tr.register dataset=precip_abs6 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-08-01" increment="2 months"
+t.info type=strds dataset=precip_abs6
+
+tr.register dataset=precip_abs7 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="20 years, 3 months, 1 days, 4 hours"
+t.info type=strds dataset=precip_abs7
+# Register with different valid time again
+tr.register dataset=precip_abs7 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="99 years, 9 months, 9 days, 9 hours"
+t.info type=strds dataset=precip_abs7
+
+t.remove --v type=raster dataset=prec_1,prec_2,prec_3
+t.remove --v type=strds dataset=precip_abs1,precip_abs2,precip_abs3,precip_abs4,precip_abs5,precip_abs6,precip_abs7
+t.remove --v type=raster dataset=prec_4,prec_5,prec_6

temporal/tr.extract/tr.extract.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:	tr.extract
+# AUTHOR(S):	Soeren Gebbert
+#
+# PURPOSE:	Register raster maps in a space time raster dataset
+# COPYRIGHT:	(C) 2011 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: Register raster maps in a space time raster dataset
+#% keywords: spacetime raster dataset
+#% keywords: raster
+#% keywords: extract
+#%end
+
+#%option
+#% key: input
+#% type: string
+#% description: Name of an existing space time raster dataset
+#% required: yes
+#% multiple: yes
+#%end
+
+#%option G_OPT_DB_WHERE
+#%end
+
+#%option
+#% key: expression
+#% type: string
+#% description: The r.mapcalc expression assigned to all extracted raster maps
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% description: Name of the output space time raster dataset
+#% required: yes
+#% multiple: no
+#%end
+
+
+import grass.script as grass
+import grass.temporal as tgis
+
+############################################################################
+
+def main():
+
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    where = options["where"]
+    expression = options["expression"]
+
+    # Make sure the temporal database exists
+    tgis.create_temporal_database()
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()
+

temporal/tr.series/Makefile

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

temporal/tr.series/test.tr.series.sh

@@ -0,0 +1,32 @@
+# This is a test to register and unregister raster maps in
+# space time raster dataset.
+# The raster maps will be registered in different space time raster
+# datasets
+
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster with r.mapcalc and create two space time raster datasets
+# with relative and absolute time
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+r.mapcalc --o expr="prec_1 = rand(0, 550)"
+r.mapcalc --o expr="prec_2 = rand(0, 450)"
+r.mapcalc --o expr="prec_3 = rand(0, 320)"
+r.mapcalc --o expr="prec_4 = rand(0, 510)"
+r.mapcalc --o expr="prec_5 = rand(0, 300)"
+r.mapcalc --o expr="prec_6 = rand(0, 650)"
+
+# The first @test
+# We create the space time raster datasets and register the raster maps with absolute time interval
+
+t.create --v --o type=strds temporaltype=absolute dataset=precip_abs gran="1 months" title="A test" descr="A test"
+
+tr.register --v dataset=precip_abs maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="1 months"
+
+tr.series --o input=precip_abs method=average output=prec_average
+tr.series --o input=precip_abs method=maximum output=prec_max
+tr.series --o input=precip_abs method=quantile output=prec_quant quantile=0.9
+
+#t.remove --v type=raster dataset=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+#t.remove --v type=strds dataset=precip_abs

temporal/tr.series/tr.series.py

@@ -0,0 +1,108 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:	tr.series
+# AUTHOR(S):	Soeren Gebbert
+#
+# PURPOSE:	Performe different aggregation algorithms from r.series on all raster maps of a space time raster dataset
+# COPYRIGHT:	(C) 2011 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: Performe different aggregation algorithms from r.series on all raster maps of a space time raster dataset
+#% keywords: spacetime raster dataset
+#% keywords: raster
+#% keywords: extract
+#%end
+
+#%option
+#% key: input
+#% type: string
+#% description: Name of an existing space time raster dataset
+#% required: yes
+#% multiple: yes
+#%end
+
+#%option
+#% key: method
+#% type: string
+#% description: Aggregate operation
+#% required: yes
+#% multiple: no
+#% options: average,count,median,mode,minimum,min_raster,maximum,max_raster,stddev,range,sum,variance,diversity,slope,offset,detcoeff,quart1,quart3,perc90,quantile,skewness,kurtosis
+#% answer: average
+#%end
+
+#%option
+#% key: quantile
+#% type: string
+#% description: Quantile to calculate for method=quantile
+#% required: no
+#% multiple: no
+#% options: 0.0-1.0
+#%end
+
+#%option G_OPT_R_OUTPUT
+#%end
+
+
+import grass.script as grass
+import grass.temporal as tgis
+
+############################################################################
+
+def main():
+
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    quantile = options["quantile"]
+    method = options["method"]
+
+    # Make sure the temporal database exists
+    tgis.create_temporal_database()
+
+    if input.find("@") >= 0:
+        id = input
+    else:
+        mapset =  grass.gisenv()["MAPSET"]
+        id = input + "@" + mapset
+
+    sp = tgis.space_time_raster_dataset(id)
+
+    if sp.is_in_db() == False:
+        grass.fatal("Dataset <%s> not found in temporal database" % (id))
+
+    sp.select()
+
+    rows = sp.get_registered_maps()
+
+    if rows:
+        inputs = ""
+
+        count = 0
+        for row in rows:
+            if count == 0:
+                inputs += row["id"]
+            else:
+                inputs += "," + row["id"]
+            count += 1
+
+        print inputs
+
+        if grass.overwrite() == True:
+            grass.run_command("r.series", input=inputs, output=output, overwrite=True, method=method,  quantile=quantile)
+        else:
+            grass.run_command("r.series", input=inputs, output=output, overwrite=False, method=method,  quantile=quantile)
+
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()
+