1#!/usr/local/bin/python3.8 2 3############################################################################ 4# 5# MODULE: i.spectral 6# AUTHOR(S): Markus Neteler, 18. August 1998 7# Converted to Python by Glynn Clements 8# PURPOSE: Displays spectral response at user specified locations in 9# group or raster images 10# COPYRIGHT: (C) 1999-2013 by the GRASS Development Team 11# 12# This program is free software under the GNU General Public 13# License (>=v2). Read the file COPYING that comes with GRASS 14# for details. 15# 16############################################################################# 17# 18# this script needs gnuplot for pretty rendering 19# TODO: use PyPlot like the wxGUI Profiling tool 20# 21# written by Markus Neteler 18. August 1998 22# neteler geog.uni-hannover.de 23# 24# bugfix: 25. Nov.98/20. Jan. 1999 25# 3 March 2006: Added multiple images and group support by Francesco Pirotti - CIRGEO 26# 27 28#%Module 29#% description: Displays spectral response at user specified locations in group or images. 30#% keyword: imagery 31#% keyword: querying 32#% keyword: raster 33#% keyword: multispectral 34#%End 35#%option G_OPT_I_GROUP 36#% required : no 37#% guisection: Input 38#%end 39#%option G_OPT_I_SUBGROUP 40#% required : no 41#% guisection: Input 42#%end 43#%option G_OPT_R_INPUTS 44#% key: raster 45#% required : no 46#% guisection: Input 47#%end 48#%option G_OPT_M_COORDS 49#% multiple: yes 50#% required: yes 51#% guisection: Input 52#%end 53#%option G_OPT_F_OUTPUT 54#% key: output 55#% description: Name for output image (or text file for -t) 56#% guisection: Output 57#% required : no 58#%end 59#%Option 60#% key: format 61#% type: string 62#% description: Graphics format for output file 63#% options: png,eps,svg 64#% answer: png 65#% multiple: no 66#% guisection: Output 67#%End 68#%flag 69#% key: c 70#% description: Show sampling coordinates instead of numbering in the legend 71#%end 72#% flag 73#% key: g 74#% description: Use gnuplot for display 75#%end 76#% flag 77#% key: t 78#% description: output to text file 79#%end 80 81import os 82import atexit 83from grass.script.utils import try_rmdir 84from grass.script import core as gcore 85 86 87def cleanup(): 88 try_rmdir(tmp_dir) 89 90 91def write2textf(what, output): 92 outf = open(output, 'w') 93 i = 0 94 for row in enumerate(what): 95 i = i + 1 96 outf.write("%d, %s\n" % (i, row)) 97 outf.close() 98 99 100def draw_gnuplot(what, xlabels, output, img_format, coord_legend): 101 xrange = 0 102 103 for i, row in enumerate(what): 104 outfile = os.path.join(tmp_dir, 'data_%d' % i) 105 outf = open(outfile, 'w') 106 xrange = max(xrange, len(row) - 2) 107 for j, val in enumerate(row[3:]): 108 outf.write("%d %s\n" % (j + 1, val)) 109 outf.close() 110 111 # build gnuplot script 112 lines = [] 113 if output: 114 if img_format == 'png': 115 term_opts = "png truecolor large size 825,550" 116 elif img_format == 'eps': 117 term_opts = "postscript eps color solid size 6,4" 118 elif img_format == 'svg': 119 term_opts = "svg size 825,550 dynamic solid" 120 else: 121 gcore.fatal(_("Programmer error (%s)") % img_format) 122 123 lines += [ 124 "set term " + term_opts, 125 "set output '%s'" % output 126 ] 127 128 lines += [ 129 "set xtics (%s)" % xlabels, 130 "set grid", 131 "set title 'Spectral signatures'", 132 "set xrange [0.5 : %d - 0.5]" % xrange, 133 "set noclabel", 134 "set xlabel 'Bands'", 135 "set xtics rotate by -40", 136 "set ylabel 'DN Value'", 137 "set style data lines" 138 ] 139 140 cmd = [] 141 for i, row in enumerate(what): 142 if not coord_legend: 143 title = 'Pick ' + str(i + 1) 144 else: 145 title = str(tuple(row[0:2])) 146 147 x_datafile = os.path.join(tmp_dir, 'data_%d' % i) 148 cmd.append(" '%s' title '%s'" % (x_datafile, title)) 149 150 cmd = ','.join(cmd) 151 cmd = ' '.join(['plot', cmd, "with linespoints pt 779"]) 152 lines.append(cmd) 153 154 plotfile = os.path.join(tmp_dir, 'spectrum.gnuplot') 155 plotf = open(plotfile, 'w') 156 for line in lines: 157 plotf.write(line + '\n') 158 plotf.close() 159 160 if output: 161 gcore.call(['gnuplot', plotfile]) 162 else: 163 gcore.call(['gnuplot', '-persist', plotfile]) 164 165 166def draw_linegraph(what): 167 yfiles = [] 168 169 xfile = os.path.join(tmp_dir, 'data_x') 170 171 xf = open(xfile, 'w') 172 for j, val in enumerate(what[0][3:]): 173 xf.write("%d\n" % (j + 1)) 174 xf.close() 175 176 for i, row in enumerate(what): 177 yfile = os.path.join(tmp_dir, 'data_y_%d' % i) 178 yf = open(yfile, 'w') 179 for j, val in enumerate(row[3:]): 180 yf.write("%s\n" % val) 181 yf.close() 182 yfiles.append(yfile) 183 184 sienna = '#%02x%02x%02x' % (160, 82, 45) 185 coral = '#%02x%02x%02x' % (255, 127, 80) 186 gp_colors = ['red', 'green', 'blue', 'magenta', 'cyan', sienna, 'orange', 187 coral] 188 189 colors = gp_colors 190 while len(what) > len(colors): 191 colors += gp_colors 192 colors = colors[0:len(what)] 193 194 gcore.run_command('d.linegraph', x_file=xfile, y_file=yfiles, 195 y_color=colors, title=_("Spectral signatures"), 196 x_title=_("Bands"), y_title=_("DN Value")) 197 198 199def main(): 200 group = options['group'] 201 subgroup = options['subgroup'] 202 raster = options['raster'] 203 output = options['output'] 204 coords = options['coordinates'] 205 img_fmt = options['format'] 206 coord_legend = flags['c'] 207 gnuplot = flags['g'] 208 textfile = flags['t'] 209 210 global tmp_dir 211 tmp_dir = gcore.tempdir() 212 213 if not group and not raster: 214 gcore.fatal(_("Either group= or raster= is required")) 215 216 if group and raster: 217 gcore.fatal(_("group= and raster= are mutually exclusive")) 218 219 # -t needs an output filename 220 if textfile and not output: 221 gcore.fatal(_("Writing to text file requires output=filename")) 222 223 # check if gnuplot is present 224 if gnuplot and not gcore.find_program('gnuplot', '-V'): 225 gcore.fatal(_("gnuplot required, please install first")) 226 227 # get data from group listing and set the x-axis labels 228 if group: 229 # Parse the group list output 230 if subgroup: 231 s = gcore.read_command('i.group', flags='g', group=group, subgroup=subgroup, quiet=True) 232 else: 233 s = gcore.read_command('i.group', flags='g', group=group, quiet=True) 234 rastermaps = s.splitlines() 235 else: 236 # get data from list of files and set the x-axis labels 237 rastermaps = raster.split(',') 238 239 xlabels = ["'%s' %d" % (n, i + 1) for i, n in enumerate(rastermaps)] 240 xlabels = ','.join(xlabels) 241 242 # get y-data for gnuplot-data file 243 what = [] 244 s = gcore.read_command('r.what', map=rastermaps, coordinates=coords, 245 null='0', quiet=True) 246 if len(s) == 0: 247 gcore.fatal(_('No data returned from query')) 248 249 for l in s.splitlines(): 250 f = l.split('|') 251 for i, v in enumerate(f): 252 if v in ['', '*']: 253 f[i] = 0 254 else: 255 f[i] = float(v) 256 what.append(f) 257 258 # build data files 259 if gnuplot: 260 draw_gnuplot(what, xlabels, output, img_fmt, coord_legend) 261 elif textfile: 262 write2textf(what, output) 263 else: 264 draw_linegraph(what) 265 266if __name__ == "__main__": 267 options, flags = gcore.parser() 268 atexit.register(cleanup) 269 main() 270