👏

Sentinel-1のSAR画像をmask処理で切りとる(筆ポリゴンで切り抜き)

2023/07/27に公開

準備

データはopen access hubなどから適宜ダウンロードしてきてください

今回は4つの

  • S1A_IW_GRDH_1SDV

で始まるプロダクトを使います

gdalが必要でインストールが面倒なのでdockerのgdalイメージを利用します

version: "3"
services:
  gdal-env:
    build: .
    ports:
      - "8888:8888"
    volumes:
      - ./:/root/
    command: tail -f /dev/null
FROM osgeo/gdal:ubuntu-full-3.5.0

RUN apt-get update -y && apt-get upgrade -y
RUN apt-get install build-essential \
    wget git  -y

RUN apt-get install libssl-dev libbz2-dev libreadline-dev \
    libsqlite3-dev python3-tk tk-dev python-tk \
    libfreetype6-dev libffi-dev liblzma-dev -y

ENV HOME  /root

WORKDIR /root

コンテナに入ったあと以下のコマンドでnotebookを起動するとgdalインストール済みのnotebookを起動できる

apt update
apt install -y python3-pip
pip install jupyterlab
jupyter-lab --ip 0.0.0.0 --allow-root --no-browser -b localhost

前処理

ファイルの保存先と、取得するデータの形式を指定

sentinel_1_dir_list = [
    '../sentinel-1/S1A_IW_GRDH_1SDV_20230408T211724_20230408T211749_048010_05C544_CFF6.SAFE/',
    '../sentinel-1/S1A_IW_GRDH_1SDV_20230502T211725_20230502T211750_048360_05D106_7A7C.SAFE/',
    '../sentinel-1/S1A_IW_GRDH_1SDV_20230607T211726_20230607T211751_048885_05E0F8_697F.SAFE/',
    '../sentinel-1/S1A_IW_GRDH_1SDV_20230701T211728_20230701T211753_049235_05EB9D_72F3.SAFE/'
]

file_type_list = ["measurement/*[v][v]*.tiff", "measurement/*[v][h]*.tiff"]

自前のポリゴンを作成する

bbox = box(minx=130.21971622772818, miny=33.071356339303456, maxx=130.50652945686997, maxy=33.30119764666647)
geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs='EPSG:4326')
geo.crs

関数の用意

def warp_img(data_dir):
    for dir_path in data_dir:
        for file_type in file_type_list:
            path_c = glob.glob(dir_path + file_type)
            path_b = path_c[0]
            print(f"warp " + dir_path + file_type)
            gdal.Warp(path_b, path_b, dstSRS="EPSG:4326")


def crop_sentinel_1(data_dir, polygon):
    save_dir = './cropped_sentinel_1/'
    cropped_img_dir = Path(save_dir)
    if cropped_img_dir.exists():
        shutil.rmtree(save_dir)
    
    cropped_img_dir.mkdir()
    
    for dir_path in data_dir:
        for file_type in file_type_list:
            path_c = glob.glob(dir_path + file_type)
            path_b = path_c[0]
    
            date_info = path_b.split('/')[-1].split('-')[4][:8]
            data_type = path_b.split('/')[-1].split('-')[3]
            cropped_file_name = date_info[0:4] + '_' + date_info[4:8] + '_' + data_type
            print(f"making {cropped_file_name}")
            
            epsg = 'epsg:4326'
            nReserve_proj = polygon.to_crs({'init': epsg})
            
            with rasterio.open(path_b) as src:
                out_image, out_transform = rasterio.mask.mask(src, nReserve_proj.geometry,crop=True)
                out_meta = src.meta.copy()
                out_meta.update({"driver": "GTiff",
                                 "height": out_image.shape[1],
                                 "width": out_image.shape[2],
                                 "transform": out_transform})
                
            with rasterio.open(f'./cropped_sentinel_1/{cropped_file_name}.tif', "w", **out_meta) as dest:
                dest.write(out_image)

def display_cropped_img(img_dir, polygon):
    crop_sentinel_1(sentinel_1_dir_list, polygon)
    for img_path in sorted(Path(img_dir).glob('**/*')):
        plt.figure(figsize=(5, 5))
        img = rasterio.open(img_path).read()
        plt.imshow(np.clip(img[0], 0, 1000)/1000)
        plt.title(str(img_path).split('/')[1].split('.')[0])
	
# warp_img(sentinel_1_dir_list)

初回のみwarp_img関数を実行してSRSをEPSG:4326に変換する(結構時間がかかる)

これでポリゴンと画像のSRSが揃ったのでレイヤーを重ねてmask処理ができるようになりました

これで可視化のための関数をつかって確認してみます
(簡単のためにnp.clip(img[0], 0, 1000)/1000で数字を丸めていますが、本来だったらもう少し丁寧な正規化をした方が良いかと思います。今回は切り取ったものを見るだけなのでこうしてます)

display_cropped_img('./cropped_sentinel_1/', geo)

次にこちらの記事にある筆ポリゴンのデータを使って切り抜きをしてみます

https://zenn.dev/fusic/articles/6807ccf192eda0

i = 24
sample = fude_polygon_copy_sorted_desc[i:i+1]
display_cropped_img('./cropped_sentinel_1/', sample)

とまあこんな感じでsentinel-1のmask処理をやってみました。
そのままやろうと思うとポリゴンで切れないという事例になりますよねw

知ってさえいればwarpでSRS(CRS)を揃えてあげれば良いのですが、初学の時はわからないですよね
自分もこれでかなり足踏みすることになりました。

この情報が誰かの役に立てば幸いです

Fusic 技術ブログ

Discussion