Numpy处理中国站点数据32766缺测值的方法

herenso

今天用到numpy来处理中国站点数据集，其中有很多记为32766的缺测值。以前一直用NCL，可以用var@_FillValue = 32766来非常方便地定义缺测值（毕竟是专门用于处理气象数据的），但今天查了半天numpy的文档，也没太弄明白numpy处理缺测值的方式。

结合搜索引擎和stackoverflow，我大概总结出自己的一套处理32766缺测值的方法，总结如下（下文涉及的Python语句均假设已经先import numpy as np）：

• 使用genfromtxt函数从文本文件中读取数据，genfromtxt函数中有参数可以设定缺测值，但这里我们不进行设置，原因最后说。
• 假设我们读取到一个两维numpy数组stn_prcp，第一维为不同站点，第二维为年份。我们要计算每个站点的气候平均值，则可以通过
  

  1. stn_prcp_clim = np.sum(stn_prcp, axis = 2)

  stn_prcp_clim = np.sum(stn_prcp, axis = 2)
  来计算。但这里的问题是，我们还没有设置缺测值，于是32766也会参与计算，所以我们通过
  

  1. stn_prcp[stn_prcp > 30000] = np.nan

  stn_prcp[stn_prcp > 30000] = np.nan
  来将stn_prcp数组中所有大于30000的异常值设定为缺测值（np.nan）。
• 这里，stn_prcp > 30000会生成一个与stn_prcp相同维度的布尔值序列，stn_prcp[stn_prcp > 30000]则会将stn_prcp中对应着True的元素设置为np.nan。
• 要说明的是，numpy中有两种缺测值的定义，一种是NA，一种是IGNORE，np.nan是前者，表示该数确定存在，只是不知道是多少，因此一个序列中只要有一个数是nan，则该序列的sum、mean等都无法准确求得。
• 因此，回到第2步，我们将大于30000的值设定为nan，会使含有缺测值的站点其气候平均值也为缺测nan，而非加入32766的错误结果。
• 下面问题是，如果我们要将生成的stn_prcp_clim输出，不想输出nan，而是还想将缺测的气候值输出为32766，可以利用isnan函数：
  

  1. stn_prcp_clim[np.ispan(stn_prcp_clim)] = 32766

  stn_prcp_clim[np.ispan(stn_prcp_clim)] = 32766
  。
• 同样的道理，np.isnan(stn_prcp_clim)会根据每个元素是否是np.nan生成一个布尔值序列，从而将stn_prcp_clim中的所有nan改为32766。
• 注意，由于np.nan的特殊性，这里只能用isnan来判断，而不能用stn_prcp_clim == np.nan。
  

　　以上就是我总结的缺测数据处理方法。另外补充一点，关于为什么不用genfromtxt来设定缺测值。在读取文本格式的数据时，使用genfromtxt函数。通过参数missing_values和usemask配合，即可设置缺测值。但要注意的是，numpy认为缺测值应该是类似于“nan”的字符串，而非用某一特定数值来表示，因此要写为

1. missing_values='32766'

missing_values='32766'

　　但即使这样，numpy仍按照IGNORE的原则处理读取到的数据，即忽略缺测值，将剩下的数计算sum或者mean。另外还有一点，就是如果你在genfromtxt中设置了names，会改变得到的数组类型，stn_prcp中每一行的类型会由numpy.ndarray变成numpy.void，此时missing_values=32766又可以不加引号了……

　　总体上感觉numpy这个包有时候解决问题的思路不是很清晰，至少文档非常杂乱，看了半天不得要领。当然Python的科学计算这部分我很少接触，希望有经验的朋友指教！

阿阿飞飞

有两种解决办法：
1、含有nan的数组运算时，使用专用函数，如np.nansum，np.nanmean这些便会忽略nan
stn_prcp[stn_prcp>30000]=np.nan
np.nansum(stn_prcp, axis = 2)
2、使用mask函数：
stn_prcp=np.ma.masked_greater_equal(stn_prcp,30000)
np.sum(stn_prcp, axis = 2)
如此即忽略了缺省值

kongfeng0824

谢谢分享经验

ginanony

{:eb502:}{:eb502:}{:eb502:}
感谢分享

Arsenal

试试np.ma.masked_equal

亽小笨蜑■ɑ

在NCL中如果缺测值是空白，怎么处理呢 谢谢

fengkp

很好的经验分享！

fengkp

很好的经验分享！

wanghp

不错，感谢分享

candyzh

{:eb502:}  不错的经验分享