[GRASS-SVN] r50196 - in grass/trunk/temporal: . tr3.extract

svn_grass at osgeo.org svn_grass at osgeo.org
Sat Jan 14 20:32:43 EST 2012


Author: huhabla
Date: 2012-01-14 17:32:43 -0800 (Sat, 14 Jan 2012)
New Revision: 50196

Added:
   grass/trunk/temporal/tr3.extract/
   grass/trunk/temporal/tr3.extract/Makefile
   grass/trunk/temporal/tr3.extract/test.tr3.extract.sh
   grass/trunk/temporal/tr3.extract/tr3.extract.html
   grass/trunk/temporal/tr3.extract/tr3.extract.py
Log:
New space time raster3d module to extract subsets of raster3d maps


Added: grass/trunk/temporal/tr3.extract/Makefile
===================================================================
--- grass/trunk/temporal/tr3.extract/Makefile	                        (rev 0)
+++ grass/trunk/temporal/tr3.extract/Makefile	2012-01-15 01:32:43 UTC (rev 50196)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = tr3.extract
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)

Added: grass/trunk/temporal/tr3.extract/test.tr3.extract.sh
===================================================================
--- grass/trunk/temporal/tr3.extract/test.tr3.extract.sh	                        (rev 0)
+++ grass/trunk/temporal/tr3.extract/test.tr3.extract.sh	2012-01-15 01:32:43 UTC (rev 50196)
@@ -0,0 +1,28 @@
+# Test the extraction of a subset of a space time raster input
+
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster with r.mapcalc and create two space time raster inputs
+# with relative and absolute time
+# 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
+
+r3.mapcalc --o expr="prec_1 = rand(0, 550)"
+r3.mapcalc --o expr="prec_2 = rand(0, 450)"
+r3.mapcalc --o expr="prec_3 = rand(0, 320)"
+r3.mapcalc --o expr="prec_4 = rand(0, 510)"
+r3.mapcalc --o expr="prec_5 = rand(0, 300)"
+r3.mapcalc --o expr="prec_6 = rand(0, 650)"
+
+t.create --o type=str3ds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
+tr3.register -i input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="3 months"
+
+# The first @test
+# We create the space time raster inputs and register the raster maps with absolute time interval
+
+tr3.extract --o input=precip_abs1 output=precip_abs2 where="start_time > '2001-06-01'" expression=" if(precip_abs1 > 400, precip_abs1, null())" base=new_prec
+
+t.info type=str3ds input=precip_abs2
+
+t.remove type=rast3d input=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.remove type=str3ds input=precip_abs1,precip_abs2

