オーストラリア　森林火災の衛星画像のタイムラプスを作ってみた．

## 1. はじめに．
　これまで，欧州の人工衛星Sentinelシリーズが提供する無料の最新の観測画像の入手方法や，入手した画像の応用編としてタイムラプスの作成方法を紹介しました．
　
　[無料で最新の衛星画像を入手する方法．](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)
  [人工衛星（Sentinel-2)の観測画像をAPIを使って自動取得してみた．](https://qiita.com/nigo1973/items/e905959b6167ee1074d9)
  
  今回は，前回の紹介した[無料で人工衛星画像のタイムラプスの作り方．](https://qiita.com/nigo1973/items/28c873a21fc155a511bf)を修正し，１月単位ではなく観測毎のタイムプラスのgifアニメーションの作成方法を紹介します．このとき，一般的な人工衛星の観測画像ではなく，その応用編として火災にフォーカスした赤外波長の観測画像を用いた画像処理についても紹介します．
　ここで作成したGifアニメーションの例（解像度を落としています）がこちらです．

![Australia_wildfire_new2 (2).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/269d6a70-3337-c168-6d20-ec83a33270db.gif)

 また，ここで紹介する記事と関係するものを宙畑にもアップしています．よろしければ，一緒に見てください．

　[【コード付き】Sentinel-2の衛星画像から指定した範囲のタイムラプスを作成する](https://sorabatake.jp/10366/)

前回の記事と同じく， gitにも今回のコードを参考にアップしています．個人での利用の範囲で楽しんでください．

　[RS_Sentinel2_timelapse_wildfire](https://github.com/shirokawakita/RS_Sentinel2_timelapse_wildfire)

その後，豪雨によってオーストラリアの火災は鎮火された，と報告がありました．ただ，豪雨によって洪水になったのはなんとも．

[豪シドニー、３０年ぶりの豪雨で洪水　森林火災は鎮火](https://www.cnn.co.jp/world/35149229.html)
 
##2.火災箇所の人工衛星の観測画像のタイムプラスのつくりかた．

###2.1 人工衛星Sentinel-2について
Sentinelシリーズは，欧州宇宙機関が開発・運用する地球観測シリーズで，すでに１から５まで打上がり運用されています．
今回は，デジカメ画像と同じく光学観測画像を撮像するSentinel-2の衛星画像の取得から，画像処理を紹介します．

![sentinel2.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/45195a8b-2260-fe4a-29e6-2eb924adfe6f.jpeg)
Sentinel-2 ©ESA 2000-2018.

Sentinel-2について，その観測計画や観測画像の詳細については，[無料で最新の衛星画像を入手する方法．](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)をご参考にしてください．

Sentinel-2は１２の波長の異なる波長の画像を撮像しています．

![スクリーンショット 2020-02-15 17.45.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/d7389bf3-9b76-181b-03a2-a8326289a24e.png)
(Credit: European Space Agency)

一般的な衛星画像（True color)であれば，バンドの４，３，２をR,G,Bとすることで作ることができます．

今回は，森林火災にフォーカスし，赤外波長の観測画像をもちいて，火災の可視化を試してみます．
今回は，以下の記事を参考にしています．

[How to use open source satellite data for your investigative reportingHow to use open source satellite data for your investigative reporting](https://towardsdatascience.com/how-to-use-open-source-satellite-data-for-your-investigative-reporting-d662cb1f9f90)

Copernicusが提供する衛星画像をブラウザにて閲覧するアプリケーションの”Eo browser”でも異なる波長の観測画像を用いた処理ができますので，試してみてください．

![スクリーンショット 2020-02-15 17.52.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/56533f6f-b1d1-b60f-9fcd-f12e87e6be77.png)
[オーストラリアの森林火災のSentinel-2観測画像（２０１９年１２月３１日撮像）](https://apps.sentinel-hub.com/eo-browser/?lat=-34.7885&lng=150.3009&zoom=10&time=2019-12-31&preset=CUSTOM&datasource=Sentinel-2%20L2A&layers=B12,B11,B8A&evalscript=cmV0dXJuIFtCMTIqMi41LEIxMSoyLjUsQjhBKjIuNV07)

###2.2 観測毎のタイムプラス画像の作成方法
では，Sentinel-2の観測画像のタイムプラスを作っていきます．
Sentinel-２は，世界中の多くの地域において，[同地域を５日毎に観測しています](https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/revisit-coverage)．
![スクリーンショット 2020-02-15 17.55.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/69fbc0cc-a7b7-ac55-7c0b-7700efd57ae7.png)
(credit:European Space Agency)

そこで，前回のタイムプラスのコードにてしている観測期間に関する部分を以下に修正しました．

```python
#L2Aデータ取得のため，２０１９年１月以降を開始日とすること．
Begin_date = '20191201'
End_date = '2020210'

m = 5 # 5日間毎の画像を取得する．

#開始日と終了日の日差の数を求める．
import datetime
begin_data_s = datetime.datetime.strptime(Begin_date, '%Y%m%d') 
end_data_s = datetime.datetime.strptime(End_date, '%Y%m%d')
d_day = end_data_s - begin_data_s
n =d_day.days//m

#Sentinel-2の観測画像の取得．
for i in range(n):
    if i < 1:
        Begin_date = begin_data_s
        End_date = begin_data_s + datetime.timedelta(days=m)
        begin= Begin_date.strftime('%Y%m%d')
        end = End_date.strftime('%Y%m%d')
    else:
        Begin_date = End_date
        End_date = Begin_date + datetime.timedelta(days=m)
        begin= Begin_date.strftime('%Y%m%d')
        end = End_date.strftime('%Y%m%d')
    
    Sentinel2_get() #Sentinel-2の観測画像をAPIで取得する関数．

```

###2.3 赤外波長の観測画像の作り方．

取得したSentinel-2の画像データから，赤外波長の観測画像を用いた火災にフォーカスした解析画像を作成します．

```python
    path = str(product_title) + '.SAFE/GRANULE'
    files = os.listdir(path)
    pathA = str(product_title) + '.SAFE/GRANULE/' + str(files[0])
    files2 = os.listdir(pathA)
    pathB = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R20m'
    files3 = os.listdir(pathB)
    
    
    path_b12 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R20m/' +str(files3[0][0:23] +'B12_20m.jp2')
    path_b11 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R20m/' +str(files3[0][0:23] +'B11_20m.jp2')
    path_b8 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R20m/' +str(files3[0][0:23] +'B8A_20m.jp2')
    
    # Band４，３，２をR,G,Bとして開く．
    b8 = rio.open(path_b8)
    b12 = rio.open(path_b12)
    b11 = rio.open(path_b11)

```
上記にて取得した観測画像から，　バンド１２，１１と８Aの観測画像を指定します．
このとき，２．１にも書きましたが，これらの観測画像の分解能は２０mになるため，そちらのアドレスをしています．

そして，取得したB１２，B１１，B8Aの観測画像を，それぞれR, G, Bとしてカラー合成を行います．

```python
    #RGBカラー合成（GeoTiffファイルの出力）
    with rio.open(str(object_name) +'.tiff','w',driver='Gtiff', width=b12.width, height=b12.height, 
              count=3,crs=b12.crs,transform=b12.transform, dtype=b12.dtypes[0]) as rgb:
        rgb.write(b12.read(1),1) 
        rgb.write(b11.read(1),2) 
        rgb.write(b8.read(1),3) 
        rgb.close()

```
ここで作成した画像はtiffファイルですので，作成後にQGISのようなGISソフトで解析したいときは，この画像をつかいます．
その後，関心域のみのtiff画像，　tiff画像のjpeg画像への変換をへて，以下のような画像を作成します．

![20191231Masked_Australia_wildfire_new2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/c8fec2c3-63b9-cfdc-dd10-f680554179a5.jpeg)

関心域の設定が正方形でなかったため，黒い部分が周囲にできています．
そこに観測日の記載がかぶってしまったので，ちょっと見づらいですね．　観測日の記載箇所については，コードにて位置をかえるなど，試してみてください．

希望期間のそれぞれの観測画像の取得からカラー合成処理を行った後，観測日が古い画像の順番にgifアニメーションを作成します．

##3. おわりに
 今回は，前回の応用編として，Sentinel-2の赤外波長の観測画像を用いて，　森林火災にフォーカスしたタイムラプスのアニメーションを作成しました．　赤外以外にも，他の異なる波長をもちいることで，森林や水源にフォーカスした画像をつくることもできます．
 [課題に応じて変幻自在？ 衛星データをブレンドして見えるモノ・コト #マンガでわかる衛星データ](https://sorabatake.jp/5192/)

上記の記事の画像はどういった波長の組み合わせでつくるのか，　WebアプリケーションのEo browserにも例がありますので，　そちらも参考にしてみてください．

つぎは，SAR衛星であるSentinel-１の観測画像の取得，およびその解析例について紹介してみようと思います．

ご意見や質問などありましたら，コメントいだけると嬉しいです．

## 追記（２０２０．２．１７）
　新型コロナウィルスの感染拡大が続く中国武漢に新設された病院の建設模様を人工衛星の観察画像を求めました．

![china_wuhan_hospital_building.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/329e4a62-8228-a51d-e0f2-39aae1ab0b36.gif)

参考にしたのは以下の記事です．特に，New York Timesの写真は，その表現力の高さ，訴える力を感じました．
新設された病院によって感染された方々の治療が進み，感染拡大に対して大きな効果がでることを期待します．

[Wuhan Reacts To Coronavirus, The Taal Volcano Erupts And More](https://www.planet.com/pulse/wuhan-reacts-to-coronavirus-the-taal-volcano-erupts-and-more/)
[China Pledged to Build a New Hospital in 10 Days. It’s Close.](https://www.nytimes.com/2020/02/03/world/asia/coronavirus-wuhan-hospital.html)
## 参考記事
　[無料で最新の衛星画像を入手する方法．](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)
  [人工衛星（Sentinel-2)の観測画像をAPIを使って自動取得してみた．](https://qiita.com/nigo1973/items/e905959b6167ee1074d9)
  [【コード付き】Sentinel-2の衛星画像から指定した範囲のタイムラプスを作成する](https://sorabatake.jp/10366/)
  [How to use open source satellite data for your investigative reportingHow to use open source satellite data for your investigative reporting](https://towardsdatascience.com/how-to-use-open-source-satellite-data-for-your-investigative-reporting-d662cb1f9f90)
  [Revisit and Coverage of Sentinel-2](https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/revisit-coverage)