Fixed map removal in t.remove

Added multiple strds observation to t.vect.observe.strds
Added flag to t.rast.series


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

temporal/t.rast.neighbors/t.rast.neighbors.py

@@ -110,6 +110,7 @@ def main():
                                    finish_=False, size=int(size),
                                    method=method, overwrite=overwrite,
                                    quiet=True)
+
     # The module queue for parallel execution
     process_queue = pymod.ParallelModuleQueue(int(nprocs))
 

temporal/t.rast.series/t.rast.series.py

@@ -55,6 +55,12 @@
 #% description: Do not assign the space time raster dataset start and end time to the output map
 #%end
 
+#%flag
+#% key: n
+#% description: Propagate NULLs
+#%end
+
+
 import grass.script as grass
 import grass.temporal as tgis
 
@@ -70,6 +76,7 @@ def main():
     order = options["order"]
     where = options["where"]
     add_time = flags["t"]
+    nulls = flags["n"]
 
     # Make sure the temporal database exists
     tgis.init()
@@ -89,7 +96,11 @@ def main():
 
         file.close()
 
-        ret = grass.run_command("r.series", flags="z", file=filename,
+        flag = "z"
+        if nulls:
+            flag += "n"
+
+        ret = grass.run_command("r.series", flags=flag, file=filename,
                                 output=output, overwrite=grass.overwrite(),
                                 method=method)
 

temporal/t.remove/t.remove.py

@@ -120,7 +120,11 @@ def main():
             name_list = []
             for map in maps:
                 map.select(dbif)
-                name_list.append(map.get_name())
+                # We may have multiple layer for a single map, hence we need
+                # to avoid multiple deletation of the same map,
+                # but the database entries are still present and must be removed
+                if map.get_name() not in name_list:
+                    name_list.append(map.get_name())
                 map_statement += map.delete(dbif=dbif, execute=False)
 
                 count += 1
@@ -134,7 +138,12 @@ def main():
             if map_statement:
                 dbif.execute_transaction(map_statement)
             if name_list:
-                remove(rast=name_list, run_=True)
+                if type == "strds":
+                    remove(rast=name_list, run_=True)
+                if type == "stvds":
+                    remove(vect=name_list, run_=True)
+                if type == "str3ds":
+                    remove(rast3d=name_list, run_=True)
 
         statement += sp.delete(dbif=dbif, execute=False)
 

temporal/t.vect.observe.strds/t.vect.observe.strds.py

@@ -23,7 +23,7 @@
 #%option G_OPT_V_INPUT
 #%end
 
-#%option G_OPT_STRDS_INPUT
+#%option G_OPT_STRDS_INPUTS
 #% key: strds
 #%end
 
@@ -32,18 +32,17 @@
 
 #%option G_OPT_V_OUTPUT
 #% key: vector_output
-#% description: Name of the new created vector map that stores the sampled values
+#% description: Name of the new created vector map that stores the sampled values in different layers
 #%end
 
 #%option
-#% key: column
+#% key: columns
 #% type: string
-#% description: Name of the vector column to be created and to store sampled raster values, otherwise the space time raster name is used as column name
-#% required: no
-#% multiple: no
+#% description: Names of the vector columns to be created and to store sampled raster values, one name for each STRDS
+#% required: yes
+#% multiple: yes
 #%end
 
-
 #%option G_OPT_DB_WHERE
 #%end
 
@@ -57,6 +56,19 @@ import grass.script.raster as raster
 
 ############################################################################
 
+class Sample(object):
+    def __init__(self, start = None, end = None, raster_names = None):
+        self.start = start
+        self.end = end
+        if raster_names != None:
+            self.raster_names = raster_names
+        else:
+            self.raster_names = []
+
+    def __str__(self):
+        return "Start: %s\nEnd: %s\nNames: %s\n"%(str(self.start), str(self.end), str(self.raster_names))
+            
+############################################################################
 
 def main():
 
@@ -66,7 +78,7 @@ def main():
     vector_output = options["vector_output"]
     strds = options["strds"]
     where = options["where"]
-    column = options["column"]
+    columns = options["columns"]
     tempwhere = options["t_where"]
 
     if where == "" or where == " " or where == "\n":
@@ -74,6 +86,13 @@ def main():
 
     overwrite = grass.overwrite()
 
+    # Check the number of sample strds and the number of columns
+    strds_names = strds.split(",")
+    column_names = columns.split(",")
+
+    if len(strds_names) != len(column_names):
+        grass.fatal(_("The number of columns must be equal to the number of space time raster datasets"))
+
     # Make sure the temporal database exists
     tgis.init()
     # We need a database interface
@@ -82,27 +101,70 @@ def main():
 
     mapset = grass.gisenv()["MAPSET"]
 
-    strds_sp = tgis.open_old_space_time_dataset(strds, "strds", dbif)
+    out_sp = tgis.check_new_space_time_dataset(output, "stvds", dbif, overwrite)
 
+    samples = []
 
-    title = _("Observaion of space time raster dataset <%s>") % (strds_sp.get_id())
-    description= _("Observation of space time raster dataset <%s>"
-                                " with vector map <%s>") % (strds_sp.get_id(),
-                                                            input)
+    first_strds = tgis.open_old_space_time_dataset(strds_names[0], "strds", dbif)
 
-    out_sp = tgis.check_new_space_time_dataset(output, "stvds", dbif,
-                                               overwrite)
+    # 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)
 
-    # Select the raster maps
-    rows = strds_sp.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") %
+                        out_sp.get_id())
 
