如何用Python制作学术动图？（数据+代码）

0.前言

2019年年初，某厂年会的一个视频火爆全网，里面说出了无数职场人士的心声：

干活的累死累活,到头来干不过写PPT的！

也有网友表示：写好PPT和做好PPT在职场上就是一种能力，一份好的PPT是内容好加视觉美观。

在平时的科研过程中，我们经常会输出一些二维的平面图，二维平面图反映某个变量在二维场景下的分布情况。

在我们常规的PPT展示中，二维平面图经常会出现在内容中，但是单单一张二维图，往往会显得比较单调。要写好PPT，首先还是要有好的内容，如何有效并酷炫的展示我们的平面二维结果？先尝试让它动起来吧。

1.准备工作

人生苦短，我用Python。

之前的推送中有提到过Matlab除了不会生孩子，什么都会！把Matlab替换成Python这句话同样成立。

1.1 用到的数据

水下滑翔机（Underwater Glider）通过浮力调节实现在海水中的升降，并借助固定翼的水动力实现水下滑翔运动。预设轨迹后，其可自动通过携带CTD等传感器测量航线上数据，使用卫星通讯返回观测数据、并接受遥控指令，是当前海洋水体自主观测的新型装备，也是实现深海大范围、长时序海洋观测与探测的有效技术手段。

我们先下载Glider观测到的一个剖面的温度数据，同时下载观测位置处的水深情况为后续数据处理做准备。分别打开以下链接后会下载两个文件，其中csv文件是温度数据，nc文件是水深数据。

1.2 用到的Python库

from matplotlib.transforms import Bbox, TransformedBbox,  blended_transform_factory
from mpl_toolkits.axes_grid1.inset_locator import BboxPatch, BboxConnector, BboxConnectorPatch
import matplotlib.pyplot as plt
from math import radians, cos, sin, asin, sqrt
import pandas as pd
import numpy as np
import netCDF4 as nc
import imageio

2. 动起手来

2.1 数据读取

nc文件使用netCDF4库来读取文件中的变量信息，csv文件中的信息使用pandas来读取。

topo为所在区域的水深情况，lon, lat, date, depth, temp分别为Glider观测到的站点经度、纬度、观测时间、深度和温度。

file2read = r'otn200_sci_water_temp_live.csv'
ncfile = r'usgsCeSrtm30v1_6ebb_eec1_d277.nc'
data = pd.read_csv(file2read)
date = data.iloc[:, 0]
lat = data.iloc[:, 1]
lon = data.iloc[:, 2]
depth = data.iloc[:, 3]
temp = data.iloc[:, 4]
etopo2 = nc.Dataset(ncfile)
latDepth = etopo2.variables['latitude'][:]
lonDepth = etopo2.variables['longitude'][:]
topo = etopo2.variables['topo'][:]

2.2计算各观测点之间的距离

因为我们要画Glider观测到的断面的数据情况，对于这个平面图来说，x方向为断面上的点距离起点的距离，y方向是水深，平面上展示的为变量的值。先得到断面上各点之间的距离。

def haversine(lon1, lat1, lon2, lat2):  # 经度1，纬度1，经度2，纬度2 （十进制度数）    
"""    using haversin to Calculate the great circle distance between two points on the earth (specified in decimal degrees)    """ 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])    
    # haversine    
    dlon = lon2 - lon1    
    dlat = lat2 - lat1    
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2    
    c = 2 * asin(sqrt(a))    
    r = 6371    
    return c * r
# Compute distance along transect
dist = np.zeros((np.size(lon)))
for i in range(1, np.size(lon)):
    dist[i] = dist[i - 1] + haversine(lon[i - 1], lat[i - 1], lon[i], lat[i])

2.3得到断面上各个点的水深情况

bathy = np.zeros((np.size(lon)))
for i in range(np.size(lon)):
    cost_func = np.abs(lonDepth - lon[i])
    xmin = np.where(cost_func == np.min(cost_func))[0]
    cost_func = np.abs(latDepth - lat[i])
    ymin = np.where(cost_func == np.min(cost_func))[0]
    bathy[i] = -topo[ymin, xmin]

2.4 局部放大图像

局部放大图像并将两个图像连接起来。主要使用到以下两个函数：

def zoom_effect(ax1, ax2, xmin, xmax, **kwargs):
    """
    ax1 : the main axes
    ax1 : the zoomed axes
    (xmin,xmax) : the limits of the colored area in both plot axes.
    connect ax1 & ax2. The x-range of (xmin, xmax) in both axes will
    be marked.  The keywords parameters will be used ti create
    patches.
    Source: http://matplotlib.org/dev/users/annotations_guide.html#zoom-effect-between-axes
    """
    trans1 = blended_transform_factory(ax1.transData, ax1.transAxes)
    trans2 = blended_transform_factory(ax2.transData, ax2.transAxes)
    bbox = Bbox.from_extents(xmin, 0, xmax, 1)
    mybbox1 = TransformedBbox(bbox, trans1)
    mybbox2 = TransformedBbox(bbox, trans2)
    prop_patches = kwargs.copy()
    prop_patches["ec"] = "r"
    prop_patches["alpha"] = None
    prop_patches["facecolor"] = 'none'
    prop_patches["linewidth"] = 2
    c1, c2, bbox_patch1, bbox_patch2, p = \
        connect_bbox(mybbox1, mybbox2,
                     loc1a=3, loc2a=2, loc1b=4, loc2b=1,
                     prop_lines=kwargs, prop_patches=prop_patches)
    ax1.add_patch(bbox_patch1)
    ax2.add_patch(bbox_patch2)
    ax2.add_patch(c1)
    ax2.add_patch(c2)
    ax2.add_patch(p)
    return c1, c2, bbox_patch1, bbox_patch2, p

