How to make a step forcing¶

Here we explain how to make a sharp step forcing.

Imports¶

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import openscm_units

import continuous_timeseries as ct

import matplotlib.pyplot as plt import numpy as np import openscm_units import continuous_timeseries as ct

Handy pint aliases¶

In [2]:

UR = openscm_units.unit_registry
Q = UR.Quantity

UR = openscm_units.unit_registry Q = UR.Quantity

Set up matplotlib to work with pint¶

For details, see the pint docs (stable docs, last version that we checked at the time of writing) or our docs on unit-aware plotting.

In [3]:

UR.setup_matplotlib(enable=True)

UR.setup_matplotlib(enable=True)

Introducing the problem¶

Imagine you want to run an experiment with a step forcing. For example, forcing before 1850 is 0 W / m² while forcing from 1850 (inclusive) onwards is equal to 2 W / m².

The naive solution¶

One way you could achieve this is to simply use a timeseries with a really fine time axis.

In [4]:

time_fine = Q(np.linspace(1750, 2100, int(1e5)), "yr")
time_fine

time_fine = Q(np.linspace(1750, 2100, int(1e5)), "yr") time_fine

Out[4]:

Magnitude	[1750.0 1750.0035000350003 1750.0070000700007 ... 2099.9929999299993 2099.9964999649997 2100.0]
Units	yr

In [5]:

forcing_fine = np.zeros_like(time_fine) * UR.Unit("W / m^2")
transition_year = 1850
forcing_fine[time_fine.m >= transition_year] = Q(2.0, "W / m^2")
forcing_fine

forcing_fine = np.zeros_like(time_fine) * UR.Unit("W / m^2") transition_year = 1850 forcing_fine[time_fine.m >= transition_year] = Q(2.0, "W / m^2") forcing_fine

Out[5]:

Magnitude	[0.0 0.0 0.0 ... 2.0 2.0 2.0]
Units	watt/meter2

This works, but if you zoom in far enough, you see the discrete steps i.e. the edge is no longer sharp. In addition, you have to carry around an array with heaps of time points and whatever API you pass this to has to know how to deal with a fine array of points.

In [6]:

fig, axes = plt.subplots(ncols=4, sharey=True, figsize=(12, 4))

for i, xlim in enumerate(
    ((1750, 2100), (1849, 1851), (1849.9, 1850.1), (1849.99, 1850.01))
):
    axes[i].plot(time_fine, forcing_fine)
    axes[i].set_xlim(xlim)
    axes[i].grid()

fig.tight_layout()

fig, axes = plt.subplots(ncols=4, sharey=True, figsize=(12, 4)) for i, xlim in enumerate( ((1750, 2100), (1849, 1851), (1849.9, 1850.1), (1849.99, 1850.01)) ): axes[i].plot(time_fine, forcing_fine) axes[i].set_xlim(xlim) axes[i].grid() fig.tight_layout()

The Continuous Timeseries solution¶

With continuous timeseries, we can achieve the same thing with far fewer data points.

In [7]:

ts = ct.Timeseries.from_arrays(
    y=Q([0.0, 2.0, 2.0], "W / m^2"),
    x=Q([1750, 1850, 3000], "yr"),
    interpolation=ct.InterpolationOption.PiecewiseConstantPreviousLeftClosed,
    name="step_forcing",
)
ts

ts = ct.Timeseries.from_arrays( y=Q([0.0, 2.0, 2.0], "W / m^2"), x=Q([1750, 1850, 3000], "yr"), interpolation=ct.InterpolationOption.PiecewiseConstantPreviousLeftClosed, name="step_forcing", ) ts

Out[7]:

continuous_timeseries.Timeseries

time_axis	bounds Magnitude [1750 1850 3000] Units yr
timeseries_continuous	name step_forcing time_units yr values_units watt/meter2 function PPolyPiecewiseConstantPreviousLeftClosed(x=array([1750, 1850, 3000]), y=array([0., 2., 2.])) domain (1750 yr, 3000 yr)

This solution also maintains the sharp edge, no matter how far we zoom in.

In [8]:

fig, axes = plt.subplots(ncols=4, sharey=True, figsize=(12, 4))

for i, xlim in enumerate(
    ((1750, 2100), (1849, 1851), (1849.9, 1850.1), (1849.99, 1850.01))
):
    ts.plot(
        ax=axes[i],
        # Make sure that the plot resolution is high enough
        # to show what is actually going on.
        continuous_plot_kwargs=dict(res_increase=int(3e6)),
    )
    axes[i].set_xlim(xlim)
    axes[i].grid()

fig.tight_layout()

