New module to fill gaps in space time raster datasets using raster to raster interpolation of r.series.interp.

Keys for temporal sorting of dataset lists by start or end time.
Many new unit tests about sorting and temporal topology building.
Temporal topology computation.


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

lib/python/array.py

@@ -7,6 +7,8 @@ Functions to use GRASS rasters with NumPy.
 Usage:
 
 @code
+from grass.script import array as garray
+
     map = 'elevation'
     x = garray.array()
     x.read(map)

lib/python/temporal/abstract_dataset.py

@@ -74,6 +74,15 @@ class abstract_dataset(object):
 	    self.relative_time.print_self()
 	self.spatial_extent.print_self()
 	self.metadata.print_self()
+	
+    def set_id(self, ident):
+	self.base.set_id(ident)
+	if self.is_time_absolute():
+	    self.absolute_time.set_id(ident)
+        if self.is_time_relative():
+	    self.relative_time.set_id(ident)
+	self.spatial_extent.set_id(ident)
+	self.metadata.set_id(ident)
 
     def get_id(self):
 	"""!Return the unique identifier of the dataset"""
@@ -325,3 +334,82 @@ class abstract_dataset(object):
 	    return self.relative_time.temporal_relation(map.relative_time)
     	return None
 
+###############################################################################
+	
+class abstract_dataset_comparison_key_start_time(object):
+    """!This comparison key can be used to sort lists of abstract datasets by start time
+    
+        Example:
+        
+        # Return all maps in a space time raster dataset as map objects
+	map_list = strds.get_registered_maps_as_objects()
+	
+	# Sort the maps in the list by start time
+	sorted_map_list = sorted(map_list, key=abstract_dataset_comparison_key_start_time)
+    """
+    def __init__(self, obj, *args):
+	self.obj = obj
+    def __lt__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return startA < startB
+    def __gt__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return startA > startB
+    def __eq__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return startA == startB
+    def __le__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return startA <= startB
+    def __ge__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return startA >= startB
+    def __ne__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return startA != startB
+	
+###############################################################################
+	
+class abstract_dataset_comparison_key_end_time(object):
+    """!This comparison key can be used to sort lists of abstract datasets by end time
+    
+        Example:
+        
+        # Return all maps in a space time raster dataset as map objects
+	map_list = strds.get_registered_maps_as_objects()
+	
+	# Sort the maps in the list by end time
+	sorted_map_list = sorted(map_list, key=abstract_dataset_comparison_key_end_time)
+    """
+    def __init__(self, obj, *args):
+	self.obj = obj
+    def __lt__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return endA < endB
+    def __gt__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return endA > endB
+    def __eq__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return endA == endB
+    def __le__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return endA <= endB
+    def __ge__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return endA >= endB
+    def __ne__(self, other):
+	startA, endA = self.obj.get_valid_time()
+	startB, endB = other.obj.get_valid_time()
+	return endA != endB

lib/python/temporal/abstract_map_dataset.py

@@ -28,7 +28,7 @@ class temporal_map_relations(abstract_dataset):
     
        This object will be set up by temporal topology creation methods.
        
-       If correctly initialize the calls get_next() and get_prev() let the user walk temporally forward
+       If correctly initialize the calls next() and prev() let the user walk temporally forward
        and backward in time.
        
        The following temporal relations with access methods are supported:
@@ -51,7 +51,7 @@ class temporal_map_relations(abstract_dataset):
            for _map in dlist:
                _map.print_info()
                
-           start = start.get_next()       
+           start = start.next()       
     """
 
     def __init__(self):
@@ -81,7 +81,7 @@ class temporal_map_relations(abstract_dataset):
 	   temporally located AFTER the start time of this map, but temporally 
 	   near than other maps of the same dataset.
 	   
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	"""
         self._temporal_topology["NEXT"] = _map
 
