準備
データは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)
次にこちらの記事にある筆ポリゴンのデータを使って切り抜きをしてみます
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)を揃えてあげれば良いのですが、初学の時はわからないですよね
自分もこれでかなり足踏みすることになりました。
この情報が誰かの役に立てば幸いです
Discussion