[GRASS-SVN] r64121 - grass-addons/grass7/raster/r.eb

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Jan 13 02:45:37 PST 2015


Author: pvanbosgeo
Date: 2015-01-13 02:45:37 -0800 (Tue, 13 Jan 2015)
New Revision: 64121

Added:
   grass-addons/grass7/raster/r.eb/r.meb.html
   grass-addons/grass7/raster/r.eb/r.meb.py
   grass-addons/grass7/raster/r.eb/r_meb_concept.png
Removed:
   grass-addons/grass7/raster/r.eb/r.eb.html
   grass-addons/grass7/raster/r.eb/r.eb.py
   grass-addons/grass7/raster/r.eb/r_eb_concept.png
Modified:
   grass-addons/grass7/raster/r.eb/Makefile
Log:
Noticed that there is a module i.eb. To avoid confusion, changed the name to r.meb (multivariate environmental bias)

Modified: grass-addons/grass7/raster/r.eb/Makefile
===================================================================
--- grass-addons/grass7/raster/r.eb/Makefile	2015-01-13 10:38:18 UTC (rev 64120)
+++ grass-addons/grass7/raster/r.eb/Makefile	2015-01-13 10:45:37 UTC (rev 64121)
@@ -1,6 +1,6 @@
 MODULE_TOPDIR = ../..
 
-PGM = r.eb
+PGM = r.meb
 
 include $(MODULE_TOPDIR)/include/Make/Script.make
 

Deleted: grass-addons/grass7/raster/r.eb/r.eb.html
===================================================================
--- grass-addons/grass7/raster/r.eb/r.eb.html	2015-01-13 10:38:18 UTC (rev 64120)
+++ grass-addons/grass7/raster/r.eb/r.eb.html	2015-01-13 10:45:37 UTC (rev 64121)
@@ -1,54 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<p>The environmental bias (EB) takes the medium conditions in an 
-area <em>N</em> and computes how much conditions in a subset of <em>N
-</em> (<em>S</em>) deviate from these medium conditions. This can for 
-example be used to see how well conditions in the protected areas of 
-a country represent conditions in the whole country.
-
-<p>The messure is based on the Multivariate Environmental Similarity 
-(<em>MES</em>) surface, which was proposed by Elith et al (2010). To compute 
-the MES First the similarity of a point <em>P</em> to the conditions 
-in <em>N</em> with respect to variable <em>V</em> is computed. The 
-similarity is expressed as the deviation from the median of <em>V</em> 
-in <em>P</em> to those in <em>N</em>. This is done for all variables of 
-interest (V<sub>1</sub>, V<sub>2</sub>, ...V<sub>j</sub>).
-
-<p>In the original equation proposed by Elith et al (2010) the final 
-<em>MES</em> of <em>P</em> is computed as the minimum of the 
-similarity values (<em>IES<sub>minimum</sub></em>) of the individual 
-variables (<em>V<sub>j</sub></em>) in <em>P</em>. In <em>r.eb</em> 
-there is the option to use the mean (<em>IES<sub>mean</sub></em>) or 
-median (<em>IES<sub>median</sub></em>) of the <em>IES</em> instead. 
-For the EB these may be a better choice as they take into account 
-the similarity along all environmnetal axes and not only the one 
-that deviates most. 
-
-<p>The <em>EB</em> is computed as the absolute difference of the median 
-of the <em>MES</em> in the whole target area (MES<sub>n</sub>) and the 
-median of the <em>MES</em> in the subset (<em>MES<sub>s</sub></em>), 
-divided by the median absolute deviation 
-(<em><a href="http://en.wikipedia.org/wiki/Median_absolute_deviation">MAD</a></em>
-) of the <em>MES</em> in <em>N</em>. It is also possible to compute 
-the <em>EB</em> based on the individual variables in which case the 
-<em>IES</em> instead of the <em>MES</em> is used. 
-
-<p><img src="r_eb_concept.png">
-
-<h2>REFERENCES</h2> 
-
-<p>Elith, J., Kearney, M., and Phillips, S. 
-2010. The art of modelling range-shifting species. Methods in 
-Ecology and Evolution 1:330-342.
-<p>van Breugel et al. Environmental gap analysis to prioritize 
-conservation efforts in eastern Africa. Submitted to Plos One
-(todo: add full reference when accepted)
-<h2>SEE ALSO</h2>
-
-<em><a href="r.mess.html">r.mess</a></em>
-
-<h2>AUTHOR</h2>
-
-Paulo van Breugel, paulo at ecodiv.org
-
-<p><em>Last changed: $Date$</em>

