Skip to content

Commit 3eb015b

Browse files
shoyershoyer
committedFeb 26, 2015
Merge pull request #337 from shoyer/cleanup
Cleanup (mostly documentation)
·
v2025.03.0v0.4rc1
2 parents 7ea9de8 + ca14eb7 commit 3eb015b

13 files changed

+400
-87
lines changed
 

‎doc/api.rst

+2-1
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,7 @@ Dataset contents
6565
Dataset.copy
6666
Dataset.merge
6767
Dataset.rename
68-
Dataset.drop_vars
68+
Dataset.drop
6969
Dataset.set_coords
7070
Dataset.reset_coords
7171

@@ -164,6 +164,7 @@ DataArray contents
164164
:toctree: generated/
165165

166166
DataArray.rename
167+
DataArray.drop
167168
DataArray.reset_coords
168169
DataArray.copy
169170

‎doc/combining.rst

+6-6
Original file line numberDiff line numberDiff line change
@@ -70,8 +70,8 @@ necessary.
7070

7171
.. _merge:
7272

73-
Merge and ``Dataset.__init__``
74-
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
73+
Merge
74+
~~~~~
7575

7676
To combine variables and coordinates between multiple Datasets, you can use the
7777
:py:meth:`~xray.Dataset.merge` and :py:meth:`~xray.Dataset.update` methods.
@@ -102,8 +102,8 @@ used in the :py:class:`~xray.Dataset` constructor:
102102
103103
.. _update:
104104

105-
Update and ``__setitem__``
106-
~~~~~~~~~~~~~~~~~~~~~~~~~~
105+
Update
106+
~~~~~~
107107

108108
In contrast, update modifies a dataset in-place without checking for conflicts,
109109
and will overwrite any existing variables with new values:
@@ -159,8 +159,8 @@ are constant along those new dimensions:
159159

160160
.. ipython:: python
161161
162-
left = Dataset(coords={'x': 0})
163-
right = Dataset({'x': [0, 0, 0]})
162+
left = xray.Dataset(coords={'x': 0})
163+
right = xray.Dataset({'x': [0, 0, 0]})
164164
left.broadcast_equals(right)
165165
166166
Like pandas objects, two xray objects are still equal or identical if they have

‎doc/examples.rst

+1
Original file line numberDiff line numberDiff line change
@@ -6,3 +6,4 @@ Examples
66

77
examples/quick-overview
88
examples/weather-data
9+
examples/monthly-means

‎doc/examples/monthly-means.rst

