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("全部完成")
10个月宝宝每天需要喝多少奶粉？