In [1]:
# NBVAL_IGNORE_OUTPUT
from os.path import join
from datetime import datetime
import pandas as pd
from pymagicc import MAGICC6
from pymagicc.io import MAGICCData
from pymagicc.definitions import (
convert_magicc7_to_openscm_variables,
PART_OF_SCENFILE_WITH_EMISSIONS_CODE_1
)
import expectexception
import matplotlib.pyplot as plt
plt.style.use('bmh')
%matplotlib inline
In [2]:
DATA_DIR = join("..", "tests", "test_data")
DATA_MAGICC6_RUN = join("..", "pymagicc", "MAGICC6", "run")
We start by reading in our raw data. As the formats we're reading don't have any easily readable metadata, we set it using the columns argument.
Note: The model column refers to the model which generated the scenario (typically an integrated assessment model), with the climate_model column being saved for storing the climate model which was used to do climate projections.
In [3]:
rcp3pd = MAGICCData(
join(DATA_DIR, "RCP3PD_EMISSIONS.DAT"),
columns={
"model": "IMAGE",
"scenario": "RCP3PD"
}
)
rcp3pd.head()
Out[3]:
This data came from MAGICC6.3.09, in order to set it up for a new run, we have to strip out that information.
In [4]:
rcp3pd.set_meta("unspecified", "climate_model")
rcp3pd.head()
Out[4]:
In [5]:
rcp85 = MAGICCData(
join(DATA_MAGICC6_RUN, "RCP85.SCEN"),
columns={
"model": "MESSAGE",
"scenario": "RCP85"
}
)
rcp85.set_meta("unspecified", "climate_model")
rcp85.head()
Out[5]:
To join the data, we first have to do a little bit of preparation. Specifically, we must make sure that the metadata for the timeseries we want to prepare is the same which means re-writing the model and scenario columns.
In [6]:
my_scen_hist = rcp3pd.copy()
my_scen_hist.set_meta("idealised", "model")
my_scen_hist.set_meta("custom", "scenario")
my_scen_future = rcp85.copy()
my_scen_future.set_meta("idealised", "model")
my_scen_future.set_meta("custom", "scenario")
my_scen_future.tail()
Out[6]:
To then join the data, we can simply filter the data and then append one onto the other.
In [7]:
# NBVAL_IGNORE_OUTPUT
my_scen_hist.filter(year=range(1, 2006)).head()
Out[7]:
In [8]:
# NBVAL_IGNORE_OUTPUT
rcp3pd_up_to_2005 = my_scen_hist.filter(year=range(1, 2006))
rcp85_2006_onwards = my_scen_future.filter(year=range(2006, 30000))
# TODO: fix openscm method so it checks for duplicate time points, not duplicate index...
my_scen = rcp3pd_up_to_2005.append(rcp85_2006_onwards, duplicate_msg="warn")
my_scen["time"].head()
Out[8]:
In [9]:
# NBVAL_IGNORE_OUTPUT
my_scen["time"].tail()
Out[9]:
Note that appending in this way keeps all of the available data. This means that we still have regional data over the future.
In [10]:
# NBVAL_IGNORE_OUTPUT
my_scen.filter(region="World", keep=False).head()
Out[10]:
We can remove this by simply filtering our scenario once again.
In [11]:
# NBVAL_IGNORE_OUTPUT
my_scen_world = my_scen.filter(region="World")
my_scen_world.head()
Out[11]:
We can check the results as shown.
In [12]:
var_to_plot = "Emissions|CO2|MAGICC Fossil and Industrial"
ax = rcp3pd.filter(
variable=var_to_plot,
region="World"
).line_plot(x="time", label="RCP3PD", lw=6.0, figsize=(16, 9))
rcp85.filter(
variable=var_to_plot,
region="World"
).line_plot(x="time", label="RCP85", lw=6.0, ax=ax)
my_scen_world.filter(
variable=var_to_plot
).line_plot(x="time", label="my_scen", lw=6.0, linestyle=":", ax=ax)
ax.set_xlim([datetime(2000, 1, 1), datetime(2011, 1, 1)])
ax.set_ylim([6.5, 9]);
If we put our scenario into a MAGICCData object, we can run it using Pymagicc.
Note that Pymagicc automatically only uses the variables which are in SCEN files.
In [13]:
with MAGICC6() as magicc:
results = magicc.run(my_scen_world)
In [14]:
results.filter(
variable="Surface Temperature",
region="World"
).line_plot(x="time", figsize=(16, 9));
If we wish, we can also MAGICCData's resample method to join our timeseries together with a linear interpolation over a given period. This allows us to use some base scenario, overwrite with a multitude of other scenarios but ensure they all have a common historical period and some roughly smooth transition. (At the moment the resample and interpolate methods are the only options, a far more complete solution would be to use the Aneris package but this would also be far more difficult.)
In [15]:
rcp3pd_up_to_2005 = my_scen_hist.filter(year=range(1, 2006))
rcp85_2015_onwards = my_scen_future.filter(year=range(2015, 30000))
# TODO: fix openscm method so it checks for duplicate time points, not duplicate index...
# TODO: work out why resample and interpolate are so slow...
target_times = [datetime(y, 1, 1) for y in range(2000, 2501)]
linear_scen = rcp3pd_up_to_2005.append(
rcp85_2015_onwards, duplicate_msg="warn"
).filter(region="World")#.interpolate(target_times)
In [16]:
var_to_plot = "Emissions|CO2|MAGICC Fossil and Industrial"
ax = rcp3pd.filter(
variable=var_to_plot,
region="World"
).line_plot(x="time", label="RCP3PD", lw=6.0, figsize=(16, 9))
rcp85.filter(
variable=var_to_plot,
region="World"
).line_plot(x="time", label="RCP85", lw=6.0, ax=ax)
linear_scen.filter(
variable=var_to_plot
).line_plot(x="time", label="my_scen", lw=6.0, linestyle=":", ax=ax)
ax.set_xlim([datetime(2000, 1, 1), datetime(2031, 1, 1)])
ax.set_ylim([6.5, 15]);
In [17]:
with MAGICC6() as magicc:
results = magicc.run(linear_scen)
In [18]:
results.filter(
variable="Surface Temperature",
region="World"
).line_plot(x="time", figsize=(16, 9));
If you try to merge data which has no overlapping indices, you'll end up with np.nan in your output. For example, if we set our "todo" column to "N/A" in base (for an explanation of what this "todo" column is, see file conventions).
As a user, you have to check carefully that the merging has actually happened as intended (pull requests with convenience functions that can help automate these checks are most welcome!!). One easy way to do this is to simply see if any of the data in the output is null.
In [19]:
rcp85_2015_onwards_world = my_scen_future.filter(
year=range(2015, 30000),
region="World",
)
rcp3pd_up_to_2005_world = my_scen_hist.filter(
year=range(1, 2006),
region="World",
)
null_output = rcp3pd_up_to_2005_world.append(rcp85_2015_onwards_world, duplicate_msg="warn").timeseries()
null_output.isnull().any().any()
Out[19]:
There are still null values, we can see where they are with something like this.
In [20]:
null_output[null_output.isnull().any(axis=1)]
Out[20]:
This makes it clear that all of the non-SCEN gases have np.nan values into the future, which makes sense as the SCEN file does not define them.
If we drop these variables first, i.e. do the merge carefully, none of the null values appear.
In [21]:
rcp85_2015_onwards_world = my_scen_future.filter(year=range(2015, 30000), region="World")
rcp3pd_up_to_2005_world_careful = my_scen_hist.filter(
year=range(1, 2006),
region="World",
variable=rcp85_2015_onwards_world["variable"],
)
no_nulls = rcp3pd_up_to_2005_world_careful.append(rcp85_2015_onwards_world, duplicate_msg="warn")
no_nulls.timeseries().isnull().any().any()
Out[21]:
As a side effect, running MAGICC with this scenario does not produce any warning about ignoring data.
In [22]:
with MAGICC6() as magicc:
results = magicc.run(linear_scen)