Python | 计算位涡平流项

写在前面

最近忙着复习、考试…都没怎么空敲代码，还得再准备一周考试。。。等考完试再慢慢更新了，今天先来浅更一个简单但是使用的python code

• 在做动力机制分析时，我们常常需要借助收支方程来诊断不同过程的贡献，其中最常见的一项就包括水平平流项，如下所示，其中var表示某一个变量，V表示水平风场。

  − V ⋅ ∇ v a r -V\cdot\mathbf{\nabla}var −V⋅∇var

  以位涡的水平平流项为例，展开为

  − ( u ∂ p v ∂ x + v ∂ p v ∂ y ) -(u\frac{\partial pv}{\partial x}+v\frac{\partial pv}{\partial y}) −(u∂x∂pv​+v∂y∂pv​)

  其物理解释为：

  • 位涡受背景气流的调控作用。在存在背景气流的情况下，这个位涡信号受到多少平流的贡献

    对于这种诊断量的计算，平流项，散度项，都可以通过metpy来计算。之前介绍过一次，因为metpy为大大简便了计算

    

    但是，如果这样还是不放心应该怎么办呢，这里可以基于numpy的方式自己来手动计算平流项，然后比较两种方法的差异。下面我就从两种计算方式以及检验方式来计算pv的水平平流项

    首先来导入基本的库和数据，我这里用的wrf输出的风场以及位涡pv,同时，为了减少计算量，对于数据进行区域和高度层的截取

    • 经纬度可以直接从nc数据中获取，我下面是之前的代码直接拿过来用了，还是以wrf*out数据来提取的lon&lat

      import numpy as np
      import pandas as pd
      import xarray as xr
      import matplotlib.pyplot as plt
      import cartopy.crs as ccrs
      import glob
      from matplotlib.colors import ListedColormap 
      from wrf import getvar,pvo,interplevel,destagger,latlon_coords
      from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
      from netCDF4 import Dataset
      import matplotlib as mpl   
      from metpy.units import units
      import metpy.constants as constants
      import metpy.calc as mpcalc
      import cmaps
      # units.define('PVU = 1e-6 m^2 s^-1 K kg^-1 = pvu')
      import cartopy.feature as cfeature
      f_w   = r'L:\JAS\wind.nc'
      f_pv  = r'L:\JAS\pv.nc'
      def  cal_dxdy(file):
          ncfile = Dataset(file)
          P = getvar(ncfile, "pressure")
          lats, lons = latlon_coords(P)
          lon = lons[0]
          lon[lon'projection': ccrs.PlateCarree()})
      proj = ccrs.PlateCarree()
      # plot isotherms
      cs = ax.contour(lon, lat, pv[-1]*10**4, 8, 
                  colors='k',
                    linewidths=0.2)
      ax.coastlines('10m',linewidths=0.5,facecolor='w',edgecolor='k')
      box = [100,121,0,20]
      ax.set_extent(box,crs=ccrs.PlateCarree())
      ax.set_xticks(np.arange(box[0],box[1], 10), crs=proj)
      ax.set_yticks(np.arange(box[2], box[3], 5), crs=proj)
      ax.xaxis.set_major_formatter(LongitudeFormatter())
      ax.yaxis.set_major_formatter(LatitudeFormatter())
      # plot tmperature advection and convert units to Kelvin per 3 hours
      cf = ax.contourf(lon, lat,   tadv[10]*10**4, 
                       np.linspace(-4, 4, 21), extend='both',
                       cmap='RdYlGn', alpha=0.75)
      plt.colorbar(cf, pad=0.05, aspect=20,shrink=0.75)
      ax.set_title('PV Advection Calculation')
      plt.show()
      key: value
                         for key, value in {'u': u, 'v': v, 'w': w}.items()
                         if value is not None}
          return_only_horizontal = {key: value
                                    for key, value in {'u': 'df/dx', 'v': 'df/dy'}.items()
                                    if key in wind_vector}
          gradient_vector = ()
          # Calculate horizontal components of gradient, if needed
          if return_only_horizontal:
              gradient_vector = geospatial_gradient(scalar, dx=dx, dy=dy, x_dim=x_dim, y_dim=y_dim,
                                                    parallel_scale=parallel_scale,
                                                    meridional_scale=meridional_scale,
                                                    return_only=return_only_horizontal.values())
          # Calculate vertical component of gradient, if needed
          if 'w' in wind_vector:
              gradient_vector = (*gradient_vector,
                                 first_derivative(scalar, axis=vertical_dim, delta=dz))
          return -sum(
              wind * gradient
              for wind, gradient in zip(wind_vector.values(), gradient_vector)
          )
      'projection': ccrs.PlateCarree()})
      # 绘制三个不同的图
      plot_contour(ax1, tadv[0] * 10**4, 'Metpy', lon, lat, levels=np.linspace(-1, 1, 21))
      plot_contour(ax2, padv[0] * 10**4, 'Numpy', lon, lat, levels=np.linspace(-1, 1, 21))
      plot_contour(ax3, (np.array(tadv[0]) - padv[0]) * 10**4, 'difference',
                   lon, lat, 
                   levels=11)
      # 显示图形
      plt.tight_layout()
      plt.show()