事象発生日：2025-12-23
記事公開日：2025-12-23
アクセス数：19
衛星データのダウンロードって大変ですよね．特に時系列解析をしたいときなどに，数年分の衛星データをダウンロードするのにものすごい時間がかかる，ということもしばしば．
ということで，必要なバンドのみ，かつ，領域のみ切り取ってダウンロードする手法を備忘録的に残します．
これは FOSS4G Advent Calendar 2025 の 12 月 23 日の記事です．
FOSS4Gとは，Free and Open Source Software for Geospatial の略で，地理空間情報に関わるフリー＆オープンソースソフトウェアの総称です．
私が FOSS4G を扱うようになったのは 2 年前（その経緯や最近取り組んでいることは，『人工衛星の作り方ではなく使い方と向き合った２年間』を参照）で，去年はこんな記事を書かせていただきました．
今年は，時系列解析などのために大量の衛星データを効率的にダウンロードする Tips を紹介しようと思います．
衛星データ解析入門については，ぴっかりん氏の Zenn Books 『光学衛星画像解析の１歩目を踏み出してみよう！― QGIS編 ―』を参照するのが良いと思います．
その「Chapter 05 Webサイトで衛星画像を探してみよう、見てみよう」で述べられているように，無料で使える衛星データはたくさんあり，例えば光学衛星 Sentinel-2 やレーダー衛星 (SAR 衛星) Sentinel-1 などは，Copernicus Browser から GUI 経由で入手できます．
GUI でのデータダウンロードは直感的操作で便利ですが，大量のデータの取得には向いていなかったり，また，1 シーン（大きな衛星画像を配信用に区切られた画像データ）ごとにまとめられる zip ファイルには，多くのバンドの画像（Sentinel-1 では VV 偏波と VH 偏波の 2 種類，Sentinel-2 では，B1 ～ B12 の 13 バンド）や校正情報，メタデータなど，様々なデータが含まれているうえに，1 シーンの範囲は下図のようにそれなりに広大なため，1 ダウンロードのファイルサイズも大きいです．
1 シーン全域のすべてのデータを，時系列解析する全期間ダウンロードするのは大変骨の折れる作業です．
さらに，これらのオープン衛星データは北米やヨーロッパのサーバーに保存されていることが多く，ただでさえサイズの大きなデータを海底ケーブル経由の国際通信でダウンロードしなければならず，想像を超える時間がかかってしまうのです．
一方で，データ解析をするときの実際は，ほしいデータ（バンド）や領域はごく一部であったりすることがしばしばです．
ということなので，バンドと領域を絞って効率的にダウンロードしていきましょう．
そのためのスクリプトは以下です．
#!/usr/bin/env python
import logging
import os
import sys
from pathlib import Path
import planetary_computer
import pystac_client
import rasterio
import rasterio.mask
import requests
OUTPUT_DIR = "./output/cog/"
COLLECTIONS = "sentinel-1-rtc"
# https://github.com/stac-api-extensions/query
QUERY = {"sat:orbit_state": {"eq": "descending"}}
SEARCH_AOI = {
"coordinates": [
[
[138.7249, 35.3591],
[138.7249, 35.3687],
[138.7365, 35.3687],
[138.7365, 35.3591],
[138.7249, 35.3591]
]
],
"type": "Polygon"
} # 富士山山頂近辺
START_DATE = "2025-12-01"
END_DATE = "2025-12-20"
# CROP_AOI = None
CROP_AOI = [820023, 3907273, 857351, 3939096] # [xmin, ymin, xmax, ymax] in UTM coordinates
def download_cropped_img(item, pol, crop_aoi):
asset = item.assets[pol]
signed_href = planetary_computer.sign(asset.href) # COG の URL に署名
filename = item.id + f"_{pol}.tif" # 出力ファイル名
file_path = Path(OUTPUT_DIR) / filename
if os.path.exists(file_path):
logging.info(f"Skip {filename} (already exists)")
return
logging.info(f"Download: {filename} ...")
try:
if crop_aoi is None:
# そのまま保存
response = requests.get(signed_href, timeout=6000.0)
response.raise_for_status()
with open(file_path, "wb") as f:
f.write(response.content)
else:
with rasterio.open(signed_href) as src:
# crop_aoi [minx, miny, maxx, maxy] を polygon geometry へ変換
minx, miny, maxx, maxy = crop_aoi
bbox_polygon = {
"type": "Polygon",
"coordinates": [
[
[minx, miny],
[maxx, miny],
[maxx, maxy],
[minx, maxy],
[minx, miny],
]
],
}
# mask を適用
out_image, out_transform = rasterio.mask.mask(src, [bbox_polygon], crop=True)
out_meta = src.meta.copy()
out_meta.update(
{
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform,
}
)
with rasterio.open(file_path, "w", **out_meta) as dest:
dest.write(out_image)
logging.info(f"Saved: {filename}")
except Exception as e:
logging.info(f"Error processing {item.id}: {e}")
def download_stac_metadata(item):
stac_metadata_ural_base = "https://planetarycomputer.microsoft.com/api/stac/v1/collections/" + COLLECTIONS + "/items/"
url = stac_metadata_ural_base + item.id
filename = item.id + ".json"
file_path = Path(OUTPUT_DIR) / filename
if os.path.exists(file_path):
logging.info(f"Skip {filename} (already exists)")
return
logging.info(f"Save meta data: {filename} ...")
try:
response = requests.get(url, timeout=600.0)
except Exception as e:
logging.info(f"Error processing {item.id}: {e}")
return
if response.status_code == 200:
with open(file_path, "wb") as f:
f.write(response.content)
logging.info(f"Saved: {filename}")
else:
logging.info(f"Error processing {item.id}: {response.status_code}")
def main():
logging.basicConfig(
level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s"
)
os.makedirs(OUTPUT_DIR, exist_ok=True)
api_key = os.getenv("PC_SDK_SUBSCRIPTION_KEY")
planetary_computer.settings.set_subscription_key(api_key)
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=COLLECTIONS,
intersects=SEARCH_AOI,
datetime=f"{START_DATE}/{END_DATE}",
query=QUERY,
)
items = search.item_collection()
logging.info(f"Items: {len(items)} items")
for item in items:
download_stac_metadata(item)
for pol in ["vv", "vh"]:
if pol not in item.assets:
continue
download_cropped_img(item, pol, CROP_AOI)
if __name__ == "__main__":
main()
sys.exit(0)
大事な部分を見ていきます．
時空間データ（主に衛星データ）のカタログを表現・検索するためのオープン標準として STAC (SpatioTemporal Asset Catalog) というものがあります．
現在，多くのデータ提供者がこの標準に基づいてデータを公開しており，STAC API を利用することで，API 経由で衛星データを検索し，実データを取得することが可能です．
ここでは，私が良く利用している Microsoft Planetary Computer を例としてサンプルコードを掲載しています．
STAC によるデータ検索などについては，STAC 公式チュートリアルや，Microsoft Planetary Computer の Sentinel-1 データのサンプルコードをご覧いただくとして，サンプルコードでは具体例として Sentinel 1 Radiometrically Terrain Corrected (RTC) データを対象に，撮影日 が 2025/12/01 から 2025/12/20 の期間における富士山山頂付近のデータを検索し，取得しようとしています．
COLLECTIONS = "sentinel-1-rtc"
# https://github.com/stac-api-extensions/query
QUERY = {"sat:orbit_state": {"eq": "descending"}}
SEARCH_AOI = {
"coordinates": [
[
[138.7249, 35.3591],
[138.7249, 35.3687],
[138.7365, 35.3687],
[138.7365, 35.3591],
[138.7249, 35.3591]
]
],
"type": "Polygon"
} # 富士山山頂近辺
START_DATE = "2025-12-01"
END_DATE = "2025-12-20"
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=COLLECTIONS,
intersects=SEARCH_AOI,
datetime=f"{START_DATE}/{END_DATE}",
query=QUERY,
)
items = search.item_collection()
logging.info(f"Items: {len(items)} items")
STAC では，実際のデータファイルは asset と呼ばれる単位で管理されます．
Sentinel などの衛星データでは，各バンド（SAR では偏波）ごとの画像ファイルが個別の asset として提供されるのが一般的です．
Sentinel-1 では偏波は VV と VH の 2 種類が利用可能あり，サンプルコードのように必要な偏波（バンド）のみを選択して画像データを取得することができます．
for item in items:
download_stac_metadata(item)
for pol in ["vv", "vh"]:
if pol not in item.assets:
continue
download_cropped_img(item, pol, CROP_AOI)
def download_cropped_img(item, pol, crop_aoi):
asset = item.assets[pol]
STAC を通じて配信される衛星データは， COG (Cloud Optimized GeoTIFF) と呼ばれる形式で配信されることが一般的です．
COG については空畑の解説記事などを参照いただくとして，Cloud Optimized と名乗っているように，HTTP の Range Request を用いることで，大きな画像ファイルのうちの必要な領域のデータのみに選択的にアクセスできるよう最適化されています．
このように必要な部分のみへ選択的にアクセスできるという特性は， rasterio などのライブラリで画像ファイルを部分的に読み込んで扱えることを意味しており，サンプルコードでは rasterio.mask によって指定した領域のみを切り抜いてダウンロードしています．
ここでの注意点は主に 2 つあります．
1 つめは rasterio.mask.mask の crop パラメタ（デフォルト値は False）で，これを True にしないと選択範囲外がマスクされるだけで，画像のサイズは変わりません（範囲外は NoData で埋められるだけ）．
2 つめは座標系であり，切り抜き処理は asset として取得されるファイルの座標系に基づいて行う必要があります．
Sentinel-1 は UTM 座標系で配信されており， CROP_AOI で指定する切り抜き範囲についても UTM 座標系で定義してます．
CROP_AOI = [820023, 3907273, 857351, 3939096] # [xmin, ymin, xmax, ymax] in UTM coordinates
with rasterio.open(signed_href) as src:
# crop_aoi [minx, miny, maxx, maxy] を polygon geometry へ変換
minx, miny, maxx, maxy = crop_aoi
bbox_polygon = {
"type": "Polygon",
"coordinates": [
[
[minx, miny],
[maxx, miny],
[maxx, maxy],
[minx, maxy],
[minx, miny],
]
],
}
# mask を適用
out_image, out_transform = rasterio.mask.mask(src, [bbox_polygon], crop=True)
out_meta = src.meta.copy()
out_meta.update(
{
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform,
}
)
with rasterio.open(file_path, "w", **out_meta) as dest:
dest.write(out_image)
サンプルコードでは，以下の行でダウンロード領域を切り抜くかどうかを指定できるようになっています．
# CROP_AOI = None CROP_AOI = [820023, 3907273, 857351, 3939096] # [xmin, ymin, xmax, ymax] in UTM coordinates
実際にダウンロードした結果は次のようになり，意図通りの動作になっていることがわかります．
実行時間（ダウンロード時間）について確認するために，サンプルコードの実行ログを掲載します（STAC Search の結果， 3 シーンがヒットしており，それぞれについてメタデータおよび VV，VH のデータをダウンロードしています．）．
通常のように富士山が含まれる 1 シーンをまるごとダウンロードしようとすると，1 画像につき 40 分ほどかかっていますが，富士山周辺のみを切り抜いた場合，約 2 分でダウンロードできています．
なお，2 つめのシーン (S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc) の VV，VH について，前者は接続エラーでダウンロードが途中で失敗し，後者は後述する理由で失敗していますね．
$ python ./download_sentinel-1.py
2025-12-21 05:50:00,973 - INFO - Items: 3 items
2025-12-21 05:50:00,973 - INFO - Save meta data: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc.json ...
2025-12-21 05:50:01,811 - INFO - Saved: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc.json
2025-12-21 05:50:01,811 - INFO - Download: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vv.tif ...
2025-12-21 06:29:54,696 - INFO - Saved: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vv.tif
2025-12-21 06:29:54,928 - INFO - Download: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vh.tif ...
2025-12-21 07:09:06,643 - INFO - Saved: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vh.tif
2025-12-21 07:09:06,875 - INFO - Save meta data: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc.json ...
2025-12-21 07:09:07,781 - INFO - Saved: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc.json
2025-12-21 07:09:07,781 - INFO - Download: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc_vv.tif ...
2025-12-21 07:35:56,339 - INFO - Error processing S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc: ('Connection broken: IncompleteRead([1280720896](tel:1280720896) bytes read, 577416664 more expected)', IncompleteRead([1280720896](tel:1280720896) bytes read, 577416664 more expected))
2025-12-21 07:35:56,340 - INFO - Download: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc_vh.tif ...
2025-12-21 09:26:52,149 - INFO - Error processing S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc: HTTPSConnectionPool(host='[sentinel1euwestrtc.blob.core.windows.net](http://sentinel1euwestrtc.blob.core.windows.net)', port=443): Read timed out.
2025-12-21 09:26:52,150 - INFO - Save meta data: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc.json ...
2025-12-21 09:26:53,025 - INFO - Saved: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc.json
2025-12-21 09:26:53,025 - INFO - Download: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vv.tif ...
2025-12-21 10:04:21,592 - INFO - Saved: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vv.tif
2025-12-21 10:04:21,818 - INFO - Download: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vh.tif ...
2025-12-21 10:41:02,698 - INFO - Saved: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vh.tif
$ python ./download_sentinel-1.py 2025-12-21 22:04:11,423 - INFO - Items: 3 items 2025-12-21 22:04:11,424 - INFO - Save meta data: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc.json ... 2025-12-21 22:04:12,270 - INFO - Saved: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc.json 2025-12-21 22:04:12,270 - INFO - Download: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vv.tif ... 2025-12-21 22:05:37,142 - INFO - Saved: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vv.tif 2025-12-21 22:05:37,143 - INFO - Download: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vh.tif ... 2025-12-21 22:07:21,609 - INFO - Saved: S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc_vh.tif 2025-12-21 22:07:21,610 - INFO - Save meta data: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc.json ... 2025-12-21 22:07:22,720 - INFO - Saved: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc.json 2025-12-21 22:07:22,721 - INFO - Download: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc_vv.tif ... 2025-12-21 22:07:23,440 - INFO - Error processing S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc: Input shapes do not overlap raster. 2025-12-21 22:07:23,441 - INFO - Download: S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc_vh.tif ... 2025-12-21 22:07:23,916 - INFO - Error processing S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc: Input shapes do not overlap raster. 2025-12-21 22:07:23,917 - INFO - Save meta data: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc.json ... 2025-12-21 22:07:24,820 - INFO - Saved: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc.json 2025-12-21 22:07:24,821 - INFO - Download: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vv.tif ... 2025-12-21 22:09:53,307 - INFO - Saved: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vv.tif 2025-12-21 22:09:53,307 - INFO - Download: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vh.tif ... 2025-12-21 22:13:13,348 - INFO - Saved: S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc_vh.tif
後者で Input shapes do not overlap raster. というエラーとともに画像の切り抜きに失敗している理由について，軽く触れておきます．
切り抜きが成功しているシーン (S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc と S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc) と失敗しているシーン (S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc) それぞれについて，シーンの領域は下図の領域でした．
どちらのシーンも SEARCH_AOI で指定した富士山山頂周辺が含まれており意図した動作です．
地理空間情報に馴染みのある勘の鋭い方々はすでにお気づきだと思いますが，富士山周辺は UTM 座標系の Zone 53 と 54 の境目に位置しています（UTM 座標系では，地球表面を平面に投影するときの歪みを低減するために，子午線で区切られた Zone ごとに異なる座標系を定義している）．
つまり，S1A_IW_GRDH_1SDV_20251214T205139_20251214T205204_062316_07CDC4_rtc と S1A_IW_GRDH_1SDV_20251202T205140_20251202T205205_062141_07C6EC_rtc は UTM zone 53N の座標系で配信されている画像ファイルであり，S1A_IW_GRDH_1SDV_20251209T204343_20251209T204408_062243_07CAEF_rtc は UTM Zone 54N の座標系で配信されている画像ファイルです．
このため，CROP_AOI を UTM zone 53N で指定していた本サンプルコードの場合，後者のシーンでは CROP_AOI の領域が意図した領域を示せず，切り抜きに失敗したわけです（座標系めんどくいさい～～～）．
というわけで，衛星データをしゅっとダウンロードする Tips の紹介でした．
とはいえ，例えば機械学習のモデルを作りたいときなど，膨大な領域のデータを取得したい，なんてこともよくあります．
そういうときには，例えば衛星データ配信サーバーに地理的に近いサーバー（eg, AWS リージョンなど）に仮想環境をたててしまって，そこで処理を行う，なんてこともあります．
衛星データをはじめとする地理空間情報の扱いは，データサイズとの戦いであったりもするので，皆さんの Tips などがありましたら，ぜひ教えて下さいね！
名前
Email (※公開されることはありません)
コメント