Deleted: grass-addons/grass7/raster/r.eb/r.eb.py
===================================================================
--- grass-addons/grass7/raster/r.eb/r.eb.py	2015-01-13 10:38:18 UTC (rev 64120)
+++ grass-addons/grass7/raster/r.eb/r.eb.py	2015-01-13 10:45:37 UTC (rev 64121)
@@ -1,382 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-########################################################################
-#
-# MODULE:       r.eb
-# AUTHOR(S):    Paulo van Breugel <p.vanbreugel AT gmail.com>
-# PURPOSE:      Compute the envirionmental bias (EB). If A is an areas within a
-#               larger region B, the EB represents how much envirionmental
-#               conditions in A deviate from median conditions in B. The first
-#               step is to compute the multi-envirionmental similarity (MES)
-#               for B, using all raster cells in B as reference points. The MES
-#               of a raster cell thus represent how much conditions deviate
-#               from median conditions in B. The EB is then computed as the
-#               absolute difference of the median of MES values in A (MESa)
-#               and median of MES values in B (MESb), divided by the median of
-#               the absolute deviations of MESb from the median of MESb (MAD)
-#
-# COPYRIGHT: (C) 2015 Paulo van Breugel
-#            http://ecodiv.org
-#            http://pvanb.wordpress.com/
-#
-#            This program is free software under the GNU General Public
-#            License (>=v2). Read the file COPYING that comes with GRASS
-#            for details.
-#
-########################################################################
-#
-#%Module
-#% description: Compute the environmental bias
-#%End
-
-#%option
-#% key: env
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Raster map(s) of environmental conditions
-#% key_desc: names
-#% required: yes
-#% multiple: yes
-#% guisection: Input
-#%end
-
-#%option
-#% key: ref
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Area for which EB should be computed (binary map with 1 and 0)
-#% key_desc: names
-#% required: yes
-#% guisection: Input
-#%end
-
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Output MES layer (and root for IES layers if kept)
-#% key_desc: names
-#% required: no
-#% multiple: no
-#% guisection: Output
-#%end
-
-#%option
-#% key: file
-#% type: string
-#% gisprompt: new
-#% description: Name of output text file (csv format)
-#% key_desc: name
-#% required: no
-#% guisection: Output
-#%end
-
-#%flag
-#% key: i
-#% description: Compute EB individual variables (IES)
-#% guisection: Output
-#%end
-
-#%flag
-#% key: m
-#% description: Use minimum values of IES layers to compute MES (MES_min)
-#% guisection: Output
-#%end
-
-#%flag
-#% key: n
-#% description: Use mean values of IES layers to compute MES (MES_av)
-#% guisection: Output
-#%end
-
-#%flag
-#% key: o
-#% description: Use median values of IES layers to compute MES (MES_med)
-#% guisection: Output
-#%end
-
-#%option
-#% key: digits
-#% type: integer
-#% description: Precision of your input layers values
-#% key_desc: string
-#% answer: 5
-#%end
-
-#----------------------------------------------------------------------------
-# Standard
-#----------------------------------------------------------------------------
-
-# import libraries
-import os
-import sys
-import csv
-import numpy as np
-import uuid
-import operator
-import atexit
-import tempfile
-import string
-import grass.script as grass
-
-if not os.environ.has_key("GISBASE"):
-    grass.message("You must be in GRASS GIS to run this program.")
-    sys.exit(1)
-
-# create set to store names of temporary maps to be deleted upon exit
-clean_rast = set()
-
-def cleanup():
-    for rast in clean_rast:
-        grass.run_command("g.remove", type="rast", name=rast, quiet=True)
-
-#----------------------------------------------------------------------------
-# Functions
-#----------------------------------------------------------------------------
-
-# Create temporary name
-def tmpname(name):
-    tmpf = name + "_" + str(uuid.uuid4())
-    tmpf = string.replace(tmpf, '-', '_')
-    clean_rast.add(tmpf)
-    return tmpf
-
-# Create temporary filename
-def CreateFileName(outputfile):
-    flname = outputfile
-    k = 0
-    while os.path.isfile(flname):
-        k = k + 1
-        fn = flname.split('.')
-        if len(fn) == 1:
-            flname = fn[0] + "_" + str(k)
-        else:
-            flname = fn[0] + "_" + str(k) + "." + fn[1]
-    return flname
-
-def CheckLayer(envlay):
-    for chl in xrange(len(envlay)):
-        ffile = grass.find_file(envlay[chl], element = 'cell')
-        if ffile['fullname'] == '':
-            grass.fatal("The layer " + envlay[chl] + " does not exist.")
-
-# Write color rule file
-def defcol(mapname):
-    # Color table
-    tmpcol = tempfile.mkstemp()
-    text_file = open(tmpcol[1], "w")
-    text_file.write("0% 244:109:67\n")
-    text_file.write("50 255:255:210\n")
-    text_file.write("100% 50:136:189\n")
-    text_file.close()
-    # Assign color table
-    grass.run_command("r.colors", flags="n", quiet=True, map=mapname, rules=tmpcol[1])
-    os.remove(tmpcol[1])
-
-# Compute EB for input file (simlay = similarity, reflay = reference layer)
-def EB(simlay, reflay):
-
-    # layer name
-    vn = simlay.split("@")[0]
-
-    # Median and mad for whole region (within current mask)
-    tmpf4 = tmpname('reb4')
-    d = grass.read_command("r.quantile", quiet=True, input=simlay, percentiles="50")
-    d = d.split(":")
-    d = float(string.replace(d[2], "\n", ""))
-    grass.mapcalc("$tmpf4 = abs($map - $d)",
-                  map=simlay,
-                  tmpf4=tmpf4,
-                  d=d, quiet=True)
-    mad = grass.read_command("r.quantile", quiet=True, input=tmpf4, percentiles="50")
-    mad = mad.split(":")
-    mad = float(string.replace(mad[2], "\n", ""))
-    grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf4)
-
-    # Median and mad for reference layer
-    tmpf5 = tmpname('reb5')
-    clean_rast.add(tmpf5)
-    grass.run_command("r.null", map=reflay, setnull=0)
-    chlay = grass.parse_command("r.info", flags="r", map=reflay)
-    if chlay['max'] != '1' or chlay['min'] != '1':
-        grass.warning('please note that your layer is not a binary layer. Continuing assuming that all values other than 0 mark the area of interest')
-        grass.mapcalc("$tmpf5 = if($reflay != 0, 1, null())",
-                      tmpf5=tmpf5, reflay=reflay, quiet=True)
-    else:
-        grass.run_command("g.copy", raster=(reflay,tmpf5), quiet=True)
-    grass.mapcalc("$tmpf5 = $tmpf5 * $simlay", tmpf5=tmpf5,
-                  simlay=simlay, quiet=True, overwrite=True)
-    e = grass.read_command("r.quantile", quiet=True, input=tmpf5, percentiles="50")
-    e = e.split(":")
-    e = float(string.replace(e[2], "\n", ""))
-    EBstat = abs(d - e) / mad
-
-    # Print results to screen and return results
-    grass.info("Median " + vn + " (all region) = " + str(d))
-    grass.info("Median " + vn + " (ref. area) = " + str(e))
-    grass.info("MAD " + vn + " (all region) = " + str(mad))
-    grass.info("EB = " + str(EBstat))
-
-    # Clean up and return data
-    grass.run_command("g.remove", flags="f", type="raster", name=tmpf5, quiet=True)
-    return (mad, d, e, EBstat)
-
-
-def main():
-
-    #----------------------------------------------------------------------------
-    # Variables
-    #----------------------------------------------------------------------------
-
-#Test
-#options = {"env":"bio_1 at climate,bio_2 at climate,bio_3 at climate", "file":"test.txt", "ref":"PAs2", "output":"AA1", "digits":"5"}
-#flags = {"m":True, "n":True, "o":True, "i":True}
-
-    # variables
-    ipl = options['env']
-    ipl = ipl.split(',')
-    CheckLayer(ipl)
-    ipn = [z.split('@')[0] for z in ipl]
-    ipn = [x.lower() for x in ipn]
-    out = options['output']
-    if out == '':
-        tmpf0 = tmpname('reb0')
-    else:
-        tmpf0 = out
-    filename = options['file']
-    if filename != '':
-        filename = CreateFileName(filename)
-    ref = options['ref']
-    flag_m = flags['m']
-    flag_n = flags['n']
-    flag_o = flags['o']
-    flag_i = flags['i']
-    digits = int(options['digits'])
-    digits2 = pow(10, digits)
-
-    #----------------------------------------------------------------------------
-    # Compute MES
-    #----------------------------------------------------------------------------
-
-    ipi = []
-    for j in xrange(len(ipl)):
-
-        # Calculate the frequency distribution
-        tmpf1 = tmpname('reb1')
-        grass.mapcalc("$tmpf1 = int($dignum * $inplay)",
-                      tmpf1=tmpf1,
-                      inplay=ipl[j],
-                      dignum=digits2,
-                      quiet=True)
-        p = grass.pipe_command('r.stats', quiet=True, flags='cn', input=tmpf1, sort='asc', sep=';')
-        stval = {}
-        for line in p.stdout:
-            [val, count] = line.strip(os.linesep).split(';')
-            stval[float(val)] = float(count)
-        p.wait()
-        sstval = sorted(stval.items(), key=operator.itemgetter(0))
-        sstval = np.matrix(sstval)
-        a = np.cumsum(np.array(sstval), axis=0)
-        b = np.sum(np.array(sstval), axis=0)
-        c = a[:,1] / b[1] * 100
-
-        # Create recode rules
-        e1 = np.min(np.array(sstval), axis=0)[0] - 99999
-        e2 = np.max(np.array(sstval), axis=0)[0] - 99999
-        a1 = np.hstack([(e1), np.array(sstval.T[0])[0, :]])
-        a2 = np.hstack([np.array(sstval.T[0])[0,:] -1, (e2)])
-        b1 = np.hstack([(0), c])
-
-        tmprule = tempfile.mkstemp()
-        text_file = open(tmprule[1], "w")
-        for k in np.arange(0, len(b1.T)):
-            rtmp = str(int(a1[k])) + ":" + str(int(a2[k])) + ":" + str(b1[k])
-            text_file.write(rtmp + "\n")
-        text_file.close()
-
-        # Create the recode layer and calculate the IES
-        tmpf2 = tmpname('reb2')
-        grass.run_command("r.recode", input=tmpf1, output=tmpf2, rules=tmprule[1])
-
-        tmpf3 = tmpname('reb3')
-        z1 = tmpf3 + " = if(" + tmpf2 + "<=50, 2*float(" + tmpf2 + ")"
-        z3 = ", if(" + tmpf2 + "<100, 2*(100-float(" + tmpf2 + "))))"
-        calcc = z1 + z3
-        grass.mapcalc(calcc, quiet=True)
-        grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf2)
-        grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf1)
-        #os.remove(tmprule[1])
-        ipi.append(tmpf3)
-
-    #-----------------------------------------------------------------------
-    # Calculate EB statistics
-    #-----------------------------------------------------------------------
-
-    if flag_i or flag_m or flag_n or flag_o:
-
-        # EB MES
-        if flag_m:
-            nmn = tmpf0 + "_MES_mean"
-            grass.run_command("r.series", quiet=True, output=nmn, input=tuple(ipi), method="average")
-            defcol(nmn)
-            ebm = EB(simlay=nmn, reflay=ref)
-
-        if flag_n:
-            nmn = tmpf0 + "_MES_median"
-            grass.run_command("r.series", quiet=True, output=nmn, input=tuple(ipi), method="median")
-            defcol(nmn)
-            ebn = EB(simlay=nmn, reflay=ref)
-
-        if flag_o:
-            nmn = tmpf0 + "_MES_minimum"
-            grass.run_command("r.series", quiet=True, output=nmn, input=tuple(ipi), method="minimum")
-            defcol(nmn)
-            ebo = EB(simlay=nmn, reflay=ref)
-
-        # EB individual layers
-        if flag_i:
-            ebi = {}
-            for mm in xrange(len(ipi)):
-                nmn = tmpf0 + "_" + ipn[mm]
-                grass.run_command("g.rename", quiet=True, raster=(ipi[mm], nmn))
-                defcol(nmn)
-                value = EB(simlay=nmn, reflay=ref)
-                ebi[nmn] = value
-        else:
-            grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=ipi)
-    else:
-        grass.fatal("You need to select at least one of IES/MES_med/MES_av/MES_min options")
-
-    if filename != '':
-        with open(filename, 'wb') as csvfile:
-            fieldnames = ['variable', 'median_region', 'median_reference', 'mad', 'eb']
-            writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
-            writer.writeheader()
-            if flag_m:
-                writer.writerow({'variable':'MES_mean', 'median_region':ebm[0],
-                                'median_reference':ebm[1],'mad':ebm[2],'eb':ebm[3]})
-            if flag_n:
-                writer.writerow({'variable':'MES_mean','median_region':ebm[0],
-                    'median_reference':ebm[1], 'mad':ebm[2],'eb':ebm[3]})
-            if flag_n:
-                writer.writerow({'variable':'MES_median', 'median_region':ebn[0],
-                    'median_reference':ebn[1], 'mad':ebn[2],'eb':ebn[3]})
-            if flag_o:
-                writer.writerow({'variable':'MES_minimum','median_region':ebo[0],
-                    'median_reference':ebo[1], 'mad':ebo[2],'eb':ebo[3]})
-            if flag_i:
-                mykeys = ebi.keys()
-                for vari in mykeys:
-                    ebj = ebi[vari]
-                    writer.writerow({'variable':vari,'median_region':ebj[0],
-                        'median_reference':ebj[1], 'mad':ebj[2],'eb':ebj[3]})
-        grass.info("\nThe results are written to " + filename + "\n")
-        grass.info("\n-------------------------------------------\n")
-
-if __name__ == "__main__":
-    options, flags = grass.parser()
-    atexit.register(cleanup)
-    sys.exit(main())
-
-

