MODISの地表反射率（8日間データ）から 小縮尺グローバル衛星画像タイル（ZL2-8、ラスタタイル）を作る

はじめに

最近、小縮尺のベクトルデータを扱っていましたので、小縮尺で使える背景用衛星画像を作ろうと思い、少し作業に取り組みました。ここにその記録をメモしておきます。NASAのオープンデータから最新の衛星画像タイルを作成できますので、興味がある方は試してみてください。下図の背景画像を作りました。

また、この記事はSmart Maps Advent Calendar 2022の7日目に登録したいと思います。長くなってしまい、あまりスマートさを感じさせないかもしれませんが、ご容赦いただければ幸いです。

MODISを選んだ理由

小縮尺の衛星画像を作ろうとしたとき、空間解像度を考えると広範囲を低解像度でカバーするMODISがいいと思います。回帰日数も早いのでデータも新鮮なものがたくさんあります。

地理院地図で使っている写真レイヤでも小縮尺の衛星画像はMODISデータです。（参考： https://maps.gsi.go.jp/development/ichiran.html#modis )　このことから、MODISを使ってみようと思いました。

私の環境

• Windows 10 Enterprise
• Powershell
• Google chrome
• QGIS 3.16.6
• Docker desktop
  • unvt/nanban
    • gdal version 3.4.1
  • unvt/rgbify-node
    • gdal version 3.0.4
      なお、今回の作業レポジトリはこちらです：　https://github.com/ubukawa/worldimage

作業

Step 1. MODISのプロダクトを確認

まず、MODISのページにアクセスして、プロダクトを確認します。Land Product の Surface reflectanceは使えそうです。
https://modis.gsfc.nasa.gov/about/

バンドを確認します。可視光のバンドもきちんとあるので、RGBを4-3-1で組み合わせれば、True colorの衛星画像合成ができそうです。
https://modis.gsfc.nasa.gov/data/dataprod/mod09.php

データをもう少し選びます。MODISは２台の衛星（TerraとAqua）についているセンサなのですが、私の用途ではどちらでもいいです（TerraとかAquaって素敵な名前ですねぇ。）。雲の影響などをなるべく減らすため、8日の画像を選びます。

In the 8-day product, each surface reflectance pixel contains the best possible L2G observation during an 8-day period as selected on the basis of high observation coverage, low view angle, the absence of clouds or cloud shadow, and aerosol loading.　(MODISページの説明より、一部抄)

データシリーズ名は MOD09A1 か MYD09A1 シリーズにすればよいです。 最初に出ている Terraの方で、MOD09A1 にしようと思います。500メートル解像度で1シーンが1,200km×1,200kmです（2,400×2,400ピクセル）。

https://modis.gsfc.nasa.gov/data/dataprod/mod09.php

このMOD09A1シリーズをLP DAACのページで調べると、バージョンはv061とv006がありますが、補正に改善を行ったようなので新しいほうv061を使うことにします。詳細な説明は https://lpdaac.usgs.gov/products/mod09a1v061/ で見られます。

Step 2. データのダウンロード

2-1. EarthExplorerでのチェック

まずはいつものとおり EarthExplorer でデータの雰囲気を見てみます。NASA LPDAAC Collectionsの下にありました。下の図ではV.6をチェックしていますが、V6.1があるのでそちらを選びましょう。

ダウンロードするには、EarthExplorer（USGS）のログインだけでなく、NASAのデータダウンロードのためのログインも必要だそうです。検索すると時間が結構かかるので、プロダクトの種類に加えて観測日でフィルターを付けておくといいです。

最近3か月分くらいにしたら結構早く結果がでました。300近いファイルがあるので、EarthExplorerでなくNASAのLPDAACから直接ダウンロードしたいと思います。

2-2. NASAからのダウンロード

LP DAAC関係データのダウンロードページからデータをダウンロードします。以前の私の記事『LP DAACデータ（NASADEM） を curl で一括ダウンロード』で書いた方法と同じ要領でできましたが以下にメモしておきます。

