記事作成時から、Tellusからデータを検索・取得するAPIが変更になっております。該当箇所のコードについては、以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/access/traveler_api_20220310_
firstpart.html

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

本記事の内容

Tellusを使って、衛星画像から特定植物の抽出をした際の処理をまとめました。
衛星画像解析の考え方は以下の記事で解説しています。記事を読んでいただくとより理解が深まると思いますのでぜひご覧ください。

Tellusにて、画像解析する際の全体像は以下の通りです。

事前準備
(1) Tellusの開発アカウントを取得
(2) Tellusのマイページにて、apiアクセストークンを発行
(3) 抽出対象の植物分布データ(サンプル)を用意
(4) Tellus OS上でタイル座標を確認
(5) Tellus OS上で衛星画像を検索し、product idを確認
(6) 衛星画像内で植物がいる場所を型取ったファイルをTellus上にアップロード

衛星画像から特定植物の抽出手順
(1) Tellus開発環境上で、衛星画像を取り込む
(2) Tellus上のgeojsonファイルを取得し、geojsonファイルから色の範囲を取得
(3) 取得した色の範囲に存在する色を白色に変更

以降、それぞれの詳細について説明していきます。

事前準備

⑴ Tellusの開発アカウントを取得

1. Tellusのログイン画面からTellus OSにログインします。

2. 右上部の自身のユーザ名をクリック後、ダッシュボードから『開発環境の申し込み』をクリック。入力フォームの内容を埋め送信ボタンをクリックします。4,5日待つと登録したメールアドレス宛に開発環境の利用可能通知が届きます。※空き状況によってかかる時間は異なります。

3. ダッシュボードの開発環境利用状況を確認し、トークン情報を確認します。このトークンは開発環境にログインする際に必要なります。

⑵ Tellusのダッシュボードにて、APIアクセストークンを発行

1. 開発環境内のAPIアクセスにアクセスし、任意のAPIキー名を決めてアクセストークンを作成します。このアクセストークンはTellusのAPIを使用する際に必要になります。

⑶ 抽出対象の植物分布データ(サンプル)を用意

本記事ではアカエゾマツの抽出を行う為アカエゾマツの分布がわかるデータを使います。

1. アカエゾマツの分布図がわかるサンプルを探します。google chrome上で「アカエゾマツ 森林」と検索すると下図のアカエゾマツの分布図が見つかりました。

⑷ Tellus OS上でタイル座標を確認

上記のサンプルが示す部分と一致する場所をTellus OS上で確認します。

1. Tellus OSにログイン後、上部の?をクリックし、調べたい住所を入力します。

2. 画面上部の格子状のアイコンをクリック後、タイル座標をクリックします。今回はLandsat-8という衛星の衛星画像を使う為、Landsat-8の最大ズーム率であるz = 12の時の、x, yの値を確認します。サンプルの示す部分と被っているタイル(青枠で囲まれた格子)は、x = 3691, y = 1495です。

⑸ Tellus OS上で衛星画像を検索し、product idを確認

1. 左側のマイライブラリをクリック後、マイライブラリの中のLandsat-8の左側にある目のマークをクリックします。表示している画面上に存在するLandsat-8のシーンが表示されます。

2. 右下のタイムスケールで時期を入力します。今回は春と秋の衛星画像を取得したいので、春は 2018/5/1を指定し、秋は 2018/10/1を指定しました。

表示スケールを1dayとし、欲しい月のシーンがあるかを確認します。
シーンがある場合は、タイムスケール上に赤い丸で表示されます。

3. マイライブラリの詳細から画像をクリック後、該当の衛星データをクリックしproductidを確認します。本記事では、春の画像のproductidとして「LC08_L1TP_105030_20180512_20180517_01_T1」を利用し、秋の画像のproductidとして「LC08_L1TP_105030_20181003_20181010_01_T1」を利用します。

4. 念のため、2.で選択した画像にx= 3691, y=1495が含まれているかを確認します。

⑹衛星画像内で植物がいる場所を型取りしたファイルを読み込み

