CMIP降尺度前期处理1：将下载降水历史数据进行nc合并【20250413】

import os
import numpy as np
import xarray as xr
import rasterio
from datetime import datetime
import pandas as pd
import glob

# 输入路径
input_dir = "J:/data/PRE_shed3"
output_nc = "J:/data/precip_1981-2000_1.nc"  # 输出文件路径
# ----------------------------------------
# 关键修改1: 定义时间解析函数（支持文件名中的多位数期数）
# ----------------------------------------
# 文件命名规则示例: precip_198101.tif → 1981年第1期
def parse_time_from_filename(fname):
    # 提取文件名中的年份和期数（例如："precip_198123.tif" → year=1981, period=23）
    parts = os.path.basename(fname).split("_")[1].split(".")[0]
    year = int(parts[:4])
    period = int(parts[4:])  # 直接提取剩余部分作为期数
    
    # 定义每期的起始日期（假设为每年1月1日开始，每16天一期）
    start_date = datetime(year, 1, 1)
    time = start_date + pd.DateOffset(days=(period-1)*16)  # 需验证时间间隔
    return time

# 获取所有TIFF文件路径并按时间排序
all_files = sorted(glob.glob(os.path.join(input_dir, "precip_*.tif")))

# 逐个读取文件并构建DataArray列表
data_arrays = []
for file in all_files:

    fname = os.path.basename(file)
    time = parse_time_from_filename(fname)
    
    # 读取TIFF数据
    with rasterio.open(file) as src:
        left, bottom, right, top = src.bounds
        # 生成经度（横向）和纬度（纵向）坐标
        transform = src.transform
        lon = transform.c + np.arange(src.width) * transform.a  # 经度：从左到右递增
        lat = transform.f + np.arange(src.height) * transform.e  # 纬度：从上到下递减（需反转） 
        lat = lat[::-1]  # 反转纬度方向
        precip = src.read(1) [::-1, :]  # 读取第一个波段
        
    # 创建DataArray（假设经纬度信息已知或从TIFF读取）
    da = xr.DataArray(
        precip[np.newaxis, ...],  # 添加时间维度
        dims=["time", "lat", "lon"],
        coords={"time": [time],"lat": lat, "lon": lon},
           # "lat": np.arange(precip.shape[0])*transform.e + transform.c,  # 根据实际坐标调整
           # "lon": np.arange(precip.shape[1])*transform.a + transform.f
        name="precip"
    )
    data_arrays.append(da)

# 合并所有时间步
ds = xr.concat(data_arrays, dim="time")

# ----------------------------------------
# 关键修改3: 统一坐标系和元数据
# ----------------------------------------
# 添加坐标系属性（假设为WGS84）
# 添加元数据
ds.lat.attrs = {"units": "degrees_north", "long_name": "latitude"}
ds.lon.attrs = {"units": "degrees_east", "long_name": "longitude"}
ds.attrs = {
    "title": "准噶尔盆地降水数据 (1981-2000)",
    "units": "mm/day",
    "time_resolution": "16-day composite",
    "source": "TIFF files from PRE_shed3"
}

# 保存为NetCDF
# 直接保存，不使用压缩
ds.to_netcdf(output_nc)
print(f"数据已保存至 {output_nc}")

ArcGIS对比：

说明空间完全吻合，可以进行下一步了！

合并了1981-2000年，再合并2001-2020年【修订代码年份】，然后合并2001-2014年，2015-2020年作为验证。