少し探すと、 https://e4ftl01.cr.usgs.gov/MOLT/ にMODISのダウンロードページがありました。

プロダクトを選ぶと日付が選べます。8日間でよい値を選ぶデータなので8日ごとの日付が入っていますね。

いつでもいいのですが、北半球の秋、南半球の春で、なるべく最近の2022年9月14日のデータにしようかと思います。

MOD09A1.061/2022.09.14 のフォルダをみていると、観測したデータは323ファイルあるようです。hdfファイルになっていて、複数のバンドが１つのファイルになっていますね。

まず、IDとパスワードを使って.netrcファイルを作ります。テキスト形式で扱うので、環境によってはパスワードの扱いに注意しましょう。なお、私はDockerで疑似Linux環境にしてやりました。curl.exeにすればPowershellでもできるのではないかと思うのですが、うまくできなかったのでdockerを使いました。

touch ~/.netrc | chmod og-rw ~/.netrc | echo machine urs.earthdata.nasa.gov >> ~/.netrc
vi　~/.netrc　

これでできたnetrcファイルをviで編集して、NASA LP DAAC のIDとパスワードを入れます。

できたらcurlでダウンロードするためのシェルスクリプトを作ります。今回の作業レポジトリのここにおいておきます。

download.sh

curl -O -b ~/.urs_cookies -c ~/.urs_cookies -L -n https://e4ftl01.cr.usgs.gov/MOLT/MOD09A1.061/2022.09.14/MOD09A1.A2022257.h00v08.061.2022266045506.hdf;
curl -O -b ~/.urs_cookies -c ~/.urs_cookies -L -n https://e4ftl01.cr.usgs.gov/MOLT/MOD09A1.061/2022.09.14/MOD09A1.A2022257.h00v09.061.2022266045516.hdf;
curl -O -b ~/.urs_cookies -c ~/.urs_cookies -L -n https://e4ftl01.cr.usgs.gov/MOLT/MOD09A1.061/2022.09.14/MOD09A1.A2022257.h00v10.061.2022266045243.hdf;
...（続き省略）

そして実行するとダウンロードが始まりました。

./download.sh

２時間半くらいかかりましたが、全部のファイルがダウンロードできました。19GBと少しくらいのサイズでした。

Step 3. ファイルを見て処理の方法を考える

3-1. メタデータの点検

gdalinfoでファイルをチェックしてみます。メタデータを確認します。

13の画像データが入っているようですね。合成に必要なバンド４，３，１はそれぞれサブセット４，３，１でよさそうです。座標が整数なので投影法を持っているのかは若干心配です。

3-2. TIFFに展開する実験

一つのファイルを試みにtiffファイルに展開してみます。-sdsオプションをつけてサブセットを個別のファイルに展開します。13個のTIFFファイルができます。私の技術力では、HDFから必要なサブセットだけ選んで出力するコマンドを知らないので（可能かどうかもしらない）、13個それぞれ展開してから不要なものを消すような方法でやっていこうと思います。展開にそんなに時間がかかるわけではないですし。

3-3. ファイルサイズを考える

TIFFにすると一つのバンドの１つのファイルが11,276KBでした。16bit（=2byte）のデータで、2,400×2,400ピクセルなので、計算上は11,520,000 byte = 11,250 KBなので、残りの26KBくらいが投影法などのメタデータなんだと思います。

例えば300ファイルでRGB画像を作るとすると、　
　約11（MB/バンド・シーン) × 3(バンド) × 300 (シーン) ≒ 9,900 MB
でだいたい10GBです。コマンドラインで扱う分には大丈夫でしょうが、QGISで色調補正してタイル化するのはちょっと大変かもしれません。レイヤの合成処理の前に8bitにできるといいなと思います。

3-4. 位置情報の有無チェック

