衛星から得られた夜間光データ画像から主要都市を特定してみた

#はじめに

1. Google earth engineから衛星データ(夜間光)をダウンロード
2. numpyやOpenCVで扱える形に加工
3. 得られたデータからmeanshiftを用いて、光が集中している部分を特定

といったことをやりました。

#Google Earth Engineとは
色々な方が解説している内容ではありますが、
GEE(Google Earth Engine)とは、フリー(少なくとも商用でなければ)で
衛星画像のアクセスできるプラットフォームです。

javascriptもしくはpythonによるプラットフォームが提供されているようですが、
各種データの取扱いを楽に行いたいため、ここではpythonを用いました。
なお、先例に倣い、各種コードをGoogle Colab上で動かしています。

#実装

事前準備

Google Earth Engine 公式サイトの右上にあるSign upから登録を済ませておきましょう.
団体(Institution)を入力する部分があるのですが、
国信大とかでも通るようなので、厳密な審査があるわけではないようです。

初期設定

まずは、Colab上で次の操作を実行して
衛星画像を入れるディレクトリを用意しておきましょう。

!mkdir drive/MyDrive/eisei

次に、GEEを利用できるようにします。
と入っても、以下のセルを実行した後に、
コンソールから出力される指示に従って貰えれば結構です。

import ee

ee.Authenticate()
ee.Initialize()

イメージのダウンロード

###関数の準備

def download_image(scale = 10000, longitude_range = [128, 148.43], latitude_range=[29.97, 46.12], fig_dir = 'eisei', fig_name = 'test'):
  # 衛星名を指定
  satellite = "NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG"

  # バンド名を指定
  band = 'avg_rad'
  
  # 期間を指定
  from_date='2020-04-01'
  to_date='2020-4-30'

  ## エリアを指定(日本エリアを緯度・経度を指定)
  geometry = ee.Geometry.Rectangle([longitude_range[0], latitude_range[0], longitude_range[1], latitude_range[1]])

  # 指定した条件でGEEからデータをロード
  dataset = ee.ImageCollection(satellite).filter(
    ee.Filter.date(from_date, to_date)).filter(
    ee.Filter.geometry(geometry)).select(band)
  
  ## リスト形式へ変換
  data_list = dataset.toList(dataset.size().getInfo())
  
  # 0番目の画像を取得
  image = ee.Image(data_list.get(0))
  
  # Gdriveへ保存
  task = ee.batch.Export.image.toDrive(**{
    'image': image,# 対象データの指定
    'description': fig_name,# ファイル名の指定
    'folder':fig_dir,# Google Drive(MyDrive)のフォルダ名
    'scale': scale,# 解像度の指定
    'region': geometry.getInfo()['coordinates']# 上記で指定した対象エリア
  })

  task.start()

  return task

はい、ぶっちゃけ、コメント含めてTakumi Bannai氏のコードのほぼパクリです(記事はこちら)。
ただ、解像度(scale)、緯度経度の範囲(latitude_range, longitude_range)、
ディレクトリ名(fig_dir)、画像名(fig_name)については、
後々変更していきたいのでこのような形にしました。
ちなみにここでは「解像度」と書かれていますが、
この値を小さくすればするほど、画像が鮮明になる代わりに、画像のサイズが重くなります。

###画像のダウンロード
先ほどの関数を呼び出すことで、作ったディレクトリに衛星からの画像(.tif形式)が出力されます。
ここでは、関東周辺(test_jp0)とほぼアジア全体(test_wide)の２つを用意しています。

longitude_range_jp0 = [136.5, 140.5]
latitude_range_jp0 =[34.5, 38.5]

longitude_range_wide = [110, 150]
latitude_range_wide =[10, 50]

task1 = download_image(scale = 1000, longitude_range = longitude_range_jp0, latitude_range=latitude_range_jp0, fig_name = 'test_jp0')
task2 = download_image(scale = 4000, longitude_range = longitude_range_wide, latitude_range=latitude_range_wide, fig_name = 'test_wide')

なお、Google Colab(python上)では、
あくまで操作のリクエストをGEE側に投げているだけなので、
Colab上でセルの実行が完了したからと言って、
画像がダウンロードされるというわけではありません。
タスクの状況を確認したい場合にはtask.active()もしくはtask.status()から見ることができるようです。

#.active()でGEE側の状況を確認
print(task1.active())
#True 未完了
#False 終了済

#.status()なら詳細な情報を表示可能
print(task1.status())

