代码收藏家技术教程 2025-01-30

【代码分享】Python根据站点位置批量提取nc文件数据

ERA5数据是常用于研究气候变化的数据库，提供了逐时，逐日，逐月不同时间尺度的数据，包含0.1°及0.25°空间分辨率，下载的文件为nc（Netcdf）数据格式。用arcmap，创建netcdf栅格图层能够打开，但在提取多站点对应的数值时处理效率较低。用Python可实现nc文件的批量读取与数据处理。下这边，以0.25°空间分辨的逐日地面2米温度为例，来介绍ERA5数据的下载及处理，并附带数据链接。

数据连接：

百度网盘：

通过网盘分享的文件：ERA5_export.rar
链接: https://pan.baidu.com/s/1XorpZXXSiCAddxkKbSZDhw 提取码: 5vfq
–来自百度网盘超级会员v2的分享

ERA5气象数据下载：

ERA5气象数据下载链接：Catalogue — Climate Data Store (copernicus.eu)

1.找到所需数据集

本次选择的数据集为ERA5 post-processed daily statistics on single levels from 1940 to present（ERA5 post-processed daily statistics on single levels from 1940 to present (copernicus.eu)）

2.选择变量、时间范围、空间范围及数据格式

范围覆盖全国，也可根据自己的研究区去填写。

站点数据提取：

本次案例需要提取的站点及温度数据如下。站点数据以区站号为标识符，每个站提取的文件以此命名。温度nc文件包含1960-1962年逐日数据，单位为K。

1.导入所需的库，填写相关的文件路径。

2.先查看nc文件信息，包括变量名、时间及经纬度名称。Spyder可以运行局部代码，还是很方便的。

3.运行当前单元，查看变量信息，主要的单位，根据自己的需求决定是否需要进行单位转换。

以下为单个站点的查询方式，包含双线性插值提取及提取站点最近栅格值两种方式，红框处的经纬度名要与先前查询的nc文件变量的经纬度名一致。后续逐站点提取即可。

完整代码展示：

import os
import xarray as xr
import pandas as pd
from tqdm import tqdm  
#根据文件夹路径创建文件夹
def creat_dir(filename):
    if not os.path.exists(filename): 
        os.makedirs(filename)
#获取某个时间的字符串格式 '1961-06-15T12:00:00.000000000' --> 1961
def get_timestr(time,form):                   
    time = pd.to_datetime(str(time))
    timestr = time.strftime(form)
    return timestr
#%%
outputdir = r'F:\download\ERA5_daily\ERA5_export\result' #结果保存文件夹 
inputdir = r'F:\download\ERA5_daily\ERA5_export\nc_data'  #需提取nc文件所在文件夹
sid_path = r"F:\download\ERA5_daily\ERA5_export\sids_table.xlsx" #需提取站点所在表
sid_table = pd.read_excel(sid_path)
id_name,slon_name,slat_name = '区站号','Lon','Lat' #Excel表中 站点标识，经度名，纬度名

file_list = os.listdir(inputdir)
file_name = file_list[0]
file_path = os.path.join(inputdir,file_name)
data = xr.open_dataset(file_path)
print(data)

#%%
variable_name = 't2m' #变量名
lon_name,lat_name = 'longitude','latitude' #nc文件中经纬度名称
time_name = 'valid_time' #时间维度名
variable_data = data[variable_name ]
print(variable_data) #查看变量单位
#%%
for sid_index in tqdm(range(sid_table.shape[0])):
    # sid_index = 0
    #获取单站点信息
    loni,lati = sid_table[slon_name][sid_index],sid_table[slat_name][sid_index]
    id_namei = sid_table[id_name][sid_index]
    #逐文件提取数据
    time_list,var_values = [],[]
    for file_name in file_list:
        # file_name = file_list[0]
        file_path = os.path.join(inputdir,file_name)
        data = xr.open_dataset(file_path)
        variable_data = data[variable_name]
        #单位转换
        variable_data = variable_data - 273.15 
        #查询数据
        #方法一：双线性插值提取
        sid_data = variable_data.interp(longitude=loni,latitude=lati)
        #方法二：提取最近的栅格值
        # sid_data = variable_data.sel(longitude =loni,latitude=lati,method='nearest') 
        #save
        time_list.extend(sid_data[time_name].values)
        var_values.extend(sid_data.values) 
    #save to table    
    sid_value_table = pd.DataFrame({'time':time_list,variable_name:var_values})
    table_path = os.path.join(outputdir,f'{variable_name}_id_{id_namei}.xlsx')
    sid_value_table.to_excel(table_path,index=False)