-    if not rows:
-        dbif.close()
-        grass.fatal(_("Space time raster dataset <%s> is empty") %
-                    out_sp.get_id())
+        for row in rows:
+            start = row["start_time"]
+            end = row["end_time"]
+            raster_maps = [row["name"] + "@" + row["mapset"],]
 
-    num_rows = len(rows)
+            s = Sample(start, end, raster_maps)
+            samples.append(s)
+    else:
+        # Multiple space time raster datasets
+        for name in strds_names[1:]:
+            dataset = tgis.open_old_space_time_dataset(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)
+                samples.append(s)
+
+    num_samples = len(samples)
 
     # Get the layer and database connections of the input vector
     vector_db = grass.vector.vector_db(input)
@@ -114,7 +176,7 @@ def main():
     else:
         layers = ""
     first = True
-    for layer in range(num_rows):
+    for layer in range(num_samples):
         layer += 1
         # Skip existing layer
         if vector_db and layer in vector_db and \
@@ -136,11 +198,15 @@ def main():
         grass.fatal(_("Unable to create new layers for vector map <%s>")
                     % (vectmap))
 
+    title = _("Observaion of space time raster dataset(s) <%s>") % (strds)
+    description= _("Observation of space time raster dataset(s) <%s>"
+                   " with vector map <%s>") % (strds, input)
+
     # Create the output space time vector dataset
     out_sp = tgis.open_new_space_time_dataset(output, "stvds",
-                                              strds_sp.get_temporal_type(),
+                                              first_strds.get_temporal_type(),
                                               title, description,
-                                              strds_sp.get_semantic_type(),
+                                              first_strds.get_semantic_type(),
                                               dbif, overwrite)
 
     dummy = out_sp.get_new_map_instance(None)
@@ -148,64 +214,72 @@ def main():
     # Sample the space time raster dataset with the vector
     # map at specific layer with v.what.rast
     count = 1
-    for row in rows:
-        start = row["start_time"]
-        end = row["end_time"]
-        rastmap = row["name"] + "@" + row["mapset"]
-
-        if column:
-            col_name = column
-        else:
-            # Create a new column with name of the
-            # sampled space time raster dataset
-            col_name = row["name"]
-
-        coltype = "DOUBLE PRECISION"
-        # Get raster map type
-        rasterinfo = raster.raster_info(rastmap)
-        if rasterinfo["datatype"] == "CELL":
-            coltype = "INT"
+    for sample in samples:
+        raster_names = sample.raster_names
+
+        if len(raster_names) != len(column_names):
+            grass.fatal(_("The number of raster maps in a granule must "
+                          "be equal to the number of column names"))
+
+        # Create the columns creation string
+        columns_string = ""
+        for name, column in zip(raster_names, column_names):
+            # The column is by default double precision
+            coltype = "DOUBLE PRECISION"
+            # Get raster map type
+            raster_map = tgis.RasterDataset(name)
+            raster_map.load()
+            if raster_map.metadata.get_datatype() == "CELL":
+                coltype = "INT"
+
+            tmp_string = "%s %s,"%(column, coltype)
+            columns_string += tmp_string
+
+        # Remove last comma
+        columns_string = columns_string[0:len(columns_string) - 1]
 
         # Try to add a column
         if vector_db and count in vector_db and vector_db[count]["table"]:
             ret = grass.run_command("v.db.addcolumn", map=vectmap,
-                                    layer=count,
-                                    column="%s %s" % (col_name, coltype),
+                                    layer=count, column=columns_string,
                                     overwrite=overwrite)
             if ret != 0:
                 dbif.close()
                 grass.fatal(_("Unable to add column %s to vector map <%s> "
-                              "with layer %i") % (col_name, vectmap, count))
+                              "with layer %i") % (columns_string, vectmap, count))
         else:
             # Try to add a new table
             grass.message("Add table to layer %i" % (count))
             ret = grass.run_command("v.db.addtable", map=vectmap, layer=count,
-                                    columns="%s %s" % (col_name, coltype),
-                                    overwrite=overwrite)
+                                    columns=columns_string, overwrite=overwrite)
             if ret != 0:
                 dbif.close()
                 grass.fatal(_("Unable to add table to vector map "
                               "<%s> with layer %i") % (vectmap, count))
 
