地図タイル画像を比較する
はじめに
本記事ではTellusのJupyter上で取得した衛星データからGeoJSONで定義した領域を抜き出す（クリッピングする）方法を紹介します。
抜き出す領域を指定する
本記事では、画像の一部を切り出して解析する方法について説明をしていきます。どのような場合に任意の領域を抜き出す必要があるでしょうか。例えば農業なら、以下のようなケースです。
・水田のA地点とB地点で稲の発育を比較し、品種に適した土壌条件を分析する。
・ある畑の生育状況を過去の同じ地点と比較し、その年の収穫高を予想する。
抜き出す領域を指定するには、位置情報のベクトルデータが必要です。今回は最も簡単な例として、Tellus OS上で地図や衛星写真を眺めながら、抜き出したい領域の図形を作成し、それを利用します。最初に抜き出したい衛星データを選択しましょう。例ではAVNIR-2のデータを使って東京の晴海周辺を調べてみます。画面中央に竹芝ふ頭を表示して、左上の「マイライブラリ」をクリックします。
開いたパネルからTellus OS上で扱える衛星データ群が列挙されていることが確認できます。初期設定では全てがオフになっていますので、目のアイコンを押して欲しい画像を取得する状態に設定します。
OS上で対象地点へとズームします。Tellus OSは、ズームに応じて最適な場所と思われるシーンデータを返します。そのため特にデータの選択をする必要なく、対象地点を含む画像が入手できます。どのようなシーンデータが含まれているのかを知りたい場合は、「詳細」をクリックしましょう。
上のように現在のズームや位置から判断してAVNIRのシーンデータが返されています。
次に抜き出したい領域の図形を作成します。今回は晴海周辺の埋立地を抜き出してみます。画面上部のアイコンから台形型のアイコンをクリックしポリゴンを作成します。画面の中心に見えている埋立地を囲ってみましょう。
きちんと図形を閉じたら完成で、できた図形は画面左の「マイライブラリ」、「シェイプ」から確認できます。今回は「polygon-1」という名前で保存されています。
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)
OS上で作成した図形（この例では「polygon-1」）を取得するには以下のコードを実行します。
APIで図形を取得する方法について詳しくは【ゼロからのTellusの使い方】Jupyter Notebookで空間データ(GeoJSON)を保存/取得するをご覧ください。
def get_shape_geojson_by_name(name):
URL = "https://api.tellusxdp.com/v2/geojson-files/"
headers = {
"Authorization": "Bearer " + TOKEN
}
r = requests.get(URL , headers=headers)
shapes = r.json()
shape = next((shape for shape in shapes["results"] if shape["original_name"]==name) ,None)
r = requests.get(shape["url"])
return r.json()
shapename1 = "polygon-1"
shape1 = get_shape_geojson_by_name(shapename1)
print(shape1)
クリッピング
取得した図形で衛星画像をクリッピングしてみましょう。図形の各点は経度緯度で表されていますが、これを画像のピクセル座標に変換する必要があります。そして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)
警告が出る場合がありますが問題ありません。この1と0で構成された図形画像と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)
2地点比較
衛星画像を抜き出すことができるようなりました。これを解析に利用してみましょう。例えば、水田のとある2地点でどちらが稲の発育が良好かをNDVIの値で比べてみます。
※画像内の植物の割合をNDVIから算出する方法については【ゼロからのTellusの使い方】Jupyter NotebookでAVNIR-2/ALOSの光学画像から緑地割合算出をご覧ください。
以下のように比較したい地点を選び、2地点の図形を作成してください。同様の場所を使う場合は、GeoJSONファイルを作成して、40.73440677988237,140.59194660134384の座標をポイントして下さい。青森県にポイントがされたのがわかると思います。開発環境を用いて、同様にポイントを作成することもできます。
インポートされたファイル名はaomori_point.geojsonになります。
対象地へとズームします。右下にある時間スケールを1monthに変えると、該当する場所に対して存在する月別のデータが示されます。今回は8月のデータを選択しましょう(上記と同様に詳細情報を調べます。シーンIDはALAV2A134782770であることがわかります)。最後に比較する2地点の図形を多角形ツールで作成します（右上がA地点、左上がB地点）。
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))
緑の濃い（植物の繁殖度合いが高い）部分を絞り込み、領域内での割合を計算。A地点では69.66%、B地点では、54.00%となります。この解析例ではB地点と比べてA地点のほうが稲の発育が良いと考えられます。
時間での変化を捉える
ある地点の生育状況を過去と比較することもできます。以下のように自分で比較したい地点を選び、2地点の図形を作成してください。場所は（40.8034148344062,140.32504320144656）付近 となります。二箇所の図形を作成した方法を踏襲します。GeoJSONから緯度経度を用いて場所を特定してください。先ほど作成した図形の北西側にポイントが作成されたのが確認できます。
該当する場所を拡大します。下の図で描いたような図形をTellus上で作成してください。
tile3 = {"x": 14578, "y": 6154, "z": 14} #取得したい地点のタイル座標
sceneid3_1 = "ALAV2A083582770" #取得したいシーンID
scene3_1 = get_AVNIR_SCENEID(sceneid3_1)
sceneid3_2 = "ALAV2A137262780" #取得したいシーン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)))
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月のほうが稲の発育が悪かったと考えられます（2007年が約10%に対して、2008年では約3%）。そのため指定した範囲内では、2008年は2007年よりも収穫高も少なかったと予想できます。
以上が、TellusのJupyterLab上で取得した衛星データから、GeoJSONで定義した領域を抜き出す方法でした。
形状のGeoJSONさえ用意できれば、さらに複雑な形状のクリッピングも可能です。また、コードを少し変えればGeoJSONでなく、Shapefileなど別フォーマットのベクトルデータでも同じことができます。
実際の解析では、シーンや地点の選択（雲の有無、観測方向、観測時刻）、閾値の調整（緑の数値の範囲）など、解析ごとに試行錯誤が必要となりますが、これを参考に衛星データと空間データを組み合わせたデータ解析に挑戦してみてください。
今回使用したスクリプトです。
import os, json, requests, math
import numpy as np
from skimage import io
from io import BytesIO
from skimage.draw import polygon
get_ipython().run_line_magic('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)
def get_shape_geojson_by_name(name):
URL = "https://api.tellusxdp.com/v2/geojson-files/"
headers = {
"Authorization": "Bearer " + TOKEN
}
r = requests.get(URL , headers=headers)
shape = next((shape for shape in shapes["results"] if shape["original_name"]==name) ,None)
r = requests.get(shape["url"])
return r.json()
shapename1 = "polygon-1"
shape1 = get_shape_geojson_by_name(shapename1)
print(shape1)
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)
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)
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)
URL = "https://api.tellusxdp.com/v2/geojson-files"
def post(data):
headers = {
"Authorization": "Bearer " + TOKEN,
"Content-Type": "application/json"
}
r = requests.post(URL, headers=headers, json=data)
return r.json()
geojson = {
"type": "Point",
"coordinates": [140.59194660134384, 40.73440677988237]
}
data = {
"originalName": "test",
"transparency": 0.5,
"windowNumber": 1,
"data": geojson
}
print(post(data))
bbox2 = get_tile_bbox(tile2["x"], tile2["y"], tile2["z"])
pxsize2 = {"width": img2.shape[1], "height": img2.shape[0]}
shapename2_1 = "polygon-4" #自分で作った図形の名前に変更する
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))
geojson = {
"type": "Point",
"coordinates": [140.32504320144656, 40.8034148344062]
}
data = {
"originalName": "test",
"transparency": 0.5,
"windowNumber": 1,
"data": geojson
}
print(post(data))
tile3 = {"x": 14578, "y": 6154, "z": 14} #取得したい地点のタイル座標
sceneid3_1 = "ALAV2A083582770" #取得したいシーンID
scene3_1 = get_AVNIR_SCENEID(sceneid3_1)
sceneid3_2 = "ALAV2A137262780" #取得したいシーン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-5" #自分で作った図形の名前に変更する
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))