據(jù)：從HDF到NDVI/ET成品圖全流程)
遙感數(shù)據(jù)處理這件事最怕的不是算法難而是重復(fù)勞動(dòng)。我最早接觸MODIS數(shù)據(jù)的時(shí)候一個(gè)研究區(qū)、一年12個(gè)月、4種產(chǎn)品光是HDF文件的拼接、投影、裁剪、格式轉(zhuǎn)換就能耗掉整整兩天中間還得盯著MRT那個(gè)老界面反復(fù)點(diǎn)選參數(shù)稍不留神就報(bào)錯(cuò)重來。后來被逼得沒辦法硬是把整套流程用Python腳本串了起來配合ArcGIS做后處理和出圖現(xiàn)在同樣的工作量壓縮到半小時(shí)以內(nèi)而且可以批量跑、無人值守。這篇內(nèi)容就是把這套流程完整拆開講清楚從HDF原始文件到NDVI、ET成品圖每一步為什么這么做、參數(shù)怎么定、坑在哪里都會(huì)講到。適合做生態(tài)遙感、農(nóng)業(yè)監(jiān)測(cè)、水文分析方向的同學(xué)也適合任何需要批量處理柵格數(shù)據(jù)的GIS從業(yè)者參考。1. 為什么放棄MRT轉(zhuǎn)向腳本化處理1.1 MRT的便利與它的天花板MRTMODIS Reprojection Tool在早期確實(shí)是處理MODIS數(shù)據(jù)的標(biāo)配工具圖形界面友好支持批量導(dǎo)入HDF、選擇波段、設(shè)定投影和輸出格式。對(duì)于偶爾處理幾個(gè)文件的人來說它夠用。但問題在于一旦數(shù)據(jù)量上去MRT的短板就暴露得非常明顯。首先是批處理能力弱。MRT雖然提供了命令行版本但參數(shù)文件.prm的編寫并不直觀尤其是涉及多個(gè)產(chǎn)品、多個(gè)波段組合的時(shí)候每個(gè)組合都要單獨(dú)寫一個(gè)prm文件管理起來很混亂。其次是投影和重采樣選項(xiàng)有限MRT內(nèi)置的投影類型雖然覆蓋了常見的幾種但如果你想用自定義的Albers等積投影或者特定的地理坐標(biāo)系配置起來就很別扭。再就是運(yùn)行穩(wěn)定性MRT在處理大區(qū)域、長(zhǎng)時(shí)間序列數(shù)據(jù)時(shí)偶爾會(huì)卡死或者輸出異常值而且錯(cuò)誤提示非常模糊排查起來費(fèi)時(shí)費(fèi)力。我印象最深的一次是處理青藏高原區(qū)域連續(xù)8年的MOD13A2 NDVI數(shù)據(jù)MRT跑到第三年的時(shí)候突然報(bào)了一個(gè)Error in processing tile的提示沒有任何具體信息。后來逐個(gè)文件排查才發(fā)現(xiàn)是某一個(gè)HDF文件的第13波段存在數(shù)據(jù)異常MRT直接崩潰。如果換成腳本處理這種問題可以在代碼里加異常捕獲跳過問題文件繼續(xù)跑不會(huì)影響整體進(jìn)度。1.2 腳本化處理的三個(gè)核心優(yōu)勢(shì)轉(zhuǎn)向Python腳本之后最直接的感受是可控性完全不一樣了。具體來說腳本化處理帶來三個(gè)層面的提升第一是流程可復(fù)用。一套腳本寫好之后換一個(gè)研究區(qū)、換一個(gè)時(shí)間段只需要改幾個(gè)參數(shù)就能直接跑不需要重新配置工具。這對(duì)于需要反復(fù)做實(shí)驗(yàn)、調(diào)整參數(shù)的研究場(chǎng)景來說效率提升是數(shù)量級(jí)的。第二是異?？勺匪?。代碼里可以加日志記錄、異常捕獲、中間結(jié)果檢查每一步都有據(jù)可查。哪個(gè)文件處理失敗了、失敗原因是什么、跳過了哪些數(shù)據(jù)全都清清楚楚。這比MRT那種黑箱式報(bào)錯(cuò)要友好太多。第三是與其他工具鏈無縫銜接。Python處理完的中間結(jié)果可以直接喂給ArcGIS做空間分析或者用GDAL、Rasterio做進(jìn)一步處理整個(gè)鏈條是打通的。而MRT的輸出往往還需要額外的格式轉(zhuǎn)換步驟才能進(jìn)入下一步分析。提示如果你的數(shù)據(jù)量在10個(gè)文件以內(nèi)MRT確實(shí)夠用。但只要超過這個(gè)量級(jí)或者需要反復(fù)處理不同區(qū)域的數(shù)據(jù)腳本化就是必然選擇。1.3 整體技術(shù)路線概覽這套流程的整體思路是這樣的HDF原始文件 → 波段提取與拼接 → 投影轉(zhuǎn)換與重采樣 → 研究區(qū)裁剪 → NDVI/ET計(jì)算 → 成品圖輸出。其中前四步用Python腳本完成后兩步根據(jù)具體需求可以在Python里做也可以導(dǎo)入ArcGIS做后處理和制圖。具體用到的工具和庫(kù)包括環(huán)節(jié)工具/庫(kù)作用HDF讀取PyHDF / GDAL讀取MODIS HDF4格式數(shù)據(jù)波段提取GDAL提取指定波段并轉(zhuǎn)為GeoTIFF拼接GDAL Warp / Rasterio多景影像拼接投影轉(zhuǎn)換GDAL Warp重投影到目標(biāo)坐標(biāo)系裁剪Rasterio / ArcPy按研究區(qū)邊界裁剪NDVI計(jì)算NumPy / Rasterio波段運(yùn)算ET計(jì)算NumPy / Rasterio基于SEBAL或其他模型制圖輸出ArcGIS / Matplotlib成品圖渲染與導(dǎo)出這條路線的好處是每一環(huán)都可以獨(dú)立調(diào)試出了問題容易定位。而且Python生態(tài)里的這些庫(kù)都是成熟穩(wěn)定的社區(qū)支持也好遇到問題基本都能找到解決方案。2. 環(huán)境搭建與依賴庫(kù)的取舍2.1 Python環(huán)境的選擇處理MODIS數(shù)據(jù)對(duì)Python環(huán)境的要求其實(shí)不算高但有幾個(gè)庫(kù)的安裝比較講究。我建議用Anaconda來管理環(huán)境因?yàn)镚DAL、Rasterio這些庫(kù)在Windows下用pip安裝經(jīng)常出問題而conda源里有預(yù)編譯好的版本省去很多麻煩。創(chuàng)建一個(gè)獨(dú)立的環(huán)境conda create -n modis python3.9 conda activate modis為什么選Python 3.9而不是更新的版本因?yàn)镚DAL和PyHDF在某些新版本Python上的兼容性還不穩(wěn)定3.9是目前最穩(wěn)妥的選擇。我試過3.11PyHDF編譯報(bào)錯(cuò)折騰了半天還是退回3.9。2.2 核心依賴庫(kù)安裝核心依賴庫(kù)包括以下幾個(gè)conda install -c conda-forge gdal rasterio pyhdf numpy matplotlib這里重點(diǎn)說一下PyHDF和GDAL的分工。PyHDF專門用來讀取HDF4格式的文件MODIS的很多產(chǎn)品比如MOD13、MOD16都是HDF4格式用PyHDF讀取最直接。GDAL雖然也能讀HDF但對(duì)HDF4的支持需要額外編譯有時(shí)候會(huì)出問題。所以我的做法是用PyHDF讀數(shù)據(jù)和元信息用GDAL做投影轉(zhuǎn)換和格式輸出。Rasterio是基于GDAL的Python封裝API更友好適合做裁剪和波段運(yùn)算。NumPy不用說了做數(shù)組運(yùn)算的基礎(chǔ)。注意如果你在Windows下用conda安裝GDAL失敗可以試試先安裝conda install -c conda-forge gdal3.4指定版本或者用mamba替代conda解決依賴沖突的能力更強(qiáng)。2.3 ArcGIS的角色定位ArcGIS在這套流程里主要承擔(dān)兩個(gè)角色后處理和制圖輸出。后處理包括一些Python腳本不太方便做的操作比如柵格計(jì)算器的復(fù)雜表達(dá)式、鄰域分析、重分類等。制圖輸出則是ArcGIS的強(qiáng)項(xiàng)做出來的圖專業(yè)、美觀適合直接放到論文或報(bào)告里。如果你用的是ArcGIS Pro它自帶的Python環(huán)境是ArcGIS Pro的conda環(huán)境可以直接調(diào)用ArcPy。但要注意ArcGIS Pro的Python環(huán)境和你自己創(chuàng)建的conda環(huán)境是分開的如果要在ArcPy里調(diào)用GDAL需要在ArcGIS Pro的環(huán)境里也安裝GDAL。我的建議是分開處理Python腳本用自己的環(huán)境跑輸出GeoTIFF后再導(dǎo)入ArcGIS做后續(xù)處理避免環(huán)境沖突。3. HDF文件的批量讀取與波段提取3.1 MODIS HDF文件的內(nèi)部結(jié)構(gòu)MODIS的HDF文件本質(zhì)上是一個(gè)容器里面包含了多個(gè)數(shù)據(jù)集SDSScientific Data Sets。以MOD13A2NDVI產(chǎn)品為例一個(gè)HDF文件里包含NDVI歸一化植被指數(shù)EVI增強(qiáng)植被指數(shù)VI_Quality質(zhì)量波段pixel_reliability像元可靠性sur_refl_b01到sur_refl_b07七個(gè)地表反射率波段角度信息太陽(yáng)天頂角、觀測(cè)天頂角等每個(gè)數(shù)據(jù)集都有自己的維度和數(shù)據(jù)類型。NDVI通常是int16類型需要乘以縮放因子scale factor才能得到真實(shí)值。MOD13A2的NDVI縮放因子是0.0001也就是說原始值10000對(duì)應(yīng)真實(shí)NDVI值1.0。用PyHDF查看文件結(jié)構(gòu)from pyhdf.SD import SD, SDC hdf SD(MOD13A2.A2020001.h26v05.006.2020015000000.hdf, SDC.READ) datasets hdf.datasets() for name, info in datasets.items(): print(f數(shù)據(jù)集: {name}, 維度: {info[1]}, 類型: {info[3]})這段代碼會(huì)列出文件里所有數(shù)據(jù)集的名稱、維度和數(shù)據(jù)類型。第一次處理某個(gè)產(chǎn)品的時(shí)候建議先跑一下這個(gè)確認(rèn)波段名稱和結(jié)構(gòu)因?yàn)椴煌a(chǎn)品的波段命名規(guī)則不一樣。3.2 批量讀取與波段提取腳本批量處理的核心邏輯是遍歷文件夾下所有HDF文件 → 逐個(gè)讀取指定波段 → 輸出為GeoTIFF。這里有個(gè)關(guān)鍵點(diǎn)MODIS的HDF文件本身帶有地理定位信息但PyHDF讀取出來的只是純數(shù)組沒有地理坐標(biāo)。所以需要額外讀取經(jīng)緯度信息或者用GDAL來讀取帶地理信息的波段。我的做法是用GDAL直接讀取HDF中的子數(shù)據(jù)集這樣輸出的GeoTIFF自帶地理坐標(biāo)import os import glob from osgeo import gdal def extract_band(hdf_path, band_name, output_path): 從HDF文件中提取指定波段并輸出為GeoTIFF # 構(gòu)建子數(shù)據(jù)集路徑 subdataset fHDF4_EOS:EOS_GRID:{hdf_path}:MODIS_Grid_16DAY_1km_VI:{band_name} # 打開子數(shù)據(jù)集 ds gdal.Open(subdataset) if ds is None: print(f無法打開: {subdataset}) return False # 輸出為GeoTIFF gdal.Translate(output_path, ds, formatGTiff) ds None return True # 批量處理 input_dir rD:\MODIS\MOD13A2 output_dir rD:\MODIS\NDVI_tif os.makedirs(output_dir, exist_okTrue) hdf_files glob.glob(os.path.join(input_dir, *.hdf)) for hdf_file in hdf_files: filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) print(f已處理: {filename})這里的關(guān)鍵是子數(shù)據(jù)集路徑的構(gòu)建。HDF4_EOS:EOS_GRID是GDAL讀取HDF4的固定前綴后面的MODIS_Grid_16DAY_1km_VI是網(wǎng)格名稱不同產(chǎn)品的網(wǎng)格名稱不一樣。MOD13A2是MODIS_Grid_16DAY_1km_VIMOD16A2ET產(chǎn)品是MODIS_Grid_8Day_1km_LSTE或者M(jìn)ODIS_Grid_8Day_1km_ET具體要看產(chǎn)品文檔。提示如果不確定網(wǎng)格名稱可以用gdalinfo命令查看HDF文件的信息里面會(huì)列出所有子數(shù)據(jù)集的完整路徑。3.3 縮放因子與無效值處理提取出來的NDVI數(shù)據(jù)是int16類型值域通常在-2000到10000之間。要得到真實(shí)的NDVI值需要乘以縮放因子0.0001。同時(shí)MODIS數(shù)據(jù)有專門的無效值填充fill value通常是-3000需要處理掉。import numpy as np import rasterio def apply_scale_factor(input_tif, output_tif, scale0.0001, fill_value-3000): 應(yīng)用縮放因子并處理無效值 with rasterio.open(input_tif) as src: data src.read(1).astype(np.float32) profile src.profile.copy() # 處理無效值 data[data fill_value] np.nan # 應(yīng)用縮放因子 data data * scale # 更新數(shù)據(jù)類型 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(data, 1)這一步看起來簡(jiǎn)單但不做縮放的話后續(xù)所有分析都是錯(cuò)的。我見過不少初學(xué)者直接拿原始int16值去做時(shí)序分析結(jié)果NDVI值域變成-2000到10000完全沒法用。所以這個(gè)步驟一定要養(yǎng)成習(xí)慣。4. 拼接、投影轉(zhuǎn)換與研究區(qū)裁剪4.1 多景影像的拼接邏輯MODIS數(shù)據(jù)是按瓦片tile組織的全球被劃分為36×18個(gè)瓦片每個(gè)瓦片覆蓋約1100km×1100km。如果你的研究區(qū)跨越多個(gè)瓦片就需要先拼接。比如青藏高原區(qū)域通常涉及h25v05、h26v05、h25v06、h26v06四個(gè)瓦片。拼接用GDAL的Warp工具最方便from osgeo import gdal def mosaic_tiles(input_files, output_path): 拼接多個(gè)瓦片 options gdal.WarpOptions( formatGTiff, resampleAlggdal.GRA_NearestNeighbour, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_path, input_files, **options)這里重采樣方法的選擇很關(guān)鍵。對(duì)于NDVI這種連續(xù)值數(shù)據(jù)理論上用雙線性插值GRA_Bilinear更平滑但會(huì)引入不屬于原始數(shù)據(jù)的新值。對(duì)于分類數(shù)據(jù)或者質(zhì)量波段必須用最近鄰GRA_NearestNeighbour。我的習(xí)慣是NDVI用雙線性質(zhì)量波段用最近鄰這樣既保證平滑又不破壞分類信息。4.2 投影轉(zhuǎn)換的參數(shù)設(shè)定MODIS原始數(shù)據(jù)用的是正弦投影Sinusoidal這種投影在低緯度地區(qū)變形較小但在高緯度地區(qū)面積變形嚴(yán)重。做區(qū)域分析的時(shí)候通常需要轉(zhuǎn)換到更適合的投影比如Albers等積投影或者UTM投影。以中國(guó)區(qū)域?yàn)槔Ｓ玫腁lbers投影參數(shù)是albers_proj projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs def reproject(input_tif, output_tif, dst_proj): 投影轉(zhuǎn)換 options gdal.WarpOptions( formatGTiff, dstSRSdst_proj, resampleAlggdal.GRA_Bilinear, xRes1000, yRes1000, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_tif, input_tif, **options)分辨率的選擇要根據(jù)原始數(shù)據(jù)來定。MOD13A2的空間分辨率是1km所以重采樣到1000m是合理的。如果你重采樣到500m雖然文件變大了但信息量并沒有增加反而引入了插值誤差。這一點(diǎn)很多人容易忽略覺得分辨率越高越好其實(shí)不是。4.3 按研究區(qū)邊界精確裁剪裁剪需要兩部分?jǐn)?shù)據(jù)待裁剪的柵格和研究區(qū)邊界矢量。用Rasterio的mask功能可以很方便地實(shí)現(xiàn)import rasterio from rasterio.mask import mask import geopandas as gpd def clip_raster(input_tif, boundary_shp, output_tif): 按矢量邊界裁剪柵格 # 讀取邊界 gdf gpd.read_file(boundary_shp) geometries gdf.geometry.values with rasterio.open(input_tif) as src: # 裁剪 out_image, out_transform mask(src, geometries, cropTrue) out_meta src.meta.copy() # 更新元數(shù)據(jù) out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(output_tif, w, **out_meta) as dst: dst.write(out_image)這里有個(gè)常見的坑如果研究區(qū)邊界是經(jīng)緯度坐標(biāo)WGS84而柵格是投影坐標(biāo)Albers直接裁剪會(huì)報(bào)錯(cuò)或者結(jié)果為空。解決辦法是先統(tǒng)一坐標(biāo)系把邊界矢量轉(zhuǎn)換到和柵格一樣的投影下gdf gdf.to_crs(src.crs)這一步看起來簡(jiǎn)單但實(shí)際工作中經(jīng)常忘記導(dǎo)致裁剪結(jié)果為空排查半天才發(fā)現(xiàn)是坐標(biāo)系不匹配。注意裁剪的時(shí)候cropTrue會(huì)把柵格裁剪到邊界的最小外接矩形邊界外的像元會(huì)被設(shè)為nodata。如果你需要精確到邊界形狀還需要額外做掩膜處理。5. NDVI與ET成品的計(jì)算與后處理5.1 NDVI的時(shí)序合成與最大值合成法NDVI數(shù)據(jù)通常是16天合成的一年有23期。做年度分析的時(shí)候常用的方法是最大值合成法MVCMaximum Value Composite即取一年中每個(gè)像元的最大NDVI值這樣可以最大程度消除云污染和大氣影響。import numpy as np import rasterio import glob def mvc_composite(ndvi_files, output_path): 最大值合成 # 讀取所有期數(shù)據(jù) arrays [] for f in sorted(ndvi_files): with rasterio.open(f) as src: arrays.append(src.read(1)) profile src.profile.copy() # 堆疊并取最大值 stack np.stack(arrays, axis0) max_ndvi np.nanmax(stack, axis0) # 輸出 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: dst.write(max_ndvi.astype(np.float32), 1)MVC的邏輯很直觀一年中NDVI最高的那期大概率是植被生長(zhǎng)最旺盛、云污染最少的時(shí)候。這個(gè)方法雖然簡(jiǎn)單但在植被動(dòng)態(tài)監(jiān)測(cè)中非常有效也是文獻(xiàn)里最常用的方法之一。5.2 ET產(chǎn)品的處理要點(diǎn)MODIS的ET產(chǎn)品MOD16A2和NDVI產(chǎn)品在處理上有幾個(gè)不同點(diǎn)需要注意第一是數(shù)據(jù)單位。MOD16A2的ET數(shù)據(jù)單位是kg/m2/8day也就是每8天的蒸散量。要換算成日均ET需要除以8。如果要換算成mm由于水的密度是1kg/L1kg/m2等于1mm所以數(shù)值上是一樣的。第二是質(zhì)量波段。MOD16A2有專門的質(zhì)量控制波段ET_QC需要根據(jù)QC值篩選有效數(shù)據(jù)。QC值是一個(gè)8位整數(shù)不同的位代表不同的質(zhì)量信息。通常的做法是只保留QC值小于某個(gè)閾值的像元。第三是無效值處理。MOD16A2的填充值是32767需要單獨(dú)處理。def process_et(et_tif, qc_tif, output_tif): 處理ET數(shù)據(jù)篩選有效像元 with rasterio.open(et_tif) as src: et src.read(1).astype(np.float32) profile src.profile.copy() with rasterio.open(qc_tif) as src: qc src.read(1) # 處理填充值 et[et 32767] np.nan # 根據(jù)QC篩選保留QC小于64的像元即質(zhì)量較好的 et[qc 64] np.nan # 轉(zhuǎn)換為日均ET除以8 et et / 8.0 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(et, 1)5.3 在ArcGIS中做后處理與制圖Python處理完的GeoTIFF導(dǎo)入ArcGIS后可以做進(jìn)一步的后處理和制圖。幾個(gè)常用的操作柵格計(jì)算器做歸一化、重分類、閾值提取等。比如把NDVI小于0.1的像元設(shè)為0非植被大于0.1的保留。符號(hào)化NDVI通常用綠-黃-紅的漸變色帶ET用藍(lán)-綠-黃的色帶。ArcGIS的色帶庫(kù)里有現(xiàn)成的也可以自定義。布局出圖添加指北針、比例尺、圖例、經(jīng)緯網(wǎng)格導(dǎo)出為高分辨率PNG或PDF。這一步是ArcGIS的強(qiáng)項(xiàng)做出來的圖比Matplotlib要專業(yè)得多。提示如果要在ArcGIS里批量出圖可以用ArcPy寫腳本結(jié)合地圖文檔.mxd模板自動(dòng)替換數(shù)據(jù)源并導(dǎo)出。這在做多期對(duì)比圖的時(shí)候特別有用。6. 批量處理中的踩坑記錄與性能優(yōu)化6.1 內(nèi)存溢出與分塊處理處理大區(qū)域、長(zhǎng)時(shí)間序列數(shù)據(jù)的時(shí)候最容易遇到的問題就是內(nèi)存溢出。尤其是做MVC合成的時(shí)候如果把23期1km分辨率的全國(guó)數(shù)據(jù)全部讀進(jìn)內(nèi)存光是一個(gè)波段就要占幾個(gè)GB加上中間變量很容易爆內(nèi)存。解決辦法是分塊處理。Rasterio支持按窗口讀取數(shù)據(jù)def mvc_composite_block(ndvi_files, output_path, block_size1024): 分塊最大值合成 with rasterio.open(ndvi_files[0]) as src: profile src.profile.copy() height, width src.height, src.width profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: for i in range(0, height, block_size): for j in range(0, width, block_size): # 計(jì)算當(dāng)前塊的范圍 win_height min(block_size, height - i) win_width min(block_size, width - j) window rasterio.windows.Window(j, i, win_width, win_height) # 讀取所有期的當(dāng)前塊 blocks [] for f in ndvi_files: with rasterio.open(f) as src: blocks.append(src.read(1, windowwindow)) # 取最大值 stack np.stack(blocks, axis0) max_block np.nanmax(stack, axis0) # 寫入 dst.write(max_block.astype(np.float32), 1, windowwindow)分塊處理雖然代碼復(fù)雜一點(diǎn)但內(nèi)存占用可以控制在幾百M(fèi)B以內(nèi)處理全國(guó)數(shù)據(jù)也沒問題。塊大小的選擇要根據(jù)你的內(nèi)存來定一般1024×1024或者2048×2048比較合適。6.2 坐標(biāo)系不匹配導(dǎo)致的裁剪失敗這個(gè)坑我在前面提過但值得再?gòu)?qiáng)調(diào)一次。裁剪失敗最常見的原因就是坐標(biāo)系不匹配。癥狀是裁剪結(jié)果為空或者只有一小塊。排查方法是用gdalinfo查看柵格的坐標(biāo)系用ogrinfo查看矢量的坐標(biāo)系確認(rèn)兩者是否一致如果不一致用gdf.to_crs()轉(zhuǎn)換矢量或者用gdal.Warp轉(zhuǎn)換柵格。我個(gè)人的習(xí)慣是統(tǒng)一用投影坐標(biāo)系因?yàn)橥队白鴺?biāo)系下的距離和面積計(jì)算更準(zhǔn)確。6.3 并行處理加速批量任務(wù)如果數(shù)據(jù)量很大單線程處理太慢可以用Python的multiprocessing庫(kù)做并行。比如批量提取波段的時(shí)候每個(gè)文件獨(dú)立處理天然適合并行from multiprocessing import Pool import os def process_single_file(hdf_file): 處理單個(gè)文件 filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) return filename if __name__ __main__: hdf_files glob.glob(os.path.join(input_dir, *.hdf)) with Pool(processes4) as pool: results pool.map(process_single_file, hdf_files) print(f完成 {len(results)} 個(gè)文件)進(jìn)程數(shù)根據(jù)你的CPU核心數(shù)來定一般設(shè)為核心數(shù)的一半到全部。注意GDAL在多進(jìn)程下有時(shí)候會(huì)有問題如果遇到崩潰可以改用ThreadPool或者降低進(jìn)程數(shù)。6.4 常見錯(cuò)誤與快速排查表錯(cuò)誤現(xiàn)象可能原因解決辦法裁剪結(jié)果為空坐標(biāo)系不匹配統(tǒng)一矢量與柵格坐標(biāo)系NDVI值域異常未應(yīng)用縮放因子乘以0.0001內(nèi)存溢出一次性讀取全部數(shù)據(jù)分塊處理HDF讀取失敗網(wǎng)格名稱錯(cuò)誤用gdalinfo查看子數(shù)據(jù)集路徑投影轉(zhuǎn)換后變形重采樣方法不當(dāng)連續(xù)值用雙線性分類值用最近鄰并行處理崩潰GDAL多進(jìn)程沖突降低進(jìn)程數(shù)或改用線程這張表是我在實(shí)際工作中總結(jié)出來的基本上覆蓋了80%以上的常見問題。遇到報(bào)錯(cuò)的時(shí)候先對(duì)照這張表排查能省不少時(shí)間。7. 從腳本到工作流我的實(shí)際使用體會(huì)整套流程跑通之后我把它封裝成了一個(gè)配置文件驅(qū)動(dòng)的工作流。所有的輸入路徑、輸出路徑、投影參數(shù)、分辨率、縮放因子都寫在一個(gè)YAML文件里主腳本讀取配置后自動(dòng)執(zhí)行。這樣換研究區(qū)的時(shí)候只需要改配置文件不用動(dòng)代碼。# config.yaml input_dir: D:/MODIS/MOD13A2 output_dir: D:/MODIS/output product: MOD13A2 band: NDVI scale_factor: 0.0001 fill_value: -3000 target_proj: projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm resolution: 1000 boundary: D:/data/study_area.shp這種配置驅(qū)動(dòng)的方式最大的好處是可追溯。半年后回頭看某個(gè)結(jié)果是怎么生成的翻出配置文件就一目了然不用去回憶當(dāng)時(shí)用了什么參數(shù)。另外一個(gè)小技巧是中間結(jié)果保留。很多人為了省磁盤空間處理完就把中間文件刪了。我的建議是至少保留波段提取后的GeoTIFF因?yàn)楹罄m(xù)如果要調(diào)整投影或者裁剪范圍可以直接從這一步重新跑不用從頭再來。磁盤空間現(xiàn)在很便宜但重新處理的時(shí)間成本很高。最后說一個(gè)關(guān)于ArcGIS和Python配合的心得。我現(xiàn)在的習(xí)慣是Python做批量和計(jì)算ArcGIS做交互和出圖。Python腳本負(fù)責(zé)把原始數(shù)據(jù)變成干凈的、帶地理坐標(biāo)的GeoTIFF然后導(dǎo)入ArcGIS做符號(hào)化、布局、導(dǎo)出。兩者各司其職效率最高。如果硬要用ArcPy做所有事情反而會(huì)因?yàn)锳rcGIS的Python環(huán)境限制而束手束腳。這套流程我用了三年多從最初的幾十個(gè)文件到現(xiàn)在每年處理上萬(wàn)景數(shù)據(jù)穩(wěn)定性一直很好。核心邏輯沒有變過只是根據(jù)具體需求做了些微調(diào)。如果你也在做類似的批量柵格處理希望這些經(jīng)驗(yàn)?zāi)軒湍闵僮咝澛贰? Python批量處理MODIS遙感數(shù)據(jù)：從HDF到NDVI/ET成品圖全流程
遙感數(shù)據(jù)處理這件事最怕的不是算法難而是重復(fù)勞動(dòng)。我最早接觸MODIS數(shù)據(jù)的時(shí)候一個(gè)研究區(qū)、一年12個(gè)月、4種產(chǎn)品光是HDF文件的拼接、投影、裁剪、格式轉(zhuǎn)換就能耗掉整整兩天中間還得盯著MRT那個(gè)老界面反復(fù)點(diǎn)選參數(shù)稍不留神就報(bào)錯(cuò)重來。后來被逼得沒辦法硬是把整套流程用Python腳本串了起來配合ArcGIS做后處理和出圖現(xiàn)在同樣的工作量壓縮到半小時(shí)以內(nèi)而且可以批量跑、無人值守。這篇內(nèi)容就是把這套流程完整拆開講清楚從HDF原始文件到NDVI、ET成品圖每一步為什么這么做、參數(shù)怎么定、坑在哪里都會(huì)講到。適合做生態(tài)遙感、農(nóng)業(yè)監(jiān)測(cè)、水文分析方向的同學(xué)也適合任何需要批量處理柵格數(shù)據(jù)的GIS從業(yè)者參考。1. 為什么放棄MRT轉(zhuǎn)向腳本化處理1.1 MRT的便利與它的天花板MRTMODIS Reprojection Tool在早期確實(shí)是處理MODIS數(shù)據(jù)的標(biāo)配工具圖形界面友好支持批量導(dǎo)入HDF、選擇波段、設(shè)定投影和輸出格式。對(duì)于偶爾處理幾個(gè)文件的人來說它夠用。但問題在于一旦數(shù)據(jù)量上去MRT的短板就暴露得非常明顯。首先是批處理能力弱。MRT雖然提供了命令行版本但參數(shù)文件.prm的編寫并不直觀尤其是涉及多個(gè)產(chǎn)品、多個(gè)波段組合的時(shí)候每個(gè)組合都要單獨(dú)寫一個(gè)prm文件管理起來很混亂。其次是投影和重采樣選項(xiàng)有限MRT內(nèi)置的投影類型雖然覆蓋了常見的幾種但如果你想用自定義的Albers等積投影或者特定的地理坐標(biāo)系配置起來就很別扭。再就是運(yùn)行穩(wěn)定性MRT在處理大區(qū)域、長(zhǎng)時(shí)間序列數(shù)據(jù)時(shí)偶爾會(huì)卡死或者輸出異常值而且錯(cuò)誤提示非常模糊排查起來費(fèi)時(shí)費(fèi)力。我印象最深的一次是處理青藏高原區(qū)域連續(xù)8年的MOD13A2 NDVI數(shù)據(jù)MRT跑到第三年的時(shí)候突然報(bào)了一個(gè)Error in processing tile的提示沒有任何具體信息。后來逐個(gè)文件排查才發(fā)現(xiàn)是某一個(gè)HDF文件的第13波段存在數(shù)據(jù)異常MRT直接崩潰。如果換成腳本處理這種問題可以在代碼里加異常捕獲跳過問題文件繼續(xù)跑不會(huì)影響整體進(jìn)度。1.2 腳本化處理的三個(gè)核心優(yōu)勢(shì)轉(zhuǎn)向Python腳本之后最直接的感受是可控性完全不一樣了。具體來說腳本化處理帶來三個(gè)層面的提升第一是流程可復(fù)用。一套腳本寫好之后換一個(gè)研究區(qū)、換一個(gè)時(shí)間段只需要改幾個(gè)參數(shù)就能直接跑不需要重新配置工具。這對(duì)于需要反復(fù)做實(shí)驗(yàn)、調(diào)整參數(shù)的研究場(chǎng)景來說效率提升是數(shù)量級(jí)的。第二是異?？勺匪?。代碼里可以加日志記錄、異常捕獲、中間結(jié)果檢查每一步都有據(jù)可查。哪個(gè)文件處理失敗了、失敗原因是什么、跳過了哪些數(shù)據(jù)全都清清楚楚。這比MRT那種黑箱式報(bào)錯(cuò)要友好太多。第三是與其他工具鏈無縫銜接。Python處理完的中間結(jié)果可以直接喂給ArcGIS做空間分析或者用GDAL、Rasterio做進(jìn)一步處理整個(gè)鏈條是打通的。而MRT的輸出往往還需要額外的格式轉(zhuǎn)換步驟才能進(jìn)入下一步分析。提示如果你的數(shù)據(jù)量在10個(gè)文件以內(nèi)MRT確實(shí)夠用。但只要超過這個(gè)量級(jí)或者需要反復(fù)處理不同區(qū)域的數(shù)據(jù)腳本化就是必然選擇。1.3 整體技術(shù)路線概覽這套流程的整體思路是這樣的HDF原始文件 → 波段提取與拼接 → 投影轉(zhuǎn)換與重采樣 → 研究區(qū)裁剪 → NDVI/ET計(jì)算 → 成品圖輸出。其中前四步用Python腳本完成后兩步根據(jù)具體需求可以在Python里做也可以導(dǎo)入ArcGIS做后處理和制圖。具體用到的工具和庫(kù)包括環(huán)節(jié)工具/庫(kù)作用HDF讀取PyHDF / GDAL讀取MODIS HDF4格式數(shù)據(jù)波段提取GDAL提取指定波段并轉(zhuǎn)為GeoTIFF拼接GDAL Warp / Rasterio多景影像拼接投影轉(zhuǎn)換GDAL Warp重投影到目標(biāo)坐標(biāo)系裁剪Rasterio / ArcPy按研究區(qū)邊界裁剪NDVI計(jì)算NumPy / Rasterio波段運(yùn)算ET計(jì)算NumPy / Rasterio基于SEBAL或其他模型制圖輸出ArcGIS / Matplotlib成品圖渲染與導(dǎo)出這條路線的好處是每一環(huán)都可以獨(dú)立調(diào)試出了問題容易定位。而且Python生態(tài)里的這些庫(kù)都是成熟穩(wěn)定的社區(qū)支持也好遇到問題基本都能找到解決方案。2. 環(huán)境搭建與依賴庫(kù)的取舍2.1 Python環(huán)境的選擇處理MODIS數(shù)據(jù)對(duì)Python環(huán)境的要求其實(shí)不算高但有幾個(gè)庫(kù)的安裝比較講究。我建議用Anaconda來管理環(huán)境因?yàn)镚DAL、Rasterio這些庫(kù)在Windows下用pip安裝經(jīng)常出問題而conda源里有預(yù)編譯好的版本省去很多麻煩。創(chuàng)建一個(gè)獨(dú)立的環(huán)境conda create -n modis python3.9 conda activate modis為什么選Python 3.9而不是更新的版本因?yàn)镚DAL和PyHDF在某些新版本Python上的兼容性還不穩(wěn)定3.9是目前最穩(wěn)妥的選擇。我試過3.11PyHDF編譯報(bào)錯(cuò)折騰了半天還是退回3.9。2.2 核心依賴庫(kù)安裝核心依賴庫(kù)包括以下幾個(gè)conda install -c conda-forge gdal rasterio pyhdf numpy matplotlib這里重點(diǎn)說一下PyHDF和GDAL的分工。PyHDF專門用來讀取HDF4格式的文件MODIS的很多產(chǎn)品比如MOD13、MOD16都是HDF4格式用PyHDF讀取最直接。GDAL雖然也能讀HDF但對(duì)HDF4的支持需要額外編譯有時(shí)候會(huì)出問題。所以我的做法是用PyHDF讀數(shù)據(jù)和元信息用GDAL做投影轉(zhuǎn)換和格式輸出。Rasterio是基于GDAL的Python封裝API更友好適合做裁剪和波段運(yùn)算。NumPy不用說了做數(shù)組運(yùn)算的基礎(chǔ)。注意如果你在Windows下用conda安裝GDAL失敗可以試試先安裝conda install -c conda-forge gdal3.4指定版本或者用mamba替代conda解決依賴沖突的能力更強(qiáng)。2.3 ArcGIS的角色定位ArcGIS在這套流程里主要承擔(dān)兩個(gè)角色后處理和制圖輸出。后處理包括一些Python腳本不太方便做的操作比如柵格計(jì)算器的復(fù)雜表達(dá)式、鄰域分析、重分類等。制圖輸出則是ArcGIS的強(qiáng)項(xiàng)做出來的圖專業(yè)、美觀適合直接放到論文或報(bào)告里。如果你用的是ArcGIS Pro它自帶的Python環(huán)境是ArcGIS Pro的conda環(huán)境可以直接調(diào)用ArcPy。但要注意ArcGIS Pro的Python環(huán)境和你自己創(chuàng)建的conda環(huán)境是分開的如果要在ArcPy里調(diào)用GDAL需要在ArcGIS Pro的環(huán)境里也安裝GDAL。我的建議是分開處理Python腳本用自己的環(huán)境跑輸出GeoTIFF后再導(dǎo)入ArcGIS做后續(xù)處理避免環(huán)境沖突。3. HDF文件的批量讀取與波段提取3.1 MODIS HDF文件的內(nèi)部結(jié)構(gòu)MODIS的HDF文件本質(zhì)上是一個(gè)容器里面包含了多個(gè)數(shù)據(jù)集SDSScientific Data Sets。以MOD13A2NDVI產(chǎn)品為例一個(gè)HDF文件里包含NDVI歸一化植被指數(shù)EVI增強(qiáng)植被指數(shù)VI_Quality質(zhì)量波段pixel_reliability像元可靠性sur_refl_b01到sur_refl_b07七個(gè)地表反射率波段角度信息太陽(yáng)天頂角、觀測(cè)天頂角等每個(gè)數(shù)據(jù)集都有自己的維度和數(shù)據(jù)類型。NDVI通常是int16類型需要乘以縮放因子scale factor才能得到真實(shí)值。MOD13A2的NDVI縮放因子是0.0001也就是說原始值10000對(duì)應(yīng)真實(shí)NDVI值1.0。用PyHDF查看文件結(jié)構(gòu)from pyhdf.SD import SD, SDC hdf SD(MOD13A2.A2020001.h26v05.006.2020015000000.hdf, SDC.READ) datasets hdf.datasets() for name, info in datasets.items(): print(f數(shù)據(jù)集: {name}, 維度: {info[1]}, 類型: {info[3]})這段代碼會(huì)列出文件里所有數(shù)據(jù)集的名稱、維度和數(shù)據(jù)類型。第一次處理某個(gè)產(chǎn)品的時(shí)候建議先跑一下這個(gè)確認(rèn)波段名稱和結(jié)構(gòu)因?yàn)椴煌a(chǎn)品的波段命名規(guī)則不一樣。3.2 批量讀取與波段提取腳本批量處理的核心邏輯是遍歷文件夾下所有HDF文件 → 逐個(gè)讀取指定波段 → 輸出為GeoTIFF。這里有個(gè)關(guān)鍵點(diǎn)MODIS的HDF文件本身帶有地理定位信息但PyHDF讀取出來的只是純數(shù)組沒有地理坐標(biāo)。所以需要額外讀取經(jīng)緯度信息或者用GDAL來讀取帶地理信息的波段。我的做法是用GDAL直接讀取HDF中的子數(shù)據(jù)集這樣輸出的GeoTIFF自帶地理坐標(biāo)import os import glob from osgeo import gdal def extract_band(hdf_path, band_name, output_path): 從HDF文件中提取指定波段并輸出為GeoTIFF # 構(gòu)建子數(shù)據(jù)集路徑 subdataset fHDF4_EOS:EOS_GRID:{hdf_path}:MODIS_Grid_16DAY_1km_VI:{band_name} # 打開子數(shù)據(jù)集 ds gdal.Open(subdataset) if ds is None: print(f無法打開: {subdataset}) return False # 輸出為GeoTIFF gdal.Translate(output_path, ds, formatGTiff) ds None return True # 批量處理 input_dir rD:\MODIS\MOD13A2 output_dir rD:\MODIS\NDVI_tif os.makedirs(output_dir, exist_okTrue) hdf_files glob.glob(os.path.join(input_dir, *.hdf)) for hdf_file in hdf_files: filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) print(f已處理: {filename})這里的關(guān)鍵是子數(shù)據(jù)集路徑的構(gòu)建。HDF4_EOS:EOS_GRID是GDAL讀取HDF4的固定前綴后面的MODIS_Grid_16DAY_1km_VI是網(wǎng)格名稱不同產(chǎn)品的網(wǎng)格名稱不一樣。MOD13A2是MODIS_Grid_16DAY_1km_VIMOD16A2ET產(chǎn)品是MODIS_Grid_8Day_1km_LSTE或者M(jìn)ODIS_Grid_8Day_1km_ET具體要看產(chǎn)品文檔。提示如果不確定網(wǎng)格名稱可以用gdalinfo命令查看HDF文件的信息里面會(huì)列出所有子數(shù)據(jù)集的完整路徑。3.3 縮放因子與無效值處理提取出來的NDVI數(shù)據(jù)是int16類型值域通常在-2000到10000之間。要得到真實(shí)的NDVI值需要乘以縮放因子0.0001。同時(shí)MODIS數(shù)據(jù)有專門的無效值填充fill value通常是-3000需要處理掉。import numpy as np import rasterio def apply_scale_factor(input_tif, output_tif, scale0.0001, fill_value-3000): 應(yīng)用縮放因子并處理無效值 with rasterio.open(input_tif) as src: data src.read(1).astype(np.float32) profile src.profile.copy() # 處理無效值 data[data fill_value] np.nan # 應(yīng)用縮放因子 data data * scale # 更新數(shù)據(jù)類型 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(data, 1)這一步看起來簡(jiǎn)單但不做縮放的話后續(xù)所有分析都是錯(cuò)的。我見過不少初學(xué)者直接拿原始int16值去做時(shí)序分析結(jié)果NDVI值域變成-2000到10000完全沒法用。所以這個(gè)步驟一定要養(yǎng)成習(xí)慣。4. 拼接、投影轉(zhuǎn)換與研究區(qū)裁剪4.1 多景影像的拼接邏輯MODIS數(shù)據(jù)是按瓦片tile組織的全球被劃分為36×18個(gè)瓦片每個(gè)瓦片覆蓋約1100km×1100km。如果你的研究區(qū)跨越多個(gè)瓦片就需要先拼接。比如青藏高原區(qū)域通常涉及h25v05、h26v05、h25v06、h26v06四個(gè)瓦片。拼接用GDAL的Warp工具最方便from osgeo import gdal def mosaic_tiles(input_files, output_path): 拼接多個(gè)瓦片 options gdal.WarpOptions( formatGTiff, resampleAlggdal.GRA_NearestNeighbour, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_path, input_files, **options)這里重采樣方法的選擇很關(guān)鍵。對(duì)于NDVI這種連續(xù)值數(shù)據(jù)理論上用雙線性插值GRA_Bilinear更平滑但會(huì)引入不屬于原始數(shù)據(jù)的新值。對(duì)于分類數(shù)據(jù)或者質(zhì)量波段必須用最近鄰GRA_NearestNeighbour。我的習(xí)慣是NDVI用雙線性質(zhì)量波段用最近鄰這樣既保證平滑又不破壞分類信息。4.2 投影轉(zhuǎn)換的參數(shù)設(shè)定MODIS原始數(shù)據(jù)用的是正弦投影Sinusoidal這種投影在低緯度地區(qū)變形較小但在高緯度地區(qū)面積變形嚴(yán)重。做區(qū)域分析的時(shí)候通常需要轉(zhuǎn)換到更適合的投影比如Albers等積投影或者UTM投影。以中國(guó)區(qū)域?yàn)槔Ｓ玫腁lbers投影參數(shù)是albers_proj projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs def reproject(input_tif, output_tif, dst_proj): 投影轉(zhuǎn)換 options gdal.WarpOptions( formatGTiff, dstSRSdst_proj, resampleAlggdal.GRA_Bilinear, xRes1000, yRes1000, srcNodata-3000, dstNodata-3000 ) gdal.Warp(output_tif, input_tif, **options)分辨率的選擇要根據(jù)原始數(shù)據(jù)來定。MOD13A2的空間分辨率是1km所以重采樣到1000m是合理的。如果你重采樣到500m雖然文件變大了但信息量并沒有增加反而引入了插值誤差。這一點(diǎn)很多人容易忽略覺得分辨率越高越好其實(shí)不是。4.3 按研究區(qū)邊界精確裁剪裁剪需要兩部分?jǐn)?shù)據(jù)待裁剪的柵格和研究區(qū)邊界矢量。用Rasterio的mask功能可以很方便地實(shí)現(xiàn)import rasterio from rasterio.mask import mask import geopandas as gpd def clip_raster(input_tif, boundary_shp, output_tif): 按矢量邊界裁剪柵格 # 讀取邊界 gdf gpd.read_file(boundary_shp) geometries gdf.geometry.values with rasterio.open(input_tif) as src: # 裁剪 out_image, out_transform mask(src, geometries, cropTrue) out_meta src.meta.copy() # 更新元數(shù)據(jù) out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(output_tif, w, **out_meta) as dst: dst.write(out_image)這里有個(gè)常見的坑如果研究區(qū)邊界是經(jīng)緯度坐標(biāo)WGS84而柵格是投影坐標(biāo)Albers直接裁剪會(huì)報(bào)錯(cuò)或者結(jié)果為空。解決辦法是先統(tǒng)一坐標(biāo)系把邊界矢量轉(zhuǎn)換到和柵格一樣的投影下gdf gdf.to_crs(src.crs)這一步看起來簡(jiǎn)單但實(shí)際工作中經(jīng)常忘記導(dǎo)致裁剪結(jié)果為空排查半天才發(fā)現(xiàn)是坐標(biāo)系不匹配。注意裁剪的時(shí)候cropTrue會(huì)把柵格裁剪到邊界的最小外接矩形邊界外的像元會(huì)被設(shè)為nodata。如果你需要精確到邊界形狀還需要額外做掩膜處理。5. NDVI與ET成品的計(jì)算與后處理5.1 NDVI的時(shí)序合成與最大值合成法NDVI數(shù)據(jù)通常是16天合成的一年有23期。做年度分析的時(shí)候常用的方法是最大值合成法MVCMaximum Value Composite即取一年中每個(gè)像元的最大NDVI值這樣可以最大程度消除云污染和大氣影響。import numpy as np import rasterio import glob def mvc_composite(ndvi_files, output_path): 最大值合成 # 讀取所有期數(shù)據(jù) arrays [] for f in sorted(ndvi_files): with rasterio.open(f) as src: arrays.append(src.read(1)) profile src.profile.copy() # 堆疊并取最大值 stack np.stack(arrays, axis0) max_ndvi np.nanmax(stack, axis0) # 輸出 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: dst.write(max_ndvi.astype(np.float32), 1)MVC的邏輯很直觀一年中NDVI最高的那期大概率是植被生長(zhǎng)最旺盛、云污染最少的時(shí)候。這個(gè)方法雖然簡(jiǎn)單但在植被動(dòng)態(tài)監(jiān)測(cè)中非常有效也是文獻(xiàn)里最常用的方法之一。5.2 ET產(chǎn)品的處理要點(diǎn)MODIS的ET產(chǎn)品MOD16A2和NDVI產(chǎn)品在處理上有幾個(gè)不同點(diǎn)需要注意第一是數(shù)據(jù)單位。MOD16A2的ET數(shù)據(jù)單位是kg/m2/8day也就是每8天的蒸散量。要換算成日均ET需要除以8。如果要換算成mm由于水的密度是1kg/L1kg/m2等于1mm所以數(shù)值上是一樣的。第二是質(zhì)量波段。MOD16A2有專門的質(zhì)量控制波段ET_QC需要根據(jù)QC值篩選有效數(shù)據(jù)。QC值是一個(gè)8位整數(shù)不同的位代表不同的質(zhì)量信息。通常的做法是只保留QC值小于某個(gè)閾值的像元。第三是無效值處理。MOD16A2的填充值是32767需要單獨(dú)處理。def process_et(et_tif, qc_tif, output_tif): 處理ET數(shù)據(jù)篩選有效像元 with rasterio.open(et_tif) as src: et src.read(1).astype(np.float32) profile src.profile.copy() with rasterio.open(qc_tif) as src: qc src.read(1) # 處理填充值 et[et 32767] np.nan # 根據(jù)QC篩選保留QC小于64的像元即質(zhì)量較好的 et[qc 64] np.nan # 轉(zhuǎn)換為日均ET除以8 et et / 8.0 profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_tif, w, **profile) as dst: dst.write(et, 1)5.3 在ArcGIS中做后處理與制圖Python處理完的GeoTIFF導(dǎo)入ArcGIS后可以做進(jìn)一步的后處理和制圖。幾個(gè)常用的操作柵格計(jì)算器做歸一化、重分類、閾值提取等。比如把NDVI小于0.1的像元設(shè)為0非植被大于0.1的保留。符號(hào)化NDVI通常用綠-黃-紅的漸變色帶ET用藍(lán)-綠-黃的色帶。ArcGIS的色帶庫(kù)里有現(xiàn)成的也可以自定義。布局出圖添加指北針、比例尺、圖例、經(jīng)緯網(wǎng)格導(dǎo)出為高分辨率PNG或PDF。這一步是ArcGIS的強(qiáng)項(xiàng)做出來的圖比Matplotlib要專業(yè)得多。提示如果要在ArcGIS里批量出圖可以用ArcPy寫腳本結(jié)合地圖文檔.mxd模板自動(dòng)替換數(shù)據(jù)源并導(dǎo)出。這在做多期對(duì)比圖的時(shí)候特別有用。6. 批量處理中的踩坑記錄與性能優(yōu)化6.1 內(nèi)存溢出與分塊處理處理大區(qū)域、長(zhǎng)時(shí)間序列數(shù)據(jù)的時(shí)候最容易遇到的問題就是內(nèi)存溢出。尤其是做MVC合成的時(shí)候如果把23期1km分辨率的全國(guó)數(shù)據(jù)全部讀進(jìn)內(nèi)存光是一個(gè)波段就要占幾個(gè)GB加上中間變量很容易爆內(nèi)存。解決辦法是分塊處理。Rasterio支持按窗口讀取數(shù)據(jù)def mvc_composite_block(ndvi_files, output_path, block_size1024): 分塊最大值合成 with rasterio.open(ndvi_files[0]) as src: profile src.profile.copy() height, width src.height, src.width profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: for i in range(0, height, block_size): for j in range(0, width, block_size): # 計(jì)算當(dāng)前塊的范圍 win_height min(block_size, height - i) win_width min(block_size, width - j) window rasterio.windows.Window(j, i, win_width, win_height) # 讀取所有期的當(dāng)前塊 blocks [] for f in ndvi_files: with rasterio.open(f) as src: blocks.append(src.read(1, windowwindow)) # 取最大值 stack np.stack(blocks, axis0) max_block np.nanmax(stack, axis0) # 寫入 dst.write(max_block.astype(np.float32), 1, windowwindow)分塊處理雖然代碼復(fù)雜一點(diǎn)但內(nèi)存占用可以控制在幾百M(fèi)B以內(nèi)處理全國(guó)數(shù)據(jù)也沒問題。塊大小的選擇要根據(jù)你的內(nèi)存來定一般1024×1024或者2048×2048比較合適。6.2 坐標(biāo)系不匹配導(dǎo)致的裁剪失敗這個(gè)坑我在前面提過但值得再?gòu)?qiáng)調(diào)一次。裁剪失敗最常見的原因就是坐標(biāo)系不匹配。癥狀是裁剪結(jié)果為空或者只有一小塊。排查方法是用gdalinfo查看柵格的坐標(biāo)系用ogrinfo查看矢量的坐標(biāo)系確認(rèn)兩者是否一致如果不一致用gdf.to_crs()轉(zhuǎn)換矢量或者用gdal.Warp轉(zhuǎn)換柵格。我個(gè)人的習(xí)慣是統(tǒng)一用投影坐標(biāo)系因?yàn)橥队白鴺?biāo)系下的距離和面積計(jì)算更準(zhǔn)確。6.3 并行處理加速批量任務(wù)如果數(shù)據(jù)量很大單線程處理太慢可以用Python的multiprocessing庫(kù)做并行。比如批量提取波段的時(shí)候每個(gè)文件獨(dú)立處理天然適合并行from multiprocessing import Pool import os def process_single_file(hdf_file): 處理單個(gè)文件 filename os.path.basename(hdf_file).replace(.hdf, .tif) output_path os.path.join(output_dir, filename) extract_band(hdf_file, NDVI, output_path) return filename if __name__ __main__: hdf_files glob.glob(os.path.join(input_dir, *.hdf)) with Pool(processes4) as pool: results pool.map(process_single_file, hdf_files) print(f完成 {len(results)} 個(gè)文件)進(jìn)程數(shù)根據(jù)你的CPU核心數(shù)來定一般設(shè)為核心數(shù)的一半到全部。注意GDAL在多進(jìn)程下有時(shí)候會(huì)有問題如果遇到崩潰可以改用ThreadPool或者降低進(jìn)程數(shù)。6.4 常見錯(cuò)誤與快速排查表錯(cuò)誤現(xiàn)象可能原因解決辦法裁剪結(jié)果為空坐標(biāo)系不匹配統(tǒng)一矢量與柵格坐標(biāo)系NDVI值域異常未應(yīng)用縮放因子乘以0.0001內(nèi)存溢出一次性讀取全部數(shù)據(jù)分塊處理HDF讀取失敗網(wǎng)格名稱錯(cuò)誤用gdalinfo查看子數(shù)據(jù)集路徑投影轉(zhuǎn)換后變形重采樣方法不當(dāng)連續(xù)值用雙線性分類值用最近鄰并行處理崩潰GDAL多進(jìn)程沖突降低進(jìn)程數(shù)或改用線程這張表是我在實(shí)際工作中總結(jié)出來的基本上覆蓋了80%以上的常見問題。遇到報(bào)錯(cuò)的時(shí)候先對(duì)照這張表排查能省不少時(shí)間。7. 從腳本到工作流我的實(shí)際使用體會(huì)整套流程跑通之后我把它封裝成了一個(gè)配置文件驅(qū)動(dòng)的工作流。所有的輸入路徑、輸出路徑、投影參數(shù)、分辨率、縮放因子都寫在一個(gè)YAML文件里主腳本讀取配置后自動(dòng)執(zhí)行。這樣換研究區(qū)的時(shí)候只需要改配置文件不用動(dòng)代碼。# config.yaml input_dir: D:/MODIS/MOD13A2 output_dir: D:/MODIS/output product: MOD13A2 band: NDVI scale_factor: 0.0001 fill_value: -3000 target_proj: projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm resolution: 1000 boundary: D:/data/study_area.shp這種配置驅(qū)動(dòng)的方式最大的好處是可追溯。半年后回頭看某個(gè)結(jié)果是怎么生成的翻出配置文件就一目了然不用去回憶當(dāng)時(shí)用了什么參數(shù)。另外一個(gè)小技巧是中間結(jié)果保留。很多人為了省磁盤空間處理完就把中間文件刪了。我的建議是至少保留波段提取后的GeoTIFF因?yàn)楹罄m(xù)如果要調(diào)整投影或者裁剪范圍可以直接從這一步重新跑不用從頭再來。磁盤空間現(xiàn)在很便宜但重新處理的時(shí)間成本很高。最后說一個(gè)關(guān)于ArcGIS和Python配合的心得。我現(xiàn)在的習(xí)慣是Python做批量和計(jì)算ArcGIS做交互和出圖。Python腳本負(fù)責(zé)把原始數(shù)據(jù)變成干凈的、帶地理坐標(biāo)的GeoTIFF然后導(dǎo)入ArcGIS做符號(hào)化、布局、導(dǎo)出。兩者各司其職效率最高。如果硬要用ArcPy做所有事情反而會(huì)因?yàn)锳rcGIS的Python環(huán)境限制而束手束腳。這套流程我用了三年多從最初的幾十個(gè)文件到現(xiàn)在每年處理上萬(wàn)景數(shù)據(jù)穩(wěn)定性一直很好。核心邏輯沒有變過只是根據(jù)具體需求做了些微調(diào)。如果你也在做類似的批量柵格處理希望這些經(jīng)驗(yàn)?zāi)軒湍闵僮咝澛贰?