[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