@@ -92,7 +92,7 @@ class temporal_map_relations(abstract_dataset):
 	   temporally located BEFORE the start time of this map, but temporally 
 	   near than other maps of the same dataset.
 	   
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	"""
         self._temporal_topology["PREV"] = _map
 
@@ -119,7 +119,7 @@ class temporal_map_relations(abstract_dataset):
     def append_equivalent(self, _map): 
 	"""!Append a map with equivalent temporal extent as this map
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	"""
         if not self._temporal_topology.has_key("EQUAL"):
             self._temporal_topology["EQUAL"] = []
@@ -137,7 +137,7 @@ class temporal_map_relations(abstract_dataset):
     def append_overlaps(self, _map): 
 	"""!Append a map that this map temporally overlaps
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	"""
         if not self._temporal_topology.has_key("OVERLAPS"):
             self._temporal_topology["OVERLAPS"] = []
@@ -155,7 +155,7 @@ class temporal_map_relations(abstract_dataset):
     def append_overlapped(self, _map): 
 	"""!Append a map that this map temporally overlapped
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	""" 
         if not self._temporal_topology.has_key("OVERLAPPED"):
             self._temporal_topology["OVERLAPPED"] = []
@@ -173,7 +173,7 @@ class temporal_map_relations(abstract_dataset):
     def append_follows(self, _map): 
 	"""!Append a map that this map temporally follows
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	""" 
         if not self._temporal_topology.has_key("FOLLOWS"):
             self._temporal_topology["FOLLOWS"] = []
@@ -191,7 +191,7 @@ class temporal_map_relations(abstract_dataset):
     def append_precedes(self, _map): 
 	"""!Append a map that this map temporally precedes
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	""" 
         if not self._temporal_topology.has_key("PRECEDES"):
             self._temporal_topology["PRECEDES"] = []
@@ -210,7 +210,7 @@ class temporal_map_relations(abstract_dataset):
 	"""!Append a map that this map is temporally located during
 	   This includes temporal relationships starts and finishes
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	""" 
         if not self._temporal_topology.has_key("DURING"):
             self._temporal_topology["DURING"] = []
@@ -230,7 +230,7 @@ class temporal_map_relations(abstract_dataset):
 	"""!Append a map that this map temporally contains
 	   This includes temporal relationships started and finished
 	
-	   @param _map This object should be of type abstract_map_dataset or derived classes
+	   @param _map: This object should be of type abstract_map_dataset or derived classes
 	""" 
         if not self._temporal_topology.has_key("CONTAINS"):
             self._temporal_topology["CONTAINS"] = []

lib/python/temporal/abstract_space_time_dataset.py

@@ -246,7 +246,7 @@ class abstract_space_time_dataset(abstract_dataset):
         """
 
         if maps == None:
-            maps = get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
+            maps = self.get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
             
         gaps = 0
 
@@ -259,7 +259,7 @@ class abstract_space_time_dataset(abstract_dataset):
 
         return gaps
         
-    def print_temporal_relation_matrix(self, maps=None, dbif=None):
+    def print_temporal_relationships(self, maps=None, dbif=None):
         """!Print the temporal relation matrix of all registered maps to stdout
 
            The temporal relation matrix includes the temporal relations between
@@ -270,25 +270,9 @@ class abstract_space_time_dataset(abstract_dataset):
         """
         
         if maps == None:
-            maps = get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
+            maps = self.get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
 	    
-	print_temporal_relations(maps, maps)
-
-    def get_temporal_relation_matrix(self, maps=None, dbif=None):
-        """!Return the temporal relation matrix of all registered maps as list of lists
-
-           The map list must be ordered by start time
-
-           The temporal relation matrix includes the temporal relations between
-           all registered maps. The relations are strings stored in a list of lists.
-           
-           @param maps: a ordered by start_time list of map objects
-           @param dbif: The database interface to be used
-        """
-        if maps == None:
-            maps = get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
-
-        return get_temporal_relation_matrix(maps, maps)
+	print_temporal_topology_relationships(maps, maps)
 
     def count_temporal_relations(self, maps=None, dbif=None):
         """!Count the temporal relations between the registered maps.
@@ -302,9 +286,9 @@ class abstract_space_time_dataset(abstract_dataset):
         """
         
         if maps == None:
-            maps = get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
+            maps = self.get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
 
-        return count_temporal_relations(maps, maps)
+        return count_temporal_topology_relationships(maps, maps)
 
     def check_temporal_topology(self, maps=None, dbif=None):
         """!Check the temporal topology
@@ -335,9 +319,9 @@ class abstract_space_time_dataset(abstract_dataset):
            @return True if topology is correct
         """
         if maps == None:
-            maps = get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
+            maps = self.get_registered_maps_as_objects(where=None, order="start_time", dbif=dbif)
 
-        relations = self.count_temporal_relations(maps)
+        relations = count_temporal_topology_relationships(maps, maps)
 
         map_time = self.get_map_time()
 

lib/python/temporal/base.py

@@ -438,6 +438,16 @@ class dataset_base(sql_database_interface):
         """
 	self.ident = ident
 	self.D["id"] = ident
+	
+        if ident != None:
+            if ident.find("@") >= 0:
+                name, mapset = ident.split("@")
+		self.set_mapset(mapset)
+            if name.find(":") >= 0:
+                name, layer = ident.split(":")
+		self.set_layer(layer)
+	    
+	    self.set_name(name)
 
     def set_name(self, name):
 	"""!Set the name of the dataset

lib/python/temporal/space_time_datasets.py

