Skip to content

Commit

Permalink
Update r.meb.py
Browse files Browse the repository at this point in the history
bugfixes
  • Loading branch information
ecodiv committed Feb 19, 2022
1 parent ba56135 commit 4962982
Showing 1 changed file with 30 additions and 29 deletions.
59 changes: 30 additions & 29 deletions src/raster/r.meb/r.meb.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
#!/usr/bin/env python


########################################################################
#
# MODULE: r.meb
Expand All @@ -19,7 +18,7 @@
# 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) 2014-2019 by Paulo van Breugel and the GRASS Development Team
# COPYRIGHT: (C) 2014-2022 by Paulo van Breugel and the GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
Expand Down Expand Up @@ -113,20 +112,16 @@
# import libraries
import os
import sys
import csv
import numpy as np
import uuid
import operator
import atexit
import uuid
import tempfile
import string
import operator
from subprocess import PIPE
import csv
import numpy as np
import grass.script as gs
from grass.pygrass.modules import Module

# for Python 3 compatibility
try:
xrange
except NameError:
xrange = range

# Rules
COLORS_MES = """\
Expand Down Expand Up @@ -158,14 +153,14 @@ def tmpname(prefix):
Use only for raster maps.
"""
tmpf = prefix + str(uuid.uuid4())
tmpf = string.replace(tmpf, "-", "_")
tmpf = tmpf.replace("-", "_")
CLEAN_RAST.append(tmpf)
return tmpf


def raster_exists(envlay):
"""Check if the raster map exists, call GRASS fatal otherwise"""
for chl in xrange(len(envlay)):
for chl in range(len(envlay)):
ffile = gs.find_file(envlay[chl], element="cell")
if not ffile["fullname"]:
gs.fatal(_("The layer {} does not exist").format(envlay[chl]))
Expand All @@ -179,11 +174,11 @@ def EB(simlay, reflay):
CLEAN_RAST.append(tmpf4)
d = gs.read_command("r.quantile", quiet=True, input=simlay, percentiles="50")
d = d.split(":")
d = float(string.replace(d[2], "\n", ""))
d = float(d[2].replace("\n", ""))
gs.mapcalc("$tmpf4 = abs($map - $d)", map=simlay, tmpf4=tmpf4, d=d, quiet=True)
mad = gs.read_command("r.quantile", quiet=True, input=tmpf4, percentiles="50")
mad = mad.split(":")
mad = float(string.replace(mad[2], "\n", ""))
mad = float(mad[2].replace("\n", ""))
gs.run_command("g.remove", quiet=True, flags="f", type="raster", name=tmpf4)

# Median and mad for reference layer
Expand All @@ -198,7 +193,7 @@ def EB(simlay, reflay):
)
e = gs.read_command("r.quantile", quiet=True, input=tmpf5, percentiles="50")
e = e.split(":")
e = float(string.replace(e[2], "\n", ""))
e = float(e[2].replace("\n", ""))
EBstat = abs(d - e) / mad

# Print results to screen and return results
Expand Down Expand Up @@ -255,8 +250,8 @@ def main(options, flags):
)

# Text for history in metadata
opt2 = dict((k, v) for k, v in options.iteritems() if v)
hist = " ".join("{!s}={!r}".format(k, v) for (k, v) in opt2.iteritems())
opt2 = dict((k, v) for k, v in options.items() if v)
hist = " ".join("{!s}={!r}".format(k, v) for (k, v) in opt2.items())
hist = "r.meb {}".format(hist)
unused, tmphist = tempfile.mkstemp()
text_file = open(tmphist, "w")
Expand All @@ -273,7 +268,7 @@ def main(options, flags):
gs.run_command("g.copy", quiet=True, raster=(ref, tmpref0))

ipi = []
for j in xrange(len(ipl)):
for j in range(len(ipl)):
# Calculate the frequency distribution
tmpf1 = tmpname("reb1")
CLEAN_RAST.append(tmpf1)
Expand All @@ -288,14 +283,20 @@ def main(options, flags):
dignum=digits2,
quiet=True,
)
p = gs.pipe_command(
"r.stats", quiet=True, flags="cn", input=tmpf1, sort="asc", sep=";"
)
stats_out = Module(
"r.stats",
flags="cn",
input=tmpf1,
sort="asc",
separator=";",
stdout_=PIPE,
).outputs.stdout
stval = {}
for line in p.stdout:
[val, count] = line.strip(os.linesep).split(";")
stats_outlines = stats_out.replace("\r", "").split("\n")
stats_outlines = [_f for _f in stats_outlines if _f]
for z in stats_outlines:
[val, count] = z.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)
Expand Down Expand Up @@ -332,7 +333,6 @@ def main(options, flags):
"g.remove", quiet=True, flags="f", type="raster", name=(tmpf2, tmpf1)
)
os.close(fd2)
os.remove(tmprule)
ipi.append(tmpf3)

# ----------------------------------------------------------------------
Expand Down Expand Up @@ -410,7 +410,7 @@ def main(options, flags):
# EB individual layers
if flag_i:
ebi = {}
for mm in xrange(len(ipi)):
for mm in range(len(ipi)):
nmn = "{}_{}".format(tmpf0, ipn[mm])
if not out:
CLEAN_RAST.append(nmn)
Expand All @@ -429,11 +429,12 @@ def main(options, flags):
description="Environmental similarity (ES) for " "{}".format(ipn[mm]),
loadhistory=tmphist,
)

else:
gs.run_command("g.remove", quiet=True, flags="f", type="raster", name=ipi)

if filename:
with open(filename, "wb") as csvfile:
with open(filename, "wt") as csvfile:
fieldnames = ["variable", "median_region", "median_reference", "mad", "eb"]
writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
writer.writeheader()
Expand Down

0 comments on commit 4962982

Please sign in to comment.