Landsat-8の衛星データを取得する

はじめに

Landsat-8は広域撮影を目的とした光学センサで撮影した解像度15mの白黒画像と30mのカラー画像です。Landsat-8は地球全体の中分解能マルチスペクトル画像だけでなく、熱赤外画像の取得・保存をしており、そのデータは誰でも自由に利用することができます。APIの使い方はリファレンスにまとめられています。 ※画像を利用する際はのデータポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/api/v1/landsat8/{min_lat}/{min_lon}/{max_lat}/{max_lon}.png
Name Input Description
min_latnumber(query) 必須 最小緯度(-90~90)
min_lonnumber(query) 必須 最小経度(-180~180)
max_latnumber(query) 必須 最大緯度(-90~90)
max_lonnumber(query) 必須 最大経度(-180~180)
cloudnumber(query) 任意 雲量

最初に衛星シーンデータ(衛星画像のメタデータ)を取得する関数を定義します。APIに必要な値は、ASNAROと同じく、最小の経度緯度と最大の経度緯度(min_lat, min_lon, max_lat, max_lon)です。この範囲にあるシーンを検索し、データはjsonフォーマットで取得しています。

import os, json, requests, math, glob, re
import numpy as np
import pandas as pd
from skimage import io
from io import BytesIO
import dateutil.parser
from datetime import datetime
from datetime import timezone
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "自分のトークンを貼る"
def get_landsat_scene(min_lat, min_lon, max_lat, max_lon):
    """Landsat8のメタデータを取得します"""
    url = "https://gisapi.tellusxdp.com/api/v1/landsat8/scene" \
        + "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return r.json()

シーン検索

今回は種子島宇宙センターがある鹿児島周辺のシーンを検索します。

# 鹿児島周辺
scenes = get_landsat_scene(30.448673679287566, 130.78125, 30.90222470517145, 130.78125)
print(len(scenes))

58箇所のシーンデータが取得できました。種子島は30.6096° N, 130.9789° Eのあたりに位置しています。種子島の座標情報を基にして、シーンデータに含まれる座標情報を取捨選択します。そのため、シーンデータから座標情報を抽出する関数を定義しましょう。

def ext_latlon(scene_data):
    max_lats = []
    max_lons = []
    min_lats = []
    min_lons = []
    range_latlon = []
    for i in range(len(scene_data)):
        max_lats = scene_data[i]['max_lat']
        max_lons = scene_data[i]['max_lon']
        min_lats = scene_data[i]['min_lat']
        min_lons = scene_data[i]['min_lon']
        range_latlon.append([max_lats, max_lons, min_lats, min_lons])
    return np.array(range_latlon)

配列をデータフレームに変換します。

range_latlons = ext_latlon(scenes)
column_names = ['max_lat', 'max_lon', 'min_lat','min_lon']
df = pd.DataFrame(data = range_latlons, columns = column_names, dtype='float')
print(df)

devLandsat_api1_20200220_5.png このデータから種子島を含む範囲の行だけを抽出します。

df_ext = (df[(df['max_lat'] < 32.0) & (df['max_lon'] < 132.0) & (df['min_lat'] > 29.0) & (df['min_lon'] > 129.0)])
df_ext

devLandsat_api1_20200220_6.png やや広い範囲で条件設定しましたが、58シーンから15シーンのみが抽出されました。元のシーンデータから条件にあったシーンを抽出するためにデータフレームからインデックス番号を取得しましょう。

index_li = list(df_ext.index)
print(index_li)

[2, 10, 12, 14, 16, 22, 24, 30, 36, 40, 41, 52, 53, 55, 56]のように結果が返ってきます。

ext_scenes = []
for i in index_li:
    ext_scenes.append(scenes[i])
print(ext_scenes)

devLandsat_api1_20200220_7.png インデックスを基にして、シーンからデータが抽出されました。選択したシーンデータに種子島が含まれているのか、サムネイル画像一覧を描写して確かめてみましょう。

ext_thumbs = [imgs['thumbs_url'] for imgs in ext_scenes]
print(ext_thumbs)

devLandsat_api1_20200220_8.png 続いて、サムネイル画像を格子状に表示する関数を定義します。

def get_tile_images(thumbs, r_num, c_num, figsize_x=12, figsize_y=12):
        """サムネイル画像を指定した引数に応じて行列状に表示"""
        f, ax_list = plt.subplots(r_num, c_num, figsize=(figsize_x,figsize_y))
        for row_num, ax_row in enumerate(ax_list):
            for col_num, ax in enumerate(ax_row):
                if len(thumbs) < r_num*c_num:
                    len_shortage = r_num*c_num - len(thumbs) # 行列の不足分を算出
                    count = row_num * c_num + col_num
                    if count < len(thumbs):
                        ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
                    else:
                        for i in range(len_shortage):
                            blank = np.zeros([100,100,3],dtype=np.uint8)
                            blank.fill(255)
                            ax.imshow(blank)
                else:
                    ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
        return plt.show()