@@ -27,6 +27,7 @@ import grass.lib.gis as libgis
 import grass.lib.raster as libraster
 import grass.lib.vector as libvector
 import grass.lib.raster3d as libraster3d
+import grass.script.array as garray
 
 from datetime_math import *
 from abstract_map_dataset import *
@@ -73,11 +74,30 @@ class raster_dataset(abstract_map_dataset):
         """Return the two dimensional spatial relation"""
         
         return self.spatial_extent.spatial_relation_2d(dataset.spatial_extent)
+        
+    def get_np_array(self):
+	"""Return this raster map as memmap numpy style array to access the raster
+	   values in numpy style without loading the whole map in the RAM. 
+	   
+	   In case this raster map does exists in the grass spatial database, the map
+	   will be exported using r.out.bin to a temporary location and assigned to the
+	   memmap object that is returned by this function.
+	   
+	   In case the raster map does not exists, an empty temporary binary file will be created
+	   and assigned to the memap object.
+	   
+	   You need to call the write function to write the memmap array back into grass.
+	"""
+	
+	a = garray.array()
+	
+	if self.map_exists():
+	    a.read(self.get_map_id())
+	
+	return a
 	
     def reset(self, ident):
 	"""!Reset the internal structure and set the identifier"""
-	self.ident = ident
-
 	self.base = raster_base(ident=ident)
 	self.absolute_time = raster_absolute_time(ident=ident)
 	self.relative_time = raster_relative_time(ident=ident)
@@ -210,9 +230,6 @@ class raster_dataset(abstract_map_dataset):
 
 
         # Fill base information
-
-        self.base.set_name(self.ident.split("@")[0])
-        self.base.set_mapset(self.ident.split("@")[1])
         self.base.set_creator(str(getpass.getuser()))
 
         # Get the data from an existing raster map
@@ -289,8 +306,6 @@ class raster3d_dataset(abstract_map_dataset):
         
     def reset(self, ident):
 	"""!Reset the internal structure and set the identifier"""
-	self.ident = ident
-
 	self.base = raster3d_base(ident=ident)
 	self.absolute_time = raster3d_absolute_time(ident=ident)
 	self.relative_time = raster3d_relative_time(ident=ident)
@@ -421,8 +436,6 @@ class raster3d_dataset(abstract_map_dataset):
         """!Load all info from an existing raster3d map into the internal structure"""
 
         # Fill base information
-        self.base.set_name(self.ident.split("@")[0])
-        self.base.set_mapset(self.ident.split("@")[1])
         self.base.set_creator(str(getpass.getuser()))
 
         # Fill spatial extent
@@ -502,8 +515,6 @@ class vector_dataset(abstract_map_dataset):
 	
     def reset(self, ident):
 	"""!Reset the internal structure and set the identifier"""
-	self.ident = ident
-
 	self.base = vector_base(ident=ident)
 	self.absolute_time = vector_absolute_time(ident=ident)
 	self.relative_time = vector_relative_time(ident=ident)
@@ -648,12 +659,6 @@ class vector_dataset(abstract_map_dataset):
 
 
         # Fill base information
-	if self.ident.find(":") >= 0:
-	    self.base.set_name(self.ident.split("@")[0].split(":")[0])
-	    self.base.set_layer(self.ident.split("@")[0].split(":")[1])
-	else:
-	    self.base.set_name(self.ident.split("@")[0])
-        self.base.set_mapset(self.ident.split("@")[1])
         self.base.set_creator(str(getpass.getuser()))
 
         # Get the data from an existing raster map
@@ -719,13 +724,7 @@ class space_time_raster_dataset(abstract_space_time_dataset):
     def reset(self, ident):
 
 	"""!Reset the internal structure and set the identifier"""
-	self.ident = ident
-
 	self.base = strds_base(ident=ident)
-
-        if ident != None:
-            self.base.set_name(self.ident.split("@")[0])
-            self.base.set_mapset(self.ident.split("@")[1])
         self.base.set_creator(str(getpass.getuser()))
         self.absolute_time = strds_absolute_time(ident=ident)
         self.relative_time = strds_relative_time(ident=ident)
@@ -779,13 +778,7 @@ class space_time_raster3d_dataset(abstract_space_time_dataset):
     def reset(self, ident):
 
 	"""!Reset the internal structure and set the identifier"""
-	self.ident = ident
-
 	self.base = str3ds_base(ident=ident)
-
-        if ident != None:
-            self.base.set_name(self.ident.split("@")[0])
-            self.base.set_mapset(self.ident.split("@")[1])
         self.base.set_creator(str(getpass.getuser()))
         self.absolute_time = str3ds_absolute_time(ident=ident)
         self.relative_time = str3ds_relative_time(ident=ident)
