【入門】Pythonによる人工衛星データ解析（Google Colab環境）

はじめに　〜衛星データとは〜

• 人工衛星データとは、人工衛星を利用した“リモートセンシング”によって取得されたデータを指します。
• これまで人工衛星データは専門ツールや大容量データ処理基盤が必要なため、利用できる組織は大学機関や一部の専門機関が限られていましたが、昨今のオープンソース・ライブラリの普及やデータ処理基盤のクラウド利用により、一般組織でも気軽に人工衛星データを扱える外部環境が整ってきました。
• 衛星データを利用することで、これまで取得することができなかった様々な場所・時間・対象の状態をビッグデータで解析することが期待できます。
• そこで本記事では、どの様にデータを扱うのかを、衛星データ解析の専門ツールを利用せず（最も身近なツールの一つであるpythonを利用）、誰でも気軽に試すために無償で利用方法を紹介していきたいと思います。
• また、今回はビジネスや社会実装に利用イメージが沸きやすい衛星データセットを選定してみました。応用できそうな場面を想像しながら、読んで頂ければと思います。

対象データ

• 今回の記事では、夜間光データを例に衛星データセットの利用方法を紹介します。 Image and Data processing by NOAA's National Geophysical Data Center. DMSP data collected by the US Air Force Weather Agency.

夜間光データ（NightTime Lights）とは

• 人工衛星が観測しているデータの一つに夜間光データというものがあります。
• 簡単に言うと、夜間の街の光の量をセンシングし続けているデータです。
• 夜間光データは、経済活動（GDPやエネルギー使用量）との相関が報告されており、経済活動の代理変数として実証研究が進んでいます。
• 特に都市単位での経済活動量の測定や正確な経済指標の測定が難しい国への応用が検討されています。

データ概要

• 夜間光データの計測はいくつかの人工衛星で測定が行われているのですが、今回は夜間光データの研究では老舗データセットであるDMSP-OLSというデータを利用します。
• スペック
  • 範囲：グローバル
  • 解像度：1kmメッシュ
  • 頻度：年次
  • 期間：1992年 ~ 2013年
  • データ量：圧縮時 ~300MB、解凍時 ~3GB
  • フォーマット：GeoTIFF
  • データダウンロード：https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html
• 備考
  • 前処理の異なる様々な夜間光データタイプが存在するが、今回は研究でもよく利用されている"Average Visible, Stable Lights, & Cloud Free Coverages"を利用

利用イメージ

• 人工衛星データはダウンロード時は全球（グローバル）サイズやリージョンサイズで巨大です。
• そのため、解析時は必要なエリアのみを切り出して利用することが一般的です。
• ダウンロードしてきたオリジナルデータ（ラスターデータ）と切り出し範囲を指定するデータ（ベクターデータ）に２つが必要になります。（どちらもWEBでダウンロードできます）

pythonでの基本操作

• 下記のライブラリを利用します（本記事では必要最低限の機能のみ紹介します）
  • データの読み込み: rasterio
  • 対象エリアの切り出し: geopandas
  • 数値計算：numpy
  • 可視化：matplotlib

解析環境

• Google Colaboratory
  • Google Colabの基本的な使い方（起動方法など）は他記事を参照ください

作業の流れ

• 1. データの取得
• 2. データの読み込み
• 3. 必要なエリアの抽出
• 4. データ可視化・解析　〜国ごとの夜間光データを可視化・比較してみる〜

1. データの取得

• DMSP-OLSのdata downloadからダウンロードできます
• アクティブのリンクをクリックするとダウンロードが開始します
• ただし、注釈で書いてある通り、ダウンロード時の圧縮ファイルは300MB程度ですが、解凍すると1ファイルあたり3GB程度になりますので、一般的なノートPCで作業するのはやや厳しいです
• そこで、データのダウンロードから解凍までの一連の処理をGoogle Colab上で行います

• （高スペックPCをお持ちの恵まれた方はポチポチとダウンロード＆解凍...で対応しても問題ないです）

• はじめに必要なライブラリをGoogle Colab上でインストールしておきます

# 必要なライブラリのインストール
!pip install sh
! pip install rasterio
!pip install geopandas

• 次にダウンロードしてきたデータセットを保存するフォルダ（ディレクトリ）を作成します

import os
from sh import wget, gunzip, mv
import tarfile

# データ保存用のディレクトリ'data'を作成
if not os.path.exists('data'):
    os.mkdir('data')

• wgetコマンドを利用して、圧縮データをダウンロードしてきます
• ダウンロードしてきたファイルを解凍します
• ファイルは.gzと.zipで二重圧縮されおりややこしいですが、最終的にtifファイルが解凍されます

