2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

衛星データから空港周辺の緑地率を出し、バードストライク(飛行機が鳥と衝突すること)率との相関があるのか、調査してみました。

結果として、衛星データから緑地率を調査することはできたものの、バードストライク率との相関を出すことについてはまだまだ考える余地があるようです。

バードストライクが航空会社に与える影響含め、空港周辺の緑地率を算出する方法、また、調査の結果をご紹介します。

※本記事は宙畑メンバーが気になったヒト・モノ・コトを衛星画像から探す不定期連載「宇宙データ使ってみた-Space Data Utilization-」の第15弾です。まだまだ修行中の身のため至らない点があるかと思いますがご容赦・アドバイスいただけますと幸いです！

(1)バードストライクとは？ 概要と影響

バードストライクとは、航空機が鳥と衝突し、窓ガラスのひびや胴体の損傷、また、エンジン部に鳥が入り込むことで、運航への支障をきたしてしまう可能性の事故です。

2018年のバードストライク発生件数は1434件で、そのうち引き返しや離陸の停止は24件。日本国内の離着陸回数が2,632,166件であることを考えると、微小な影響のようにも思えますが、バードストライクを要因とした航空機墜落による死亡事故も世界では発生しています。

現状の対策としては、空港内の目視巡視や専用レーダーを導入しているケースもあるようですが、決して事故が0になっているということはありません。

(2)周りに植物が多い空港ほどバードストライク率高いという仮説

では、バードストライクがいつ、どのような場所で起こるかに違いはあるのでしょうか。

上のスライドはバードストライクが離陸から飛行中、到着までのどこで発生しているかを表すデータです。

どうやら飛行中ではなく、離発着時の事故が最も多いようです。空港周辺の鳥の生息状況、つまり、「鳥が生息する木々がある」がバードストライク率に関係していると考えられます。

そこで、今回はひとつの仮説「空港周辺に植物が多いところほどバードストライク率が高い」を立て、調査を進めました。

(3)検出方法と使用するデータの説明

では、空港周辺の植物の状況をどのように調査するのか。そこで利用するのが衛星データです。

衛星データでは地上約500km以上の宇宙空間より、地上の状況を定期的に観測することができます。そして、衛星データからは人間の目で見える可視光以外にも様々な光波長をセンサで捉えることができ、その一つの特性として、地上の植物の状況を把握することが可能です。

詳しくは「光の波長って何？ なぜ人工衛星は人間の目に見えないものが見えるのか」をご覧ください。

今回は、Tellusに搭載されているAVNIR-2という光学センサのデータを活用して、国内28の拠点空港(※)周辺の緑地率を調べてました。
※空港法第４条第１項各号に掲げる空港（成田国際空港、東京国際空港、中部国際空港、 関西国際空港、大阪国際空港並びに国際航空輸送網又は国内航空輸送網の拠点となる空港）

(4)Tellusの衛星データを利用した空港周辺の緑地率の算出方法

それでは実際に、空港周辺の緑地率を確認していきましょう。

まずは必要なライブラリをインストールしておきます。

import os, json, requests, math
from pathlib import Path
from skimage import io
from io import BytesIO
from skimage.draw import polygon
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

次に空港周辺の植物活性度のデータを取得する記述をします。

植物活性度はpresetのndviを利用し、空港周辺のデータを取得するのは「【ゼロからのTellusの使い方】大きなタイル画像を取得するには？地図タイルを扱う上で知っておきたいTips4選」で紹介しているように、(空港の入った)特定のタイルを中心にした周り一つずつのタイルを連結させて大きなタイル画像を取得します。

TOKEN = '自分のトークンを入れる'
DOMAIN = 'tellusxdp.com'

def get_AVNIR2_brend_image(option):
    url = 'https://gisapi.{}/av2ori/{}/{}/{}/{}/{}.png'.format(DOMAIN, option['rspid'], option['product_id'], option['zoom'], option['xtile'], option['ytile']) 
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }

    r = requests.get( url, params=option['payload'], headers=headers )
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return io.imread(BytesIO(r.content))  


def get_AVNIR2_series_image(rspid, product_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1, payload={}):
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            option = {
                'rspid': rspid,
                'product_id': product_id,
                'zoom': zoom,
                'xtile': topleft_x + x,
                'ytile': topleft_y + y,
                'payload': payload
            }
            row.append(get_AVNIR2_brend_image(option))
        img.append(np.hstack(row))
    return  np.vstack(img)

try:
    img = get_AVNIR2_series_image(
        'D083P3', 
        'ALAV2A035302940', 14, 14146, 6607, 3, 3, {'preset': 'ndvi'})
except ValueError as e:
    print(e)
# 熊本空港
io.imshow(img)

無事、熊本空港周辺の植物活性度を確認する大きなタイル画像を取得することができました。今回はzoomレベルを14(1枚のタイルにちょうど空港がまるっと収まるくらいのzoomレベル)に統一して取得していきます。

