model_scripts/ampang.py at master · EnvModellingGroup/model_scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
#!/usr/bin/env python
# Copyright: This document has been placed in the public domain.
"""
Yannick Copin's pythoncode for Taylor diagram (Taylor, 2001) changed to tidal amplitude & phase diagram by Takai Tsubono (2021/11/11)
Note: If you have found these software useful for your research, I would
appreciate an acknowledgment.

Tsubono, T., 2022. Diagram statistically displaying model performance for tides or quasi-periodic oscillations, Deep Sea Research I, Vol. 180, 103686,
https://www.sciencedirect.com/science/article/pii/S0967063721002235?via%3Dihub

"""
__version__ = "modified from Time-stamp: <2018-12-06 11:43:41 ycopin>"
__author__ = "modified from Yannick Copin's code "
import numpy as np
import matplotlib.pyplot as plt
class AmpPhsDiagram(object):
#" AmpPhs diagram.
#    Plot model constants of constituents and correlation to reference (data)
#    sample in a single-quadrant polar plot, r=amplitude and
#    theta= phaselag          .  "
    def __init__(self, reff=1.0, fig=None, rect=111, label='_', amprange=[0, 1.5],ampthck=[0,1.5,0.25],phsrange=[-30,30],phsthck=[-30,30,15],amptitle='amplitude',phstitle='phase lag',rd_fmt="%.1f"):
        """Set up amplitude and phase reffered to Taylor diagram axes, i.e. polar
        plot, using `mpl_toolkits.axisartist.floating_axes`.
        Parameters:
        * reff: reference standard deviation to be compared to
        * fig: input Figure or None
        * rect: subplot definition
        * label: reference label
        * amprange & phsrange: range angurar and radius
        * ampthck & phsthick : thick def. for angular and radius
        """
        from matplotlib.projections import PolarAxes
        import mpl_toolkits.axisartist.floating_axes as FA
        import mpl_toolkits.axisartist.grid_finder as GF
        self.reff = reff            # Reference standard deviation
        self.amprange = amprange
        self.phsrange = phsrange
        idv = np.int_( ( ampthck[1]+0.0000001 - ampthck[0] ) / ampthck[2] )+1
        jdv = np.int_( ( phsthck[1]+0.0000001 - phsthck[0] ) / phsthck[2] )+1

        tr = PolarAxes.PolarTransform()
        degree_ticks = lambda d: (d*np.pi/180, "%d$^\\circ$"%(d))
#       angle_ticks = map(degree_ticks, np.linspace(phsthck[0],phsthck[1],jdv))
#       grid_locator1 = GF.FixedLocator([v for v, s in angle_ticks])
        angle_ticks = dict(map(degree_ticks, np.linspace(phsthck[0],phsthck[1],jdv)))
        grid_locator1 = GF.FixedLocator( angle_ticks.keys() )
        tick_formatter1 = GF.DictFormatter(angle_ticks)
#       angle_ticks = map(degree_ticks, np.linspace(phsthck[0],phsthck[1],jdv))
#       tick_formatter1 = GF.DictFormatter(dict(angle_ticks))

#       STDgrid = np.round(np.linspace(ampthck[0],ampthck[1]+0.0000,idv),1)
#       STDgrid = np.arange(ampthck[0],ampthck[1]+0.0000,ampthck[2])
        STDgrid = np.round(np.linspace(ampthck[0],ampthck[1]+0.0001,idv),3)
        grid_locator2 = GF.FixedLocator(STDgrid)
#       radial_func  = lambda d: (d, "%.1f"%(d))
        radial_func  = lambda d: (d, rd_fmt%(d))
        radial_ticks = dict(map(radial_func, STDgrid))
        tick_formatter2 = GF.DictFormatter(radial_ticks)
#       tick_formatter2 = GF.DictFormatter(dict(zip(STDgrid,map(str,STDgrid))))

        phspi=phsrange/180*np.pi
        gh = FA.GridHelperCurveLinear(tr,
             extremes=(phspi[0], phspi[1], amprange[0],amprange[1]),
             grid_locator1=grid_locator1,
             grid_locator2=grid_locator2,
             tick_formatter1=tick_formatter1,
             tick_formatter2=tick_formatter2)
#                               extremes=(-np.pi/4.*1, np.pi/4.*1, smin, smax),

