Q-Q-plotDemo.py

#
# example script showing how to use opalsHisto to generate a quantile-quantile plot
#
# script parameters: <inputfile> [1|0]
#   <inputfile>: odm or rasterfile
#   1|0:         flag if plot window should be shown
#
import sys
import opals
from opals import Histo
import matplotlib.pylab as pl
import numpy
from scipy.stats import norm
 
def generate_plot(filename, showPlot = False, plotFile = 'QQplot.png'):
    # used probabilities for quantile-quantile plot
    probabilities = [0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99]
 
    # compute quantile values based on normal distribution
    dist = norm()
    quNorm = [ dist.ppf(p) for p in  probabilities]
    quNormArr = numpy.array(quNorm)
 
    # compute histogram including quantiles from
    his=Histo.Histo()
    his.commons.screenLogLevel = opals.Types.LogLevel.error
    his.inFile = [filename]
    his.probabilities = probabilities
    his.exactComputation = True  # to get exact quantile values
    his.plotFile = "off"       # do not create a svg file
    his.run()
    stats = his.histogram[0]
 
    # get quantiles of data
    quantiles = stats.getQuantiles()
    alpha, qu = zip(*quantiles)
    qu = list(qu)  # convert tuple to list
    quArr = numpy.array(qu)
 
    # compute trend line
    z = numpy.polyfit(quNormArr, quArr, 1)
    p = numpy.poly1d(z)
 
    # generate q-q plot
    pl.plot(quNormArr, quArr, 'bo-', ms = 5, mec = 'b', mfc = 'w')   # plot quantiles
    pl.plot(quNormArr, p(quNormArr), color="purple", linewidth=2, linestyle="--")   # plot trend line
    pl.title('Q-Q Plot')
    pl.xlabel('Quantiles of standard normal distribution')
    pl.ylabel('Quantiles of residuals')
    pl.grid(True)
 
    # build legend text (insert line feed after 4 numbers)
    text = 'Quantiles:'
    lf_after = 4  # line feed after 4 numbers
    for idx in range(0,len(probabilities),lf_after):
        text += '\n' + ', '.join(map(str,probabilities[idx:idx+lf_after]))
    pl.text(0.8, min(qu), text, backgroundcolor = 'w')
 
    # show and save plot
    if showPlot:
        pl.show()
    if plotFile:
        pl.savefig(plotFile, format='png')  # store plot as png file
 
 
if __name__ == "__main__":
    inputFile = None
    showPlot = False
 
    if len(sys.argv) == 1:
        print("ODM/Raster file misssing")
        sys.exit(-1)
    else:
        inputFile = sys.argv[1]
 
    if len(sys.argv) > 2:
        showPlot = bool(int(sys.argv[2]))
 
    generate_plot(inputFile, showPlot)