Zenn
🗺️

tif形式の気象データをnumpy配列に変換する

2020/10/23に公開
Python
QGIS
numpy
tech

やりたいこと

  • tif形式で格納されている気象データをnumpy.arrayに変換しpythonでかんたんに利用できるようにしたい！

そもそもtif形式の気象データ とは？

世界全体の気温や降水量などの数値が格納されているデータ．ラスタデータと呼ばれる．
世界地図をグリッド（1マスを1度，0.5度など様々な粒度がある）に区切り，その1マス毎に気象データを格納されている．

例：worldclim

WorldClimの[Historical climate data]
1970-2000年までの1−12月の月平均数値を記録したもの．最低/最高/平均気温，降水量などのデータが10m~30sまでの粒度で利用できる．

とりあえず触ってみる．

今回はとりあえずminimum temperature-5minutesをダウンロードしてみる．

ちなみに

5minutesとは，1辺が緯度経度で言う1/12度の大きさ(1hで1度)の正方形グリッドで世界地図を区切ったもの．
世界地図は緯度が90~-90，経度が180~-180の180360の長方形で表すことができる．今回のデータは1/12度で1マスなので，縦2160横4320マスである．
以下のQiitaの記事に詳しい記載があったので，詳細はそちらを.とてもわかり易い．

ダウンロードして展開すると，以下のような[wc2.1_5m_tmin_XX.tif]と名前のついた12のtifファイルが確認できる．それぞれが1−12月の30年分の月最低気温データになっている．
スクリーンショット 2020-10-23 14.25.50.png

なお今回の粒度だと1ファイルあたりのサイズは11MB程度と可愛いが，これが30s(現在の10倍の粒度)とかになるとファイルサイズも可愛くなくなってくるので注意してほしい．

QGISを使ってデータを取り出す

QGISとはオープンソースGISソフトで，地理情報システムの一種だ．今回のような気象データや地形データなど，様々な地理データを扱うことに向いている．今回はtifファイルをnumpyに変換することが目的なので詳しくは触れない意が，興味がある方は以下を参照していただきたい．

かんたんな使い方

  • QGISをダウンロードし起動，CMD＋Nまたはメニューボタンで新規プロジェクトを開く.

  • 左側にあるレイヤウィンドウに先程のtifファイルをドラッグアンドドロップする．

  • 世界地図が出た！色の濃淡がデータの数値の大小を表している．
    スクリーンショット 2020-10-23 14.34.46.png

  • 右クリックしtifファイルの情報を確認できる．ファイルの幅と高さを確認し，自分の求めているスケールと一致しているか確認する．

ここからは実際にラスタデータをnumpy配列に変換していく．

QGIS内のPythonコンソールを使っていく．

メニューバーにあるPythonアイコンをクリックするとIpythonコンソールが起動される．
これ以降の操作はラスターデータの読み込み · GIS実習オープン教材を参考にした．

numpyに変換する．

はじめに必要なライブラリーをインポートする．


from osgeo import gdal
import numpy as np

次に読み込みたいtifファイルのパスを与え，開く．


uri = 'YOUR_PASS/wc2.1_5m_tmin/wc2.1_5m_tmin_01.tif'
src = gdal.Open(uri, gdal.GA_ReadOnly)
res = src.GetRasterBand(1)
res_arr = res.ReadAsArray()

これでres_arrにnumpy配列が格納されたはずだ．
確認してみよう．


>>> res_arr
array([[-3.4000000e+38, -3.4000000e+38, -3.4000000e+38, ...,
        -3.4000000e+38, -3.4000000e+38, -3.4000000e+38],
       [-3.4000000e+38, -3.4000000e+38, -3.4000000e+38, ...,
        -3.4000000e+38, -3.4000000e+38, -3.4000000e+38],
       [-3.4000000e+38, -3.4000000e+38, -3.4000000e+38, ...,
        -3.4000000e+38, -3.4000000e+38, -3.4000000e+38],
       ...,
       [-1.9660000e+01, -2.1010000e+01, -2.1010000e+01, ...,
        -2.2410000e+01, -2.2410000e+01, -1.9708000e+01],
       [-2.2403999e+01, -2.3939999e+01, -2.3939999e+01, ...,
        -2.5539999e+01, -2.5539999e+01, -2.2451000e+01],
       [-1.6927999e+01, -1.8080000e+01, -1.8080000e+01, ...,
        -1.9279999e+01, -1.9279999e+01, -1.6963999e+01]], dtype=float32)

問題なさそうだ．ちなみに値が-3.4000000e+38となっている箇所は，陸地がなく，データが入っていないNAのマスだ．
念の為arrayのshapeも確認しておく．


>>> res_arr.shape
(2160, 4320)

先程プロパティで確認したサイズと一致している．
np.saveで保存すればnumpyarrayとして利用できる．

おまけ：12ヶ月の平均データを取る．

ここまでくればあとは簡単で，01~12までのtifを読み込んで平均をとってやれば
目的のデータは入手できるはずだ．


dict = {}
from osgeo import gdal
import numpy as np
for a in range(1,13):
    uri = '/YOUR_PASS/wc2.1_5m_tmin/wc2.1_5m_tmin_'+str(a).zfill(2)+'.tif'
    print(uri)
    src = gdal.Open(uri, gdal.GA_ReadOnly)
    res = src.GetRasterBand(1)
    res_arr = red.ReadAsArray()
    dict[a] = res_arr
res = np.stack(dict.values(),axis=1)
data_m = np.mean(res, axis=1)
np.save('/YOUR_PASS/wc2.1_5m_tmin/mean.npy', data_m)

これで完成．

以上です．はじめ検索した際にやり方がまとまったサイトを見つけられなかったのでまとめておきました．

Discussion