get_tile_images(ext_thumbs,5,3,16,16)

devLandsat_api1_20200220_1.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

サムネイルに種子島と屋久島が含まれていることがわかります。島に雲がかかっているのを除けば、どの画像を選んでも種子島宇宙センターが映り込む地図タイル画像を取得できそうです。

地図タイル画像の取得

Landsat-8の地図タイル画像を取得するAPIにはproductIdを指定する必要はありません。一方で、波長合成画像を取得するAPIにはproductIdの指定が必要です。地図タイル画像はシーンを選ぶことができませんが、より高解像度の画像を取得することができます(Landsat-8の波長合成画像ではズーム率は12までしか選択できません)。 シーンを設定した波長合成画像を取得する前に、地図タイル画像で宇宙センターの位置を確認してみましょう。ASNARO-1ではentityIdを含める必要がありましたが、こちらは(zoom, xtile, ytile)を値として入れるだけになっています。https://gisapi.tellusxdp.com/landsat8/{}/{}/{}.pngの{}に各々の値が入ります。

# Landsatではタイルデータを取得する場合には、シーンidは指定しない
def get_landsat_image(zoom, xtile, ytile):
    """取得した画像のタイルを定義し、さらに拡大して表示"""
    url = "https://gisapi.tellusxdp.com/landsat8/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }    
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))
tanegashima_south = get_landsat_image(13, 7076, 3369)
io.imshow(tanegashima_south)

devLandsat_api1_20200220_4.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

画像の左隅が宇宙センターが位置する場所になります。

Name Description
path Landsat8のPATH番号 必須項目 デフォルトなし
row Landsat8のROW番号 必須項目 デフォルトなし
productid Landsat8のシーン番号 必須項目 デフォルトなし
z タイルマップのZ座標 必須項目 デフォルトなし ( 6 ~ 12 )
x タイルマップのX座標 必須項目 デフォルトなし
y タイルマップのY座標 必須項目 デフォルトなし
preset true / false / natural / ndvi / urban / agli / atmos / healthy / landw / naturalwoatom / shortir / vegitation のいずれかを指定 デフォルトなし
opacity 透過度 0で完全に透過、1で不透過 デフォルト1.00 ( 0.00 ~ 1.00 )
r RGBカラーのRに入るバンドの番号を指定 デフォルトr,g,b = 4,3,2 ( 1,2,3,4,5,6,7,9,10,11 )
g RGBカラーのGに入るバンドの番号を指定 デフォルトr,g,b = 4,3,2 ( 1,2,3,4,5,6,7,9,10,11 )
b RGBカラーのBに入るバンドの番号を指定 デフォルトr,g,b = 4,3,2 ( 1,2,3,4,5,6,7,9,10,11 )
rdepth Rの濃さを指定 デフォルト1.0 ( 0.00 ~ 2.00 )
gdepth Gの濃さを指定 デフォルト1.0 ( 0.00 ~ 2.00 )
bdepth Bの濃さを指定 デフォルト1.0 ( 0.00 ~ 2.00 )

最後に波長合成画像を取得しましょう。波長合成画像のAPIはASNARO-1のものより多くの変数を割り当てる必要があります。https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.pngのカッコ部分に入る値は、それぞれ、(path, row, productId, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth)となります。関数を定義して、デフォルトの値を割り当てておけば、全ての変数に値を代入する必要はありません。波長合成画像APIではさらにpresetが必要になります。下で定義した関数では、初期値でNoneが与えられています。また、RGBにはTrue color画像を取得できる初期値が与えられています。

def get_landsat_blend(path,row,productId,zoom, xtile, ytile, opacity=1, r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
    url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
    (path, row, productId, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth)
    if preset is not None:
        url += '&preset='+preset
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

抽出したシーンの後ろから三番目のデータを選択しました。今回は二つの画像を水平方向につなぎ合わせた画像を出力します。

low_cloud_scene = ext_scenes[-3]
low_cloud_scene_img1 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3537, 1684)
low_cloud_scene_img2 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3538, 1684)
io.imshow(np.hstack((low_cloud_scene_img1,low_cloud_scene_img2)))

devLandsat_api1_20200220_3.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

画像の真ん中下辺りに宇宙センターの位置する場所があります。続いてFalse color画像を取得してみましょう。

low_cloud_scene_false1 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3537, 1684, preset = "false")
low_cloud_scene_false2 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3538, 1684, preset = "false")
io.imshow(np.hstack((low_cloud_scene_false1,low_cloud_scene_false2)))

devLandsat_api1_20200220_2.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

波長合成画像ではpresetを用いる方法もありますが、開発者が自分で好きなバンドを当てはめた画像を取得することもできます。例えば、RGB値に衛星のバンドを直接割り当ててFalse color画像を作成してみましょう。