QGISで確認したら、位置情報も持っていました。MODISのSinusoidal Tile Gridの雰囲気が見えますね。

3-5. ストレッチを考える

16bit符号付の整数値と反射率の関係を見てみると、値を10,000で割ると反射率になっているようです。この方式は各ファイルで共通なようなので、マージするときは一括のストレッチングだけを考えればよさそうです。

処理としては、理論上０～１に存在しうる反射率の値が16ビット符号付き整数データの中では0～10,000に分布しています。最終的にRGB画像なら各バンドは8bitなので、これらのデータもgdal_calcで計算して、0~250の分布くらいの8bit整数にできるといいなと思います。8bitなら16bitの半分のデータサイズにできるかなと思います。多分、--calc="(A<=0)*0 + (A>0)*(A<=10000)*A/400 + (A>10000)*250"という感じかなと思います。

試しに１シーンでストレッチをやってみたら、0～1でなくて 0～0.35の反射率（データでいうと0～3500）を0-255の色に移すくらいがちょうどよかったです。14で割るくらいがいいかもしれません。calc="(A<=0)*0 + (A>0)*(A<=3500)*A/14 + (A>3500)*250"という感じでやってからマージすればいいかと思います。

Step 4. HDFからTIFFにエキスポートし、不要なサブレイヤを消す。

ダウンロードしてきたHDFファイルを処理していきます。まずはgdal_translateを使って、サブセットごとにTIFFファイルに出力します。バンド４，３，１以外は不要なので消します。これをbashのコマンドラインで一括でやります。

hdf2tiff.sh

