-
Notifications
You must be signed in to change notification settings - Fork 0
/
fig2.era5.precip%from1991-2020.py
131 lines (93 loc) · 3.57 KB
/
fig2.era5.precip%from1991-2020.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
123
124
125
126
127
128
129
130
131
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import pandas as pd
'''
Script summary: WY 2021 sesaonal departure of precipitation from the 1991-2020 climatological means. Autumn 2021 (OND), winter 2022 (JFM), spring 2022 (AMJ), and summer 2022 (JAS).
Author: Siiri Bigalke
Last edited: October 13, 2022
'''
# ==========================
''' Load ERA5 Data '''
# ==========================
pwd = '/work2/Reanalysis/ERA5/ERA5_monthly/monolevel/'+ \
'ERA5.mon.tp.1950-2022.nc'
ds = xr.open_dataset(pwd).sel(
latitude=slice(90, 60),
longitude=slice(0, 360))['tp']
ds = ds.mean('expver') # As of 10.13.2022 september data is still ERA5T
lon = ds.coords['longitude']
lat = ds.coords['latitude']
X, Y = np.meshgrid(lon, lat)
# Weight by cosine of latitude
weights = np.cos(np.deg2rad(ds.latitude))
weights.name = 'weights'
wds = ds.weighted(weights)
ds = wds.mean(('latitude', 'longitude'))
# Create new coordinate for water year
water_year = (ds.time.dt.month >= 10) + ds.time.dt.year
ds.coords['water_year'] = water_year
# =========================
# Trend Analysis
# =========================
def seasonal_trend(ds, months, name):
s = ds[ds.time.dt.month.isin(months)] # Select seasonal month range
if name == 'WY':
ds = s.groupby('water_year').mean(['time'])
avg = ds.sel(water_year = slice('1991', '2020')).mean()
depart = (ds - avg) / avg # Find percent departure from 1991-2020 average
else:
ds = s.groupby('time.year').mean(['time'])
avg = ds.sel(year = slice('1991', '2020')).mean()
depart = (ds - avg) / avg # Find percent departure from 1991-2020 average
#return(ds*120000) # Use for raw values of mm/decade
return(depart * 100 + 100)
ond = seasonal_trend(ds, [10,11,12], 'ond').sel(year = slice(1950,2021))
jfm = seasonal_trend(ds, [1, 2, 3], 'jfm').sel(year = slice(1951,2022))
amj = seasonal_trend(ds, [4, 5, 6], 'amj').sel(year = slice(1951,2022))
jas = seasonal_trend(ds, [7, 8, 9], 'jas').sel(year = slice(1951,2022))
wy = seasonal_trend(ds, [np.arange(1, 13)], 'WY').sel(water_year=slice(1951,2022))
# Check length
# ------------
seasons = [ond, jfm, amj, jas, wy]
labs = ['OND', 'JFM', 'AMJ', 'JAS', 'WY']
for sdx, sea in enumerate(seasons):
print(labs[sdx])
print(len(sea))
date = np.arange(1951, 2023)
df = pd.DataFrame({'year': date,
'OND': ond,
'JFM': jfm,
'AMJ': amj,
'JAS': jas,
'WY': wy})
df.to_csv('v4.ERA5.precipitation.%from1991-2020.csv') # v4 includes September 2022 data
exit()
# ===========================
# Plot Data
# ===========================
# (figure not included in 2022 Arctic Report Card)
fig, ax = plt.subplots(figsize=(10, 6))
seasons = [ond, jfm, amj, jas]#, year]
labs = ['OND', 'JFM', 'AMJ', 'JAS']#, 'Annual']
colors = ['brown', 'blue', 'pink', 'orange']#, 'black']
x = np.arange(1951, 2023)
for i, season in enumerate(seasons):
ax.plot(x, season,
label = labs[i],
color = colors[i])
ax.plot(x, wy,
label = 'WY',
color = 'black',
linewidth = 3)
ax.legend(loc='center', bbox_to_anchor =(1.05, 0.5))
ax.axhline(y = 100, color = 'black')
ax.spines.right.set_visible(False)
ax.spines.top.set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_ylabel('Percent of 1991-2020 average')
plt.grid(linestyle = '--')
#plt.savefig('weighted.arctic.precip.departure1991-2020.png', dpi=500)
plt.show()
exit()