用于识别低空急流的Python算法(气象学)

Question

低空急流是地面附近的局部风速最大值.Baas (2009)定义了一个低空急流如下：

大气最低500米处风速剖面的最低最大值，至少比下一个最低风速高2米/S和25% (.)如果高于该最小风速的风速增加小于1m/S，则忽略最小风速，然后再降至低于最小风速的值。当没有最小值时，将射流上方风速剖面的最低值作为最小值。

我有一个二维的numpy阵列，它包含了几个测量高度和时间的风速值。我已经编写了一个算法来确定每个时间步骤的风廓线是否符合这个标准。但是，我认为我的代码可以改进。现在我使用的是嵌套循环和很多if语句。另外，我需要捕获一个索引错误，因为我需要检查数组中的“next”值是低还是高。

我希望从更有经验的程序员那里得到一些反馈。这是我的密码：

def findLLJ(wspd_matrix, time, hgt):
  '''Find low-level jets according to Peter Baas' algorithm
     input: wind speed matrix with time as first axis and height as second axis
            time and height arrays
     output: jet index array, jet height array, jet speed array
     use: jetidx,jethgt,jetspd = findLLJ(wsp_matrix,time,height)'''
  #Allocate output arrays:
  jet = np.zeros(time.shape,dtype=bool)
  jetspd = np.zeros(time.shape)
  jethgt = np.zeros(time.shape)
  #Shape of input array:
  ntime,nlev = wspd_matrix.shape
  #Consistency check:
  if not (ntime == len(time) and nlev == len(hgt)):
    raise Exception("Shapes of input arrays are not consistent")
  #Loop over all wind profiles at given times
  for t,prof in enumerate(wspd_matrix):
    maxflg = False
    minflg = False
    maxspd = 0.
    #Loop over all vertical points in the profile
    for i,x in enumerate(prof):
      #Find out whether vertical point is a maximum
      if np.any(x>(prof[i:]+2.5)) and np.any(x>1.25*prof[i:]) and x>maxspd:
        maxflg = True
        maxind = i
        maxspd = x
        print 'Maximum found at time',t,'and level',i,'with speed',x,'m/s'
      #Find out whether vertical point is a minimum
      try:
        if x<0.8*maxspd and x<maxspd-2.5 and x<prof[i+1] and not prof[i+2]<x:
          minflg = True
          minind = i
          minspd = x
          print 'Minimum found at time',t,'and level',i,'with speed',x,'m/s'
      except IndexError:
        if i == nlev-2 and x<0.8*maxspd and x<maxspd-2.5 and x<(prof[nlev-1]-1):
          minflg = True
          minind = i
          minspd = x
          print 'Minimum found at time',t,'and level',i,'with speed',x,'m/s'   
        else:
          minflg = True
          minind = nlev-1
          minspd = prof[nlev-1] # this should be changed to minumum of the profile, but that might return NaN, need to find an update for that.
          print 'No minumum found, using lowest value from profile'
      #Determine wheter a LLJ occured
      if maxflg and minflg:
        print 'Low-level jet found at time',t,'and height',hgt[maxind],'m with speed',maxspd,'m/s and minimum at',hgt[minind],'m with speed',minspd,'m/s'
        jet[t] = True
        jetspd[t] = maxspd
        jethgt[t] = hgt[maxind]
        break
      else:
        continue
  if maxflg and not minflg:
    print 'probably something wrong with NaNs'
  #End of loop, next profile...
  return jet,jethgt,jetspd

Answer 1

命名变量

• 如果wspd_matrix被称为wind_speed_matrix，那么它将更容易阅读。
• 如果jetspd被称为jet_speed，那么它将更容易阅读。
• 如果jethgt被称为jet_height，那么它将更容易阅读。

删除代码

中的数值。

0.8和2.5经常出现在代码中。对于一个没有领域知识的人来说，很难知道他们代表什么。如果你因为某种原因不得不改变它们，这是很容易出错的。为什么不只是在函数的开头命名它们(例如max_multiplier = 0.8)？

冗余条件

if x<0.8*maxspd and x<maxspd-2.5

在代码中出现两次。它代表了一个明确的数学功能：

if x<min(0.8*maxspd , maxspd-2.5)

更容易读懂(并在脑海中描绘)

冗余码

这三种指纹需要做些什么：

print 'Maximum found at time',t,'and level',i,'with speed',x,'m/s'
print 'Minimum found at time',t,'and level',i,'with speed',x,'m/s' 
print 'Minimum found at time',t,'and level',i,'with speed',x,'m/s' 

也许是

print_local_optimum_found(type_optimum, at_time, at_level, at_speed)

检查输入

if maxflg and not minflg:
    print 'probably something wrong with NaNs'

为什么要等到程序结束才意识到发生了什么错误呢？如果你知道NAs会带来麻烦，就检查一下之前是否有NAs存在。

算法本身

现在很难建议对同一函数中的每一段代码进行改进。如果您修复了所有这些，这些改进将变得更加可读性！

Answer 2

我的建议，在RUser4512概述的同一条线上，将开始专注于更好地命名变量。我将为prof、nlev和ntime添加更好的命名，这样就更容易知道它们包含了什么

稍后，Point类可能很方便，它有一个函数来确定点是否是最大值，而另一个函数则用来知道它的最小值

您可以减少主要代码块：

if np.any(x>(prof[i:]+2.5)) and np.any(x>1.25*prof[i:]) and x>maxspd:
    maxflg = True
    maxind = i
    maxspd = x
    print 'Maximum found at time',t,'and level',i,'with speed',x,'m/s'
  #Find out whether vertical point is a minimum
  try:
    if x<0.8*maxspd and x<maxspd-2.5 and x<prof[i+1] and not prof[i+2]<x:
      minflg = True
      minind = i
      minspd = x
      print 'Minimum found at time',t,'and level',i,'with speed',x,'m/s'
  except IndexError:
    if i == nlev-2 and x<0.8*maxspd and x<maxspd-2.5 and x<(prof[nlev-1]-1):
      minflg = True
      minind = i
      minspd = x
      print 'Minimum found at time',t,'and level',i,'with speed',x,'m/s'   
    else:
      minflg = True
      minind = nlev-1
      minspd = prof[nlev-1] # this should be changed to minumum of the profile, but that might return NaN, need to find an update for that.
      print 'No minumum found, using lowest value from profile'
  #Determine wheter a LLJ occured
  if maxflg and minflg:
    print 'Low-level jet found at time',t,'and height',hgt[maxind],'m with speed',maxspd,'m/s and minimum at',hgt[minind],'m with speed',minspd,'m/s'

这样的事情：

if point.is_maximum() and point.is_minimum():
    print 'Low-level jet found at time',t,'and height',hgt[point.maxind],'m with speed',point.maxspd,'m/s and minimum at',hgt[point.minind],'m with speed',point.minspd,'m/s'

Point.is_maximum()的示例：

def is_maximum(self):
    if np.any(self.x>(self.prof[self.i:]+2.5)) and np.any(self.x>1.25*self.prof[self.i:]) and self.x>self.maxspd:
    self.maxflg = True
    self.maxind = i
    self.maxspd = x
    print 'Maximum found at time',self.prof.t,'and level',self.i,'with speed',self.x,'m/s'

注意，您必须使用一些参数初始化Point对象，并且创建Profile对象也可能是明智的