@@ -833,13 +826,7 @@ class space_time_vector_dataset(abstract_space_time_dataset):
     def reset(self, ident):
 
 	"""!Reset the internal structure and set the identifier"""
-	self.ident = ident
-
 	self.base = stvds_base(ident=ident)
-
-        if ident != None:
-            self.base.set_name(self.ident.split("@")[0])
-            self.base.set_mapset(self.ident.split("@")[1])
         self.base.set_creator(str(getpass.getuser()))
         self.absolute_time = stvds_absolute_time(ident=ident)
         self.relative_time = stvds_relative_time(ident=ident)

lib/python/temporal/spatial_extent.py

@@ -23,7 +23,7 @@ for details.
 from base import *
 
 class spatial_extent(sql_database_interface):
-    """!This is the spatial extent base class for all maps and spacetime datasets"""
+    """!This is the spatial extent base class for all maps and space time datasets"""
     def __init__(self, table=None, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None, proj="XY"):
 
 	sql_database_interface.__init__(self, table, ident)
@@ -61,16 +61,16 @@ class spatial_extent(sql_database_interface):
                 E -= 360.0
                 W -= 360.0
                 
-        if(self.get_north() < S):
+        if(self.get_north() <= S):
             return False
         
-        if(self.get_south() > N):
+        if(self.get_south() >= N):
             return False
             
-        if self.get_east() < W:
+        if self.get_east() <= W:
             return False
         
-        if self.get_west() > E:
+        if self.get_west() >= E:
             return False
         
         return True
@@ -92,10 +92,10 @@ class spatial_extent(sql_database_interface):
         T = extent.get_top()
         B = extent.get_bottom()
         
-        if self.get_top() < B:
+        if self.get_top() <= B:
             return False
         
-        if self.get_bottom() > T:
+        if self.get_bottom() >= T:
             return False
         
         return True
@@ -678,7 +678,7 @@ class spatial_extent(sql_database_interface):
 	return True
 
     def meet(self,extent):
-	""" Check if self and extent touch meet other in three dimensions"""
+	""" Check if self and extent meet other in three dimensions"""
 	eN = extent.get_north()
         eS = extent.get_south()
         eE = extent.get_east()
@@ -727,7 +727,7 @@ class spatial_extent(sql_database_interface):
 	    edge = "B"
 	    edge_count += 1	
 	
-	# Meet a a single edge only
+	# Meet a single edge only
 	if edge_count != 1:
 	    return False
 	

lib/python/temporal/temporal_relationships.py

@@ -23,8 +23,10 @@ for details.
 from abstract_map_dataset import *
 from datetime_math import *
 
+###############################################################################
+
 class temporal_topology_builder(object):
-    """!This class is designed to build the temporal topology based on a list of maps
+    """!This class is designed to build the temporal topology based on a lists of maps
     
 	Example:
 	
@@ -34,7 +36,7 @@ class temporal_topology_builder(object):
 	
 	# Now lets build the temporal topology of the maps in the list
 	tb = temporal_topology_builder()
-	tb.build(maps, True)
+	tb.build(maps)
 	
 	for _map in tb:
 	    _map.print_temporal_topology_info()
@@ -63,6 +65,7 @@ class temporal_topology_builder(object):
     def _reset(self):
         self._store = {}
         self._first = None
+        self._temporal_iteratable = False
 
     def _set_first(self, first):
         self._first = first
