python - 计算卫星数据中 u 和 v 风分量的旋度 - Python

Question

我不确定如何对卫星数据中风的 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 也不正确。如果有人对在卫星数据中进行这些类型的计算有好主意，我将不胜感激。谢谢!!

Answer 1

正如 @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] 的引用。当您遇到 广播 错误时，请检查各种参数的形状。

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

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

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

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