'''
{'attempt': 1,
 'creation_timestamp_ms': 1615546592726,
 'description': 'test_jp',
 'destination_uris': ['https://drive.google.com/#folders/1qqpdlnujTBWPj2klgvlldWobPm6YdfAg'],
 'id': 'TYEGXMHLCGDYLDZ5O3AGNTXT',
 'name': 'projects/earthengine-legacy/operations/TYEGXMHLCGDYLDZ5O3AGNTXT',
 'start_timestamp_ms': 1615546611779,
 'state': 'COMPLETED', #終了時ならCOMPLETED, ダウンロード中ならRUNIINGなど
 'task_type': 'EXPORT_IMAGE',
 'update_timestamp_ms': 1615546655234}
'''

また、task.active()はFalseになっているのにもかかわらず、
ドライブを見ても、画像がダウンロードされていないといったエラーが起きている場合には、
task.status()を確認することで、エラーコードを見ることができます。
###matplotlibを利用した画像の表示
あらかじめライブラリrasterioを環境にインストール並びにインポートしたうえで、

!pip install rasterio #Colabには入っていないので、pipで入れておく
import rasterio
import matplotlib.pyplot as plt

次のようにすれば簡単に中身を表示できるようです。(ここもパクリ)

path = '/content/drive/MyDrive/{}/{}.tif'.format(fig_dir, fig_name)
with rasterio.open(path) as src:
    arr = src.read()

plt.imshow(arr[0])

この際、少し見づらいと感じたら、
plt.imshow()時に変数(vmin, vmax)を指定することで解決できます。

path = '/content/drive/MyDrive/{}/{}.tif'.format(fig_dir, fig_name)
with rasterio.open(path) as src:
    arr = src.read()

plt.imshow(arr[0], vmin = 0, vmax = 10)

###衛星とその画像を指定するパラメータについて(補足)

今回は衛星名(satelliteに対応)とバンド名(bandに対応)を指定していますが、
これを適切なものに変えることで、
色々なデータ画像(夜間光に限らず)(ページはこちら)を取得することができるようです。
夜間光に限定すれば、使えるデータとして次のようなものがあります。

(2021年3月現在、GEEサイトより)
ここで、今回利用した衛星(一番右)について開くとこのようなページが出てきます。

ここにある「Earth Engine Snippet」で示される文字列を入れてあげれば、
衛星を指定することができるようです。
次にバンド名(衛星からどんなデータを取りたいのか)については、
少し下にスクロールして、ページ下部にあるBandsタブをクリックすることで一覧が表示されます。

ここでは、「avg_rad」と「cf_cvg」の二つが指定できることが見て取れます。
興味があれば試してみてください。

##ダウンロードしたデータから主要都市(光が集まっている部分の中心)を表示
まず、次のような5つの関数を用意します。

• combine_to_array：画像データをnumpy arrayに変換する作業（ほぼパクリ）
• meanshift_core：meanshiftにより光っている画素をクラスタリング
• draw_plane_image：ただの白黒画像を出力
• get_cluster_center_locations_dummy：画像内の位置を緯度経度にそのまま変換
• draw_marked_image：光っている場所+集合の中心+それに相当する緯度経度を表示

from sklearn.cluster import MeanShift, estimate_bandwidth

#画像データをアレーに変換
def combine_to_array(fig_dir, fig_name, index = 0):
  path = '/content/drive/MyDrive/{}/{}.tif'.format(fig_dir, fig_name)
  with rasterio.open(path) as src:
    arr = src.read()
  
  plt.imshow(arr[0], vmin = 0, vmax = 10)
  return arr[index]

#クラスタリングの実行およびその準備
def meanshift_core(arr, vmin,bandwidth_c):
  coordinates = []

  for i in range(len(arr)):
    for j in range(len(arr[0])):
      v = arr[i][j]
      if v > vmin:
        #左下を原点にしたうえで,[経度(横), 緯度(縦)]の順に直す
        coordinates.append([j,len(arr) - i])
  
  X = np.array(coordinates)

  #参考用の"bandwidth"を取得(ここら辺は適当)
  bandwidth = estimate_bandwidth(X, quantile=0.2)

  #ミーンシフトの実行
  clustering = MeanShift(bandwidth=bandwidth/bandwidth_c).fit(X)

  return clustering

