- 积分
- 976
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-8-23
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
1、数据为国家卫星中心下载的HDF5是全球标称数据,2288*2288的格点2、每个格点对应一个一个经纬度。对应的经纬度数据文件可在卫星中心下载
3、读取hdf5数据使用h5py库
4、cmap为使用ncl的色标,论坛里有,自行查找
- import h5py
- import cmaps
- import numpy as np
- import matplotlib.pyplot as plt
- from mpl_toolkits.basemap import Basemap
- # read hdf5 file
- f = h5py.File('data/FY2E_TBB_IR1_NOM_20170912_1030.hdf','r')
- # 查看文件信息
- # print(list(f.keys()))
- data = f.get('FY2E TBB Hourly Product')
- # print(list(data.attrs.items()))
- data = np.array(data)
- '''
- 处理缺值,数据里有效数据的最大值为340,最小为160
- print(list(data.attrs.items()))即可查看
- '''
- data = np.ma.masked_outside(data,160,340)
- info = f.get('NomFileInfo')[0]
- '''
- 卫星的中心经度 info[4]
- 中心纬度为0
- '''
- lonCenter = info[4]
- # read lonlat msg
- loc = np.fromfile('data/NOM_ITG_2288_2288(0E0N)_LE.dat',dtype=np.float32,count=-1)
- latlon = loc.reshape(2,2288,2288)
- lon = np.array(latlon[0]+lonCenter)
- lat = np.array(latlon[1])
- '''
- 经纬度信息里300为缺值
- '''
- lon = np.ma.masked_where(lon>300,lon)
- lat = np.ma.masked_values(lat,300)
- # 画图
- fig = plt.figure(figsize=(18.5, 10.5))
- m=Basemap(projection='cyl',llcrnrlat=10,llcrnrlon=90,urcrnrlat=35,urcrnrlon=140)
- m.drawcoastlines(color='blue')
- m.drawstates(color='blue')
- m.drawcountries(color='blue')
- x,y = m(lon, lat) # 将lats / lons转换为地图投影坐标
- cf = m.contourf(x,y,data,levels=np.linspace(np.min(data),np.max(data),400),cmap=cmaps.MPL_RdGy )
- cbar = m.colorbar(cf ,location='right',size='3%' ,pad='2%')
- cbar.set_label('Brightness Temperature ( K )', fontsize=12, fontweight='bold')
- m.drawmeridians(np.arange(90, 140, 5),labels=[0, 0, 0, 1])
- m.drawparallels(np.arange(10, 35, 5),labels=[1, 0, 0, 0])
- plt.title('Brightness Temperature(20170912_1030) ',fontweight='bold', fontsize=14)
- plt.savefig('test.png')
- plt.show()
复制代码
|
评分
-
查看全部评分
|