[GRASS-SVN] r66954 - in grass-addons/grass7/raster: . r.series.filter
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 27 18:24:21 PST 2015
Author: DmitryKolesov
Date: 2015-11-27 18:24:21 -0800 (Fri, 27 Nov 2015)
New Revision: 66954
Added:
grass-addons/grass7/raster/r.series.filter/
grass-addons/grass7/raster/r.series.filter/Makefile
grass-addons/grass7/raster/r.series.filter/flt.png
grass-addons/grass7/raster/r.series.filter/r.series.filter.html
grass-addons/grass7/raster/r.series.filter/r.series.filter.py
grass-addons/grass7/raster/r.series.filter/signal.png
Log:
r.series.filter module
Added: grass-addons/grass7/raster/r.series.filter/Makefile
===================================================================
--- grass-addons/grass7/raster/r.series.filter/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.series.filter/Makefile 2015-11-28 02:24:21 UTC (rev 66954)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR =../..
+
+PGM = r.series.filter
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.series.filter/flt.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.series.filter/flt.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.series.filter/r.series.filter.html
===================================================================
--- grass-addons/grass7/raster/r.series.filter/r.series.filter.html (rev 0)
+++ grass-addons/grass7/raster/r.series.filter/r.series.filter.html 2015-11-28 02:24:21 UTC (rev 66954)
@@ -0,0 +1,66 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.series.filter</em> is a module to filter raster time series <em>X</em>
+in time domain.
+<p>
+<em>input</em> Raster names of equally spaced time series <em>X</em>.
+<p>
+<em>result_prefix</em> Prefix for raster names of filterd <em>X</em>.
+<em>method</em> filtering method. Savitzky-Golay filter is only implemented now.
+<em>winsize</em> The length of the filter window. <em>winsize</em> must be a positive odd integer.
+<em>order</em> The order of the polynomial used to fit the samples. <em>order</em> must be less than <em>winsize</em>.
+
+<h2>NOTES</h2>
+
+<em>X</em> must be equally spaced time series. If the series isn't equally
+spaced, insert NULL raster maps into <em>X</em>.
+
+
+<h2>EXAMPLES</h2>
+Create test data: <em>X = sin(t) + E</em>,
+where <em>X</em> is raster time series, <em>E</em> is a error term.
+<pre>
+for T in $(seq -w 0 10 360)
+do
+ name="test_raster"$T
+ r.mapcalc -s "$name = sin($T) + rand(-0.3, 0.3)"
+done
+</pre>
+<p>
+Create smooth raster series using Savitzky-Golay method:
+<div class="code"><pre>
+maps=$(g.list rast patt="test_*" sep=,)
+r.series.filter input=$maps result_prefix="flt." method=savgol winsize=9 order=2 --o
+</pre></div>
+
+<p>
+Look at the result (plot the curves for a pixel):
+<div class="code"><pre>
+maps=$(g.list rast patt="test_*" sep=,)
+fmaps=$(g.list rast patt="flt.*" sep=,)
+
+eval $(g.region -cg)
+i.spectral -g raster=$maps coor=$center_easting,$center_northing out=signal.png
+i.spectral -g raster=$fmaps coor=$center_easting,$center_northing out=flt.png
+</pre></div>
+
+<center>
+<img src="signal.png" border="0" alt="Values of the time series in the central point">
+<img src="flt.png" border="0" alt="Filtered values of the time series in the central point">
+</center>
+
+
+
+
+<h2>SEE ALSO</h2>
+
+<em><a href="addons/r.mregression.series.html">r.regression.series</a></em>,
+<em><a href="r.series.html">r.series</a></em>,
+<em><a href="r.regression.line.html">r.regression.line</a></em>,
+<em><a href="g.list.html">g.list</a></em>,
+
+<h2>AUTHOR</h2>
+
+Dmitry Kolesov
+
+<p><i>Last changed: $Date: 2015-11-20 17:36:50 +0300 (Пт., 20 нояб. 2015) $</i>
Added: grass-addons/grass7/raster/r.series.filter/r.series.filter.py
===================================================================
--- grass-addons/grass7/raster/r.series.filter/r.series.filter.py (rev 0)
+++ grass-addons/grass7/raster/r.series.filter/r.series.filter.py 2015-11-28 02:24:21 UTC (rev 66954)
@@ -0,0 +1,190 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+############################################################################
+#
+# MODULE: r.series.filter
+#
+# AUTHOR(S): Dmitry Kolesov <kolesov.dm at gmail.com>
+#
+# PURPOSE: Perform filtering of raster time series X
+#
+# COPYRIGHT: (C) 2015 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+#%Module
+#% description: Perform filtering of raster time series X (in time domain)
+#% overwrite: yes
+#%End
+#%option
+#% key: input
+#% type: string
+#% gisprompt: list of raster names
+#% description: Raster names of equally spaced time series.
+#% required : yes
+#% multiple: yes
+#%end
+#%option
+#% key: result_prefix
+#% type: string
+#% gisprompt: prefix of result raster names
+#% description: Prefix for raster names of filtered X(t)
+#% required : yes
+#% multiple: no
+#%end
+#%option
+#% key: method
+#% type: string
+#% required : no
+#% multiple: no
+#% answer: savgol
+#% description: Used method
+#% descriptions: savgol; Savitzky–Golay filter
+#%end
+#%option
+#% key: winsize
+#% type: integer
+#% answer: 5
+#% required: no
+#% multiple: no
+#% description: Length of running window for the filter
+#%end
+#%option
+#% key: order
+#% type: integer
+#% answer: 2
+#% required: no
+#% multiple: no
+#% description: Order of the Savitzky-Golay filter
+#%end
+
+
+
+import os
+import sys
+
+import numpy as np
+from scipy.signal import savgol_filter
+
+if "GISBASE" not in os.environ:
+ sys.stderr.write("You must be in GRASS GIS to run this program.\n")
+ sys.exit(1)
+
+import grass.script as grass
+from grass.pygrass import raster
+from grass.pygrass.raster.buffer import Buffer
+from grass.exceptions import OpenError
+from grass.pygrass.gis.region import Region
+
+CNULL = -2147483648 # null value for CELL maps
+FNULL = np.nan # null value for FCELL and DCELL maps
+
+
+def init_rasters(names):
+ """Get list of raster names,
+ return array of the rasters
+ """
+ rasters = []
+ for name in names:
+ r = raster.RasterSegment(name)
+ rasters.append(r)
+ return rasters
+
+
+def open_rasters(raster_list, write=False):
+ for r in raster_list:
+ try:
+ if write:
+ if r.exist():
+ r.open('w', 'DCELL', overwrite=grass.overwrite())
+ else:
+ r.open('w', 'DCELL')
+ else:
+ r.open()
+ except OpenError:
+ grass.error("Can't open raster %s" % (r.name, ))
+ sys.exit(1)
+
+def close_rasters(raster_list):
+ for r in raster_list:
+ if r.is_open():
+ r.close()
+
+
+def _filter(row_data, winsize, order):
+ result = np.empty(row_data.shape)
+ _, cols = row_data.shape
+ for i in range(cols):
+ arr = _fill_nulls(row_data[:, i])
+ arr = savgol_filter(arr, winsize, order)
+ result[:, i] = arr
+
+ return result
+
+def _fill_nulls(arr):
+ """Fill no-data values in arr
+ Return np.array with filled data
+ """
+ nz = lambda z: z.nonzero()[0]
+ nans = np.isnan(arr)
+ arr[nans] = np.interp(nz(nans), nz(~nans), arr[~nans])
+
+ return arr
+
+def filter(names, winsize, order, prefix):
+ inputs = init_rasters(names)
+ output_names = [prefix + name for name in names]
+ outputs = init_rasters(output_names)
+ try:
+ open_rasters(outputs, write=True)
+ open_rasters(inputs)
+
+ reg = Region()
+ for i in range(reg.rows):
+ row_data = np.array([r.get_row(i) for r in inputs])
+ filtered_rows = _filter(row_data, winsize, order)
+ # import ipdb; ipdb.set_trace()
+ for map_num in range(len(outputs)):
+ map = outputs[map_num]
+ row = filtered_rows[map_num]
+ buf = Buffer(row.shape, map.mtype, row)
+ map.put_row(i, buf)
+ finally:
+ close_rasters(outputs)
+ close_rasters(inputs)
+
+def main(options, flags):
+ xnames = options['input']
+ xnames = xnames.split(',')
+
+ winsize = options['winsize']
+ winsize = int(winsize)
+
+ order = options['order']
+ order = int(order)
+
+ res_prefix = options['result_prefix']
+
+ N = len(xnames)
+ if N < winsize:
+ grass.error("The used running window size is to big. Decrease the paramether or add more rasters to the series.")
+ sys.exit(1)
+
+ _, rem = divmod(winsize, 2)
+ if rem == 0:
+ grass.error("Window length must be odd.")
+ sys.exit(1)
+
+ if order >= winsize:
+ grass.error("Order of the filter must be less than window length")
+
+ filter(xnames, winsize, order, res_prefix)
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main(options, flags)
Added: grass-addons/grass7/raster/r.series.filter/signal.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.series.filter/signal.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
More information about the grass-commit
mailing list