はじめに

　衛星画像を処理する際には、解像度の違いが問題となります。例えば異なる解像度（ピクセル数）の衛星プロダクトをNumpy arrayとして読み込んだ場合、それらは要素数異なる行列として読み込まれるため、行列についての多くの演算を行う事ができません。そのため、複数の衛星画像プロダクトを用いるのに先立ち解像度や投影を統一しておく必要があります。この課題を解決する為の処理コードを書きましたので、簡単にご紹介したいと思います。

コードができる事

①基準となる画像に合わせてリサンプリング及びクリップ処理を行う
②任意の領域に入力画像をクリップする
③任意の解像度へリサンプリングする

前提条件

　以下の状態にフォルダが構成されていることが前提となっています。例としましては、Target folderにプロダクトの名前、Sub folderにある年、もしくは月の画像を入れてある様な状態です。投影法は主に使われるWGS84とUTMに対応しています。

コードと使い方

ClipResample.py

import rasterio 
from rasterio.mask import mask
from shapely.geometry import box
import geopandas as gpd
from fiona.crs import from_epsg
import json
import pycrs
import numpy as np
from rasterio.enums import Resampling
import os, glob, shutil
from tqdm import tqdm

class ClipResampleRaster:
    
    def __init__(self,pathToTargetFolder,pathToRefImage=None):
        self.pathToTargetFolder=pathToTargetFolder 
        if pathToRefImage==None:
            pass
        else:
            self.pathToRefImage=pathToRefImage
    
    def boundaryBox_ref(self):
        ref_image=rasterio.open(self.pathToRefImage)
        bounds=ref_image.bounds
        bbox = box(bounds[0], bounds[1], bounds[2], bounds[3])
        geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
        geo = geo.to_crs(crs=ref_image.crs)
        coords=[json.loads(geo.to_json())['features'][0]['geometry']]
        return coords
    
    def outputRasterMeta_ref(self):
        ref_image=rasterio.open(self.pathToRefImage)
        self.height=ref_image.meta["height"]
        self.width=ref_image.meta["width"]
        self.count=ref_image.meta["count"]
        self.crs=ref_image.meta["crs"]
        self.transform=ref_image.meta["transform"]
        
    def boundaryBox(self,boxCoordinates=None):
        pathToIndividualFolder=os.listdir(self.pathToTargetFolder)[0]
        fileName=os.listdir(self.pathToTargetFolder+f"/{pathToIndividualFolder}")[0]
        path=self.pathToTargetFolder+f"/{pathToIndividualFolder}"+f"/{fileName}"
        src_image=rasterio.open(path) 
        
        if boxCoordinates==None:
            boxCoordinates=list(src_image.bounds)
            bbox=box(boxCoordinates[0],boxCoordinates[1],boxCoordinates[2],boxCoordinates[3])
        else:
            bbox=box(boxCoordinates[0],boxCoordinates[1],boxCoordinates[2],boxCoordinates[3])
        geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=src_image.crs)
        coords=[json.loads(geo.to_json())['features'][0]['geometry']]
        self.mapCenter=((boxCoordinates[0]+boxCoordinates[2])/2,(boxCoordinates[1]+boxCoordinates[3])/2)
        
        out_img, out_transform = rasterio.mask.mask(src_image, coords, crop=True,filled=False)
        self.width=out_img.shape[2]
        self.height=out_img.shape[1]
        self.transform=out_transform
        return coords
    
    def outputRasterMeta(self,resolutionByMeter,resolutionByDegree=None):
        pathToIndividualFolder=os.listdir(self.pathToTargetFolder)[0]
        fileName=os.listdir(self.pathToTargetFolder+f"/{pathToIndividualFolder}")[0]
        path=self.pathToTargetFolder+f"/{pathToIndividualFolder}"+f"/{fileName}"
        src_image=rasterio.open(path)
        if (src_image.meta["crs"]==from_epsg(4326))==True:
            d=resolutionByMeter/1000
            r=6371
            lat=np.radians(self.mapCenter[1])
            resolution=abs(np.degrees(2*np.arcsin(np.sin(d/(2*r)))))
        elif (resolutionByDegree!=None)==True:
            resolution=resolutionByDegree
        else:
            resolution=resolutionByMeter
        
        src_affine=src_image.meta["transform"]
        src_width=src_image.meta["width"]
        src_height=src_image.meta["height"]
        target_transform,target_width,target_height=rasterio.warp.aligned_target(self.transform,self.width,self.height,resolution)
        
        self.height=target_height
        self.width=target_width
        self.count=src_image.meta["count"]
        self.crs=src_image.meta["crs"]
        self.transform=target_transform
        
    def getTargetFiles(self):
        targetFolderNames=os.listdir(self.pathToTargetFolder)
        targetFiles={}
        for name in targetFolderNames:
            targetFiles[name]=os.listdir(self.pathToTargetFolder+f"/{name}")
        self.targetFiles=targetFiles
    
    def resampleRaster(self,path,fileName,resamplingMethod):
        if resamplingMethod=="nearest":
            resamplingMethod_func=Resampling.nearest
        elif resamplingMethod=="bilinear":
            resamplingMethod_func=Resampling.bilinear
        else:
            print("Please specify the resampling method either nearest or bilinear.")
            
        name=fileName.replace(".tif","")
        with rasterio.open(path+f"/{name}_clipped.tif") as dataset:
            dataset=dataset.read(out_shape=(self.height,self.width),resampling=resamplingMethod_func)
        return dataset
    
    def writeClippedRaster(self,path,fileName):
        data=rasterio.open(path+f"/{fileName}")
        out_img, out_transform = rasterio.mask.mask(data, coords, crop=True,filled=False)
        name=fileName.replace(".tif","")
        dtype=out_img.dtype
        with rasterio.open(path+f"/{name}_clipped.tif",
                           'w',
                           driver='GTiff',
                           width=out_img.shape[2],
                           height=out_img.shape[1],
                           count=data.count,
                           crs=data.crs,
                           transform=out_transform,
                           dtype=dtype) as output:
            for i in range(data.count):
                output.write(out_img[i],i+1)
            output.close()
        data.close()
    
    def writeResampledRaster(self,path,fileName,raster):
        dtype=raster.dtype
        name=fileName.replace(".tif","")
        with rasterio.open(path+f"/{name}_resampled.tif",
                           'w',
                           driver='GTiff',
                           width=self.width,
                           height=self.height,
                           count=self.count,
                           crs=self.crs,
                           transform=self.transform,
                           dtype=dtype) as output:
            for i in range(self.count):
                output.write(raster[i],i+1)
            output.close()
    
    def writeOutputRaster(self,resamplingMethod,suffixOfSavedFileFolder):
        targetFiles=self.targetFiles
        writeClippedRaster=self. writeClippedRaster
        resampleRaster=self.resampleRaster
        writeResampledRaster=self.writeResampledRaster
        for folder in tqdm(targetFiles.keys()):
            path=self.pathToTargetFolder+f"/{folder}"
            path_to_savedFileFolder=self.pathToTargetFolder+f"/{folder}_{suffixOfSavedFileFolder}"
            os.mkdir(path_to_savedFileFolder)
            for fileName in targetFiles[f"{folder}"]:
                writeClippedRaster(path,fileName)
                dataset=resampleRaster(path,fileName,resamplingMethod)
                writeResampledRaster(path_to_savedFileFolder,fileName,dataset)
    
    def cleanUp(self):
        targetFiles=self.targetFiles
        for folder in tqdm(targetFiles.keys()):
            path=self.pathToTargetFolder+f"/{folder}"
            shutil.rmtree(path)
        print("The clean up process is completed!")

