[GRASS-SVN] r64122 - in grass-addons/grass7/raster: . r.meb
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 13 02:57:07 PST 2015
Author: pvanbosgeo
Date: 2015-01-13 02:57:07 -0800 (Tue, 13 Jan 2015)
New Revision: 64122
Added:
grass-addons/grass7/raster/r.meb/
grass-addons/grass7/raster/r.meb/Makefile
grass-addons/grass7/raster/r.meb/r.meb.html
grass-addons/grass7/raster/r.meb/r.meb.py
grass-addons/grass7/raster/r.meb/r_meb_concept.png
Removed:
grass-addons/grass7/raster/r.meb/r.eb.html
grass-addons/grass7/raster/r.meb/r.eb.py
grass-addons/grass7/raster/r.meb/r_eb_concept.png
Log:
Changed name to r.meb to avoid confusion with i.eb set of functions
Copied: grass-addons/grass7/raster/r.meb/Makefile (from rev 64121, grass-addons/grass7/raster/r.eb/Makefile)
===================================================================
--- grass-addons/grass7/raster/r.meb/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.meb/Makefile 2015-01-13 10:57:07 UTC (rev 64122)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.meb
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Deleted: grass-addons/grass7/raster/r.meb/r.eb.html
===================================================================
--- grass-addons/grass7/raster/r.eb/r.eb.html 2015-01-12 11:16:12 UTC (rev 64073)
+++ grass-addons/grass7/raster/r.meb/r.eb.html 2015-01-13 10:57:07 UTC (rev 64122)
@@ -1,51 +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. <br><br> <img src=
-"r_eb_concept.png" alt="EB concept"> <br><br> <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.meb/r.eb.py
===================================================================
--- grass-addons/grass7/raster/r.eb/r.eb.py 2015-01-12 11:16:12 UTC (rev 64073)
+++ grass-addons/grass7/raster/r.meb/r.eb.py 2015-01-13 10:57:07 UTC (rev 64122)
@@ -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.meb/r.meb.html (from rev 64121, grass-addons/grass7/raster/r.eb/r.meb.html)
===================================================================
--- grass-addons/grass7/raster/r.meb/r.meb.html (rev 0)
+++ grass-addons/grass7/raster/r.meb/r.meb.html 2015-01-13 10:57:07 UTC (rev 64122)
@@ -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.meb/r.meb.py (from rev 64121, grass-addons/grass7/raster/r.eb/r.meb.py)
===================================================================
--- grass-addons/grass7/raster/r.meb/r.meb.py (rev 0)
+++ grass-addons/grass7/raster/r.meb/r.meb.py 2015-01-13 10:57:07 UTC (rev 64122)
@@ -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.meb/r_eb_concept.png
===================================================================
(Binary files differ)
Copied: grass-addons/grass7/raster/r.meb/r_meb_concept.png (from rev 64121, grass-addons/grass7/raster/r.eb/r_meb_concept.png)
===================================================================
(Binary files differ)
More information about the grass-commit
mailing list