使用EC再分析画形势场

image.png

代码如下：

#读ec数据并画图hgt,wind,rh
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.util import add_cyclic_point
import matplotlib as mpl
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from cartopy.io.shapereader import Reader
import os
import numpy as np 
import pdb
import datetime,time
import pandas as pd
st=datetime.datetime(1900,1,1,0,0)
start=datetime.datetime.now()

from matplotlib.font_manager import FontProperties
from matplotlib.transforms import offset_copy
font = FontProperties(fname='/mnt/d/work/数据/SimHei.ttf', size=16)

SHP = u'/mnt/d/work/数据/china_shp'
nc = Dataset(u'/mnt/d/work_no_backup/20200828/ECDATA/20200827china.nc')
lons = nc.variables['longitude'][:]
lats = nc.variables['latitude'][:]
t=np.array(nc.variables['time'][:])
olon,olat=np.meshgrid(lons,lats)


timelist=[]
for i in range(0,len(t)):
    timelist.append(st+datetime.timedelta(hours=int(t[i])+8))
    print(timelist[i])

    lev=np.array(nc.variables['level'][:])

    for l in [21,25,30,33]:
        print(lev[l])


        hgt=nc.variables['z'][i,l,:,:]*0.01
        rh=nc.variables['r'][i,l,:,:]
        temp=nc.variables['t'][i,l,:,:]
        u=nc.variables['u'][i,l,:,:]
        v=nc.variables['v'][i,l,:,:]

        proj= ccrs.PlateCarree() 
        fig = plt.figure(figsize=(10,8),dpi=300)  
        ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
        ax.add_geometries(Reader(os.path.join(SHP, 'cnhimap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.3)
#        ax.add_geometries(Reader(os.path.join(SHP_world, '10m_coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.3)
        ax.set_xticks([60,65,70,75,80,85,90,95,100,105,110,115,120,125,130,135,140,145])#需要显示的经度，一般可用np.arange
        ax.set_yticks([10,15,20,25,30,35,40,45,50,55,60])#需要显示的纬度
        ax.xaxis.set_major_formatter(LongitudeFormatter())#将横坐标转换为经度格式
        ax.yaxis.set_major_formatter(LatitudeFormatter())#将纵坐标转换为纬度格式
        ax.tick_params(axis='both',labelsize=8,direction='in',length=2.75,width=0.55,right=True,top=True)#修改刻度样式
        ax.grid(linewidth=0.4, color='k', alpha=0.45, linestyle='--')#开启网格线
        extent=[108, 127, 22,40]
        ax.set_extent(extent,ccrs.PlateCarree())
        #高度场
        levels=np.arange(pd.Series(hgt.flatten()).dropna().round(0).min(),pd.Series(hgt.flatten()).dropna().round(0).max(),1)
        hgt_line = ax.contour(olon,olat,hgt,levels=levels,colors=['blue'],linewidths=0.8,transform=ccrs.PlateCarree())
        ax.clabel(
                hgt_line,  # Typically best results when labelling line contours.
                colors=['black'],
                manual=False,  # Automatic placement vs manual placement.
                inline=True,  # Cut the line where the label will be placed.
                fmt=' {:.0f} '.format,  # Labes as integers, with some extra space.
                fontsize=4.5,
                inline_spacing=0,
                )
        #温度场
        levels=np.arange(pd.Series(temp.flatten()).dropna().round(0).min(),pd.Series(temp.flatten()).dropna().round(0).max(),2)
        temp_line = ax.contour(olon,olat,temp,levels=levels,colors=['red'],linestyles='--',linewidths=0.8,transform=ccrs.PlateCarree())
        ax.clabel(
                temp_line,  # Typically best results when labelling line contours.
                colors=['black'],
                manual=False,  # Automatic placement vs manual placement.
                inline=True,  # Cut the line where the label will be placed.
                fmt=' {:.0f} '.format,  # Labes as integers, with some extra space.
                fontsize=5.5,
                inline_spacing=0,
                )
        #相对湿度
        levels=np.arange(70,110,5)
        rh_fill = ax.contourf(olon,olat,rh,levels=levels,cmap='Greens',transform=ccrs.PlateCarree())
        position=fig.add_axes([0.92,0.16,0.02,0.68])
        cb=fig.colorbar(rh_fill,cax=position,shrink=0.3,pad=0.01)

        ax.text(0.09, 1.01, str(timelist[i]),horizontalalignment='center',verticalalignment='center',transform=ax.transAxes)
        ax.text(0.9, 1.01,str(lev[l])+'hPa,H,T,Wind,RH',horizontalalignment='center',verticalalignment='center',transform=ax.transAxes)


        barbs_inspace=2
        ax.barbs(olon[::barbs_inspace,::barbs_inspace],olat[::barbs_inspace,::barbs_inspace],
        u[::barbs_inspace,::barbs_inspace],v[::barbs_inspace,::barbs_inspace],
        linewidth=0.45,
        barb_increments={'half':2,'full':4,'flag':20},
        sizes=dict(spacing=0.1,width=0.05,emptybarb=0.02,height=0.35),length=3.5,zorder=5)
        geodetic_transform = ccrs.Geodetic()._as_mpl_transform(ax)
        text_transform = offset_copy(geodetic_transform, units='dots', x=-25)
        ax.text(121.4283,29.0079,u'台州机场',
                horizontalalignment='center',
                verticalalignment='top',
                transform=text_transform,
                fontproperties=font,
                fontsize=10)
        ax.scatter(121.4283,28.5579,marker='*',c='red',s=40.0,transform=ccrs.PlateCarree())
        plt.savefig('./pic/形势场改/'+str(timelist[i])+'_'+str(lev[l])+'.png')
        plt.close()