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

一つ前の記事で、クラウドに対応した衛星データのファイル形式「COG」について紹介しました。

COGとはCloud Optimized Geotiffの略で、名前の通り、「クラウドに最適化された」＋「GeoTIFF」という意味です。「GeoTIFF」は、位置情報のついたTIFF（画像ファイル）という意味で、衛星データで最近よく用いられているファイル形式です。

これまで撮影した画像（シーン）単位でダウンロードしないと読み取れなかった情報を、範囲やピクセルを指定して、欲しいところだけダウンロード（またはネットワーク越しにアクセス）できるようにしたというのがCOGの特徴です。

Tellusでは、2021年10月26日にリリースしたver.3.0からCOG形式でのデータ提供に対応しています。

本記事では、Pythonを使ってTellusのCOG形式の衛星データを取得する方法を具体的に紹介します。

取得したいエリアのポリゴンを用意する

先に述べた通り、COGは衛星データの欲しい部分だけを取得できることが利点です。
そこでまずは「欲しい部分」を指定するためのポリゴンを準備します。

今回は例として、千代田区の形にデータを取得します。
国土数値情報ダウンロードサービスから、千代田区のデータを取得して簡単のために点を間引いています。

#COGを取得するポリゴン
chiyoda = {
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          139.77286583728858,
          35.70370213544527
        ],
        [
          139.7730391708211,
          35.70509546858108
        ],
        [
          139.77013861260696,
          35.70535158650654
        ],
        [
          139.7698949682782,
          35.703366865488476
        ],
        [
          139.76677299607695,
          35.70264602739229
        ],
        [
          139.76640105176716,
          35.7017824138274
        ],
        [
          139.76828411252188,
          35.699700748590715
        ],
        [
          139.74906406075786,
          35.70312438758049
        ],
        [
          139.744478780971,
          35.70230827890401
        ],
        [
          139.7377947731104,
          35.69344108134089
        ],
        [
          139.73117527901593,
          35.68914246855013
        ],
        [
          139.730001362473,
          35.68554191942411
        ],
        [
          139.73120160846747,
          35.68197610841071
        ],
        [
          139.73383379998995,
          35.678878027598955
        ],
        [
          139.73677522699518,
          35.678873027368354
        ],
        [
          139.73906131083095,
          35.67402553225259
        ],
        [
          139.74747936448114,
          35.669656360452734
        ],
        [
          139.74951603123986,
          35.670830802199816
        ],
        [
          139.7543168875144,
          35.66990497353578
        ],
        [
          139.76368416299374,
          35.67251272148496
        ],
        [
          139.77032228384644,
          35.67841050514073
        ],
        [
          139.77090066123662,
          35.684508531843335
        ],
        [
          139.76967569467763,
          35.68691158600092
        ],
        [
          139.77043208837028,
          35.68929797302309
        ],
        [
          139.78258856081345,
          35.694717694166286
        ],
        [
          139.78266080425283,
          35.697855441375566
        ],
        [
          139.77997747098357,
          35.69807797321994
        ],
        [
          139.78039360618004,
          35.700331306444696
        ],
        [
          139.77317053209447,
          35.70187880226496
        ],
        [
          139.77286583728858,
          35.70370213544527
        ]
      ]
    ]
  },
  "bbox": [
    139.730001362473,
    35.669656360452734,
    139.78266080425283,
    35.70535158650654
  ]
}

ライブラリの準備

衛星データを扱う時によく使うライブラリをインストール、インポートしておきます。

!pip install rasterio shapely geopandas scikit-image

import rasterio
from rasterio.mask import mask
from shapely.geometry import shape, Polygon
import geopandas as gpd

Tellusの衛星データのダウンロードURLの取得

Tellusから衛星データをダウンロードしてくるためのURLを取得します。
Tellusでは衛星ごとにデータポリシーが定められており、Tellus外にダウンロードできるものとできないものがあります。

2022年2月現在、Tellus外へのダウンロードが可能なデータはAVNIR-2、PALSAR、ASTER GDEMの3種類です。

ダウンロードURLはAPIで取得する方法もありますが、今回はブラウザ上から取得します。
Tellusのアカウント作成後、Travelerの「衛星データを見る」画面に移動し、ダウンロードしたい衛星データを検索します。

今回は以下のシーンを用いることにしました。
https://www.tellusxdp.com/traveler/myitem/dataset/ea71ef6e-9569-49fc-be16-ba98d876fb73/data/1772c62e-096a-4db4-8926-95a5476913a0

衛星データは、以下のように1シーンの中に複数のファイルが存在しています。

それぞれがダウンロードまたは、URLを表示できるようになっているので、この中から画像データ（79.51MiBとなっている容量の大きいもの）を一つ選んで、URLを表示させます。

1時間有効なダウンロードURLが発行されるので、コピーします。

# ダウンロードURL
url = "ダウンロードURLを貼り付け"

衛星データを必要なエリアだけトリミングして取得する

冒頭で作ったポリゴンデータを元に、衛星データの必要な部分だけをくり抜いて取得します。

# ポリゴンを含む最小限の範囲だけを取得の上、トリミングして返す
def get_masked_image(image_href, poly):
    data = rasterio.open(image_href)
    masked, mask_transform = mask(dataset=data, shapes=poly.geometry.to_crs(crs=data.crs.data), crop=True, all_touched=True)
    return masked[0,:,:], mask_transform

このget_masked_imageの中でファイルの投影座標系に変換するのにgeopandasを使用するため、ポリゴンをgeopandasのオブジェクトに変換しておく必要があります。

mask_polygon = gpd.GeoDataFrame(crs={'init': 'epsg:4326'}, geometry=[shape(chiyoda["geometry"])])

衛星データからポリゴンデータでマスクされたデータを取得します。

masked, _ = get_masked_image(url, mask_polygon)

取得したデータを表示する

実際にちゃんとマスクされているのか、表示させてみます。

from skimage import io
io.imshow(masked)

以下のように千代田区の形に切り抜かれたデータが表示されます。

試しに千代田区だけでなく衛星データのシーン全体を呼び出してみましょう。シーン全体をダウンロードするのには時間がかかるので注意してください。

io.imshow(rasterio.open(url).read(1))

千代田区を含む、かなり広域なエリアが取得されることが分かります。

COGが如何に素早く必要なデータだけを取得できるかが体感いただけたかと思います。

まとめ

本記事ではCOGデータの特性を活かして、必要な範囲だけのデータを取得する方法をご紹介しました。

衛星データは1シーンのデータ量が大きく、全部を取ってくるのはなかなか時間もストレージも必要になります。COG形式を上手く活用して、素早く衛星データ解析を行っていただければと思います。