计算卫星数据中 u 和 v 风分量的卷曲 - Python答案

【问题标题】：calculating the curl of u and v wind components in satellite data - Python计算卫星数据中 u 和 v 风分量的卷曲 - Python
【发布时间】：2015-12-18 13:14:51
【问题描述】：

我不确定如何在卫星数据中获取风的 u 和 v 分量的导数。我以为我可以这样使用 numpy.gradient ：

    from netCDF4 import Dataset      
    import numpy as np      
    import matplotlib.pyplot as plt 

    GridSat = Dataset('analysis_20040713_v11l30flk.nc4','r',format='NETCDF4')
    missing_data = -9999.0
    lat = GridSat.variables['lat']   
    lat = lat[:]     
    lat[np.where(lat==missing_data)] = np.nan  
    lat[np.where(lat > 90.0)] = np.nan     

    lon = GridSat.variables['lon']   
    lon = lon[:]                
    lon[np.where(lon==missing_data)] = np.nan


    uwind_data = GridSat.variables['uwnd']  
    uwind = GridSat.variables['uwnd'][:]
    uwind_sf = uwind_data.scale_factor   
    uwind_ao = uwind_data.add_offset
    miss_uwind = uwind_data.missing_value

    uwind[np.where(uwind==miss_uwind)] = np.nan    


    vwind_data = GridSat.variables['vwnd']  
    vwind = GridSat.variables['vwnd'][:]
    vwind_sf = vwind_data.scale_factor    
    vwind_ao = vwind_data.add_offset
    miss_vwind = vwind_data.missing_value

    vwind[np.where(vwind==miss_vwind)] = np.nan  


    uwind = uwind[2,:,:]
    vwind = vwind[2,:,:]  

    dx = 28400.0 # meters calculated from the 0.25 degree spatial gridding 
    dy = 28400.0 # meters calculated from the 0.25 degree spatial gridding 

    dv_dx, dv_dy = np.gradient(vwind, [dx,dy])
    du_dx, du_dy = np.gradient(uwind, [dx,dy])


    File "<ipython-input-229-c6a5d5b09224>", line 1, in <module>
     np.gradient(vwind, [dx,dy])

    File "/Users/anaconda/lib/python2.7/site-packages/nump/lib/function_base.py", line 1040, in gradient
out /= dx[axis]

    ValueError: operands could not be broadcast together with shapes (628,1440) (2,) (628,1440)

老实说，我不确定如何计算具有 (0.25x0.25) 度间距的卫星数据的中心差异。我也不认为我的 dx 和 dy 是正确的。如果有人对在卫星数据中进行这些类型的计算有一个好主意，我将不胜感激。谢谢！！

【问题讨论】：

这是一个有趣的问题。计算一组离散点的旋度与计算连续场有些不同。有几种方法可以解决：1）如果您想要“更平滑”的数据输出或有噪声输入信息，则取向量在 (i,j) 和八个（或更多）处的叉积的平均值) 周围的网格方块按与感兴趣的网格方块的距离加权。 2) 最简单的方法是仅在 4 个相邻的方格上使用第一种方法，但对角线特征的分辨率会较低。
3) 有一些数值方法可用于获取离散数据并将其映射为 n 阶连续函数，这将是“最佳”方法，但也是迄今为止计算量最大的方法实施起来很复杂。
numpy.gradient 的文档在这方面有点狡猾，但正确的称呼是：gradient(vwind, dx, dy)。 IE。函数签名是gradient(f, *varargs, **kwargs)，其中varargs 是一个由“splat”或“unpack”运算符扩展的列表。
我删除了 curl 标签，因为它指的是 Linux 工具，而不是数学概念。

标签： python numpy vector signal-processing weather

【解决方案1】：

下面的代码可以在matlab风数据集上运行，文件wind.mat在

http://bioinformatics.intec.ugent.be/MotifSuite/INCLUSive_for_users/CPU_64/Matlab_Compiler_Runtime/v79/toolbox/matlab/demos/

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import scipy.io as sio

