卫星数据Python中u、v风分量的卷曲计算

Question

我不知道如何获取卫星数据中风的u和v分量的导数。我想我可以用数字梯度以这种方式：

    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也不正确。如果有人能在卫星数据中处理这些类型的计算，我将非常感激。谢谢您！！在

Answer 1

如前所述，存在一个必须实现某种离散curl算子的问题。这大概是大气物理学中的一个常规问题，所以你可以查阅一下教科书。在

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

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

s这里有一个平滑因子，您应该使用它；如果数据非常精确，s=0会产生最好的结果；如果它们有大量的分散，您将需要一些平滑。现在可以直接计算旋度：

^{pr2}$

编辑： 上面的表达式没有给出旋度，但基本思想是正确的。在

Answer 2

正如@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和更高的数据。基本上它在计算

^{pr2}$

对于内部值，1项步骤用于边界值。在

我将把从梯度（或不从）导出curl的问题留给其他人。在