#!/usr/bin/env python3 ############################################################################ # # 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()