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