1#!/usr/local/bin/python3.8
2
3############################################################################
4#
5# MODULE:       i.oif
6# AUTHOR(S):    Markus Neteler 21.July 1998
7#               Updated for GRASS 5.7 by Michael Barton 2004/04/05
8#               Converted to Python by Glynn Clements
9#               Customised by Nikos Alexandris 04. August 2013
10#               Customised by Luca Delucchi Vienna Code Sprint 2014
11# PURPOSE:      calculates the Optimum Index factor of all band combinations
12#               for LANDSAT TM 1,2,3,4,5,7
13# COPYRIGHT:    (C) 1999,2008 by the GRASS Development Team
14#
15#               This program is free software under the GNU General Public
16#               License (>=v2). Read the file COPYING that comes with GRASS
17#               for details.
18#
19# Ref.: Jensen: Introductory digital image processing 1996, p.98
20#############################################################################
21
22#% Module
23#% description: Calculates Optimum-Index-Factor table for spectral bands
24#% keyword: imagery
25#% keyword: multispectral
26#% keyword: statistics
27#% End
28#% option G_OPT_R_INPUTS
29#% end
30#% option G_OPT_F_OUTPUT
31#% description: Name for output file (if omitted or "-" output to stdout)
32#% required: no
33#% end
34#% Flag
35#% key: g
36#% description: Print in shell script style
37#% End
38#% Flag
39#% key: s
40#% description: Process bands serially (default: run in parallel)
41#% End
42
43import sys
44import os
45from grass.script.utils import parse_key_val
46from grass.script import core as grass
47
48
49def oifcalc(sdev, corr, k1, k2, k3):
50    grass.debug(_("Calculating OIF for combination: %s, %s, %s" % (k1, k2,
51                                                                   k3)), 1)
52    # calculate SUM of Stddeviations:
53    ssdev = [sdev[k1], sdev[k2], sdev[k3]]
54    numer = sum(ssdev)
55
56    # calculate SUM of absolute(Correlation values):
57    scorr = [corr[k1, k2], corr[k1, k3], corr[k2, k3]]
58    denom = sum(map(abs, scorr))
59
60    # Calculate OIF index:
61    #     Divide (SUM of Stddeviations) and (SUM of Correlation)
62    return numer / denom
63
64
65def perms(bands):
66    n = len(bands)
67    for i in range(0, n - 2):
68        for j in range(i + 1, n - 1):
69            for k in range(j + 1, n):
70                yield (bands[i], bands[j], bands[k])
71
72
73def main():
74    shell = flags['g']
75    serial = flags['s']
76    bands = options['input'].split(',')
77
78    if len(bands) < 4:
79        grass.fatal(_("At least four input maps required"))
80
81    output = options['output']
82    # calculate the Stddev for TM bands
83    grass.message(_("Calculating standard deviations for all bands..."))
84    stddev = {}
85
86    if serial:
87        for band in bands:
88            grass.verbose("band %d" % band)
89            s = grass.read_command('r.univar', flags='g', map=band)
90            kv = parse_key_val(s)
91            stddev[band] = float(kv['stddev'])
92    else:
93        # run all bands in parallel
94        if "WORKERS" in os.environ:
95            workers = int(os.environ["WORKERS"])
96        else:
97            workers = len(bands)
98        proc = {}
99        pout = {}
100
101        # spawn jobs in the background
102        n = 0
103        for band in bands:
104            proc[band] = grass.pipe_command('r.univar', flags='g', map=band)
105            if n % workers is 0:
106                # wait for the ones launched so far to finish
107                for bandp in bands[:n]:
108                    if not proc[bandp].stdout.closed:
109                        pout[bandp] = proc[bandp].communicate()[0]
110                    proc[bandp].wait()
111            n = n + 1
112
113        # wait for jobs to finish, collect the output
114        for band in bands:
115            if not proc[band].stdout.closed:
116                pout[band] = proc[band].communicate()[0]
117            proc[band].wait()
118
119    # parse the results
120        for band in bands:
121            kv = parse_key_val(pout[band])
122            stddev[band] = float(kv['stddev'])
123
124    grass.message(_("Calculating Correlation Matrix..."))
125    correlation = {}
126    s = grass.read_command('r.covar', flags='r', map=[band for band in bands],
127                           quiet=True)
128
129    # We need to skip the first line, since r.covar prints the number of values
130    lines = s.splitlines()
131    for i, row in zip(bands, lines[1:]):
132        for j, cell in zip(bands, row.split(' ')):
133            correlation[i, j] = float(cell)
134
135    # Calculate all combinations
136    grass.message(_("Calculating OIF for all band combinations..."))
137
138    oif = []
139    for p in perms(bands):
140        oif.append((oifcalc(stddev, correlation, *p), p))
141    oif.sort(reverse=True)
142
143    grass.verbose(_("The Optimum Index Factor analysis result "
144                    "(best combination shown first):"))
145
146    if shell:
147        fmt = "%s,%s,%s:%.4f\n"
148    else:
149        fmt = "%s, %s, %s:  %.4f\n"
150
151    if not output or output == '-':
152        for v, p in oif:
153            sys.stdout.write(fmt % (p + (v,)))
154    else:
155        outf = open(output, 'w')
156        for v, p in oif:
157            outf.write(fmt % (p + (v,)))
158        outf.close()
159
160if __name__ == "__main__":
161    options, flags = grass.parser()
162    main()
163