书画网站建设方案策划丰县做淘宝网站

之前的一篇文章介绍了如何在终端/命令行窗口使用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.configfilename = '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.axesax.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的图吧！（红、蓝、黑三线是自己画的，用于标明研究区！）