-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcomposite_cldmaps.py
109 lines (79 loc) · 4.05 KB
/
composite_cldmaps.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
import numpy as np
import numpy.ma as ma
import iris as iris
import iris.plot as iplt
import iris.quickplot as qplt
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import mycmaps as mc
import scipy.stats as stats
import sys
import troposave as ta
# Load in RHsfc, RHtropos (or RH300hPa?)
# calculate the regression between the timeseries
# RHsfc = alpha * RHtropo + epsilon
ncfile_path = '/home/nicholat/project/pacemaker/ncfiles/'
# import the ACCESS data using iris
cld = iris.load_cube(ncfile_path + 'cld.thlev.4ysl.m48.nc')
cloud_clim = iris.load_cube(ncfile_path + 'cld.thlev.4ysl.nc')
cld.coord('Hybrid height').standard_name = 'atmosphere_hybrid_height_coordinate'
cloud_clim.coord('Hybrid height').standard_name = 'atmosphere_hybrid_height_coordinate'
# Define regions
cld.coord('latitude').guess_bounds()
cld.coord('longitude').guess_bounds()
cloud_clim.coord('latitude').guess_bounds()
cloud_clim.coord('longitude').guess_bounds()
# define variance for cloud
cloud_clim_std = cloud_clim.collapsed('time',iris.analysis.STD_DEV)
mons = 6
lag = 0
max_i = 35 + lag
max_f = max_i + mons
min_i = 11 + lag
min_f = min_i + mons
high = iris.Constraint(atmosphere_hybrid_height_coordinate = lambda h: 5000 <= h <= 15000)
low = iris.Constraint(atmosphere_hybrid_height_coordinate = lambda h: 0 <= h <= 5000)
full = iris.Constraint(atmosphere_hybrid_height_coordinate = lambda h: 0 <= h <= 10000)
cld_max = cld[max_i:max_f,::].collapsed('time',iris.analysis.MEAN)
cld_min = cld[min_i:min_f,::].collapsed('time',iris.analysis.MEAN)
cld_high_max = cld_max.extract(high).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cld_high_min = cld_min.extract(high).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cld_low_max = cld_max.extract(low).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cld_low_min = cld_min.extract(low).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cld_full_max = cld_max.extract(full).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cld_full_min = cld_min.extract(full).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cloud_std = cloud_clim_std.extract(full).collapsed('atmosphere_hybrid_height_coordinate',iris.analysis.MEAN)
cld_by_std = cld_full_max/(cloud_std)
rhmm = 0.02
plt.figure(1)
qplt.pcmeshclf(cld_by_std,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite())
plt.title('Normalized Cloud response to max positive forcing, '+str(lag)+'-'+str(mons+lag)+' months')
# plt.savefig('figures/comp_cldfull_max.png')
sys.exit()
plt.figure(3)
qplt.pcmeshclf(cld_high_max,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite())
plt.title('High Cloud response to max positive forcing. '+str(lag)+'-'+str(mons+lag)+' months')
plt.savefig('figures/comp_cldhigh_max.png')
# plt.figure(4)
# qplt.pcmeshclf(cld_high_min,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite_r())
# plt.title('High Cloud response to min negative forcing. '+str(lag)+'-'+str(mons+lag)+' months')
# plt.savefig('figures/comp_cldhigh_min.png')
plt.figure(5)
qplt.pcmeshclf(cld_low_max,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite())
plt.title('Low Cloud response to max positive forcing, 700hPa. '+str(lag)+'-'+str(mons+lag)+' months')
plt.savefig('figures/comp_cldlow_max.png')
# plt.figure(6)
# qplt.pcmeshclf(cld_low_min,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite_r())
# plt.title('Low Cloud response to min negative forcing, 700hPa. '+str(lag)+'-'+str(mons+lag)+' months')
# plt.savefig('figures/comp_cldlow_min.png')
plt.figure(7)
qplt.pcmeshclf(cld_full_max,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite())
plt.title('All Cloud response to max positive forcing, 1000hPa. '+str(lag)+'-'+str(mons+lag)+' months')
plt.savefig('figures/comp_cldfull_max.png')
# plt.figure(8)
# qplt.pcmeshclf(cld_full_min,vmin=-rhmm,vmax=rhmm,cmap=mc.jetwhite_r())
# plt.title('All Cloud response to min negative forcing, 1000hPa'+str(mons)+'-'+str(mons+lag)+' months')
# plt.savefig('figures/comp_cldfull_min.png')
# plt.figure(3)
# qplt.pcmeshclf(RH300max - RH300min,vmin=-1,vmax=1,cmap=mc.jetwhite())
# plt.savefig('figures/reg_RH_sfc_RH_sfc.png')