Copied: grass-addons/grass7/raster/r.eb/r.meb.html (from rev 64090, grass-addons/grass7/raster/r.eb/r.eb.html)
===================================================================
--- grass-addons/grass7/raster/r.eb/r.meb.html	                        (rev 0)
+++ grass-addons/grass7/raster/r.eb/r.meb.html	2015-01-13 10:45:37 UTC (rev 64121)
@@ -0,0 +1,58 @@
+<h2>DESCRIPTION</h2>
+
+<p>The Multivariate Environmental Bias (MEB) takes the medium 
+conditions in an area <em>N</em> and computes how much conditions in 
+a subset of <em>N </em> (<em>S</em>) deviate from these medium 
+conditions. This can for example be used to see how well conditions 
+in the protected areas of a country represent conditions in the 
+whole country.
+
+<p>The messure is based on the Multivariate Environmental Similarity 
+(<em>MES</em>) surface, which was proposed by Elith et al (2010). To compute 
+the MES First the similarity of a point <em>P</em> to the conditions 
+in <em>N</em> with respect to variable <em>V</em> is computed. The 
+similarity is expressed as the deviation from the median of <em>V</em> 
+in <em>P</em> to those in <em>N</em>. This is done for all variables of 
+interest (V<sub>1</sub>, V<sub>2</sub>, ...V<sub>j</sub>).
+
+<p>In the original equation proposed by Elith et al (2010) the final 
+<em>MES</em> of <em>P</em> is computed as the minimum of the 
+similarity values (<em>IES<sub>minimum</sub></em>) of the individual 
+variables (<em>V<sub>j</sub></em>) in <em>P</em>. In <em>r.eb</em> 
+there is the option to use the mean (<em>IES<sub>mean</sub></em>) or 
+median (<em>IES<sub>median</sub></em>) of the <em>IES</em> instead. 
+For the MEB these may be a better choice as they take into account 
+the similarity along all environmnetal axes and not only the one 
+that deviates most. 
+
+<p>The <em>MEB</em> is computed as the absolute difference of the median 
+of the <em>MES</em> in the whole target area (MES<sub>n</sub>) and the 
+median of the <em>MES</em> in the subset (<em>MES<sub>s</sub></em>), 
+divided by the median absolute deviation 
+(<em><a href="http://en.wikipedia.org/wiki/Median_absolute_deviation">MAD</a></em>
+) of the <em>MES</em> in <em>N</em>. It is also possible to compute 
+the bias based on the individual variables (IEB) in which case the 
+<em>IES</em> instead of the <em>MES</em> is used. 
+
+<p><img src="r_eb_concept.png">
+
+<h2>NOTE</h2>
+Input variables are expected to be or to represent continuous variables. 
+
+<h2>REFERENCES</h2> 
+
+<p>Elith, J., Kearney, M., and Phillips, S. 
+2010. The art of modelling range-shifting species. Methods in 
+Ecology and Evolution 1:330-342.
+<p>van Breugel et al. Environmental gap analysis to prioritize 
+conservation efforts in eastern Africa. Submitted to Plos One
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.mess.html">r.mess</a></em>
+
+<h2>AUTHOR</h2>
+
+Paulo van Breugel, paulo at ecodiv.org
+
+<p><em>Last changed: $Date$</em>