・ClipResampleRaster.boudaryBox_ref()
基準となる画像（Reference image）を用いる場合、その領域を計算する
・ClipResampleRaster.outputRasterMeta_ref()
基準となる画像に合わせる為のメタデータを作成する
・ClipResampleRaster.boudaryBox(boxCoordinates=None)
特定範囲で切り抜きたい場合には、関心領域のboxCoordinatesに緯度経度（長方形）を与える事で、適切な形式へと変換を行った境界の緯度経度の情報を計算する。UTMの場合には、緯度経度を自動的にUTM座標へと変換する。デフォルトでは切り抜きを行わない。
・ClipResampleRaster.outputRasterMeta(resolutionByMeter,resolutionByDegree=None)
基準となる画像を用いない場合に、出力画像に必要なメタデータを計算する。解像度をm単位で与えれば、WGS84の場合にはinverse haversine formulaによって対象領域中央におけるdegree単位での解像度へと変換する。UTMならばそのまま用いる。resolutionByDegreeにより、解像度をdegreeで与える事もできる。
・ClipResampleRaster.getTargetFiles()
処理を行う対象のファイルのリストを取得する。
・ClipResampleRaster.writeOutputRaster(resamplingMethod,suffixOfSavedFileFolder)
出力画像を書き込む。resamplingMethodでnearest neighborもしくはbilinearを指定できる。suffixOfSavedFileFolderには、処理後の画像が入るファイル名を指定できる。
・ClipResampleRaster.cleanUp()
実行することで、処理前の画像や処理で生じる中間ファイル（_clippedの付くデータ等）を削除する。

使い方１
基準となる画像（"MCD12Q1.A2001001.IGBP.Buryat.geotiff.tif"）に合わせて、ターゲット画像の入ったフォルダ（"/content/drive/MyDrive/Colab Notebooks/satellite_data/MOD11A2.2001-2020.LST"）内にある全てのデータでclipとresamplingを適用する。

convertRaster_to_referennce.py

# ターゲット画像の入ったフォルダへのパス
targetFolderPath="/content/drive/MyDrive/Colab Notebooks/satellite_data/MOD11A2.2001-2020.LST" 
# 基準となる画像へのパス
referenceImagePath="/content/drive/MyDrive/Colab Notebooks/satellite_data/MCD12Q1/MCD12Q1.A2001001.IGBP.Buryat.geotiff.tif" 

database=ClipResampleRaster(targetFolderPath,referenceImagePath) 
coords=database.boundaryBox_ref() #基準となる画像の境界の緯度経度を取得
database.getTargetFiles() #ターゲット画像のリストを取得
database.outputRasterMeta_ref() #基準となる画像からメタデータを取得
database.writeOutputRaster("nearest","outputFileFolder") #出力画像を書き込む

実行前のフォルダ構造

実行後のフォルダ構造（outputFileFolderには、基準となる画像に合わせてClipとResamplingが適用された画像が入っている）

使い方2
任意の解像度に変更する（Clipは行わない）

ChangeResolution.py

targetFolderPath="(Target Folderへのパス)"
database=ClipResampleRaster(targetFolderPath)
coords=database.boundaryBox() 
database.getTargetFiles()
database.outputRasterMeta(2000) #解像度を2000 mに設定
database.writeOutputRaster("nearest","outputFileFolder")

使い方3
解像度を変更し、対象領域も同時にClipする

ChangeResolutionAndROI.py

bounds=[100,52,105,55] #対象領域を長方形の境界の緯度経度で与える
ROI=[xmin, ymin, xmax, ymax]

targetFolderPath="（Target folderへのパス）"
database=ClipResampleRaster(targetFolderPath)
coords=database.boundaryBox(bounds) 
database.getTargetFiles()
database.outputRasterMeta(2000)
database.writeOutputRaster("nearest","outputFileFolder")