Geografic-Information-System/scripts/d.polar/d.polar.py

#!/usr/bin/env python
############################################################################
#
# MODULE:       d.polar
# AUTHOR(S):    Markus Neteler. neteler itc.it
#               algorithm + EPS output by Bruno Caprile
#               d.graph plotting code by Hamish Bowman
#               Converted to Python by Glynn Clements
# PURPOSE:      Draws polar diagram of angle map. The outer circle considers
#               all cells in the map. If one or many of them are NULL (no data),
#               the figure will not reach the outer circle. The vector inside
#               indicates the prevalent direction.
# COPYRIGHT:    (C) 2006,2008 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.
#
#############################################################################

#%Module
#%  description: Draws polar diagram of angle map such as aspect or flow directions
#%  keywords: display, diagram
#%End
#%option
#% key: map
#% type: string
#% key_desc: name
#% gisprompt: old,cell,raster
#% description: Name of raster angle map
#% required : yes
#%End
#%option
#% key: undef
#% type: double
#% description: Pixel value to be interpreted as undefined (different from NULL)
#% required : no
#%End
#%option
#% key: eps
#% type: string
#% gisprompt: new_file,file,output
#% description: Name of optional EPS output file
#% required : no
#%end
#%flag
#% key: x
#% description: Plot using Xgraph
#%end

import sys
import os
import string
import types
import math
import atexit
import glob
import shutil
import grass

def cleanup():
    grass.try_remove(tmp)
    for f in glob.glob(tmp + '_*'):
	grass.try_remove(f)

def plot_xgraph():
    newline = ['\n']
    p = grass.Popen(['xgraph'], stdin = grass.PIPE)
    for point in sine_cosine_replic + newline + outercircle + newline + vector:
	if isinstance(point, types.TupleType):
	    p.stdin.write("%f %f\n" % point)
	else:
	    p.stdin.write(point + '\n')
    p.stdin.close()
    p.wait()

def plot_dgraph():
    # use d.info and d.frame to create a square frame in the center of the window.
    s = grass.read_command('d.info', flags = 'd')
    f = s.split()
    frame_width = float(f[1])
    frame_height = float(f[2])

    # take shorter side as length of frame side
    min_side = min(frame_width, frame_height)
    dx = (frame_width  - min_side) / 2
    dy = (frame_height - min_side) / 2

    fl = dx
    fr = frame_width - dx
    ft = dy
    fb = frame_height - dy

    tenv = os.environ.copy()
    tenv['GRASS_FRAME'] = '%f,%f,%f,%f' % (ft, fb, fl, fr)

    # polyline calculations
    ring = 0.95
    scaleval = ring * totalvalidnumber / totalnumber

    sine_cosine_replic_normalized = [
	((scaleval * p[0]/maxradius + 1)*50,
	 (scaleval * p[1]/maxradius + 1)*50)
	for p in sine_cosine_replic
	if isinstance(p, types.TupleType)]

    # create circle
    circle = [(50*(1 + ring * math.sin(math.radians(i))),
	       50*(1 + ring * math.cos(math.radians(i))))
	      for i in range(0, 361)]

    # trend vector
    vect = [((scaleval * p[0]/maxradius + 1)*50,
	     (scaleval * p[1]/maxradius + 1)*50)
	    for p in vector[1:]]

    # Possible TODOs:
    # To fill data area with color, use BOTH d.graph's polyline and polygon commands.
    #  Using polygon alone gives a jagged boundary due to sampling technique (not a bug).

    # plot it!
    lines = [
	# draw circle
	#   mandatory when drawing proportional to non-square frame
	"color 180:255:180",
	"polyline"] + circle + [
	# draw axes
        "color 180:180:180",
        "width 0",
        "move 0 50",
        "draw 100 50",
        "move 50 0",
        "draw 50 100",

	# draw the goods
        "color red",
        "width 0",
        "polyline"] + sine_cosine_replic_normalized + [
	# draw vector
        "color blue",
        "width 3",
        "polyline"] + vect + [

	# draw compass text
        "color black",
        "width 2",
        "size 3 3",
        "move 51 97",
        "text N",
        "move 51 1",
        "text S",
        "move 1 51",
        "text W",
        "move 97 51",
        "text E",

	# draw legend text
        "width 0",
        "size 2",
        "color 0:180:0",
        "move 1.5 96.5",
        "text All data (incl. NULLs)",
        "color red",
        "move 1.5 93.5",
        "text Real data angles",
        "color blue",
        "move 1.5 90.5",
        "text Avg. direction"
	]

    p = grass.feed_command('d.graph', env = tenv)
    for point in lines:
	if isinstance(point, types.TupleType):
	    p.stdin.write("%f %f\n" % point)
	else:
	    p.stdin.write(point + '\n')
    p.stdin.close()
    p.wait()

