はじめに

本記事では、衛星データの扱いとは切っても切り離せない地理空間処理するうえで便利なライブラリ「GDAL」 について紹介します。

GDALは画像やポリゴンの地理的な位置座標を確認、編集したりしてくれるライブラリです。

実際にラスター(画像形式)のデータの使用方法や例をみながら使っていきましょう。これらの加工ができるようになると自分で衛星データの解析やモデリングが自在にできるようになるのでぜひ極めたいところですね。

座標系の事前知識については以下の記事が詳しいです。

サンプルで扱う地理座標付きデータは以下に保管しています。
Github

練習する際は使ってみてください。

よく使うコマンド

まずは、よく使う以下5つのコマンドについて、使用例と合わせて紹介します。

1.GdalInfo：地理空間情報やデータ形式を表示します。
2.GdalWarp：画像を切ったり、座標変換をしたりできます。
3.Gdal Plygonize：画像をgeojson や shp ファイルなどのポリゴン(ベクター形式)に変更します。
4.Gdal Rasterize：今度はポリゴンからtiff ファイルなどのラスタ(画像形式)に変換します。
5.Gdal Transrate：解像度変更やデータ形式などの変更をします。

1.GdalInfo

https://gdal.org/programs/gdalinfo.html
扱うデータについて、地理空間やデータ形式などの情報を表示します。

使用方法

$ gdalinfo 詳細が知りたいファイル名

使用例①

座標系や欠損データの扱いについて出力します。サイズや緯度経度などの地理空間情報を知りたい時に使用しましょう。

$ gdalinfo sample/sorabatake_geocode_image.tif

Driver: GTiff/GeoTIFF
Files: sample/sorabatake_geocode_image.tif
Size is 353, 173
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - 85°S to 85°N"],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (15556490.661657888442278,4257688.364789782091975)
Pixel Size = (7.399196882220372,-7.399196882220372)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Band 1 Block=353x23 Type=Byte, ColorInterp=Palette
  NoData Value=251

使用例②

出力結果をJson形式で保管します。地理空間の内部情報を出力して見える形にしたい時に使えます。

$ gdalinfo -json sample/sorabatake_geocode_image.tif > output/sorabatake_geocode_image.json

2.GdalWarp

https://gdal.org/programs/gdalwarp.html
画像を切り抜きや座標の変換ができます。

使用方法

$ gdalwarp -te 欲しい画像の座標 -t_srs 元の座標系 -te_srs 出力座標系 入力ファイル 出力ファイル

使用例①

皇居周辺だけ切り抜く場合など、画像を地理情報を元にくり抜きたい時に使用できます。

$ gdalwarp -te 139.746543 35.677909 139.756337 35.686451 -t_srs EPSG:4326 -te_srs EPSG:4326 sample/sorabatake_geocode_image.tif output/sorabatake_cut.tif

因みにこの時のWGS84(EPSG:4326)の座標では
右下が (139.7465430000000026,35.6779089999999997)
左上が (139.7563370000000020,35.6864509999999981)
の位置です。

使用例②

座標を変換する場合など、投影座標系を変えたい時に使いましょう。微妙ですがリサンプリングされて変化しています。

$ gdalwarp -t_srs EPSG:4329 -te_srs EPSG:4326 sample/sorabatake_geocode_image.tif output/sorabatake_EPSG4329.tif

使用例③

欠損データがある場合などは、その場所の値に0を入れて可視化することができます。

$ gdalwarp -t_srs EPSG:4326 -te_srs EPSG:4326 -srcnodata 0 sample/sorabatake_geocode_image.tif output/sorabatake_nan0.tif

黒い部分が 0 の値になっています

3.Gdal Plygonize

https://gdal.org/programs/gdal_polygonize.html
画像をgeojson や shp ファイルなどのポリゴン(ベクター形式)に変更します。

使用方法

$ gdal_polygonize.py -f 出力形式 -nomask マスクも出力するか 入力ファイル 出力ファイル

使用例①

特定の値のみのベクトルデータに変換したい時に使用しましょう。画像の値がある部分をポリゴンの geojson として出力します。

$ gdal_polygonize.py -f GeoJson sample/sorabatake_geocode_image.tif output/sorabatake_polygon.geojson

使用例②

画像全体での全ての値をベクトルデータにしたい時に使いましょう。欠損データが入っていない部分も含めてポリゴンとして出力されます。

$ gdal_polygonize.py -f GeoJson -nomask sample/sorabatake_geocode_image.tif output/sorabatake_polygon_nomask.geojson

4.Gdal Rasterize

https://gdal.org/programs/gdal_rasterize.html
今度はポリゴンからtiff ファイルなどのラスタ(画像形式)に変換します

使用方法

$ gdal_rasterize -burn 画像に代入する値 -te 四隅の座標 -a_srs 出力座標系 -ot 型 -ts 解像度 -a_nodata 欠損値の値 入力ファイル 出力ファイル

使用例①

ポリゴンから光学衛星のような画像を取得したい時に使用しましょう。ポリゴンファイルから3チャンネルで256 x 128の解像度のRGB画像を作成します。

$ gdal_rasterize -burn 255 -burn 255 -burn 255 -a_srs EPSG:3857 -ot Byte -a_nodata 0 -ts 256 128 output/sorabatake_polygon.geojson output/sorabatake_rasterize.tif

