[GRASS-SVN] r67589 - grass-addons/grass7/raster/r.series.diversity
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jan 15 03:06:44 PST 2016
Author: pvanbosgeo
Date: 2016-01-15 03:06:44 -0800 (Fri, 15 Jan 2016)
New Revision: 67589
Modified:
grass-addons/grass7/raster/r.series.diversity/r.series.diversity.html
grass-addons/grass7/raster/r.series.diversity/r.series.diversity.py
Log:
r.series.diversity addon: rewrote some expressions in the script to avoid "in-place modifications in r.mapcalc. Changed computation to avoid large number of temporary layers. Corrections in description of indices in manual page.
Modified: grass-addons/grass7/raster/r.series.diversity/r.series.diversity.html
===================================================================
--- grass-addons/grass7/raster/r.series.diversity/r.series.diversity.html 2016-01-15 11:06:39 UTC (rev 67588)
+++ grass-addons/grass7/raster/r.series.diversity/r.series.diversity.html 2016-01-15 11:06:44 UTC (rev 67589)
@@ -76,21 +76,21 @@
equivalent to <i>-1 * 1 / exp(R2)</i>, with <i>R2</i> the renyi
index for <i>alpha=2</i>. With this index, 0 represents infinite
diversity and 1, no diversity. This is counterintuitive behavior for
-a diversity index. An alternative is the inverse Simpson index is
-defined as <i>ID = 1 / D)</i>. The index represents the probability
-that two individuals randomly selected from a sample will belong to
-different species. The value ranges between 0 and 1, with greater
-values representing greater sample diversity. The name of the output
-layer is composed of the basename + invsimpson.
+a diversity index. An alternative is the inverse Simpson index,
+which is defined as <i>ID = 1 / D)</i>. The lowest value of this
+index is 1 and represent a community containing only one species
+(Gini, 1912; Simpson, 1949). The higher the value, the greater the
+diversity. The maximum value is the number of species in the sample.
+The name of the output layer is composed of the basename + invsimpson.
-<h4>Gini-Simpson index</h4>
+<h4>Gini-Simpson index (Simpson's index of diversity)</h4>
An alternative way to overcome the problem of the counter-intuitive
-nature of Simpson's Index is to use <i>1 - D)</i>. The lowest value
-of this index is 1 and represent a community containing only one
-species (Gini, 1912; Simpson, 1949). The higher the value, the
-greater the diversity. The maximum value is the number of species in
-the sample. The name of the output layer is composed of the basename
+nature of Simpson's Index is to use <i>1 - D)</i>. The index
+represents the probability that two individuals randomly selected
+from a sample will belong to different species. The value ranges
+between 0 and 1, with greater values representing greater sample
+diversity. The name of the output layer is composed of the basename
+ ginisimpson.
<h2>NOTES</h2>
@@ -102,6 +102,10 @@
below). These functions requires one input layer and compute the
diversity using a moving window.
+<p>Currently when working with very large raster layers and many
+input layers, computations can take a long time. On the todo list:
+find a way to reduce this.
+
<h2>EXAMPLES</h2>
Suppose we have five layers, each representing number of
Modified: grass-addons/grass7/raster/r.series.diversity/r.series.diversity.py
===================================================================
--- grass-addons/grass7/raster/r.series.diversity/r.series.diversity.py 2016-01-15 11:06:39 UTC (rev 67588)
+++ grass-addons/grass7/raster/r.series.diversity/r.series.diversity.py 2016-01-15 11:06:44 UTC (rev 67589)
@@ -106,6 +106,12 @@
#% guisection: Indices
#%end
+#%flag
+#% key: t
+#% description: Total counts
+#% guisection: Indices
+#%end
+
#%rules
#% required: -r,-s,-h,-e,-p,-g,-n
#%end
@@ -162,6 +168,7 @@
#--------------------------------------------------------------------------
# Layers
+ grass.info("Preparing...")
OUT = options['output']
IN = options['input']
IN = IN.split(',')
@@ -175,6 +182,7 @@
flag_p = flags['p']
flag_g = flags['g']
flag_n = flags['n']
+ flag_t = flags['t']
if options['alpha']:
Q = map(float, options['alpha'].split(','))
else:
@@ -204,46 +212,50 @@
#--------------------------------------------------------------------------
# Renyi entropy
#--------------------------------------------------------------------------
- tmpt = tmpname("sht")
- clean_rast.add(tmpt)
- grass.info("Computing the sum across all input layers")
- grass.run_command("r.series", quiet=True, output=tmpt,
- input=IN, method="sum")
+ tmp_1 = tmpname("sht")
+ clean_rast.add(tmp_1)
+ grass.info("Computing the sum across all input layers (this may take a while)")
+ grass.run_command("r.series", quiet=True, output=tmp_1, input=IN, method="sum")
for n in xrange(len(Q)):
grass.info(_("Computing alpha = {n}").format(n=Q[n]))
Qn = str(Q[n])
Qn = Qn.replace('.', '_')
- out_renyi = OUT + "_Renyi_" + Qn
+ renyi = OUT + "_Renyi_" + Qn
if Q[n] == 1:
# If alpha = 1
- grass.mapcalc("$tmpa = 0", tmpa=out_renyi, quiet=True)
+ # TODO See email 14-01-16 about avoiding loop below
+ grass.mapcalc("$renyi = 0", renyi=renyi, quiet=True)
for i in xrange(len(IN)):
- grass.mapcalc("$tmpta = $tmpta - (($inl/$tmpt) * log(($inl/$tmpt)))",
- tmpta=out_renyi,
- inl=IN[i],
- tmpt=tmpt,
- overwrite=True,
- quiet=True)
+ grass.info(_("Computing map {j} from {n} maps").format(j=i+1, n=len(IN)))
+ tmp_2 = tmpname("sht")
+ clean_rast.add(tmp_2)
+ grass.mapcalc("$tmp_2 = if($inl==0, $renyi, $renyi - \
+ (($inl/$tmp_1) * log(($inl/$tmp_1))))",
+ renyi=renyi, tmp_2=tmp_2,
+ inl=IN[i], tmp_1=tmp_1, quiet=True)
+ grass.run_command("g.rename", raster='%s,%s' % (tmp_2,renyi),
+ overwrite=True, quiet=True)
else:
# If alpha != 1
- tmptb = tmpname("sht")
- clean_rast.add(tmptb)
- grass.mapcalc("$tmpb = 0", tmpb=tmptb, quiet=True)
+ tmp_3 = tmpname("sht")
+ clean_rast.add(tmp_3)
+ tmp_4 = tmpname("sht")
+ clean_rast.add(tmp_4)
+ grass.mapcalc("$tmp_3 = 0", tmp_3=tmp_3, quiet=True)
for i in xrange(len(IN)):
- grass.mapcalc("$tmptb = if($inl==0,$tmptb+0, $tmptb+pow($inl/$tmpt,$alpha))",
- tmptb=tmptb,
- tmpt=tmpt,
- inl=IN[i],
- alpha=Q[n],
- quiet=True,
- overwrite=True)
- grass.mapcalc("$outl = (1/(1-$alpha)) * log($tmptb)",
- outl=out_renyi,
- tmptb=tmptb,
- alpha=Q[n],
- quiet=True)
+ grass.info(_("Computing map {j} from {n} maps").format(j=i+1, n=len(IN)))
+ grass.mapcalc("$tmp_4 = if($inl==0,$tmp_3+0, $tmp_3 + \
+ pow($inl/$tmp_1,$alpha))",
+ tmp_3=tmp_3, tmp_4=tmp_4,
+ tmp_1=tmp_1, inl=IN[i],
+ alpha=Q[n], quiet=True)
+ grass.run_command("g.rename", raster='%s,%s' % (tmp_4,tmp_3),
+ overwrite=True, quiet=True)
+ grass.mapcalc("$outl = (1/(1-$alpha)) * log($tmp_3)",
+ outl=renyi, tmp_3=tmp_3,
+ alpha=Q[n], quiet=True)
grass.run_command("g.remove", type="raster",
- name=tmptb, flags="f", quiet=True)
+ name=tmp_3, flags="f", quiet=True)
#--------------------------------------------------------------------------
# Species richness
@@ -336,10 +348,16 @@
grass.run_command("g.remove", flags="f", type="raster",
name=in_div, quiet=True)
- # Clean up temporary files
- # grass.run_command("g.remove", type="raster", name=tmpt,
- # flags="f", quiet=True)
+ #--------------------------------------------------------------------------
+ # Total count (base unit, like individuals)
+ #--------------------------------------------------------------------------
+ if flag_t:
+ rast = OUT + "_count"
+ grass.run_command("g.rename", raster=(tmp_1,rast), quiet=True)
+ else:
+ grass.run_command("g.remove", type="raster", name=tmp_1, flags="f", quiet=True)
+
if __name__ == "__main__":
options, flags = grass.parser()
atexit.register(cleanup)
More information about the grass-commit
mailing list