いよいよ植物の割合を算出する記述を書きます。まず、RGBからHSVに色順を入れ替えます。この記述を加えることで、色をひとつなぎのものとして扱うことができます。

※HSV変換については「【ゼロからのTellusの使い方】Jupyter NotebookでAVNIR-2/ALOSの光学画像から緑地割合算出」で解説をしております

次は、自分のお好みの色味(今回では緑色)以外の色を黒色にします。HSVに変換したことで、衛星画像の色を数値化させてお好みの操作を加えることができるのです。
※緑色と指定した範囲について、本来厳密にやるべき作業ですが、今回は厳密な検証を行ったものではありません。

そして、黒でマスキングされたデータから、緑色の割合を算出する記述をします。今回は取得するタイル画像が人の目にも分かりやすいよう、再度RGBに戻す記述も加えています。

import cv2
def masking(img, hsv_min, hsv_max):
    # OpenCVで扱うために色順をRGB(RGBA)からBGRに入れ替える
    cv_bgr = cv2.cvtColor(img, cv2.COLOR_RGBA2BGR) 
    # BGRからHSVに変換
    cv_hsv = cv2.cvtColor(cv_bgr, cv2.COLOR_BGR2HSV)
    
    # 指定した範囲の色以外を黒に変換する
    mask = cv2.inRange(cv_hsv, hsv_min, hsv_max)
    masked_img = cv2.bitwise_and(cv_bgr, cv_bgr, mask=mask)

    # 指定した範囲の色の割合
    rate = len(np.where(mask.ravel()!=0)[0])/len(mask.ravel())

    # 得られた画像の色順をBGRからRGBに戻す
    return cv2.cvtColor(masked_img, cv2.COLOR_BGR2RGB), rate


green_min = np.array([25,0,0]) 
green_max = np.array([75,255,255]) 

masked_img1, rate = masking(img, green_min, green_max)
io.imshow(masked_img1)
print(rate)

そして出てきたのがこちら。

きちんと緑色以外はマスキングされてますね。画像上の数字が熊本空港周辺の緑地率。熊本空港周辺の緑地率は「0.7359178331163194」、つまり、約73.6%だと分かりました。

この作業をあと27回繰り返して、拠点空港の緑地データを一覧にまとめたのがこちらです。
※今回はAVNIR2のデータの中でも、(緑が茂っているであろう)5月から10月の期間でタイル画像に雲がかかっていないデータを1シーンずつ選んでいます

(5)相関性を確認する

では、空港周辺の緑地率が揃ったところで、バードストライク率との相関性を見てみましょう。2018年のバードストライク率(※)と合わせて数字を並べてみたのがこちらの表。
※バードストライク率は、離着陸回数10,000あたりのバードストライク数になります。

データの散布図は以下のようになりました。すでに負の相関があるようで、どうやら最初に立てた仮説が間違っているようです。

相関係数、P値を求めたところ、以下のようになりました。

相関係数：-0.2476165351
P値：0.2039348279

仮説は大きく外れ、負の相関がある上に、P値も統計的に有意とは言えない値となってしまいました。この結果が出てあらためて調べてみたところ、バードストライクの要因となるのは海に生息する鳥が多いようです。※なぜ最初に気づけなかったのか甚だ疑問です苦笑

そこで、取得したタイルに海が入っている画像と入っていない画像を分けて相関を出したのが以下の値です。

【海が入っている】

相関係数：-0.1601917078
P値：0.5684743732

【海が入っていない(内陸)】

相関係数：0.4998903329
P値：0.0819419944

取得した衛星画像に海が入っている場合の相関については、P値が0.568と、相関があるとはいえないデータに……。海が入っていない、内陸の空港周辺のデータのみについては、正の相関が見られました。ただし、これもP値が0.08と、0.05に近い数字ではあるものの、まだまだ統計的に有意であるとは言いづらい結果となってしまいました。

(6)まとめ

今回、衛星データを利用して空港周辺の緑地率を算出するまではできたものの、バードストライク率との相関性までは仮説通りの結果を出すまでには至りませんでした。

内陸の空港については、P値に課題があったものの、今回調査の対象としていない54の地方管理空港も調査するほか、例えばLandsatといった他の衛星データも確認する、海からの距離を算出する、zoomレベルを変更してより広範囲を見るようにする、植物の種類(林なのか、畑なのか)を見分ける、緑色の閾値を変更する……など、複数の調整を加えることで、相関があると言える可能性も見えたように思います。

海に近い空港についても、海に対しての離着陸の向きや、温度、風の向きや強さといった周辺の状況を踏まえたデータと合わせることで、バードストライクが起こりやすい空港を把握できるかもしれません。

そして、今回取り組んだ空港周辺の状況とバードストライク率の関係性調査の精度が高まれば、今後新しい空港を建設する際の参考になるかもしれません。

衛星データを活用した新しい施設の建設候補地探しは今後活用が期待される事例のひとつ。ぜひ衛星データを使った〇〇との相関を皆様も調査してみませんか？