Python使用ccplot绘制CALIPSO L1B后向散射
之前的一篇文章介绍了如何在终端/命令行窗口使用ccplot绘制后向散射图，只需要一行简单的代码即可，但···是···，这样不好二创啊，所以还是回归python吧~
好在ccplot官网提供了代码，下面简单讲解一下，嘻嘻。
首先导入必要的包，成功安装ccplot后这些包都会有的，就算不安装ccplot，理清了思路应该也是可以实现的（应该吧哈哈）！
x1和x2是横坐标的索引，下述代码也就是说要取前1000个点进行绘图。
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from ccplot.hdf import HDF
from ccplot.algorithms import interp2d_12
import ccplot.utils
import ccplot.config
filename = 'CAL_LID_L1-ValStage1-V3-01.2007-06-12T03-42-18ZN.hdf'
name = 'Total_Attenuated_Backscatter_532'
label = 'Total Attenuated Backscatter 532nm (km$^{-1}$ sr$^{-1}$)'
colormap = os.path.join(ccplot.config.sharepath, 'cmap', 'calipso-backscatter.cmap')
x1 = 0
x2 = 1000
h1 = 0 # km
h2 = 20 # km
nz = 500 # Number of pixels in the vertical.
下述代码就是打开文件，获取时间、高度以及后向散射数据，这里可能会报错：TypeError: expected bytes, str found，意思就是说传入的期望类型是bytes，但实际上是str，解决方法就是在每个str后面加上.enconde()就解决了。
获取完数据后根据x和height构建一个均匀网格，并将后向散射数据插值进去，也就是interp2d_12这个方法的功能。
if __name__ == '__main__':
with HDF(filename) as product:
# Import datasets.
time = product['Profile_UTC_Time'][x1:x2, 0]
height = product['metadata']['Lidar_Data_Altitudes']
dataset = product[name][x1:x2]
# Convert time to datetime.
time = np.array([ccplot.utils.calipso_time2dt(t) for t in time])
# Mask missing values.
dataset = np.ma.masked_equal(dataset, -9999)
# Interpolate data on a regular grid.
X = np.arange(x1, x2, dtype=np.float32)
Z, null = np.meshgrid(height, X)
data = interp2d_12(
dataset[::],
X.astype(np.float32),
Z.astype(np.float32),
x1, x2, x2 - x1,
h2, h1, nz,
)
之后的代码就是绘图了，也是比较容易理解的。
# Import colormap.
cmap = ccplot.utils.cmap(colormap)
cm = mpl.colors.ListedColormap(cmap['colors']/255.0)
cm.set_under(cmap['under']/255.0)
cm.set_over(cmap['over']/255.0)
cm.set_bad(cmap['bad']/255.0)
norm = mpl.colors.BoundaryNorm(cmap['bounds'], cm.N)
# Plot figure.
fig = plt.figure(figsize=(12, 6))
TIME_FORMAT = '%e %b %Y %H:%M:%S UTC'
im = plt.imshow(
data.T,
extent=(mpl.dates.date2num(time[0]), mpl.dates.date2num(time[-1]), h1, h2),
cmap=cm,
norm=norm,
aspect='auto',
interpolation='nearest',
)
ax = im.axes
ax.set_title('CALIPSO %s - %s' % (
time[0].strftime(TIME_FORMAT),
time[-1].strftime(TIME_FORMAT)
))
ax.set_xlabel('Time')
ax.set_ylabel('Altitude (km)')
ax.xaxis.set_major_locator(mpl.dates.AutoDateLocator())
ax.xaxis.set_major_formatter(mpl.dates.DateFormatter('%H:%M:%S'))
cbar = plt.colorbar(
extend='both',
use_gridspec=True
)
cbar.set_label(label)
fig.tight_layout()
plt.savefig('calipso-plot.png')
plt.show()以上就是绘制后向散射的代码，还是比较简单的。但是···上述代码以横坐标索引（0，1000）进行绘图的，并不是以经纬度，这样就不是很方便，因此要进行修改。
简单说一下思路，我们在原来的基础上，在获取经纬度数据，这样就可以对时间进行掩膜处理，处理好的时间范围也就是卫星经过研究区的时间范围，这样（0，len（time））的索引既是研究区域的索引！问题解决！
if __name__ == '__main__':
with HDF(filename.encode()) as product:
# Import datasets.
time = product['Profile_UTC_Time'.encode()][:, 0]
lat = product['Latitude'.encode()][:,0]
lon = product['Longitude'.encode()][:,0]
height = product['metadata'.encode()]['Lidar_Data_Altitudes'.encode()]
tab = product[name.encode()][:]
mask = (lon >= 72.5) & (lon <= 97.5) & (lat >= 32.5) & (lat <= 50)
lat = lat[mask]
lon = lon[mask]
time = time[mask]
tab = tab[mask, :]
print(time.shape,tab.shape)
# Mask missing values.
dataset = np.ma.masked_equal(tab, -9999)
# Interpolate data on a regular grid.
X = np.arange(0,len(time),dtype=np.float32)
Z,null = np.meshgrid(height, X)
data = interp2d_12(
dataset[::],
X.astype(np.float32),
Z.astype(np.float32),
0, len(time),len(time),
h2, h1, nz,
)最后附上一张综合了L1B和VFM的图吧！（红、蓝、黑三线是自己画的，用于标明研究区！）
