利用 pandas 和 xarray 整理气象站点数据

平时用 xarray 库在处理 nc 格式的数据非常方便，但偶尔还是要用到一些站点数据来辅助分析，而站点数据一般都是用文本文件存储的，比如下图这种格式，从外到内的坐标依次是：年、月、站点、日

这种格式与CSV格式还有点不同，CSV格式是字段间用相同的符号隔开，而图中的文件可能是用 Fortran 写的，每个字段的长度固定为30个字符，此外，其中有不少特征值比如30XXX代表缺测/微量的情况，用Fortran处理也有不小的麻烦。

用Python处理这种文本列表就需要用上 pandas 库了， xarray 库就是基于 pandas 的，虽然天天在用 xarray ，但是这还是第一次正儿八经用 pandas 处理数据，就当做一次学习的过程啦。

一、 目标和步骤

将上图示例的文件处理为（站点，时间）坐标的 nc 格式数据，方便以后直接读取，主要有以下几个步骤：

• 将文本文件读取为 DataFrame 并将无效值替换为 Nan
• 将时间信息处理为 pandas 可用的时间坐标
• 将 DataFrame 进一步转换为 Dataset 并补充经纬度、站点名称信息

目标如图所示

二、 具体处理

1. 文件读取与预处理

• 导入所需的库

import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt

• 定义处理过程中的函数：
  • 处理时间坐标，利用 datetime 将整形的年、月、日转换为 pandas 的时间戳

def YMD_todatetime(ds):
  # 读取年月日数据，转换为Timestape，由于本质上还是遍历所有行，因此这个步骤最费时间
    import pandas as pd
    from datetime import datetime
    time = datetime(  # datetime 只接收整形参数，返回一个datetime类型的日期
        ds['年'].astype(int), ds['月'].astype(int), ds['日'].astype(int)
    )
    return pd.to_datetime(time)

• 具体的处理，包括特征值替换、插入日期列（利用 apply 函数逐行处理，这一步很费时间，暂时也没想到更快的方法），精度转换

def PreProcess(df_t):
  # 每读取一个文本文件做一步预处理
    df_t.loc[df_t['20-20时降水量'] >= 29999, '20-20时降水量'] = np.nan  # 替换掉所有特征值
    df_t.insert( # 插入日期列，此时并不以此为索引
        1, 'Date',df_t.iloc[:, 1:4].apply(YMD_todatetime, axis=1)
    ) 
    df_t.drop(columns=['年', '月', '日'], inplace=True, errors='raise')
    
    cols = ['平均本站气压', '平均风速', '平均气温', '日照时数',
            '平均水汽压', '20-20时降水量']  # 需要转换单位精度的变量 *0.1
    df_t[cols] = df_t[cols]/10.  # 转换精度
    return df_t

• 循环读取文件并处理 注意：
  • 不是用 pd.read_csv 而是用 pd.read_table 读取，选项sep='\s+'表示字段间至少有一个空格，\s 代表空白字符，+ 表示前面的字符至少重复一次（具体查看正则表达式的用法）
  • na_values 选项将把指定的值替换为 Nan
  • parse_dates=False 防止将某些字符解析为日期

StaDir = './Station/' # 文件路径，自定义
year = list(range(2012, 2014)) # 提取年份
usecols = ['区站号',   '年', '月', '日', '平均本站气压', '平均风速',
           '平均气温', '日照时数', '平均水汽压', '平均相对湿度', '20-20时降水量']  # 需要的变量
na_values = [32700, 32744, 32766]  # 分别代表 微量、空白、缺测，读取时替换为Nan

df = pd.DataFrame() # 先建立一个空表，然后append进去
for yr in year:
    print(yr)
    for i in [1, 2]:
        df_t = pd.read_table(
            StaDir+'SURF_CLI_CHN_MUL_DAY_1396725598869_' + f'{yr}_{i}.txt',
            sep='\s+',  parse_dates=False, na_values=na_values,
            engine='python', usecols=usecols, encoding='utf-8')
        df_t = PreProcess(df_t)  # 中间处理
        df = df.append(df_t, ignore_index=True,
                       verify_integrity=False, sort=None)
