梅雨明けが例年になく早く明けたと思えば、雨と夏の暑さが交互に訪れて寝苦しい今日この頃、皆様いかがお過ごしでしょうか。今回の記事では、ちょっと気が早いですが、今秋に向けて衛星データを使った紅葉予測ができないかを検証してみました。

本記事で紹介する方法はまだまだ未完成のものにはなりますが、前例の少ない衛星データ活用案を実用化するための挑戦過程として、参考になるかと思います。また、今回は解析が若干長くなってしまうため、前編(本記事)と後編の２つの 記事として公開予定です。

1. 企画のきっかけ

最近、衛星データの産業活用はかなり盛んになってきたように思います。ただ、衛星データとSNSのデータを掛け合わせるような活用例はまだ多くは見られないのではないでしょうか。今回はそのような「衛星データ×SNS」という切り口から衛星データの活用例を探ってみました。

そこで思い当たったのが、YAMAPのリアルタイム紅葉モニターです。株式会社ヤマップは、登山地図GPSアプリ「YAMAP」を提供し、遊んだ山の情報や写真を共有するSNS機能もあります。なかでも、上述のリアルタイム紅葉モニターは、登山者がアプリ上のSNSに投稿した登山道中の写真を解析し、山の紅葉具合をリアルタイムでモニターするサービスです。

YAMAPは、リアルタイム紅葉モニターの他、「リアルタイム桜モニター」を公開しており、2022年6月には環境省と協力した取り組み「ライチョウモニター』」を公開しています。

これらの登山者さんのSNSデータ(登山道の写真)と衛星データを組み合わせることで、宇宙から紅葉時期の予測が行えないかと考えました。

2. 取っ掛かりをつかむ

実際にとりかかろうとした際、当初想定していた以上に難易度が高いことに気づきました。

例えば、衛星データに映っているゴルフ場や飛行機を検出するといったタスクであれば、機械学習ですでに標準的な手法(e.g.YOLO)が数多くあるので、それらを適用すればある程度の精度で実行可能です。

一方で、衛星データを使って紅葉予測をしようと思っても、標準的な手法は確立しておらず、何から始めればいいのか検討もつきません……。そこで、まずは紅葉予測に関する先行研究を探してみました。具体的には、”remote sensing”, “Fall foliage”, “satellite”などのキーワードで論文を検索してみました。本記事の最後にいくつか代表的なものを参考文献として載せておいたので、興味のある方はご覧になって下さい。

先行研究を調べたところ、多くの研究ではNDVI (Normalized Difference Vegetation Index)などの植生指数を用いて紅葉解析を行っているようです。しかしながら、論文ではアクセスの難しい衛星データを使用していたり、複雑なデータ処理をしていることが多く、完全な再現は困難でした。そこで、まず、今回の記事ではSentinel-2から得られたNDVIを用いて簡易的な紅葉予測ができるかどうかを検証してみることにしました。

3. 地上データの準備

衛星データを用いて紅葉予測を行っても、衛星データだけからはその結果が正しいかどうかが分かりません。そのため、予測モデルの調整や検証のために地上データが必要となります[1]。今回はそのような地上データとして、YAMAPのリアルタイム紅葉モニターの位置情報も分かる画像データを利用しました。

4. 衛星データ(Sentinel-2)の準備

さて、次に紅葉予測に用いる衛星データです。衛星データに関しては、最近AWS Open DataにてSnetinel-2 Cloud-Optimized GeoTIFF (COG)データセットが公開されたので、こちらを使用します。こちらのデータセットはSentinel-2 Public DatasetのシーンすべてがCOGに変換されて提供されています。詳細はこちらの記事をご覧ください。

まずは緯度、経度、日時範囲を指定したときに、該当するSentinel-2Aのシーンを検索するための関数を定義します。

def get_catalog(lon, lat, distance_km, dates):
    """
    sentinel-2衛星のデータを検索し、検索結果を返す関数。
    ・入力
      - lon: 経度
      - lat: 緯度
      - distance_km: どれほどの空間的な範囲(km)を見るか 
      - date: 検索する日時の範囲

    ・出力
      - catalog: シーンの検索結果
    """
    R_earth = 6377.397 # km # 赤道半径
    eps = (distance_km/R_earth)*(2.0*90/math.pi) # km単位で指定した距離を緯度経度の単位に変換する。
    bbox = [lon-eps, lat-eps, lon+eps, lat+eps]
   #print(eps)
    URL='https://earth-search.aws.element84.com/v0'
    try:
        results = satsearch.Search.search(url=URL,
                                          collections=['sentinel-s2-l2a-cogs'], # note collection='sentinel-s2-l2a-cogs' doesn't work
                                          datetime=dates,    
                                          bbox=bbox,    
                                          sort=[''])    
        print('%s items have been found!' % results.found())
        items = results.items()
        catalog = intake_stac.catalog.StacItemCollection(items)
    except KeyError:
        print("Error")

    return catalog

