최근 옆 데이터를 사용한 위성 제품의 결손 데이터에 대한 보충

개요

안녕하세요.위성 데이터로 제작된 제품은 각양각색이지만, 때로는 관심 분야가 포함되지 않는다.이번에는 그 중 하나로 붉은 숲 속의 해면 온도 부족 문제를 해결해 봤다.

정황

적수림은 육지와 해역 한가운데 위치해 있어 MODIS 위성의 해수면 온도 제품 중 일부가 덮여 있지 않다.이를 위해 최근 해수면 온도 정보를 옆에서 끌어온 것에 따라 대상 구역 전체에 대한 해수면 온도 정보를 제작한다는 지침을 마련했다.1300km*1600km 정도의 넓은 영역을 처리했지만 최근 인접 픽셀 검색은 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 값의 상호 보완 결과