df['区站号'] = df['区站号'].astype(int) 
df.columns = ['StaNum', 'time', 'prec', 'pres',
              'wind', 'temp', 'vpres', 'rh', 'sunhr']  # 更改变量名
df.to_hdf('Station_test.hdf', key='df', mode='w') # 保存变量

得到的DateFrame如图所示

Dataframe信息

2. 转换为 nc 文件

到此为止，上面得到的文件已经可以用于基本的分析了，直接筛选站点、指定日期即可。

但是我自己还是习惯了直接用 xarray 处理文件，因此还是做了进一步处理。

首先读取站点的地理资料，比如下图这种格式

• 变量读取

df = pd.read_hdf('Station_test.hdf')

def LatLng_Rad2Dec(x): # 度分格式转为十进制
    d = int(str(x)[:-2])
    m = int(str(x)[-2:])
    decNum = d + m/60.0
    return decNum

stainfo = pd.read_excel(StaDir + 'Stations.xlsx')
stainfo['区站号'] = stainfo['区站号'] .astype(int).values
stainfo.index = stainfo['区站号']
# ind = stainfo['省份'] == '西藏'
ind = stainfo.index # 目标索引
stas = stainfo.loc[ind, '区站号'].values
name = stainfo.loc[ind, '台站名称'].apply(lambda x: x.replace(' ', ''))  # 去除台站名称中的空格
lat = stainfo.loc[ind, '纬度'].apply(LatLng_Rad2Dec)  # 转换为十进制小数
lon = stainfo.loc[ind, '经度'].apply(LatLng_Rad2Dec)
elev = stainfo.loc[ind, '海拔']/10.
prov = stainfo.loc[ind, '省份']

• nc 文件合并，沿着站点合并，取并集，个别站点缺少的时间坐标自动填充，变量填充为 Nan

ds_merge = xr.Dataset(
    data_vars={},
    coords={'station': (['station'], np.empty(shape=0, dtype=np.int64)),
            'time': (['time'], np.empty(shape=0, dtype=np.datetime64))}
) # 思路和上面读取dataframe一样，先建立一个空DataSet
n = 0
for s in stas: # 遍历每一个站点
    n = n+1
    print(f'\r{n}', end=' ')
    df_s = df[df['StaNum'] == s] # 站号为 s 的站点dataframe
    df_s.index = df_s['time']
    df_s = df_s.sort_index()
    ds_1 = xr.Dataset(data_vars=df_s.iloc[:, 1:])
    ds_1 = xr.concat([ds_1], pd.Index([s], name='station'))  # 增加一个站点维度
    ds_merge = xr.merge([ds_merge, ds_1])  # 沿站点合并
# 为Dataset补充地理信息
ds_merge['name'] = (('station'), name)
ds_merge['lon'] = (('station'), lon)
ds_merge['lat'] = (('station'), lat)
ds_merge['elev'] = (('station'), elev)
ds_merge['prov'] = (('station'), prov)

ds_merge.to_netcdf('Station_test.nc')

至此，文本格式的站点数据就转化成了便于读取和分析的 nc 数据了，结构如开头那张目标示意图所示。

三、 数据处理实例

1. 2012年夏季平均气温的空间分布

此例所用数据即上面生成的数据

ds = xr.open_dataset('Station_test.nc')
temp = ds['temp'].groupby('time.season').mean('time') # 季节平均
plt.scatter(ds.lon, ds.lat, s=6, c= temp.sel(season='JJA'), cmap='jet')

结果如下

2012夏季气温

2. 西藏站点平均的逐月风速距平序列

本文不含此例数据

ds = xr.open_dataset(
    '/Data/China753-1979-2014.nc')
indp = np.where(ds['prov'] == '西藏')[0]  # 筛选西藏的站点
TibetWind = ds['wind'][indp, :].mean('station')\
    .resample(time='M').mean()
TibetWindAnom = TibetWind.groupby(
    'time.month') - TibetWind.groupby('time.month').mean()  # 计算距平
TibetWindAnom.plot()  # 绘制距平序列
plt.axhline(0, c='k')
TibetWindAnom.rolling(dim={'time': 36}).mean().plot()  # 36个月滑动平均

结果如下

西藏风速距平

示例数据和代码

链接：https://pan.baidu.com/s/1yNYIIyg02kTyPw9HDqwddQ 提取码：tfuy