temporal modules: New module that utilizes r.what for STRDS sampling

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

temporal/t.rast.what/Makefile

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

temporal/t.rast.what/t.rast.what.html

@@ -0,0 +1,102 @@
+<h2>DESCRIPTION</h2>
+
+<em>t.rast.what</em> is designed to sample space time raster datasets
+at specific coordinates using <a href="r.what.html">r.what</a> internally. 
+The output of <a href="r.what.html">r.what</a> 
+is transformed in different output layouts.
+The output layouts can be specified using the option <em>layout</em>.
+Please have a look at the output layout demonstration in the example.
+<p>
+This module is designed to run several instances of r.what to sample
+parts of a space time raster dataset in parallel. Several intermediate
+text files will be produces that are merged into a single file at the
+and of the processing.
+<p>
+Coordinates can be provided using the <em>coordinates</em> option
+or as vector map using the <em>points</em> option.
+<p>
+An output file must be provided with the <em>output</em> option.
+    
+<h2>EXAMPLE</h2>
+
+In this example we sample a space time raster dataset that contains
+4 raster map layer with 3 random vector points. The generated output
+of all layout options are demonstrated. 
+
+<div class="code"><pre>
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+# Generate data
+r.mapcalc expr="a_1 = 1" -s
+r.mapcalc expr="a_2 = 2" -s
+r.mapcalc expr="a_3 = 3" -s
+r.mapcalc expr="a_4 = 4" -s
+
+v.random output=points n=3
+
+t.create type=strds output=A title="A test" descr="A test"
+
+t.register -i type=raster input=A maps=a_1,a_2,a_3,a_4 \
+    start='1990-01-01' increment="1 month"
+
+# Now we create 3 outputs of the same data using different layouts
+
+# Output in "row" order
+t.rast.what strds=A points=points output=result.txt layout=row -n
+
+cat result.txt
+
+x|y|start|end|value
+107.9761951|14.4780173|1990-01-01 00:00:00|1990-02-01 00:00:00|1
+107.9761951|14.4780173|1990-02-01 00:00:00|1990-03-01 00:00:00|2
+107.9761951|14.4780173|1990-03-01 00:00:00|1990-04-01 00:00:00|3
+107.9761951|14.4780173|1990-04-01 00:00:00|1990-05-01 00:00:00|4
+48.1672585|75.5186798|1990-01-01 00:00:00|1990-02-01 00:00:00|1
+48.1672585|75.5186798|1990-02-01 00:00:00|1990-03-01 00:00:00|2
+48.1672585|75.5186798|1990-03-01 00:00:00|1990-04-01 00:00:00|3
+48.1672585|75.5186798|1990-04-01 00:00:00|1990-05-01 00:00:00|4
+114.3944372|36.2390227|1990-01-01 00:00:00|1990-02-01 00:00:00|1
+114.3944372|36.2390227|1990-02-01 00:00:00|1990-03-01 00:00:00|2
+114.3944372|36.2390227|1990-03-01 00:00:00|1990-04-01 00:00:00|3
+114.3944372|36.2390227|1990-04-01 00:00:00|1990-05-01 00:00:00|4
+
+# Output in column order
+t.rast.what strds=A points=points output=result.txt layout=col -n
+
+cat result.txt
+
+star|end|107.9761951;14.4780173|107.9761951;14.4780173|107.9761951;14.4780173
+1990-01-01 00:00:00|1990-02-01 00:00:00|1|1|1
+1990-02-01 00:00:00|1990-03-01 00:00:00|2|2|2
+1990-03-01 00:00:00|1990-04-01 00:00:00|3|3|3
+1990-04-01 00:00:00|1990-05-01 00:00:00|4|4|4
+
+# Output in one time series per row order
+t.rast.what strds=A points=points output=result.txt layout=timerow -n
+
+cat result.txt
+
+x|y|1990-01-01 00:00:00;1990-02-01 00:00:00|1990-02-01 00:00:00;1990-03-01 00:00:00|1990-03-01 00:00:00;1990-04-01 00:00:00|1990-04-01 00:00:00;1990-05-01 00:00:00
+107.976195095579|14.4780172666265|1|2|3|4
+48.1672585088805|75.5186797967386|1|2|3|4
+114.394437197952|36.2390227089844|1|2|3|4
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.what.html">r.what</a> ,
+<a href="r.neighbors.html">r.neighbors</a>,
+<a href="t.rast.aggregate.ds.html">t.rast.aggregate.ds</a>,
+<a href="t.rast.extract.html">t.rast.extract</a>,
+<a href="t.info.html">t.info</a>,
+<a href="g.region.html">g.region</a>,
+<a href="r.mask.html">r.mask</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+S&ouml;ren Gebbert, Th&uuml;nen Institute of Climate-Smart Agriculture
+
+<p><i>Last changed: $Date: 2014-12-03 20:39:26 +0100 (Mi, 03. Dez 2014) $</i>