def plot_eps(psout):
    # EPS output (by M.Neteler and Bruno Caprile, ITC-irst)
    grass.message("Generating %s ..." % psout)

    outerradius = maxradius
    epsscale = 0.1
    frameallowance = 1.1
    halfframe = 3000
    center = (halfframe, halfframe)
    scale = halfframe / (outerradius * frameallowance)

    diagramlinewidth = halfframe / 400
    axeslinewidth = halfframe / 500
    axesfontsize = halfframe / 16
    diagramfontsize = halfframe / 20
    halfframe_2 = halfframe * 2

    averagedirectioncolor = 1 #(blue)
    diagramcolor = 4 #(red)
    circlecolor = 2 #(green)
    axescolor = 0 #(black)

    northjustification = 2
    eastjustification = 6
    southjustification = 8
    westjustification = 8

    northxshift = 1.02 * halfframe 
    northyshift = 1.98 * halfframe 
    eastxshift  = 1.98 * halfframe  
    eastyshift  = 1.02 * halfframe  
    southxshift = 1.02 * halfframe 
    southyshift = 0.02 * halfframe 
    westxshift  = 0.01 * halfframe  
    westyshift  = 1.02 * halfframe  

    alldatastring = "all data (null included)"
    realdatastring = "real data angles"
    averagedirectionstring = "avg. direction"

    legendsx = 1.95 * halfframe
    alldatalegendy = 1.95 * halfframe
    realdatalegendy = 1.90 * halfframe
    averagedirectionlegendy = 1.85 * halfframe

    ##########
    outf = file(psout, 'w')

    prolog = os.path.join(os.environ['GISBASE'], 'etc', 'd.polar', 'ps_defs.eps')
    inf = file(prolog)
    shutil.copyfileobj(inf, outf)
    inf.close()

    t = string.Template("""
$EPSSCALE $EPSSCALE scale                           %% EPS-SCALE EPS-SCALE scale
%%
%% drawing axes
%%

col0                                    %% colAXES-COLOR
$AXESLINEWIDTH setlinewidth                          %% AXES-LINEWIDTH setlinewidth
[] 0 setdash
newpath
 $HALFFRAME     0.0 moveto                  %% HALF-FRAME 0.0 moveto
 $HALFFRAME  $HALFFRAME_2 lineto                  %% HALF-FRAME (2 * HALF-FRAME) lineto
    0.0  $HALFFRAME moveto                  %% 0.0 HALF-FRAME moveto
 $HALFFRAME_2  $HALFFRAME lineto                  %% (2 * HALF-FRAME) HALF-FRAME lineto
stroke

%%
%% drawing outer circle
%%

col2                                    %% colCIRCLE-COLOR
$DIAGRAMFONTSIZE /Times-Roman choose-font            %% DIAGRAM-FONTSIZE /Times-Roman choose-font
$DIAGRAMLINEWIDTH setlinewidth                          %% DIAGRAM-LINEWIDTH setlinewidth
[] 0 setdash
newpath
                                        %% coordinates of rescaled, translated outer circle follow
                                        %% first point moveto, then lineto
""")
    s = t.substitute(
	AXESLINEWIDTH = axeslinewidth,
	DIAGRAMFONTSIZE = diagramfontsize,
	DIAGRAMLINEWIDTH = diagramlinewidth,
	EPSSCALE = epsscale,
	HALFFRAME = halfframe,
	HALFFRAME_2 = halfframe_2
	)
    outf.write(s)

    sublength = len(outercircle) - 2
    (x, y) = outercircle[1]
    outf.write("%.2f %.2f moveto\n" % (x * scale + halfframe, y * scale + halfframe))
    for x, y in outercircle[2:]:
        outf.write("%.2f %.2f lineto\n" % (x * scale + halfframe, y * scale + halfframe))

    t = string.Template("""
stroke

%%
%% axis titles
%%

col0                                    %% colAXES-COLOR
$AXESFONTSIZE /Times-Roman choose-font            %% AXES-FONTSIZE /Times-Roman choose-font
(N) $NORTHXSHIFT $NORTHYSHIFT $NORTHJUSTIFICATION just-string         %% NORTH-X-SHIFT NORTH-Y-SHIFT NORTH-JUSTIFICATION just-string
(E) $EASTXSHIFT $EASTYSHIFT $EASTJUSTIFICATION just-string         %% EAST-X-SHIFT EAST-Y-SHIFT EAST-JUSTIFICATION just-string
(S) $SOUTHXSHIFT $SOUTHYSHIFT $SOUTHJUSTIFICATION just-string           %% SOUTH-X-SHIFT SOUTH-Y-SHIFT SOUTH-JUSTIFICATION just-string
(W) $WESTXSHIFT $WESTYSHIFT $WESTJUSTIFICATION just-string           %% WEST-X-SHIFT WEST-Y-SHIFT WEST-JUSTIFICATION just-string
$DIAGRAMFONTSIZE /Times-Roman choose-font            %% DIAGRAM-FONTSIZE /Times-Roman choose-font


%%
%% drawing real data diagram
%%

col4                                    %% colDIAGRAM-COLOR
$DIAGRAMLINEWIDTH setlinewidth                          %% DIAGRAM-LINEWIDTH setlinewidth
[] 0 setdash
newpath
                                        %% coordinates of rescaled, translated diagram follow
                                        %% first point moveto, then lineto
""")
    s = t.substitute(
	AXESFONTSIZE = axesfontsize,
	DIAGRAMFONTSIZE = diagramfontsize,
	DIAGRAMLINEWIDTH = diagramlinewidth,
	EASTJUSTIFICATION = eastjustification,
	EASTXSHIFT = eastxshift,
	EASTYSHIFT = eastyshift,
	NORTHJUSTIFICATION = northjustification,
	NORTHXSHIFT = northxshift,
	NORTHYSHIFT = northyshift,
	SOUTHJUSTIFICATION = southjustification,
	SOUTHXSHIFT = southxshift,
	SOUTHYSHIFT = southyshift,
	WESTJUSTIFICATION = westjustification,
	WESTXSHIFT = westxshift,
	WESTYSHIFT = westyshift
	)
    outf.write(s)

    sublength = len(sine_cosine_replic) - 2
    (x, y) = sine_cosine_replic[1]
    outf.write("%.2f %.2f moveto\n" % (x * scale + halfframe, y * scale + halfframe))
    for x, y in sine_cosine_replic[2:]:
        outf.write("%.2f %.2f lineto\n" % (x * scale + halfframe, y * scale + halfframe))

    t = string.Template("""
stroke
%%
%% drawing average direction
%%

col1                                    %% colAVERAGE-DIRECTION-COLOR
$DIAGRAMLINEWIDTH setlinewidth                          %% DIAGRAM-LINEWIDTH setlinewidth
[] 0 setdash
newpath
                                        %% coordinates of rescaled, translated average direction follow
                                        %% first point moveto, second lineto
""")
    s = t.substitute(DIAGRAMLINEWIDTH = diagramlinewidth)
    outf.write(s)

    sublength = len(vector) - 2
    (x, y) = vector[1]
    outf.write("%.2f %.2f moveto\n" % (x * scale + halfframe, y * scale + halfframe))
    for x, y in vector[2:]:
        outf.write("%.2f %.2f lineto\n" % (x * scale + halfframe, y * scale + halfframe))

    t = string.Template("""
stroke

%%
%% drawing legends
%%

col2                                    %% colCIRCLE-COLOR
%% Line below: (ALL-DATA-STRING) LEGENDS-X ALL-DATA-LEGEND-Y 4 just-string
($ALLDATASTRING) $LEGENDSX $ALLDATALEGENDY 4 just-string

col4                                    %% colDIAGRAM-COLOR
%% Line below: (REAL-DATA-STRING) LEGENDS-X REAL-DATA-LEGEND-Y 4 just-string
($REALDATASTRING) $LEGENDSX $REALDATALEGENDY 4 just-string

col1                                    %% colAVERAGE-DIRECTION-COLOR
%% Line below: (AVERAGE-DIRECTION-STRING) LEGENDS-X AVERAGE-DIRECTION-LEGEND-Y 4 just-string
($AVERAGEDIRECTIONSTRING) $LEGENDSX $AVERAGEDIRECTIONLEGENDY 4 just-string
""")
    s = t.substitute(
	ALLDATALEGENDY = alldatalegendy,
	ALLDATASTRING = alldatastring,
	AVERAGEDIRECTIONLEGENDY = averagedirectionlegendy,
	AVERAGEDIRECTIONSTRING = averagedirectionstring,
	LEGENDSX = legendsx,
	REALDATALEGENDY = realdatalegendy,
	REALDATASTRING = realdatastring
	)
    outf.write(s)

    outf.close()

    grass.message("Done.")

