[GRASS-SVN] r55471 - in grass/trunk/temporal: . t.rast.aggregate t.rast.neighbors
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Mar 20 11:11:45 PDT 2013
Author: huhabla
Date: 2013-03-20 11:11:45 -0700 (Wed, 20 Mar 2013)
New Revision: 55471
Added:
grass/trunk/temporal/t.rast.neighbors/
grass/trunk/temporal/t.rast.neighbors/Makefile
grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.html
grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.py
grass/trunk/temporal/t.rast.neighbors/test.t.rast.neighbors.sh
Modified:
grass/trunk/temporal/Makefile
grass/trunk/temporal/t.rast.aggregate/t.rast.aggregate.html
Log:
Added new module to perform neighborhood computations on STRDS
Modified: grass/trunk/temporal/Makefile
===================================================================
--- grass/trunk/temporal/Makefile 2013-03-20 17:53:11 UTC (rev 55470)
+++ grass/trunk/temporal/Makefile 2013-03-20 18:11:45 UTC (rev 55471)
@@ -17,6 +17,7 @@
t.rast.univar \
t.rast.list \
t.rast.mapcalc \
+ t.rast.neighbors \
t.rast.series \
t.rast.export \
t.rast.out.vtk \
Modified: grass/trunk/temporal/t.rast.aggregate/t.rast.aggregate.html
===================================================================
--- grass/trunk/temporal/t.rast.aggregate/t.rast.aggregate.html 2013-03-20 17:53:11 UTC (rev 55470)
+++ grass/trunk/temporal/t.rast.aggregate/t.rast.aggregate.html 2013-03-20 18:11:45 UTC (rev 55471)
@@ -59,7 +59,7 @@
description="Test dataset with daily precipitation"
t.register -i type=rast input=precipitation_daily \
- file=map_list.txt start=20-08-2012 increment="1 days"
+ file=map_list.txt start="2012-08-20" increment="1 days"
t.info type=strds input=precipitation_daily
Added: grass/trunk/temporal/t.rast.neighbors/Makefile
===================================================================
--- grass/trunk/temporal/t.rast.neighbors/Makefile (rev 0)
+++ grass/trunk/temporal/t.rast.neighbors/Makefile 2013-03-20 18:11:45 UTC (rev 55471)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast.neighbors
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)
Added: grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.html
===================================================================
--- grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.html (rev 0)
+++ grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.html 2013-03-20 18:11:45 UTC (rev 55471)
@@ -0,0 +1,178 @@
+<h2>DESCRIPTION</h2>
+
+<em>t.rast.neighbors</em> performs <a href="r.neighbors.html">r.neighbors</a> computations
+on the maps of a space time raster dataset (STRDS). This module supports a subset of options
+that are available in <a href="r.neighbors.html">r.neighbors</a>.
+The size of the neighborhood and the aggregation method can be chosen.
+<p>
+The user must provide an input and an output space time raster dataset and the basename
+of the resulting raster maps. The resulting STRDS will have
+the same temporal resolution as the input dataset.
+All maps will be processed using the current region settings.
+<p>
+The user can select a subset of the input space time raster dataset for processing
+using a SQL WHERE statement. The number of CPU's to be used for parallel processing can be specified
+with the <em>nprocs</em> option, to speedup the computation on multicore system.
+
+
+
+<h2>EXAMPLE</h2>
+
+In this example we create 7 raster maps that will be registered in a single space time
+raster dataset named <em>precipitation_daily</em> using a daily temporal granularity.
+The names of the raster maps are stored in a text file that is used for raster map registration.
+<p>
+A neighborhood average smoothing with a window size of 5 cells will be performed on
+the space time raster dataset <em>precipitation_daily</em> resulting in the output space time raster dataset
+<em>precipitation_daily_smooth</em>. The base name of the new generated raster maps is <em>prec_smooth</em>.
+The temporal resolution of the resulting STRDS is identical to the input STRDS.
+
+
+<div class="code"><pre>
+
+MAPS="map_1 map_2 map_3 map_4 map_5 map_6 map_7"
+
+count=1
+for map in ${MAPS} ; do
+ export GRASS_RND_SEED=${count}
+ r.mapcalc --o expr="${map} = rand(0, 550)"
+ echo ${map} >> map_list.txt
+ count=$((count+1))
+done
+
+t.create --o type=strds temporaltype=absolute \
+ output=precipitation_daily \
+ title="Daily precipitation" \
+ description="Test dataset with daily precipitation"
+
+t.register -i --o type=rast input=precipitation_daily \
+ file=map_list.txt start="2012-08-20" increment="1 days"
+
+t.info type=strds input=precipitation_daily
+
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ | |
+ +-------------------- Basic information -------------------------------------+
+ | Id: ........................ precipitation_daily at PERMANENT
+ | Name: ...................... precipitation_daily
+ | Mapset: .................... PERMANENT
+ | Creator: ................... soeren
+ | Creation time: ............. 2013-03-20 18:47:25.867658
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2012-08-20 00:00:00
+ | End time:................... 2012-08-27 00:00:00
+ | Granularity:................ 1 days
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | Number of registered maps:.. 7
+ | Title:
+ | Daily precipitation
+ | Description:
+ | Test dataset with daily precipitation
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 0.0
+ | Minimum value max:.......... 8.0
+ | Maximum value min:.......... 533.0
+ | Maximum value max:.......... 549.0
+ | Raster register table:...... precipitation_daily_PERMANENT_raster_register
+ +----------------------------------------------------------------------------+
+
+
+t.rast.list input=precipitation_daily columns=name,start_time,min,max
+
+map_1 2012-08-20 00:00:00 2.0 549.0
+map_2 2012-08-21 00:00:00 4.0 549.0
+map_3 2012-08-22 00:00:00 8.0 548.0
+map_4 2012-08-23 00:00:00 3.0 549.0
+map_5 2012-08-24 00:00:00 0.0 533.0
+map_6 2012-08-25 00:00:00 2.0 543.0
+map_7 2012-08-26 00:00:00 1.0 539.0
+
+
+t.rast.neighbors --o input=precipitation_daily \
+ output=precipitation_daily_smooth base=prec_smooth \
+ size=5 method=average nprocs=3
+
+t.info type=strds input=precipitation_daily_smooth
+
+ +-------------------- Space Time Raster Dataset -----------------------------+
+ | |
+ +-------------------- Basic information -------------------------------------+
+ | Id: ........................ precipitation_daily_smooth at PERMANENT
+ | Name: ...................... precipitation_daily_smooth
+ | Mapset: .................... PERMANENT
+ | Creator: ................... soeren
+ | Creation time: ............. 2013-03-20 18:47:28.697990
+ | Temporal type: ............. absolute
+ | Semantic type:.............. mean
+ +-------------------- Absolute time -----------------------------------------+
+ | Start time:................. 2012-08-20 00:00:00
+ | End time:................... 2012-08-27 00:00:00
+ | Granularity:................ 1 days
+ | Temporal type of maps:...... interval
+ +-------------------- Spatial extent ----------------------------------------+
+ | North:...................... 80.0
+ | South:...................... 0.0
+ | East:.. .................... 120.0
+ | West:....................... 0.0
+ | Top:........................ 0.0
+ | Bottom:..................... 0.0
+ +-------------------- Metadata information ----------------------------------+
+ | Number of registered maps:.. 7
+ | Title:
+ | Daily precipitation
+ | Description:
+ | Test dataset with daily precipitation
+ | North-South resolution min:. 10.0
+ | North-South resolution max:. 10.0
+ | East-west resolution min:... 10.0
+ | East-west resolution max:... 10.0
+ | Minimum value min:.......... 140.166667
+ | Minimum value max:.......... 215.166667
+ | Maximum value min:.......... 293.3
+ | Maximum value max:.......... 373.388889
+ | Raster register table:...... precipitation_daily_smooth_PERMANENT_raster_register
+ +----------------------------------------------------------------------------+
+
+
+t.rast.list input=precipitation_daily_smooth columns=name,start_time,min,max
+
+prec_smooth_1 2012-08-20 00:00:00 172.388889 360.566667
+prec_smooth_2 2012-08-21 00:00:00 215.166667 322.833333
+prec_smooth_3 2012-08-22 00:00:00 140.166667 373.388889
+prec_smooth_4 2012-08-23 00:00:00 184.3 356.583333
+prec_smooth_5 2012-08-24 00:00:00 163.966667 293.3
+prec_smooth_6 2012-08-25 00:00:00 204.98 324.611111
+prec_smooth_7 2012-08-26 00:00:00 154.25 340.366667
+
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.neighbors.html">r.neighbors</a><br>
+<a href="t.rast.aggregate.ds.html">t.rast.aggregate.ds</a><br>
+<a href="t.rast.extract.html">t.rast.extract</a><br>
+<a href="t.info.html">t.info</a><br>
+<a href="g.region.html">g.region</a><br>
+<a href="r.mask.html">r.mask</a><br>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Sören Gebbert
+
+<p><i>Last changed: $Date: 2012-08-29 14:55:21 +0200 (Mi, 29. Aug 2012) $</i>
Added: grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.py
===================================================================
--- grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.py (rev 0)
+++ grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.py 2013-03-20 18:11:45 UTC (rev 55471)
@@ -0,0 +1,268 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE: t.rast.neighbors
+# AUTHOR(S): Soeren Gebbert
+#
+# PURPOSE: Performs a neighborhood analysis for each map in a space time raster dataset.
+# 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: Performs a neighborhood analysis for each map in a space time raster dataset.
+#% keywords: temporal
+#% keywords: aggregation
+#%end
+
+#%option G_OPT_STRDS_INPUT
+#%end
+
+#%option G_OPT_T_WHERE
+#%end
+
+#%option
+#% key: size
+#% type: integer
+#% description: Neighborhood size
+#% required: no
+#% multiple: no
+#% answer: 3
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% description: Name of the output space time raster dataset
+#% required: yes
+#% multiple: no
+#%end
+
+#%option
+#% key: method
+#% type: string
+#% description: Aggregate operation to be performed on the raster maps
+#% required: yes
+#% multiple: no
+#% options: average,median,mode,minimum,maximum,range,stddev,sum,count,variance,diversity,interspersion,quart1,quart3,perc90,quantile
+#% answer: average
+#%end
+
+
+#%option G_OPT_R_BASE
+#%end
+
+#%option
+#% key: nprocs
+#% type: integer
+#% description: Number of r.neighbor processes to run in parallel
+#% required: no
+#% multiple: no
+#% answer: 1
+#%end
+
+#%flag
+#% key: n
+#% description: Register Null maps
+#%end
+
+from multiprocessing import Process
+import grass.script as grass
+import grass.temporal as tgis
+
+############################################################################
+
+
+def main():
+
+ # Get the options
+ input = options["input"]
+ output = options["output"]
+ where = options["where"]
+ size = options["size"]
+ base = options["base"]
+ register_null = flags["n"]
+ method = options["method"]
+ nprocs = options["nprocs"]
+
+ # Make sure the temporal database exists
+ tgis.init()
+ # We need a database interface
+ dbif = tgis.SQLDatabaseInterfaceConnection()
+ dbif.connect()
+
+ mapset = grass.gisenv()["MAPSET"]
+
+ if input.find("@") >= 0:
+ id = input
+ else:
+ id = input + "@" + mapset
+
+ sp = tgis.SpaceTimeRasterDataset(id)
+
+ if sp.is_in_db() == False:
+ dbif.close()
+ grass.fatal(_("Space time %s dataset <%s> not found") % (
+ sp.get_new_map_instance(None).get_type(), id))
+
+ sp.select(dbif)
+ dummy = sp.get_new_map_instance(None)
+
+ if output.find("@") >= 0:
+ out_id = output
+ else:
+ out_id = output + "@" + mapset
+
+ # The new space time raster dataset
+ new_sp = tgis.SpaceTimeRasterDataset(out_id)
+ if new_sp.is_in_db(dbif):
+ if not grass.overwrite():
+ dbif.close()
+ grass.fatal(_("Space time raster dataset <%s> is already in the "
+ "database, use overwrite flag to overwrite") % out_id)
+
+ rows = sp.get_registered_maps("id,start_time", where, "start_time", dbif)
+
+ if not rows:
+ dbif.close()
+ grass.fatal(_("Space time raster dataset <%s> is empty") % out_id)
+
+ count = 0
+ proc_count = 0
+ proc_list = []
+
+ num_rows = len(rows)
+ new_maps = {}
+
+ for row in rows:
+ count += 1
+
+ grass.percent(count, num_rows, 1)
+
+ map_name = "%s_%i" % (base, count)
+ map_id = dummy.build_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():
+ # Remove the existing temporal database entry
+ new_map.delete(dbif)
+ new_map = sp.get_new_map_instance(map_id)
+ else:
+ grass.error(_("Map <%s> is already in temporal database,"
+ " use overwrite flag to overwrite") %
+ (new_map.get_map_id()))
+ continue
+
+ proc_list.append(Process(target=run_neighbors,
+ args=(row["id"],map_name,method,size)))
+
+ proc_list[proc_count].start()
+ proc_count += 1
+
+ # Join processes if the maximum number of processes are
+ # reached or the end of the loop is reached
+ if proc_count == nprocs or proc_count == num_rows:
+ proc_count = 0
+ exitcodes = 0
+ for proc in proc_list:
+ proc.join()
+ exitcodes += proc.exitcode
+
+ if exitcodes != 0:
+ dbif.close()
+ grass.fatal(_("Error while computation"))
+
+ # Empty process list
+ proc_list = []
+
+ # Store the new maps
+ new_maps[row["id"]] = new_map
+
+ grass.percent(0, num_rows, 1)
+
+ # Insert the new space time dataset
+ if new_sp.is_in_db(dbif):
+ if grass.overwrite():
+ new_sp.delete(dbif)
+ new_sp = sp.get_new_instance(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)
+
+ # collect empty maps to remove them
+ empty_maps = []
+
+ # Register the maps in the database
+ count = 0
+ for row in rows:
+ count += 1
+
+ grass.percent(count, num_rows, 1)
+ # Register the new maps
+ if row["id"] in new_maps:
+ new_map = new_maps[row["id"]]
+ # Read the raster map data
+ new_map.load()
+
+ # In case of a empty map continue, do not register empty maps
+ if new_map.metadata.get_min() is None and \
+ new_map.metadata.get_max() is None:
+ if not register_null:
+ empty_maps.append(new_map)
+ continue
+
+ old_map = sp.get_new_map_instance(row["id"])
+ old_map.select(dbif)
+
+ # 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, unit = old_map.get_relative_time()
+ new_map.set_relative_time(start, end, unit)
+
+ # Insert map in temporal database
+ new_map.insert(dbif)
+ new_sp.register_map(new_map, dbif)
+
+ # Update the spatio-temporal extent and the metadata table entries
+ new_sp.update_from_registered_maps(dbif)
+ grass.percent(num_rows, num_rows, 1)
+
+ # Remove empty maps
+ if len(empty_maps) > 0:
+ names = ""
+ count = 0
+ for map in empty_maps:
+ if count == 0:
+ names += "%s" % (map.get_name())
+ else:
+ names += ",%s" % (map.get_name())
+ count += 1
+
+ grass.run_command("g.remove", rast=names, quiet=True)
+
+
+ dbif.close()
+
+def run_neighbors(input, output, method, size):
+ """Helper function to run r.neighbors in parallel"""
+ return grass.run_command("r.neighbors", input=input, output=output,
+ method=method, size=size, overwrite=grass.overwrite(),
+ quiet=True)
+
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
Property changes on: grass/trunk/temporal/t.rast.neighbors/t.rast.neighbors.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass/trunk/temporal/t.rast.neighbors/test.t.rast.neighbors.sh
===================================================================
--- grass/trunk/temporal/t.rast.neighbors/test.t.rast.neighbors.sh (rev 0)
+++ grass/trunk/temporal/t.rast.neighbors/test.t.rast.neighbors.sh 2013-03-20 18:11:45 UTC (rev 55471)
@@ -0,0 +1,37 @@
+#!/bin/sh
+# Space time raster dataset neighborhood operations
+# 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 -p
+# Generate data
+r.mapcalc --o expr="prec_1 = rand(0, 550)"
+r.mapcalc --o expr="prec_2 = rand(0, 450)"
+r.mapcalc --o expr="prec_3 = rand(0, 320)"
+r.mapcalc --o expr="prec_4 = rand(0, 510)"
+r.mapcalc --o expr="prec_5 = rand(0, 300)"
+r.mapcalc --o expr="prec_6 = rand(0, 650)"
+
+t.create --o type=strds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
+t.register -i --o type=rast input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-15 12:05:45" increment="14 days"
+
+# The first @test
+
+t.rast.neighbors --o input=precip_abs1 output=precip_abs2 base=prec_avg \
+ size=3 method=average nprocs=1
+t.info type=strds input=precip_abs2
+t.rast.neighbors --o input=precip_abs1 output=precip_abs2 base=prec_avg \
+ size=5 method=average nprocs=2
+t.info type=strds input=precip_abs2
+t.rast.neighbors --o input=precip_abs1 output=precip_abs2 base=prec_avg \
+ size=7 method=average nprocs=3
+t.info type=strds input=precip_abs2
+t.rast.neighbors --o input=precip_abs1 output=precip_abs2 base=prec_avg \
+ size=7 method=max nprocs=3
+t.info type=strds input=precip_abs2
+t.rast.neighbors --o input=precip_abs1 output=precip_abs2 base=prec_avg \
+ size=7 method=min nprocs=3
+t.info type=strds input=precip_abs2
+
+t.unregister type=rast maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.remove type=strds input=precip_abs1,precip_abs2
Property changes on: grass/trunk/temporal/t.rast.neighbors/test.t.rast.neighbors.sh
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list