b_low_cloud_scene_false1 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3537, 1684,r = 5, g = 4, b = 3)
b_low_cloud_scene_false2 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3538, 1684, r = 5, g = 4, b = 3)
io.imshow(np.hstack((b_low_cloud_scene_false1,b_low_cloud_scene_false2)))

devLandsat_api1_20200220_2.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

同様の画像が取得できました。このように波長合成画像取得APIでは任意に衛星のバンドをRGB値に当てはめることができるため、目的に応じた衛星画像を取得することが可能となります。 波長合成画像を使いこなして、衛星データで様々な解析に応用してみてください。

今回使用したスクリプトです。

import os, json, requests, math, glob, re
import numpy as np
import pandas as pd
from skimage import io
from io import BytesIO
import dateutil.parser
from datetime import datetime
from datetime import timezone
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここにトークンを貼る"
def get_landsat_scene(min_lat, min_lon, max_lat, max_lon):
    """Landsat8のメタデータを取得します"""
    url = "https://gisapi.tellusxdp.com/api/v1/landsat8/scene" \
        + "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return r.json()
# 鹿児島周辺
scenes = get_landsat_scene(30.448673679287566, 130.78125, 30.90222470517145, 130.78125)
# {\min_lat\:40.49,\min_lon\:140.502,\max_lat\:40.744,\max_lon\:140.797}
print(len(scenes))
def ext_latlon(scene_data):
    max_lats = []
    max_lons = []
    min_lats = []
    min_lons = []
    range_latlon = []
    for i in range(len(scene_data)):
        max_lats = scene_data[i]['max_lat']
        max_lons = scene_data[i]['max_lon']
        min_lats = scene_data[i]['min_lat']
        min_lons = scene_data[i]['min_lon']
        range_latlon.append([max_lats, max_lons, min_lats, min_lons])
    return np.array(range_latlon)
range_latlons = ext_latlon(scenes)
column_names = ['max_lat', 'max_lon', 'min_lat','min_lon']
df = pd.DataFrame(data = range_latlons, columns = column_names, dtype='float')
print(df)
df_ext = (df[(df['max_lat'] < 32.0) & (df['max_lon'] < 132.0) & (df['min_lat'] > 29.0) & (df['min_lon'] > 129.0)])
df_ext
index_li = list(df_ext.index)
print(index_li)
xt_scenes = []
for i in index_li:
    ext_scenes.append(scenes[i])
print(ext_scenes)
ext_thumbs = [imgs['thumbs_url'] for imgs in ext_scenes]
print(ext_thumbs)
def get_tile_images(thumbs, r_num, c_num, figsize_x=12, figsize_y=12):
        """サムネイル画像を指定した引数に応じて行列状に表示"""
        f, ax_list = plt.subplots(r_num, c_num, figsize=(figsize_x,figsize_y))
        for row_num, ax_row in enumerate(ax_list):
            for col_num, ax in enumerate(ax_row):
                if len(thumbs) < r_num*c_num:
                    len_shortage = r_num*c_num - len(thumbs) # 行列の不足分を算出
                    count = row_num * c_num + col_num
                    if count < len(thumbs):
                        ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
                    else:
                        for i in range(len_shortage):
                            blank = np.zeros([100,100,3],dtype=np.uint8)
                            blank.fill(255)
                            ax.imshow(blank)
                else:
                    ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
        return plt.show()
get_tile_images(ext_thumbs,5,3,16,16)
# Landsatではタイルデータを取得する場合には、シーンidは指定しない
def get_landsat_image(zoom, xtile, ytile):
    """取得した画像のタイルを定義し、さらに拡大して表示"""
    url = "https://gisapi.tellusxdp.com/landsat8/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }    
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))
tanegashima_south = get_landsat_image(13, 7076, 3369)
io.imshow(tanegashima_south)
def get_landsat_blend(path,row,productId,zoom, xtile, ytile, opacity=1, r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
    url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
    (path, row, productId, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth)
    if preset is not None:
        url += '&preset='+preset
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))
low_cloud_scene = ext_scenes[-3]
low_cloud_scene_img1 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3537, 1684)
low_cloud_scene_img2 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3538, 1684)
io.imshow(np.hstack((low_cloud_scene_img1,low_cloud_scene_img2)))
low_cloud_scene_false1 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3537, 1684, preset = "false")
low_cloud_scene_false2 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3538, 1684, preset = "false")
io.imshow(np.hstack((low_cloud_scene_false1,low_cloud_scene_false2)))
b_low_cloud_scene_false1 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3537, 1684,r = 5, g = 4, b = 3)
b_low_cloud_scene_false2 = get_landsat_blend(low_cloud_scene['path'],low_cloud_scene['row'],\
                                         low_cloud_scene['productId'],12, 3538, 1684, r = 5, g = 4, b = 3)
io.imshow(np.hstack((b_low_cloud_scene_false1,b_low_cloud_scene_false2)))