def curl(x,y,z,u,v,w):
    dx = x[0,:,0]
    dy = y[:,0,0]
    dz = z[0,0,:]

    dummy, dFx_dy, dFx_dz = np.gradient (u, dx, dy, dz, axis=[1,0,2])
    dFy_dx, dummy, dFy_dz = np.gradient (v, dx, dy, dz, axis=[1,0,2])
    dFz_dx, dFz_dy, dummy = np.gradient (w, dx, dy, dz, axis=[1,0,2])

    rot_x = dFz_dy - dFy_dz
    rot_y = dFx_dz - dFz_dx
    rot_z = dFy_dx - dFx_dy

    l = np.sqrt(np.power(u,2.0) + np.power(v,2.0) + np.power(w,2.0));

    m1 = np.multiply(rot_x,u)
    m2 = np.multiply(rot_y,v)
    m3 = np.multiply(rot_z,w)

    tmp1 = (m1 + m2 + m3)
    tmp2 = np.multiply(l,2.0)

    av = np.divide(tmp1, tmp2)

    return rot_x, rot_y, rot_z, av

mat = sio.loadmat('wind.mat')
x = mat['x']; y = mat['y']; z = mat['z']
u = mat['u']; v = mat['v']; w = mat['w']

rot_x, rot_y, rot_z, av = curl(x,y,z,u,v,w)


# plot a small area of the wind
i=5;j=7;k=8;S = 3
x1 = x[i-S:i+S, j-S:j+S, k-S:k+S]; 
y1 = y[i-S:i+S, j-S:j+S, k-S:k+S]; 
z1 = z[i-S:i+S, j-S:j+S, k-S:k+S];
u1 = u[i-S:i+S, j-S:j+S, k-S:k+S]; 
v1 = v[i-S:i+S, j-S:j+S, k-S:k+S]; 
w1 = w[i-S:i+S, j-S:j+S, k-S:k+S];

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.view_init(elev=47, azim=-145)

ax.quiver(x1, y1, z1, u1, v1, w1, length=0.05, color = 'black')

i=5;j=7;k=8;
x0=x[i,j,k]
y0=y[i,j,k]
z0=z[i,j,k]
cx0=rot_x[i,j,k]
cy0=rot_y[i,j,k]
cz0=rot_z[i,j,k]
ax.quiver(x0, y0, z0, 0, cy0, cz0, length=1.0, color = 'blue')

plt.show()

【讨论】：

【解决方案2】：

如前所述，存在必须实现某种离散 curl 运算符的问题。这大概是大气物理学中的一个常规问题，因此您可以查看有关这方面的教科书。

另一种方法可能是将样条拟合到数据中，以便您可以使用连续操作。例如

bspl = scipy.interpolate.SmoothBivariateSpline(x,y,z,s=0)

s 这是一个你应该使用的平滑因子；如果数据非常精确s=0 给出最好的结果；如果他们有很大的分散，你会想要一些平滑。现在你可以直接计算卷曲：

curl = bspl.integral(x0,x1,y0,y1) / ((x1-x0)*(y1-y0))

编辑：上面的表达式没有给出curl，但是基本思路是合理的。

【讨论】：

【解决方案3】：

正如@moarningsun 评论的那样，更改您调用np.gradient 的方式应该更正ValueError

dv_dx, dv_dy = np.gradient(vwind, dx,dy)
du_dx, du_dy = np.gradient(uwind, dx,dy)

您如何从文件中获得vwind 并不是特别重要，尤其是因为我们无权访问该文件。 vwind 的形状会很有用，尽管我们可以从错误消息中猜到这一点。错误中对(2,) 数组的引用是[dx,dy]。当您收到broadcasting 错误时，请检查各种参数的形状。

np.gradient 代码很简单，只是因为它可以处理 1、2、3d 和更高的数据而变得复杂。基本上它会进行类似的计算

(z[:,2:]-z[:,:-2])/2
(z[2:,:]-z[:-2,:])/2

用于内部值，1 项用于边界值。

我将把从渐变（或不）推导出curl 的问题留给其他人。

【讨论】：