3
3
Calculating Seasonal Averages from Timeseries of Monthly Means
4
4
==============================================================
5
5
6
- Author: `Joe Hamman <http://www.hydro.washington.edu/~jhamman/ >`_
6
+ Author: `Joe Hamman <http://uw-hydro.github.io/current_member/joe_hamman/ >`_
7
+
8
+ The data for this example can be found in the `xray-data <https://github.com/xray/xray-data >`_ repository. This example is also available in an IPython Notebook that is available `here <https://github.com/xray/xray/tree/master/examples/xray_seasonal_means.ipynb >`_.
7
9
8
10
Suppose we have a netCDF or xray Dataset of monthly mean data and we
9
11
want to calculate the seasonal average. To do this properly, we need to
@@ -12,6 +14,7 @@ different number of days.
12
14
13
15
.. code :: python
14
16
17
+ % matplotlib inline
15
18
import numpy as np
16
19
import pandas as pd
17
20
import xray
@@ -24,9 +27,9 @@ different number of days.
24
27
25
28
.. parsed-literal ::
26
29
27
- numpy version : 1.9.1
28
- pandas version : 0.15 .2
29
- xray version : 0.4rc1-20-g52bbca3
30
+ numpy version : 1.9.2
31
+ pandas version : 0.16 .2
32
+ xray version : 0.5.1
30
33
31
34
32
35
Some calendar information so we can support any netCDF calendar.
@@ -42,6 +45,7 @@ Some calendar information so we can support any netCDF calendar.
42
45
' all_leap' : [0 , 31 , 29 , 31 , 30 , 31 , 30 , 31 , 31 , 30 , 31 , 30 , 31 ],
43
46
' 366_day' : [0 , 31 , 29 , 31 , 30 , 31 , 30 , 31 , 31 , 30 , 31 , 30 , 31 ],
44
47
' 360_day' : [0 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 ]}
48
+
45
49
A few calendar functions to determine the number of days in each month
46
50
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
47
51
@@ -80,59 +84,39 @@ the ``calendar.month_range`` function.
80
84
if leap_year(year, calendar = calendar):
81
85
month_length[i] += 1
82
86
return month_length
87
+
83
88
Open the ``Dataset ``
84
89
^^^^^^^^^^^^^^^^^^^^
85
90
86
91
.. code :: python
87
92
88
- monthly_mean_file = ' /raid2/jhamman/projects/RASM/data/processed/R1002RBRxaaa01a/lnd/monthly_mean_timeseries/R1002RBRxaaa01a.vic.hmm.197909-201212 .nc'
93
+ monthly_mean_file = ' RASM_example_data .nc'
89
94
ds = xray.open_dataset(monthly_mean_file, decode_coords = False )
90
- ds.attrs[' history' ] = ' ' # get rid of the history attribute because its obnoxiously long
91
95
print (ds)
92
96
97
+
93
98
.. parsed-literal ::
94
99
95
100
<xray.Dataset>
96
- Dimensions: (depth: 3, time: 400 , x: 275, y: 205)
101
+ Dimensions: ( time: 36 , x: 275, y: 205)
97
102
Coordinates:
98
- * time (time) datetime64[ns] 1979-09-16T12:00:00 1979-10-17 ...
99
- * depth (depth) int64 0 1 2
100
- * x (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ...
101
- * y (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ...
103
+ * time (time) datetime64[ns] 1980-09-16T12:00:00 1980-10-17 ...
104
+ * x (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...
105
+ * y (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...
102
106
Data variables:
103
- Precipitation (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
104
- Evap (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
105
- Runoff (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
106
- Baseflow (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
107
- Soilw (time, depth, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
108
- Swq (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
109
- Swd (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
110
- Swnet (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
111
- Lwnet (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
112
- Lwin (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
113
- Netrad (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
114
- Swin (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
115
- Latht (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
116
- Senht (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
117
- Grdht (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
118
- Albedo (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
119
- Radt (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
120
- Surft (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
121
- Relhum (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
122
- Tair (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
123
- Tsoil (time, depth, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
124
- Wind (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
107
+ Tair (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan ...
125
108
Attributes:
126
109
title: /workspace/jhamman/processed/R1002RBRxaaa01a/lnd/temp/R1002RBRxaaa01a.vic.ha.1979-09-01.nc
127
110
institution: U.W.
128
111
source: RACM R1002RBRxaaa01a
129
112
output_frequency: daily
130
113
output_mode: averaged
131
114
convention: CF-1.4
132
- history:
133
115
references: Based on the initial model of Liang et al., 1994, JGR, 99, 14,415- 14,429.
134
116
comment: Output from the Variable Infiltration Capacity (VIC) model.
135
117
nco_openmp_thread_number: 1
118
+ NCO: 4.3.7
119
+ history: history deleted for brevity
136
120
137
121
138
122
Now for the heavy lifting:
@@ -148,8 +132,10 @@ allong the time dimension.
148
132
.. code :: python
149
133
150
134
# Make a DataArray with the number of days in each month, size = len(time)
151
- month_length = xray.DataArray(get_dpm(ds.time.to_index(), calendar = ' noleap' ),
135
+ month_length = xray.DataArray(get_dpm(ds.time.to_index(),
136
+ calendar = ' noleap' ),
152
137
coords = [ds.time], name = ' month_length' )
138
+
153
139
# Calculate the weights by grouping by 'time.season'
154
140
weights = month_length.groupby(' time.season' ) / month_length.groupby(' time.season' ).sum()
155
141
@@ -158,49 +144,30 @@ allong the time dimension.
158
144
159
145
# Calculate the weighted average
160
146
ds_weighted = (ds * weights).groupby(' time.season' ).sum(dim = ' time' )
147
+
161
148
.. code :: python
162
149
163
150
print (ds_weighted)
164
151
152
+
165
153
.. parsed-literal ::
166
154
167
155
<xray.Dataset>
168
- Dimensions: (depth: 3, season: 4, x: 275, y: 205)
156
+ Dimensions: ( season: 4, x: 275, y: 205)
169
157
Coordinates:
170
- * y (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ...
171
- * x (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ...
172
- * depth (depth) int64 0 1 2
173
- * season (season) object 'DJF' 'JJA' 'MAM' 'SON'
158
+ * y (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...
159
+ * x (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...
160
+ * season (season) object 'DJF' 'JJA' 'MAM' 'SON'
174
161
Data variables:
175
- Baseflow (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
176
- Tsoil (season, depth, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
177
- Wind (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
178
- Swin (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
179
- Swq (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
180
- Netrad (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
181
- Albedo (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
182
- Evap (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
183
- Swd (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
184
- Radt (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
185
- Lwin (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
186
- Relhum (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
187
- Soilw (season, depth, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
188
- Lwnet (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
189
- Senht (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
190
- Surft (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
191
- Latht (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
192
- Runoff (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
193
- Tair (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
194
- Grdht (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
195
- Swnet (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
196
- Precipitation (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
162
+ Tair (season, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...
197
163
198
164
199
165
.. code :: python
200
166
201
167
# only used for comparisons
202
168
ds_unweighted = ds.groupby(' time.season' ).mean(' time' )
203
169
ds_diff = ds_weighted - ds_unweighted
170
+
204
171
.. code :: python
205
172
206
173
# Quick plot to show the results
0 commit comments