New module for spatio-temporal sampling of raster maps with vector points

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

temporal/t.remove/t.remove.py

@@ -51,6 +51,10 @@ import grass.temporal as tgis
 
 ############################################################################
 
+#
+# TODO: Add recursive and physical removement of maps and datasets
+#
+
 def main():
     
     # Get the options

temporal/tr.aggregate.ds/tr.aggregate.ds.py

@@ -175,7 +175,7 @@ def main():
 
     if not rows:
             dbif.close()
-            grass.fatal(_("Aggregation dataset <%s> is empty") % out_id)
+            grass.fatal(_("Aggregation dataset <%s> is empty") % id)
 
     count = 0
     for row in rows:

temporal/tv.what.rast/Makefile

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

temporal/tv.what.rast/test.tv.what.rast.sh

@@ -0,0 +1,52 @@
+# Test the temporal and spatial sampling of raster maps by vector point maps
+# 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.0"
+r.mapcalc --o expr="prec_2 = 200.0"
+r.mapcalc --o expr="prec_3 = 300"
+r.mapcalc --o expr="prec_4 = 400"
+r.mapcalc --o expr="prec_5 = 500.0"
+r.mapcalc --o expr="prec_6 = 600.0"
+
+v.random --o -z output=soil_1 n=5 zmin=0 zmax=100 column=heigh seed=1
+v.random --o -z output=soil_2 n=5 zmin=0 zmax=100 column=height seed=2
+v.random --o -z output=soil_3 n=5 zmin=0 zmax=100 column=height seed=3
+
+n1=`g.tempfile pid=1 -d` 
+
+cat > $n1 << EOF
+soil_1|2001-01-01|2001-04-01
+soil_2|2001-05-01|2001-07-01
+soil_3|2001-08-01|2001-12-01
+EOF
+
+t.create --o type=stvds temporaltype=absolute output=soil_abs1 title="A test" descr="A test"
+tv.register input=soil_abs1 file="$n1" start=file end=file
+
+t.create --o type=strds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
+tr.register -i input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-03-01 00:00:00" increment="1 months"
+
+# The @test
+
+tv.what.rast --v input=soil_abs1 strds=precip_abs1 sampling=start,during column=map_vals
+v.db.select map=soil_1
+v.db.select map=soil_2
+v.db.select map=soil_3
+
+tv.what.rast --v input=soil_abs1 strds=precip_abs1 sampling=start,during
+v.db.select map=soil_1
+v.db.select map=soil_2
+v.db.select map=soil_3
+
+# @postprocess
+t.remove type=vect input=soil_1,soil_2,soil_3
+t.remove type=stvds input=soil_abs1
+
+t.remove type=rast input=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.remove type=strds input=precip_abs1
+
+g.remove rast=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+g.remove vect=soil_1,soil_2,soil_3

temporal/tv.what.rast/tv.what.rast.py

