LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@wwaaoonn

【JAXA Eearth API】PythonでGeoTIFFから生データを抽出する2つの方法

Last updated at Posted at 2025-09-07

背景

人工衛星の観測データを取得する方法の中でも、特に有力なものが JAXA Earth API である。
大まかな傾向を把握するだけなら、公式ドキュメント中の使用例にあるように，任意の地点の画像や記述統計量の取得で十分なことも多い。

Figure_1.png

一方で、研究や実務においては「特定の1地点ごとの変化量の把握」や「他の指標との併用分析」を行うことが多い。そのためには、可視化や記述統計量として加工される前の生データが必要となる。しかし、公式ドキュメントでは生データの出力方法についてほとんど言及されておらず、抽出や加工の工程が課題となる。

そこで本記事では、この問題に対応するため2種類の方法を検討する。「方法1: raster.imgを直接DataFrame化」の方が簡単にデータを取得できるが、より正確な位置情報を付与したい場合には「方法2: GeoTIFFを解析」も利用してもよい。

取得したいデータの形状

本記事では、以下のようなデータを受け取ることを目的とする。

列名 データ型 中身
y float64 緯度
x float64 経度
value float64 習得した値(植生指数、地表面温度など）

方法1: raster.imgを直接DataFrame化

get_images()→ImageProcessの戻り値に 生の NumPy 配列が格納されおり、公式サンプルにも「生配列を取り出すにはimg.raster.img[0]」と明記されている。こちらはこの NumPy 配列を利用する方法である。

ただし、ip.raster.img[0]の中身は以下のような4次元配列となっている。データそのものは2階層目および3階層目に格納されているため、扱いやすくするためnumpy.sqeeze()などを利用して次元数を減らす。

import numpy as np

# すでに ip = je.ImageProcess(...) まで済んでいる前提

def first_2d(a):
    """配列aを必ず(Y, X)の2次元に落とす。先頭スライス優先。"""
    a = np.asarray(a)
    a = np.squeeze(a)           # サイズ1の軸は潰す
    while a.ndim > 2:           # まだ3D以上なら先頭スライスを取る
        a = a[0]
    if a.ndim != 2:
        raise ValueError(f"Expected 2D, got {a.ndim}D shape={a.shape}")
    return a

arr2 = first_2d(ip.raster.img)

print("変換前の形:", ip.raster.img.shape)
print("変換後の形:", arr2.shape)

# 出力例
# 変換前の形: (1, 168, 383, 1)
# 変換後の形: (168, 383)

また、ip.raster.latlimおよびip.raster.lonlimには、最小(最大)緯度および経度が2次元配列として格納されている。これと先ほど取得した生データが格納されている配列の形状を利用して、各データに対応する緯度経度を作成する。

def lim_pair(lim):
    """latlim/lonlim から (min, max) を安全に取得"""
    b = np.array(lim, dtype=float)
    b = np.squeeze(b)
    if b.ndim == 1 and b.size == 2:
        return float(b[0]), float(b[1])
    b = b.reshape(-1, 2)
    return float(b[0, 0]), float(b[0, 1])

latmin, latmax = lim_pair(ip.raster.latlim)
lonmin, lonmax = lim_pair(ip.raster.lonlim)

nrows, ncols = arr2.shape
dlon = (lonmax - lonmin) / ncols
dlat = (latmax - latmin) / nrows
lons = lonmin + dlon * (0.5 + np.arange(ncols))
lats = latmax - dlat * (0.5 + np.arange(nrows))  # 北→南に並ぶ

print("緯度の配列:", np.repeat(lats, ncols))
print("経度の配列:"np.tile(lons, nrows))

# 出力例
# 緯度の配列: [41.5375 41.5375 41.5375 ... 40.2125 40.2125 40.2125]
# 経度の配列: [139.8625 139.8875 139.9125 ... 141.6375 141.6625 141.6875]

後は、生データと緯度・経度をDataFrameに変換すればよい。

df = pd.DataFrame({
"lat": np.repeat(lats, ncols),
"lon": np.tile(lons, nrows),
"temperature": arr2.ravel()
})
df = df.dropna(subset=["temperature"])

print(df.head())

# 出力例
"""
         lat       lon   value
102  45.5375  141.8875  0.4562
103  45.5375  141.9125  0.4562
104  45.5375  141.9375  0.4562
372  45.5375  148.6375  0.1016
"""

方法2: GeoTIFFを解析

JAXA Earth は COG (Cloud Optimized GeoTIFF) と STAC 形式で配信されている。ImageCollection オブジェクトの stac_* プロパティには STAC の url/json が格納されており、そこから GeoTIFF のアセット URL を取得できる。これを rioxarray で読み込めば、座標情報付きの配列を直接扱えるため、DataFrame 化が容易である。

data.stac_band.url は 2次元配列として URL を保持しているため、対象の文字列を取り出して利用する。

import rioxarray as rxr

# すでに ic = ImageCollection(...).filter_*().select(...) まで済んでいる前提

item_url = data.stac_band.url
href = item_url[0][0]

# GeoTIFFを直接読み込む
da = rxr.open_rasterio(href, masked=True).squeeze(drop=True).rio.write_crs(4326, inplace=False)

このとき、GeoTIFF全体のタイル（例: 東経90°〜180°, 北緯0°〜90°）が読み込まれる。特定の地域だけを切り出したい場合は、クリップ用ポリゴンを別途指定する必要がある。

import geopandas as gpd

# 任意のクリップ用ポリゴンを読み込む
geom = gpd.read_file(geoj_path).geometry[i]

# ポリゴンでクリップ
da_clip = da.rio.clip([geom], drop=True)

GeoTIFFの値はスケーリングが必要な場合がある。例えば NDVI では 0.00001 倍することで実際の値に変換できる。

# DataFrame化とスケーリング
df = da_clip.to_dataframe().reset_index()
df = df.dropna(subset=["NDVI"])
df["NDVI"] = df["NDVI"] * 0.00001

print(df.head())

# 出力例
"""
           y         x  spatial_ref     NDVI
373  45.5375  148.6625            0  0.13172
374  45.5375  148.6875            0  0.13530
375  45.5375  148.7125            0  0.12922
376  45.5375  148.7375            0  0.14475
"""
データ種別ごとのスケーリング係数一覧
データ種別 元の格納値の特徴 スケーリング係数 単位 / 実際の範囲 変換例
NDVI (SGLI, MODIS, VIIRS など) 整数 (int16) ÷10000 -1.0 ～ +1.0 （無次元） 10074 → 10074 / 10000 = 1.0074
EVI (Enhanced Vegetation Index) 整数 (int16) ÷10000 -1.0 ～ +1.0 （無次元） 3456 → 0.3456
反射率 (Surface Reflectance) 整数 (int16) ÷10000 0 ～ 1.0 （無次元、時に >1.0 は補間値） 5678 → 0.5678
LST (Land Surface Temperature, MODIS MYD11 系) 整数 (int16) ×0.02 Kelvin (K) 15000 → 300 K → 26.85 °C
海面水温 (SST, MODIS/VIIRS) 整数 (int16) ×0.01 Kelvin (K) 30000 → 300.00 K
植生量指標 (LAI, FPAR) 整数 (int16) ÷1000 LAI: 0–10, FPAR: 0–1 250 → 0.25
アルベド (Albedo) 整数 (int16) ÷1000 0–1.0 345 → 0.345
降水量 (GSMaP, GPMなど) 整数 (float32もあり) ×0.1 (またはそのまま) mm/h, mm/day など 125 → 12.5 mm/h
標高 (DEM, ALOS AW3D30, SRTM) 整数 (int16/float32) スケーリングなし メートル (m) 250 → 250 m
土壌水分量 (SMAP, AMSR2 など) 整数 (int16) ÷1000 体積含水率 (m³/m³) 245 → 0.245

まとめ

  • JAXA Earth API から取得できるデータは、大きく分けて 2つの扱い方 がある。
    • 方法1: ip.raster.img の NumPy 配列を直接利用する。緯度・経度は自前で生成する。
    • 方法2: GeoTIFF（COG）を rioxarray で解析する。座標情報が埋め込まれているので、そのまま DataFrame 化できる。
  • 試行的・軽量に確認する場合は方法1、座標を活かして本格的な解析を行う場合は方法2を用いるとよい
項目 方法1: raster.img直接 方法2: GeoTIFF解析
データの取得経路 API の戻り値 (ImageProcess) STAC から GeoTIFF を取得
必要な処理 NumPy 配列から次元削減＋緯度経度の生成 rioxarrayで読み込み＋必要ならクリップ
座標情報の扱い 自前で計算して付与する必要あり GeoTIFFに含まれる座標を利用可能
スケーリング データ種類によって要調整 データ種類によって要調整
利点 軽量・依存ライブラリ少ない 座標情報付きで解析が簡単
向いている用途 迅速に小規模な実験を行う場合 正確な座標付き解析やGISとの連携

補足：各方法のサンプルコード

以下のコードは、それぞれの方法で「2018/01/01の各都道府県のNDVI（正規化植生指数）」を取得するサンプルコードである。なお、ポリゴンデータは以下のサイトから取得し、resource/prefectures.geojsonとして保存している。

方法1:raster.imgを直接DataFrame化

from jaxa.earth import je
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

# Set query parameters
dlim = ["2024-01-01T00:00:00","2024-01-10T23:59:59"]
ppu  = 180

# Set information of collection,band
collection = "JAXA.G-Portal_GCOM-C.SGLI_standard.L3-NDVI.daytime.v3_global_half-monthly" 
band = "NDVI"

def first_2d(a):
    """配列aを必ず(Y, X)の2次元に落とす。先頭スライス優先。"""
    import numpy as np
    a = np.asarray(a)
    a = np.squeeze(a)           # サイズ1の軸は潰す
    while a.ndim > 2:           # まだ3D以上なら先頭スライスを取る
        a = a[0]
    if a.ndim != 2:
        raise ValueError(f"Expected 2D, got {a.ndim}D shape={a.shape}")
    return a

def lim_pair(lim):
    """latlim/lonlim から (min, max) を安全に取得"""
    import numpy as np
    b = np.array(lim, dtype=float)
    b = np.squeeze(b)
    if b.ndim == 1 and b.size == 2:
        return float(b[0]), float(b[1])
    b = b.reshape(-1, 2)
    return float(b[0, 0]), float(b[0, 1])

# Get feature collection data
geoj_path = "resource/prefectures.geojson"
geoj = je.FeatureCollection().read(geoj_path).select()

for i in range(len(geoj)):
    # Get an image
    data = je.ImageCollection(collection=collection,ssl_verify=True)\
            .filter_date(dlim=dlim)\
            .filter_resolution(ppu=ppu)\
            .filter_bounds(geoj=geoj[i])\
            .select(band=band)\
            .get_images()
    ip = je.ImageProcess(data)
    print(ip)

    # 画像と緯度経度レンジを“安全に”取り出す
    print(ip.raster.img.shape)
    arr2 = first_2d(ip.raster.img)            # ←ここがポイント（常に2Dに）
    latmin, latmax = lim_pair(ip.raster.latlim)
    lonmin, lonmax = lim_pair(ip.raster.lonlim)

    # 2D前提で座標ベクトルを作成
    nrows, ncols = arr2.shape
    dlon = (lonmax - lonmin) / ncols
    dlat = (latmax - latmin) / nrows
    lons = lonmin + dlon * (0.5 + np.arange(ncols))
    lats = latmax - dlat * (0.5 + np.arange(nrows))  # 北→南に並ぶ

    df = pd.DataFrame({
        "lat": np.repeat(lats, ncols),
        "lon": np.tile(lons, nrows),
        "value": arr2.ravel()
    })
    df = df.dropna(subset=["temperature"])

    print(df.head())

方法2:GeoTIFFを解析

from jaxa.earth import je
import pandas as pd
import numpy as np
import rioxarray as rxr
import geopandas as gpd

# Set query parameters
dlim = ["2018-01-01T00:00:00", "2018-01-10T23:59:59"]
ppu = 180

# Set information of collection, band
collection = "JAXA.G-Portal_GCOM-C.SGLI_standard.L3-NDVI.daytime.v3_global_half-monthly"
band = "NDVI"

# Get feature collection data
geoj_path = "resource/prefectures.geojson"
geoj = je.FeatureCollection().read(geoj_path).select()

geoj = geoj[:2]
print(len(geoj))

for i in range(len(geoj)):
    # Get an image
    data = (
        je.ImageCollection(collection=collection, ssl_verify=True)
        .filter_date(dlim=dlim)
        .filter_resolution(ppu=ppu)
        .filter_bounds(geoj=geoj[i])
        .select(band=band)
        .get_images()
    )

    item_url = data.stac_band.url
    href = item_url[0][0]
    print(f"access to {href}")

    # 直接GeoTIFFを読む
    da = rxr.open_rasterio(href, masked=True).squeeze(drop=True).rio.write_crs(4326, inplace=False)

    # GeoJSONから i 番目の都道府県ポリゴンを取得
    geom = gpd.read_file(geoj_path).geometry[i]

    # ポリゴンでクリップ
    da_clip = da.rio.clip([geom], drop=True)

    # DataFrame化（← clip 後に変換するのが重要）
    df = da_clip.to_dataframe(name=band).reset_index()

    # データの加工
    df = df.dropna(subset=["NDVI"])
    df["NDVI"] = df["NDVI"] * 0.00001  # スケーリング
    print(df.head())
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?