def connect_bbox(bbox1, bbox2,
                 loc1a, loc2a, loc1b, loc2b,
                 prop_lines, prop_patches=None):
    # Fuctions for zoom effect ****************************************************
    # Source: http://matplotlib.org/dev/users/annotations_guide.html#zoom-effect-between-axes
    if prop_patches is None:
        prop_patches = prop_lines.copy()
        prop_patches["alpha"] = prop_patches.get("alpha", 1) * 0.2
    c1 = BboxConnector(bbox1, bbox2, loc1=loc1a, loc2=loc2a, **prop_lines)
    c1.set_clip_on(False)
    c2 = BboxConnector(bbox1, bbox2, loc1=loc1b, loc2=loc2b, **prop_lines)
    c2.set_clip_on(False)
    bbox_patch1 = BboxPatch(bbox1, **prop_patches)
    bbox_patch2 = BboxPatch(bbox2, **prop_patches)
    p = BboxConnectorPatch(bbox1, bbox2,
                           # loc1a=3, loc2a=2, loc1b=4, loc2b=1,
                           loc1a=loc1a, loc2a=loc2a, loc1b=loc1b, loc2b=loc2b,
                           **prop_patches)
    p.set_clip_on(False)
    return c1, c2, bbox_patch1, bbox_patch2, p

2.5 生成多张图片

要生成动图，首先要准备多张图片来成为动图每一帧的素材，在下面这段代码中，nframes表示生成多少幅图，这个相当于将平面图在x方向分割成nframes个部分，在平面图的下方新建一个图层，然后将分割出来的部分再展示出来，实现局部放大效果。

# Make plots
nframes = 20
overlap = 0.95
window = np.floor(max(dist) - min(dist)) / (nframes - (nframes * overlap) + overlap)
xmin = 0
xmax = xmin + window
cmp = plt.cm.get_cmap('jet', 16)
fig1 = plt.figure()
for i in range(0, nframes):
    ax1 = plt.subplot(211)
    plt.fill_between(dist, bathy, 1000, color='k')
    plt.scatter(dist, np.asarray(depth), s=15, c=temp,
                cmap = cmp, marker='o', edgecolor='none')
    plt.ylim((-0.5, max(depth) + 5))
    ax1.set_ylim(ax1.get_ylim()[::-1])
    cbar = plt.colorbar(orientation='vertical', extend='both')
    cbar.ax.set_ylabel('Temperature ($^\circ$C)')
    plt.title('OTN Glider transect')
    plt.ylabel('Depth (m)')
    ax1.set_xlim(min(dist), max(dist))
    ax2 = plt.subplot(212)
    plt.fill_between(dist, bathy, 1000, color='k')
    plt.scatter(dist, depth, s=15, c=temp,
                cmap= cmp, marker='o', edgecolor='none')
    plt.ylim((-0.5, max(depth) + 5))
    ax2.set_ylim(ax2.get_ylim()[::-1])
    plt.ylabel('Depth (m)')
    plt.xlabel('Distance along transect (km)')
    ax2.set_xlim(xmin, xmax)
    zoom_effect(ax1, ax2, xmin, xmax)
    file2write = 'glider_' + '%02d'%i + '.png'
    plt.savefig(file2write, bbox_inches='tight')
    plt.close()
    xmin = xmax - np.floor(window * overlap)
    xmax = xmin + window

2.6 生成gif文件

这里使用到一个imageio的库，首先读取上文中生成的静态图到列表，作为GIF的每一帧；接着设置输入（静态图)、输出(动态图）和一些必要参数，我们这里设置每一帧间隔时间为0.3秒，默认是1秒，可以通过调整dutation来设置逐帧之间的间隔时间。之后调用imageio.mimsave函数去保存结果。

imagelist = ['glider_'+'%02d'%i +'.png' for i in range(0,20)]
frames = [imageio.imread(_) for _ in imagelist]
outname = 'gliders.gif'
imageio.mimsave(outname, frames, 'GIF', duration=0.3)

3.小结

一图胜千言，就用本篇生成的动图来结束这篇推送吧。

和前文的动图似乎有点不一样？可以从公众号推送中找找如何实现。

最后，感谢阅读。

关注作者公众号 海洋纪 ，后台回复 Python动图 获取数据及本文代码。