python如何绘制登陆时的卫星云图(TBB)
作者:小朱小朱绝不服输
这篇文章主要介绍了python如何绘制登陆时的卫星云图(TBB),具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教
根据前面的文章python处理卫星云图获取亮温值可以从hdf文件中获取某时刻某经纬度的亮温值。
然后就是把这些读取的亮温值可视化,绘制成TBB的卫星云图。
先展示一下结果:
1.读取hdf文件
详见上篇文章,这里不再赘述,再贴一下代码。
主要是从全圆盘标对称文件经纬度对照表获取经纬度和从hdf文件中获取亮温值。
# 从全圆盘标对称文件经纬度对照表获取经纬度 def getLatLonFromdat(): lonlatfile = 'F:/Satellite_Imagery/Code/NOM_ITG_2288_2288(0E0N)_LE.dat' with open(lonlatfile, 'rb') as f: #lon_fy = np.fromfile(f, count=2288 * 2288, dtype='float32') + 79 # 先存经度,根据卫星的不同加上对应的经度值 #lat_fy = np.fromfile(f, count=2288 * 2288, dtype='float32') # 再存纬度 data = np.fromfile(f, dtype='float32') data = data.reshape([2288, 2288, 2], order='F') #lon = lon_fy.reshape([2288, 2288], order='F') #lat = lat_fy.reshape([2288, 2288], order='F') lon = data[:, :, 0] + 104.5 lat = data[:, :, 1] return lon, lat
# 从hdf文件中获取亮温值 def getTBBFromhdf(): hdfFile = h5py.File('F:/IR/Satellite_Imagery/IR_data/利奇马/FY2G_FDI_ALL_NOM_20190811_1200.hdf', 'r') db1 = hdfFile['/CALChannelIR1'] hw1 = hdfFile['/NOMChannelIR1'] # db2 = hdfFile['/CALChannelIR2'] # hw2 = hdfFile['/NOMChannelIR2'] # db3 = hdfFile['/CALChannelIR3'] # hw3 = hdfFile['/NOMChannelIR3'] # db4 = hdfFile['/CALChannelIR4'] # hw4 = hdfFile['/NOMChannelIR4'] infoh = hdfFile['/NomFileInfo'] # 查看卫星的经纬度 lat_hdf = infoh[0][3] lon_hdf = infoh[0][4] # print(lat_hdf) # print(lon_hdf) hw = hw1[()] db = db1[()] # 获取定标表的值 tb = np.zeros(shape=(2288, 2288)) # 2288*2288的图像每个具体的亮温值 for i in range(2288): for j in range(2288): if hw[i][j] == 65535 or hw[i][j] == 65534: tb[i][j] = 0 else: a = hw[i][j] tb[i][j] = db[0][a] tb = tb.T return tb
2. 画图
fig = plt.figure(figsize=(8, 6)) m = Basemap(projection='cyl', llcrnrlat=10, llcrnrlon=110, urcrnrlat=40, urcrnrlon=140) # 使用Basemap绘制地图,这里可以读取对应的地图shp文件。 m.drawcoastlines(color='black') m.drawstates(color='black') m.drawcountries(color='black') x, y = m(lon, lat) # 将lats / lons转换为地图投影坐标 # 绘制轮廓图 # 这里data就是计算的亮温值,x,y就是经纬度投影的坐标 cf = m.contourf(x, y, data, levels=np.linspace(180, 301, 400), cmap='jet') cbar = m.colorbar(cf, location='right', size='5%', pad='2%') font = {'family': 'serif', 'color': 'darkred', 'weight': 'normal', 'size': 16, } cbar.set_label('Brightness Temperature ( K )', fontdict=font) m.drawmeridians(np.arange(110, 140, 5), labels=[0, 0, 0, 1]) m.drawparallels(np.arange(10, 40, 5), labels=[1, 0, 0, 0]) # 将最佳路径集上的经纬度映射到地图上,再把该点绘制在地图上 best_lon, best_lat = m(best_lon, best_lat) m.plot(best_lon, best_lat, 'o', color='fuchsia', ms=5) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.title('Brightness Temperature(2019081112) ', fontdict=font) plt.savefig('test.png') plt.show()
注:读取地图文件时,可以使用readshapefile函数。
以2019081112时刻为例,对应地图10-40,110-140的区域,读取对应的hdf文件,绘制该点的卫星云图的亮温值。
总结
这是找了好久,调整过的结果。
以上仅为个人经验,希望能给大家一个参考,也希望大家多多支持脚本之家。