@@ -84,8 +87,117 @@ class temporal_topology_builder(object):
 	   @return The map with the earliest start time
 	"""
 	return self._first
+
+    def _build_internal_iteratable(self, maps):
+	"""!Build an iteratable temporal topology structure for all maps in the list and store the maps internally
+	
+	   Basically the "next" and "prev" relations will be set in the temporal topology structure of each map
+	   The maps will be added to the object, so they can be accessed using the iterator of this class
+	   
+	   @param maps: A sorted (by start_time)list of abstract_dataset objects with initiated temporal extent
+	"""
+	self._build_iteratable(maps)
+
+	for _map in maps:
+	    self._insert(_map)
+	
+	# Detect the first map
+	self._detect_first()
+	
+    def _build_iteratable(self, maps):
+	"""!Build an iteratable temporal topology structure for all maps in the list
+	
+	   Basically the "next" and "prev" relations will be set in the temporal topology structure of each map.
+	   
+	   @param maps: A sorted (by start_time)list of abstract_dataset objects with initiated temporal extent
+	"""
+	for i in xrange(len(maps)):
+	    offset = i + 1
+	    for j in xrange(offset, len(maps)):		
+		# Get the temporal relationship
+		relation = maps[j].temporal_relation(maps[i])
+		
+		# Build the next reference
+		if relation != "equivalent" and relation != "started":
+		    maps[i].set_next(maps[j])
+		    break
+		
+	for _map in maps:
+	    _next = _map.next()
+	    if _next:
+		_next.set_prev(_map)
+	    _map.set_temporal_topology_build_true()
+	
+    def build2(self, mapsA, mapsB):
+	"""!Build the temporal topology structure between two ordered lists of maps
+	
+	   This method builds the temporal topology from mapsA to mapsB and vice verse.
+	   The temporal topology structure of each map, defined in class temporal_map_relations,
+	   will be reseted and rebuild for mapsA and mapsB. 
+	   
+	   After building the temporal topology the modified map objects of mapsA can be accessed 
+	   in the same way as a dictionary using there id. The implemented iterator assures 
+	   the chronological iteration over the mapsA.	    
+	   
+	   @param mapsA: A sorted list (by start_time) of abstract_dataset objects with initiated temporal extent
+	   @param mapsB: A sorted list (by start_time) of abstract_dataset objects with initiated temporal extent
+	"""
+	
+	if mapsA == mapsB:
+	    self.build(mapsA, True)
+	    return
+	
+	for _map in mapsA:
+	    _map.reset_temporal_topology()
+	    
+	for _map in mapsB:
+	    _map.reset_temporal_topology()
+	
+	for i in xrange(len(mapsA)):
+	    for j in xrange(len(mapsB)):
+		
+		# Get the temporal relationship
+		relation = mapsB[j].temporal_relation(mapsA[i])
+		
+		if relation == "before":
+		    continue
+			    
+		if relation == "equivalent":
+		    mapsB[j].append_equivalent(mapsA[i])
+		    mapsA[i].append_equivalent(mapsB[j])
+		elif relation == "follows":
+		    mapsB[j].append_follows(mapsA[i])
+		    mapsA[i].append_precedes(mapsB[j])
+		elif relation == "precedes":
+		    mapsB[j].append_precedes(mapsA[i])
+		    mapsA[i].append_follows(mapsB[j])
+		elif relation == "during" or relation == "starts" or relation == "finishes":
+		    mapsB[j].append_during(mapsA[i])
+		    mapsA[i].append_contains(mapsB[j])
+		elif relation == "contains" or relation == "started" or relation == "finished":
+		    mapsB[j].append_contains(mapsA[i])
+		    mapsA[i].append_during(mapsB[j])
+		elif relation == "overlaps":
+		    mapsB[j].append_overlaps(mapsA[i])
+		    mapsA[i].append_overlapped(mapsB[j])
+		elif relation == "overlapped":
+		    mapsB[j].append_overlapped(mapsA[i])
+		    mapsA[i].append_overlaps(mapsB[j])
+
+		# Break if the next map follows and the over-next maps is after
+		if relation == "follows":
+		    if j < len(mapsB) - 1:
+			relation = mapsB[j + 1].temporal_relation(mapsA[i])
+			if relation == "after":
+			    break
+		# Break if the the next map is after
+		if relation == "after":
+		    break 
 	
-    def build(self, maps, is_sorted = False):
+	self._build_internal_iteratable(mapsA)
+	self._build_iteratable(mapsB)
+			    
+    def build(self, maps):
 	"""!Build the temporal topology structure
 	
 	   This method builds the temporal topology based on all maps in the provided map list.
@@ -94,58 +206,19 @@ class temporal_topology_builder(object):
 	   
 	   After building the temporal topology the modified map objects can be accessed 
 	   in the same way as a dictionary using there id. The implemented iterator assures 
-	   the chronological iteration over the maps.
+	   the chronological iteration over the maps.	   
 	   
