In the second part of this two-part post we’ll learn to use different tools and techniques to visualize and save animations (first part here). All the code discussed here is available on a GitHub repository, here.
This part focuses on the moviepy Python library, and all the neat things one can do with it. There actually are some nice tutorials for when we have a continuous function t -> f(t) to work with (see here). Instead, we are often working with data structures that are indexed on time in a discrete way.
Moviepy could be used from any data source dependent on time, including netCDF data such as the one manipulated by VisIt in the first part of this post. But in this second part, we are instead going to focus on how to draw time -dependent trajectories to make sense of nonlinear dynamical systems, then animate them in GIF. I will use the well-known shallow lake problem, and go through a first example with detailed explanation of the code. Then I’ll finish with a second example showing trajectories.
Part I: using state trajectories to understand the concept of stable equilibria
The shallow lake problem is a classic problem in the management of coupled human and natural system. Some human (e.g. agriculture) produce phosphorus that eventually end up in water bodies such as lakes. Too much phosphorus in lake causes a processus called eutrophication which usually destroys lakes’ diverse ecosystems (no more fish) and lower water quality. A major problem with that is that eutrophication is difficult or even sometimes impossible to reverse: lowering phosphorus inputs to what they were pre-eutrophication simply won’t work. Simple nonlinear dynamics, first proposed by Carpenter et al. in 1999 (see here) describe the relationship between phosphorus inputs (L) and concentration (P). The first part of the code (uploaded to GitHub as
import attractors import numpy as np from moviepy.video.io.bindings import mplfig_to_npimage from moviepy.video.VideoClip import DataVideoClip import matplotlib.pyplot as plt import matplotlib.lines as mlines # Lake parameters b = 0.65 q = 4 # One step dynamic (P increment rate) # arguments are current state x and lake parameters b,q and input l def Dynamics(x, b, q, l): dp = (x ** q) / (1 + x ** q) - b * x + l return dp
Where the first 6 lines contain the usual library imports. Note that I am importing an auxiliary Python function “attractors” to enable me to plot the attractors (see
attractors.pyon the GitHub repository). The function “Dynamics” correspond to the evolution of P given L and lake parameters b and q, also given in this bit of code. Then we introduce the time parameters:
# Time parameters dt = 0.01 # time step T = 40 # final horizon nt = int(T/dt+1E-6) # number of time steps
To illustrate that lake phosphorus dynamics depend not only on the phosphorus inputs L but also on initial phosphorus levels, we are going to plot P trajectories for different constant values of L, and three cases regarding the initial P. We first introduce these initial phosphorus levels, and the different input levels, then declare the arrays in which we’ll store the different trajectories
# Initial phosphorus levels pmin = 0 pmed = 1 pmax = 2.5 # Inputs levels l = np.arange(0.001,0.401,0.005) # Store trajectories low_p = np.zeros([len(l),nt+1]) # Correspond to pmin med_p = np.zeros([len(l),nt+1]) # Correspond to pmed high_p = np.zeros([len(l),nt+1]) # Correspond to pmax
Once that is done, we can use the attractor import to plot the equilibria of the lake problem. This is a bit of code that is the GitHub repository associated to this post, but that I am not going to comment on further here.
After that we can generate the trajectories for P with constant L, and store them to the appropriate arrays:
# Generating the data: trajectories def trajectory(b,q,p0,l,dt,T): # Declare outputs time = np.arange(0,T+dt,dt) traj = np.zeros(len(time)) # Initialize traj traj = p0 # Fill traj with values for i in range(1,len(traj)): traj[i] = traj[i-1] + dt * Dynamics(traj[i-1],b,q,l) return traj # Get them! for i in range(len(l)): low_p[i,:] = trajectory(b,q,pmin,l[i],dt,T) med_p[i, :] = trajectory(b, q, pmed, l[i], dt, T) high_p[i,:] = trajectory(b,q,pmax,l[i],dt,T)
Now we are getting to the interesting part of making the plots for the animation. We need to declare a figure that all the frames in our animation will use (we don’t want the axes to wobble around). For that we use matplotlib / pyplot libraries:
# Draw animated figure fig, ax = plt.subplots(1) ax.set_xlabel('Phosphorus inputs L') ax.set_ylabel('Phosphorus concentration P') ax.set_xlim(0,l[-1]) ax.set_ylim(0,pmax) line_low, = ax.plot(l,low_p[:,0],'.', label='State, P(0)=0') line_med, = ax.plot(l,med_p[:,0],'.', label='State, P(0)=1') line_high, = ax.plot(l,high_p[:, 0], '.', label='State, P(0)=2.5')
Once that is done, the last things we need to do before calling the core moviepy functions are to 1) define the parameters that manage time, and 2) have a function that makes frames for the instant that is being called.
For 1), we need to be careful because we are juggling with different notions of time, a) time in the dynamics, b) the index of each instant in the dynamics (i.e., in the data, the arrays where we stored the trajectories), and c) time in the animation. We may also want to have a pause at the beginning or at the end of the GIF, rather than watch with tired eyes as the animation is ruthlessly starting again before we realized what the hell happened. So here is how I declared all of this:
# Parameters of the animation initial_delay = 0.5 # in seconds, delay where image is fixed before the animation final_delay = 0.5 # in seconds, time interval where image is fixed at end of animation time_interval = 0.25 # interval of time between two snapshots in the dynamics (time unit or non-dimensional) fps = 20 # number of frames per second on the GIF # Translation in the data structure data_interval = int(time_interval/dt) # interval between two snapshots in the data structure t_initial = -initial_delay*fps*data_interval t_final = final_delay*fps*data_interval time = np.arange(t_initial,low_p.shape+t_final,data_interval) # time in the data structure
Now for 2), the function to make the frames resets the parts of the plot that change for different time indexes(“t” below is the index in the data). If we don’t do that, the plot will keep the previous plotted elements, and will grow messier at the animation goes on.
# Making frames def make_frame(t): t = int(t) if t<0: return make_frame(0) elif t>nt: return make_frame(nt) else: line_low.set_ydata(low_p[:,t]) line_med.set_ydata(med_p[:,t]) line_high.set_ydata(high_p[:, t]) ax.set_title(' Lake attractors, and dynamics at t=' + str(int(t*dt)), loc='left', x=0.2) if t > 0.25*nt: alpha = (t-0.25*nt) / (1.5*nt) lakeAttBase(eqList, 0.001, alpha=alpha) plt.legend(handles=[stable, unstable], loc=2) return mplfig_to_npimage(fig)
In the above mplfig_to_npimage(fig) is a moviepy function that turns a figure into a frame of our GIF. Now we just have to call the function to do frames using the data, and to turn it into a GIF:
# Animating animation = DataVideoClip(time,make_frame,fps=fps) animation.write_gif("lake_attractors.gif",fps=fps)
Where the moviepy function DataVideoClip takes as arguments the sequences of indexes defined by the vector “time” defined in the parameters of the animation, the “make_frame” routine we defined, and the number of frame per second we want to output. The last lines integrates each frame to the GIF that is plotted below:
Each point on the plot represent a different world (different constant input level, different initial phosphorus concentration), and the animation shows how these states converge towards an stable equilibriun point. The nonlinear lake dynamics make the initial concentration important towards knowing if the final concentration is low (lower set of stable equilibria), or if the lake is in a eutrophic state (upper set of stable equilibria).
Part II: plotting trajectories in the 2-D plane
Many trajectories can be plotted at the same time to understand the behavior of attractors, and visualize system dynamics for fixed human-controlled parameters (here, the phosphorus inputs L). Alternatively, if one changes the policy, trajectories evolve depending on both L and P. This redefines how trajectories are defined.
I did a similar bit of code to show how one could plot trajectories in the 2D plane. It is also uploaded on the GitHub repository (under
movie2.py), and is similar in its structure to the code above. The definition of the trajectories and where to store them change. We define trajectories where inputs are lowered at a constant rate, with a minimum input of 0.08. For three different initial states, that gives us the following animation that illustrates how the system’s nonlinearity leads to very different trajectories even though the starting positions are close and the management policy, identical:
This could easily be extended to trajectories in higher dimensional planes, with and without sets of equilibria to guide our eyes.