1. 項(xiàng)目概述與核心價(jià)值最近在做一個(gè)農(nóng)業(yè)氣象分析的項(xiàng)目需要用到長(zhǎng)時(shí)間序列的全球降水?dāng)?shù)據(jù)。CHIRPSClimate Hazards Group InfraRed Precipitation with Station data數(shù)據(jù)因其高分辨率0.05°和長(zhǎng)時(shí)序1981年至今而成為我的首選。但官網(wǎng)手動(dòng)下載幾十甚至上百個(gè)壓縮包再逐個(gè)解壓、裁剪到研究區(qū)這個(gè)流程繁瑣到讓人頭皮發(fā)麻。用Python把整個(gè)過程自動(dòng)化就成了一個(gè)非常實(shí)際且高效的需求。這不僅僅是寫個(gè)下載腳本那么簡(jiǎn)單它涉及到網(wǎng)絡(luò)請(qǐng)求的穩(wěn)定性處理、大文件的分塊下載、特定壓縮格式的解壓、以及基于地理坐標(biāo)的柵格數(shù)據(jù)批量處理是一個(gè)典型的數(shù)據(jù)工程與地理信息科學(xué)GIS結(jié)合的實(shí)戰(zhàn)場(chǎng)景。無論你是從事氣候變化研究、農(nóng)業(yè)模型構(gòu)建還是水資源評(píng)估這套自動(dòng)化流程都能幫你把寶貴的時(shí)間從重復(fù)勞動(dòng)中解放出來聚焦在真正的數(shù)據(jù)分析上。2. 技術(shù)棧選型與思路拆解面對(duì)“批量下載CHIRPS氣象數(shù)據(jù)并完成解壓裁剪”這個(gè)目標(biāo)我們需要一個(gè)清晰的技術(shù)實(shí)現(xiàn)路徑。整個(gè)流程可以拆解為四個(gè)核心環(huán)節(jié)數(shù)據(jù)發(fā)現(xiàn)與鏈接獲取、穩(wěn)定批量下載、特定格式解壓、以及空間裁剪。每個(gè)環(huán)節(jié)都有多種技術(shù)方案我的選型基于穩(wěn)定性、效率以及可維護(hù)性。2.1 核心工具庫選擇網(wǎng)絡(luò)請(qǐng)求與下載requests與tqdmrequests庫是Python中進(jìn)行HTTP通信的事實(shí)標(biāo)準(zhǔn)其API簡(jiǎn)潔且功能強(qiáng)大。對(duì)于批量下載尤其是大文件直接使用requests.get()的streamTrue參數(shù)進(jìn)行流式下載是必須的這樣可以避免將整個(gè)文件加載到內(nèi)存中。配合tqdm庫我們可以為每個(gè)下載任務(wù)添加一個(gè)美觀的進(jìn)度條實(shí)時(shí)監(jiān)控下載速度和預(yù)計(jì)剩余時(shí)間這對(duì)于處理動(dòng)輒幾百M(fèi)B的GeoTIFF文件至關(guān)重要。壓縮文件處理zipfile與gzipCHIRPS數(shù)據(jù)通常提供兩種壓縮格式.zip和.gz。對(duì)于.zip文件Python標(biāo)準(zhǔn)庫中的zipfile是唯一且最佳選擇它可以方便地解壓?jiǎn)蝹€(gè)或多個(gè)文件。對(duì)于.gz文件一種使用gzip壓縮的TAR包常見后綴為.tar.gz或.tgz我們需要結(jié)合tarfile和gzip庫。這里有一個(gè)關(guān)鍵點(diǎn)tarfile庫本身支持直接打開.tar.gz文件并自動(dòng)處理gzip解壓因此我們通常直接使用tarfile.open(moder:gz)而無需顯式調(diào)用gzip庫。柵格數(shù)據(jù)處理rasterio與geopandas這是整個(gè)流程的技術(shù)核心。rasterio是專門為讀寫柵格數(shù)據(jù)如GeoTIFF而設(shè)計(jì)的庫其API設(shè)計(jì)深受GDAL影響但更加Pythonic。我們將用它來打開下載的CHIRPS柵格文件、讀取空間參考信息、執(zhí)行裁剪操作以及寫入新的裁剪后文件。geopandas則用于處理我們的研究區(qū)矢量邊界通常是Shapefile或GeoJSON格式它可以方便地讀取矢量數(shù)據(jù)、進(jìn)行坐標(biāo)參考系統(tǒng)CRS的統(tǒng)一并將幾何對(duì)象轉(zhuǎn)換為rasterio可用的掩膜形狀。rasterio的mask函數(shù)是實(shí)現(xiàn)裁剪的關(guān)鍵。2.2 流程架構(gòu)設(shè)計(jì)整個(gè)自動(dòng)化腳本的骨架設(shè)計(jì)如下它遵循“獲取-處理-輸出”的線性流程但每個(gè)模塊內(nèi)部都有完善的錯(cuò)誤處理和日志記錄。輸入與配置用戶提供研究區(qū)的矢量邊界文件路徑、目標(biāo)時(shí)間范圍起始年-月、以及數(shù)據(jù)保存的根目錄。鏈接生成與列表構(gòu)建根據(jù)CHIRPS數(shù)據(jù)服務(wù)器的公開URL命名規(guī)則程序化生成指定時(shí)間范圍內(nèi)所有數(shù)據(jù)文件的下載鏈接列表。CHIRPS的URL通常有固定模式例如https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p05/YYYY/chirps-v2.0.YYYY.MM.DD.tif.gz穩(wěn)定批量下載遍歷鏈接列表。對(duì)每個(gè)鏈接檢查本地是否已存在該文件避免重復(fù)下載。使用requests進(jìn)行流式下載并利用tqdm顯示進(jìn)度。實(shí)現(xiàn)重試機(jī)制如使用tenacity庫或自定義while循環(huán)應(yīng)對(duì)網(wǎng)絡(luò)波動(dòng)。將文件保存到本地的臨時(shí)或原始數(shù)據(jù)目錄。智能解壓處理遍歷下載的壓縮文件。根據(jù)文件后綴.zip或.gz調(diào)用相應(yīng)的解壓函數(shù)。將解壓后的.tif文件統(tǒng)一存放到一個(gè)指定目錄如/unzipped_tifs。刪除原始的壓縮文件以節(jié)省空間可選建議在確認(rèn)解壓成功后再進(jìn)行。批量空間裁剪使用geopandas讀取研究區(qū)矢量邊界并獲取其邊界框bounds和幾何形狀geometry。確保矢量邊界的CRS與CHIRPS數(shù)據(jù)通常是WGS84EPSG:4326一致如果不一致進(jìn)行坐標(biāo)轉(zhuǎn)換。遍歷解壓后的所有.tif文件。對(duì)每個(gè)文件用rasterio打開使用rasterio.mask.mask函數(shù)以上一步獲取的幾何形狀為掩膜執(zhí)行裁剪操作。mask函數(shù)會(huì)返回裁剪后的柵格數(shù)組、變換信息以及標(biāo)簽。將裁剪后的數(shù)組寫入新的GeoTIFF文件保存到輸出目錄如/clipped_tifs并繼承原始數(shù)據(jù)的CRS等信息。日志與狀態(tài)管理在整個(gè)過程中使用Python的logging模塊記錄關(guān)鍵步驟、成功信息和錯(cuò)誤警告便于后期排查問題。注意CHIRPS數(shù)據(jù)是公開的科研數(shù)據(jù)在使用時(shí)請(qǐng)遵守其數(shù)據(jù)使用政策通常要求注明出處。批量下載時(shí)請(qǐng)務(wù)必設(shè)置合理的請(qǐng)求間隔如time.sleep(1)避免對(duì)服務(wù)器造成過大壓力體現(xiàn)良好的網(wǎng)絡(luò)公民素養(yǎng)。3. 核心模塊實(shí)現(xiàn)與代碼詳解理論講清楚了接下來我們進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)把每個(gè)模塊的代碼掰開揉碎講明白。我會(huì)提供可直接運(yùn)行的函數(shù)代碼并解釋關(guān)鍵參數(shù)和設(shè)計(jì)考量。3.1 數(shù)據(jù)鏈接生成模塊CHIRPS數(shù)據(jù)的URL結(jié)構(gòu)相對(duì)規(guī)整。我們需要根據(jù)用戶輸入的日期范圍批量生成這些鏈接。import requests from datetime import datetime, timedelta import logging logging.basicConfig(levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s) def generate_chirps_links(start_date, end_date, base_url_template): 生成指定日期范圍內(nèi)的CHIRPS數(shù)據(jù)下載鏈接列表。 參數(shù): start_date (datetime): 開始日期。 end_date (datetime): 結(jié)束日期。 base_url_template (str): URL模板使用{year}, {month:02d}, {day:02d}作為占位符。 示例: “https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p05/{year}/chirps-v2.0.{year}.{month:02d}.{day:02d}.tif.gz” 返回: list: 下載鏈接列表。 links [] current_date start_date delta timedelta(days1) while current_date end_date: url base_url_template.format( yearcurrent_date.year, monthcurrent_date.month, daycurrent_date.day ) links.append(url) current_date delta logging.debug(fGenerated link for {current_date - delta}: {url}) logging.info(fGenerated {len(links)} links from {start_date.strftime(%Y-%m-%d)} to {end_date.strftime(%Y-%m-%d)}.) return links # 使用示例 if __name__ __main__: base_template https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p05/{year}/chirps-v2.0.{year}.{month:02d}.{day:02d}.tif.gz start datetime(2023, 6, 1) end datetime(2023, 6, 10) url_list generate_chirps_links(start, end, base_template) print(url_list[:3]) # 打印前三個(gè)鏈接作為示例關(guān)鍵點(diǎn)解析日期遍歷使用datetime和timedelta可以安全、準(zhǔn)確地遍歷每一天避免手動(dòng)計(jì)算月份和閏年等問題。URL模板將URL定義為模板字符串使用str.format()方法進(jìn)行格式化代碼清晰且易于修改。{month:02d}確保月份總是兩位數(shù)字這是CHIRPS文件名所要求的。日志記錄使用logging模塊而非print可以靈活控制輸出級(jí)別如調(diào)試時(shí)用DEBUG運(yùn)行時(shí)用INFO信息也更結(jié)構(gòu)化。3.2 帶重試與進(jìn)度顯示的下載模塊這是最容易出錯(cuò)的環(huán)節(jié)。網(wǎng)絡(luò)不穩(wěn)定、服務(wù)器忙、甚至本地磁盤空間不足都可能導(dǎo)致下載失敗。一個(gè)健壯的下載函數(shù)必須包含重試機(jī)制和斷點(diǎn)續(xù)傳的考量雖然這里簡(jiǎn)化了但思路很重要。import os from pathlib import Path import requests from tqdm import tqdm import time def download_file_with_retry(url, save_path, max_retries3, chunk_size8192): 下載文件支持重試和進(jìn)度顯示。 參數(shù): url (str): 文件下載鏈接。 save_path (str or Path): 本地保存路徑。 max_retries (int): 最大重試次數(shù)。 chunk_size (int): 下載流塊大小單位字節(jié)。 返回: bool: 下載是否成功。 Path(save_path).parent.mkdir(parentsTrue, exist_okTrue) # 確保目錄存在 for attempt in range(max_retries): try: # 發(fā)起HEAD請(qǐng)求獲取文件大小用于進(jìn)度條 resp_head requests.head(url, timeout10, allow_redirectsTrue) total_size int(resp_head.headers.get(content-length, 0)) # 流式下載 response requests.get(url, streamTrue, timeout30) response.raise_for_status() # 如果狀態(tài)碼不是200拋出HTTPError # 初始化進(jìn)度條 progress_bar tqdm(totaltotal_size, unitiB, unit_scaleTrue, descPath(save_path).name, leaveFalse) with open(save_path, wb) as f: for chunk in response.iter_content(chunk_sizechunk_size): if chunk: size f.write(chunk) progress_bar.update(size) progress_bar.close() # 簡(jiǎn)單校驗(yàn)如果知道大小可以檢查本地文件大小 if total_size ! 0 and os.path.getsize(save_path) ! total_size: logging.warning(fFile size mismatch for {url}. Retrying...) os.remove(save_path) raise IOError(Downloaded file size does not match expected size.) logging.info(fSuccessfully downloaded: {save_path}) return True except (requests.exceptions.RequestException, IOError) as e: logging.warning(fAttempt {attempt 1} failed for {url}: {e}) if attempt max_retries - 1: wait_time 2 ** attempt # 指數(shù)退避策略 logging.info(fWaiting {wait_time} seconds before retry...) time.sleep(wait_time) else: logging.error(fFailed to download {url} after {max_retries} attempts.) if os.path.exists(save_path): os.remove(save_path) # 刪除不完整的文件 return False return False def batch_download(url_list, download_dir): 批量下載文件列表。 參數(shù): url_list (list): 下載鏈接列表。 download_dir (str or Path): 下載文件存儲(chǔ)目錄。 download_dir Path(download_dir) download_dir.mkdir(parentsTrue, exist_okTrue) successful_downloads [] for url in tqdm(url_list, descOverall Download Progress): # 從URL中提取文件名 filename url.split(/)[-1] save_path download_dir / filename # 檢查文件是否已存在且完整這里簡(jiǎn)化處理僅檢查存在性 if save_path.exists(): logging.info(fFile already exists, skipping: {save_path}) successful_downloads.append(save_path) continue if download_file_with_retry(url, save_path): successful_downloads.append(save_path) else: logging.error(fSkipping {url} due to download failure.) logging.info(fBatch download finished. {len(successful_downloads)}/{len(url_list)} files succeeded.) return successful_downloads實(shí)操心得chunk_size的選擇8192字節(jié)8KB是一個(gè)比較平衡的值。太小會(huì)增加循環(huán)次數(shù)和系統(tǒng)調(diào)用開銷太大會(huì)占用更多內(nèi)存。對(duì)于網(wǎng)絡(luò)環(huán)境好、文件大的情況可以嘗試增加到3276832KB或6553664KB。指數(shù)退避重試等待時(shí)間采用指數(shù)增長(zhǎng)2 ** attempt這是處理臨時(shí)性網(wǎng)絡(luò)故障的經(jīng)典策略避免在服務(wù)器恢復(fù)前進(jìn)行“狂轟濫炸”。文件存在性檢查這里只做了簡(jiǎn)單的存在性檢查。更嚴(yán)謹(jǐn)?shù)淖龇ㄊ怯涗浺严螺d文件的MD5或文件大小下次運(yùn)行時(shí)進(jìn)行比對(duì)實(shí)現(xiàn)“斷點(diǎn)續(xù)傳”和“完整性校驗(yàn)”。對(duì)于CHIRPS由于其文件命名包含日期且數(shù)據(jù)更新不頻繁簡(jiǎn)單檢查通常足夠。tqdm的leave參數(shù)設(shè)置為False可以讓單個(gè)文件的進(jìn)度條在完成后消失保持整體進(jìn)度條的整潔。如果你希望查看每個(gè)文件的具體下載情況可以設(shè)為True。3.3 多格式解壓處理模塊下載下來的文件可能是.zip或.gz我們需要一個(gè)能自動(dòng)識(shí)別并處理的解壓函數(shù)。import zipfile import tarfile import gzip import shutil from pathlib import Path def extract_file(compressed_path, extract_to_dir): 解壓文件支持 .zip 和 .gz/.tar.gz 格式。 參數(shù): compressed_path (str or Path): 壓縮文件路徑。 extract_to_dir (str or Path): 解壓目標(biāo)目錄。 返回: Path or None: 解壓出的主要文件路徑通常是.tif失敗則返回None。 compressed_path Path(compressed_path) extract_to_dir Path(extract_to_dir) extract_to_dir.mkdir(parentsTrue, exist_okTrue) extracted_file None try: if compressed_path.suffix .zip: with zipfile.ZipFile(compressed_path, r) as zip_ref: # 假設(shè)zip內(nèi)只有一個(gè)主要文件如.tif for file_info in zip_ref.infolist(): if file_info.filename.endswith(.tif): zip_ref.extract(file_info, extract_to_dir) extracted_file extract_to_dir / file_info.filename logging.info(fExtracted {file_info.filename} from ZIP.) break # 找到第一個(gè)tif文件就退出 elif compressed_path.suffix in [.gz, .tgz]: # 處理 .tar.gz 或 .tgz if .tar in compressed_path.suffixes: # 是 .tar.gz with tarfile.open(compressed_path, r:gz) as tar_ref: members tar_ref.getmembers() for member in members: if member.name.endswith(.tif): tar_ref.extract(member, extract_to_dir) extracted_file extract_to_dir / member.name logging.info(fExtracted {member.name} from TAR.GZ.) break else: # 是純 .gz (非tar包)CHIRPS的.gz通常是tar包但這里做兼容處理 logging.warning(fFile {compressed_path} is a pure .gz, which is uncommon for CHIRPS. Trying gzip decompression.) output_path extract_to_dir / compressed_path.stem # 去除 .gz 后綴 with gzip.open(compressed_path, rb) as f_in: with open(output_path, wb) as f_out: shutil.copyfileobj(f_in, f_out) extracted_file output_path logging.info(fDecompressed pure GZIP to {output_path}.) else: logging.error(fUnsupported compression format for file: {compressed_path}) return None # 驗(yàn)證解壓文件是否存在 if extracted_file and extracted_file.exists(): logging.info(fExtraction successful: {extracted_file}) # 可選刪除原壓縮包以節(jié)省空間 # compressed_path.unlink() # logging.info(fDeleted original archive: {compressed_path}) return extracted_file else: logging.error(fFailed to find extracted .tif file in {compressed_path}) return None except (zipfile.BadZipFile, tarfile.ReadError, OSError) as e: logging.error(fFailed to extract {compressed_path}: {e}) return None def batch_extract(compressed_file_list, extract_dir): 批量解壓文件。 參數(shù): compressed_file_list (list): 壓縮文件路徑列表。 extract_dir (str or Path): 解壓目標(biāo)根目錄。 extract_dir Path(extract_dir) extracted_files [] for comp_file in tqdm(compressed_file_list, descExtracting files): result extract_file(comp_file, extract_dir) if result: extracted_files.append(result) logging.info(fBatch extraction finished. {len(extracted_files)} files extracted.) return extracted_files踩過的坑壓縮包內(nèi)結(jié)構(gòu)不是所有壓縮包解壓后文件都在根目錄。有些可能包含一層以日期命名的文件夾。上述代碼假設(shè).tif文件直接在壓縮包根目錄或第一層。更健壯的做法是遞歸查找所有.tif文件或者根據(jù)CHIRPS的固定命名模式如chirps-v2.0.YYYY.MM.DD.tif來定位。.gz與.tar.gz務(wù)必區(qū)分純gzip壓縮文件.gz和先用tar打包再用gzip壓縮的文件.tar.gz或.tgz。CHIRPS數(shù)據(jù)通常是后者。用tarfile處理.tar.gz是最方便的方式。代碼中做了判斷但實(shí)際使用中應(yīng)確認(rèn)數(shù)據(jù)格式。權(quán)限與路徑在Windows系統(tǒng)上解壓出的文件可能因?yàn)閠ar包內(nèi)記錄的權(quán)限信息而導(dǎo)致權(quán)限錯(cuò)誤。通常用tar_ref.extract(member, extract_to_dir)沒問題但如果遇到問題可以嘗試使用tar_ref.extractall(extract_to_dir)并設(shè)置filterdata參數(shù)Python 3.12或手動(dòng)處理權(quán)限。3.4 基于矢量邊界的批量柵格裁剪模塊這是地理空間處理的核心。我們將使用rasterio的mask函數(shù)它需要一個(gè)幾何圖形列表作為輸入。import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np from pathlib import Path def clip_raster_with_shapefile(raster_path, shapefile_path, output_path, cropTrue, all_touchedFalse): 使用矢量面文件裁剪柵格。 參數(shù): raster_path (str or Path): 輸入柵格文件路徑。 shapefile_path (str or Path): 裁剪用的矢量面文件路徑支持.shp或.geojson等。 output_path (str or Path): 輸出裁剪后柵格文件路徑。 crop (bool): 是否將輸出柵格的范圍裁剪到掩膜幾何圖形的最小外接矩形。默認(rèn)為True可以減小文件體積。 all_touched (bool): 是否裁剪所有被幾何圖形接觸到的像素。默認(rèn)為False只裁剪中心點(diǎn)在圖形內(nèi)的像素。 對(duì)于低分辨率數(shù)據(jù)或需要更精確邊界時(shí)可設(shè)為True。 返回: bool: 裁剪是否成功。 try: # 1. 讀取矢量邊界 gdf gpd.read_file(shapefile_path) # 確保矢量數(shù)據(jù)是地理坐標(biāo)系WGS84與CHIRPS數(shù)據(jù)一致 if gdf.crs is None: logging.warning(fShapefile {shapefile_path} has no CRS. Assuming WGS84 (EPSG:4326).) gdf.set_crs(EPSG:4326, inplaceTrue) elif gdf.crs.to_epsg() ! 4326: logging.info(fReprojecting shapefile from {gdf.crs} to WGS84 (EPSG:4326).) gdf gdf.to_crs(EPSG:4326) # 將所有幾何圖形合并為一個(gè)如果有多邊形 if len(gdf) 1: geometry gdf.unary_union else: geometry gdf.geometry.iloc[0] # mask函數(shù)需要geojson-like的字典列表 geoms [geometry.__geo_interface__] # 2. 打開柵格并執(zhí)行裁剪 with rasterio.open(raster_path) as src: # 執(zhí)行裁剪操作 out_image, out_transform mask(src, geoms, cropcrop, all_touchedall_touched, nodatasrc.nodata) # 獲取裁剪后的元數(shù)據(jù) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: src.nodata }) # 3. 寫入輸出文件 with rasterio.open(output_path, w, **out_meta) as dest: dest.write(out_image) logging.info(fSuccessfully clipped: {output_path}) return True except Exception as e: logging.error(fFailed to clip {raster_path} with {shapefile_path}: {e}) return False def batch_clip(raster_file_list, shapefile_path, output_dir, clip_kwargsNone): 批量裁剪柵格文件。 參數(shù): raster_file_list (list): 待裁剪的柵格文件路徑列表。 shapefile_path (str or Path): 矢量邊界文件路徑。 output_dir (str or Path): 輸出目錄。 clip_kwargs (dict, optional): 傳遞給clip_raster_with_shapefile函數(shù)的其他參數(shù)。 返回: list: 成功裁剪的輸出文件路徑列表。 if clip_kwargs is None: clip_kwargs {} output_dir Path(output_dir) output_dir.mkdir(parentsTrue, exist_okTrue) clipped_files [] for raster_file in tqdm(raster_file_list, descClipping rasters): raster_path Path(raster_file) # 構(gòu)建輸出文件名例如chirps-v2.0.2023.06.01_clipped.tif output_filename raster_path.stem _clipped.tif output_path output_dir / output_filename if output_path.exists(): logging.info(fClipped file already exists, skipping: {output_path}) clipped_files.append(output_path) continue if clip_raster_with_shapefile(raster_file, shapefile_path, output_path, **clip_kwargs): clipped_files.append(output_path) logging.info(fBatch clipping finished. {len(clipped_files)}/{len(raster_file_list)} files processed.) return clipped_files關(guān)鍵技術(shù)細(xì)節(jié)CRS一致性這是空間分析中最常見的錯(cuò)誤來源。CHIRPS數(shù)據(jù)的CRS是WGS84地理坐標(biāo)系EPSG:4326。你的研究區(qū)矢量邊界也必須轉(zhuǎn)換或確認(rèn)為此坐標(biāo)系否則裁剪會(huì)錯(cuò)位甚至失敗。geopandas的to_crs方法可以方便地進(jìn)行轉(zhuǎn)換。crop參數(shù)設(shè)置為True時(shí)rasterio.mask.mask函數(shù)不僅將圖形外的像素設(shè)為nodata還會(huì)將輸出圖像的范圍邊界框縮小到掩膜幾何圖形的最小外接矩形。這能顯著減小輸出文件的大小通常是推薦的做法。all_touched參數(shù)這個(gè)參數(shù)決定了像素被保留或丟棄的規(guī)則。False默認(rèn)使用“中心點(diǎn)規(guī)則”即像素中心點(diǎn)在圖形內(nèi)才被保留。True使用“全部接觸規(guī)則”只要像素被圖形接觸到即使中心點(diǎn)在外面就被保留。對(duì)于低分辨率數(shù)據(jù)或需要精確邊界時(shí)如小流域設(shè)為True可能更合適但會(huì)導(dǎo)致邊界更“粗糙”。nodata值裁剪后圖形區(qū)域外的像素會(huì)被賦予nodata值。我們繼承了原始柵格的nodata值通常是-9999或nan確保數(shù)據(jù)的一致性。在后續(xù)分析中需要使用如numpy.nan或rasterio的掩膜數(shù)組來處理這些無效值。4. 完整流程串聯(lián)與配置管理將上述模塊組合起來并添加一些工程化的配置管理就構(gòu)成了一個(gè)完整的自動(dòng)化腳本。我們使用一個(gè)配置文件如config.yaml或config.json來管理所有路徑和參數(shù)這樣更清晰也便于復(fù)用和分享。config.yaml示例# config.yaml project: name: CHIRPS_Data_Pipeline paths: shapefile: ./data/boundary/study_area.shp # 研究區(qū)矢量邊界 download_dir: ./data/raw_downloads # 原始?jí)嚎s文件下載目錄 extract_dir: ./data/unzipped_tifs # 解壓后的TIFF文件目錄 output_dir: ./data/clipped_tifs # 裁剪后的TIFF文件輸出目錄 log_file: ./logs/pipeline.log # 日志文件路徑 download: base_url_template: https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p05/{year}/chirps-v2.0.{year}.{month:02d}.{day:02d}.tif.gz start_date: 2023-06-01 end_date: 2023-06-30 max_retries: 5 chunk_size_kb: 8192 # 單位KB clip: crop: true all_touched: false主程序main.py# main.py import yaml from datetime import datetime from pathlib import Path import logging from modules.download import generate_chirps_links, batch_download from modules.extract import batch_extract from modules.clip import batch_clip def setup_logging(log_file): 配置日志記錄到文件和控制臺(tái) Path(log_file).parent.mkdir(parentsTrue, exist_okTrue) logging.basicConfig( levellogging.INFO, format%(asctime)s - %(name)s - %(levelname)s - %(message)s, handlers[ logging.FileHandler(log_file, encodingutf-8), logging.StreamHandler() ] ) def main(config_pathconfig.yaml): # 1. 加載配置 with open(config_path, r, encodingutf-8) as f: config yaml.safe_load(f) paths config[paths] dl_cfg config[download] clip_cfg config[clip] # 2. 設(shè)置日志 setup_logging(paths[log_file]) logger logging.getLogger(__name__) logger.info(*50) logger.info(fStarting CHIRPS Data Processing Pipeline: {config[project][name]}) logger.info(*50) # 3. 生成下載鏈接 start_date datetime.strptime(dl_cfg[start_date], %Y-%m-%d) end_date datetime.strptime(dl_cfg[end_date], %Y-%m-%d) url_list generate_chirps_links(start_date, end_date, dl_cfg[base_url_template]) if not url_list: logger.error(No download links generated. Check date range and URL template.) return # 4. 批量下載 logger.info(Step 1: Batch Downloading...) downloaded_files batch_download(url_list, paths[download_dir]) if not downloaded_files: logger.error(Download failed or no files downloaded. Exiting.) return # 5. 批量解壓 logger.info(Step 2: Batch Extracting...) extracted_files batch_extract(downloaded_files, paths[extract_dir]) if not extracted_files: logger.error(Extraction failed or no files extracted. Exiting.) return # 6. 批量裁剪 logger.info(Step 3: Batch Clipping...) clipped_files batch_clip( extracted_files, paths[shapefile], paths[output_dir], clip_kwargs{crop: clip_cfg[crop], all_touched: clip_cfg[all_touched]} ) # 7. 完成 logger.info(*50) logger.info(Pipeline finished.) logger.info(f Downloaded: {len(downloaded_files)} files) logger.info(f Extracted: {len(extracted_files)} files) logger.info(f Clipped: {len(clipped_files)} files) logger.info(f Output directory: {paths[output_dir]}) logger.info(*50) if __name__ __main__: main()這個(gè)主程序?qū)⒏鱾€(gè)模塊串聯(lián)起來并通過配置文件管理所有參數(shù)使得整個(gè)流程清晰、可配置、易于維護(hù)。你可以通過修改config.yaml文件來處理不同的時(shí)間范圍、研究區(qū)域和輸出設(shè)置。5. 常見問題與排查技巧實(shí)錄在實(shí)際運(yùn)行中你幾乎一定會(huì)遇到各種問題。下面是我在多次運(yùn)行類似腳本后總結(jié)的“排坑指南”。5.1 網(wǎng)絡(luò)與下載問題問題1下載速度極慢或頻繁中斷。排查首先檢查網(wǎng)絡(luò)連接。其次CHIRPS數(shù)據(jù)服務(wù)器位于國(guó)外國(guó)內(nèi)直接訪問可能不穩(wěn)定。解決增加重試次數(shù)和等待時(shí)間將max_retries提高到5或10并適當(dāng)增加指數(shù)退避的基數(shù)。使用會(huì)話Sessionrequests.Session()可以復(fù)用TCP連接對(duì)批量下載同一服務(wù)器的文件有性能提升。考慮代理如果條件允許且合規(guī)配置網(wǎng)絡(luò)代理可能改善連接穩(wěn)定性。注意此處僅作技術(shù)可能性探討具體實(shí)施需符合當(dāng)?shù)胤煞ㄒ?guī)和網(wǎng)絡(luò)使用政策。分時(shí)段運(yùn)行在網(wǎng)絡(luò)相對(duì)空閑的時(shí)段如凌晨運(yùn)行腳本。問題2HTTP 403 Forbidden 或 404 Not Found 錯(cuò)誤。排查URL鏈接錯(cuò)誤或服務(wù)器上文件不存在/權(quán)限變更。解決手動(dòng)驗(yàn)證URL在瀏覽器中打開一個(gè)生成的鏈接看是否能正常訪問或下載。檢查日期范圍確認(rèn)你請(qǐng)求的日期在CHIRPS數(shù)據(jù)發(fā)布范圍內(nèi)。CHIRPS數(shù)據(jù)通常有幾天到一周的延遲。檢查URL模板CHIRPS的URL結(jié)構(gòu)偶爾會(huì)微調(diào)。去官網(wǎng)查看最新數(shù)據(jù)的下載鏈接核對(duì)模板。5.2 解壓與文件處理問題問題3zipfile.BadZipFile或tarfile.ReadError錯(cuò)誤。排查下載的文件不完整或已損壞。解決刪除損壞文件并重新下載在download_file_with_retry函數(shù)中我們已經(jīng)在最終失敗后刪除了不完整文件。你可以手動(dòng)刪除報(bào)錯(cuò)的那個(gè)壓縮包重新運(yùn)行腳本。增加下載完整性校驗(yàn)更嚴(yán)格的做法是在下載函數(shù)中計(jì)算下載文件的MD5或SHA256哈希值與服務(wù)器提供的如果有或已知的正確值進(jìn)行比對(duì)。問題4解壓后找不到.tif文件。排查壓縮包內(nèi)的文件結(jié)構(gòu)可能與預(yù)期不符或者解壓到了子文件夾里。解決修改extract_file函數(shù)將只尋找第一個(gè).tif文件的邏輯改為遞歸搜索整個(gè)解壓目錄并返回找到的所有.tif文件列表。手動(dòng)檢查用解壓軟件手動(dòng)打開一個(gè)出問題的壓縮包觀察其內(nèi)部結(jié)構(gòu)然后調(diào)整代碼中的提取路徑邏輯。5.3 空間裁剪與GIS相關(guān)問題問題5裁剪后的圖像全是nodata值或者范圍完全不對(duì)。排查幾乎肯定是坐標(biāo)參考系統(tǒng)CRS不匹配。解決打印CRS信息在裁剪函數(shù)開頭添加代碼打印輸入柵格和矢量數(shù)據(jù)的CRS。with rasterio.open(raster_path) as src: print(fRaster CRS: {src.crs}) print(fVector CRS: {gdf.crs})強(qiáng)制統(tǒng)一CRS確保在裁剪前矢量數(shù)據(jù)通過gdf.to_crs(EPSG:4326)轉(zhuǎn)換到了WGS84EPSG:4326。如果柵格不是4326CHIRPS是則需要先將矢量轉(zhuǎn)換到柵格的CRS。檢查幾何圖形有效性有時(shí)矢量數(shù)據(jù)可能存在自相交等無效幾何圖形?？梢杂胓df.geometry.is_valid檢查并用gdf.geometry.buffer(0)進(jìn)行簡(jiǎn)單修復(fù)。問題6裁剪過程內(nèi)存占用過高甚至導(dǎo)致程序崩潰。排查研究區(qū)過大或原始CHIRPS數(shù)據(jù)分辨率高全球0.05°一次性讀入內(nèi)存的數(shù)組很大。解決使用rasterio.windows對(duì)于非常大的裁剪任務(wù)可以考慮使用rasterio.windows.Window進(jìn)行分塊讀取和寫入但這會(huì)顯著增加代碼復(fù)雜度。升級(jí)硬件增加系統(tǒng)內(nèi)存是最直接的方法?？紤]使用專業(yè)GIS軟件對(duì)于超大規(guī)模的批量裁剪像GDAL命令行工具gdalwarp或gdal_translate配合-cutline可能更內(nèi)存高效。你可以用Python的subprocess模塊來調(diào)用這些命令。問題7輸出文件體積異常大。排查數(shù)據(jù)類型和壓縮選項(xiàng)。解決在clip_raster_with_shapefile函數(shù)的out_meta更新部分添加壓縮選項(xiàng)。out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: src.nodata, compress: lzw, # 或 deflate, packbits predictor: 2, # 對(duì)于浮點(diǎn)型數(shù)據(jù)配合lzw或deflate壓縮效果更好 tiled: True, # 創(chuàng)建分塊存儲(chǔ)有利于大數(shù)據(jù)讀取 blockxsize: 256, blockysize: 256 })lzw是一種無損壓縮通常能減少30%-70%的文件大小且大多數(shù)GIS軟件都支持讀取。5.4 性能優(yōu)化建議當(dāng)處理長(zhǎng)時(shí)間序列如一整年甚至十年的數(shù)據(jù)時(shí)效率變得很重要。并行下載可以使用concurrent.futures.ThreadPoolExecutor來并發(fā)下載多個(gè)文件。但務(wù)必謹(jǐn)慎對(duì)同一服務(wù)器發(fā)起過多并發(fā)連接是不禮貌的可能導(dǎo)致你的IP被暫時(shí)封鎖。建議將最大并發(fā)數(shù)限制在3-5個(gè)。并行裁剪裁剪是CPU密集型任務(wù)適合使用multiprocessing.Pool進(jìn)行多進(jìn)程并行。每個(gè)進(jìn)程處理一個(gè)柵格文件可以充分利用多核CPU。注意將矢量數(shù)據(jù)讀入每個(gè)子進(jìn)程。緩存矢量數(shù)據(jù)在批量裁剪函數(shù)batch_clip中每次循環(huán)都讀取一次矢量邊界文件是低效的。應(yīng)該在循環(huán)外讀取一次然后將幾何對(duì)象傳遞給每個(gè)裁剪任務(wù)。最后記得良好的日志記錄是你的最佳幫手。確保日志級(jí)別設(shè)置為INFO并輸出到文件這樣即使程序在后臺(tái)運(yùn)行你也能隨時(shí)查看進(jìn)度和定位錯(cuò)誤。這套自動(dòng)化流程一旦跑通以后任何需要CHIRPS數(shù)據(jù)的項(xiàng)目你只需要修改配置文件中的日期和邊界然后泡杯咖啡等待結(jié)果即可。