@@ -0,0 +1,192 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:	tv.what.rast
+# AUTHOR(S):	Soeren Gebbert
+#
+# PURPOSE:	Uploads raster map values at spatial and temporal positions of vector points to the tables. The maps are registered in space time datasets
+# COPYRIGHT:	(C) 2011 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 raster map values at spatial and temporal positions of vector points to the tables. The maps are registered in space time datasets
+#% keywords: space time vector dataset
+#% keywords: space time raster dataset
+#% keywords: sampling
+#%end
+
+#%option
+#% key: input
+#% type: string
+#% description: Name of a space time vector dataset
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: strds
+#% type: string
+#% description: The space time raster dataset to use
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: column
+#% type: string
+#% description: Name of the vector column to be created and to store sampled raster values. The use of a column name forces tv.what.rast to sample only values from the first map found in an interval. Otherwise the raster map names are used as column names.
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: where
+#% type: string
+#% description: A where statement for attribute selection related to the input dataset without "WHERE" e.g: "cat < 200"
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: tempwhere
+#% type: string
+#% description: A where statement for temporal selection related to the input dataset without "WHERE" e.g: "start_time < '2001-01-01' AND end_time > '2001-01-01'"
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: sampling
+#% type: string
+#% description: The method to be used for temporal sampling
+#% required: no
+#% multiple: yes
+#% options: start,during,overlap,contain,equal
+#% answer: start
+#%end
+
+
+import grass.script as grass
+import grass.temporal as tgis
+import grass.script.raster as raster
+
+############################################################################
+
+def main():
+
+    # Get the options
+    input = options["input"]
+    strds = options["strds"]
+    where = options["where"]
+    column = options["column"]
+    tempwhere = options["tempwhere"]
+    sampling = options["sampling"]
+
+    if where == "" or where == " " or where == "\n":
+        where = None
+
+    # Make sure the temporal database exists
+    tgis.create_temporal_database()
+    # We need a database interface
+    dbif = tgis.sql_database_interface()
+    dbif.connect()
+   
+    mapset =  grass.gisenv()["MAPSET"]
+
+    if input.find("@") >= 0:
+        id = input
+    else:
+        id = input + "@" + mapset
+
+    sp = tgis.space_time_vector_dataset(id)
+    
+    if sp.is_in_db() == False:
+        dbif.close()
+        grass.fatal(_("Dataset <%s> not found in temporal database") % (id))
+
+    sp.select(dbif)
+
+    if strds.find("@") >= 0:
+        strds_id = strds
+    else:
+        strds_id = strds + "@" + mapset
+
+    strds_sp = tgis.space_time_raster_dataset(strds_id)
+    
+    if strds_sp.is_in_db() == False:
+        dbif.close()
+        grass.fatal(_("Dataset <%s> not found in temporal database") % (id))
+
+    strds_sp.select(dbif)
+
+    if strds_sp.get_temporal_type() != sp.get_temporal_type():
+        dbif.close()
+        grass.fatal(_("Input and aggregation dataset must have the same temporal type"))
+
+    # Check if intervals are present
+    if sp.get_temporal_type() == "absolute":
+        map_time = strds_sp.absolute_time.get_map_time()
+    else:
+        map_time = strds_sp.relative_time.get_map_time()
+    
+    if map_time != "interval":
+        dbif.close()
+        grass.fatal(_("All registered maps of the space time vector dataset must have time intervals"))
+
+    rows = sp.get_registered_maps("id,name,start_time,end_time", tempwhere, "start_time", dbif)
+
+    if not rows:
+        dbif.close()
+        grass.fatal(_("Space time vector dataset <%s> is empty") % sg.get_id())
+
+    # Sample the raster dataset with the vector dataset and run v.what.rast
+    for row in rows:
+        start = row["start_time"]
+        end = row["end_time"]
+        vectmap = row["id"]
+
+        raster_maps = tgis.collect_map_names(strds_sp, dbif, start, end, sampling)
+
+        if raster_maps:
+            # Collect the names of the raster maps
+            for rastermap in raster_maps:
+                
+                if column:
+                    col_name = column
+                else:
+                    # Create a new column with the SQL compliant name of the sampled raster map
+                    col_name = rastermap.split("@")[0].replace(".", "_")
+
+                coltype = "DOUBLE PRECISION"
+                # Get raster type
+                rasterinfo = raster.raster_info(rastermap)
+                if rasterinfo["datatype"] == "CELL":
+                    coltype = "INT"
+
+                ret = grass.run_command("v.db.addcolumn", map=vectmap, column="%s %s" % (col_name, coltype), overwrite=grass.overwrite())
+                if ret != 0:
+                    dbif.close()
+                    grass.fatal(_("Unable to add column %s to vector map <%s>")%(col_name, vectmap))
+
+                # Call v.what.rast
+                ret = grass.run_command("v.what.rast", map=vectmap, raster=rastermap, column=col_name, where=where)
+                if ret != 0:
+                    dbif.close()
+                    grass.fatal(_("Unable to run v.what.rast for vector map <%s> and raster map <%s>")%(vectmap, rastermap))
+
+                # Use the first map in case a column names was provided
+                if column:
+                    break
+
+    dbif.close()
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()
+