【问题标题】:Calculating wind divergence of u and v using Python, np.gradient使用 Python,np.gradient 计算 u 和 v 的风散度
【发布时间】:2015-11-27 20:01:33
【问题描述】:

我对 Python 非常陌生,目前正在尝试复制我以前使用 Grad 的绘图等。我想使用来自 netCDF 气候模型文件的 u 和 v 风场(仅按特定湿度 q 缩放)计算每个网格框的散度。

从无休止的搜索中,我知道我需要使用 np.gradient 和 np.sum 的某种组合,但找不到正确的组合。我只知道要“手工”完成,计算将是 divg = dqu/dx + dqv/dy 我知道下面的内容是错误的,但这是迄今为止我所拥有的最好的......

nc = Dataset(ifile)
q = np.array(nc.variables['hus'][0,:,:])
u = np.array(nc.variables['ua'][0,:,:])
v = np.array(nc.variables['va'][0,:,:])
lon=nc.variables['lon'][:]
lat=nc.variables['lat'][:]

qu = q*u
qv = q*v  

dqu/dx, dqu/dy = np.gradient(qu, [dx, dy])
dqv/dx, dqv/dy = np.gradient(qv, [dx, dy])

divg = np.sum(dqu/dx, dqv/dy)

这会给出错误“语法错误:无法分配给运算符”。

任何帮助将不胜感激。

【问题讨论】:

  • 附带说明,gradient 返回d_row, d_column。这取决于您的输入数组,但这通常是 dy, dx,而不是 dx, dy

标签: python numpy netcdf


【解决方案1】:

如果您的网格是高斯网格并且文件中的风名称是“u”和“v”,您也可以直接使用 cdo 计算散度:

cdo uv2dv in.nc out.nc

更多详情请见https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-6850002.13.2

【讨论】:

  • 虽然这个答案可能对 OP 有帮助,但它没有回答为什么代码会抛出语法错误并且与 python 无关的问题
  • 如果有人想知道如何计算涡度和散度,他们很可能会找到这个页面,并且可能不知道他们可以从命令行轻松地做到这一点......这取决于是否考虑堆栈-exchange 仅作为针对特定问题或更广泛的信息发布特定调试建议的一种方式。我认为后者更有用。其他人已经认为它有用,你却认为它太笼统了……这就是投票系统的用途。
  • 这个 cdo 命令只计算风矢量场的散度。如果一个向量场为例如海洋或冰的速度,需要计算它们的散度吗?
  • 这个函数是相当硬连线的 - 在这里看到这个线程code.mpimet.mpg.de/boards/2/topics/4647 但我想你可能会通过用 nco 重命名 advective 字段来捏造它?不幸的是,CDO 没有更通用的通量字段运算符...
【解决方案2】:

尝试删除 [dx, dy]。

   [dqu_dx, dqu_dy] = np.gradient(qu)
   [dqv_dx, dqv_dy] = np.gradient(qv)

还要指出您是否正在重新创建情节。 numpy 的梯度在 1.82 和 1.9 之间变化。这对在 python 中重新创建 matlab 图有影响,因为 1.82 是 matlab 方法。我不确定这与 Grad 有什么关系。这是两者的措辞。

1.82

"梯度是使用内部的中心差异计算的 和边界处的第一个差异。因此返回的梯度有 与输入数组的形状相同。”

1.9

“梯度是使用内部的二阶精确中心差异边界处的一阶差异或二阶精确的单边(前向或后向)差异计算的 >。因此返回的梯度与输入数组具有相同的形状。”

1.82的梯度函数来了。

