[GRASS-SVN] r66870 - in grass-addons/grass7/raster: . r.mregression.series r.series.decompose
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Nov 19 23:43:28 PST 2015
Author: DmitryKolesov
Date: 2015-11-19 23:43:28 -0800 (Thu, 19 Nov 2015)
New Revision: 66870
Added:
grass-addons/grass7/raster/r.series.decompose/
grass-addons/grass7/raster/r.series.decompose/Makefile
grass-addons/grass7/raster/r.series.decompose/r.series.decompose.html
grass-addons/grass7/raster/r.series.decompose/r.series.decompose.py
Modified:
grass-addons/grass7/raster/r.mregression.series/r.mregression.series.py
Log:
Add r.series.decompose (draft version, refactoring is required)
Modified: grass-addons/grass7/raster/r.mregression.series/r.mregression.series.py
===================================================================
--- grass-addons/grass7/raster/r.mregression.series/r.mregression.series.py 2015-11-20 06:01:25 UTC (rev 66869)
+++ grass-addons/grass7/raster/r.mregression.series/r.mregression.series.py 2015-11-20 07:43:28 UTC (rev 66870)
@@ -7,7 +7,7 @@
#
# AUTHOR(S): Dmitry Kolesov <kolesov.dm at gmail.com>
#
-# PURPOSE: multiple regression between two time series
+# PURPOSE: multiple regression between several time series
#
# COPYRIGHT: (C) 2015 by the GRASS Development Team
#
@@ -18,7 +18,7 @@
#############################################################################
#%Module
-#% description: Program calculates multiple regression between two series: Y = a1*X1 + ... + an*Xn
+#% description: Program calculates multiple regression between time series: Y = b1*X1 + ... + bn*Xn
#% overwrite: yes
#%End
#%option
@@ -244,7 +244,6 @@
b.put(r, c, coefs[i])
finally:
self.close_rasters()
- print 'ols'
def main(options, flags):
Added: grass-addons/grass7/raster/r.series.decompose/Makefile
===================================================================
--- grass-addons/grass7/raster/r.series.decompose/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.series.decompose/Makefile 2015-11-20 07:43:28 UTC (rev 66870)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR =../..
+
+PGM = r.series.decompose
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.series.decompose/r.series.decompose.html
===================================================================
--- grass-addons/grass7/raster/r.series.decompose/r.series.decompose.html (rev 0)
+++ grass-addons/grass7/raster/r.series.decompose/r.series.decompose.html 2015-11-20 07:43:28 UTC (rev 66870)
@@ -0,0 +1,148 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.series.decompose</em> is a module to calculate decomposintion of signal X.
+<div class="code"><pre>
+X(t) = B0 + B1*t + B2*sin(B1*t) + B3 * cos(B1*t) + ... + B{n-1}*sin(Bk*t) + Bn * cos(Bk*t) + e
+</pre></div>
+<p>
+where <em>X</em> is a raster time series, <em>t</em> is time (<em>t</em>
+in <em>[0, pi]</em>), <em>sin(Fi*t)</em> and <em>cos(Fi*t)</em> are time
+variables; <em>Fi</em> are user specifed frequencies; <em>e</em> is a error.
+<p>
+The module used r.mregression.series to find the regression coefficients
+<em>Bi</em>, then it produces the fitted rasters series <em>X(t)</em>
+using the coefficients.
+<p>
+So the module makes each output cell value a function of the values
+assigned to the corresponding cells in the time variable raster map series
+and the rasters of the coefficients.
+
+<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(t)</em>
+<p>
+<em>coef_prefix</em> Prefix for names of result raster
+ (rasters of coefficients)
+<p>
+<em>timevar_prefix</em> Prefix for names of result raster
+ (rasters of time variables)
+<p>
+<em>freq</em> List of frequencies for sin and cos functions
+<p>
+
+
+<h2>NOTES</h2>
+
+<em>X</em> must be equally spaced time serie. If the serie isn't equally
+spaced, insert NULL raster maps into <em>X</em>.
+<p>
+The list of inputs for each cell (including NULLs) is passed to the
+regression function. The functin compute the parameters over the
+non-NULL values, producing a NULL result only if there aren't enough
+non-NULL values for computing.
+
+The regression coefficients <em>Bi</em> are stored in raster maps.
+They can be used for construct more detail time series via the equation:
+<div class="code"><pre>
+X(t) = B0 + B1*t + B2*sin(B1*t) + B3 * cos(B1*t) + ... + B{n-1}*sin(Bk*t) + Bn * cos(Bk*t) + e
+</pre></div>
+<p>
+To do that the user have to create time variables (<em>t</em>,
+<em>sin(Fi*t)</em> and <em>cos(Fi*t)</em>) at desired time <em>T0</em>
+and then use r.mapcalc to produce the <em>X(T0)</em>.
+
+<p>
+The maximum number of raster maps to be processed is limited by the
+operating system. For example, both the hard and soft limits are
+typically 1024. The soft limit can be changed with e.g. <tt>ulimit -n
+1500</tt> (UNIX-based operating systems) but not higher than the hard
+limit. If it is too low, you can as superuser add an entry in
+
+<div class="code"><pre>
+/etc/security/limits.conf
+# <domain> <type> <item> <value>
+your_username hard nofile 1500
+</pre></div>
+
+This would raise the hard limit to 1500 file. Be warned that more
+files open need more RAM.
+
+
+<h2>EXAMPLES</h2>
+Suppose we have time series of MODIS NDVI data (from 01 jan to 27 dec):
+<pre>
+> g.mlist rast pattern="mod*", separator=','
+mod2003_01_01,mod2003_01_09,...,mod2003_12_27
+</pre>
+<p>
+We use one year data, so we suppose the there is a half of sinusoid signal
+in the data (NDVI values icrease, then decrease usualy). So 01 jan is t0==0,
+27 dec is tn==2*pi, there is a frequence 0.5 in the data (and there are more
+frequencies, for example 1.0 and 1.5).
+<p>
+Decompose the signal:
+<pre>
+> maps = $(g.list rast pattern="mod*", separator=',')
+> r.series.decompose input=$maps coef_prefix="coef." \
+ timevar_prefix="dec." result_pref="res." \
+ freq=0.5,1.0,1.5
+</pre>
+<p>
+The command creates rasters of the coefficiens <em>coef.*</em>:
+<pre>
+coef.const
+coef.time
+coef.sin_fr0.5
+coef.cos_fr0.5
+coef.sin_fr1.0
+coef.cos_fr1.0
+coef.cos_fr1.5
+coef.sin_fr1.5
+</pre>
+and rasters of fitted NDVI <em>res.*</em>:
+<pre>
+res.mod2003_01_01
+res.mod2003_01_09
+...
+</pre>
+<p>
+To compute NDVI for 03 jan we need: (1) find time <em>T</em> for 03 jan
+(2) create time variables for 02 jan.
+<p>
+The length (in days) of the NDVI time series is 362, 03 jan
+is the third day of the series, so <em>T</em> = <em>3 * (2*pi/362)</em>
+radians. But r.mapcalc uses degrees for <em>sin()</em> and <em>cos()</em>
+functions. So <em>T</em> = <em>3 * 360/362</em> degrees.
+<p>
+Create time variables:
+<pre>
+r.mapcalc "T = 3.0 * 360.0/362.0"
+r.mapcalc "sin0.5 = sin(0.5*3.0*360.0/362.0)"
+r.mapcalc "cos0.5 = cos(0.5*3.0*360.0/362.0)"
+r.mapcalc "sin1 = sin(3.0*360.0/362.0)"
+r.mapcalc "cos1 = cos(3.0*360.0/362.0)"
+r.mapcalc "sin1.5 = sin(1.5*3.0*360.0/362.0)"
+r.mapcalc "cos1.5 = cos(1.5*3.0*360.0/362.0)"
+</pre>
+<p>
+Create NDVI for 03 jan:
+<pre>
+r.mapcals "ndvi03jan = coef.const + coef.time*T +\
+ coef.sin_fr0.5*sin0.4 + coef.cos_fr0.5*cos0.5 +\
+ coef.sin_fr1.0*sin1 + coef.cos_fr1.0*cos1 +\
+ coef.sin_fr1.5*sin1.5 + coef.cos_fr1.5*cos1.5"
+</pre>
+
+<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: 2014-10-29 02:11:08 +0300 (Ср., 29 окт. 2014) $</i>
Added: grass-addons/grass7/raster/r.series.decompose/r.series.decompose.py
===================================================================
--- grass-addons/grass7/raster/r.series.decompose/r.series.decompose.py (rev 0)
+++ grass-addons/grass7/raster/r.series.decompose/r.series.decompose.py 2015-11-20 07:43:28 UTC (rev 66870)
@@ -0,0 +1,238 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+############################################################################
+#
+# MODULE: r.series.decompose
+#
+# AUTHOR(S): Dmitry Kolesov <kolesov.dm at gmail.com>
+#
+# PURPOSE: Perform decomposition of time series X :
+# X(t) = b0 + b1 *t + sin(2*pi*t/N) + cos(2*pi*t/N) + ... + sin(2*pi*t/k*N) + cos(2*pi*t/k*N)
+#
+# 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: Program calculates decomposition of time series X
+#% 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 filterd X(t)
+#% required : yes
+#% multiple: no
+#%end
+#%option
+#% key: coef_prefix
+#% type: string
+#% gisprompt: prefix for raster names of decomposition coefficients
+#% description: Prefix for names of result raster (rasters of coefficients)
+#% required : yes
+#% multiple: no
+#%end
+#%option
+#% key: timevar_prefix
+#% type: string
+#% gisprompt: prefix for raster names of time variables
+#% description: Prefix for names of result raster (rasters of time variables)
+#% required : yes
+#% multiple: no
+#%end
+#%option
+#% key: freq
+#% type: double
+#% gisprompt: list of used frequiences
+#% description: List of frequencies for sin and cos functions
+#% required : yes
+#% multiple: yes
+#%end
+
+
+
+
+import os
+import sys
+import csv
+import uuid
+
+from math import pi, degrees
+
+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
+
+
+def get_time(N, t, deg=True):
+ """Returns point number t from linspace [0, 2*pi]
+ if deg == True, return degrees
+ """
+
+ step = 2 * pi / N
+ x = step * t
+ if deg:
+ x = degrees(x)
+ return x
+
+def _freq_to_name(freq):
+ return "_fr%s" % (freq)
+
+def _time_to_name(time):
+ return "_t%04d" % (time, )
+
+
+def _generate_time(N, prefix):
+ """Generate (constant) rasters of time line
+ t0, t1, t2, ..., t_{N-1})
+ where t0 = 0, t_{N-1} = 2*pi
+
+ N -- count of result rasters
+ prefix -- prefix for the raster names
+
+ return list of names for the result rasters
+ """
+ assert N > 1
+ names = dict()
+ for i in range(N):
+ output = prefix + _time_to_name(i)
+ t = get_time(N, i, deg=True)
+ grass.mapcalc("${out} = ${t}", out=output, t=t,
+ quiet=True, overwrite=grass.overwrite())
+ names[i] = output
+
+ return names
+
+
+def _generate_harmonics(time_names, freq, prefix):
+ """Generate (constant) rasters of harmonics
+ t0, t1, t2, ..., t_{N-1})
+ where t0 = 0, t_{N-1} = 2*pi
+
+ freq -- frequencies for sin() and cos() rasters
+ prefix -- prefix for the raster names
+
+ return list of names for the result rasters
+ """
+ names = dict()
+ for i in time_names:
+ harm = dict()
+ t = get_time(len(time_names), i)
+ for f in freq:
+ sin_output = prefix + 'sin' + _time_to_name(i) + _freq_to_name(f)
+ grass.mapcalc("${out} = sin(${f} * ${t})", out=sin_output, t=t, f=f,
+ quiet=True, overwrite=grass.overwrite())
+
+ cos_output = prefix + 'cos' + _time_to_name(i) + _freq_to_name(f)
+ grass.mapcalc("${out} = cos(${f} * ${t})", out=cos_output, t=t, f=f,
+ quiet=True, overwrite=grass.overwrite())
+ harm[f] = dict(sin=sin_output, cos=cos_output)
+
+ names[i] = harm
+
+ return names
+
+def _generate_const(prefix):
+ output = prefix + 'const'
+ grass.mapcalc("${out} = 1.0", out=output,
+ quiet=True, overwrite=grass.overwrite())
+ return output
+
+def generate_vars(time_count, freq, prefix):
+ """Generate time_count sets of variables.
+ """
+ const_name = _generate_const(prefix)
+ time_names = _generate_time(time_count, prefix)
+ harm_names = _generate_harmonics(time_names, freq, prefix)
+
+ return const_name, time_names, harm_names
+
+def _generate_sample_descr(fileobj, freq, xnames, const_name, time_names, harm_names):
+ """Generate settings file for r.mregression.series
+
+ """
+ freq_names = []
+ for f in freq:
+ freq_names.append('sin' + _freq_to_name(f))
+ freq_names.append('cos' + _freq_to_name(f))
+ header = ['x', 'const', 'time'] + freq_names
+
+ writer = csv.writer(fileobj, delimiter=',')
+ writer.writerow(header)
+
+ size = len(xnames)
+ for i in range(size):
+ row = [xnames[i], const_name, time_names[i]]
+ f_name = harm_names[i]
+ for f in freq:
+ sin_name = f_name[f]['sin']
+ cos_name = f_name[f]['cos']
+ row.append(sin_name)
+ row.append(cos_name)
+ writer.writerow(row)
+
+def regression(settings_name, coef_prefix):
+ grass.run_command('r.mregression.series',
+ samples=settings_name, result_prefix=coef_prefix,
+ overwrite=grass.overwrite())
+
+
+def inverse_transform(settings_name, coef_prefix, result_prefix='res.'):
+ reader = csv.reader(open(settings_name), delimiter=',')
+ header = reader.next()
+
+ data_names = [coef_prefix + name for name in header[1:]]
+
+ for row in reader:
+ s = "%s%s = " % (result_prefix, row[0])
+ sums = []
+ for i in range(len(data_names)):
+ sums.append("%s*%s" % (data_names[i], row[i+1]))
+ s += ' + '.join(sums)
+
+ grass.mapcalc(s, overwrite=grass.overwrite(), quite=True)
+
+def main(options, flags):
+ xnames = options['input']
+ coef_pref = options['coef_prefix']
+ timevar_pref = options['timevar_prefix']
+ result_pref = options['result_prefix']
+ freq = options['freq']
+ freq = [float(f) for f in freq.split(',')]
+
+ xnames = xnames.split(',')
+
+ N = len(xnames)
+ if len(freq) >= (N-1)/2:
+ grass.error("Count of used harmonics is to large. Reduce the paramether.")
+ sys.exit(1)
+
+ const_name, time_names, harm_names = generate_vars(N, freq, timevar_pref)
+
+ settings_name = uuid.uuid4().hex
+ settings = open(settings_name, 'w')
+ _generate_sample_descr(settings, freq, xnames, const_name, time_names, harm_names)
+ settings.close()
+ regression(settings_name, coef_pref)
+ inverse_transform(settings_name, coef_pref, result_pref)
+ os.unlink(settings_name)
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main(options, flags)
More information about the grass-commit
mailing list