%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
precip = xr.open_dataset('http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ghcngridded/precip.mon.total.nc')
precip
Let's look at how the distribution of Tropical precipitation rates have changed over time. Let's group the data by decade, and plot a distribution of the point by point precipitation rates within that region.
First, let's subset the data by space.
precip_space_subset=precip['precip'].sel(lat=slice(20 -20)) #go in ame order as the data when slicing
precip_space_subset
precip_time_subset=precip_space_subset.sel(time=slice('1900-01-01', '1909-12-31'))
precip_time_subset
#print(precip_time_subset['time'])
OK, now that we know how to create a subset, let's do the statistical analysis of this field. First, let's create a histogram of all the values (note we are not doing any spatial or time averaging, just subsetting), and we will set up the bins in the histogram to have 31 values, 30 mm wide between 0 and 900 mm. These are the bin edges. The plt.hist function will count the values between the bin edges, and display them as a plot. Then we will create a PDF by enabling normed=True, and a CDF by normed=True and cumulative=True, respectively. Note that the returned value has the values, and the bin edges.
hist_1910=xr.plot.hist(precip_time_subset,bins=np.arange(0,2000,30))
hist_1910
pdf_1910=xr.plot.hist(precip_time_subset,bins=30*np.arange(31),normed=True)
pdf_1910
cdf_1910=xr.plot.hist(precip_time_subset, bins=30*np.arange(31), normed=True, cumulative=True)
cdf_1910
Note that the distribution is far from "normal" - the left tail of the distribution is "heavier" than the right tail. This physically corresponds with a higher probability of light rain than heavy rain in the observations which makes sense. Examine how the PDF and CDF looks in this case corresponding with this distribution.
Now, let's process all of the decades. I'll set up a for loop to process through. We will use pandas handy built in date functionality which is nice for selecting date periods automatically. Decades can be selected using 10AS as the frequency string, and we'll select 11 decades starting in 1900.
import pandas as pd
time_start = pd.date_range('1900-01-01', freq='10AS', periods=11)
time_end = pd.date_range('1909-12-31', freq='10A', periods=11)
print(time_start)
print(time_end)
Now, it is just a matter of looping over the decades, storing the data into arrays, and then plotting. We will set up a numpy array to store the data from the histograms. Instead of plt.hist i'll use np.histogram, which is very similar, because I don't want to plot the data yet.
ndecades=len(time_start)
#bins for histograms
bins=30*np.arange(31)
#What is going on here? Test!
bin_centers=bins[0:-1]+np.diff(bins)/2.
nbins=len(bins)
#need to subtract 1 from nbins because bins are the edges
#and are 1 more value than the histogram
histogram_decades=np.zeros((nbins-1,ndecades))
for i in np.arange(ndecades):
data_this_decade=precip_space_subset.sel(time=slice(time_start[i], time_end[i]))
histogram_decades[:,i]=xr.plot.hist(data_this_decade,bins=bins)[0]
histogram_decades
plt.figure()
plt.plot(bin_centers,histogram_decades)
plt.xlabel('Precipitation (mm/mon)')
plt.ylabel('Histogram counts')
plt.legend(time_start,fontsize=8)
plt.show()
Interestingly, it looks like there are more lications with stations in the mid-20th century than now in the Tropics. This has been associated with a decline in the observational network in Africa.
That takes care of the histogram, but what about the PDF?
The definition of a PDF is the sum of the histogram divided by the total of the histogram, which yields the fraction or probability of obtaining a data in that range out of the population of datapoints. We can divide the histogram by the sum on axis=0 (which will sum each decadal pdf) to get the PDF.
pdf_decades=histogram_decades/np.sum(histogram_decades,axis=0)
plt.figure()
plt.plot(bin_centers,pdf_decades)
plt.xlabel('Precipitation (mm/mon)')
plt.ylabel('PDF')
plt.legend(time_start,fontsize=8)
plt.show()
To get the CDF, we can use the cumsum object, which will sum the PDF over a specified axis, in this case we want axis 0, which are the precipitation rates.
cdf_decades=pdf_decades.cumsum(axis=0)
plt.figure()
plt.plot(bin_centers,cdf_decades)
plt.xlabel('Precipitation (mm/mon)')
plt.xlabel('CDF')
plt.legend(time_start,fontsize=8,loc='lower right')
plt.show()
OK, we see that the statistical distribution of observed Tropical rain rates is quite similar over the decades. We have different numbers of stations, but the distribution of rain rates at those stations hasn't changed much.
For you homework, you need to do the statistics when subsetting over space. It is quite similar to what we did here.