[GRASS-SVN] r68900 - grass-addons/grass7/raster/r.exdet
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jul 9 05:00:58 PDT 2016
Author: pvanbosgeo
Date: 2016-07-09 05:00:58 -0700 (Sat, 09 Jul 2016)
New Revision: 68900
Modified:
grass-addons/grass7/raster/r.exdet/r.exdet.html
grass-addons/grass7/raster/r.exdet/r.exdet.py
Log:
r.exdet changes to make code more compliant with pep8 guidelines, write metadata to layers
Modified: grass-addons/grass7/raster/r.exdet/r.exdet.html
===================================================================
--- grass-addons/grass7/raster/r.exdet/r.exdet.html 2016-07-09 12:00:49 UTC (rev 68899)
+++ grass-addons/grass7/raster/r.exdet/r.exdet.html 2016-07-09 12:00:58 UTC (rev 68900)
@@ -80,15 +80,10 @@
<h2>CITATION</h2>
-When using this tool, please use the citation to the original paper describing the method in your publications or other derived products:
+When using this tool, please use the citation to the original paper describing the method in your publications or other derived products.
<ul><li>Mesgaran, M.B., Cousens, R.D. and Webber, B.L. (2014) Here be dragons: a tool for quantifying novelty due to covariate range and correlation change when projecting species distribution models. Diversity & Distributions, 20: 1147-1159, DOI: 10.1111/ddi.12209.</li></ul>
-<p>If, in addition, you want to cite this addon, you can use:
-
-<ul><li>van Breugel, P. (2015) r.exdet. http://grasswiki.osgeo.org/wiki/AddOns/GRASS_7</li></ul>
-
-
<h2>REFERENCES</h2>
[1] Mesgaran, M.B., Cousens, R.D. and Webber, B.L. (2014) Here be dragons: a tool for quantifying novelty due to covariate range and correlation change when projecting species distribution models. Diversity & Distributions, 20: 1147-1159, DOI: 10.1111/ddi.12209.
Modified: grass-addons/grass7/raster/r.exdet/r.exdet.py
===================================================================
--- grass-addons/grass7/raster/r.exdet/r.exdet.py 2016-07-09 12:00:49 UTC (rev 68899)
+++ grass-addons/grass7/raster/r.exdet/r.exdet.py 2016-07-09 12:00:58 UTC (rev 68900)
@@ -12,31 +12,23 @@
# combinations between covariates (NT2).
# Based on Mesgaran et al.2014 [1]
#
-# COPYRIGHT: (C) 2015 Paulo van Breugel
-# http://ecodiv.org
-# http://pvanb.wordpress.com/
-# NOTE: Requires Numpy >= 1.8, GRASS >=7.0
+# COPYRIGHT: (C) 1997-2016 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.
##############################################################################
-# [1] Mesgaran, M.B., Cousens, R.D. & Webber, B.L. (2014) Here be dragons: a
-# tool for quantifying novelty due to covariate range and correlation change
-# when projecting species distribution models. Diversity & Distributions,
-# 20: 1147–1159, DOI: 10.1111/ddi.12209.
-
#%module
-#% description: Detection and quantification of novel uni- and multi-variate environments
+#% description: Quantification of novel uni- and multi-variate environments
+#% keyword: similarity
+#% keyword: multivariate
#% keyword: raster
#% keyword: modelling
#%end
-#%option
+#%option G_OPT_R_INPUTS
#% key: reference
-#% type: string
-#% gisprompt: old,cell,raster
#% description: Reference environmental raster layers
#% key_desc: raster
#% required: yes
@@ -44,10 +36,8 @@
#% guisection: Input
#%end
-#%option
+#%option G_OPT_R_INPUTS
#% key: projection
-#% type: string
-#% gisprompt: old,cell,raster
#% description: Projected environmental raster layers
#% key_desc: raster
#% required: yes
@@ -56,17 +46,6 @@
#%end
#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Output layer
-#% key_desc: raster
-#% required: yes
-#% multiple: no
-#% guisection: Output
-#%end
-
-#%option
#% key: region
#% type: string
#% gisprompt: new,region
@@ -77,6 +56,19 @@
#% guisection: Input
#%end
+#%rules
+#%exclusive: projection,region
+#%end
+
+#%option G_OPT_R_BASENAME_OUTPUT
+#% key: output
+#% description: Root name of the output layers
+#% key_desc: raster
+#% required: yes
+#% multiple: no
+#% guisection: Output
+#%end
+
#%flag
#% key: o
#% description: Most influential covariate (MIC) for NT1
@@ -95,110 +87,77 @@
#% description: Keep layer mahalanobis distance?
#%end
-#%flag
-#% key: m
-#% description: Restore mask?
-#%end
-
-#%flag
-#% key: r
-#% description: Restore region?
-#%end
-
# Dependencies
# Numpy >= 1.8, Scipy, grass >= 7.0
-#----------------------------------------------------------------------------
-# Standard
-#----------------------------------------------------------------------------
-
# import libraries
import os
import sys
import atexit
import numpy as np
-import grass.script as grass
+import grass.script as gs
import grass.script.array as garray
import tempfile
import uuid
import string
-if not os.environ.has_key("GISBASE"):
- grass.message("You must be in GRASS GIS to run this program.")
- sys.exit(1)
+COLORS_EXDET = """\
+0% 255:0:0
+0 255:200:200
+0.000001 200:255:200
+1.000001 200:200:255
+100% 0:0:255
+"""
-# 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
-#----------------------------------------------------------------------------
+CLEAN_LAY = []
-def tmpname(name):
- tmpf = name + "_" + str(uuid.uuid4())
- tmpf = string.replace(tmpf, '-', '_')
- clean_rast.add(tmpf)
- return tmpf
-def tmpregion():
- regionname = tmpname("mahal")
- k = 0
- ffile = grass.find_file(regionname, element = 'region',
- mapset=grass.gisenv()['MAPSET'])
- while ffile['fullname'] != '':
- regionname = regionname + k
- k = k + 1
- ffile = grass.find_file(regionname, element = 'region',
- mapset=grass.gisenv()['MAPSET'])
- grass.region(save=regionname)
- grass.run_command("g.region", save=regionname)
- return regionname
+def cleanup():
+ """Remove temporary maps specified in the global list"""
+ cleanrast = list(reversed(CLEAN_LAY))
+ for rast in cleanrast:
+ gs.run_command("g.remove", flags="f", type="all",
+ name=rast, quiet=True)
-def tmpmask():
- maskname = tmpname("mask")
- k = 0
- ffile = grass.find_file(maskname, element = 'cell',
- mapset=grass.gisenv()['MAPSET'])
- while ffile['fullname'] != '':
- maskname = maskname + k
- k = k + 1
- ffile = grass.find_file(maskname, element = 'cell',
- mapset=grass.gisenv()['MAPSET'])
- return maskname
def checkmask():
- ffile = grass.find_file("MASK", element = 'cell',
- mapset=grass.gisenv()['MAPSET'])
+ """Check if there is a MASK set"""
+ ffile = gs.find_file("MASK", element="CELL",
+ mapset=gs.gisenv()['MAPSET'])
mask_presence = ffile['fullname'] != ''
return mask_presence
-def icovar(input):
- tmpcov = tempfile.mkstemp()[1]
- text_file = open(tmpcov, "w")
- text_file.write(grass.read_command("r.covar", quiet=True, map=input))
- text_file.close()
- covar = np.loadtxt(tmpcov, skiprows=1)
- VI = np.linalg.inv(covar)
- os.remove(tmpcov)
- return VI
+def tmpname(prefix):
+ """Generate a tmp name which contains prefix
+ Store the name in the global list.
+ Use only for raster maps.
+ """
+ tmpf = prefix + str(uuid.uuid4())
+ tmpf = string.replace(tmpf, '-', '_')
+ CLEAN_LAY.append(tmpf)
+ return tmpf
+
+
def mahal(v, m, VI):
- delta = v - m[:,None,None]
- mahdist = np.sum(np.sum(delta[None,:,:,:] * VI[:,:,None,None],axis=1) * delta, axis=0)
+ """Compute the mahalanobis distance over reference layers"""
+ delta = v - m[:, None, None]
+ mahdist = np.sum(np.sum(delta[None, :, :, :] *
+ VI[:, :, None, None], axis=1) * delta, axis=0)
stat_mahal = garray.array()
stat_mahal[...] = mahdist
return stat_mahal
-def main():
- #----------------------------------------------------------------------------
- # Variables
- #----------------------------------------------------------------------------
+def main(options, flags):
+ gisbase = os.getenv('GISBASE')
+ if not gisbase:
+ gs.fatal(_('$GISBASE not defined'))
+ return 0
+
+ # Variables
ref = options['reference']
ref = ref.split(',')
pro = options['projection']
@@ -207,26 +166,33 @@
opn = [x.lower() for x in opn]
out = options['output']
region = options['region']
- flag_r = flags['r']
- flag_m = flags['m']
flag_d = flags['d']
flag_o = flags['o']
- # flag_p = flags['p'] # not implemented yet
+ #TODO: flag_p = flags['p']
- if ref == pro and checkmask() == False and region == '':
- grass.warning("1) reference and projected layers are the same, 2) you \
- did not define a mask, 3) and you did not define a region to delimit \
- the region for projection. Was that intentional?")
-
+ if region:
+ ffile = gs.find_file(region, element="region")
+ if not ffile:
+ gs.fatal(_("the region {} does not exist").format(region))
+ if ref == pro and checkmask() is False and not region:
+ gs.fatal(_("Reference and projected layers are the same, and "
+ "there are no mask or region defined. This means"
+ "that there is difference between the reference and "
+ "projection areas"))
if len(ref) != len(pro):
- grass.fatal("the number of reference and test layers should be the same")
+ gs.fatal(_("be the same. Your provided %d reference and %d"
+ "projection variables") % (len(ref), len(pro)))
# Create covariance table
- VI = icovar(input=ref)
+ tmpcov = tempfile.mkstemp()[1]
+ text_file = open(tmpcov, "w")
+ text_file.write(gs.read_command("r.covar", quiet=True, map=ref))
+ text_file.close()
+ covar = np.loadtxt(tmpcov, skiprows=1)
+ VI = np.linalg.inv(covar)
+ os.remove(tmpcov)
- #----------------------------------------------------------------------------
# Import reference data & compute univar stats per reference layer
- #----------------------------------------------------------------------------
s = len(ref)
dat_ref = stat_mean = stat_min = stat_max = None
layer = garray.array()
@@ -245,50 +211,27 @@
stat_min[i] = np.nanmin(layer)
stat_mean[i] = np.nanmean(layer)
stat_max[i] = np.nanmax(layer)
- dat_ref[i,:,:] = layer
+ dat_ref[i, :, :] = layer
del layer
- #----------------------------------------------------------------------------
- # Compute mahalanobis over reference layers
- #----------------------------------------------------------------------------
-
- # Compute mahalanobis over full set of layers
+ # Compute mahalanobis over full set of reference layers
mahal_ref = mahal(v=dat_ref, m=stat_mean, VI=VI)
mahal_ref_max = max(mahal_ref[np.isfinite(mahal_ref)])
del mahal_ref
- #----------------------------------------------------------------------------
- # Compute mahalanobis over projected layers
- #----------------------------------------------------------------------------
-
- # Remove mask (and optionally backup the MASK)
+ # Remove mask and set new region based on user-defined region or
+ # otherwise based on projection layers
if checkmask():
- if flag_m:
- mask_backup = tmpmask()
- grass.mapcalc("$mask = MASK", mask=mask_backup, quiet=True)
- grass.run_command("r.mask", flags="r")
+ gs.run_command("r.mask", flags="r")
+ if region:
+ gs.run_command("g.region", region=region)
+ gs.info(_("The region has set to the region {}").format(region))
+ if options['projection']:
+ gs.run_command("g.region", raster=pro[0])
+ # TODO: only set region to pro[0] when different from current region
+ gs.info(_("The region has set to match the proj raster layers"))
- # Make backup of region if required
- if flag_r:
- region_backup = tmpregion()
-
- # Set new region based on user-defined region or otherwise based on
- # projected layers
- if region != '':
- ffile = grass.find_file(region, element = 'region',
- mapset=grass.gisenv()['MAPSET'])
- if ffile != '':
- grass.run_command("g.region", region=region)
- grass.info("Region is set to the region. Please note that the function" + region)
- grass.info("does not check if projected layers fall within this region")
- else:
- grass.fatal("the region you selected does not exist")
- else:
- grass.run_command("g.region", raster=pro[0])
- if flag_r:
- grass.info("The original region was saved as " + region_backup)
-
- # Import projected layers
+ # Import projected layers in numpy array
s = len(pro)
dat_pro = None
layer = garray.array()
@@ -298,86 +241,61 @@
r, c = layer.shape
if dat_pro is None:
dat_pro = np.empty((s, r, c), dtype=np.double)
- dat_pro[i,:,:] = layer
+ dat_pro[i, :, :] = layer
del layer
# Compute mahalanobis distance
mahal_pro = mahal(v=dat_pro, m=stat_mean, VI=VI)
if flag_d:
- mapname = out + "_mahalanobis"
+ mapname = "{}_mahalanobis".format(out)
mahal_pro.write(mapname=mapname)
- grass.info("Mahalanobis distance map saved: " + mapname)
+ gs.info(_("Mahalanobis distance map saved: {}").format(mapname))
- #----------------------------------------------------------------------------
# Compute NT2
- #----------------------------------------------------------------------------
nt2map = garray.array()
nt2map[...] = mahal_pro / mahal_ref_max
- nt2 = out + "_NT2"
- clean_rast.add(nt2)
+ nt2 = "{}__NT2".format(out)
+ CLEAN_LAY.append(nt2)
nt2map.write(mapname=nt2)
- #----------------------------------------------------------------------------
# Compute NT1
- #----------------------------------------------------------------------------
-
tmplay = tmpname(out)
mnames = [None] * len(ref)
for i in xrange(len(ref)):
- tmpout = out + "_" + opn[i]
- clean_rast.add(tmpout)
- grass.mapcalc("$Dij = min(($prolay - $refmin), ($refmax - $prolay), 0) / ($refmax - $refmin)",
- Dij=tmpout,
- prolay=pro[i],
- refmin=stat_min[i],
- refmax=stat_max[i],
- quiet=True)
+ tmpout = "{}_{}".format(out, opn[i])
+ CLEAN_LAY.append(tmpout)
+ gs.mapcalc("$Dij = min(($prolay - $refmin), ($refmax - $prolay),"
+ "0) / ($refmax - $refmin)", Dij=tmpout, prolay=pro[i],
+ refmin=stat_min[i], refmax=stat_max[i], quiet=True)
mnames[i] = tmpout
- nt1 = out + "_NT1"
- grass.run_command("r.series", quiet=True, input=mnames, output=nt1, method="sum")
+ nt1 = "{}_NT1".format(out)
+ gs.run_command("r.series", quiet=True, input=mnames, output=nt1,
+ method="sum")
- #----------------------------------------------------------------------------
# Compute exdet novelty map
- #----------------------------------------------------------------------------
+ gs.mapcalc("$final = if($nt1<0,$nt1,$nt2)", final=out, nt1=nt1,
+ nt2=nt2, quiet=True, overwrite=True)
- grass.mapcalc("$final = if($nt1<0,$nt1,$nt2)",
- final=out,
- nt1=nt1,
- nt2=nt2,
- quiet=True, overwrite=True)
-
# Color table
- tmpcol = tempfile.mkstemp()
- text_file = open(tmpcol[1], "w")
- text_file.write("0% 255:0:0\n")
- text_file.write("0 255:200:200\n")
- text_file.write("0.000001 200:255:200\n")
- text_file.write("1 50:210:150\n")
- text_file.write("1.000001 200:200:255\n")
- text_file.write("100% 0:0:255\n")
- text_file.close()
- grass.run_command("r.colors", quiet=True, map=out, rules=tmpcol[1])
- os.remove(tmpcol[1])
+ gs.write_command("r.colors", map=out, rules='-', stdin=COLORS_MES,
+ quiet=True)
-#----------------------------------------------------------------------------
# Compute most influential covariate (MIC) metric for NT1
- #----------------------------------------------------------------------------
-
if flag_o:
- mod = out + "_MIC1"
- grass.run_command("r.series", quiet=True, output=mod, input=mnames, method="min_raster")
+ mod = "{}_MIC1".format(out)
+ gs.run_command("r.series", quiet=True, output=mod, input=mnames,
+ method="min_raster")
tmpcat = tempfile.mkstemp()
text_file = open(tmpcat[1], "w")
for cats in xrange(len(opn)):
- text_file.write(str(cats) + ":" + opn[cats] + "\n")
+ text_file.write("{}:{}\n".format(cats, opn[cats]))
text_file.close()
- grass.run_command("r.category", quiet=True, map=mod, rules=tmpcat[1], separator=":")
+ gs.run_command("r.category", quiet=True, map=mod, rules=tmpcat[1],
+ separator=":")
os.remove(tmpcat[1])
- #----------------------------------------------------------------------------
# Compute most influential covariate (MIC) metric for NT2
- # Not implemented yet
- #----------------------------------------------------------------------------
+ # TODO: Not implemented yet
# if flag_p:
# laylist = []
@@ -395,10 +313,10 @@
# # Combine ICp layers
# icpname = out + "_MIC2"
# icptemp = tmpname(name="icp")
-# grass.run_command("r.series", quiet=True, output=icptemp,
+# gs.run_command("r.series", quiet=True, output=icptemp,
# input=laylist, method="max_raster")
# # ICp values apply only to parts where NT2 > 1 and NT1 = 0
-# grass.mapcalc("$icpname = if($nt1<0,-888,if($nt2>1,$icptemp,-999))",
+# gs.mapcalc("$icpname = if($nt1<0,-888,if($nt2>1,$icptemp,-999))",
# nt1=nt1, nt2=nt2, icpname=icpname, icptemp=icptemp,
# quiet=True)
# # Write categories
@@ -409,34 +327,33 @@
# for cats in xrange(len(opn)):
# text_file.write(str(cats) + ":" + opn[cats] + "\n")
# text_file.close()
-# grass.run_command("r.category", quiet=True, map=icpname,
+# gs.run_command("r.category", quiet=True, map=icpname,
# rules=tmpcat[1], separator=":")
# os.remove(tmpcat[1])
-# grass.run_command("g.remove", quiet=True, type="raster", name=icptemp, flags="f")
+# gs.run_command("g.remove", quiet=True, type="raster", name=icptemp, flags="f")
#=======================================================================
## Clean up
#=======================================================================
- grass.run_command("g.remove", type="raster", flags="f",
+ gs.run_command("g.remove", type="raster", flags="f",
quiet=True, pattern=tmplay + "*")
- grass.run_command("g.remove", quiet=True, type="raster",
+ gs.run_command("g.remove", quiet=True, type="raster",
name=mnames, flags="f")
if flag_m:
- grass.run_command("g.rename", raster=[mask_backup,"MASK"])
- grass.info("Note, the original mask was restored")
+ gs.run_command("g.rename", raster=[mask_backup,"MASK"])
+ gs.info("Note, the original mask was restored")
- grass.info("Done....")
+ gs.info("Done....")
if flag_r:
- grass.info("Original region is restored")
+ gs.info("Original region is restored")
else:
- grass.info("Region is set to match the projected/test layers")
+ gs.info("Region is set to match the projected/test layers")
if __name__ == "__main__":
- options, flags = grass.parser()
atexit.register(cleanup)
- sys.exit(main())
+ sys.exit(main(*gs.parser()))
More information about the grass-commit
mailing list