fig, axes = plt.subplots(ncols=4, sharey=True, figsize=(12, 4)) for i, xlim in enumerate( ((1750, 2100), (1849, 1851), (1849.9, 1850.1), (1849.99, 1850.01)) ): ts.plot( ax=axes[i], # Make sure that the plot resolution is high enough # to show what is actually going on. continuous_plot_kwargs=dict(res_increase=int(3e6)), ) axes[i].set_xlim(xlim) axes[i].grid() fig.tight_layout()

In [9]:

# No matter how fine we make the time axis,
# the sharp transition is maintained.
ts.interpolate(Q([1850 - 1e-9, 1850.0, 1850 + 1e-9], "yr")).discrete

# No matter how fine we make the time axis, # the sharp transition is maintained. ts.interpolate(Q([1850 - 1e-9, 1850.0, 1850 + 1e-9], "yr")).discrete

Out[9]:

continuous_timeseries.TimeseriesDiscrete

name	step_forcing
time_axis	bounds Magnitude [1849.999999999 1850.0 1850.000000001] Units yr
values_at_bounds	values Magnitude [0.0 2.0 2.0] Units watt/meter2

A few other notes¶

The way ts is created above, at 1850, the value is 2.0 W / m². If you want the value at 1850 to be zero instead, you can use PPolyPiecewiseConstantPreviousLeftOpen interpolation as shown below.

In [10]:

ts_open = ct.Timeseries.from_arrays(
    x=Q([1750, 1850, 3000], "yr"),
    y=Q([0.0, 2.0, 2.0], "W / m^2"),
    interpolation=ct.InterpolationOption.PiecewiseConstantPreviousLeftOpen,
    name="step_forcing",
)
ts_open.plot()

ts_open = ct.Timeseries.from_arrays( x=Q([1750, 1850, 3000], "yr"), y=Q([0.0, 2.0, 2.0], "W / m^2"), interpolation=ct.InterpolationOption.PiecewiseConstantPreviousLeftOpen, name="step_forcing", ) ts_open.plot()

Out[10]:

<Axes: xlabel='yr', ylabel='watt/meter²'>

In [11]:

# Notice here that the 1850 value is 0.0 W / m^2, not 2.0 W / m^2
ts_open.discrete

# Notice here that the 1850 value is 0.0 W / m^2, not 2.0 W / m^2 ts_open.discrete

Out[11]:

continuous_timeseries.TimeseriesDiscrete

name	step_forcing
time_axis	bounds Magnitude [1750 1850 3000] Units yr
values_at_bounds	values Magnitude [0.0 0.0 2.0] Units watt/meter2

You can also achieve something similar using PPolyPiecewiseConstantNextLeftClosed and PPolyPiecewiseConstantNextLeftOpen. The process is slightly different but the outcome is the largely the same, with the main difference being in interpolation (for a further overview of the differences, see the docs of InterpolationOption).

In [12]:

ts_next_left_closed = ct.Timeseries.from_arrays(
    x=Q([1750, 1850, 3000], "yr"),
    y=Q([0.0, 0.0, 2.0], "W / m^2"),
    interpolation=ct.InterpolationOption.PiecewiseConstantNextLeftClosed,
    name="next_left_closed",
)
ts_next_left_closed.plot()

ts_next_left_closed = ct.Timeseries.from_arrays( x=Q([1750, 1850, 3000], "yr"), y=Q([0.0, 0.0, 2.0], "W / m^2"), interpolation=ct.InterpolationOption.PiecewiseConstantNextLeftClosed, name="next_left_closed", ) ts_next_left_closed.plot()

Out[12]:

<Axes: xlabel='yr', ylabel='watt/meter²'>

In [13]:

ts_next_left_closed.discrete

ts_next_left_closed.discrete

Out[13]:

continuous_timeseries.TimeseriesDiscrete

name	next_left_closed
time_axis	bounds Magnitude [1750 1850 3000] Units yr
values_at_bounds	values Magnitude [0.0 2.0 2.0] Units watt/meter2

In [14]:

ts_next_left_open = ct.Timeseries.from_arrays(
    x=Q([1750, 1850, 3000], "yr"),
    y=Q([0.0, 0.0, 2.0], "W / m^2"),
    interpolation=ct.InterpolationOption.PiecewiseConstantNextLeftOpen,
    name="next_left_open",
)
ts_next_left_open.plot()

ts_next_left_open = ct.Timeseries.from_arrays( x=Q([1750, 1850, 3000], "yr"), y=Q([0.0, 0.0, 2.0], "W / m^2"), interpolation=ct.InterpolationOption.PiecewiseConstantNextLeftOpen, name="next_left_open", ) ts_next_left_open.plot()

Out[14]:

<Axes: xlabel='yr', ylabel='watt/meter²'>

In [15]:

ts_next_left_open.discrete

ts_next_left_open.discrete

Out[15]:

continuous_timeseries.TimeseriesDiscrete

name	next_left_open
time_axis	bounds Magnitude [1750 1850 3000] Units yr
values_at_bounds	values Magnitude [0.0 0.0 2.0] Units watt/meter2