请问python可以实现水汽通量散度计算吗

slddbs

卡在这很久，找到资料不多，利用metpy库的函数divergence不知为何dx，dy步长一直报错

edwardli

##引用部分
import numpy as np                           # ┐  数  基础
from datetime import datetime                # ┝  学  时间处理
import xarray as xr                          # |  工  数据结构
from scipy.ndimage import gaussian_filter    # ┘  具  平滑工具
import matplotlib.pyplot as plt            # ┐
import cartopy.crs as ccrs                 # ┝ 绘图工具
import cartopy.feature as cfeature         # ┘
import metpy.calc as mpcalc          # ┐  MetPy
from metpy.units import units        # ┘专业工具

##数据准备
ds=xr.open_dataset('data/Metpy/fnl_20160719_00_00.grib2',engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa'}})

# 850hPa数据选择
r = ds['r']
t = ds['t']

# h850 = hght.sel(isobaricInhPa=850).data
u850 = uwnd.sel(isobaricInhPa=850).data * units('m/s')
v850 = vwnd.sel(isobaricInhPa=850).data * units('m/s')
r850 = r.sel(isobaricInhPa=850).data * units.percent
t850 = t.sel(isobaricInhPa=850).data * units.kelvin
dew  = mpcalc.dewpoint_from_relative_humidity(t850,r850)
q = mpcalc.specific_humidity_from_dewpoint(850*units.mbar, dew)

lats = ds.latitude.data
lons = ds.longitude.data
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)  

g = 9.8* units.meter / (units.second * units.second)
qdiv = mpcalc.divergence(q*u850/g ,q*v850/g ,dx=dx,dy=dy, x_dim=-1, y_dim=-2)

就是他了。提示：metpy版本1.0.1

edwardli

这个很简单，找一下公式，写几行代码就行。
还可以去找GrADS、C、FORTRAN语言的脚本、程序，改成Python即可。

还有一个办法是理解水汽通量散度和散度的关系：
metpy.calc.divergence(u, v, *, dx=None, dy=None, x_dim=- 1, y_dim=- 2)是散度
[size=13.125px]将其中的u替换为  q*u/g  v替换为   q*v/g，计算出的应该就是水汽通量散度。


============================
关于dx，dy正常的写法应该是
lats = ds.latitude.data
lons = ds.longitude.data
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)  

edwardli

slddbs 发表于 2021-1-14 15:43
divergence
metpy.calc.divergence(u, v, *, dx=None, dy=None, x_dim=- 1, y_dim=- 2)[source]
Calcul ...

首先，pint.Quantity不是数据格式，而是一种数据结构，就和pandas.dataframe、ndarray、xarray.dataset、netcdf4.dataset一样，是某个库/包/模块存储数据的结构体。只是每种结构有各自优势和特点。
其次，就如你提到的如果是xarray.dataset就不需要dx、dy，就是因为这个结构体不光存储了经纬度格点上的数据，也存储了经纬度网格的描述信息，就可以被metpy这些默认获取（范围、分辨率）并计算。这就是xarray.dataset这种数据结构的优势。如果用普通二维数组ndarray存风场，那么计算散度的时候，也需要自己手动“凑”一个lons、lats、dx、dy
最后，看你的回复意思是用xarray并且不带dxdy运行成功了，那就可以了，写程序以实现目标为主，而不是以用很多路很多方法实现。或者将示例数据放上来，大家可以进行具体操作，更为直接有效率。（在我的计算温度平流、散度的脚本中，使用dx，dy没有报过错……）

slddbs

edwardli 发表于 2021-1-15 10:45
首先，pint.Quantity不是数据格式，而是一种数据结构，就和pandas.dataframe、ndarray、xarray.dataset、 ...

恩，恩，感谢，感谢，

slddbs

edwardli 发表于 2021-1-13 14:31
这个很简单，找一下公式，写几行代码就行。
还可以去找GrADS、C、FORTRAN语言的脚本、程序，改成Python即 ...

感谢大佬的指点
运行过程中出现raise TypeError(
TypeError: too many positional arguments
我把GFS_20101026_1200.nc数据下载下来，然后利用https://www.kesci.com/mw/project/5f278156d278b1002c2486b3中计算例子，同样出现以上错误。
但仅用metpy.calc.divergence(u, v），其中u,v为xarray.DataArray时则可以运行通过
现在对于lons 和lats各种数据格式均运行，仍然存在问题喜欢帮忙分析分析啊

slddbs

edwardli 发表于 2021-1-13 14:31
这个很简单，找一下公式，写几行代码就行。
还可以去找GrADS、C、FORTRAN语言的脚本、程序，改成Python即 ...

divergence
metpy.calc.divergence(u, v, *, dx=None, dy=None, x_dim=- 1, y_dim=- 2)[source]
Calculate the horizontal divergence of a vector.

Parameters
u ((…, M, N) xarray.DataArray or pint.Quantity) – x component of the vector

v ((…, M, N) xarray.DataArray or pint.Quantity) – y component of the vector

dx (pint.Quantity, optional) – The grid spacing(s) in the x-direction. If an array, there should be one item less than the size of u along the applicable axis. Optional if xarray.DataArray with latitude/longitude coordinates used as input. Keyword-only argument.

dy (pint.Quantity, optional) – The grid spacing(s) in the y-direction. If an array, there should be one item less than the size of u along the applicable axis. Optional if xarray.DataArray with latitude/longitude coordinates used as input. Keyword-only argument
观察其使用方法，好像当为xarray.DataArray格式时，dx,dy是不需要参数
那么，请问pint.Quantity是个什么数据格式

yqcxm

不错，学习了先，

花开半夏

用的ERA5资料，也出现TypeError: too many positional arguments这个问题，博主现在解决了吗

xiaoyun11070

lats = ds.latitude.data
lons = ds.longitude.data
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)  

没有理解该怎么写，用的nc数据

xiaoyun11070

最终问题解决了么