<jupyter_start><jupyter_text>Sentinel-2の観測画像をpythonのjupyter noteboook(lab)のプラットフォームに取得し，火災検知に特化した衛星画像の５日間毎の時系列変化（Timelapse）のGifアニメーションを作成 目的：pythonのjupyter labにて，希望する地域を含むのSentinel-2の観測画像ファイルを自動でダウンロードする．その画像より赤外波長領域の画像作成し，その画像より関心域の時系列のGIFアニメーションを作成する．  関連するモジュールのインストールおよびインポート SentinelのAPIを含め，必要なモジュールを総合開発環境(python）にインストールする．<jupyter_code>!pip install sentinelsat
!pip install rasterio
!pip install folium
!apt install gdal-bin python-gdal python3-gdal 
!apt install python3-rtree 
!pip install git+git://github.com/geopandas/geopandas.git
!pip install descartes<jupyter_output><empty_output><jupyter_text>geopandasを使うための関連モジュールもインストールする．<jupyter_code>!pip install shapely
!pip install fiona
!pip install six
!pip install pyproj
!pip install descartes
!pip install geopandas<jupyter_output><empty_output><jupyter_text>処理・解析に必要なパッケージをインポートする．<jupyter_code>import folium
import os
import datetime

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt 
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from shapely.geometry import MultiPolygon, Polygon
import rasterio as rio
from rasterio.plot import show
import rasterio.mask
import fiona

import zipfile
import shutil
import glob

from osgeo import gdal

from PIL import Image, ImageDraw, ImageFont<jupyter_output><empty_output><jupyter_text>関心域の緯度経度情報の取得と地図ポリゴンの作成 海外の関心域の緯度・経度情報を取得するために，米国キーン州立大学の以下のWEBアプリケーションを利用します． 関心域（四隅）を右クリックで指定し，ポリゴンを作成します．枠を作成したあと，”Close Shape”をクリックし形を整えます． 関心域を新たに取得する場合は，”Reset”をクリックしてポリゴンを消去する．注意：基本的に画像を拡大後，左回りで関心域（例えば，東京など）を検索し，緯度・経度を取得しています．右回りで探したときは，経度情報が逆転しますのでご注意ください．<jupyter_code>#関心領域のポリゴン情報の取得．
from IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://www.keene.edu/campus/maps/tool/" frameborder="0"></iframe>')<jupyter_output>/opt/conda/lib/python3.7/site-packages/IPython/core/display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")<jupyter_text>作成したポリゴンの緯度・経度情報をコピーし，以下に貼り付けます．<jupyter_code>#上記にて取得した地理情報をコピペ．
AREA =  [
         [
        -210.4512048,
        -35.2522692
      ],
      [
        -210.4512048,
        -35.2477833
      ],
      [
        -210.4525244,
        -35.8537092
      ],
      [
        -209.5915246,
        -35.8444562
      ],
      [
        -209.5799482,
        -35.2564395
      ],
      [
        -210.4512048,
        -35.2522692
      ]
    ]<jupyter_output><empty_output><jupyter_text>経度情報が３６０°異なるため（左回り検索の場合），以下の処理で標準に修正します．<jupyter_code>for i in range(len(AREA)):
    AREA[i][0] = AREA[i][0] +360<jupyter_output><empty_output><jupyter_text>関心域の地理情報ポリゴンの作成<jupyter_code>from geojson import Polygon

# no hole within polygon
m=Polygon([AREA])<jupyter_output><empty_output><jupyter_text>今後使用する任意のファイル名をセットする．<jupyter_code>object_name = 'Australia_wildfire_new3'
import json
with open(str(object_name) +'.geojson', 'w') as f:
    json.dump(m, f)
footprint_geojson = geojson_to_wkt(read_geojson(str(object_name) +'.geojson'))<jupyter_output><empty_output><jupyter_text>foliumの地図閲覧機能を用いて関心域を確認する．<jupyter_code>m = folium.Map([(AREA[0][1]+AREA[len(AREA)-1][1])/2,(AREA[0][0]+AREA[len(AREA)-1][0])/2], zoom_start=10)

folium.GeoJson(str(object_name) +'.geojson').add_to(m)
m<jupyter_output><empty_output><jupyter_text>フォントファイルのダウンロード及び設定準備 作成した画像ファイルに，衛星画像の取得日やクレジットを記載します．このとき，フォントファイルがない場合はダウンロードし準備する必要があります。 いくつか無料でフォントを提供されていますが，今回は以下のサイトよりダウンロードします．[M+ FONTS](http://mplus-fonts.osdn.jp/about-en.htmldownload-1)<jupyter_code>!wget https://osdn.net/dl/mplus-fonts/mplus-TESTFLIGHT-063a.tar.xz
!xz -dc mplus-TESTFLIGHT-*.tar.xz | tar xf -
fontfile = "./mplus-TESTFLIGHT-063a/mplus-1c-bold.ttf" #downloadしたfontファイルの保存先を指定します．<jupyter_output><empty_output><jupyter_text>Sentinel-2の衛星画像取得，およびGIFアニメーション作成の関数の準備． コペルニクスハブのアカウントは以下のサイトなどを参考に取得してください．[無料で最新の衛星画像を入手する方法．](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)<jupyter_code>user = '********' #ご自身のCopernicus hubのアカウント情報
password = '********' #ご自身のCopernicus hubのアカウント情報

api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')<jupyter_output><empty_output><jupyter_text>一般的な人工衛星の観測画像（True Color）とは異なり，火災検知を見やすくするためにSentinel-2の赤外波長の画像（B12, B11, B8A）を使用します．<jupyter_code>def Sentinel2_get_init(i):
    products = api.query(footprint_geojson,
                     date = (Begin_date, End_date), #画像取得月の画像取得
                     platformname = 'Sentinel-2',
                     processinglevel = 'Level-2A', #Leve-1C
                     cloudcoverpercentage = (0,100)) #被雲率（０％〜５０％）
    
    products_gdf = api.to_geodataframe(products)
    
    #画像取得月の被雲率の小さい画像からソートする．
    products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
    title = products_gdf_sorted.iloc[i]["title"]
    print(str(title))
    
    uuid = products_gdf_sorted.iloc[i]["uuid"]
    product_title = products_gdf_sorted.iloc[i]["title"]
    
    #画像の観測日を確認．
    date = products_gdf_sorted.iloc[i]["ingestiondate"].strftime('%Y-%m-%d')
    print(date)
    
    #Sentinel-2 data download
    api.download(uuid)
    file_name = str(product_title) +'.zip'
    
    
    #ダウロードファイルの解凍．
    with zipfile.ZipFile(file_name) as zf:
        zf.extractall()
    
    #フォルダの画像ファイルのアドレスを取得．
    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)
    
    #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()
    
    #ポリゴン情報の取得．
    nReserve_geo = gpd.read_file(str(object_name) +'.geojson')
    
    #取得画像のEPSGを取得
    epsg = b12.crs
    
    nReserve_proj = nReserve_geo.to_crs({'init': str(epsg)})
    
    #マスクTiffファイルの一時置き場
    os.makedirs('./Image_tiff', exist_ok=True)

    #カラー合成画像より関心域を抽出．
    with rio.open(str(object_name) +'.tiff') as src:
        out_image, out_transform = rio.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('./Image_tiff/' +'Masked_' +str(object_name) +'.tif', "w", **out_meta) as dest:
        dest.write(out_image)
    
    #抽出画像のjpeg処理．
    scale = '-scale 0 255 0 15'
    options_list = [
        '-ot Byte',
        '-of JPEG',
        scale
    ] 
    
    options_string = " ".join(options_list)
    
    #ディレクトリの作成 
    os.makedirs('./Image_jpeg_'+str(object_name), exist_ok=True)
    
    #jpeg画像の保存．
    gdal.Translate('./Image_jpeg_'+str(object_name) +'/' + str(date) + 'Masked_' +str(object_name) +'.jpg',
                   './Image_tiff/Masked_' +str(object_name) +'.tif',
                   options=options_string)
    
    #画像への撮像日を記載
    img = Image.open('./Image_jpeg_'+str(object_name) +'/' + str(date) + 'Masked_' +str(object_name) +'.jpg')


    x = img.size[0]/100 #日付の記載位置の設定
    y = img.size[1]/100 #日付の記載位置の設定
    fs = img.size[0]/50 #日付のフォントサイズの設定
    fs1 = int(fs)
    obj_draw = ImageDraw.Draw(img)
    #obj_font = ImageFont.truetype("./font/mplus-1c-bold.ttf", fs1)
    obj_font = ImageFont.truetype(fontfile, fs1)
    obj_draw.text((x, y), str(date), fill=(255, 255, 255), font=obj_font)
    obj_draw.text((img.size[0]/2, img.size[1]-y - img.size[1]/20 ), 'produced from ESA remote sensing data', fill=(255, 255, 255), font=obj_font)
    

    img.save('./Image_jpeg_'+str(object_name) +'/' + str(date) + 'Masked_' +str(object_name) +'.jpg')
    
    #Copernicus hubからのダウンロードファイル，および解凍フォルダの削除
    shutil.rmtree( str(product_title) + '.SAFE')
    os.remove(str(product_title) +'.zip')
    
    return
def Sentinel2_get():
    products = api.query(footprint_geojson,
                     date = (begin, end), #取得希望期間の入力
                     platformname = 'Sentinel-2',
                     processinglevel = 'Level-2A', #Leve-1C
                     cloudcoverpercentage = (0,100)) #被雲率（０％〜５０％）
    
    products_gdf = api.to_geodataframe(products)
    products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
    
    #同一シーンの画像を取得するため，placenumberを固定する．
    products_gdf_sorted = products_gdf_sorted[products_gdf_sorted["title"].str.contains(placenumber)]
    title = products_gdf_sorted.iloc[0]["title"]
    print(str(title))
    
    uuid = products_gdf_sorted.iloc[0]["uuid"]
    product_title = products_gdf_sorted.iloc[0]["title"]
    
  
    date = products_gdf_sorted.iloc[0]["ingestiondate"].strftime('%Y%m%d')
    date1 = products_gdf_sorted.iloc[0]["ingestiondate"].strftime('%Y-%m-%d')
    print(date)
    
    api.download(uuid)
    file_name = str(product_title) +'.zip'
    
    with zipfile.ZipFile(file_name) as zf:
        zf.extractall()
    
    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)
    
    #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()
    
    os.makedirs('./Image_tiff', exist_ok=True)
    
    nReserve_geo = gpd.read_file(str(object_name) +'.geojson')
    epsg = b12.crs
    
    nReserve_proj = nReserve_geo.to_crs({'init': str(epsg)})

    with rio.open(str(object_name) +'.tiff') as src:
        out_image, out_transform = rio.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('./Image_tiff/' +'Masked_' +str(object_name) +'.tif', "w", **out_meta) as dest:
        dest.write(out_image)
    
    from osgeo import gdal

    scale = '-scale 0 255 0 15'
    options_list = [
        '-ot Byte',
        '-of JPEG',
        scale
    ] 
    options_string = " ".join(options_list)
    
    os.makedirs('./Image_jpeg_'+str(object_name), exist_ok=True)

    gdal.Translate('./Image_jpeg_'+str(object_name) +'/' + str(date) + 'Masked_' +str(object_name) +'.jpg',
                   './Image_tiff/Masked_' +str(object_name) +'.tif',
                   options=options_string)
    
    #画像への撮像日の記載
    img = Image.open('./Image_jpeg_'+str(object_name) +'/' + str(date) + 'Masked_' +str(object_name) +'.jpg')


    x = img.size[0]/100 #日付の記載位置の設定
    y = img.size[1]/100 #日付の記載位置の設定
    fs = img.size[0]/50 #日付のフォントサイズの設定
    fs1 = int(fs)
    obj_draw = ImageDraw.Draw(img)
    obj_font = ImageFont.truetype(fontfile, fs1)
    obj_draw.text((x, y), str(date1), fill=(0, 0, 128), font=obj_font) #Red:(255, 0, 0), Blue:(0,0,255), Navy:(0,0,128),Yellow:(255,255,0), bluegreen:(0,128,128),Perple:(128,0,128)
    obj_draw.text((img.size[0]/2, img.size[1]-y - img.size[1]/20 ), 'produced from ESA remote sensing data', fill=(0, 0, 128), font=obj_font)
    img.save('./Image_jpeg_'+str(object_name) +'/' + str(date) + 'Masked_' +str(object_name) +'.jpg')

    
    #Copernicus hubからのダウンロードファイル，および解凍フォルダの削除
    shutil.rmtree( str(product_title) + '.SAFE')
    os.remove(str(product_title) +'.zip')
    
    return<jupyter_output><empty_output><jupyter_text>対象画像の確認 ここでは，選択した画像が関心域を多く含むものかどうか，タイル画像番号を確認します．<jupyter_code>#L2Aデータ取得のため，２０１９年１月以降を開始日とすること．