Copied: grass-addons/grass7/raster/r.eb/r.meb.py (from rev 64074, grass-addons/grass7/raster/r.eb/r.eb.py)
===================================================================
--- grass-addons/grass7/raster/r.eb/r.meb.py	                        (rev 0)
+++ grass-addons/grass7/raster/r.eb/r.meb.py	2015-01-13 10:45:37 UTC (rev 64121)
@@ -0,0 +1,384 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+########################################################################
+#
+# MODULE:       r.meb
+# AUTHOR(S):    Paulo van Breugel <p.vanbreugel AT gmail.com>
+# PURPOSE:      Compute the multivariate envirionmental bias (EB). If A 
+#               is an areas within a larger region B, the EB represents 
+#               how much envirionmental conditions in A deviate from 
+#               median conditions in B. The first step is to compute the 
+#               multi-envirionmental similarity (MES) for B, using all 
+#               raster cells in B as reference points. The MES of a 
+#               raster cell thus represent how much conditions deviate
+#               from median conditions in B. The EB is then computed as the
+#               absolute difference of the median of MES values in A (MESa)
+#               and median of MES values in B (MESb), divided by the median of
+#               the absolute deviations of MESb from the median of MESb (MAD)
+#
+# COPYRIGHT: (C) 2015 Paulo van Breugel
+#            http://ecodiv.org
+#            http://pvanb.wordpress.com/
+#
+#            This program is free software under the GNU General Public
+#            License (>=v2). Read the file COPYING that comes with GRASS
+#            for details.
+#
+########################################################################
+#
+#%Module
+#% description: Compute the multivariate environmental bias (MEB)
+#%End
+
+#%option
+#% key: env
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Raster map(s) of environmental conditions
+#% key_desc: names
+#% required: yes
+#% multiple: yes
+#% guisection: Input
+#%end
+
+#%option
+#% key: ref
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Area for which EB should be computed (binary map with 1 and 0)
+#% key_desc: names
+#% required: yes
+#% guisection: Input
+#%end
+
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output MES layer (and root for IES layers if kept)
+#% key_desc: names
+#% required: no
+#% multiple: no
+#% guisection: Output
+#%end
+
+#%option
+#% key: file
+#% type: string
+#% gisprompt: new
+#% description: Name of output text file (csv format)
+#% key_desc: name
+#% required: no
+#% guisection: Output
+#%end
+
+#%flag
+#% key: i
+#% description: Compute EB for individual variables
+#% guisection: Output
+#%end
+
+#%flag
+#% key: m
+#% description: Use minimum values of IES layers to compute MES (MES_min)
+#% guisection: Output
+#%end
+
+#%flag
+#% key: n
+#% description: Use mean values of IES layers to compute MES (MES_av)
+#% guisection: Output
+#%end
+
+#%flag
+#% key: o
+#% description: Use median values of IES layers to compute MES (MES_med)
+#% guisection: Output
+#%end
+
+#%option
+#% key: digits
+#% type: integer
+#% description: Precision of your input layers values
+#% key_desc: string
+#% answer: 5
+#%end
+
+#----------------------------------------------------------------------------
+# Standard
+#----------------------------------------------------------------------------
+
+# import libraries
+import os
+import sys
+import csv
+import numpy as np
+import uuid
+import operator
+import atexit
+import tempfile
+import string
+import grass.script as grass
+
+if not os.environ.has_key("GISBASE"):
+    grass.message("You must be in GRASS GIS to run this program.")
+    sys.exit(1)
+
+# create set to store names of temporary maps to be deleted upon exit
+clean_rast = set()
+
+def cleanup():
+    for rast in clean_rast:
+        grass.run_command("g.remove", type="rast", name=rast, quiet=True)
+
+#----------------------------------------------------------------------------
+# Functions
+#----------------------------------------------------------------------------
+
+# Create temporary name
+def tmpname(name):
+    tmpf = name + "_" + str(uuid.uuid4())
+    tmpf = string.replace(tmpf, '-', '_')
+    clean_rast.add(tmpf)
+    return tmpf
+
+# Create temporary filename
+def CreateFileName(outputfile):
+    flname = outputfile
+    k = 0
+    while os.path.isfile(flname):
+        k = k + 1
+        fn = flname.split('.')
+        if len(fn) == 1:
+            flname = fn[0] + "_" + str(k)
+        else:
+            flname = fn[0] + "_" + str(k) + "." + fn[1]
+    return flname
+
+def CheckLayer(envlay):
+    for chl in xrange(len(envlay)):
+        ffile = grass.find_file(envlay[chl], element = 'cell')
+        if ffile['fullname'] == '':
+            grass.fatal("The layer " + envlay[chl] + " does not exist.")
+
+# Write color rule file
+def defcol(mapname):
+    # Color table
+    tmpcol = tempfile.mkstemp()
+    text_file = open(tmpcol[1], "w")
+    text_file.write("0% 244:109:67\n")
+    text_file.write("50 255:255:210\n")
+    text_file.write("100% 50:136:189\n")
+    text_file.close()
+    # Assign color table
+    grass.run_command("r.colors", flags="n", quiet=True, map=mapname, rules=tmpcol[1])
+    os.remove(tmpcol[1])
+
+# Compute EB for input file (simlay = similarity, reflay = reference layer)
+def EB(simlay, reflay):
+
+    # layer name
+    vn = simlay.split("@")[0]
+
+    # Median and mad for whole region (within current mask)
+    tmpf4 = tmpname('reb4')
+    d = grass.read_command("r.quantile", quiet=True, input=simlay, percentiles="50")
+    d = d.split(":")
+    d = float(string.replace(d[2], "\n", ""))
+    grass.mapcalc("$tmpf4 = abs($map - $d)",
+                  map=simlay,
+                  tmpf4=tmpf4,
+                  d=d, quiet=True)
+    mad = grass.read_command("r.quantile", quiet=True, input=tmpf4, percentiles="50")
+    mad = mad.split(":")
+    mad = float(string.replace(mad[2], "\n", ""))
+    grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf4)
+
+    # Median and mad for reference layer
+    tmpf5 = tmpname('reb5')
+    clean_rast.add(tmpf5)
+    grass.run_command("r.null", map=reflay, setnull=0)
+    chlay = grass.parse_command("r.info", flags="r", map=reflay)
+    if chlay['max'] != '1' or chlay['min'] != '1':
+        grass.warning('please note that your layer is not a binary layer. Continuing assuming that all values other than 0 mark the area of interest')
+        grass.mapcalc("$tmpf5 = if($reflay != 0, 1, null())",
+                      tmpf5=tmpf5, reflay=reflay, quiet=True)
+    else:
+        grass.run_command("g.copy", raster=(reflay,tmpf5), quiet=True)
+    grass.mapcalc("$tmpf5 = $tmpf5 * $simlay", tmpf5=tmpf5,
+                  simlay=simlay, quiet=True, overwrite=True)
+    e = grass.read_command("r.quantile", quiet=True, input=tmpf5, percentiles="50")
+    e = e.split(":")
+    e = float(string.replace(e[2], "\n", ""))
+    EBstat = abs(d - e) / mad
+
+    # Print results to screen and return results
+    grass.info("Median " + vn + " (all region) = " + str(d))
+    grass.info("Median " + vn + " (ref. area) = " + str(e))
+    grass.info("MAD " + vn + " (all region) = " + str(mad))
+    grass.info("EB = " + str(EBstat))
+
+    # Clean up and return data
+    grass.run_command("g.remove", flags="f", type="raster", name=tmpf5, quiet=True)
+    return (mad, d, e, EBstat)
+
+
+def main():
+
+    #----------------------------------------------------------------------------
+    # Variables
+    #----------------------------------------------------------------------------
+
+#Test
+#options = {"env":"bio_1 at climate,bio_2 at climate,bio_3 at climate", "file":"test.txt", "ref":"PAs2", "output":"AA1", "digits":"5"}
+#flags = {"m":True, "n":True, "o":True, "i":True}
+
+    # variables
+    ipl = options['env']
+    ipl = ipl.split(',')
+    CheckLayer(ipl)
+    ipn = [z.split('@')[0] for z in ipl]
+    ipn = [x.lower() for x in ipn]
+    out = options['output']
+    if out == '':
+        tmpf0 = tmpname('reb0')
+    else:
+        tmpf0 = out
+    filename = options['file']
+    if filename != '':
+        filename = CreateFileName(filename)
+    ref = options['ref']
+    flag_m = flags['m']
+    flag_n = flags['n']
+    flag_o = flags['o']
+    flag_i = flags['i']
+    digits = int(options['digits'])
+    digits2 = pow(10, digits)
+
+    #----------------------------------------------------------------------------
+    # Compute MES
+    #----------------------------------------------------------------------------
+
+    ipi = []
+    for j in xrange(len(ipl)):
+
+        # Calculate the frequency distribution
+        tmpf1 = tmpname('reb1')
+        # todo - check if layer is integer. If so, skip step below
+        grass.mapcalc("$tmpf1 = int($dignum * $inplay)",
+                      tmpf1=tmpf1,
+                      inplay=ipl[j],
+                      dignum=digits2,
+                      quiet=True)
+        p = grass.pipe_command('r.stats', quiet=True, flags='cn', input=tmpf1, sort='asc', sep=';')
+        stval = {}
+        for line in p.stdout:
+            [val, count] = line.strip(os.linesep).split(';')
+            stval[float(val)] = float(count)
+        p.wait()
+        sstval = sorted(stval.items(), key=operator.itemgetter(0))
+        sstval = np.matrix(sstval)
+        a = np.cumsum(np.array(sstval), axis=0)
+        b = np.sum(np.array(sstval), axis=0)
+        c = a[:,1] / b[1] * 100
+
+        # Create recode rules
+        e1 = np.min(np.array(sstval), axis=0)[0] - 99999
+        e2 = np.max(np.array(sstval), axis=0)[0] - 99999
+        a1 = np.hstack([(e1), np.array(sstval.T[0])[0, :]])
+        a2 = np.hstack([np.array(sstval.T[0])[0,:] -1, (e2)])
+        b1 = np.hstack([(0), c])
+
+        tmprule = tempfile.mkstemp()
+        text_file = open(tmprule[1], "w")
+        for k in np.arange(0, len(b1.T)):
+            rtmp = str(int(a1[k])) + ":" + str(int(a2[k])) + ":" + str(b1[k])
+            text_file.write(rtmp + "\n")
+        text_file.close()
+
+        # Create the recode layer and calculate the IES
+        tmpf2 = tmpname('reb2')
+        grass.run_command("r.recode", input=tmpf1, output=tmpf2, rules=tmprule[1])
+
+        tmpf3 = tmpname('reb3')
+        z1 = tmpf3 + " = if(" + tmpf2 + "<=50, 2*float(" + tmpf2 + ")"
+        z3 = ", if(" + tmpf2 + "<100, 2*(100-float(" + tmpf2 + "))))"
+        calcc = z1 + z3
+        grass.mapcalc(calcc, quiet=True)
+        grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf2)
+        grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf1)
+        #os.remove(tmprule[1])
+        ipi.append(tmpf3)
+
+    #-----------------------------------------------------------------------
+    # Calculate EB statistics
+    #-----------------------------------------------------------------------
+
+    if flag_i or flag_m or flag_n or flag_o:
+
+        # EB MES
+        if flag_m:
+            nmn = tmpf0 + "_MES_mean"
+            grass.run_command("r.series", quiet=True, output=nmn, input=tuple(ipi), method="average")
+            defcol(nmn)
+            ebm = EB(simlay=nmn, reflay=ref)
+
+        if flag_n:
+            nmn = tmpf0 + "_MES_median"
+            grass.run_command("r.series", quiet=True, output=nmn, input=tuple(ipi), method="median")
+            defcol(nmn)
+            ebn = EB(simlay=nmn, reflay=ref)
+
+        if flag_o:
+            nmn = tmpf0 + "_MES_minimum"
+            grass.run_command("r.series", quiet=True, output=nmn, input=tuple(ipi), method="minimum")
+            defcol(nmn)
+            ebo = EB(simlay=nmn, reflay=ref)
+
+        # EB individual layers
+        if flag_i:
+            ebi = {}
+            for mm in xrange(len(ipi)):
+                nmn = tmpf0 + "_" + ipn[mm]
+                grass.run_command("g.rename", quiet=True, raster=(ipi[mm], nmn))
+                defcol(nmn)
+                value = EB(simlay=nmn, reflay=ref)
+                ebi[nmn] = value
+        else:
+            grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=ipi)
+    else:
+        grass.fatal("You need to select at least one of IES/MES_med/MES_av/MES_min options")
+
+    if filename != '':
+        with open(filename, 'wb') as csvfile:
+            fieldnames = ['variable', 'median_region', 'median_reference', 'mad', 'eb']
+            writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
+            writer.writeheader()
+            if flag_m:
+                writer.writerow({'variable':'MES_mean', 'median_region':ebm[0],
+                                'median_reference':ebm[1],'mad':ebm[2],'eb':ebm[3]})
+            if flag_n:
+                writer.writerow({'variable':'MES_mean','median_region':ebm[0],
+                    'median_reference':ebm[1], 'mad':ebm[2],'eb':ebm[3]})
+            if flag_n:
+                writer.writerow({'variable':'MES_median', 'median_region':ebn[0],
+                    'median_reference':ebn[1], 'mad':ebn[2],'eb':ebn[3]})
+            if flag_o:
+                writer.writerow({'variable':'MES_minimum','median_region':ebo[0],
+                    'median_reference':ebo[1], 'mad':ebo[2],'eb':ebo[3]})
+            if flag_i:
+                mykeys = ebi.keys()
+                for vari in mykeys:
+                    ebj = ebi[vari]
+                    writer.writerow({'variable':vari,'median_region':ebj[0],
+                        'median_reference':ebj[1], 'mad':ebj[2],'eb':ebj[3]})
+        grass.info("\nThe results are written to " + filename + "\n")
+        grass.info("\n-------------------------------------------\n")
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    atexit.register(cleanup)
+    sys.exit(main())
+
+

Deleted: grass-addons/grass7/raster/r.eb/r_eb_concept.png
===================================================================
(Binary files differ)

Copied: grass-addons/grass7/raster/r.eb/r_meb_concept.png (from rev 64073, grass-addons/grass7/raster/r.eb/r_eb_concept.png)
===================================================================
(Binary files differ)



More information about the grass-commit mailing list