記事作成時から、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のJupyter上で取得した衛星データからGeoJSONで定義した領域を抜き出す（クリッピングする）方法を紹介します。

TellusでJupyter Notebookを使う方法は、「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。

1. 抜き出す領域を指定する

宙畑では以前、衛星データを取得し解析する方法の例（【ゼロからのTellusの使い方】Jupyter NotebookでAVNIR-2/ALOSの光学画像から緑地割合算出）を紹介しました。

この記事では、取得した衛星画像を丸々1枚使って解析していますが、画像全体を利用するのではなく、一部の領域だけを抜き出して解析したくなる場合もありますよね。

例えば農業なら、以下のようなケースです。

・田んぼのA地点とB地点で稲の発育を比較し、品種に適した土壌条件を分析する。
・ある畑の生育状況を過去の同じ地点と比較し、その年の収穫高を予想する。

抜き出す領域を指定するには、位置情報のベクトルデータが必要です。

今回は最も簡単な例として、Tellus OS上で地図や衛星写真を眺めながら、抜き出したい領域の図形を作成し、それを利用します。

では、まず抜き出したい衛星データを選択しましょう。
例では、AVNIR-2のデータを使って東京の晴海周辺を抜きだしてみます。

画面中央に竹芝ふ頭を表示して、左上の「メタデータ検索」をクリックします。

開いたパネルから「光学画像（AVNIR-2）」を選択肢し「検索」をクリックすると、表示範囲を含んだ観測データの一覧が表示されます。

一覧の中から利用したい観測データを選択してOS上にデータが表示しましょう（この例ではALAV2A266502880を選択しています）。

次に抜き出したい領域の図形を作成します。
今回は晴海周辺の埋立地を抜き出してみます。

画面上部のアイコンから「多角形を追加」をクリックし、埋立地の形に沿って図形を作ります。

きちんと図形を閉じたら完成で、できた図形は画面右の「図形リスト」から確認できます。
今回は「polygon-1」という名前で保存されています。

2. Jupyter Notebook上で衛星データと抜き出す領域のGeoJSONを取得する

解析したい衛星データと作成した図形をJupyter Notebookで取得します。

シーンID（この例ではALAV2A266502880）とタイル座標からAVNIR-2のデータを取得するには以下のコードを実行します。

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

TOKEN = "ここには自分のアカウントのトークンを貼り付ける"

def get_AVNIR_SCENEID(scene_id):
    url = "https://gisapi.tellusxdp.com/api/v1/av2ori/get_scene/{}".format(scene_id) 
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get( url, headers=headers )
    return r.json()
    
def get_AVNIR_image(rspid, sceneid, zoom, xtile, ytile, opacity=1, r=3, g=2, b=1, rdepth=1, gdepth=1, bdepth=1, preset=None):
    url = "https://gisapi.tellusxdp.com/av2ori/{}/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format(rspid, sceneid, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth) 
    headers = {
        "Authorization": "Bearer " + TOKEN
    }

    url = "https://gisapi.tellusxdp.com/av2ori/{}/{}/{}/{}/{}.png".format(rspid, sceneid, zoom, xtile, ytile) 
    payload = { "preset": preset, "opacity": opacity, "r": r, "g": g, "b": b, "rdepth": rdepth, "gdepth": gdepth, "bdepth": bdepth }

    r = requests.get( url, params=payload, headers=headers )
    
    return io.imread(BytesIO(r.content))

sceneid1 = "ALAV2A266502880" #取得したいシーンID
tile1 = {"x": 7276, "y": 3226, "z": 13} #取得したい地点のタイル座標

scene1 = get_AVNIR_SCENEID(sceneid1)
img1 = get_AVNIR_image(rspid=scene1["rspId"], sceneid=sceneid1, zoom=tile1["z"], xtile=tile1["x"], ytile=tile1["y"])
io.imshow(img1)

トークンはマイページのＡＰＩアクセス設定（要ログイン）から取得してください。