temporal/t.rast.what/t.rast.what.py

@@ -0,0 +1,449 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:       t.rast.what
+# AUTHOR(S):    Soeren Gebbert
+#
+# PURPOSE:      Sample a space time raster dataset at specific coordinates 
+#               and write the output to stdout using different layouts
+#               
+# 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: Sample a space time raster dataset at specific coordinates and write the output to stdout using different layouts
+#% keyword: temporal
+#% keyword: raster
+#% keyword: sampling
+#%end
+
+#%option G_OPT_V_INPUT
+#% key: points
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#% key: strds
+#%end
+
+#%option G_OPT_F_OUTPUT
+#%end
+
+#%option G_OPT_T_WHERE
+#%end
+
+#%option G_OPT_M_COORDS
+#%end
+
+#%option G_OPT_M_NULL_VALUE
+#%end
+
+#%option G_OPT_F_SEP
+#%end
+
+#%option
+#% key: order
+#% type: string
+#% description: Sort the maps by category
+#% required: no
+#% multiple: yes
+#% options: id, name, creator, mapset, creation_time, modification_time, start_time, end_time, north, south, west, east, min, max
+#% answer: start_time
+#%end
+
+#%option
+#% key: layout
+#% type: string
+#% description: The layout of the ouput. One point per row (row), one point per column (col), all timsteps in one row (timerow)
+#% required: no
+#% multiple: no
+#% options: row, col, timerow
+#% answer: row
+#%end
+
+#%option
+#% key: nprocs
+#% type: integer
+#% description: Number of r.neighbor processes to run in parallel
+#% required: no
+#% multiple: no
+#% answer: 1
+#%end
+
+#%flag
+#% key: n
+#% description: Output header row
+#%end
+
+#%flag
+#% key: f
+#% description: Show the category labels of the grid cell(s)
+#%end
+
+#%flag
+#% key: r
+#% description: Output color values as RRR:GGG:BBB
+#%end
+
+#%flag
+#% key: i
+#% description: Output integer category values, not cell values
+#%end
+
+import copy
+import grass.script as grass
+import grass.temporal as tgis
+import grass.pygrass.modules as pymod
+
+
+############################################################################
+
+def main(options, flags):
+
+    # Get the options
+    points = options["points"]
+    strds = options["strds"]
+    output = options["output"]
+    where = options["where"]
+    order = options["order"]
+    layout = options["layout"]
+    coordinates = options["coordinates"]
+    null_value = options["null_value"]
+    separator = options["separator"]
+    
+    nprocs = int(options["nprocs"])
+    write_header = flags["n"]
+    output_cat_label = flags["f"]
+    output_color = flags["r"]
+    output_cat = flags["i"]
+    
+    overwrite = grass.overwrite()
+    
+    if coordinates and points:
+        grass.error(_("Options coordinates and points are mutually exclusive"))
+
+    # Make sure the temporal database exists
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+
+    sp = tgis.open_old_stds(strds, "strds", dbif)
+    maps = sp.get_registered_maps_as_objects(where=where, order=order, 
+                                             dbif=dbif)
+    dbif.close()
+
+    if not maps:
+        grass.warning(_("Space time raster dataset <%s> is empty") % sp.get_id())
+        return
+
+    # Setup separator
+    if separator == "pipe":
+        separator = "|"
+    if separator == "comma":
+        separator = ","
+    if separator == "space":
+        separator = " "
+    if separator == "tab":
+        separator = "\t"
+    if separator == "newline":
+        separator = "\n"
+
+    # Setup flags
+    flags = ""
+    if output_cat_label is True:
+        flags += "f"
+    if output_color is True:
+        flags += "r"
+    if output_cat is True:
+        flags += "f"
+
+    # Configure the r.what module
+    if points:
+        r_what = pymod.Module("r.what", map="dummy",
+                              output="dummy", run_=False,
+                              separator=separator, points=points,
+                              overwrite=overwrite, flags=flags,
+                              quiet=True)
+    elif coordinates:
+        r_what = pymod.Module("r.what", map="dummy",
+                              output="dummy", run_=False,
+                              separator=separator, 
+                              coordinates=coordinates,
+                              overwrite=overwrite, flags=flags,
+                              quiet=True)
+    else:
+        grass.error(_("Please specify points or coordinates"))
+
+    # The module queue for parallel execution
+    process_queue = pymod.ParallelModuleQueue(int(nprocs))
+    num_maps = len(maps)
+    
+    # 400 Maps is the absolute maximum in r.what
+    # We need to determie the number of maps that can be processed
+    # in parallel
+
+    # First estimate the number of maps per process. We use 400 maps
+    # simultaniously as maximum for a single process
+
+    num_loops = int(num_maps / (400 * nprocs))
+    remaining_maps = num_maps % (400 * nprocs)
+
+    if num_loops == 0:
+        num_loops = 1
+        remaining_maps = 0
+
+    # Compute the number of maps for each process
+    maps_per_loop = int((num_maps - remaining_maps) / num_loops)
+    maps_per_process = int(maps_per_loop / nprocs)
+    remaining_maps_per_loop = maps_per_loop % nprocs
+
+    # We put the output files in an ordered list
+    output_files = []
+    output_time_list = []
+
+    count = 0
+    for loop in range(num_loops):
+        file_name = "out_%i"%(loop)
+        count = process_loop(nprocs, maps, file_name, count, maps_per_process, 
+                             remaining_maps_per_loop, output_files, 
+                             output_time_list, r_what, process_queue)
+    
+    process_queue.wait()
+    
+    print "Remaining maps", remaining_maps, count, len(maps)
+    if remaining_maps > 0:
+        # Use a single process if less then 400 maps
+        if remaining_maps <= 100:
+            print "Single remain process"
+            mod = copy.deepcopy(r_what)
+            mod(map=map_names, output=file_name)
+            process_queue.put(mod)
+        else:
+            maps_per_process = int(remaining_maps / nprocs)
+            remaining_maps_per_loop = remaining_maps % nprocs
+            
+            file_name = "out_remain"
+            process_loop(nprocs, maps, file_name, count, maps_per_process, 
+                         remaining_maps_per_loop, output_files, 
+                         output_time_list, r_what, process_queue)
+
+    # Wait for unfinished processes
+    process_queue.wait()
+    
+    # Out the output files in the correct order together
+    if layout == "row":
+        one_point_per_row_output(separator, output_files, output_time_list,
+                                 output, write_header)
+    elif layout == "col":
+        one_point_per_col_output(separator, output_files, output_time_list,
+                                 output, write_header)
+    else:
+        one_point_per_timerow_output(separator, output_files, output_time_list,
+                                     output, write_header)
+
+############################################################################
+
+def one_point_per_row_output(separator, output_files, output_time_list,
+                             output, write_header):
+    """Write one point per row
+       output is of type: x,y,start,end,value
+    """
+    # open the output file for writing
+    out_file = open(output, "w")
+    if write_header is True:
+        out_file.write("x%(sep)sy%(sep)sstart%(sep)send%(sep)svalue\n"\
+                       %({"sep":separator}))
+    
+    for count in range(len(output_files)):
+        file_name = output_files[count]
+        map_list = output_time_list[count]
+        print "Transform", file_name
+        in_file = open(file_name, "r")
+        for line in in_file:
+            line = line.split(separator)
+            x = line[0]
+            y = line[1]
+            # We ignore the site name
+            values = line[3:]
+            for i in range(len(values)):
+                start, end = map_list[i].get_temporal_extent_as_tuple()
+                coor_string = "%(x)10.7f%(sep)s%(y)10.7f%(sep)s"\
+                               %({"x":float(x),"y":float(y),"sep":separator})
+                time_string = "%(start)s%(sep)s%(end)s%(sep)s%(val)s\n"\
+                               %({"start":str(start), "end":str(end),
+                                  "val":(values[i].strip()),"sep":separator})
+
+                out_file.write(coor_string + time_string)
+        
+        in_file.close()
+    out_file.close()
+        
+############################################################################
+
+def one_point_per_col_output(separator, output_files, output_time_list,
+                             output, write_header):
+    """Write one point per col
+       output is of type: 
+       start,end,point_1 value,point_2 value,...,point_n value
+       
+       Each row represents a single raster map, hence a single time stamp
+    """
+    # open the output file for writing
+    out_file = open(output, "w")
+        
+    first = True
+    for count in range(len(output_files)):
+        file_name = output_files[count]
+        map_list = output_time_list[count]
+        print "Transform", file_name
+        in_file = open(file_name, "r")
+        lines = in_file.readlines()
+        
+        matrix = []
+        for line in lines:
+            matrix.append(line.split(separator))
+            
+        num_cols = len(matrix[0])
+        
+        if first is True:
+            if write_header is True:
+                out_file.write("star%(sep)send"%({"sep":separator}))
+                for line in lines:
+                    x = matrix[0][0]
+                    y = matrix[0][1]
+                    out_file.write("%(sep)s%(x)10.7f;%(y)10.7f"\
+                                   %({"sep":separator,
+                                      "x":float(x), 
+                                      "y":float(y)}))
+                out_file.write("\n")
+
+        first = False
+
+        for col in xrange(num_cols - 3):
+            start, end = output_time_list[count][col].get_temporal_extent_as_tuple()
+            time_string = "%(start)s%(sep)s%(end)s"\
+                               %({"start":str(start), "end":str(end),
+                                  "sep":separator})
+            out_file.write(time_string)
+            for row in xrange(len(matrix)):
+                value = matrix[row][col + 3]
+                out_file.write("%(sep)s%(value)s"\
+                                   %({"sep":separator,
+                                      "value":value.strip()}))
+            out_file.write("\n")
+
+        in_file.close()
+    out_file.close()
+
+############################################################################
+
+def one_point_per_timerow_output(separator, output_files, output_time_list,
+                             output, write_header):
+    """Use the original layout of the r.waht output and print instead of 
+       the raster names, the time stamps as header
+       
+       One point per line for all time stamps:
+        x|y|1991-01-01 00:00:00;1991-01-02 00:00:00|1991-01-02 00:00:00;1991-01-03 00:00:00|1991-01-03 00:00:00;1991-01-04 00:00:00|1991-01-04 00:00:00;1991-01-05 00:00:00
+        3730731.49590371|5642483.51236521|6|8|7|7
+        3581249.04638104|5634411.97526282|5|8|7|7
+    """
+    out_file = open(output, "w")
+
+    matrix = []
+    header = ""
+
+    first = True
+    for count in range(len(output_files)):
+        file_name = output_files[count]
+        map_list = output_time_list[count]
+        print "Transform orig", file_name
+        in_file = open(file_name, "r")
+
+        if write_header:
+            if first is True:
+                header = "x%(sep)sy"%({"sep":separator})
+            for map in map_list:
+                start, end = map.get_temporal_extent_as_tuple()
+                time_string = "%(sep)s%(start)s;%(end)s"\
+                              %({"start":str(start), "end":str(end),
+                                 "sep":separator})
+                header += time_string
+
+        lines = in_file.readlines()
+
+        for i in xrange(len(lines)):
+            cols = lines[i].split(separator)
+
+            if first is True:
+                matrix.append(cols[:2])
+
+            matrix[i] = matrix[i] + cols[3:]
+
+        first = False
+
+        in_file.close()
+
+    out_file.write(header + "\n")
+
+    for row in matrix:
+        first = True
+        for col in row:
+            value = col.strip()
+            
+            if first is False:
+                out_file.write("%s"%(separator))
+            out_file.write(value)
+            
+            first = False
+
+        out_file.write("\n")
+
+    out_file.close()
+ 
+############################################################################
+
+def process_loop(nprocs, maps, file_name, count, maps_per_process, 
+                 remaining_maps_per_loop, output_files, 
+                 output_time_list, r_what, process_queue):
+    """Call r.what in parallel subprocesses"""
+    first = True
+    for process in range(nprocs):
+        num = maps_per_process
+        # Add the remaining maps to the first process
+        if first is True:
+            num += remaining_maps_per_loop
+        first = False
+
+        # Temporary output file
+        final_file_name = file_name + "_%i"%(process)
+        output_files.append(final_file_name)
+        
+        map_names = []
+        map_list = []
+        for i in range(num):
+            map = maps[count]
+            map_names.append(map.get_id())
+            map_list.append(map)
+            count += 1
+
+        output_time_list.append(map_list)
+
+        print "Process", process, "Maps", len(map_names)
+        mod = copy.deepcopy(r_what)
+        mod(map=map_names, output=final_file_name)
+        #print(mod.get_bash())
+        process_queue.put(mod)
+
+    return count
+    
+############################################################################
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main(options, flags)