Geografic-Information-System/temporal/t.vect.observe.strds/t.vect.observe.strds.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
############################################################################
#
# MODULE:       t.vect.observe.strds
# AUTHOR(S):    Soeren Gebbert
#
# PURPOSE:      Observe specific locations in a space time raster dataset over a period of time using vector points
# COPYRIGHT:    (C) 2011-2017 by the GRASS Development Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#############################################################################

# %module
# % description: Observes specific locations in a space time raster dataset over a period of time using vector points.
# % keyword: temporal
# % keyword: sampling
# % keyword: vector
# % keyword: time
# %end

# %option G_OPT_V_INPUT
# %end

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

# %option G_OPT_STVDS_OUTPUT
# %end

# %option G_OPT_V_OUTPUT
# % key: vector_output
# % description: Name of the new created vector map that stores the sampled values in different layers
# %end

# %option
# % key: columns
# % type: string
# % 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

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

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


class Sample(object):
    def __init__(self, start=None, end=None, raster_names=None):
        self.start = start
        self.end = end
        if raster_names is not 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():
    # lazy imports
    import grass.temporal as tgis

    # Get the options
    input = options["input"]
    output = options["output"]
    vector_output = options["vector_output"]
    strds = options["strds"]
    where = options["where"]
    columns = options["columns"]

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

    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
    dbif = tgis.SQLDatabaseInterfaceConnection()
    dbif.connect()

    mapset = grass.gisenv()["MAPSET"]

    out_sp = tgis.check_new_stds(output, "stvds", dbif, overwrite)

    samples = []

    first_strds = tgis.open_old_stds(strds_names[0], "strds", dbif)

    # Single space time raster dataset
    if len(strds_names) == 1:
        rows = first_strds.get_registered_maps(
            columns="name,mapset,start_time,end_time", order="start_time", dbif=dbif
        )

        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"],
            ]

            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_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,
        )

        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)
                samples.append(s)

    num_samples = len(samples)

    # Get the layer and database connections of the input vector
    vector_db = grass.vector.vector_db(input)

    # We copy the vector table and create the new layers
    if vector_db:
        # Use the first layer to copy the categories from
        layers = "1,"
    else:
        layers = ""
    first = True
    for layer in range(num_samples):
        layer += 1
        # Skip existing layer
        if vector_db and layer in vector_db and vector_db[layer]["layer"] == layer:
            continue
        if first:
            layers += "%i" % (layer)
            first = False
        else:
            layers += ",%i" % (layer)

    vectmap = vector_output

    # We create a new vector map using the categories of the original map
    try:
        grass.run_command(
            "v.category",
            input=input,
            layer=layers,
            output=vectmap,
            option="transfer",
            overwrite=overwrite,
        )
    except CalledModuleError:
        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_stds(
        output,
        "stvds",
        first_strds.get_temporal_type(),
        title,
        description,
        first_strds.get_semantic_type(),
        dbif,
        overwrite,
    )

    dummy = out_sp.get_new_map_instance(None)

    # Sample the space time raster dataset with the vector
    # map at specific layer with v.what.rast
    count = 1
    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"]:
            try:
                grass.run_command(
                    "v.db.addcolumn",
                    map=vectmap,
                    layer=count,
                    column=columns_string,
                    overwrite=overwrite,
                )
            except CalledModuleError:
                dbif.close()
                grass.fatal(
                    _("Unable to add column %s to vector map <%s> " "with layer %i")
                    % (columns_string, vectmap, count)
                )
        else:
            # Try to add a new table
            grass.message("Add table to layer %i" % (count))
            try:
                grass.run_command(
                    "v.db.addtable",
                    map=vectmap,
                    layer=count,
                    columns=columns_string,
                    overwrite=overwrite,
                )
            except CalledModuleError:
                dbif.close()
                grass.fatal(
                    _("Unable to add table to vector map " "<%s> with layer %i")
                    % (vectmap, count)
                )

        # Call v.what.rast for each raster map
        for name, column in zip(raster_names, column_names):
            try:
                grass.run_command(
                    "v.what.rast",
                    map=vectmap,
                    layer=count,
                    raster=name,
                    column=column,
                    where=where,
                )
            except CalledModuleError:
                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, first_strds.get_relative_time_unit())

        if vect.is_in_db(dbif):
            vect.update_all(dbif)
        else:
            vect.insert(dbif)

        out_sp.register_map(vect, dbif)
        count += 1

    out_sp.update_from_registered_maps(dbif)
    dbif.close()


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