AVNIR-2のシーンID検索APIとシーン取得APIについて詳しくはリファレンスを確認してください。

OS上で作成した図形（この例では「polygon-1」）を取得するには以下のコードを実行します。

APIで図形を取得する方法について詳しくは【ゼロからのTellusの使い方】Jupyter Notebookで空間データ(GeoJSON)を保存/取得するをご覧ください。

def get_shape_geojson_by_name(name):
    URL = "https://api.tellusxdp.com/v1/shapes"
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(URL , headers=headers)
    shapes = r.json()
    shape = next((shape for shape in shapes["items"] if shape["name"]==name) ,None)
    return json.loads(shape["data"])


shapename1 = "polygon-1"
shape1 = get_shape_geojson_by_name(shapename1)
print(shape1)

3.　クリッピング

では取得した図形で衛星画像をクリッピングしてみましょう。

図形の各点は経度緯度で表されていますが、これを画像のピクセル座標に直してやる必要があります。
そしてskimageという画像ライブラリのpolygonを使い、領域内は1、それ以外は0の配列を作ります。

def get_tile_bbox(xtile, ytile, zoom):
    """
    タイル座標からバウンダリボックスを取得
    https://tools.ietf.org/html/rfc7946#section-5
    """
    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])


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

points1 = shape1["features"][0]["geometry"]["coordinates"][0]
bbox1 = get_tile_bbox(tile1["x"], tile1["y"], tile1["z"])
pxsize1 = {"width": img1.shape[1], "height": img1.shape[0]}

mask1 = get_polygon_image(points1, bbox1,pxsize1)

できた配列を画像として表示すると次のようになります。

io.imshow(mask1)

警告が出る場合がありますが問題ありません。

そしてこの図形画像とnumpyのchooseメソッドを使って衛星画像をクリッピングします。

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)

cliped1 = clip(img1, mask1)
io.imshow(cliped1)

埋立地の形状に衛星画像が切り抜かれました

4.　2地点比較

衛星画像を抜き出すことができるようになったので、これを解析に利用してみましょう。

例えば、田んぼのとある2地点でどちらが稲の発育がいいか、NDVIの値で比べてみます。

※画像内の植物の割合をNDVIから算出する方法については【ゼロからのTellusの使い方】Jupyter NotebookでAVNIR-2/ALOSの光学画像から緑地割合算出をご覧ください。

以下のように自分で比較したい地点を選び、2地点の図形を作成してください。

sceneid2 = "ALAV2A134782770" #取得したいシーンID
tile2 = {"x": 14590, "y": 6158, "z": 14} #取得したい地点のタイル座標

scene2 = get_AVNIR_SCENEID(sceneid2)
img2 = get_AVNIR_image(rspid=scene2["rspId"], sceneid=sceneid2, zoom=tile2["z"], xtile=tile2["x"], ytile=tile2["y"], preset="ndvi")
io.imshow(img2)

bbox2 = get_tile_bbox(tile2["x"], tile2["y"], tile2["z"])
pxsize2 = {"width": img2.shape[1], "height": img2.shape[0]}

shapename2_1 = "polygon-2" #自分で作った図形の名前に変更する
shape2_1 = get_shape_geojson_by_name(shapename2_1)
points2_1 = shape2_1["features"][0]["geometry"]["coordinates"][0]
mask2_1 = get_polygon_image(points2_1, bbox2, pxsize2)

shapename2_2 = "polygon-3" #自分で作った図形の名前に変更する
shape2_2 = get_shape_geojson_by_name(shapename2_2)
points2_2 = shape2_2["features"][0]["geometry"]["coordinates"][0]
mask2_2 = get_polygon_image(points2_2, bbox2, pxsize2)


cliped2_1 = clip(img2, mask2_1)
cliped2_2 = clip(img2, mask2_2)
io.imshow(np.hstack((cliped2_1, cliped2_2)))

import cv2

