Python¿ÉÊÓ»¯ | Öг߶ȶÔÁ÷ϵͳ·´ÉäÂʽØÃæ

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib as mpl
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.units import units
from metpy.interpolate import cross_section, log_interpolate_1d

import warnings
warnings.filterwarnings('ignore')

def yuv_rainbow_24(nc):
    path1 = np.linspace(0.8*np.pi, 1.8*np.pi, nc)
    path2 = np.linspace(-0.33*np.pi, 0.33*np.pi, nc)

    y = np.concatenate([np.linspace(0.3, 0.85, nc*2//5),
                        np.linspace(0.9, 0.0, nc - nc*2//5)])
    u = 0.40*np.sin(path1)
    v = 0.55*np.sin(path2) + 0.1

    rgb_from_yuv = np.array([[1, 0, 1.13983],
                             [1, -0.39465, -0.58060],
                             [1, 2.03211, 0]])
    cmap_dict = {'blue': [], 'green': [], 'red': []}
    for i in range(len(y)):
        yuv = np.array([y[i], u[i], v[i]])
        rgb = rgb_from_yuv.dot(yuv)
        red_tuple = (i/(len(y)-1), rgb[0], rgb[0])
        green_tuple = (i/(len(y)-1), rgb[1], rgb[1])
        blue_tuple = (i/(len(y)-1), rgb[2], rgb[2])
        cmap_dict['blue'].append(blue_tuple)
        cmap_dict['red'].append(red_tuple)
        cmap_dict['green'].append(green_tuple)

    return cmap_dict

cmap = colors.LinearSegmentedColormap('homeyer_rainbow', yuv_rainbow_24(15), mpl.rcParams['image.lut'])

data = xr.open_dataset('/home/mw/input/data3672/wrf_for_metpy_demo.nc').metpy.parse_cf()
data

data['dbz'].max(data['dbz'].metpy.vertical.name).plot(cmap=cmap, vmin=0,
                                                      vmax=60, figsize=(14,6.),
                                                      aspect='equal')

<matplotlib.collections.QuadMesh at 0x7f2c15ca67d0>

start, end = ((48.423, -99.771),(43.476, -88.490))

cross = cross_section(data, start, end, 1000)
cross

fig = plt.figure(1, figsize=(14., 6.))
ax = plt.axes()
pres_axis = np.arange(1000, 80, -10) * units.hPa

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    dbz_interp = log_interpolate_1d(pres_axis, cross['pressure'], cross['dbz'].values)

refl_contour = ax.contourf(cross['XLONG'], pres_axis, dbz_interp,
                           levels=np.arange(0, 75, 5), cmap=cmap)
refl_colorbar = fig.colorbar(refl_contour)

ax.set_yscale('symlog')
ax.set_yticklabels(np.arange(1000, 50, -100))
ax.set_ylim(pres_axis.max(), pres_axis.min())
ax.set_yticks(np.arange(1000, 50, -100))

data_crs = data['dbz'].metpy.cartopy_crs
ax_inset = fig.add_axes([0.072, 0.63, 0.25, 0.25], projection=data_crs)


ax_inset.contourf(data['west_east'], data['south_north'], data['dbz'].max('bottom_top'),
                  levels=np.arange(0, 75, 5), cmap=cmap, zorder=5)


endpoints = data_crs.transform_points(ccrs.Geodetic(),
                                      *np.vstack([start, end]).transpose()[::-1])
ax_inset.scatter(endpoints[:, 0], endpoints[:, 1], c='k', zorder=10)
ax_inset.plot(cross['west_east'], cross['south_north'], c='k', zorder=10)

ax_inset.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'], zorder=1)
ax_inset.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='#9f9f68', zorder=2)
ax_inset.add_feature(cfeature.LAKES.with_scale('50m'), facecolor=cfeature.COLORS['water'],
               edgecolor='#9f9f68', zorder=2)

ax_inset.set_title('')
ax.set_title('WRF Simulated Reflectivity - 2015-07-13 05:00Z Forecast\n'
             'Initialized 2015-07-12 12:00Z, Thompson Microphysics, MYJ PBL\n'
             'Inset: Composite Reflectivity')
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Longitude (degrees east)')
refl_colorbar.set_label('Z (dbZ)')

plt.show()