Python | 1901—2025年全球CRU气象数据下载与栅格提取

https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.10/cruts.2604091129.v4.10/

import osimport requests# 可选变量：# pre = 降水# tmp = 平均气温# tmn = 最低气温# tmx = 最高气温# dtr = 昼夜温差# pet = 潜在蒸散发# vap = 水汽压# wet = 湿日频率# cld = 云量# frs = 霜冻日频率var = input("请输入变量（如 pre）：").strip()start_year = int(input("请输入开始年份："))end_year = int(input("请输入结束年份："))base = "https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.10/cruts.2604091129.v4.10"save_dir = rf"E:\CRU\{var}"os.makedirs(save_dir, exist_ok=True)# CRU TS v4.10 的年份分段periods = [    (1901, 1910), (1911, 1920), (1921, 1930),    (1931, 1940), (1941, 1950), (1951, 1960),    (1961, 1970), (1971, 1980), (1981, 1990),    (1991, 2000), (2001, 2010), (2011, 2020),    (2021, 2025)]for y1, y2 in periods:    # 只下载与目标年份有重叠的数据    if y2 < start_year or y1 > end_year:        continue    name = f"cru_ts4.10.{y1}.{y2}.{var}.dat.nc.gz"    url = f"{base}/{var}/{name}"    path = os.path.join(save_dir, name)    if os.path.exists(path):        print("已存在，跳过：", name)        continue    print("正在下载：", name)    with requests.get(url, stream=True) as r:        r.raise_for_status()        with open(path, "wb") as f:            for chunk in r.iter_content(1024 * 1024):                if chunk:                    f.write(chunk)print("下载完成：", save_dir)

import osimport gzipimport shutilimport tempfileimport numpy as npimport rasteriofrom netCDF4 import Dataset, num2datefrom rasterio.transform import from_originfolder = r"E:\CRU\pre"monthly_dir = os.path.join(folder, "monthly")yearly_dir = os.path.join(folder, "yearly")os.makedirs(monthly_dir, exist_ok=True)os.makedirs(yearly_dir, exist_ok=True)year_data = {}files = [f for f in os.listdir(folder)         if f.endswith(".nc") or f.endswith(".nc.gz")]for file in files:    path = os.path.join(folder, file)    temp_nc = None    # 如果是 .nc.gz，临时解压    if file.endswith(".gz"):        temp = tempfile.NamedTemporaryFile(suffix=".nc", delete=False)        temp.close()        temp_nc = temp.name        with gzip.open(path, "rb") as f_in, open(temp_nc, "wb") as f_out:            shutil.copyfileobj(f_in, f_out)        nc_path = temp_nc    else:        nc_path = path    with Dataset(nc_path) as nc:        lon = nc.variables["lon"][:]        lat = nc.variables["lat"][:]        time = nc.variables["time"]        pre = nc.variables["pre"]        dates = num2date(            time[:],            units=time.units,            calendar=getattr(time, "calendar", "standard")        )        dx = abs(lon[1] - lon[0])        dy = abs(lat[1] - lat[0])        transform = from_origin(            lon.min() - dx / 2,            lat.max() + dy / 2,            dx, dy        )        for i, date in enumerate(dates):            data = pre[i, :, :]            if np.ma.isMaskedArray(data):                data = data.filled(np.nan)            data = np.array(data, dtype="float32")            # CRU纬度为南→北，GeoTIFF需北→南            if lat[0] < lat[-1]:                data = np.flipud(data)            # ---------- 月度栅格 ----------            monthly_file = os.path.join(                monthly_dir,                f"pre_{date.year}{date.month:02d}.tif"            )            out_data = np.where(np.isnan(data), -9999, data)            with rasterio.open(                monthly_file,                "w",                driver="GTiff",                height=data.shape[0],                width=data.shape[1],                count=1,                dtype="float32",                crs="EPSG:4326",                transform=transform,                nodata=-9999,                compress="lzw"            ) as dst:                dst.write(out_data, 1)            print("月度：", monthly_file)            # 保存用于年度求和            year_data.setdefault(date.year, []).append(data)    if temp_nc:        os.remove(temp_nc)# ---------- 年度降水 ----------for year, months in sorted(year_data.items()):    # 只处理完整的12个月    if len(months) != 12:        print(f"跳过 {year} 年：只有 {len(months)} 个月")        continue    annual = np.nansum(np.stack(months), axis=0)    # 所有月份均无数据的位置仍设为无数据    all_nan = np.all(np.isnan(np.stack(months)), axis=0)    annual[all_nan] = np.nan    annual_file = os.path.join(        yearly_dir,        f"pre_{year}.tif"    )    out_data = np.where(np.isnan(annual), -9999, annual).astype("float32")    with rasterio.open(        annual_file,        "w",        driver="GTiff",        height=annual.shape[0],        width=annual.shape[1],        count=1,        dtype="float32",        crs="EPSG:4326",        transform=transform,        nodata=-9999,        compress="lzw"    ) as dst:        dst.write(out_data, 1)    print("年度：", annual_file)print("全部完成")