echo "start";date;for f in 0_download/*.hdf; do gdal_translate ${f} -sds 1_tiff/`basename ${f} .hdf`.tif;rm 1_tiff/`basename ${f} .hdf`_02.tif;rm 1_tiff/`basename ${f} .hdf`_05.tif;rm 1_tiff/`basename ${f} .hdf`_06.tif;rm 1_tiff/`basename ${f} .hdf`_07.tif;rm 1_tiff/`basename ${f} .hdf`_08.tif;rm 1_tiff/`basename ${f} .hdf`_09.tif;rm 1_tiff/`basename ${f} .hdf`_10.tif;rm 1_tiff/`basename ${f} .hdf`_11.tif;rm 1_tiff/`basename ${f} .hdf`_12.tif;rm 1_tiff/`basename ${f} .hdf`_13.tif; done; echo end;date 

実行

./hdf2tiff.sh

動き出してからだいたい2時間半くらいで処理が終わりました。

Step 5. 各バンドを8bitで0~250までの分布に変換

gdal_calc.pyを使って、先ほど出力したTIFFファイルを8bitに変えます。0-3,500の幅（反射率でいうと0-0.35）を8bitの0-250に移し替えるので、1/14を掛け算してあげます（0以下と3,500以上はそれぞれ0と250になるようにします。）。typeオプションでByteを指定することで8bitにできます。（gdalのドキュメントだと8Intになっていました。私のgdalのバージョンが少し古いからかもしれません。）

calc8bit.sh

echo "start";date;for f in test/*.tif; do gdal_calc.py -A ${f} --calc="(A<=0)*0 + (A>0)*(A<=3500)*A/14 + (A>3500)*250" --type=Byte --outfile 2_8bit/`basename ${f}`; done; echo end;date 

この計算は30分くらいで終わりました。

Step 6. マージ

色合いのストレッチはすでに終わっているので、あとは単純にバンドごとにマージして、そのあとにRGb合成をします。

6-1. バンドごとのマージ

gdal_merge.py -o 3_merge/out-b.tif 2_8bit/*_01.tif;
gdal_merge.py -o 3_merge/out-g.tif 2_8bit/*_03.tif;
gdal_merge.py -o 3_merge/out-r.tif 2_8bit/*_04.tif;

※TIFFにマージするとき、結局空白タイルもデータサイズがかさむので、18×36=648シーン分のデータを覚悟しないといけません。1シーンはだいたい5MBですから、3GBと少しのサイズになると思います。(→結局1バンドは3.6 GBくらい、処理時間はそれぞれ3時間半くらいかかりました。)

https://modis-land.gsfc.nasa.gov/MODLAND_grid.html より　↓

6-2. ひとつのRGBファイルにマージ

gdal_merge.py -o 3_merge/testRGB2.tif -co PHOTOMETRIC=RGB -separate 3_merge/out-b.tif 3_merge/out-r.tif 3_merge/out-g.tif

Step 7. ラスタタイル化

最後にラスタタイル化をします。投影変換と、もしかしたら色調の微修正の手間があるかもしれないので、QGISで開いてQGISのツールで処理してみたいと思います。（10GB近いファイルですができるかな？？）
本当は、gdal2tilesを使おうと思ったのですが、インストールしてなかったのでQGISでやってしまいます。

7-1. QGISでファイルを開く

10GB以上のファイルを開けました。ちょっと感動しました。画像はやや暗めですね。

しかし、QGISではプロジェクトの投影変換にラスタの投影変換がうまくついていかなかったので、コマンドラインにもどってgdalwarpを使います。

7-2. 投影変換（gdalwarp）

まず試しに座標の範囲だけ指定して変換しました。

デフォルトだと12,865×12,865ピクセルになってしまったのですが、ZL8のグローバルなピクセル数（65,536×65,536）を確保したいので、解像度も指定して投影変換をしたいと思います。結構時間がかかりそうです。

gdalwarp -t_srs EPSG:3857 -te -20037508.34 -20048966.1 20037508.34 20048966.1 -ts 65536 65536 3_merge/RGB.tif rgb-3857.tif

7-3. ラスタタイル化

投影変換した画像をQGISで開いて、そのままツールのGenerate XYZ Tilesを使ってタイル化しました。gdal2tilesがあればコマンドラインでそのままでもいいと思いますが、QGISでも簡単にできあｍす。
Qiitaの執筆を間に合わせるため、前のステップで作った 12,865 ×　12,873 ピクセル のものを使いました。このため ZLは2～5にしておきました。解像度の高い画像が出力されたらZL2～8にしたいと思います。
特に気を付けることはないですが、Extentについては投影変換の時にウェブメルカトルの最大と最小を指定しているので、画像から取得すると簡単だと思います。

Runのボタンを押すとタイル作成が完了します。

Step 8. 結果の確認

ということで、小縮尺の地図に重ね合わせられる画像タイルができました。GitHubページで見てみます。
例： https://ubukawa.github.io/worldimage/map.html#3/28.98/49.81

画像が全体に暗いので、ストレッチの下端は調整した方がよいかもしれません。砂漠の明るさを考えると上限は今くらいでいいのかなと思います。おいおいやっていこうと思いますが、今回の実験はタイルができたのでOKとしたいと思います。

ベクトルタイルにも重ねることができました。

まとめ

今回は、衛星画像を使って、小縮尺でグローバルな画像タイルを作りました。既存のサービスもありますが、データサイズもそんなに大きくならないので、オープンデータから最新の画像タイルを作れると面白いのではないかなと思います。

謝辞

MODIS画像を提供してくれるNASA LP DAACに感謝します。QGISとGDALツールにお世話になりました。開発者の方に感謝します。ありがとうございました。

参考

NASAのMODISページ： https://modis.gsfc.nasa.gov/
LP DAACのページ： https://lpdaac.usgs.gov/resources/e-learning/how-access-lp-daac-data-command-line/
GDALのページ： https://gdal.org/
私の前の記事（LP DAACからのダウンロード）：　https://qiita.com/T-ubu/items/8f955f21a5603a3483b3
私の前の記事（Landsat画像のタイル化）：　https://qiita.com/T-ubu/items/d31e767e5d624b1fc715