#!/usr/bin/env python ############################################################################ # # MODULE: i.spectral # AUTHOR(S): Markus Neteler, 18. August 1998 # Converted to Python by Glynn Clements # PURPOSE: displays spectral response at user specified locations in # group or images # COPYRIGHT: (C) 1999 by the GRASS Development Team # # This program is free software under the GNU General Public # License (>=v2). Read the file COPYING that comes with GRASS # for details. # ############################################################################# # written by Markus Neteler 18. August 1998 # neteler geog.uni-hannover.de # # bugfix: 25. Nov.98/20. Jan. 1999 # 3 March 2006: Added multiple images and group support by Francesco Pirotti - CIRGEO # this script needs gnuplot #%Module #% description: Displays spectral response at user specified locations in group or images. #% keywords: imagery #% keywords: querying #% keywords: raster #% keywords: multispectral #%End #%option G_OPT_I_GROUP #% required : no #%end #%option G_OPT_R_INPUTS #% key: raster #% required : no #%end #%option G_OPT_F_OUTPUT #% key: output #% description: Name for output PNG image #% required : no #%end #%option G_OPT_M_COORDS #% multiple: yes #% required: yes #%end #%flag #%key: c #%description: Label with coordinates instead of numbering #%end #%flag #%key: g #%description: Use gnuplot for display #%end import os import atexit from grass.script import core as grass def cleanup(): grass.try_rmdir(tmp_dir) def draw_gnuplot(what, xlabels, output, label): xrange = 0 for i, row in enumerate(what): outfile = os.path.join(tmp_dir, 'data_%d' % i) outf = open(outfile, 'w') xrange = max(xrange, len(row) - 2) for j, val in enumerate(row[3:]): outf.write("%d %s\n" % (j + 1, val)) outf.close() # build gnuplot script lines = [] if output: lines += [ "set term png large", "set output '%s'" % output ] lines += [ "set xtics (%s)" % xlabels, "set grid", "set title 'Spectral signatures'", "set xrange [0:%d]" % xrange, "set noclabel", "set xlabel 'Bands'", "set ylabel 'DN Value'", "set style data lines" ] cmd = [] for i, row in enumerate(what): if not label: title = str(i + 1) else: title = str(tuple(row[0:2])) cmd.append("'data_%d' title '%s'" % (i, title)) cmd = ','.join(cmd) cmd = ' '.join(['plot', cmd, "with linespoints pt 779"]) lines.append(cmd) plotfile = os.path.join(tmp_dir, 'spectrum.gnuplot') plotf = open(plotfile, 'w') for line in lines: plotf.write(line + '\n') plotf.close() if output: grass.call(['gnuplot', plotfile]) else: grass.call(['gnuplot', '-persist', plotfile]) def draw_linegraph(what): yfiles = [] xfile = os.path.join(tmp_dir, 'data_x') xf = open(xfile, 'w') for j, val in enumerate(what[0][3:]): xf.write("%d\n" % (j + 1)) xf.close() for i, row in enumerate(what): yfile = os.path.join(tmp_dir, 'data_y_%d' % i) yf = open(yfile, 'w') for j, val in enumerate(row[3:]): yf.write("%s\n" % val) yf.close() yfiles.append(yfile) sienna = '#%02x%02x%02x' % (160, 82, 45) coral = '#%02x%02x%02x' % (255, 127, 80) gp_colors = ['red', 'green', 'blue', 'magenta', 'cyan', sienna, 'orange', coral] colors = gp_colors while len(what) > len(colors): colors += gp_colors colors = colors[0:len(what)] grass.run_command('d.linegraph', x_file=xfile, y_file=yfiles, y_color=colors, title='Spectral signatures', x_title='Bands', y_title='DN Value') def main(): group = options['group'] raster = options['raster'] output = options['output'] coords = options['coordinates'] label = flags['c'] gnuplot = flags['g'] global tmp_dir tmp_dir = grass.tempdir() if not group and not raster: grass.fatal(_("Either group= or raster= is required")) if group and raster: grass.fatal(_("group= and raster= are mutually exclusive")) #check if present if gnuplot and not grass.find_program('gnuplot'): grass.fatal(_("gnuplot required, please install first")) # get y-data for gnuplot-data file # get data from group files and set the x-axis labels if group: # ## PARSES THE GROUP FILES - gets rid of ugly header info from group #list output s = grass.read_command('i.group', flags='g', group=group, quiet=True) rastermaps = s.splitlines() else: # ## get data from list of files and set the x-axis labels rastermaps = raster.split(',') xlabels = ["'%s' %d" % (n, i + 1) for i, n in enumerate(rastermaps)] xlabels = ','.join(xlabels) what = [] s = grass.read_command('r.what', map=rastermaps, coordinates=coords, quiet=True) for l in s.splitlines(): f = l.split('|') for i, v in enumerate(f): if v in ['', '*']: f[i] = 0 else: f[i] = float(v) what.append(f) # build data files if gnuplot: draw_gnuplot(what, xlabels, output, label) else: draw_linegraph(what) if __name__ == "__main__": options, flags = grass.parser() atexit.register(cleanup) main()