def main():
    global tmp
    global sine_cosine_replic, outercircle, vector
    global totalvalidnumber, totalnumber, maxradius

    map = options['map']
    undef = options['undef']
    eps = options['eps']
    xgraph = flags['x']

    tmp = grass.tempfile()

    if eps and xgraph:
	grass.fatal("Please select only one output method")

    #### check if we have xgraph (if no EPS output requested)
    if xgraph and not grass.find_program('xgraph'):
	grass.fatal("xgraph required, please install first (www.xgraph.org)")

    #################################
    # this file contains everthing:
    rawfile = tmp + "_raw"
    rawf = file(rawfile, 'w')
    grass.run_command('r.stats', flags = '1', input = map, stdout = rawf)
    rawf.close()

    rawf = file(rawfile)
    totalnumber = 0
    for line in rawf:
	totalnumber += 1
    rawf.close()

    grass.message("Calculating statistics for polar diagram... (be patient)")

    #wipe out NULL data and undef data if defined by user
    # - generate degree binned to integer, eliminate NO DATA (NULL):
    # change 360 to 0 to close polar diagram:
    rawf = file(rawfile)
    nvals = 0
    sumcos = 0
    sumsin = 0
    freq = {}
    for line in rawf:
	line = line.rstrip('\r\n')
	if line in ['*', undef]:
	    continue
	nvals += 1
	x = float(line)
	rx = math.radians(x)
	sumcos += math.cos(rx)
	sumsin += math.sin(rx)
	ix = round(x)
	if ix == 360:
	    ix = 0
	if ix in freq:
	    freq[ix] += 1
	else:
	    freq[ix] = 1
    rawf.close()

    totalvalidnumber = nvals
    if totalvalidnumber == 0:
	grass.fatal("No data pixel found")

    #################################
    # unit vector on raw data converted to radians without no data:

    unitvector = (sumcos / nvals, sumsin / nvals)

    #################################
    # how many are there?:

    occurrences = [(math.radians(x), freq[x]) for x in freq]
    occurrences.sort()
    
    # find the maximum value
    maxradius = max([f for a, f in occurrences])

    # now do cos() sin()
    sine_cosine = [(math.cos(a) * f, math.sin(a) * f) for a, f in occurrences]

    sine_cosine_replic = ['"Real data angles'] + sine_cosine + sine_cosine[0:1]

    if eps or xgraph:
	outercircle = []
	outercircle.append('"All Data incl. NULLs')
	scale = 1.0 * totalnumber / totalvalidnumber * maxradius
	for i in range(0, 361):
	    a = math.radians(i)
	    x = math.cos(a) * scale
	    y = math.sin(a) * scale
	    outercircle.append((x, y))

    # fix vector length to become visible (x? of $MAXRADIUS):
    vector = []
    vector.append('"Avg. Direction\n')
    vector.append((0, 0))
    vector.append((unitvector[0] * maxradius, unitvector[1] * maxradius))

    ###########################################################

    # Now output:

    if eps:
	psout = grass.basename(eps, 'eps') + '.eps'
	plot_eps(psout)
    elif xgraph:
	plot_xgraph()
    else:
	plot_dgraph()

    grass.message("Average vector:")
    grass.message("direction: %.1f degrees CCW from East" % math.degrees(math.atan2(unitvector[1], unitvector[0])))
    grass.message("magnitude: %.1f percent of fullscale" % (100 * math.hypot(unitvector[0], unitvector[1])))

if __name__ == "__main__":
    options, flags = grass.parser()
    atexit.register(cleanup)
    main()