1. 项目缘起与核心价值最近在做一个区域性的水文模型项目需要用到长时间序列的降水数据。全球公开的气象数据源不少但综合考虑时空分辨率、数据连续性、易获取性和权威性CHIRPSClimate Hazards Group InfraRed Precipitation with Station data数据成了我的首选。它提供了1981年至今、0.05°分辨率约5.6公里的准全球日/月降水量数据对于区域尺度的研究非常友好。然而真到动手下载和处理时问题就来了。我需要的是中国西南某流域过去20年的日尺度数据。这意味着要面对几十甚至上百个.tif.gz压缩文件。手动去官网一个个点开、下载、解压、再裁剪到我的研究区光是想想这个流程就足以让人放弃。这不仅是效率问题更关乎可重复性——下次换一个区域或者更新数据难道还要再来一遍所以这个项目的核心价值就非常明确了用Python实现CHIRPS数据的全自动化、批量化处理流水线。从根据需求自动生成下载列表到稳定下载、校验、解压再到利用GIS库进行精确的空间裁剪最终输出可以直接用于分析的、规整的GeoTIFF文件。整个过程脚本化一次编写终身受益。这对于从事气候、水文、生态等领域研究需要处理大量栅格数据的朋友来说是一个能极大解放生产力的实用技能。2. CHIRPS数据源解析与自动化下载策略CHIRPS数据由美国加州大学圣巴巴拉分校的气候灾害中心CHC制作并维护。其官网提供了FTP和HTTP两种访问方式。对于自动化脚本我们通常选择结构更清晰、更稳定的HTTP目录。2.1 数据目录结构与命名规则理解文件命名规则是自动化的第一步。CHIRPS日数据的典型路径和文件名如下https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p06/2022/chirps-v2.0.2022.01.01.tif.gz我们来拆解一下global_daily/tifs/p06/: 表示这是全球日数据TIFF格式降水量单位是毫米/天p06。2022/: 年份目录。chirps-v2.0.2022.01.01.tif.gz: 具体的文件名。其格式为chirps-v2.0.YYYY.MM.DD.tif.gz。v2.0是版本号YYYY.MM.DD对应年、月、日。月数据、年数据以及不同空间聚合层级的数据目录结构类似只是路径和文件名中的时间精度不同。例如月数据文件可能命名为chirps-v2.0.2022.01.tif.gz。2.2 使用requests库构建稳健的下载器手动点击下载不可取我们需要用程序模拟这一行为。Python的requests库是完成HTTP请求的利器。但直接用一个简单的get下载大文件一个日数据压缩包约1-2MB会遇到几个问题网络不稳定导致中断、服务器限制、以及如何避免重复下载。首先我们需要构造目标文件的URL。这需要根据用户输入的起始结束日期动态生成。import requests from datetime import datetime, timedelta def generate_daily_urls(start_date, end_date): 生成指定日期范围内的CHIRPS日数据URL列表 base_url https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p06 urls [] current_date start_date while current_date end_date: year current_date.strftime(%Y) filename fchirps-v2.0.{current_date.strftime(%Y.%m.%d)}.tif.gz url f{base_url}/{year}/{filename} urls.append((url, filename)) # 同时返回URL和期望的文件名 current_date timedelta(days1) return urls # 示例下载2022年1月的数据 start datetime(2022, 1, 1) end datetime(2022, 1, 31) url_file_pairs generate_daily_urls(start, end)接下来是下载函数。这里有几个关键点设置请求头有些服务器会检查User-Agent模拟一个浏览器的请求头可以避免被简单的反爬机制拦截。流式下载使用streamTrue参数这样requests不会立即将整个文件加载到内存而是允许我们迭代内容这对于大文件至关重要。分块写入与进度显示将数据流分块写入本地文件同时可以计算并显示下载进度对于批量下载时的心理安慰和故障排查很有帮助。异常处理与重试网络请求充满不确定性。必须用try-except包裹并实现简单的重试逻辑。文件完整性校验下载完成后比较本地文件大小和服务器返回的Content-Length头信息进行初步校验。def download_file(url, filename, save_dir, max_retries3): 下载单个文件支持重试和进度显示 filepath os.path.join(save_dir, filename) # 如果文件已存在且大小正常可跳过下载简易去重 if os.path.exists(filepath): local_size os.path.getsize(filepath) try: head requests.head(url, timeout10) remote_size int(head.headers.get(content-length, 0)) if local_size remote_size and remote_size 0: print(f文件已存在且完整跳过下载: {filename}) return True except: pass # 如果HEAD请求失败则尝试重新下载 headers {User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64) ...} for attempt in range(max_retries): try: print(f开始下载: {filename} (尝试 {attempt 1}/{max_retries})) with requests.get(url, headersheaders, streamTrue, timeout30) as r: r.raise_for_status() # 检查HTTP请求是否成功 total_size int(r.headers.get(content-length, 0)) downloaded_size 0 with open(filepath, wb) as f: for chunk in r.iter_content(chunk_size8192): if chunk: f.write(chunk) downloaded_size len(chunk) # 显示进度条 if total_size: percent (downloaded_size / total_size) * 100 print(f\r进度: {percent:.1f}%, end) print() # 换行 # 最终校验 if total_size and os.path.getsize(filepath) ! total_size: raise IOError(f文件大小不匹配: {filename}) print(f下载成功: {filename}) return True except (requests.exceptions.RequestException, IOError) as e: print(f下载失败: {e}) if os.path.exists(filepath): os.remove(filepath) # 删除不完整的文件 if attempt max_retries - 1: print(f达到最大重试次数放弃: {filename}) return False time.sleep(2 ** attempt) # 指数退避等待 return False注意频繁、大量地请求服务器可能对其造成压力。在实际应用中建议在循环下载每个文件之间添加一个短暂的随机延时例如time.sleep(random.uniform(0.5, 1.5))以示友好也降低自身IP被临时限制的风险。3. 高效解压与格式处理.gz与.tif下载得到的是.tif.gz文件这是经过Gzip压缩的TIFF栅格数据。我们需要两步操作解压得到.tif然后才能进行后续的空间处理。3.1 使用gzip标准库进行解压Python标准库中的gzip模块可以直接处理.gz文件无需调用外部命令。这保证了脚本的跨平台性Windows/macOS/Linux均可运行。import gzip import shutil def decompress_gz(gz_path, output_dir): 解压.gz文件到指定目录 filename os.path.basename(gz_path) # 去除.gz后缀得到原始的.tif文件名 tif_filename filename[:-3] if filename.endswith(.gz) else filename tif_path os.path.join(output_dir, tif_filename) try: with gzip.open(gz_path, rb) as f_in: with open(tif_path, wb) as f_out: shutil.copyfileobj(f_in, f_out) print(f解压成功: {tif_filename}) # 可选解压后删除原.gz文件以节省空间 # os.remove(gz_path) return tif_path except Exception as e: print(f解压失败 {gz_path}: {e}) return None这里有一个实操心得对于几百个文件逐个解压是串行操作可能会比较慢。如果追求极致效率可以考虑使用Python的concurrent.futures模块实现多线程/多进程并行解压特别是当解压是CPU密集型操作时。但对于IO密集型磁盘读写且单个文件不大的情况串行和并行的差别可能不大甚至并行可能因磁盘争用而变慢。我的建议是先实现串行版本确保稳定如果确实成为瓶颈再考虑并行优化。3.2 理解GeoTIFF不仅仅是图片解压后我们得到.tif文件但它不是普通的图片。它是GeoTIFF一种包含地理参考信息的栅格数据格式。简单来说一个GeoTIFF文件包含两部分图像数据一个二维数组每个像素像元的值代表该位置的降水量单位mm。地理信息定义了数组如何映射到真实地球表面。这包括投影CRS例如WGS84地理坐标系EPSG:4326。仿射变换参数定义了像素坐标行、列如何转换为地理坐标经度、纬度。通常包含左上角坐标、像素宽度和高度等信息。NoData值用于标记无效数据或海洋区域的特殊值CHIRPS中常用-9999。这些信息都被编码在TIFF文件的标签Tags里。我们需要专门的库来正确读取这些信息并进行空间运算。4. 基于GDAL/Rasterio的精确空间裁剪裁剪就是从全球数据中“切”出我们感兴趣区域AOI的那一部分。这需要用到GIS库。在Python生态中rasterio基于GDAL是处理栅格数据的首选其API设计非常Pythonic。4.1 准备工作定义裁剪区域裁剪区域通常用一个矢量面文件来定义如Shapefile.shp或GeoJSON。这里假设我们有一个研究区的矢量边界文件basin_boundary.shp。首先我们需要用geopandas处理矢量和rasterio处理栅格来读取这个边界并获取其几何信息和坐标参考系。import geopandas as gpd import rasterio from rasterio.mask import mask import numpy as np def get_features_from_shapefile(shp_path): 从Shapefile中读取几何特征并转换为rasterio.mask所需的格式 gdf gpd.read_file(shp_path) # 确保所有几何图形都是多边形并且坐标参考系已知 if gdf.crs is None: # 如果矢量文件没有CRS需要根据实际情况指定这里假设是WGS84 gdf gdf.set_crs(EPSG:4326, allow_overrideTrue) # rasterio.mask需要geojson-like的字典列表 features [json.loads(gdf.to_json())[features][0][geometry]] # 更稳健的写法处理多个多边形的情况 # features [geom for geom in gdf.geometry] # 但rasterio.mask的geometries参数可以直接接受GeoDataFrame的geometry迭代器 return gdf.geometry, gdf.crs4.2 执行裁剪操作核心步骤是使用rasterio.mask.mask函数。这个函数会计算栅格数据与裁剪区域的重叠部分。根据边界对栅格数据进行切割。可选地进行重投影如果栅格和矢量的CRS不同。将所有边界外的像元值设置为nodata。def clip_raster_with_vector(tif_path, vector_geom, vector_crs, output_path): 使用矢量边界裁剪TIFF栅格 try: with rasterio.open(tif_path) as src: # 检查CRS是否一致如果不一致需要重投影矢量几何 if src.crs ! vector_crs: # 将矢量几何重投影到栅格的CRS vector_geom_proj vector_geom.to_crs(src.crs) geoms [geom for geom in vector_geom_proj] else: geoms [geom for geom in vector_geom] # 执行裁剪 # cropTrue会调整输出图像的范围使其恰好包含裁剪区域节省空间。 # all_touchedFalse意味着只有中心点在多边形内的像元会被保留这是默认且通常更精确的做法。 out_image, out_transform mask(src, geoms, cropTrue, all_touchedFalse) # 获取原数据的元数据profile并更新以反映裁剪后的变化 out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], # 注意numpy数组维度是波段高宽 width: out_image.shape[2], transform: out_transform, crs: src.crs }) # 写入新的裁剪后的TIFF文件 with rasterio.open(output_path, w, **out_meta) as dest: dest.write(out_image) print(f裁剪成功: {os.path.basename(output_path)}) return True except Exception as e: print(f裁剪失败 {tif_path}: {e}) return False关键参数解析all_touched: 这个参数非常重要。如果设置为True那么任何被多边形接触到的像元都会被保留。这会导致边界“膨胀”特别是对于低分辨率数据。对于像CHIRPS这种约5.6km分辨率的降水数据设置为False只保留中心点在多边形内的像元通常能获得更精确的边界但可能会损失最边缘的一点点数据。选择哪种取决于你的研究对边界精度的要求。我通常先设为False如果发现裁剪后的区域比预期小太多再考虑设为True或对矢量做一点缓冲buffer。4.3 处理过程中的常见陷阱与优化内存管理全球日数据TIFF解压后大约7-8MB一次性全部读入内存进行批量裁剪问题不大。但如果你处理的是更高分辨率或更大范围的数据需要注意内存使用。rasterio的窗口读取windowed reading/writing功能可以处理超大型栅格。NoData值一致性裁剪后边界外的区域被赋予了NoData值。你需要确保输出文件的NoData值与原始数据一致通常是-9999并在后续分析如求区域平均降水量时正确处理这些值。rasterio在mask操作时会自动处理但写入时需在out_meta中明确指定nodata: src.nodata。坐标参考系CRS匹配这是最易出错的一环。务必确保你的矢量边界文件的CRS与CHIRPS数据通常是WGS84EPSG:4326一致或者在裁剪时进行正确的重投影。上面的代码演示了如何动态处理CRS不一致的情况。输出文件组织批量处理会产生大量文件。良好的文件组织至关重要。建议按日期或处理阶段建立子文件夹例如./downloaded/,./decompressed/,./clipped/。5. 构建完整的自动化处理流水线现在我们将下载、解压、裁剪三个模块串联起来形成一个完整的流水线脚本。这个脚本应该足够灵活允许用户指定时间范围、研究区矢量文件、以及输出目录。import os import time import json from datetime import datetime, timedelta import requests import gzip import shutil import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np # 配置参数 START_DATE datetime(2022, 6, 1) END_DATE datetime(2022, 8, 31) SHAPEFILE_PATH ./data/basin_boundary.shp OUTPUT_BASE_DIR ./processed_data # 创建子目录 DOWNLOAD_DIR os.path.join(OUTPUT_BASE_DIR, downloaded) DECOMPRESS_DIR os.path.join(OUTPUT_BASE_DIR, decompressed) CLIP_DIR os.path.join(OUTPUT_BASE_DIR, clipped) for d in [DOWNLOAD_DIR, DECOMPRESS_DIR, CLIP_DIR]: os.makedirs(d, exist_okTrue) # 步骤1: 生成下载列表并下载 print( 步骤1: 生成下载列表 ) url_file_pairs generate_daily_urls(START_DATE, END_DATE) print(f共发现 {len(url_file_pairs)} 个文件待处理。) success_download [] for url, filename in url_file_pairs: if download_file(url, filename, DOWNLOAD_DIR): success_download.append(filename) print(f下载完成。成功: {len(success_download)} / {len(url_file_pairs)}) # 步骤2: 解压所有下载的.gz文件 print(\n 步骤2: 解压文件 ) tif_paths [] for filename in success_download: gz_path os.path.join(DOWNLOAD_DIR, filename) tif_path decompress_gz(gz_path, DECOMPRESS_DIR) if tif_path: tif_paths.append(tif_path) # 步骤3: 读取矢量边界 print(\n 步骤3: 加载裁剪边界 ) try: aoi_gdf gpd.read_file(SHAPEFILE_PATH) if aoi_gdf.crs is None: aoi_gdf aoi_gdf.set_crs(EPSG:4326, allow_overrideTrue) print(f边界文件加载成功CRS: {aoi_gdf.crs}) except Exception as e: print(f加载边界文件失败: {e}) exit(1) # 步骤4: 批量裁剪 print(\n 步骤4: 批量裁剪 ) success_clip 0 for tif_path in tif_paths: tif_filename os.path.basename(tif_path) output_filename tif_filename.replace(.tif, _clipped.tif) # 给裁剪后的文件加后缀 output_path os.path.join(CLIP_DIR, output_filename) if clip_raster_with_vector(tif_path, aoi_gdf.geometry, aoi_gdf.crs, output_path): success_clip 1 print(f\n 处理总结 ) print(f计划处理天数: {(END_DATE - START_DATE).days 1}) print(f成功下载: {len(success_download)}) print(f成功解压: {len(tif_paths)}) print(f成功裁剪: {success_clip}) print(f最终裁剪数据保存在: {CLIP_DIR})这个脚本提供了一个基础框架。在实际项目中你可能还需要增加更多功能比如日志记录将运行信息成功、失败、错误原因写入日志文件便于后期排查。配置文件将时间范围、路径等参数外置到JSON或YAML配置文件避免硬编码。断点续传记录已成功处理的文件列表下次运行时跳过它们。更复杂的日期处理支持下载月数据、年数据或非连续日期。并行处理使用concurrent.futures或multiprocessing库并行执行下载或裁剪任务大幅提升处理速度。6. 环境配置与依赖管理为了让这个脚本在任何机器上都能顺利运行我们需要管理好Python环境。强烈推荐使用conda或venv创建独立的虚拟环境并使用requirements.txt文件记录依赖。创建一个requirements.txt文件内容如下requests2.28.0 rasterio1.3.0 geopandas0.12.0 numpy1.21.0geopandas的安装稍微复杂因为它依赖一些C库如GDAL, GEOS, PROJ。在Windows上最简单的方法是通过conda安装因为它能自动处理这些二进制依赖。# 使用conda推荐尤其对于Windows用户 conda create -n chirps_processor python3.9 conda activate chirps_processor conda install -c conda-forge requests rasterio geopandas # 或者使用pip确保系统已安装GDAL等开发库 pip install -r requirements.txt对于Linux/macOS用户可能需要先通过系统包管理器安装GDAL等库然后再用pip安装Python包。一个我踩过的坑不同版本的rasterio和GDAL库有时存在兼容性问题。如果遇到奇怪的错误比如无法打开某些TIFF文件或投影信息读取错误可以尝试固定版本或者使用conda-forge频道安装那里的二进制兼容性通常维护得更好。例如conda install -c conda-forge rasterio1.3.4 gdal3.6.0。7. 从数据到洞察后续分析示例得到裁剪好的每日降水TIFF序列后真正的分析才刚刚开始。这里举两个简单的例子展示如何用rasterio和xarray进行后续处理。示例1计算研究区面平均日降水量import rasterio import numpy as np def calculate_basin_average_precipitation(tif_path): 读取裁剪后的TIFF计算有效像元非NoData的平均值 with rasterio.open(tif_path) as src: data src.read(1) # 读取第一个波段 nodata src.nodata if nodata is not None: # 创建掩膜排除NoData值 valid_mask data ! nodata valid_data data[valid_mask] if valid_data.size 0: avg_precip np.mean(valid_data) return avg_precip else: return None # 整个区域都是NoData理论上裁剪后不应发生 else: # 如果没有定义NoData计算整个数组的平均不推荐 return np.mean(data) # 遍历裁剪后的文件夹计算每日平均降水量 daily_avg {} for clipped_file in os.listdir(CLIP_DIR): if clipped_file.endswith(.tif): date_str clipped_file.split(_)[0] # 假设文件名以日期开头 filepath os.path.join(CLIP_DIR, clipped_file) avg calculate_basin_average_precipitation(filepath) if avg is not None: daily_avg[date_str] round(avg, 2) print(daily_avg)示例2使用xarray进行时间序列分析与可视化xarray库非常适合处理多维栅格数据如时间序列。我们可以将多日的TIFF文件堆叠成一个数据立方体。import xarray as xr import rioxarray # 使xarray能直接读取GeoTIFF import matplotlib.pyplot as plt # 读取所有裁剪后的文件并沿时间维度合并 file_pattern os.path.join(CLIP_DIR, *.tif) ds xr.open_mfdataset(file_pattern, enginerasterio, combineby_coords, chunks{time: 10}) # chunks参数用于启用Dask并行计算处理大数据时非常有用 # 计算整个区域空间维度的平均时间序列 ts_mean ds[band_data].mean(dim[x, y]) # 假设数据变量名为band_data需根据实际情况调整 # 绘制时间序列图 plt.figure(figsize(12, 4)) ts_mean.plot() plt.title(研究区平均日降水量时间序列 (2022年夏季)) plt.ylabel(降水量 (mm)) plt.grid(True, alpha0.3) plt.tight_layout() plt.savefig(./precipitation_timeseries.png, dpi150) plt.show() # 也可以计算月总量或进行其他统计分析通过这样一个从数据获取、预处理到初步分析的完整自动化流程我们不仅节省了大量重复劳动的时间更重要的是建立了一个可重复、可追溯、可扩展的数据处理框架。无论是更换研究区、延长研究时段还是接入其他类似的气象数据源如ERA5, GPM这套方法论和代码骨架都能提供有力的支持。