+260
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,260 @@
1+
Calculating Seasonal Averages from Timeseries of Monthly Means
2+
==============================================================
3+
4+
Author: `Joe Hamman <http://www.hydro.washington.edu/~jhamman/>`_
5+
6+
Suppose we have a netCDF or xray Dataset of monthly mean data and we
7+
want to calculate the seasonal average. To do this properly, we need to
8+
calculate the weighted average considering that each month has a
9+
different number of days.
10+
11+
.. code:: python
12+
13+
import numpy as np
14+
import pandas as pd
15+
import xray
16+
from netCDF4 import num2date
17+
import matplotlib.pyplot as plt
18+
19+
print("numpy version : ", np.__version__)
20+
print("pandas version : ", pd.version.version)
21+
print("xray version : ", xray.version.version)
22+
23+
.. parsed-literal::
24+
25+
numpy version : 1.9.1
26+
pandas version : 0.15.2
27+
xray version : 0.3.2
28+
29+
30+
Some calendar information so we can support any netCDF calendar.
31+
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
32+
33+
.. code:: python
34+
35+
dpm = {'noleap': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
36+
'365_day': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
37+
'standard': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
38+
'gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
39+
'proleptic_gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
40+
'all_leap': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
41+
'366_day': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
42+
'360_day': [0, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30]}
43+
A few calendar functions to determine the number of days in each month
44+
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
45+
46+
If you were just using the standard calendar, it would be easy to use
47+
the ``calendar.month_range`` function.
48+
49+
.. code:: python
50+
51+
def leap_year(year, calendar='standard'):
52+
"""Determine if year is a leap year"""
53+
leap = False
54+
if ((calendar in ['standard', 'gregorian',
55+
'proleptic_gregorian', 'julian']) and
56+
(year % 4 == 0)):
57+
leap = True
58+
if ((calendar == 'proleptic_gregorian') and
59+
(year % 100 == 0) and
60+
(year % 400 != 0)):
61+
leap = False
62+
elif ((calendar in ['standard', 'gregorian']) and
63+
(year % 100 == 0) and (year % 400 != 0) and
64+
(year < 1583)):
65+
leap = False
66+
return leap
67+
68+
def get_dpm(time, calendar='standard'):
69+
"""
70+
return a array of days per month corresponding to the months provided in `months`
71+
"""
72+
month_length = np.zeros(len(time), dtype=np.int)
73+
74+
cal_days = dpm[calendar]
75+
76+
for i, (month, year) in enumerate(zip(time.month, time.year)):
77+
month_length[i] = cal_days[month]
78+
if leap_year(year, calendar=calendar):
79+
month_length[i] += 1
80+
return month_length
81+
Open the ``Dataset``
82+
^^^^^^^^^^^^^^^^^^^^
83+
84+
.. code:: python
85+
86+
monthly_mean_file = '/raid2/jhamman/projects/RASM/data/processed/R1002RBRxaaa01a/lnd/monthly_mean_timeseries/R1002RBRxaaa01a.vic.hmm.197909-201212.nc'
87+
ds = xray.open_dataset(monthly_mean_file, decode_coords=False)
88+
ds.attrs['history'] = '' # get rid of the history attribute because its obnoxiously long
89+
print(ds)
90+
91+
.. parsed-literal::
92+
93+
<xray.Dataset>
94+
Dimensions: (depth: 3, time: 400, x: 275, y: 205)
95+
Coordinates:
96+
* time (time) datetime64[ns] 1979-09-16T12:00:00 1979-10-17 ...
97+
* depth (depth) int64 0 1 2
98+
* 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 ...
99+
* 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 ...
100+
Variables:
101+
Precipitation (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
102+
Evap (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
103+
Runoff (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
104+
Baseflow (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
105+
Soilw (time, depth, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
106+
Swq (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
107+
Swd (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
108+
Swnet (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
109+
Lwnet (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
110+
Lwin (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
111+
Netrad (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
112+
Swin (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
113+
Latht (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
114+
Senht (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
115+
Grdht (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
116+
Albedo (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
117+
Radt (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
118+
Surft (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
119+
Relhum (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
120+
Tair (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
121+
Tsoil (time, depth, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
122+
Wind (time, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan nan nan ...
123+
Attributes:
124+
title: /workspace/jhamman/processed/R1002RBRxaaa01a/lnd/temp/R1002RBRxaaa01a.vic.ha.1979-09-01.nc
125+
institution: U.W.
126+
source: RACM R1002RBRxaaa01a
127+
output_frequency: daily
128+
output_mode: averaged
129+
convention: CF-1.4
130+
history:
131+
references: Based on the initial model of Liang et al., 1994, JGR, 99, 14,415- 14,429.
132+
comment: Output from the Variable Infiltration Capacity (VIC) model.
133+
nco_openmp_thread_number: 1
134+
135+
136+
Now for the heavy lifting:
137+
^^^^^^^^^^^^^^^^^^^^^^^^^^
138+
139+
We first have to come up with the weights, - calculate the month lengths
140+
for each monthly data record - calculate weights using
141+
``groupby('time.season')``
142+
143+
Finally, we just need to multiply our weights by the ``Dataset`` and sum
144+
allong the time dimension.
145+
146+
.. code:: python
147+
148+
# Make a DataArray with the number of days in each month, size = len(time)
149+
month_length = xray.DataArray(get_dpm(ds.time.to_index(), calendar='noleap'),
150+
coords=[ds.time], name='month_length')
151+
# Calculate the weights by grouping by 'time.season'
152+
weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()
153+
154+
# Test that the sum of the weights for each season is 1.0
155+
np.testing.assert_allclose(weights.groupby('time.season').sum().values, np.ones(4))
156+
157+
# Calculate the weighted average
158+
ds_weighted = (ds * weights).groupby('time.season').sum(dim='time')
159+
.. code:: python
160+
161+
print(ds_weighted)
162+
163+
.. parsed-literal::
164+
165+
<xray.Dataset>
166+
Dimensions: (depth: 3, time.season: 4, x: 275, y: 205)
167+
Coordinates:
168+
* depth (depth) int64 0 1 2
169+
* 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 ...
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+
* time.season (time.season) int32 1 2 3 4
172+
Variables:
173+
Lwnet (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
174+
Tair (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
175+
Surft (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
176+
Senht (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
177+
Tsoil (time.season, depth, y, x) float64 nan nan nan nan nan nan nan nan nan ...
178+
Netrad (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
179+
Evap (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
180+
Latht (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
181+
Wind (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
182+
Precipitation (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
183+
Soilw (time.season, depth, y, x) float64 nan nan nan nan nan nan nan nan nan ...
184+
Relhum (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
185+
Swd (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
186+
Swnet (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
187+
Swq (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
188+
Swin (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
189+
Albedo (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
190+
Lwin (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
191+
Radt (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
192+
Runoff (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
193+
Grdht (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
194+
Baseflow (time.season, y, x) float64 nan nan nan nan nan nan nan nan nan nan nan ...
195+
196+
197+
.. code:: python
198+
199+
# only used for comparisons
200+
ds_unweighted = ds.groupby('time.season').mean('time')
201+
ds_diff = ds_weighted - ds_unweighted
202+
.. code:: python
203+
204+
# Quick plot to show the results
205+
is_null = np.isnan(ds_weighted['Tair'][0].values)
206+
207+
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(14,12))
208+
for i, season in enumerate(('DJF', 'MAM', 'JJA', 'SON')):
209+
plt.sca(axes[i, 0])
210+
plt.pcolormesh(np.ma.masked_where(is_null, ds_weighted['Tair'][i].values),
211+
vmin=-30, vmax=30, cmap='Spectral_r')
212+
plt.colorbar(extend='both')
213+
214+
plt.sca(axes[i, 1])
215+
plt.pcolormesh(np.ma.masked_where(is_null, ds_unweighted['Tair'][i].values),
216+
vmin=-30, vmax=30, cmap='Spectral_r')
217+
plt.colorbar(extend='both')
218+
219+
plt.sca(axes[i, 2])
220+
plt.pcolormesh(np.ma.masked_where(is_null, ds_diff['Tair'][i].values),
221+
vmin=-0.1, vmax=.1, cmap='RdBu_r')
222+
plt.colorbar(extend='both')
223+
for j in range(3):
224+
axes[i, j].axes.get_xaxis().set_ticklabels([])
225+
axes[i, j].axes.get_yaxis().set_ticklabels([])
226+
axes[i, j].axes.axis('tight')
227+
228+
axes[i, 0].set_ylabel(season)
229+
230+
axes[0, 0].set_title('Weighted by DPM')
231+
axes[0, 1].set_title('No Weighting')
232+
axes[0, 2].set_title('Difference')
233+
234+
plt.tight_layout()
235+
236+
fig.suptitle('Seasonal Surface Air Temperature', fontsize=16, y=1.02)
237+
238+
239+
.. image:: monthly_means_output.png
240+
241+
242+
.. code:: python
243+
244+
# Wrap it into a simple function
245+
def season_mean(ds, calendar='standard'):
246+
# Make a DataArray of season/year groups
247+
year_season = xray.DataArray(ds.time.to_index().to_period(freq='Q-NOV').to_timestamp(how='E'),
248+
coords=[ds.time], name='year_season')
249+
250+
# Make a DataArray with the number of days in each month, size = len(time)
251+
month_length = xray.DataArray(get_dpm(ds.time.to_index(), calendar=calendar),
252+
coords=[ds.time], name='month_length')
253+
# Calculate the weights by grouping by 'time.season'
254+
weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()
255+
256+
# Test that the sum of the weights for each season is 1.0
257+
np.testing.assert_allclose(weights.groupby('time.season').sum().values, np.ones(4))
258+
259+
# Calculate the weighted average
260+
return (ds * weights).groupby('time.season').sum(dim='time')

‎doc/examples/monthly_means_output.png

435 KB
Loading

0 commit comments

Comments
 (0)
Please sign in to comment.