-        # Call v.what.rast
-        ret = grass.run_command("v.what.rast", map=vectmap,
-                                layer=count, raster=rastmap,
-                                column=col_name, where=where)
-        if ret != 0:
-            dbif.close()
-            grass.fatal(_("Unable to run v.what.rast for vector map <%s> "
-                          "with layer %i and raster map <%s>") % \
-                        (vectmap, count, rastmap))
+        # Call v.what.rast for each raster map
+        for name, column in zip(raster_names, column_names):
+            ret = grass.run_command("v.what.rast", map=vectmap,
+                                    layer=count, raster=name,
+                                    column=column, where=where)
+            if ret != 0:
+                dbif.close()
+                grass.fatal(_("Unable to run v.what.rast for vector map <%s> "
+                            "with layer %i and raster map <%s>") % \
+                            (vectmap, count, str(raster_names)))
 
         vect = out_sp.get_new_map_instance(dummy.build_id(vectmap,
                                                           mapset, str(count)))
         vect.load()
-
+        
+        start = sample.start
+        end = sample.end
+        
         if out_sp.is_time_absolute():
             vect.set_absolute_time(start, end)
         else:
             vect.set_relative_time(
-                start, end, strds_sp.get_relative_time_unit())
+                start, end, first_strds.get_relative_time_unit())
 
         if vect.is_in_db(dbif):
             vect.update_all(dbif)

temporal/t.vect.observe.strds/test.t.vect.observe.strds.sh

@@ -14,12 +14,20 @@ r.mapcalc  expr="prec_4 = 400"
 r.mapcalc  expr="prec_5 = 500.0"
 r.mapcalc  expr="prec_6 = 600.0"
 
+r.mapcalc  expr="prec_7 = 400"
+r.mapcalc  expr="prec_8 = 500.0"
+r.mapcalc  expr="prec_9 = 600.0"
+
+
 v.random output=prec n=5 seed=1
 v.random -z output=test_1 column=test n=5 seed=1 
 
 t.create type=strds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
 t.register -i input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-03-01 00:00:00" increment="1 months"
 
+t.create type=strds temporaltype=absolute output=precip_abs2 title="A test" descr="A test"
+t.register -i input=precip_abs2 maps=prec_7,prec_8,prec_9 start="2001-05-01 00:00:00" increment="1 months"
+
 # The @test
 t.vect.observe.strds input=prec strds=precip_abs1 output=prec_observer vector=prec_observer column="test_val"
 v.info prec_observer
@@ -27,19 +35,15 @@ t.info type=stvds input=prec_observer
 t.vect.list input=prec_observer
 t.vect.db.select input=prec_observer
 
-t.vect.observe.strds input=test_1 strds=precip_abs1 output=test_1_observer  vector=test_1_observer
+t.vect.observe.strds columns=test1,test2,test3 input=test_1 \
+    strds=precip_abs1,precip_abs1,precip_abs2 output=test_1_observer  \
+    vector=test_1_observer
+
 v.info test_1_observer
 t.info type=stvds input=test_1_observer
 t.vect.list input=test_1_observer
-t.vect.db.select input=test_1_observer
-
+t.vect.db.select input=test_1_observer columns=cat,test1,test2,test3
 
 # @postprocess
-t.unregister type=rast maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
-t.remove type=strds input=precip_abs1
-t.remove type=stvds input=prec_observer,test_1_observer
-t.unregister type=vect maps=prec_observer:1,prec_observer:2,prec_observer:3,prec_observer:4,prec_observer:5,prec_observer:6
-t.unregister type=vect maps=test_1_observer:1,test_1_observer:2,test_1_observer:3,test_1_observer:4,test_1_observer:5,test_1_observer:6
-
-g.remove vect=prec_observer,test_1_observer
-g.remove rast=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.remove -rf type=strds input=precip_abs1,precip_abs2
+t.remove -rf type=stvds input=prec_observer,test_1_observer