Palsar-2のSAR画像を取得する

Tag

はじめに

PALSAR-2はALOS-2に搭載された、マイクロ波を地表面に照射し、地表面から反射される電波を受診することで情報を取得する、合成開口レーダ(SAR)と呼ばれるセンサです。得られるSARデータの解像度3~100mとなります。APIを利用して、Palsar-2のシーンと地図タイル画像を取得してみましょう。ASNAROで使用したAPIのように、地図タイル画像を取得するためには、メタデータであるscene_idを知る必要があります。

APIの使い方は登録商品のAPIリファレンスにまとめられています。

※画像を利用する際は、利用ポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/api/v1/palsar2/{min_lat}/{min_lon}/{max_lat}/{max_lon}.png

シーンを取得する

シーン検索APIの変数は以下のようになっています。

Name Input Description
min_lat 必須 最小緯度(-90~90)
min_lon 必須 最小経度(-180~180)
max_lat 必須 最大緯度(-90~90)
max_lon 必須 最大経度(-180~180)

シーンデータをを検索する関数を定義しましょう。

import os, json, requests, math, glob
import numpy as np
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "トークンを貼る"
def get_palser_scene(min_lat, min_lon, max_lat, max_lon):
    """PALSAR-2のメタデータを取得します"""
    url = "https://gisapi.tellusxdp.com/api/v1/palsar2/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_palser_scene(35.524375, 139.785347, 35.724375, 140.005347)
print(len(scenes))
# 霞ヶ浦周辺
img_thumbs = io.imread(scenes[-1]['thumbs_url'])
io.imshow(img_thumbs)

devPalsar_api1_20200220_2.png

Credit:Original data provided by JAXA

霞ヶ浦周辺のSAR画像が取得できました。続いてタイル画像を取得するために、シーンデータの中心座標からタイル座標を計算する関数を定義します。

def get_tile_num(lat_deg, lon_deg, zoom):
    """メタデータの中心座標からタイルのXYを定義"""
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

関数を実行してタイル座標を取得します。

ext_scene = scenes[-1]
zoom = 13
(xtile, ytile) = get_tile_num(ext_scene['clat'], ext_scene['clon'], zoom)
print(xtile, ytile)

(7289,3216)の座標値が取得できました。これで地図タイルを取得する準備ができました。

地図タイル画像を取得する

地図タイル取得APIの変数は以下のようになっています。

https://gisapi.tellusxdp.com//PALSAR-2/{scene_id}/{z}/{x}/{y}.png
Name Input Description
scene_id 必須 シーン検索APIで取得したシーンIDで呼び出すことが可能です。 scene_idは「2237752900」のような数字で示されます。
z 必須 Zoom Level 9-16
x 必須 X座標
y 必須 Y座標

同様に関数を定義します。使用するscene_idはentityIdとなります。

def get_PALSAR_image(scene_id, zoom, xtile, ytile):
    """取得した画像のタイルを定義し、さらに拡大して表示"""
    url = " https://gisapi.tellusxdp.com/PALSAR-2/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

関数を実行して、地図タイル画像を取得し描写します。

ext_tile = get_PALSAR_image(ext_scene['entityId'], zoom, xtile, ytile)
io.imshow(ext_tile)

devPalsar_api1_20200220_3.png

Credit:Original data provided by JAXA

白黒の画像を一枚切り出しただけでは、地図タイル画像がどこを示しているのか判断するのが難しいです。同じ倍率の地図タイル画像を結合して表示して、再度確認しましょう。

def get_palsar_series_image(scene_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """切り出したタイル画像を結合する"""
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            row.append(get_PALSAR_image(scene_id, zoom, topleft_x + x, topleft_y + y))
        img.append(np.hstack(row))
    return  np.vstack(img)

上で定義した関数を実行して、結合画像を描写します。求めたタイル座標値を基準にして、6枚の画像を結合しています。

io.imshow(get_palsar_series_image(ext_scene['entityId'], 13, 7289, 3216, 6, 6))
plt.savefig("SAR_kasumigaura.png")

devPalsar_api1_20200220_1.png

Credit:Original data provided by JAXA

霞ヶ浦と、川、川に沿った農地らしきものが見えてきました。上の画像のようにSARはモノクロの画像であり、光学画像のように対象物が何であるのかが明確ではありません。SARを用いて何ができるのかについては、リモートセンシング基礎(マイクロ波編)に詳しい説明があります。ご参考ください。

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

import os, json, requests, math, glob
import numpy as np
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここにトークンを貼る"
def get_palser_scene(min_lat, min_lon, max_lat, max_lon):
    """PALSAR-2のメタデータを取得します"""
    url = "https://gisapi.tellusxdp.com/api/v1/palsar2/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_palser_scene(35.524375, 139.785347, 35.724375, 140.005347)
print(len(scenes))
# 霞ヶ浦周辺
img_thumbs = io.imread(scenes[-1]['thumbs_url'])
io.imshow(img_thumbs)
def get_tile_num(lat_deg, lon_deg, zoom):
    """メタデータの中心座標からタイルのXYを定義"""
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)
ext_scene = scenes[-1]
zoom = 13
(xtile, ytile) = get_tile_num(ext_scene['clat'], ext_scene['clon'], zoom)
print(xtile, ytile)
def get_PALSAR_image(scene_id, zoom, xtile, ytile):
    """取得した画像のタイルを定義し、さらに拡大して表示"""
    url = " https://gisapi.tellusxdp.com/PALSAR-2/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))
def get_palsar_series_image(scene_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """切り出したタイル画像を結合する"""
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            row.append(get_PALSAR_image(scene_id, zoom, topleft_x + x, topleft_y + y))
        img.append(np.hstack(row))
    return  np.vstack(img)
io.imshow(get_palsar_series_image(ext_scene['entityId'], 13, 7289, 3216, 6, 6))
plt.savefig("SAR_kasumigaura.png")