# ダウンロード対象のURLを作成
target_data = 'F101992'
url = f'https://ngdc.noaa.gov/eog/data/web_data/v4composites/{target_data}.v4.tar'

# データのダウンロード
wget(url)

# 圧縮ファイルの解凍（該当ファイルのみ抽出）
with tarfile.open(f'/content/{target_data}.v4.tar') as tar:
    # stable_lights.avg_vis.tif.gzというファイル名を取得
    file = [tarinfo for tarinfo in tar.getmembers() if tarinfo.name.endswith("web.stable_lights.avg_vis.tif.gz")]
    # 対象ファイルを解凍（.gz）
    tar.extractall(path='/content/', members=[file[0]])
    # 対象ファイルを解凍（Unzip）
    gunzip(f'/content/{target_data}.v4b_web.stable_lights.avg_vis.tif.gz')
    # 対象ファイルの移動
    mv(f'/content/{target_data}.v4b_web.stable_lights.avg_vis.tif', '/content/data/')

2. データの読み込み

• 解凍したtifファイルを読み込みます
• rasterio.open('ファイルへのパス')で簡単に読み取ることができます
• 読み込んたオブジェクトread()するとデータをnumpy形式で取得することができます
• numpy形式まで変換すれば、あとはお好みで解析・可視化することができます

import rasterio
import numpy as np
import matplotlib.pyplot as plt

with rasterio.open('/content/data/F101992.v4b_web.stable_lights.avg_vis.tif') as src:
    data = src.read()# numpy形式で読み込み

# データのサイズの確認
data.shape

# データの可視化
plt.imshow(data[0])
plt.colorbar()

1992年の全球の夜間光データ

• ただし、ここではダウンロードしたままの状態（全球データ）なので、解析を行いたい単位でデータを切り出す必要があります

3. 必要なエリアの抽出

• 必要なエリアを抽出するには、切り出し範囲を指定するデータが必要になります
• 切り出し範囲を指定するデータには様々な形式が存在していますが、ここではgeojson形式のデータを利用します（他にもシェープファイルが有名です）
• こちらのサイトから国ごとの境界線のデータがダウンロードできます
• ここでもクリック＆ダウンロードでも良いのですが、wgetを利用してGoogle Colab上でダウンロードします

# ベクターファイルのダウンロード
wget('https://datahub.io/core/geo-countries/r/countries.geojson')

• 次に取得したgeojsonファイルを読み込みます
• geopandasを利用すると、geojsonをファイルを読み込めることはもちろん、pandas.dataframe様に扱うことができます
• これにより任意のエリア（国）の境界線のデータを取得します

import geopandas as gpd

# geojsonファイルの読み込み
countries = gpd.read_file('/content/countries.geojson')

# 中身の確認
countries.head()

# 国境データの可視化
countries.plot()

# 日本の国境線の抽出（通常のpandas.dataframeの操作）
countries.query('ADMIN == "Japan"')

• 次に取得した日本の国境線のデータを全球データセットに適用して、日本エリアのデータのみ抽出してみます
• 必要箇所を切り出すにはrasterio.maskというメソッドを利用します
  • rasterio.mask.mask("tifファイル"を読み込んだオブジェクト, "国境線データ"のgeometry列, crop=True)
• これによりout_image, out_transformという2つのデータが出力され、out_imageにnumpy形式で切り出したデータが格納されます（out_transformは切り出したデータの座標変換情報が格納されているのですが、ちょっとややこしいのでここでは割愛します）

import rasterio.mask

with rasterio.open('/content/data/F101992.v4b_web.stable_lights.avg_vis.tif') as src:
    out_image, out_transform = rasterio.mask.mask(src, countries.query('ADMIN == "Japan"').geometry, crop=True)

• 日本の国境線で切り出したout_imageを可視化してみます

1992年の日本の夜間光データ（首都圏の光量が大きいことが見てとれる）

4. データ可視化・解析

• ここまでで任意の国の夜間光データを取得する方法を紹介しました
• これまで紹介した方法を利用していくつかの国の夜間光データを可視化・解析してみようと思います

