Data for atmospheric CO$_2$ concentration¶
As greenhouse gas we consider only CO$_2$. But the increase of CO$_2$ in the atmosphere should correlate strongly with the emission of other greenhouse gases.
In [3]:
# data from http://www.esrl.noaa.gov/gmd/ccgg/trends/
# since 1959 direct measurements!
co2_concentration = pd.read_csv("ghg-concentrations_fig-1_.csv")
# set year as index
co2_concentration = co2_concentration.set_index("Year")
co2_concentration.tail(3)
Out[3]:
Use Pandas for computing the mean of each row. NaN values are ignored:
In [5]:
# set new mean series
co2_concentration['mean'] = co2_concentration.mean(axis=1)
In [6]:
co2_concentration.tail(3)
Out[6]:
In [7]:
plt.plot(co2_concentration.index, co2_concentration['mean'])
plt.xlim((1400, 2020))
plt.xlabel("Year")
plt.ylabel("CO$_2$ concentration / ppm")
Out[7]:
In [12]:
# convert c02 concentation to numpy
co2_year = co2_concentration.index.get_values()
co2_mean = co2_concentration['mean'].get_values()
Solar Irradiance¶
In [ ]:
lean_df = pd.read_csv("./lean2000_irradiance_only_data.csv", sep=',')
#lean_df
In [15]:
# solar irridiance
year_mask = (lean_df["YEAR"]- 0.5 >= df_air.Year.min())
scy = "11yrCYCLE+BKGRND"
#scy = "11yrCYCLE"
plt.plot(lean_df["YEAR"][year_mask],lean_df[scy][year_mask], "b*")
smo_ = np.ndarray(len(lean_df.YEAR))
import pandas.stats.moments
smoothed = pandas.stats.moments.ewma(lean_df[scy], halflife=3.)
plt.plot(lean_df["YEAR"][year_mask],smoothed[year_mask], "g-")
l = len(smoothed[year_mask])
smo_[:l] = smoothed[year_mask]
# according to a graph in http://www.climate4you.com/Sun.htm
# the irridience goes sligtly down sice 2000, to be more conservative we leave it constant
# from 2000-2014
smo_[l:] = smoothed[year_mask].max()
plt.xlabel("Year")
plt.ylabel("Solar irradiance / ?")
Out[15]:
Global air temperature¶
In [16]:
#mean = (df_air["J-D"]*0.01).mean()
temp_dev_air = -0.2 + ((df_air["J-D"]-df_air["J-D"].min())*0.01).get_values()
#temp_dev_lo = ((df_land_and_ocean["J-D"] - (df_land_and_ocean["J-D"]).min())*0.01).get_values()
In [18]:
plot_tempe()
In [26]:
plt.plot(co2_original, alo_temp, 'bo')
plt.xlabel("CO$_2$ / ppm")
plt.ylabel("$\Delta$ temperature / ${}^\circ$C")
Out[26]:
In [27]:
l = (year_mask).sum()
plt.plot(smoothed[year_mask], alo_temp[:l], 'bo')
plt.xlabel("solar irridiance / ?")
plt.ylabel("$\Delta$ temperature / ${}^\circ$C")
Out[27]:
In [28]:
import pymc3
model = pymc3.Model()
with model:
co_2_weight = pymc3.Normal("co_2_weight", mu=0, sd=1.)
irridiation_weight = pymc3.Normal("irridiaton_weight", mu=0, sd=1.)
offset = pymc3.Normal("offset", mu=0, sd=1.)
i_year = year - year[0]
i_co2_year = year - co2_year[0] # + co_2_delay
temp_mu = offset + irridiation_weight * smo[i_year] + co_2_weight * co2_
temp_sigma = pymc3.HalfNormal('sigma', sd=0.5)
temp_dev_ = pymc3.Normal("temp_dev_", mu=temp_mu, sd=temp_sigma, observed = alo_temp)
temp_pred = pymc3.Normal("temp_pred", mu=temp_mu, sd=temp_sigma, shape = alo_temp.shape )
In [30]:
with model:
# obtain starting values via MAP
start = pymc3.find_MAP()#fmin=optimize.fmin_powell)
#step = pymc3.Slice()
tr = pymc3.sample(5000, start=start)
In [32]:
plt.figure(figsize=(7, 7))
pymc3.traceplot(tr)
plt.tight_layout();
Compare the influence of irridiation and CO$_2$(Normalized Data):
In [34]:
irridiation_weight_.mean()
Out[34]:
In [35]:
co_2_weight_.mean()
Out[35]:
In [36]:
plt.figure(figsize=(12,10))
temp_pred_ = tr.get_values(temp_pred)[4000:]
for t in temp_pred_:
plt.plot(year, t, 'b.', alpha=0.02)
plt.xlabel("Year")
plt.ylabel("$\Delta T$ / ${}^\circ$C")
plt.plot(year, alo_temp, 'k.')
Out[36]:
In [37]:
plt.figure(figsize=(12,10))
for t in temp_.T:
plt.plot(year, t, 'b-', alpha=0.04)
plt.xlabel("Year")
plt.ylabel("$\Delta T$ / ${}^\circ$C")
plt.plot(year, alo_temp, 'k.')
Out[37]:
In [39]:
import sklearn
import sklearn.linear_model
lr = sklearn.linear_model.LinearRegression()
lr.fit(last_years[:,None], co2_last_years)
future_co2 = lr.predict(future[:,None])
Predicting the temperature increase¶
Assumption: The CO$_2$ level in atmosphere will (linear) increase as in the last years (from 1980 to now).
In [44]:
with pymc3.Model() as model_co2:
intercept = pymc3.Normal("intercept", mu=0, sd=100.)
l_weight = pymc3.Normal("l_weight", mu=0, sd=10.)
co2_mu = intercept + l_weight * last_years_shared
co2_sigma = pymc3.HalfNormal('co2_sigma', sd=0.2)
co2_pred = pymc3.Normal("co2_pred", mu=co2_mu, sd=co2_sigma, observed = co2_last_years)
co2_sample = pymc3.Normal("co2_sample", mu=co2_mu, sd=co2_sigma, shape = co2_last_years.shape)
In [45]:
with model_co2:
#step = pymc3.Slice()
trace_co2 = pymc3.sample(10000)
In [46]:
plt.figure(figsize=(7, 7))
pymc3.traceplot(trace_co2)
plt.tight_layout();
In [55]:
nb_samples=500
last_years_shared.set_value(future)
# Simply running PPC will use the updated values and do prediction
ppc = pymc3.sample_ppc(trace_co2[5000:], model=model_co2, samples=nb_samples)
In [56]:
co2_future = ppc["co2_pred"]
co2_future.shape
Out[56]:
In [80]:
plt.figure(figsize=(14,10))
for t in temp_[-len(temp_future_.T):].T:
plt.plot(year, t, 'b.', alpha=0.02)
for t in temp_future_.T:
plt.plot(future, t, 'r.', alpha=0.02)
plt.plot(year, alo_temp, 'k.')
plt.ylabel("$\delta T$ / C")
plt.xlabel("Year")
Out[80]:
Literature¶
- R. Nave: Modeling Human Impact on Global Warming
- Solar activity: http://www.climate4you.com/Sun.htm
- Global Warming: http://solar-center.stanford.edu/sun-on-earth/glob-warm.html