-	   
-	   @param maps: A list of abstract_dataset objects with initiated temporal extent
-	   @param is_sorted Set to True if the map list is sorted by start time, sorting will dramatically reduce computation time
+	   @param maps: A sorted list (by start_time) of abstract_dataset objects with initiated temporal extent
 	"""
 	for _map in maps:
 	    _map.reset_temporal_topology()
 	
 	for i in xrange(len(maps)):
-	    if is_sorted:
-		offset = i + 1
-		found_next = False
-	    else:
-		offset = 0
-		# Needed for "next" computation 
-		start0, end0 = maps[i].get_valid_time()
-	    
+	    offset = i + 1
 	    for j in xrange(offset, len(maps)):
 		
-		# Do not build topology of the same maps
-		if maps[i] == maps[j]:
-		    continue
-		
 		# Get the temporal relationship
 		relation = maps[j].temporal_relation(maps[i])
-		
-		# We can skip several relationships if not sorted
-		if not is_sorted:
-		    if relation == "before":
-			continue
-		    if relation == "precedes":
-			continue
-		    if relation == "overlapped":
-			continue
-		    if relation == "finished":
-			continue
-		
-		# Build the next reference
-		if is_sorted:
-		    if not found_next and relation != "equivalent" and relation != "started":
-			maps[i].set_next(maps[j])
-			found_next = True
-		else: 
-		    start2, end2 = maps[j].get_valid_time()
-		    if maps[i].next():
-			start1, end1 = maps[i].next().get_valid_time()
-			if start2 > start0 and start2 < start1:
-			    maps[i].set_next(maps[j])
-		    else:
-			if start2 > start0:
-			    maps[i].set_next(maps[j])
 			    
 		# The start time of map j is equal or later than map i
 		if relation == "equivalent":
@@ -174,23 +247,14 @@ class temporal_topology_builder(object):
 		# Break if the the next map is after
 		if relation == "after":
 		    break 
-	
-	# Build the previous pointer and store the map internally
-	for _map in maps:
-	    _next = _map.next()
-	    if _next:
-		_next.set_prev(_map)
-	    _map.set_temporal_topology_build_true()
-	    self._insert(_map)
-	
-	# Detect the first map
-	self._detect_first()
+		    
+	self._build_internal_iteratable(maps)
 	
     def __iter__(self):
-        _start = self._first
-        while _start != None:
-            yield _start
-            _start = _start.next()
+	_start = self._first
+	while _start != None:
+	    yield _start
+	    _start = _start.next()
 
     def __getitem__(self, index):
         return self._store[index.get_id()]
@@ -202,12 +266,14 @@ class temporal_topology_builder(object):
         return _map in self._store.values()
 
 
-def print_temporal_relations(maps1, maps2):
+###############################################################################
+
+def print_temporal_topology_relationships(maps1, maps2):
     """!Print the temporal relation matrix of the temporal ordered map lists maps1 and maps2
        to stdout.
 	
-	@param maps1: a ordered by start_time list of map objects with initiated temporal extent
-	@param maps2: a ordered by start_time list of map objects with initiated temporal extent
+	@param maps1: A sorted (by start_time) list of abstract_dataset objects with initiated temporal extent
+	@param maps2: A sorted (by start_time) list of abstract_dataset objects with initiated temporal extent
     """
     
     identical = False
@@ -240,60 +306,16 @@ def print_temporal_relations(maps1, maps2):
 	    if relation == "after":
 		break
 
-def get_temporal_relation_matrix(maps1, maps2):
-    """!Return the temporal relation matrix of two map lists
-
-	Booth map lists must be ordered by start time
-
-	The temporal relationship matrix includes the temporal relations between
-	the two map lists. Returned is a nested dict representing 
-	a sparse (upper right side in case maps1 == maps2) relationship matrix.
-	
-	@param maps1: A sorted (start_time) list of abstract_dataset objects with initiated temporal extent
-	@param maps2: A sorted (start_time) list of abstract_dataset objects with initiated temporal extent
-    """
-
-    matrix = {}
-    identical = False
-    
-    if maps1 == maps2:
-	identical = True
-
-    for i in range(len(maps1)):
-	if identical == True:
-	    start = i + 1
-	else:
-	    start = 0
-	    
-	row = {}
-	    
-	for j in range(start, len(maps2)):
-	    relation = maps1[j].temporal_relation(maps2[i])
-
-	    row[maps2[j].base.get_id()] = relation 
-
-	    # Break if the last map follows
-	    if relation == "follows":
-		if j < len(maps1) - 1:
-		    relation = maps1[j + 1].temporal_relation(maps2[i])
-		    if relation == "after":
-			break
-	    # Break if the the next map is after
-	    if relation == "after":
-		break
-
-	matrix[maps1[i].base.get_id()] = row
-
-    return matrix
+###############################################################################
 
-def count_temporal_relations(maps1, maps2):
+def count_temporal_topology_relationships(maps1, maps2):
     """!Count the temporal relations between the two lists of maps
 
 	The map lists must be ordered by start time. Temporal relations are counted 
 	by analyzing the sparse (upper right side in case maps1 == maps2) temporal relationships matrix.
 