Added: grass/trunk/temporal/tr3.extract/tr3.extract.html
===================================================================
Added: grass/trunk/temporal/tr3.extract/tr3.extract.py
===================================================================
--- grass/trunk/temporal/tr3.extract/tr3.extract.py	                        (rev 0)
+++ grass/trunk/temporal/tr3.extract/tr3.extract.py	2012-01-15 01:32:43 UTC (rev 50196)
@@ -0,0 +1,205 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:	tr3.extract
+# AUTHOR(S):	Soeren Gebbert
+#
+# PURPOSE:	Extract a subset of a space time raster3d dataset
+# 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: Extract a subset of a space time raster3d dataset
+#% keywords: spacetime raster3d dataset
+#% keywords: raster3d
+#% keywords: extract
+#%end
+
+#%option
+#% key: input
+#% type: string
+#% description: Name of a space time raster3d dataset
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: where
+#% type: string
+#% description: A where statement for selected listing without "WHERE" e.g: "start_time < '2001-01-01' AND end_time > '2001-01-01'"
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: expression
+#% type: string
+#% description: The r3.mapcalc expression assigned to all extracted raster3d maps
+#% required: no
+#% multiple: no
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% description: Name of the output space time raster3d dataset
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: base
+#% type: string
+#% description: Base name of the new created raster3d maps
+#% required: no
+#% multiple: no
+#%end
+
+#%flag
+#% key: n
+#% description: Register Null maps
+#%end
+
+import grass.script as grass
+import grass.temporal as tgis
+
+############################################################################
+
+def main():
+
+    # Get the options
+    input = options["input"]
+    output = options["output"]
+    where = options["where"]
+    expression = options["expression"]
+    base = options["base"]
+    register_null = flags["n"]
+
+    if expression and not base:
+        grass.fatal(_("Please specify %s=")%"base")
+
+    # Make sure the temporal database exists
+    tgis.create_temporal_database()
+    
+    mapset =  grass.gisenv()["MAPSET"]
+
+    if input.find("@") >= 0:
+        id = input
+    else:
+        id = input + "@" + mapset
+
+    sp = tgis.space_time_raster3d_dataset(id)
+    
+    if sp.is_in_db() == False:
+        grass.fatal(_("Space time raster3d dataset <%s> not found in temporal database") % (id))
+
+    dbif = tgis.sql_database_interface()
+    dbif.connect()
+
+    sp.select(dbif)
+
+    if output.find("@") >= 0:
+        out_id = output
+    else:
+        out_id = output + "@" + mapset
+
+    # The new space time raster3d dataset
+    new_sp = tgis.space_time_raster3d_dataset(out_id)
+    if new_sp.is_in_db():
+        if grass.overwrite() == True:
+            new_sp.delete(dbif)
+	    new_sp = tgis.space_time_raster3d_dataset(out_id)
+        else:
+            grass.fatal(_("Space time raster3d dataset <%s> is already in database, use overwrite flag to overwrite") % out_id)
+
+    temporal_type, semantic_type, title, description = sp.get_initial_values()
+    new_sp.set_initial_values(temporal_type, semantic_type, title, description)
+    new_sp.insert(dbif)
+
+    rows = sp.get_registered_maps("id", where, "start_time", dbif)
+
+    if rows:
+	num_rows = len(rows)
+	
+	grass.percent(0, num_rows, 1)
+	
+        count = 0
+        for row in rows:
+            count += 1
+	    
+	    grass.percent(count, num_rows, 1)
+
+            old_map = sp.get_new_map_instance(row["id"])
+            old_map.select(dbif)
+            
+            if expression:
+
+                map_name = "%s_%i" % (base, count)
+
+                expr = "%s = %s" % (map_name, expression.replace(sp.get_id(), row["id"]))
+                expr = expr.replace(sp.base.get_name(), row["id"])
+
+                map_id = map_name + "@" + mapset
+
+                new_map = sp.get_new_map_instance(map_id)
+
+                # Check if new map is in the temporal database
+                if new_map.is_in_db(dbif):
+                    if grass.overwrite() == True:
+                        # Remove the existing temporal database entry
+                        new_map.delete(dbif)
+                        new_map = sp.get_new_map_instance(map_id)
+                    else:
+                        grass.error(_("Raster3d map <%s> is already in temporal database, use overwrite flag to overwrite"))
+                        continue
+
+                grass.verbose(_("Apply r3.mapcalc expression: \"%s\"") % expr)
+
+                ret = grass.run_command("r3.mapcalc", expression=expr, overwrite=grass.overwrite(), quiet=True)
+
+                if ret != 0:
+                    grass.error(_("Error while r3.mapcalc computation, continue with next map"))
+                    break
+
+                # Read the raster3d map data
+                new_map.load()
+                
+                # In case of a null map continue, do not register null maps
+                print new_map.metadata.get_min(), new_map.metadata.get_max()
+                if new_map.metadata.get_min() == None and new_map.metadata.get_max() == None:
+                    print("Found null map")
+                    if not register_null:
+                        continue
+
+                # Set the time stamp
+                if old_map.is_time_absolute():
+                    start, end, tz = old_map.get_absolute_time()
+                    new_map.set_absolute_time(start, end, tz)
+                else:
+                    start, end = old_map.get_relative_time()
+                    new_map.set_relative_time(start, end)
+
+                # Insert map in temporal database
+                new_map.insert(dbif)
+
+                new_sp.register_map(new_map, dbif)
+            else:
+                new_sp.register_map(old_map, dbif)
+
+        # Update the spatio-temporal extent and the raster3d metadata table entries
+        new_sp.update_from_registered_maps(dbif)
+	
+	grass.percent(num_rows, num_rows, 1)
+        
+    dbif.close()
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()
+


Property changes on: grass/trunk/temporal/tr3.extract/tr3.extract.py
___________________________________________________________________
Added: svn:executable
   + *



More information about the grass-commit mailing list