Python批量绘制环流形势---era5资料

中者

from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter
import xarray as xr
from netCDF4 import Dataset
import os
import matplotlib.pyplot as plt
import numpy as np
from copy import copy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import shapely.geometry as sgeom
import xarray as xr
import gc

#给画刻度用来辅助，确定四周
def find_side(ls, side):
    """
Given a shapely LineString which is assumed to be rectangular, return the
line corresponding to a given side of the rectangle.

"""
    minx, miny, maxx, maxy = ls.bounds
    points = {'left': [(minx, miny), (minx, maxy)],
              'right': [(maxx, miny), (maxx, maxy)],
              'bottom': [(minx, miny), (maxx, miny)],
              'top': [(minx, maxy), (maxx, maxy)],}
    return sgeom.LineString(points[side])
#用来画坐标轴的刻度（包括经度和纬度）
def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
    """得到一个兰伯特正形投影的轴的刻度位置和标签."""
    outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
    axis = find_side(outline_patch, tick_location)
    n_steps = 30
    extent = ax.get_extent(ccrs.PlateCarree())
    _ticks = []
    for t in ticks:
        xy = line_constructor(t, n_steps, extent)
        proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
        xyt = proj_xyz[..., :2]
        ls = sgeom.LineString(xyt.tolist())
        locs = axis.intersection(ls)
        if not locs:
            tick = [None]
        else:
            tick = tick_extractor(locs.xy)
        _ticks.append(tick[0])
    # Remove ticks that aren't visible:
    ticklabels = copy(ticks)
    while True:
        try:
            index = _ticks.index(None)
        except ValueError:
            break
        _ticks.pop(index)
        ticklabels.pop(index)
    return _ticks, ticklabels
#设置x轴标签（用来画纬度）
def lambert_xticks(ax, ticks):
    """Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
    te = lambda xy: xy[0]
    lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
    xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
    ax.xaxis.tick_bottom()
    ax.set_xticks(xticks)
    ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])
#设置y轴标签（用来画经度）
def lambert_yticks(ax, ticks):
    """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    te = lambda xy: xy[1]
    lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
    yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
    ax.yaxis.tick_left()
    ax.set_yticks(yticks)
    ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])


#################################以上为绘制地图########################
def huanliu(file,levs,windmin,path):
for file1 in file:
    plt.rcParams['font.sans-serif']=['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    print(file1)
    ds = xr.open_dataset(file1)
    lats = ds.latitude.data
    lons = ds.longitude.data

# 选取并抓取数据
    hght = ds['z']
    uwnd = ds['u']
    vwnd = ds['v']
    hght_500 = gaussian_filter(hght.sel(level=levs).data[0], sigma=3.0)
    uwnd_500 = gaussian_filter(uwnd.sel(level=levs).data[0], sigma=3.0) * units('m/s')
    vwnd_500 = gaussian_filter(vwnd.sel(level=levs).data[0], sigma=3.0) * units('m/s')
# 使用MetPy计算填色图的风速，将单位从m / s更改为节
    sped_500 = mpcalc.wind_speed(uwnd_500, vwnd_500).to('kt')
# 创建一个干净的datetime对象以根据地势高度的时间进行绘图
    vtime = datetime.strptime(str(ds.time.data[0].astype('datetime64[ms]')),'%Y-%m-%dT%H:%M:%S.%f')
# 设置将用于绘图的投影
    mapcrs = ccrs.LambertConformal(central_longitude=70,central_latitude=40)

# 设置数据投影； 如果是经/纬度，那么PlateCarree是您想要的
    datacrs = ccrs.PlateCarree()

# 启动图形并使用正确的投影创建绘图轴
    fig = plt.figure(1, figsize=(12, 10),dpi=150)
    ax = plt.subplot(111, projection=mapcrs)
#ax = fig.add_axes([0.08, 0.05, 0.8, 0.94], projection=ccrs.LambertConformal(central_latitude=45, central_longitude=80))
    ax.set_extent([40, 100, 20, 55],crs=ccrs.PlateCarree())
    fig.canvas.draw()

#设置刻度值，画经纬网格
    xticks = [40,45,50,55,60, 65,70, 75,80, 85,90, 95, 100]
    yticks = [15,20,25,30 , 35 , 40 , 45 , 50]
    ax.gridlines(xlocs=xticks, ylocs=yticks,linestyle='--',lw=1,color='dimgrey')
    ax.tick_params(labelsize=12)
# 设置经纬网格的端点（也是用来配合画刻度的，注释掉以后刻度不能正常显示了）
    ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
    ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
#画经纬度刻度
    lambert_xticks(ax, xticks)
    lambert_yticks(ax, yticks)
    with open(r'D:\MeteoInfo_3.2.2\FY4\CN-border-La.dat') as src:
      context = src.read()
      blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
      borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
    for line in borders:
      ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k',
            transform=ccrs.Geodetic())
