[GRASS-SVN] r63987 - grass-addons/grass7/raster/r.mess
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jan 7 15:06:24 PST 2015
Author: pvanbosgeo
Date: 2015-01-07 15:06:24 -0800 (Wed, 07 Jan 2015)
New Revision: 63987
Modified:
grass-addons/grass7/raster/r.mess/r.mess.py
Log:
minor corrections
Modified: grass-addons/grass7/raster/r.mess/r.mess.py
===================================================================
--- grass-addons/grass7/raster/r.mess/r.mess.py 2015-01-07 20:53:11 UTC (rev 63986)
+++ grass-addons/grass7/raster/r.mess/r.mess.py 2015-01-07 23:06:24 UTC (rev 63987)
@@ -146,21 +146,21 @@
def cleanup():
for rast in clean_rast:
- grass.run_command("g.remove",
+ grass.run_command("g.remove",
type="rast", name = rast, quiet = True)
# main function
def main():
-
+
#----------------------------------------------------------------------------
# Variables
#----------------------------------------------------------------------------
-
+
#Test
#options = {"ref_rast":"bradypus", "ref_vect":"", "env_old":"bio1 at bradypus,bio5,bio6,bio7,bio8,bio12,bio16,bio17", "env_new":"bio1 at bradypus,bio5,bio6", "output":"MESSR", "digits":"5"}
#flags = {"m":True, "k":True, "n":False, "i":True, "k":True}
-
-
+
+
# reference layer
REF_VECT = options['ref_vect']
REF_RAST = options['ref_rast']
@@ -171,23 +171,23 @@
grass.run_command('g.message', flags="w", quiet=True,
message = 'Two reference layers defined (vector and raster). Using the raster layer!')
vtl = REF_RAST
- rtl = 'R'
+ rtl = 'R'
else:
if len(REF_RAST) > 0:
vtl = REF_RAST
- rtl = 'R'
+ rtl = 'R'
else:
vtl = REF_VECT
rtl = 'V'
-
+
# old environmental layers & variable names
- ipl = options['env_old']
+ ipl = options['env_old']
ipl = ipl.split(',')
ipn = [i.split('@')[0] for i in ipl]
ipn = [x.lower() for x in ipn]
-
+
# new environmental variables
- ipl2 = options['env_new']
+ ipl2 = options['env_new']
if ipl2 == '':
ipl_dif = False
ipl2 = ipl
@@ -196,12 +196,12 @@
ipl2 = ipl2.split(',')
if len(ipl2) != len(ipl) and len(ipl2) != 0:
grass.fatal('number of old and new environmental variables is not the same')
-
+
# output layers
opl = options['output']
opc = opl + '_MESS'
ipi = [opl + '_' + i for i in ipn]
-
+
# flags
flm = flags['m']
flk = flags['k']
@@ -209,50 +209,50 @@
il = flags['i']
flr = flags['r']
fll = flags['c']
-
+
# digits / precision
digits = int(options['digits'])
digits2 = pow(10, digits)
-
+
# Color table
tmpcol = tempfile.mkstemp()
text_file = open(tmpcol[1], "w")
text_file.write("0% 244:109:67\n")
text_file.write("0 255:255:210\n")
text_file.write("100% 50:136:189\n")
- text_file.close()
-
+ text_file.close()
+
# Check if there is a MASK
citiam = grass.find_file('MASK', element = 'cell')
-
+
#----------------------------------------------------------------------------
# Create the recode table - Reference distribution is raster
#----------------------------------------------------------------------------
-
+
if rtl=="R":
-
+
# Copy mask to temporary layer)
if citiam['fullname'] != '':
rname = "MASK" + str(uuid.uuid4())
rname = string.replace(rname, '-', '_')
grass.run_command('g.copy', quiet=True, raster = ('MASK', rname))
- clean_rast.add(rname)
-
+ clean_rast.add(rname)
+
for i in xrange(len(ipl)):
-
+
# Create temporary layer based on reference layer
tmpf0 = "rmess_tmp_" + str(uuid.uuid4())
tmpf0 = string.replace(tmpf0, '-', '_')
grass.mapcalc("$tmpf0 = $vtl * 1", vtl = vtl, tmpf0=tmpf0, quiet=True)
clean_rast.add(tmpf0)
-
+
# Remove mask
if citiam['fullname'] != '':
grass.run_command("r.mask", quiet=True, flags="r")
-
+
# Create mask based on combined MASK/input layer
grass.run_command("r.mask", quiet=True, raster=tmpf0)
-
+
# Calculate the frequency distribution
tmpf1 = "rmess_tmp2_" + str(uuid.uuid4())
tmpf1 = string.replace(tmpf1, '-', '_')
@@ -273,7 +273,7 @@
a = np.cumsum(np.array(sstval), axis=0)
b = np.sum(np.array(sstval), axis=0)
c = a[:,1] / b[1] * 100
-
+
# Restore mask and set region to new env_new
# Note: if region env_new is different from region env_old, setting
# the mask will not do anything. If there is partial overlap, the mask
@@ -285,7 +285,7 @@
grass.run_command("r.mask", quiet=True, flags="r")
if ipl_dif and not flr:
grass.run_command("g.region", quiet=True, raster=ipl2[0])
-
+
# Get min and max values for recode table (based on full map)
tmpf2 = "rmess_tmp2_" + str(uuid.uuid4())
tmpf2 = string.replace(tmpf2, '-', '_')
@@ -294,88 +294,88 @@
inplay = ipl2[i],
dignum = digits2)
d = grass.parse_command("r.univar",quiet=True, flags="g", map=tmpf2)
-
+
# Create recode rules
Dmin = int(d['min'])
Dmax = int(d['max'])
envmin = np.min(np.array(sstval), axis=0)[0]
envmax = np.max(np.array(sstval), axis=0)[0]
-
+
if Dmin < envmin: e1 = Dmin - 1
else: e1 = envmin -1
if Dmax > envmax: e2 = Dmax + 1
else: e2 = envmax + 1
-
+
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()
-
+ text_file.close()
+
# Create the recode layer and calculate the IES
tmpf3 = "rmess_tmp3_" + str(uuid.uuid4())
tmpf3 = string.replace(tmpf3, '-', '_')
grass.run_command("r.recode", input=tmpf2, output=tmpf3, rules=tmprule[1])
-
- z1 = ipi[i] + " = if(" + tmpf3 + "==0, (float(" + tmpf2 + ")-" + str(float(envmin)) + ")/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100"
+
+ z1 = ipi[i] + " = if(" + tmpf3 + "==0, (float(" + tmpf2 + ")-" + str(float(envmin)) + ")/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100"
z2 = ", if(" + tmpf3 + "<=50, 2*float(" + tmpf3 + ")"
- z3 = ", if(" + tmpf3 + "<100, 2*(100-float(" + tmpf3 + "))"
+ z3 = ", if(" + tmpf3 + "<100, 2*(100-float(" + tmpf3 + "))"
z4 = ", (" + str(float(envmax)) + "- float(" + tmpf2 + "))/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100.0)))"
- calcc = z1 + z2 + z3 + z4
+ calcc = z1 + z2 + z3 + z4
grass.mapcalc(calcc, quiet=True)
grass.run_command("r.colors", quiet=True, map=ipi[i], rules=tmpcol[1])
grass.run_command("g.remove", quiet=True, flags="f", type="raster", pattern=[tmpf0,tmpf1,tmpf2,tmpf3])
os.remove(tmprule[1])
-
+
# Recover original mask
if citiam['fullname'] != '':
grass.run_command("r.mask", quiet=True, raster=rname)
grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=rname)
-
-
+
+
#----------------------------------------------------------------------------
# Create the recode table - Reference distribution is vector
#----------------------------------------------------------------------------
-
+
if rtl=="V":
-
+
# Copy point layer and add columns for variables
tmpf0 = "rmess_tmp_" + str(uuid.uuid4())
tmpf0 = string.replace(tmpf0, '-', '_')
clean_rast.add(tmpf0)
-
+
grass.run_command("v.extract", quiet=True, flags="t", input=vtl, type="point", output=tmpf0)
- grass.run_command("v.db.addtable", quiet=True, map=tmpf0)
+ grass.run_command("v.db.addtable", quiet=True, map=tmpf0)
grass.run_command("v.db.addcolumn", quiet=True, map=tmpf0, columns="envvar double precision")
-
-
+
+
# Upload raster values and get value in python as frequency table
check_n = len(np.hstack(db.db_select(sql = "SELECT cat FROM " + tmpf0)))
for m in xrange(len(ipl)):
-
+
grass.run_command("db.execute" , quiet=True, sql = "UPDATE " + tmpf0 + " SET envvar = NULL")
grass.run_command("v.what.rast", quiet=True, map=tmpf0, layer=1, raster=ipl[m], column="envvar")
volval = np.vstack(db.db_select(sql = "SELECT envvar,count(envvar) from " + tmpf0 +
" WHERE envvar IS NOT NULL GROUP BY envvar ORDER BY envvar"))
- volval = volval.astype(np.float, copy=False)
+ volval = volval.astype(np.float, copy=False)
a = np.cumsum(volval[:,1], axis=0)
b = np.sum(volval[:,1], axis=0)
c = a / b * 100
-
+
# Check for point without values
if b < check_n:
- print("Please note that there were " + str(check_n - b) + " points without value")
- print("This is probably because they lies outside the current region or input rasters")
-
+ grass.info("Please note that there were " + str(check_n - b) + " points without value")
+ grass.info("This is probably because they lies outside the current region or input rasters")
+
# Set region to env_new layers (if different from env_old)
if ipl_dif and not flr:
- grass.run_command("g.region", quiet=True, raster=ipl2[0])
-
+ grass.run_command("g.region", quiet=True, raster=ipl2[0])
+
# Multiply env layer with dignum
tmpf2 = "rmess_tmp2_" + str(uuid.uuid4())
tmpf2 = string.replace(tmpf2, '-', '_')
@@ -384,66 +384,66 @@
inplay = ipl2[m],
dignum = digits2,
quiet=True)
-
+
# Calculate min and max values of sample points and raster layer
envmin = int(min(volval[:,0]) * digits2)
envmax = int(max(volval[:,0]) * digits2)
Drange = grass.read_command("r.info", flags="r", map=tmpf2)
Drange = Drange.split('\n')
Drange = np.hstack([i.split('=') for i in Drange])
- Dmin = int(Drange[1])
- Dmax = int(Drange[3])
-
+ Dmin = int(Drange[1])
+ Dmax = int(Drange[3])
+
if Dmin < envmin: e1 = Dmin - 1
else: e1 = envmin -1
if Dmax > envmax: e2 = Dmax + 1
else: e2 = envmax + 1
-
+
a0 = volval[:,0] * digits2
- a0 = a0.astype(np.int, copy=False)
+ a0 = a0.astype(np.int, copy=False)
a1 = np.hstack([(e1) , a0 ])
a2 = np.hstack([a0 -1, (e2) ])
b1 = np.hstack([(0), c])
-
+
tmprule = tempfile.mkstemp(suffix=ipn[m])
text_file = open(tmprule[1], "w")
for k in np.arange(0,len(b1)):
rtmp = str(int(a1[k])) + ":" + str(int(a2[k])) + ":" + str(b1[k])
text_file.write(rtmp + "\n")
- text_file.close()
-
+ text_file.close()
+
# Create the recode layer and calculate the IES
tmpf3 = "rmess_tmp3_" + str(uuid.uuid4())
tmpf3 = string.replace(tmpf3, '-', '_')
grass.run_command("r.recode", quiet=True, input=tmpf2, output=tmpf3, rules=tmprule[1])
-
- z1 = ipi[m] + " = if(" + tmpf3 + "==0, (float(" + tmpf2 + ")-" + str(float(envmin)) + ")/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100"
+
+ z1 = ipi[m] + " = if(" + tmpf3 + "==0, (float(" + tmpf2 + ")-" + str(float(envmin)) + ")/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100"
z2 = ", if(" + tmpf3 + "<=50, 2*float(" + tmpf3 + ")"
- z3 = ", if(" + tmpf3 + "<100, 2*(100-float(" + tmpf3 + "))"
+ z3 = ", if(" + tmpf3 + "<100, 2*(100-float(" + tmpf3 + "))"
z4 = ", (" + str(float(envmax)) + "- float(" + tmpf2 + "))/(" + str(float(envmax)) + "-" + str(float(envmin)) + ") * 100.0)))"
- calcc = z1 + z2 + z3 + z4
+ calcc = z1 + z2 + z3 + z4
grass.mapcalc(calcc, quiet=True)
grass.run_command("r.colors", quiet=True, map=ipi[m], rules=tmpcol[1])
-
- # Clean up
+
+ # Clean up
grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=[tmpf2,tmpf3])
os.remove(tmprule[1])
-
+
# Clean up
grass.run_command("g.remove", quiet=True, flags="f", type="vector", name=tmpf0)
-
-
+
+
#-----------------------------------------------------------------------
# Calculate MESS statistics
#-----------------------------------------------------------------------
-
+
# MESS
grass.run_command("r.series", quiet=True, output=opc, input=ipi, method="minimum")
-
+
# Area with negative MESS
if fln:
grass.mapcalc(opc + "_neg = if(" + opc + "<0,1,0)", quiet=True)
-
+
# Most dissimilar variable (MoD)
if flm:
mod = opc + "_MoD"
@@ -454,7 +454,7 @@
text_file.write(str(cats) + ":" + ipi[cats] + "\n")
text_file.close()
grass.run_command("r.category", quiet=True, map=mod, rules=tmpcat[1], separator=":")
-
+
# mean(IES), where IES < 0
if flk:
MinMes = grass.read_command("r.info", quiet=True, flags="r", map=opc)
@@ -463,7 +463,7 @@
mod = opc + "MeanNeg"
grass.run_command("r.series", quiet=True, input=ipi, output=mod, method="average", range=(MinMes,0))
grass.run_command("r.colors", quiet=True, map=mod, col="ryg")
-
+
# Number of layers with negative values
if fll:
MinMes = grass.read_command("r.info", quiet=True, flags="r", map=opc)
@@ -472,18 +472,18 @@
mod = opc + "CountNeg"
grass.run_command("r.series", flags="n", quiet=True, input=ipi, output=mod, method="count", range=(MinMes,0))
grass.run_command("r.colors", quiet=True, map=mod, col="ryg")
-
+
# Remove IES layers
if il:
grass.run_command("g.remove", quiet=True, flags="f", type="raster", name=ipi)
-
+
#=======================================================================
## Clean up
#=======================================================================
-
+
os.remove(tmpcol[1])
os.remove(tmpcat[1])
-
+
if __name__ == "__main__":
options, flags = grass.parser()
atexit.register(cleanup)
More information about the grass-commit
mailing list