使用例②

ポリゴンからSAR衛星のような画像を取得したい時に使用しましょう。
ポリゴンファイルから1チャンネルで256 x 256の解像度の32Bit画像を作成します。

$ gdal_rasterize -burn 1  -a_srs EPSG:3857 -ot Float32 -a_nodata 0 -ts 256 256 output/sorabatake_polygon.geojson output/sorabatake_rasterize_fp32.tif

5.Gdal Transrate

https://gdal.org/programs/gdal_translate.html
解像度変更やデータ形式などの変更をします。

使用方法

$ gdal_translate -of 出力形式 -a_srs 出力座標系 -ot 型 -tr 解像度 入力ファイル 出力ファイル

使用例①

一般的な画像から浮動点小数の細かい画像へ変換したい時に使用しましょう。
8Bit 画像から 32Bit画像のgeotiff に変換します。

$ gdal_translate -of GTiff -ot Float32 sample/sorabatake_geocode_image.tif output/sorabatake_fp32.tif

使用例②

異なる解像度の画像の解像度を調節したい時に使用しましょう。地理座標空間で解像度が64ピクセル四方の画像を生成します。

$ gdal_translate -tr 64 64  sample/sorabatake_geocode_image.tif output/sorabatake_resolution64.tif

以上、5つのコマンドとそれぞれの使用例を紹介しました。より詳細な設定が必要になった際は公式ドキュメントをご覧ください。

環境構築について

おすすめは Docker です。コマンド１つで誰でも環境構築できます。仮想環境なので開発環境を汚しません。手軽さでは Anaconda ですが社員規模などでの商用利用などには注意してください。

Docker

https://www.docker.com/

Docker とは仮想環境を作成したり、管理してくれるツールです。

Gdal に加えて JupyterNook 環境も設定しています。
Tellus 開発環境のようにすぐに衛星データの分析ができる方が便利だからです。

Dockerfile

FROM osgeo/gdal:ubuntu-full-3.4.3


ENV DEBIAN_FRONTEND=noninteractive
SHELL ["/bin/bash", "-c"]


RUN apt-get update && apt-get upgrade -y
RUN apt-get update --fix-missing && apt-get install -y --no-install-recommends\   
       build-essential \
       libfreetype6-dev \
       libpng-dev \
       libzmq3-dev \
       python3-dev \
       libspatialindex-dev \
       gdal-bin \
       libgdal-dev \
       libsm6 \
       pip \
       && \
   apt-get clean && \
   rm -rf /var/lib/apt/lists/*


RUN pip install -U pip
RUN pip install --upgrade pip setuptools
COPY env/requirements.txt .
RUN  pip install -r requirements.txt
ENV PYTHONPATH $PYTHONPATH:/workspace


RUN mkdir /workspace
WORKDIR /workspace


ENTRYPOINT ["jupyter-lab", "--NotebookApp.token=''", "--allow-root", "--no-browser", "--port", "8888", "--ip=0.0.0.0"]

compose.yml

services:
 gdal-python:
   build:
     context: .
     dockerfile: env/Dockerfile
   ports:
     - 8888:8888
   volumes:
     - ./src/:/workspace
   working_dir: /workspace

Docker Compose で管理することで再現も楽にしています。

$ git clone https://github.com/sorabatake/article_32245_gdal.git
$ docker compose up -d

そして、リンク: http://localhost:8888/lab にアクセスします。

Anaconda

https://www.anaconda.com/download
Anaconda は Python の仮想環境を管理してくれるサービスアプリです。

$ conda install gdal

# ミラー補完
$ conda install -c conda-forge gdal
$ conda install gdal==3.4.3

$ pip install -r env/requirements.txt
$ jupyter-lab --NotebookApp.token='' --allow-root --port 8888 --ip=0.0.0.0

逆引きリスト

●データがよくわからない
〇Gdalinfo 使用例①
●２枚の衛星画像の同じ場所を切り抜きたい
〇共通の切り抜く場所の四隅の座標を取得して Gdalwarp 使用例① をそれぞれの画像で切り抜く
●位置は同じだが解像度の違う衛星画像を同じ解像度に揃える
〇Gdal Transrate 使用例②
●ポリゴンの教師データを衛星画像と同じ画像形式にする
〇衛星の画像サイズを取得して Gdal Rasterize 使用例① , Gdal Rasterize 使用例②
●衛星画像を解析データのエリア情報として切り抜きたい
〇Gdal Plygonize 使用例①
●geojsonをtiff 画像のデータ形式に変換する
〇Gdal Rasterize 使用例①, Gdal Rasterize 使用例②
●tiff画像をgeojsonのデータ形式に変換する
〇Gdal Polygonize 使用例①

使用例と環境構築のコード

Github にて公開していますので、こちらをご覧ください。

記事の内容を実践してみて、エラーが出た場合は、こちらのGoogleフォームにエラー内容のログを全てと改修したいことをリクエストいただけますと幸いです。
⇒エラー発見&改善リクエストはこちら

まとめ

GDALをマスターして衛星画像の分析や可視化、予測モデルの検討などをこれからの宇宙活用の発展の１歩を一緒に担って行きましょう！