# 添加地缘边界以供地图参考
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
    ax.add_feature(cfeature.STATES.with_scale('50m'))

# 以节为单位绘制500-hPa填色
    clevs_500_sped = np.arange(windmin,windmin+30,4 )
    cbar_ax = fig.add_axes([0.9, 0.2, 0.02, 0.6])
    cf = ax.contourf(lons, lats, sped_500/2, clevs_500_sped, cmap=plt.cm.BuPu,transform=datacrs)
#cb=plt.colorbar(cf, orientation='vertical', pad=0.0, aspect=30)
#cb.set_label('Abs. Vorticity ($s^{-1}$)')
    plt.colorbar(cf,label="速度·m/s", cax=cbar_ax,ticks=clevs_500_sped)
# 以米为单位绘制500-hPa地势高度
    clevs_500_hght = np.arange(500, 800, 4)
    cs = ax.contour(lons, lats, hght_500/100., clevs_500_hght,linewidths=1.5, colors='black',transform=datacrs)
    plt.clabel(cs, fmt='%d',inline=1,fontsize = 14)

# 以节为单位绘制500-hPa的风羽图，以减少风羽的数量
    ax.barbs(lons, lats, uwnd_500.to('kt').m, vwnd_500.to('kt').m, pivot='middle',color='black', regrid_shape=20, transform=datacrs)

# 为该图制作一些漂亮的标题（一右一左）
#plt.title('500-hPa NAM Geopotential Heights (m), Wind Speed (kt), and Wind Barbs (kt)', loc='left')
#plt.title('Valid Time: {}'.format(vtime), loc='right')

# 调整图像并展示
    plt.subplots_adjust(bottom=0, top=1)
    timss=datetime.strftime(vtime,'%Y-%m-%d%H')
    fig.savefig("{}\\{}.png".format(path,timss))#快速预览图
    fig.savefig("{}\\{}.svg".format(path,timss))#高清矢量图
    plt.clf()
    ax.cla()
    del fig,ax
def get_files_path_list(path):
    files = []
    filesList = os.listdir(path)
    for filename in filesList:
     if os.path.splitext(filename)[1] == '.nc':
        fileAbsPath = os.path.join(path,filename)
        files.append(fileAbsPath)
    return files
##################画图设置###########################
path = r'D:\MeteoInfo_3.2.2\era5\ERA5_netCDF'#era5资料存放路径
pathsave=r'D:\MeteoInfo_3.2.2\era5\tu'#图片保存地址
#获取文件的绝对路径列表
files_path = get_files_path_list(path)
huanliu(files_path,500,20,pathsave)#files_path文件名，500代表层次，20代表急流起始，pathsave图像保存路径

cookie-o-o

ERA高度单位:units = "m**2 s**-2";需要除以重力加速度，转换为米，这样画出来就和所谓实况一致了

_丢丢_

加油！

love2100

学习一下，谢谢楼主

gaoyuan215

太复杂了。。。

中者

gaoyuan215 发表于 2021-12-2 11:35
太复杂了。。。

前边属于给兰波托投影添加坐标，其实可以去掉，实际代码不多

fantasy_v

老师好，我想请问您是否碰到过高度场数值和实况差别较大的情况，如实况588，绘图出来为572

中者

cookie-o-o 发表于 2021-12-3 07:53
ERA高度单位:units = "m**2 s**-2";需要除以重力加速度，转换为米，这样画出来就和所谓实况一致了

里面有除转换成位势十米

中者

fantasy_v 发表于 2021-12-2 18:58
老师好，我想请问您是否碰到过高度场数值和实况差别较大的情况，如实况588，绘图出来为572

重力加速度的问题，除以98试试，层次越高表现的差异越大

fantasy_v

中者 发表于 2021-12-3 11:04
重力加速度的问题，除以98试试，层次越高表现的差异越大

明白了，单位需要转为位势米，感谢您