# 一連処理を関数化
def load_ntl(target_data, area):
    # データが存在しないときのみダウンロード
    if not os.path.exists(f'/content/data/{target_data}.v4b_web.stable_lights.avg_vis.tif'):
        url = f'https://ngdc.noaa.gov/eog/data/web_data/v4composites/{target_data}.v4.tar'
        # データのダウンロード
        wget(url)
        # 圧縮ファイルの解凍（該当ファイルのみ抽出）
        with tarfile.open(f'/content/{target_data}.v4.tar') as tar:
            # stable_lights.avg_vis.tif.gzというファイル名を取得
            file = [tarinfo for tarinfo in tar.getmembers() if tarinfo.name.endswith("web.stable_lights.avg_vis.tif.gz")]
            # 対象ファイルを解凍（.gz）
            tar.extractall(path='/content/', members=[file[0]])
            # 対象ファイルを解凍（Unzip）
            gunzip(f'/content/{target_data}.v4b_web.stable_lights.avg_vis.tif.gz')
            # 対象ファイルの移動
            mv(f'/content/{target_data}.v4b_web.stable_lights.avg_vis.tif', '/content/data/')
    # TIFファイルから該当エリアのデータを抽出
    with rasterio.open(f'/content/data/{target_data}.v4b_web.stable_lights.avg_vis.tif') as src:
        out_image, out_transform = rasterio.mask.mask(src, countries.query(f'ADMIN == "{area}"').geometry, crop=True)
    return out_image

# 可視化用の関数
def show(data):
    plt.figure(figsize=(15, 5))
    plt.subplot(121)
    plt.imshow(data[0])
    plt.subplot(122)
    plt.hist(data.reshape(-1), bins=np.arange(1, 63, 1))

# 利用例（1992年の日本のデータを取得）
japan_1992 = load_ntl(target_data='F101992', area='Japan')

-　色々な国の夜間光を可視化してみます。（target_data, areaの変数を変更すれば好きな国・年のデータが確認できます）

# 日本・中国・タイ・カンボジアのデータを取得
japan_1992 = load_ntl(target_data='F101992', area='Japan')
china_1992 = load_ntl(target_data='F101992', area='China')
thailand_1992 = load_ntl(target_data='F101992', area='Thailand')
cambodia_1992 = load_ntl(target_data='F101992', area='Cambodia')

# 可視化
show(japan_1992)
show(china_1992)
show(thailand_1992)
show(cambodia_1992 )

1992年の日本の夜間光データ

1992年の中国の夜間光データ

1992年のタイの夜間光データ

1992年のカンボジアの夜間光データ

• こうしてみると、日本は1992年時点でも全体的に明るい、中国は国土が広いため都市間の差が大きい、タイはバンコクを中心に幹線道路や地方都市が点々と明るい、カンボジアは首都以外は暗い、、、といった国ごとの傾向が見て取れます
• ここでは上記のパターンの解析にとどめますが、他の年代や国々を比較しても面白いかと思います
• 研究領域はSum of NTL(Night time Light）（エリアごとの合計夜間光量）を指標して、GDPやエネルギー消費量インデックスと比較分析などが行われています
• また、国よりも細かい単位（都道府県や市区町村）も他のgeojsonファイル（シェープファイル）をダウンロードしてくれば解析を行うことができます

補足

• ここまで夜間光データを紹介してきましたが、夜間光データにもいくつかの問題点も報告されています。
• 例えば、一定機関ごとにセンサータイプが異なる（データ名のヘッダーのF10, F12等）ため、それぞれ微妙に異なるセンサーバイアスが発生しています。そのため、長期トレンドでの光量を解析する際は、それらのバイアスを除去するための前処理が必要になります（Calibrationという処理で多くの研究論文が発表されています）。
• 他にも今回のデータは夜間光量を0~63の整数値で表しており、取り扱いが簡単な反面、光量がかなり大きい地点では飽和現象（Saturation）が発生して、正しく光量を測定できてないという性質も持っています。
• また、今回はDMSP-OLSという衛星データセットを紹介しましたが、こちらの衛星ミッションは2013年で終えており、現在は後継機としてNPP VIIRSというより高精度のセンサーを持った衛星が運用されています。
• 衛星データセットも万能ではないので、利用時には先行研究を参照して、この様な諸問題を把握したうえで解析に用いる必要があります。

最後に

• 本記事ではGoogle Colab環境で衛星データセットを取り扱う方法を紹介しました。
• 今回は夜間光データセットを紹介しましたが、世の中には様々な衛星データセットがオープンデータとして公開されており、基本的には今回と同様の方法で取り扱うことができます。
• 各企業単位で収集・蓄積しているミクロ側のビッグデータと衛星データセットの様なマクロ側のビッグデータを組み合わせて解析することで、新たなビッグデータの可能性が生まれるかと思います。
• せっかくオープンデータとして公開されている多様な衛星データセットが、もっと社会のために使われていけばと思っています。