-	@param maps1: A sorted (start_time) list of abstract_dataset objects with initiated temporal extent
-	@param maps2: A sorted (start_time) list of abstract_dataset objects with initiated temporal extent
+	@param maps1: A sorted (by start_time) list of abstract_dataset objects with initiated temporal extent
+	@param maps2: A sorted (by start_time) list of abstract_dataset objects with initiated temporal extent
 	@return A dictionary with counted temporal relationships
     """
     
@@ -330,4 +352,3 @@ def count_temporal_relations(maps1, maps2):
 		break  
 
     return tcount
-

lib/python/temporal/univar_statistics.py

@@ -101,7 +101,7 @@ def print_gridded_dataset_univar_statistics(type, input, where, extended, header
 ###############################################################################
     
 def print_vector_dataset_univar_statistics(input, twhere, layer, type, column, where, extended, header, fs):
-    """!Print univariate statistics for a space time raster or raster3d dataset
+    """!Print univariate statistics for a space time vector dataset
     
        @param input The name of the space time dataset
        @param twhere A temporal database where statement

raster/Makefile

@@ -86,6 +86,7 @@ SUBDIRS = \
 	r.rescale.eq \
 	r.ros \
 	r.series \
+	r.series.interp \
 	r.slope.aspect \
 	r.solute.transport \
 	r.spread \

temporal/Makefile

@@ -20,6 +20,7 @@ SUBDIRS = \
 	t.rast.export \
 	t.rast.out.vtk \
 	t.rast.import \
+	t.rast.gapfill \
 	t.rast.extract \
 	t.rast3d.list \
 	t.rast3d.extract \

temporal/t.rast.aggregate/t.rast.aggregate.py

@@ -58,15 +58,6 @@
 #%option G_OPT_R_BASE
 #%end
 
-#%option
-#% key: nprocs
-#% type: integer
-#% description: The number of r.mapcalc processes to run in parallel
-#% required: no
-#% multiple: no
-#% answer: 2
-#%end
-
 #%flag
 #% key: n
 #% description: Register Null maps
@@ -89,7 +80,6 @@ def main():
     register_null = flags["n"]
     method = options["method"]
     sampling = options["sampling"]
-    nprocs = int(options["nprocs"])
 
     # Make sure the temporal database exists
     tgis.create_temporal_database()

temporal/t.rast.gapfill/Makefile

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

temporal/t.rast.gapfill/t.rast.gapfill.py

@@ -0,0 +1,176 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:	t.rast.gapfill
+# AUTHOR(S):	Soeren Gebbert
+#
+# PURPOSE:	Replace gaps in a space time raster dataset with interpolated raster maps.
+# COPYRIGHT:	(C) 2012 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: Replace gaps in a space time raster dataset with interpolated raster maps.
+#% keywords: temporal
+#% keywords: interpolation
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#%end
+
+#%option G_OPT_T_WHERE
+#%end
+
+#%option
+#% key: base
+#% type: string
+#% description: The base name of the interpolated raster maps. This name will be extended with a numerical prefix
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: nprocs
+#% type: integer
+#% description: The number of interpolation processes to run in parallel
+#% required: no
+#% multiple: no
+#% answer: 1
+#%end
+
+#%flag
+#% key: t
+#% description: Assign the space time raster dataset start and end time to the output map
+#%end
+
+from multiprocessing import Process
+import grass.script as grass
+import grass.temporal as tgis
+
+############################################################################
+
+    
+def main():
+
+    # Get the options
+    input = options["input"]
+    base = options["base"]
+    where = options["where"]
+    nprocs = options["nprocs"]
+
+    # Make sure the temporal database exists
+    tgis.create_temporal_database()
+
+    if input.find("@") >= 0:
+        id = input
+    else:
+        mapset =  grass.gisenv()["MAPSET"]
+        id = input + "@" + mapset
+
+    # We need a database interface
+    dbif = tgis.sql_database_interface_connection()
+    dbif.connect()
+    
+    sp = tgis.space_time_raster_dataset(id)
+
+    if sp.is_in_db(dbif) == False:
+        grass.fatal(_("Space time %s dataset <%s> not found") % (sp.get_new_map_instance(None).get_type(), id))
+
+    sp.select(dbif)
+    
+    maps = sp.get_registered_maps_as_objects_with_gaps(where, dbif)
+
+    num = len(maps)
+    
+    gap_list = []
+    
+    # Identify all gaps and create new names
+    count = 0
+    for _map in maps:
+	if _map.get_id() == None:
+	    count += 1
+	    _id = "%s_%d@%s"%(base, num + count, mapset)
+	    print _id
+	    _map.set_id(_id)
+	    if _map.map_exists() or _map.is_in_db(dbif):
+		grass.fatal(_("Map with name <%s> already exists. Please use another base name."%(_id)))
+	    
+	    gap_list.append(_map)
+	    
+    if len(gap_list) == 0:
+	grass.message(_("No gaps found"))
+	return
+
+    # Build the temporal topology
+    tb = tgis.temporal_topology_builder()
+    tb.build(maps)
+    
+    # Do some checks before computation
+    for _map in gap_list:	
+	if not _map.get_precedes() or not _map.get_follows():
+	    grass.fatal(_("Unable to determine successor and predecessor of a gap."))
+	
+	if len(_map.get_precedes()) > 1:
+	    grass.warning(_("More than one successor of the gap found. Using the first found."))
+	    
+	if len(_map.get_follows()) > 1:
+	    grass.warning(_("More than one predecessor of the gap found. Using the first found."))
+
+    # Interpolate the maps using parallel processing   
+    proc_list = []
+    proc_count = 0
+    num = len(gap_list)
+        
+    for _map in gap_list:	
+	predecessor = _map.get_follows()[0]
+	successor = _map.get_precedes()[0]
+	
+	# Build the module inputs strings
+	inputs = "%s,%s"%(predecessor.get_map_id(), successor.get_map_id())
+	dpos = "0,1"
+	output = "%s"%(_map.get_name())
+	outpos = "0.5"
+	
+	# Start several processes in parallel
+	proc_list.append(Process(target=run_interp, args=(inputs,dpos,output,outpos)))
+	proc_list[proc_count].start()
+	proc_count += 1
+	
+	if proc_count == nprocs or proc_count == num:
+	    proc_count = 0
+	    exitcodes = 0
+	    for proc in proc_list:
+		proc.join()
+		exitcodes += proc.exitcode
+		
+	    if exitcodes != 0:
+		dbif.close()
+		grass.fatal(_("Error while interpolation computation"))
+		
+	    # Empty process list
+	    proc_list = []
+
+    # Insert new interpolated maps in temporal database and dataset
+    for _map in gap_list:
+	_map.load()
+	_map.insert(dbif)
+	sp.register_map(_map, dbif)
+	
+    sp.update_from_registered_maps(dbif)
+    dbif.close()
+	
+###############################################################################
+
+def run_interp(inputs, dpos, output, outpos):
+    """Helper function to run r.series.interp in parallel"""
+    return grass.run_command("r.series.interp", input=inputs, datapos=dpos, output=output, samplingpos=outpos, overwrite=grass.overwrite(), quiet=True)
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()
+

temporal/t.rast.gapfill/test.t.rast.gapfill.sh

@@ -0,0 +1,32 @@
+# We need to set a specific region in the
+# @preprocess step of this test. 
+# 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 = 100"
+r.mapcalc --o expr="prec_2 = 300"
+r.mapcalc --o expr="prec_3 = 500"
+
+n1=`g.tempfile pid=1 -d` 
+
+cat > "${n1}" << EOF
+prec_1|2001-01-01|2001-02-01
+prec_2|2001-03-01|2001-04-01
+prec_3|2001-05-01|2001-06-01
+EOF
+
+t.create --v --o type=strds temporaltype=absolute output=precip_abs title="A test" descr="A test"
+t.register --v type=rast input=precip_abs file=${n1}
+
+# @test
+t.rast.gapfill input=precip_abs base="prec" nprocs=2
+t.info precip_abs
+
+t.info type=rast input=prec_6
+t.info type=rast input=prec_7
+
+# @postprocess
+
+t.unregister --v type=rast maps=prec_1,prec_2,prec_3,prec_6,prec_7
+t.remove --v type=strds input=precip_abs
+g.remove rast=prec_1,prec_2,prec_3,prec_6,prec_7

temporal/t.rast.mapcalc/t.rast.mapcalc.py

@@ -39,8 +39,12 @@
 #%option G_OPT_STRDS_OUTPUT
 #%end
 
-#%option G_OPT_R_BASE
+#%option
+#% key: base
+#% type: string
+#% description: The base name of the new created raster maps. This name will be extended with a numerical prefix
 #% required: yes
+#% multiple: no
 #%end
 
 #%option

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

@@ -26,7 +26,7 @@
 #%option
 #% key: method
 #% type: string
-#% description: Aggregate operation to be peformed on the raster maps
+#% description: Aggregate operation to be performed on the raster maps
 #% 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

temporal/t.rast3d.mapcalc/t.rast3d.mapcalc.py

@@ -39,10 +39,12 @@
 #%option G_OPT_STR3DS_OUTPUT
 #%end
 
-#%option G_OPT_R_BASE
+#%option
+#% key: base
+#% type: string
+#% description: The base name of the new created 3d raster maps. This name will be extended with a numerical prefix
 #% required: yes
-#% description: Name of the base raster3d map
-#% gisprompt: old,grid3,3d-raster
+#% multiple: no
 #%end
 
 #%option

temporal/t.topology/t.topology.py

@@ -69,7 +69,7 @@ def main():
     maps = sp.get_registered_maps_as_objects(where=where, order="start_time", dbif=None)
 
     if tmatrix:
-        matrix = sp.print_temporal_relation_matrix(maps)
+        sp.print_temporal_relationships(maps)
         return
 
     sp.base.print_info()