Geografic-Information-System/scripts/v.what.strds/v.what.strds.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
############################################################################
#
# MODULE:       v.what.strds
# AUTHOR(S):    Luca delucchi
#
# PURPOSE:      Uploads space time raster dataset values at positions of vector points to the table
# COPYRIGHT:    (C) 2013 by the GRASS Development Team
#
#               This program is free software under the GNU General Public
#               License (version 2). Read the file COPYING that comes with GRASS
#               for details.
#
#############################################################################

#%module
#% description: Uploads space time raster dataset values at positions of vector points to the table.
#% keyword: vector
#% keyword: temporal
#% keyword: sampling
#% keyword: position
#% keyword: querying
#% keyword: attribute table
#% keyword: time
#%end

#%option G_OPT_V_INPUT
#%end

#%option G_OPT_STRDS_INPUTS
#% key: strds
#%end

#%option G_OPT_V_OUTPUT
#% required: no
#%end

#%option G_OPT_DB_WHERE
#%end

#%option G_OPT_T_WHERE
#% key: t_where
#%end

#%flag
#% key: u
#% label: Update attribute table of input vector map
#% description: Instead of creating a new vector map update the attribute table with value(s)
#%end

import grass.script as grass
from grass.exceptions import CalledModuleError

############################################################################


class Sample(object):

    def __init__(self, start=None, end=None, raster_names=None,
                 strds_name=None, granularity=None):
        self.start = start
        self.end = end
        if raster_names is not None:
            self.raster_names = raster_names
        else:
            self.raster_names = []
        self.strds_name = strds_name
        self.granu = granularity

    def __str__(self):
        return "Start: %s\nEnd: %s\nNames: %s\n" % (str(self.start),
                                                    str(self.end),
                                                    str(self.raster_names))

    def printDay(self, date='start'):
        output = ''
        if date == 'start':
            output = str(self.start).split(' ')[0].replace('-', '_')
        elif date == 'end':
            output = str(self.end).split(' ')[0].replace('-', '_')
        else:
            grass.fatal("The values accepted by printDay in Sample are:"
                        " 'start', 'end'")
        if self.granu:
            if self.granu.find('minute') != -1 or self.granu.find('second') != -1:
                output += '_' + str(self.start).split(' ')[1].replace(':', '_')
        return output

############################################################################


def main():
    # lazy imports
    import grass.temporal as tgis
    from grass.pygrass.utils import copy as gcopy
    from grass.pygrass.messages import Messenger
    from grass.pygrass.vector import Vector
    
    # Get the options
    input = options["input"]
    output = options["output"]
    strds = options["strds"]
    where = options["where"]
    tempwhere = options["t_where"]

    if output and flags['u']:
        grass.fatal(_("Cannot combine 'output' option and 'u' flag"))
    elif not output and not flags['u']:
        grass.fatal(_("'output' option or 'u' flag must be given"))
    elif not output and flags['u']:
        grass.warning(_("Attribute table of vector {name} will be updated...").format(name=input))

    if where == "" or where == " " or where == "\n":
        where = None

    overwrite = grass.overwrite()

    quiet = True

    if grass.verbosity() > 2:
        quiet = False

    # Check the number of sample strds and the number of columns
    strds_names = strds.split(",")

    # Make sure the temporal database exists
    tgis.init()
    # We need a database interface
    dbif = tgis.SQLDatabaseInterfaceConnection()
    dbif.connect()

    samples = []

    first_strds = tgis.open_old_stds(strds_names[0], "strds", dbif)
    # Single space time raster dataset
    if len(strds_names) == 1:
        granu = first_strds.get_granularity()
        rows = first_strds.get_registered_maps("name,mapset,start_time,end_time",
                                               tempwhere, "start_time",
                                               dbif)

        if not rows:
            dbif.close()
            grass.fatal(_("Space time raster dataset <%s> is empty") %
                        first_strds.get_id())
        for row in rows:
            start = row["start_time"]
            end = row["end_time"]
            raster_maps = [row["name"] + "@" + row["mapset"], ]

            s = Sample(start, end, raster_maps, first_strds.get_name(), granu)
            samples.append(s)
    else:
        # Multiple space time raster datasets
        for name in strds_names[1:]:
            dataset = tgis.open_old_stds(name, "strds", dbif)
            if dataset.get_temporal_type() != first_strds.get_temporal_type():
                grass.fatal(_("Temporal type of space time raster "
                              "datasets must be equal\n<%(a)s> of type "
                              "%(type_a)s do not match <%(b)s> of type "
                              "%(type_b)s" % {"a": first_strds.get_id(),
                                              "type_a": first_strds.get_temporal_type(),
                                              "b": dataset.get_id(),
                                              "type_b": dataset.get_temporal_type()}))

        mapmatrizes = tgis.sample_stds_by_stds_topology("strds", "strds",
                                                        strds_names,
                                                        strds_names[0],
                                                        False, None,
                                                        "equal", False,
                                                        False)
        #TODO check granularity for multiple STRDS
        for i in range(len(mapmatrizes[0])):
            isvalid = True
            mapname_list = []
            for mapmatrix in mapmatrizes:

                entry = mapmatrix[i]

                if entry["samples"]:
                    sample = entry["samples"][0]
                    name = sample.get_id()
                    if name is None:
                        isvalid = False
                        break
                    else:
                        mapname_list.append(name)

            if isvalid:
                entry = mapmatrizes[0][i]
                map = entry["granule"]

                start, end = map.get_temporal_extent_as_tuple()
                s = Sample(start, end, mapname_list, name)
                samples.append(s)

    # Get the layer and database connections of the input vector
    if output:
        gcopy(input, output, 'vector')
    else:
        output = input

    msgr = Messenger()
    perc_curr = 0
    perc_tot = len(samples)
    pymap = Vector(output)
    try:
        pymap.open('r')
    except:
        dbif.close()
        grass.fatal(_("Unable to create vector map <%s>" % output))

    if len(pymap.dblinks) == 0:
        try:
            pymap.close()
            grass.run_command("v.db.addtable", map=output)
        except CalledModuleError:
            dbif.close()
            grass.fatal(_("Unable to add table <%s> to vector map <%s>" % output))
    if pymap.is_open():
        pymap.close()

    for sample in samples:
        raster_names = sample.raster_names
        # Call v.what.rast for each raster map

        for name in raster_names:
            coltype = "DOUBLE PRECISION"
            # Get raster map type
            raster_map = tgis.RasterDataset(name)
            raster_map.load()
            if raster_map.metadata.get_datatype() == "CELL":
                coltype = "INT"
            day = sample.printDay()
            column_name = "%s_%s" % (sample.strds_name, day)
            column_string = "%s %s" % (column_name, coltype)
            column_string.replace('.', '_')
            try:
                grass.run_command("v.db.addcolumn", map=output,
                                  column=column_string,
                                  overwrite=overwrite)
            except CalledModuleError:
                dbif.close()
                grass.fatal(_("Unable to add column %s to vector map "
                              "<%s> ") % (column_string, output))
            try:
                grass.run_command("v.what.rast", map=output, raster=name,
                                  column=column_name, where=where,
                                  quiet=quiet)
            except CalledModuleError:
                dbif.close()
                grass.fatal(_("Unable to run v.what.rast for vector map"
                              " <%s> and raster map <%s>") %
                            (output, str(raster_names)))

        msgr.percent(perc_curr, perc_tot, 1)
        perc_curr += 1

    dbif.close()
if __name__ == "__main__":
    options, flags = grass.parser()
    main()