5. 衛星データ(Sentinel-2)を用いたNDVIの時間進化の計算

次に、指定した領域(緯度、経度、領域の半径)におけるNDVIの時間進化を計算するための諸関数を定義します。

def get_mean_NDVI(item, lon, lat, distance_km, dates):
    """
    対称領域でのNDVIの平均値を計算する関数。
    ・入力
      - item: 特定のシーン画像
    　- lon: 経度
      - lat: 緯度
      - distance_km: どれほどの空間的な範囲(km)を見るか 
      - date: 検索する日時の範囲
    ・出力
      - date: 時刻
      - ndvi: 平均のNDVI 
    """
    date = item.metadata['datetime']
    if date in dates:
        return -1, -1

    try:
        da = item.B04().to_dask()
    except:
        print("Error in opening the file.")
        return -1, -1
    bands = ['nir','red']
    stack = item.stack_bands(bands)
    da = stack().to_dask()
  
    # Reorganize into xarray DataSet with common band names
    da['band'] = bands
    ds = da.to_dataset(dim='band')
    
    # Extract and plot NDVI subset
    NDVI = (ds['nir'] - ds['red']) / (ds['nir'] + ds['red']+1e-8)
    (x_min, y_min, x_max, y_max) = get_coordinates_for_region_of_interest(lon, lat, distance_km, item)
    
    #subset = NDVI.sel(y=slice(9.63e6,9.64e6), x=slice(7.55e5,7.71e5))
    subset = NDVI.isel(y=slice(y_min,y_max), x=slice(x_min, x_max))
    # subset.plot.imshow(cmap='BrBG', vmin=-1, vmax=1)

    ndvi_arr = np.array(subset)
    ndvi_arr = ndvi_arr[(ndvi_arr<1.0) & (ndvi_arr>=-1)] # 異常値を除外。
    return date, np.nanmean(ndvi_arr)

def get_NDVI_time_evolution(lon, lat, distance_km, dates, cloud_cover_thres):
    """
    指定した緯度、経度、距離範囲に対してNDVIの時間進化を計算する関数。
    ・入力
      - lon: 経度
      - lat: 緯度
      - distance_km: どれほどの空間的な範囲(km)を見るか 
      - date: 検索する日時の範囲
      - cloud_cover_thres: 雲の被覆率
    ・出力
      - dates: 時刻のリスト
      - ndvi: NDVI平均値のリスト
    """
    catalog = get_catalog(lon, lat, distance_km, dates)
    catalog_list = list(catalog)
    catalog_list_little_cloud = [catalog_list[i] for i in range(len(catalog_list)) if catalog[catalog_list[i]].metadata['eo:cloud_cover']<=cloud_cover_thres]
    
    dates = []
    ndvis = []
    # 検索結果すべてに対して、時刻と指定した範囲におけるNDVIの平均値を計算。
    for i in range(len(catalog_list_little_cloud)):
        filename = catalog_list_little_cloud[i]
        item = catalog[filename]
        date_, ndvi_ = get_mean_NDVI(item, lon, lat, distance_km, dates)
        if ndvi_==-1 or np.isnan(ndvi_):
            continue
        dates.append(date_)
        ndvis.append(ndvi_)
        
        print(date_, ndvi_)

    print("Done")
    return dates, ndvis
    

def get_coordinates_for_region_of_interest(lon, lat, distance_km, item):
    """
    (lon, lat)の周囲distance_kmの四角領域の座標を取得する。
    """
    # シーンの横幅と高さを取得。
    width = item.B04().to_dask().sizes['x']
    height = item.B04().to_dask().sizes['y']

    # 指定した領域のバウンディングボックスを取得。
    R_earth = 6377.397 # km # 赤道半径
    eps = (distance_km/R_earth)*(2.0*90/math.pi) 
    # distance_km/R_earth # km単位で指定した距離を緯度経度の単位に変換する。
        
    bb_scene = item.metadata['bbox']
    lon_min = min(bb_scene[0], bb_scene[2])
    lon_max = max(bb_scene[0], bb_scene[2])
    lat_min = min(bb_scene[1], bb_scene[3])
    lat_max = max(bb_scene[1], bb_scene[3])
    
    # print(lon_min, lon_max, lat_min, lat_max)
    # 抽出したい領域のx, y座標を取得。
    eps_width =  eps*width/(lon_max-lon_min)
    eps_height = eps*height/(lat_max-lat_min)
  
    x_min = (lon - lon_min)*width/(lon_max-lon_min)
    y_min = (lat - lat_min)*height/(lat_max-lat_min)
    x_max = x_min + eps_width
    y_max = y_min + eps_height

    return (int(x_min), int(y_min), int(x_max), int(y_max))

