最近傍のデータを使った衛星プロダクトにおける欠損データの補完について

##概要
こんにちは。衛星データから作られるプロダクトには様々なものがありますが、時には関心領域がプロダクト中に含まれていない場合があります。今回はその一例として、マングローブ林における海面温度の欠損の問題に取り組んでみました。
##状況
マングローブ林は陸域と海域のちょうど中間に存在していることから、MODIS衛星の海面温度プロダクトでは一部がカバーされていませんでした。そのため、最近傍から海面温度情報を引っ張ってくることによって、対象領域全体についての海面温度情報を作成する方針を立てました。1300 km* 1600 kmほどのかなり広い領域の処理を行いましたが、最近傍ピクセルの探索でBrute forceな手法ではなく、より効率的なKDTreeで処理を行ったためか全体の処理時間は20~30秒ほどでした。処理の過程で距離による重みづけ等が必要な場合には異なる方法やライブラリが必要となるかと思いますが、ひとまずの欠損値埋めに使用する事が可能な方法かなと思います。

##コード

SST_interpolation.py

import rasterio
import numpy as np
import scipy.spatial as ss

#MODIS画像と対象領域のファイルは同一の解像度、同一領域に調整済み
#対象領域のファイルにおいて、対象領域では値が１、それ以外の領域は値が０となっている

#ファイルパスの設定
StudySiteFilePath="#対象領域のファイルへのパス"
MODIS_SST_FilePath="#MODIS SSTプロダクトファイルへのパス"
OutputFilePath="#（出力ファイルを格納したい場所へのパス）/出力したいファイル名.tif"

#rasterioでファイルの読み込み
StudySite=rasterio.open(StudySiteFilePath)
SST_Jan=rasterio.open(MODIS_SST_FilePath)

#行列としての読み込み
StudySite_matrix=StudySite.read(1)
SST_Jan_matrix=SST_Jan.read(1)

#ゼロより大きいSSTの場所のインデックスを取得
SST_loc=(SST_Jan_matrix>0).astype(np.uint8)
list_of_SST_points=np.where(SST_loc==1)

#対象領域のインデックスを取得
list_of_targetPoints=np.where(StudySite_matrix==1)

#最近傍の値の探索を高速化する為、KDTreeを利用
#ユークリッド距離に基づいて評価
data=[(list_of_SST_points[0][i],list_of_SST_points[1][i]) for i in range(len(list_of_SST_points[0]))]
tree=ss.KDTree(data)
target_coordinates=[(list_of_targetPoints[0][i],list_of_targetPoints[1][i]) for i in range(len(list_of_targetPoints[0]))]

#値を格納する為の行列を作成
SST_interpolated=np.zeros_like(SST_loc)

#対象領域のピクセルについて、最近傍のMODIS SSTの値を取得して行列へ格納
for coords in target_coordinates:
    d,i=tree.query(coords,p=2)
    SST_interpolated[coords[0],coords[1]]=SST_Jan_matrix[data[i][0],data[i][1]]

#GeoTiff形式で保存
with rasterio.open(OutputFilePath,'w',driver='GTiff',width=StudySite.width,height=StudySite.height,count=1,crs=StudySite.crs,transform=StudySite.transform,dtype=np.float32) as output:
    output.write(SST_interpolated,1)
    output.close()

##結果
対象領域のファイルの一部（白は対象領域、黒は非対象領域）
MODISプロダクト
最近傍の値によるSST値の補完結果