1. 左部の『マイライブラリ』＞『インポート』から下図青色枠内のアイコンをクリック後、取り込みたいgeojsonファイルを選択し、Tellus上に取り込みます。

本記事で使用したアカエゾマツのgeojsonファイル

Tellus上でのgeojsonファイル詳しい取り扱いについてはこちらをご覧ください。

取り込むと下図のように表示されます。

衛星画像から特定植物の抽出手順

事前準備が完了したので、ここからは実際にTellus開発環境上で衛星画像を取得し、植物の分布を抽出します。また今回は衛星画像を使って、アカエゾマツの分布を抽出します。アカエゾマツの分布を抽出するにあたっての概念はこちらの記事をご確認ください。また末尾で今回使用したプログラムを記載します。

⑴Tellus開発環境上で、衛星画像を取り込む

1. Tellus の右上部のユーザーIDをクリックし、開発環境にアクセスします。この時tokenの入力を求められるので、開発環境利用状況のトークンを入力します。

2. 開発環境にログインし、jupyter notebookファイルを開きます。
3. 下記ライブラリをインポートし、APIトークンを変数でしておきます。

import requests
import json
import math
from skimage import io
from io import BytesIO
import numpy as np
from skimage.draw import polygon
　　　
#APIアクセストークン
api_token = ''

4. 事前準備で確認した設定値を変数に代入します。また今回使用する衛星画像は、Landsat-8が撮影したものを使用するので、sat_nameには、Landsat-8を代入します。

#ズーム率(6 ~ 12)
zoom_rate = 12
#探したい地点のタイルの座標
(xtile,ytile)=(3691, 1495)

#観測衛星
sat_name='landsat8'

#春の画像のId
spring_product_id='LC08_L1TP_105030_20180512_20180517_01_T1'

#秋の画像のId
fall_product_id='LC08_L1TP_105030_20181003_20181010_01_T1'

5. アカエゾマツの分布を抽出するには、下表のBand3, Band4, Band5を観測する必要がある為、paramsに下記の値を設定します。paramsの設定値については、get_blend_img関数を説明する際に詳しく説明します。

#観測するバンド
params={'r':5, 'g':4, 'b':3}

6. get_scene関数でproduct_idに該当するsceneを取得し、get_blend_img関数でsceneに応じたブレンド画像(今回はband5,band4,band3のブレンド画像)を取得します。

#春の画像を取得
scene = get_scene(sat_name,spring_product_id)  #該当するシーンの取得
img_spring = get_blend_img(scene,sat_name,zoom_rate,xtile,ytile,params)
io.imshow(img_spring)
img_spring.shape

取得結果