上で定義した関数を用いて、高尾山周辺のNDVIの時間進化を計算してみましょう。

# 高尾山 (引用: www.kyorikeisan.com)
lat_takao = 35.6254
lon_takao = 139.2437 

takao_date_list_2020, takao_ndvi_list_2020 = get_NDVI_time_evolution(lon_takao, lat_takao, distance_km=3, dates='2020-01-01/2020-12-31', cloud_cover_thres=10)

出力は以下のようになります。

72 items have been found!
2020-12-23 01:37:18+00:00 0.6974050338522484
2020-12-18 01:37:20+00:00 0.6318476175932995
2020-12-08 01:37:19+00:00 0.6531392720332513
2020-11-28 01:37:22+00:00 0.7416216981449081
2020-11-23 01:37:22+00:00 0.7148621431468472
2020-11-18 01:37:25+00:00 0.6410134524599524
2020-11-13 01:37:23+00:00 0.7236792606806962
2020-10-24 01:37:25+00:00 0.13458207612678422
2020-09-09 01:37:27+00:00 0.8267262952678611
2020-09-04 01:37:25+00:00 0.8048793089895871
2020-08-20 01:37:28+00:00 0.7866235166038421
2020-08-15 01:37:25+00:00 0.7934460316164677
2020-05-17 01:37:21+00:00 0.8013247270434646
2020-05-07 01:37:19+00:00 0.792317609013527
2020-05-02 01:37:25+00:00 0.7416633791494799
2020-04-02 01:37:18+00:00 0.6495288378869245
2020-03-18 01:37:20+00:00 0.6068543433629968
2020-02-02 01:37:14+00:00 0.6740520167271628
2020-01-13 01:37:16+00:00 0.6468788635874622
2020-01-03 01:37:16+00:00 0.6580539780895539
Done

同様に粟駒山、九重連山についても計算して、図にすると以下のようになりました！

秋から冬にかけてNDVIが減少していく傾向が見られますね。データを見ると、高尾山であれば10~12月(10月のデータは不明)、九重連山であれば10~11月が色づき始めるタイミングなのではないかと予測できます。

実際に、YAMAPのデータでも、以下のように実際に同時期に紅葉している写真が投稿されていると確認できました。また、それぞれの場所で1カ月後には見頃を追えて、落葉してしまっていることも分かります。

高尾山

九重連山

6.解析を進めるにあたっての問題点とNDVIの時間進化を用いた紅葉予測

NDVIと実際のデータを見比べた上で、解析するにあたって、少なくとも以下2つの問題点があると考えられます。

・いくつかのデータ点においてはNDVI値が０にとても近くなっており、おそらくこれは雲の影響と推測される。
・場所によってデータ点数が異なり、例えば栗駒山ではかなりデータ点数が少ない
※今回利用している衛星データはSentinel-2のみであり、有料の衛星データを購入することで観測頻度を増やすことができます

合わせて、以下のようなデータ以外の問題点も見えてきました。

・衛星データを用いた紅葉予測は前例が少なく、解説記事などがほとんど見当たらない。そのため英語で書かれた論文にあたるほかなく、当分野の非専門家である自分にはかなり困難に感じる部分が多い。

・それらの論文ではアクセスが難しい衛星データを使用している場合もあり、結果の再現が難しい。

上記の問題点があると把握したうえで、NDVIの時間進化を用いて紅葉予測を行います。文献[2]を読むと、まずこれらのNDVIの時間発展をシグモイド関数でフィッティングをします。次に得られたパラメータを用いてTNBI(Temprally-Normalized Brownness Index)という量を計算し、これらの値の大きさに基づいて紅葉具合の予測を行っています。

今回の記事でこれらの解析すべてを行うととても長い記事になってしまうので、解析の続きはもう一つの記事として、より詳細なモデルの作成を行い、公開する予定です。

7.参考文献

・Linking Spaceborne and Ground Observations of Autumn Foliage Senescence in Southern Québec, Canada
・Prototype for monitoring and forecasting fall foliage coloration in real time from satellite data

YAMAPの取り組みを見る

ライチョウモニター

リアルタイム紅葉モニター

おすすめの関連記事