#夜間光のみ描画
def draw_plane_image(arr, vmin):
  arr_np = np.array(arr, dtype=np.uint8)

  #値がvmin以上なら白、以下なら黒に変換
  arr_np = [[v*300 if v >= vmin else 0 for v in arr_np_] for arr_np_ in arr_np]

  #cv2で利用できる形にアレーの形状を変更
  img = np.tile(arr_np, (3, 1, 1)) 
  img = np.transpose(img, (1,2,0))

  return img

#画素数データから緯度経度に無理やり変換
def get_cluster_center_locations_dummy(clustering, img, longitude_range, latitude_range):
  centers = clustering.cluster_centers_

  img_height = img.shape[0]
  img_width = img.shape[1]

  #座標を移し替える関数
  def transform_coordinates(x, xmin, xmax, ymin, ymax):
    return ymin + (ymax - ymin) / (xmax - xmin) * (x - xmin)

  locations = [[transform_coordinates(x = center[0], xmin = 0, xmax = img_width, ymin = longitude_range[0], ymax = longitude_range[1]),
              transform_coordinates(x= center[1], xmin = 0, xmax = img_height, ymin = latitude_range[0], ymax= latitude_range[1])] for center in centers]

  return np.array(locations)

#plane_imageの上にマーカーを表示
def draw_marked_image(im, clustering, locations):
  im2 = im.copy()
  
  d = 10 #マーカー円の半径
  cluster_centers_int = np.array(clustering.cluster_centers_, dtype = np.uint32) #整数値に変換 
  for loc, cc_center_int in zip(locations, cluster_centers_int):

    #位置の指定方法に注意
    image_point = (cc_center_int[0], im.shape[0] - cc_center_int[1])

    #マーカーの描画
    im2 = cv2.circle(im2,image_point, d, (120,120,255), 2)
    
    #テキストの記述
    cv2.putText(im2, '{:.2f}, {:.2f}'.format(loc[0],loc[1]),image_point, cv2.FONT_HERSHEY_SIMPLEX, 0.3, (255, 255, 255), 
                )
  
  #numpyアレー形式に戻す
  im2 = im2.get()

  return im2

なお、transform_coordinates(get_cluster_center_locations_dummy内)では、
こんな感じのことをやっています。

次に、これらを組み合わせて夜間光が集まっている場所を出力する関数を作成します。

def get_image_and_clustering(bandwidth_c = 20, vmin = 20, longitude_range = [68.60, 148.43], latitude_range=[-29.97, 46.12],
              fig_dir = 'eisei', fig_name = 'test'):
  
  arr = combine_to_array(fig_dir, fig_name, index = 0)
  _, clustering = meanshift_core(arr, vmin, bandwidth_c)
  img = draw_plane_image(arr, vmin)
  
  locations_dummy = get_cluster_center_locations_dummy(clustering, img, longitude_range, latitude_range)
  img2 = draw_marked_image(img, clustering, locations_dummy)

  cv2_imshow(img2)

  return img, clustering

この関数を実行すると、次のような感じの画像が出力されます。

img_jp0, locations_jp0 = get_image_and_clustering(
    bandwidth_c = 5,vmin = 50, longitude_range = longitude_range_jp0, latitude_range = latitude_range_jp0, fig_name ='test_jp0')

これだけだとよくわからないので、
参考までにplt.imshow()したものも載せておきます。

光が集中している場所が取得できている上、
それらが都市となんとなく対応していることが見て取れます。
また、緯度経度について確認すると、
例えば東京の場合には、(139.75, 35.67)を指しており(見づらいですが)、
マップで検索すると、一応それっぽい場所を指しているようです。

(画像はGoogleマップ)
ただ、画像の範囲が狭いうちはまだ良いのですが、
指定する緯度、経度の範囲を広くすると状況は一変します。
今度は、アジア全体について同様に実行すると、

img_w, locations_w = get_image_and_locations(
    bandwidth_c = 20,vmin = 60, longitude_range = longitude_range_wide, latitude_range = latitude_range_wide, fig_name ='test_wide')

!

となり、点群の中心は得られているものの、緯度経度に関する情報が若干怪しくなります。
先ほどと同様に、東京っぽい場所の緯度経度を検索すると、

このように、先ほどの例と比べても大きく異なった位置を指すことがわかります。
こういった問題を補正するのが、衛星画像の幾何補正と呼ばれるものです。

ここまでのまとめ

• GEE+Google Colab を用いて取得した夜間光の衛星画像から光が集まっている位置を取得
• 光が集中している場所の抽出にはmeanshiftを利用、一応成功
• 位置に関してはそのまま移しただけではダメで別の補正が必要