【get_scene関数】に関して
Tellusのシーン取得API( https://gisapi.tellusxdp.com/api/v1/{}/get_scene/{} )を使って、シーンの取得を実施、引数にはシーンを取得するために必要な衛星名とproductidを指定しています。シーン取得APIは、Landsat-8のAPIリファレンスページで仕様を確認できます。

def get_scene(sat_name,product_id):
    """
    sceneを取得する
    Parameters
    ----------
    sat_name : str
        観測衛星名
    product_id : int 
        衛星画像のプロダクトID
    Returns
    -------
    url: json
    　　観測画像の情報
    """
    url = 'https://gisapi.tellusxdp.com/api/v1/{}/get_scene/{}'.format(
        sat_name, product_id)
    headers = {
        'Authorization': 'Bearer ' + api_token
    }
    
    r = requests.get(url, headers=headers)
    
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return r.json()

【get_blend_img関数】に関して
Tellusのバンド合成画像取得API( https://gisapi.tellusxdp.com/blend/{}/{}/{}/{}/{}/{}/{}.png )を使って、画像を取得します。引数paramsの‘r’, ‘g’, ‘b’に対して取得したいバンド値を指定することで特定バンドで合成された画像を取得することができます。バンド合成画像取得APIは、Landsat-8のAPIリファレンスページで仕様を確認できます。

def get_blend_img(scene, sat_name, zoom_rate, xtile, ytile, params={'preset':None, 'opacity':1, 'r':4, 'g':3, 'b':2, 'rdepth':1, 'gdepth':1, 'bdepth':1}):
    """
    Landsat-8のブレンドタイル画像を取得する
    Parameters
    ----------
    scene : dict
        観測対象の衛星画像
    sat_name : str
        観測衛星名
    zoom_rate : int 
        タイルのズーム率 
    xtile : int 
        タイルのX座標 
    ytile : int 
        タイルのY座標 
    params : dict 
        クエリパラメータ
    Returns
    -------
    img: ndarray
        タイル画像
    """
    url = 'https://gisapi.tellusxdp.com/blend/{}/{}/{}/{}/{}/{}/{}.png'.format(
        sat_name, scene['path'], scene['row'], scene['productId'], zoom_rate, xtile, ytile)
    headers = {
        'Authorization': 'Bearer ' + api_token
    }
    r = requests.get(url, params=params, headers=headers)
    
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return io.imread(BytesIO(r.content)) 

⑵Tellus OS上のgeojsonファイルを取得し、geojsonファイルから色の範囲を取得

1. Tellus OS上にアップロードしたgeojsonファイルを取得します。

filename = "akaezomatsu"

#自身のgeojsonファイルを取得
files = get_geojson_file_list()

for file in files["items"]:
    if file["originalName"]==filename:
        file_id = file["url"]
        
#ここで使うファイルを呼び出す
sample_polygon = get_geojson_file(file_id)

#読み込んだgeojsonファイルの改行コードを省く為、dict型に変換。
shape =  = requests.get(file_id).json()

【get_geojson_file_list関数】に関して
Tellusのファイルリスト取得API( https://api.tellusxdp.com/v1/geojson-files )を使って、Tellusにアップロードされているgeojsonファイルの一覧を取得します。ファイルリスト取得APIは、Tellus User APIリファレンスページで仕様を確認できます。

#Tellus OSにあげた全Geojsonファイルを読み込む
def get_geojson_file_list():
    url = "https://api.tellusxdp.com/v1/geojson-files/"
    headers = {
        "Authorization": "Bearer " + api_token
    }
    
    r = requests.get(url, headers=headers)
    return r.json()

【get_geojson_file関数】に関して
Tellusのファイル取得API( https://api.tellusxdp.com/v1/geojson-files/{} )を使って、Tellusにアップロードされているgeojsonファイルを取得します。引数にファイルIDを指定することで特定のファイルを取得することができます。ファイル取得APIは、Tellus User APIリファレンスページで仕様を確認できます。

# ファイルIDに従ったファイルを読み込む
def get_geojson_file(file_id):
    url = "https://api.tellusxdp.com/v1/geojson-files/{}".format(file_id)
    headers = {
        "Authorization": "Bearer " + api_token
    }
    
    r = requests.get(url, headers=headers)
    return r.json()

2. 取得したgeojsonファイルの座標を元に、クリッピングするためのクリッピング枠を作成します。

points = shape["features"][0]["geometry"]["coordinates"][0]
bbox = get_tile_bbox( xtile, ytile, zoom_rate)
pxsize = {"width": img_fall.shape[1], "height": img_fall.shape[0]}
mask = get_polygon_image(points, bbox, pxsize)
#クリッピング枠表示
io.imshow(mask[0:140,0:140])

クリッピング枠作成結果

【get_tile_bbox関数】に関して
タイル座標からバウンディングボックスを取得します。バウンディングボックスの四隅の値は緯度経度で算出します。

#バウンディングボックスを取得
def get_tile_bbox(xtile, ytile, zoom):
    """
    タイル座標からバウンディングボックスを取得
    """
    def num2deg(xtile, ytile, zoom):
        # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        
        return (lon_deg, lat_deg)
    
    right_top = num2deg(xtile + 1, ytile, zoom)
    left_bottom = num2deg(xtile, ytile + 1, zoom)
    return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])

【get_polygon_image関数】に関して
geojsonファイルの座標、バウンディングボックスの値、ピクセル値のサイズを元に、クリッピング枠を作成します。

#緯度経度をピクセル座標に変換し、クリッピング枠を作成
def get_polygon_image(points, bbox, imgsize):
    def world2pixel(bbox, size, lon, lat):
        """
        経度緯度をピクセル座標に変換
        """
        dist = ((bbox[2] - bbox[0])/size[0], (bbox[3] - bbox[1])/size[1])
        pixel = (int((lon - bbox[0]) / dist[0]), int((bbox[3] - lat) / dist[1]))
        return pixel
    
    pixels = []
    for p in points:
        pixels.append(world2pixel(bbox, (imgsize["width"], imgsize["height"]), p[0], p[1]))
    pixels = np.array(pixels)

    poly = np.ones((imgsize["height"], imgsize["width"]), dtype=np.uint8)
    (rr, cc) =  polygon(pixels[:, 1], pixels[:, 0], poly.shape)
    poly[rr, cc] = 0
    return poly

3. 2.で作成したクリッピング枠で、クリッピングします。

img_spring_cliped = clip(img_spring, mask)
io.imshow(img_spring_cliped[0:140,0:140])

クリッピング結果

【clip関数】に関して
imgにクリッピング対象を、maskにクリッピング枠を指定し、クリッピングします。

#クリッピング対象とクリッピング枠を指定しクリッピング
def clip(img, mask):
    base = img.transpose(2, 0, 1)[:3,:]
    cliped = np.choose(mask, (base, 0)).astype(np.uint8)
    return cliped.transpose(1, 2, 0)

4. クリッピングした画像の色の範囲を取得します。

spring_max = get_max_color_value(img_spring_cliped)
print('spring_max:', spring_max)
spring_min = get_min_color_value(img_spring_cliped)
print('spring_min:', spring_min)

算出結果

【get_max_color_value関数】に関して
引数imgで受け取った画像の画素値の最大値を取得するので、画素値の最大値である白色[255 255 255]を除外して算出します。

#画素値の最大値を取得
def get_max_color_value(img):
    #縦,横,画素値の3次元行列になっているので、縦,画素値の2次元行列に変換します。
    #この変換を行うことで、画素値の最大を容易に取得できるようにします。
    img_3dto1d = img.flatten() 
    img_1dto2d = img_3dto1d.reshape(img.shape[0]*img.shape[1],img.shape[2])
    
    # [255 255 255]の行列を排除。
    for i in range(len(img_1dto2d)):
        img_1dto2d[i][np.all(img_1dto2d[i] == 255)] = -1
    max_color_value = img_1dto2d.max(axis=0)
    
    return max_color_value

【get_min_color_value関数】に関して
引数imgで受け取った画像の画素値の最小値を取得するので、画素値の最小値である黒色[0 0 0]を除外して算出します。

#画素値の最小値を取得
def get_min_color_value(img):
    #縦,横,画素値の3次元行列になっているので、縦,画素値の2次元行列に変換します。
    #この変換を行うことで、画素値の最小を容易に取得できるようにします。
    img_3dto1d = img.flatten()
    img_1dto2d = img_3dto1d.reshape(img.shape[0]*img.shape[1],img.shape[2])
    # [0 0 0]の行列を排除。
    for i in range(len(img_1dto2d)):
        img_1dto2d[i][np.all(img_1dto2d[i] == 0)] = 500
    min_color_value = img_1dto2d.min(axis=0)
    
    return min_color_value

⑶取得した色の範囲に存在する色を白色に変更

1. 上記で取得した画素値の範囲を元に、範囲内に存在するピクセルを抽出します。今回は範囲内にあるもの白色、範囲外にあるものを黒色にしています。

#春
masked_img_spring = masking(img_spring,spring_max,spring_min)
io.imshow(masked_img_spring[0:140,0:140])

結果

【masking関数】に関して
抽出対象の画像、最大画素値、最小画素値を設定すると、最大最小内に存在するか否かを判別します。

def masking(img, img_max_value, img_min_value):
    masked_img = img.copy()
    end_of_comparison_target = masked_img.shape[2] - 1
    num_in_range = 0
    
    #各ピクセルの画素値が、最大値と最小値の範囲内に存在するか判定しています。
    for i in range(masked_img.shape[0]):
        for j in range(masked_img.shape[1]):
            if (masked_img[i,j,:end_of_comparison_target] = img_min_value).all():
                #範囲内なら白色を設定
                masked_img[i,j,:end_of_comparison_target] = 255
            else:
                #範囲内なら黒色を設定
                masked_img[i,j,:end_of_comparison_target] = 0
    return masked_img

この画像に先ほどのgeojsonファイルが表す外枠を重ねると、実際のアカエゾマツの保護林区域が抽出できているかを確認することができます。

重ねてみると、赤枠で囲まれたエリアとその周辺の地域が抽出されていました。

以下のソースコードでは、geojsonが表す形のアウトラインを取得する為の処理と、「geojsonのアウトライン画像」と「抽出処理後の画像」の二枚の画像を重ねる処理を記載しました。また今回の本題ではないので詳しい説明は割愛しました。

#====以下は、画像のアウトラインを抽出する関数==========
def get_outline(img):
    outline = img.copy()
    flg = 0
    mark1 = 10
    mark2 = 20
    
    for j in range(outline.shape[1]):
        for i in range(outline.shape[0]):
            if (outline[i,j] == flg) and (flg == 0):
                outline[i-1,j] = mark1
                flg = 1
            if (outline[i,j] == flg) and (flg == 1):
                outline[i,j] = mark1
                flg = 0
    flg = 0
    for i in range(outline.shape[0]):
        for j in range(outline.shape[1]):
            if not outline[i,j] == mark1:
                if (flg == 0) and (outline[i,j] == flg):
                    outline[i,j-1] = mark2
                    flg = 1
                if (flg == 1) and (outline[i,j] == flg or outline[i,j] == mark1):
                    if not outline[i,j-1] == mark1:
                        outline[i,j] = mark2
                    flg = 0

    for i in range(outline.shape[0]):
        for j in range(outline.shape[1]):
            if outline[i,j] == 0:
                outline[i,j] = 1
            if outline[i,j] >= mark1:
                outline[i,j] = 0
                
    return outline
#====================================================

#====以下は、枠線の色を変更する関数===================
def change_color(outline_img):
    for i in range(outline_img.shape[0]):
        for j in range(outline_img.shape[1]):
            if not (np.all(outline_img[i,j] == 0)):
                outline_img[i,j] = np.array([255, 0, 0])
    return outline_img
#===================================================

#====以下は、2枚の画像の重ね合わせる関数=============
def put_image(top_img, bottom_img):
    img = bottom_img.copy()
    n = top_img.shape[2] - bottom_img.shape[2]
    topimg_edge = top_img.shape[2]
    bottomimg_edge = bottom_img.shape[2]
    if n > 0:
        topimg_edge -= n
    else:
        bottomimg_edge += n
    
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            if not np.all(top_img[i,j,:topimg_edge] == np.zeros(topimg_edge, dtype=np.int)):
                img[i,j,:bottomimg_edge] = top_img[i,j,:topimg_edge]
    return img
#===================================================

#====以下は、geojsonのアウトラインを取得し、 =======
#====2枚の画像を重ねるまでの処理を記載       =======
mask_outline = get_outline(mask)
outline_img = clip(img_spring, mask_outline)
change_color_img = change_color(outline_img)
io.imshow(change_color_img[0:140,0:140])
change_color_img.shape
mask_img_spring2 = put_image(change_color_img,masked_img_spring)
io.imshow(mask_img_spring2[0:140,0:140])

あとがき

今回の記事では、以前の記事でご紹介したアカエゾマツの抽出について、Tellus上のJupyter Notebookで行う際の解析の手順を1からご紹介しました。この通りにTellusの操作を行えば、難解だと思われがちな衛星データ解析もサクサク進んで楽しめるはず……！

これからもTellusの開発環境で、衛星データ含む様々なデータでできることを徐々に増やし、ご紹介できればと思います。

最後に、今回のプログラムソースをzipファイルにまとめておりますのでぜひご活用ください。

付録：プログラム全体

本記事のプログラムソースはこちら（zipファイル）。