气象数据分析与可视化的几个小作业

NCEP/DOE Reanalysis II Monthly数据资料，绘制同一高度层的位势高度叠加温度场的图，绘制同一高度层位势高度叠加风场的图，绘制垂直速度的高度-纬度剖面图

# 引用部分
import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False #这两行需要手动设置

# 数据准备
file1 =r’F:\Desktop\python\44\hgt.mon.mean.nc’
file2 =r’F:\Desktop\python\44\uwnd.mon.mean.nc’
file3 =r’F:\Desktop\python\44\vwnd.mon.mean.nc’

ds1=xr.open_dataset(file1)
ds2=xr.open_dataset(file2)
ds3=xr.open_dataset(file3)

time=ds1[‘time’]
level=ds1[‘level’]
lon=ds1[‘lon’]
lat=ds1[‘lat’]
lons,lats=np.meshgrid(lon,lat)
hgt=ds1[‘hgt’].loc[‘1980-12-01’,850,:,:]
u=ds2[‘uwnd’].loc[‘1980-12-01’,850,:,:]
v=ds3[‘vwnd’].loc[‘1980-12-01′,850,:,:]

#绘制图形
fig = plt.figure(figsize=(10,9),dpi=500)
ax=plt.axes([0,0,1,1],projection=ccrs.PlateCarree(central_longitude=0))
proj = ccrs.PlateCarree(central_longitude=0)

#######
ac1=ax.contourf(lon,lat,hgt,cmap=’Spectral_r’,
transform=ccrs.PlateCarree())
ac2=ax.quiver(lons[::2,::2],lats[::2,::2],u[::2,::2],v[::2,::2]
,width=0.01,transform=proj)
ax.quiverkey(ac2,-0.15,0.75,0.5,label=’风速(10m/s)’,color=’r’,fontproperties={‘size’:2.5})
region=[70,140,0,65]
ax.set_extent(region,crs=proj)

#######
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=0.2, color=’k’, alpha=0.8, linestyle=’–‘)
gl.xlabels_top = False ##关闭上侧坐标显示
gl.ylabels_right = False ##关闭右侧坐标显示
gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(region[0], region[1], 1))
gl.ylocator = mticker.FixedLocator(np.arange(region[2], region[3], 1))

#########
gl.xlabel_style={‘size’:20}
gl.ylabel_style={‘size’:20}

#####
ax.coastlines(’50m’,lw=0.5)

ax.set_title(‘Placeholder’,fontsize=30)
ax.set_title(“850hPa”,y=1,loc=’left’)
ax.set_title(“1980年12月1日”,y=1,loc=’right’)

NCEP/DOE Reanalysis II的位势高度、垂直速度、风场、相对湿度等逐日资料绘制此次天气过程中的环流形势图

同上，限制经纬度到达范围内即可

FY-2B TBB数据绘制

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(r’E:\Desktop\python\data\FY2E_TBB_IR1_NOM_20180228_2030.hdf’)

# 查看文件信息
# 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()))即可查看
”’
#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(r’E:\Desktop\python\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=14, 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=16)
plt.savefig(‘test.png’)
plt.show()

好看就完事了