#       fig = plt.figure(figsize=(8,6))
        ax = FA.FloatingSubplot(fig, rect, grid_helper=gh)
        fig.add_subplot(ax)

        ax.axis["top"].set_axis_direction("bottom")  # "Angle axis"
        ax.axis["top"].toggle(ticklabels=True, label=True)
        ax.axis["top"].major_ticklabels.set_axis_direction("top")
        ax.axis["top"].label.set_axis_direction("top")
        ax.axis["top"].label.set_text(phstitle)
        ax.axis["bottom"].set_axis_direction("bottom")  # "Angle axis"
        ax.axis["bottom"].toggle(ticklabels=False, label=False)
        ax.axis["right"].set_axis_direction("left") # "X axis"
        ax.axis["right"].toggle(ticklabels=True, label=True)
        ax.axis["right"].label.set_text(amptitle)
        ax.axis["left"].set_axis_direction("bottom") # "X axis"
        ax.axis["left"].toggle(ticklabels=False, label=False)
#       ax.grid(color='0.7')
        ax.grid(color='dimgray', linewidth =1.5 )
    def add_contours(self, levs=[0.1,0.2,0.3],  **kwargs):
        piphi=np.pi/180.*self.phsrange
        rs,ts=np.meshgrid( np.linspace(self.amprange[0],self.amprange[1],num=100), np.linspace(piphi[0],piphi[1],num=100) )
        rms = np.sqrt( ( self.reff)**2 + rs**2 -2*(self.reff)*rs*np.cos(ts))
        contours = plt.contour(rs*np.cos(ts),rs*np.sin(ts),rms,levels=levs,**kwargs)
        return contours
    def bdd_contours(self, levs=[0.1,0.2,0.3],  **kwargs):
        piphi=np.pi/180.*self.phsrange
        rs,ts=np.meshgrid( np.linspace(self.amprange[0],self.amprange[1],num=100), np.linspace(piphi[0],piphi[1],num=100) )
        rms = np.sqrt( ( self.reff)**2 + rs**2 -2*(self.reff)*rs*np.cos(ts))
#       rms = np.sqrt( (rs*np.cos(ts)-self.reff.real)**2+ (rs*np.sin(ts)-self.reff.imag) **2     )
        contours = plt.contour(rs*np.cos(ts),rs*np.sin(ts),rms,levels=levs,**kwargs)
        return contours

if __name__=='__main__':
#
    cal=np.zeros((4))+[10.80-5.03j, 8.44-5.40j, 6.84-4.10j, 8.71-4.83j]
    obs=np.zeros((4))+[14.78-0.45j,12.89-0.27j,11.21-0.14j,12.93-0.00j]
    rat=cal/obs

    smin=5.0
    smax=15
    swdt=2.5
    tmin=-40.
    tmax=20.
    twdt=10.
    aaa=' M2_amplitude'
    bbb='phase lag'
    refft=np.abs(obs[3])

    arng=np.array((smin-0.5,smax+0.5))
    prng=np.array((tmin-2.0,tmax+2.0))
    athc=np.array((smin,smax,swdt))
    pthc=np.array((tmin,tmax,twdt))

    fig = plt.figure(figsize=(8,6))
    dia = AmpPhsDiagram(reff=refft,fig=fig,amprange=arng,ampthck=athc,phsrange=prng,phsthck=pthc,amptitle=aaa,phstitle=bbb)
    contours=dia.add_contours(levs=[2.5,5.0,7.5],colors='0.5')
    plt.clabel(contours, inline=1, fontsize=10,fmt="%.1f")
    #plt.scatter(1,0,marker='*',s=70)
    plt.scatter(obs[0:4].real,obs[0:4].imag,marker='s', s=40, color='orange',label='Obs.')
    plt.scatter(obs[3].real,obs[3].imag,marker='s', s=42, color='orange',edgecolor='black')
    plt.scatter(cal[0:4].real,cal[0:4].imag,marker='o', s=40, color='white',edgecolor='dodgerblue',label='Cal.')
    plt.legend(loc="lower left",ncol=3,scatterpoints=1)
    plt.savefig('test_orig.png')
    plt.show()

    smin=0.5
    smax=1.5
    swdt=0.25
    tmin=-30.
    tmax=30
    twdt=15
    aaa='ratio of M2_amplitude(cal/obs)'
    bbb='phase lag'
    refft=1.0

    arng=np.array((smin-0.10,smax+0.1))
    prng=np.array((tmin-5.00,tmax+5.0))
    athc=np.array((smin,smax,swdt))
    pthc=np.array((tmin,tmax,twdt))

    fig = plt.figure(figsize=(8,6))
    dia = AmpPhsDiagram(reff=refft,fig=fig,amprange=arng,ampthck=athc,phsrange=prng,phsthck=pthc,amptitle=aaa,phstitle=bbb)
    contours=dia.add_contours(levs=[0.1,0.2,0.3],colors='0.5')
    plt.clabel(contours, inline=1, fontsize=10,fmt="%.1f")

    plt.savefig('test_ratio.png')
    plt.show()