在線看seo網(wǎng)站,常見的網(wǎng)站名稱有哪些,南京今天的最新消息,wordpress緩存無法清除以下代碼基于rasterio讀取 / 寫入 TIFF、scipy實(shí)現(xiàn)插值#xff0c;可對(duì) DEM 中的 0 值#xff08;異常值#xff09;進(jìn)行鄰近插值或反距離加權(quán)插值#xff08;IDW#xff09;修復(fù)。步驟說明#xff1a;讀取 DEM 數(shù)據(jù)#xff0c;提取有效數(shù)據(jù)#xff08;非 0 值#xf…以下代碼基于rasterio讀取 / 寫入 TIFF、scipy實(shí)現(xiàn)插值可對(duì) DEM 中的 0 值異常值進(jìn)行鄰近插值或反距離加權(quán)插值IDW修復(fù)。步驟說明讀取 DEM 數(shù)據(jù)提取有效數(shù)據(jù)非 0 值和待插值的 0 值位置可選兩種插值方式鄰近插值快速、反距離加權(quán)插值更平滑插值修復(fù) 0 值區(qū)域后保存為新的 TIFF 文件增加數(shù)據(jù)范圍校驗(yàn)確保插值結(jié)果在合理區(qū)間。import numpy as np import rasterio from scipy.interpolate import NearestNDInterpolator, RBFInterpolator from rasterio.transform import from_origin def repair_dem_zero_values(input_tif, output_tif, methodnearest, power2): 修復(fù)DEM中0值異常值的函數(shù) :param input_tif: 輸入DEM文件路徑dem.tif :param output_tif: 輸出修復(fù)后的DEM文件路徑 :param method: 插值方法可選 nearest鄰近插值或 idw反距離加權(quán) :param power: IDW插值的冪次默認(rèn)2值越大近點(diǎn)權(quán)重越高 # 1. 讀取DEM數(shù)據(jù) with rasterio.open(input_tif) as src: dem_data src.read(1) # 讀取第一波段 profile src.profile # 獲取影像元數(shù)據(jù)投影、分辨率、變換等 transform src.transform # 地理變換參數(shù) nodata src.nodata # 原始nodata值若有 # 2. 標(biāo)記有效數(shù)據(jù)和待插值位置 # 有效數(shù)據(jù)非0值0為異常值且在合理范圍0-5500m valid_mask (dem_data ! 0) (dem_data 0) (dem_data 5500) # 待插值位置0值且在影像范圍內(nèi) interpolate_mask dem_data 0 # 檢查是否有需要插值的點(diǎn) if not np.any(interpolate_mask): print(無0值需要插值直接保存原始數(shù)據(jù)) with rasterio.open(output_tif, w, **profile) as dst: dst.write(dem_data, 1) return # 3. 提取有效數(shù)據(jù)的坐標(biāo)和值 # 獲取行列坐標(biāo)可轉(zhuǎn)換為地理坐標(biāo)此處用行列坐標(biāo)插值更高效 rows, cols np.meshgrid(np.arange(dem_data.shape[0]), np.arange(dem_data.shape[1]), indexingij) # 有效數(shù)據(jù)的行列坐標(biāo) valid_rows rows[valid_mask] valid_cols cols[valid_mask] valid_values dem_data[valid_mask] # 待插值的行列坐標(biāo) interp_rows rows[interpolate_mask] interp_cols cols[interpolate_mask] # 組合坐標(biāo)為N×2的數(shù)組適配scipy插值接口 valid_points np.column_stack((valid_rows, valid_cols)) interp_points np.column_stack((interp_rows, interp_cols)) # 4. 執(zhí)行插值 if method nearest: # 鄰近插值快速適合大區(qū)域 interpolator NearestNDInterpolator(valid_points, valid_values) interp_values interpolator(interp_points) elif method idw: # 反距離加權(quán)插值更平滑計(jì)算稍慢 # 用RBFInterpolator模擬IDW函數(shù)類型選inverse_multiquadric等價(jià)IDW interpolator RBFInterpolator( valid_points, valid_values, functioninverse_multiquadric, smoothing0, # 無平滑純IDW kernel_sizeNone ) # 計(jì)算IDW權(quán)重冪次power interp_values interpolator(interp_points, normpower) else: raise ValueError(method僅支持 nearest 或 idw) # 5. 修正插值結(jié)果確保在0-5500m范圍內(nèi) interp_values np.clip(interp_values, 0, 5500) # 6. 替換0值為插值結(jié)果 dem_repaired dem_data.copy() dem_repaired[interpolate_mask] interp_values # 7. 保存修復(fù)后的DEM with rasterio.open(output_tif, w, **profile) as dst: dst.write(dem_repaired, 1) print(f修復(fù)完成輸出文件{output_tif}) print(f插值點(diǎn)數(shù){len(interp_values)}) print(f有效數(shù)據(jù)占比{np.sum(valid_mask)/dem_data.size*100:.2f}%) # -------------------------- 調(diào)用示例 -------------------------- if __name__ __main__: # 輸入輸出路徑根據(jù)實(shí)際情況修改 input_dem dem.tif # 原始DEM文件 output_dem dem_repaired.tif # 修復(fù)后的DEM文件 # 方法1鄰近插值快速推薦先試用 repair_dem_zero_values(input_dem, output_dem, methodnearest) # 方法2反距離插值更平滑計(jì)算稍慢 # repair_dem_zero_values(input_dem, output_dem, methodidw, power2)依賴安裝執(zhí)行代碼前需安裝以下庫pip install rasterio scipy numpy關(guān)鍵參數(shù)說明methodnearest鄰近插值速度極快結(jié)果為最近有效像素的值適合快速修復(fù)idw反距離加權(quán)插值結(jié)果更平滑符合地形連續(xù)性適合對(duì)精度要求高的場景。power僅 IDW反距離的冪次默認(rèn) 2值越大近點(diǎn)對(duì)插值結(jié)果的影響越大通常取 1-3。clip強(qiáng)制插值結(jié)果在 0-5500m 范圍內(nèi)避免插值出現(xiàn)異常值。注意事項(xiàng)若 DEM 中 0 值區(qū)域過大無鄰近有效數(shù)據(jù)插值結(jié)果可能不準(zhǔn)確建議先檢查數(shù)據(jù)有效性地理坐標(biāo)插值若需基于經(jīng)緯度 / 投影坐標(biāo)插值可將rows/cols轉(zhuǎn)換為地理坐標(biāo)通過transform參數(shù)代碼中已預(yù)留地理變換接口大文件優(yōu)化若 DEM 文件超大如 GB 級(jí)可分塊處理參考rasterio的window功能避免內(nèi)存溢出。驗(yàn)證結(jié)果修復(fù)后可通過以下方式驗(yàn)證# 讀取修復(fù)后的數(shù)據(jù)檢查0值是否被替換 with rasterio.open(dem_repaired.tif) as src: data src.read(1) print(f修復(fù)后0值數(shù)量{np.sum(data 0)}) print(f數(shù)據(jù)范圍{np.min(data)} - {np.max(data)})