HadCRUT analysis

code.py

import numpy as np
import netCDF4 as nc
from scipy import stats

import pylab as plt
import seaborn as sns

rootgrp = nc.Dataset('HadCRUT.4.6.0.0.median.nc')  

t = rootgrp['time'][:] / 365 + 1850
temp = np.nanmean(rootgrp['temperature_anomaly'][:], axis=(1,2))
rate = np.diff(temp) / np.diff(t)
t_diff = t[:-1]

plt.figure()
plt.plot(t_diff, rate, label='Raw data, month resolution', alpha=0.2)
plt.xlim(1959, 2019)
plt.ylim(-0.04, 0.04)
plt.title('Reproduction of figure')
plt.xlabel('Year')
plt.ylabel('Rate (º/yr)')

poly1 = np.polyfit(t_diff, rate, 2)
poly2 = np.polyfit(t_diff[t_diff > 1959], rate[t_diff > 1959], 2)

x = np.linspace(1959, 2019, num=100)
plt.plot(x, np.polyval(poly1, x), label='Full fit')
plt.plot(x, np.polyval(poly2, x), label='Fit on 1959-')

plt.legend(loc=0)



slope, _, r, p_value, err = stats.linregress(t_diff[t_diff > 2000], rate[t_diff > 2000])
print(slope, err, p_value)
plt.figure('reg')
ax_reg = plt.subplot(211)
sns.regplot(t_diff[t_diff > 2000], rate[t_diff > 2000], label='Unbinned (monthly)')  
plt.legend(loc=0)
plt.ylim(-0.6, 0.6)
plt.title('Last trend')
plt.ylabel('Rate (º/yr)')


t = rootgrp['time'][:] / 365 + 1850
temp = np.nanmean(rootgrp['temperature_anomaly'][:], axis=(1,2))
temp, t_, _ = stats.binned_statistic(t, temp, bins=t.ptp())
rate = np.diff(temp) / np.diff(t_[:-1])
t_diff = t_[:-2]

plt.figure()
plt.plot(t_diff, rate, label='Raw data, year resolution', alpha=0.2)
plt.xlim(1959, 2019)
plt.ylim(-0.04, 0.04)
plt.title('Reproduction of figure')
plt.xlabel('Year')
plt.ylabel('Rate (º/yr)')

poly1 = np.polyfit(t_diff, rate, 2)
poly2 = np.polyfit(t_diff[t_diff > 1959], rate[t_diff > 1959], 2)

x = np.linspace(1959, 2019, num=100)
plt.plot(x, np.polyval(poly1, x), label='Full fit')
plt.plot(x, np.polyval(poly2, x), label='Fit on 1959-')

plt.legend(loc=0)


slope, _, r, p_value, err = stats.linregress(t_diff[t_diff > 2000], rate[t_diff > 2000])
print(slope, err, p_value)

plt.figure('reg')
plt.subplot(212, sharex=ax_reg)
sns.regplot(t_diff[t_diff > 2000], rate[t_diff > 2000], label='Binned (yearly)')  
plt.ylim(-0.6, 0.6)
plt.xlabel('Year')
plt.ylabel('Rate (º/yr)')
plt.legend(loc=0)
plt.tight_layout()
plt.show()