def gradient(f, *varargs):
"""
Return the gradient of an N-dimensional array.

The gradient is computed using central differences in the interior
and first differences at the boundaries. The returned gradient hence has
the same shape as the input array.

Parameters
----------
f : array_like
  An N-dimensional array containing samples of a scalar function.
`*varargs` : scalars
  0, 1, or N scalars specifying the sample distances in each direction,
  that is: `dx`, `dy`, `dz`, ... The default distance is 1.


Returns
-------
gradient : ndarray
  N arrays of the same shape as `f` giving the derivative of `f` with
  respect to each dimension.

Examples
--------
>>> x = np.array([1, 2, 4, 7, 11, 16], dtype=np.float)
>>> np.gradient(x)
array([ 1. ,  1.5,  2.5,  3.5,  4.5,  5. ])
>>> np.gradient(x, 2)
array([ 0.5 ,  0.75,  1.25,  1.75,  2.25,  2.5 ])

>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float))
[array([[ 2.,  2., -1.],
       [ 2.,  2., -1.]]),
array([[ 1. ,  2.5,  4. ],
       [ 1. ,  1. ,  1. ]])]

"""
  f = np.asanyarray(f)
  N = len(f.shape)  # number of dimensions
  n = len(varargs)
  if n == 0:
      dx = [1.0]*N
  elif n == 1:
      dx = [varargs[0]]*N
  elif n == N:
      dx = list(varargs)
  else:
      raise SyntaxError(
              "invalid number of arguments")

# use central differences on interior and first differences on endpoints

  outvals = []

# create slice objects --- initially all are [:, :, ..., :]
  slice1 = [slice(None)]*N
  slice2 = [slice(None)]*N
  slice3 = [slice(None)]*N

  otype = f.dtype.char
  if otype not in ['f', 'd', 'F', 'D', 'm', 'M']:
      otype = 'd'

# Difference of datetime64 elements results in timedelta64
  if otype == 'M' :
    # Need to use the full dtype name because it contains unit information
      otype = f.dtype.name.replace('datetime', 'timedelta')
  elif otype == 'm' :
    # Needs to keep the specific units, can't be a general unit
      otype = f.dtype

  for axis in range(N):
    # select out appropriate parts for this dimension
      out = np.empty_like(f, dtype=otype)
      slice1[axis] = slice(1, -1)
      slice2[axis] = slice(2, None)
      slice3[axis] = slice(None, -2)
    # 1D equivalent -- out[1:-1] = (f[2:] - f[:-2])/2.0
      out[slice1] = (f[slice2] - f[slice3])/2.0
      slice1[axis] = 0
      slice2[axis] = 1
      slice3[axis] = 0
    # 1D equivalent -- out[0] = (f[1] - f[0])
      out[slice1] = (f[slice2] - f[slice3])
      slice1[axis] = -1
      slice2[axis] = -1
      slice3[axis] = -2
    # 1D equivalent -- out[-1] = (f[-1] - f[-2])
      out[slice1] = (f[slice2] - f[slice3])

    # divide by step size
      outvals.append(out / dx[axis])

    # reset the slice object in this dimension to ":"
      slice1[axis] = slice(None)
      slice2[axis] = slice(None)
      slice3[axis] = slice(None)

  if N == 1:
      return outvals[0]
  else:
      return outvals

【讨论】:

    【解决方案3】:

    尝试类似:

    dqu_dx, dqu_dy = np.gradient(qu, [dx, dy])
    dqv_dx, dqv_dy = np.gradient(qv, [dx, dy])
    

    你不能分配给python中的任何操作;其中任何一个都是语法错误:

    a + b = 3
    a * b = 7
    # or, in your case:
    a / b = 9
    

    更新

    遵循Pinetwig 的评论:a/b 不是有效的标识符名称;它是一个运算符(的返回值)。

    【讨论】:

    • 澄清一下:不能在变量名中使用正斜杠,因为在 Python 中,正斜杠表示除法。
    • 谢谢 - 这解决了这个问题,但现在我收到错误 NameError: name 'dx' is not defined
    • 这是因为我没有将 dx 定义为 lon/dy 的变化作为 lat 的变化吗?
    猜你喜欢
    • 2018-09-05
    • 1970-01-01
    • 2015-12-18
    • 2012-07-11
    • 2016-10-03
    • 2021-11-25
    • 1970-01-01
    • 1970-01-01
    • 2020-05-29
    相关资源
    最近更新 更多