def filter_green(img):
    green_min = np.array([55,0,0]) 
    green_max = np.array([75,255,255])
    # 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, green_min, green_max)
    masked_img = cv2.bitwise_and(cv_bgr, cv_bgr, mask=mask)

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

filtered_img2_1 = filter_green(cliped2_1)
filtered_img2_2 = filter_green(cliped2_2)
io.imshow(np.hstack((filtered_img2_1, filtered_img2_2)))

def count_not_black(img):
    return len(np.nonzero(img.ravel())[0])

# 切り出した領域内の緑色のピクセルの割合
print("領域1: {:.2f}%".format(count_not_black(filtered_img2_1) / count_not_black(cliped2_1) * 100))
print("領域2: {:.2f}%".format(count_not_black(filtered_img2_2) / count_not_black(cliped2_2) * 100))

この解析例ではB地点と比べてA地点のほうが稲の発育が良い（A地点約71%、B地点約68%）と言えます。

5.　時系列比較

またある地点の生育状況を過去と比較することもできます。
以下のように自分で比較したい地点を選び、2地点の図形を作成してください。

tile3 = {"x": 14578, "y": 6154, "z": 14} #取得したい地点のタイル座標
sceneid3_1 = "ALAV2A137262780" #取得したいシーンID
scene3_1 = get_AVNIR_SCENEID(sceneid3_1)
sceneid3_2 = "ALAV2A083582770" #取得したいシーンID
scene3_2 = get_AVNIR_SCENEID(sceneid3_2)

img3_1 = get_AVNIR_image(rspid=scene3_1["rspId"], sceneid=sceneid3_1, zoom=tile3["z"], xtile=tile3["x"], ytile=tile3["y"], preset="ndvi")
img3_2 = get_AVNIR_image(rspid=scene3_2["rspId"], sceneid=sceneid3_2, zoom=tile3["z"], xtile=tile3["x"], ytile=tile3["y"], preset="ndvi")
io.imshow(np.hstack((img3_1, img3_2)))

bbox3 = get_tile_bbox(tile3["x"], tile3["y"], tile3["z"])
pxsize3 = {"width": img3_1.shape[1], "height": img3_1.shape[0]}

shapename3 = "polygon-4" #自分で作った図形の名前に変更する
shape3 = get_shape_geojson_by_name(shapename3)
points3 = shape3["features"][0]["geometry"]["coordinates"][0]
mask3 = get_polygon_image(points3, bbox3, pxsize3)

cliped3_1 = clip(img3_1, mask3)
cliped3_2 = clip(img3_2, mask3)
io.imshow(np.hstack((cliped3_1, cliped3_2)))

filtered_img3_1 = filter_green(cliped3_1)
filtered_img3_2 = filter_green(cliped3_2)
io.imshow(np.hstack((filtered_img3_1, filtered_img3_2)))

# 切り出した領域内の緑色のピクセルの割合
print("領域1: {:.2f}%".format(count_not_black(filtered_img3_1) / count_not_black(cliped3_1) * 100))
print("領域2: {:.2f}%".format(count_not_black(filtered_img3_2) / count_not_black(cliped3_2) * 100))

この解析例では2007年8月と比べて2008年8月のほうが稲の発育が良くない（A地点約4%、B地点約12%）、つまり2008年は2007年よりも収穫高が少ないと予想できます。

以上が、TellusのJupyter Notebook上で取得した衛星データから、GeoJSONで定義した領域を抜き出す方法でした。

形状のGeoJSONさえ用意できれば、さらに複雑な形状のクリッピングも可能です。また、コードを少し変えればGeoJSONでなく、Shapefileなど別フォーマットのベクトルデータでも同じことができます。

今回は解析例を2つご紹介しました。
実際の解析では、シーンや地点の選択（雲の有無、観測方向、観測時刻）、閾値の調整（緑の数値の範囲）など、解析ごとに試行錯誤が必要となりますが、これを参考に衛星データと空間データを組み合わせたデータ解析に挑戦してみてください。

Tellus参考記事