-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathplotta_Hayashi.py
122 lines (100 loc) · 4.21 KB
/
plotta_Hayashi.py
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
#!/usr/bin python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import os
import sys
from netCDF4 import Dataset
import subprocess
# input
threshold=float(sys.argv[1]) #4.5
#threshold=3.5 #only select years with heatwaves above this threshold
#threshold=-np.inf # select all heatwaves
area = 'Scandinavia'
T = 30 # running mean in days
latnorth=int(sys.argv[3])#80#60#70#75#80
latsouth=int(sys.argv[2])#70#50#60#55#40
#output
outfolder='/scratch/gmiloshe/ERA/Hayashi'
infolder='/scratch/gmiloshe/ERA/Hayashi'
inname='test_Hayashi_output'+area+'_'+str(T)+'_'+str(threshold)+'_'+str(latsouth)+'_'+str(latnorth)
# output
outfolder=infolder
outname='Hayashi_test'+area+'_'+str(T)+'_'+str(threshold)+'_'+str(latsouth)+'_'+str(latnorth)
###############################################################
npzfile=np.load(infolder+'/'+inname+'.npz')
locals().update(npzfile)
#normalization terms (frequency and wavenumber)
time=np.arange(0,0.5,1/ntime)
zonal=np.arange(0,nlon/2+1)
#time period
period=1/time[1:int(ntime/2)+1]
#normalization
for itime in range(int(ntime/2)):
for ilon in range(int(nlon/2)):
avepeast[itime,ilon] = avepeast[itime,ilon] *time[itime]*zonal[ilon]*ntime
avepwest[itime,ilon] = avepwest[itime,ilon] *time[itime]*zonal[ilon]*ntime
avetr[itime,ilon] = avetr[itime,ilon] *time[itime]*zonal[ilon]*ntime
avetreast[itime,ilon]=avetreast[itime,ilon] *time[itime]*zonal[ilon]*ntime
avetrwest[itime,ilon]=avetrwest[itime,ilon] *time[itime]*zonal[ilon]*ntime
avetotF[itime,ilon] = avetotF[itime,ilon] *time[itime]*zonal[ilon]*ntime
avestaF[itime,ilon] = avestaF[itime,ilon] *time[itime]*zonal[ilon]*ntime
avestaP[itime,ilon] = avestaP[itime,ilon] *time[itime]*zonal[ilon]*ntime
avetotP[itime,ilon] = avetotP[itime,ilon] *time[itime]*zonal[ilon]*ntime
#meshgrid (using period instead of frequency)
wn,per=np.meshgrid(zonal[1:12],period[0:int(ntime/2)+1])
avetotF=np.real(avetotF)
avestaF=np.real(avestaF)
avepeast=np.real(avepeast)
avepwest=np.real(avepwest)
print(np.max(avetotF[1:int(ntime/2),1:12]))
color_levels = np.arange(0,80000+8000,8000)
col_map = "viridis"
#total spectrum
plt.figure(figsize=(6, 6))
plt.contourf(per,wn,avetotF[1:int(ntime/2),1:12], color_levels, cmap=col_map, extend='both')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('Period (days)')
plt.ylabel('Wavenumber')
plt.xlim((100,2))
plt.colorbar()
plt.contour(per,wn,avetotF[1:int(ntime/2),1:12], color_levels, colors='black', linewidths=.5, extend='both')
plt.title(r'$H_T(k,\omega)$ th = '+str(threshold)+' over ' +str(latsouth)+' to '+str(latnorth))
plt.savefig(outfolder+'/'+outname+'_avetotF.png')
#stationary component
plt.figure(figsize=(6, 6))
plt.contourf(per,wn,avestaF[1:int(ntime/2),1:12], color_levels, cmap=col_map, extend='both')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('Period (days)')
plt.ylabel('Wavenumber')
plt.xlim((100,2))
plt.colorbar()
plt.contour(per,wn,avestaF[1:int(ntime/2),1:12], color_levels, colors='black', linewidths=.5, extend='both')
plt.title(r'$H_S(k,\omega)$ th = '+str(threshold)+' over ' +str(latsouth)+' to '+str(latnorth))
plt.savefig(outfolder+'/'+outname+'_avestaF.png')
#eastward propagating component
plt.figure(figsize=(6, 6))
plt.contourf(per,wn,avepeast[1:int(ntime/2),1:12], color_levels, cmap=col_map, extend='both')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('Period (days)')
plt.ylabel('Wavenumber')
plt.xlim((100,2))
plt.colorbar()
plt.contour(per,wn,avepeast[1:int(ntime/2),1:12], color_levels, colors='black', linewidths=.5, extend='both')
plt.title(r'$H_E(k,\omega)$ th = '+str(threshold)+' over ' +str(latsouth)+' to '+str(latnorth))
plt.savefig(outfolder+'/'+outname+'_avepeast.png')
#westward propagating component
plt.figure(figsize=(6, 6))
plt.contourf(per,wn,avepwest[1:int(ntime/2),1:12], color_levels, cmap=col_map, extend='both')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('Period (days)')
plt.ylabel('Wavenumber')
plt.xlim((100,2))
plt.colorbar()
plt.contour(per,wn,avepwest[1:int(ntime/2),1:12], color_levels, colors='black', linewidths=.5, extend='both')
plt.title(r'$H_W(k,\omega)$ th = '+str(threshold)+' over ' +str(latsouth)+' to '+str(latnorth))
plt.savefig(outfolder+'/'+outname+'_avepwest.png')