#-------------------------------------------------------------------------------
# Name:        intreew-drains.py
# Purpose:     Reproduceren plaatje uit Jaarverslag ICW 1966 p16
# Created:     17-05-2012
# Copyright:   (c) Grondwaterformules.nl
# Licence:     Python 2.6
#-------------------------------------------------------------------------------
#!/usr/bin/env python
from scipy import *
import pylab

def plot(k,w,Pk,Pw,lb):

    # Define colors
    orange = '#FF9900'
    blue   = '#2F64B2'

    # Create plots
    pylab.subplots_adjust(bottom=0.22)
    graph = pylab.subplot(111, aspect='equal')
    pylab.hold(True)


    # Set ticklines
    xticklines = graph.get_xticklines()
    yticklines = graph.get_yticklines()

    for line in xticklines+yticklines:
        line.set_linewidth(3)

    for line in xticklines:
        line.set_visible(False)

    # Set gridlines
    gridlines = graph.get_xgridlines()
    gridlines.extend( graph.get_ygridlines() )
    for line in gridlines:
        line.set_linestyle('-')
    pylab.grid(True,which="both",ls="-")

    # Set ticklabels
    xticklabels = graph.get_xticklabels()
    yticklabels = graph.get_yticklabels()

    for label in xticklabels+yticklabels:
        label.set_color('k')
        label.set_fontsize('small')
        label.set_visible(True)

    # Create annotation
    pylab.xlabel('doorlatendheid k in m/dag',size=9)
    pylab.ylabel('intreeweerstand w in dagen/meter',size=9)

    # Plot
    graph.loglog(k,w,'-',color=blue, lw=2.0)
    graph.plot(Pk,Pw,markeredgewidth=0, linewidth=0,marker='o', markersize=6,color=orange)

    # annotation
    bbox_props = dict(boxstyle="round, pad=0.15", fc="white", ec="white", lw=1)
    for li in range(len(lb)):
        lbl = lb[li][0]
        xl = float(lb[li][1])
        yl = lb[li][2]
        graph.annotate(lbl,xy=(xl,yl),fontsize=10.0,bbox=bbox_props)

    # Filename
    filename = "intreeweerstand-drains.png"

    # Save figure to file
    fig = pylab.figure(1)
    mydpi = 100.0
    fig.set_dpi(mydpi)
    fig.set_size_inches(5.0,5.0)
    fig.savefig(filename,dpi=mydpi, facecolor='w', edgecolor='w')

    # Show the plot on screen
    pylab.show()


if __name__ == '__main__':

    # Array P : [locatie, doorlatendheid, intreeweerstand]
    # Coordinaten zijn afgeLezen uit figuur 6, jaarverslag ICW 1966, p16
    P = array([
    ['Groningen',0.0003,245.6095],
    ['Zeeland',0.0011,160.8614],
    [r'Steenbergen 0-70cm',0.0044,34.5662],
    ['Venhuizen',0.0735,1.1429],
    [r'Steenbergen 70-80cm',0.2411,0.5358],
    ['Blokzijlzand',0.1935,0.6696],
    ['Blokzijlzand',0.3209,0.4194],
    ['Blokzijlzand',0.5093,0.1758],
    ['Stuifzand',11.3283,0.0133],
    ['Stuifzand',11.8367,0.0106],
    ['Stuifzand',10.1483,0.0095],
    ['Stuifzand',13.5064,0.0078],
    ['Metselzand',80.2514,0.0050]
    ])

    # coordinaten zijn aangepast om labels op de mooiste plek te krijgen!
    Plabels = array([
    ['Groningen',0.0004,400.0],
    ['Zeeland',0.0015,150.0],
    ['Steenbergen\n0-70cm',0.007,34.5],
    ['Venhuizen',0.08,1.85],
    ['Steenbergen\n70-80cm',0.5,0.5],
    ['Blokzijlzand',0.8,0.2],
    ['Stuifzand',8.0,0.02],
    ['Metselzand',30.0,0.008]
    ])

    Pk = P[:,1]
    Pw = P[:,2]

    k = linspace(0.0001,100,100)
    w = 10**(-log10(k)-1)

    plot(k,w,Pk,Pw,Plabels)