Begin_date = '20191201'
End_date = Begin_date[:6] +'28'#１月間の観測画像より，関心域を多く含む画像のタイル番号を取得する．<jupyter_output><empty_output><jupyter_text>画像取得の開始月の画像を確認する． 被雲率の小さい画像から確認するため，デフォルトの引数は０とする．<jupyter_code>Sentinel2_get_init(4)<jupyter_output>S2B_MSIL2A_20191226T000239_N0213_R030_T55HGA_20191226T014948
2019-12-26<jupyter_text>上記の処理により作成されたjpeg画像（object_nameのjpeg画像のディレクトリに保存されている）を確認する．関心域を多く含む画像であればOK.端部のみや，関心域の画像が少ない場合は0を１に変更して，異なる衛星画像で同処理を行いjpeg画像を確認する． 上記の画像で問題ければ，placenumber(例えば，S2A_MSIL2A_20191012T000241_N0213_R030"_T56HKD_"20191012T020312)をコピーし，以下のコードに貼り付ける．これにより，同じ位置の衛星画像のみが取得される． Sentinel-2の衛星画像のplacenumberに関する情報は以下を参考とする．[Sentinel-2 Product Types](https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/product-types) 注） エラー”SentinelAPIError: HTTP status 200 OK: Invalid query string. Check the parameters and format.”は関心域が大きく，複数シーンにまたぐときに発信される．関心域を小さくすることでエラーが回避される．  関心域のGIFアニメーションの作成．<jupyter_code>placenumber = '_T55HGA_' #placenumberはT54SUE
#L2Aデータ取得のため，２０１９年１月以降を開始日とすること．
Begin_date = '20191201'
End_date = '2020210'<jupyter_output><empty_output><jupyter_text>開始日と終了日の日差を算出．　Sentinel-2の再観測までの期間は５日間なのでデフォルトは５日とする．<jupyter_code>m = 5 # 5日間毎の画像を取得する．
#開始日と終了日の日差の数を求める．
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
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()<jupyter_output><empty_output><jupyter_text>設定した日毎の衛星画像の取得，および関心域のjpeg画像を以下のコードにて作成する． jpeg画像の古い観測画像から順番にGifアニメーションを作成する．<jupyter_code>images =[]

files = sorted(glob.glob('./Image_jpeg_'+str(object_name) +'/*.jpg'))
images = list(map(lambda file: Image.open(file), files))

images[0].save('./Image_jpeg_'+str(object_name) +'/' + str(object_name) + '.gif', save_all=True, append_images=images[1:], duration=1000, loop=0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15145831986031172, 'fraction_numerical': 0.043781931061242965, 'mean_word_length': 4.0635232481990835, 'pattern_counts': {'":': 8, '<': 70, '<?xml version=': 0, '>': 69, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '298247', 'license': 'CC0-1.0', 'n_tokens_mistral': 7921, 'n_tokens_neox': 7494, 'n_words': 770, 'path': 'RS_Python_S2L2A_fire_Timelapse_public.ipynb', 'repo_name': 'sorabatake/article_10961_Sentinel2timelapse', 'size': '15463'}