NCDC气象数据的提取与处理（四）：python批量读取、写入nc数据经纬度格点数值

Table of Contents

1.问题描述：

NCDC的站点数据处理在之前三节里已经介绍过了，但是NCDC的就那么几种数据可能不能满足日常使用，比如说辐射数据他就没有。这时候我们找到其他类型数据要和它原有数据融合，比如本例找的nc格式数据。

2.思路：

本例所用数据集是网格化的面尺度数据，之前处理的NCDC是点数据，最简单的处理方法就是把对应站点所在经纬度找到，再与nc格式数据的格点位置匹配，提取对应位置的nc数据，放入表格。

3.实现过程：

3.1格点位置匹配

def matching(lat_value,lon_value,data): # 文件编码方式[lats,lons]，data:需要输入xr.open_dataset()
    '''
    经纬度与格点位置匹配
    输入经纬度,nc数据
    输出所在格点位置对应索引值 lon_idx：经度 lat_idx：纬度
    '''
    # 获取经度和纬度的坐标值
    lons = data['lon'].values
    lats = data['lat'].values
    
    # 获取经纬度坐标的索引值
    lat_idx = np.abs(lats - lat_value).argmin()
    lon_idx = np.abs(lons - lon_value).argmin()
    
    return lat_idx,lon_idx

上述函数要求输入站点所在的经度、纬度，以及要读的nc文件

3.2写入表格

我直接放完整示例代码上来，大家按需修改。

# 读取站点经纬度。将对应格点数值写入dataframe
import xarray as xr
import netCDF4 as nc
import numpy as np
import pandas as pd
import os
import openpyxl


# 定义函数
def matching(lat_value,lon_value,data): # 文件编码方式[lats,lons]，data:需要输入xr.open_dataset()
    '''
    经纬度与格点位置匹配
    输入经纬度,nc数据
    输出所在格点位置对应索引值 lon_idx：经度 lat_idx：纬度
    '''
    # 获取经度和纬度的坐标值
    lons = data['lon'].values
    lats = data['lat'].values
    # 获取经纬度坐标的索引值
    lat_idx = np.abs(lats - lat_value).argmin()
    lon_idx = np.abs(lons - lon_value).argmin()    
    return lat_idx,lon_idx


# 输入站点列表
stations = pd.read_excel('../../Desktop/YellowRiver_station_ID.xlsx') # 站点列表位置，这句按自己需求修改***************
station_code_list = [str(x)[:6] for x in stations["STATION_ID"]]  # 没看懂的复习一下     NCDC气象数据的提取与处理（二）
LAT_list = [x for x in stations["LAT"]]
LON_list = [x for x in stations["LON"]]
station_data = zip(station_code_list,LAT_list,LON_list)
station_data = list(station_data)     # 包含每个站点的编号、纬度、经度，如coordinate[0]为站点编号


# 设置输入文件列表
in_rootdir = r"E:\China_station_Srad" # 设置输入根目录， 这句按自己需求修改********************************************
file_list = [] 
for root, dirs, files in os.walk(in_rootdir):           # 获取输入根目录下所有文件
    for file in files:                                  # 遍历所有文件名
        if os.path.splitext(file)[-1] == '.nc':         # 获取所有nc文件绝对路径列表
            file_list.append(os.path.join(root, file))   # 拼接出 绝对路径 并放入列表
            
# 坐标匹配
print("*"*15,"坐标匹配","*"*15)       # 绘制分割线
data = xr.open_dataset(file_list[0])
coord_idx_list = []
for coordinate in station_data:      # coordinate[i],i=0,1,2；依次表示站点编号、纬度、经度
    coordinate_index = matching(lat_value=coordinate[1],lon_value=coordinate[2],data=data)
    coord_idx_list.append((coordinate[0],coordinate_index))       # 形如：[('534630', (258, 418)),…]
print("*"*15,"匹配完成","*"*15)      # 绘制分割线


# 读取nc数据，写入xlsx文件
for file_name in file_list:          
    data = nc.Dataset(file_name)
    times = nc.num2date(data["time"][:],data["time"].units)       # 时间的格式转换,得到一个数组
    dt = [times[i].strftime() for i in range(data["time"].size)] # 获取时间序列

    # 创建一个新dataframe
    df = pd.DataFrame(columns=['Date','srad'])                # 这句按自己需求修改***************************************
    df['Date'] = pd.to_datetime(dt)
    df = df.set_index('Date')

    # 设置输出路径
    out_rootdir = r"E:\Srad_TEST_2"        # 设置输出根目录 ，这句按自己需求修改*****************************************
    workbook = openpyxl.Workbook()       # 创建Workbook对象
    moy = os.path.splitext(file_name)[0][-6:]                # 获取年和月moy:month of year
    print("-"*25,moy,"-"*25)             # 绘制分割线，年月
    outfile_dir = f'{out_rootdir}\\test{moy}.xlsx'
    workbook.save(outfile_dir)           # 保存空的xlsx文件

    for idx in coord_idx_list:     # idx[0]表示站点编号、idx[1][0],idx[1][1]分别表示纬度、经度
        print("*"*15,idx[0],"*"*15) # 绘制分割线，站点代号
        location = idx[1] # 获得格点坐标location，包含（纬度，经度）
        tempfile = []
        for i in range(data["time"].size):         
            tempfile.append(data["srad"][i][location])        # 获得对应格点的属性值，并存入列表
        df["srad"] = tempfile                                 # 对应格点的属性值存入dataframe
        srad_mean = df.resample('d').mean()                   # mean：取属性值的日均值
        writer = pd.ExcelWriter(outfile_dir, engine='openpyxl',mode='a') # 将excel写入对象writer
        """
        将各站点的数据表写入Excel中的
        sheet1、sheet2、sheet3……，
        sheet以站点编号命名
        """
        srad_mean.to_excel(writer, sheet_name=idx[0])  # idx[0]为站点编号
        writer.save()                   # 保存读写的内容
    wb = openpyxl.load_workbook(outfile_dir)      
    del wb['Sheet']                     # 删除新建Excel时默认的空表
    wb.save(outfile_dir)

4.运行效果

4.1打包站点信息

首先我有一个xlsx文件，名为“YellowRiver_station_ID.xlsx’，之前文章出现过，存放的是NCDC对应站点的代号和经纬度，（还有高程，这里用不到）形式如下：

我这边用了一次 list() 和 zip() 函数的套用，总之最后打包效果就是[(a1, b1,c1 ), (a2, b2,c2 ), …]，abc分别表示站点代号，经度，纬度。

4.2读取nc文件列表

我这边还有一个文件夹，名为“China_station_Srad”。里面存放着这次要读的nc文件

我这里做了一个 file_list ，把文件夹里所有nc格式文件的绝对路径存进去了

对于每一条绝对路径，做如下操作

4.3提取对应格点的nc数据

创建一个dataframe，形式如下：

4.4数据写入

将各站点的数据表写入xlsx中的sheet1、sheet2、sheet3……，sheet以站点编号命名

最后，文件名称以年月命名，形如 201812

如有错漏希望各位大佬不吝指教。

文章出处登录后可见！

已经登录？立即刷新