再解析データ

再解析データとは，比較的新しい予報モデルに対して過去の観測データを同化したものです。 数値天気予報においては，予報モデルに対して観測を同化した解析値が作成され，初期値として用いられています。 予報モデルは改良が重ねられていくため，解析値の品質はモデルの変遷と共に変化していきます。 データ同化手法も大きく進歩し，数十年前とは全く異なる手法が用いられています。 モデルやデータ同化手法の改良の影響を取り除くことにより，時間的により品質が安定したデータセットになります。 再解析とは，現業予報で一度行った解析を再びやり直すということを意味しています。

Caution

衛星データは世代間のセンサの品質の差やセンサの経年劣化が考慮されていますが，解析値は入力される観測データの種類や量にも影響を受けます。 品質を一貫させるため，大きな精度向上をもたらす衛星データを使わず従来型観測のみを同化したデータセット（JRA-55C[1]，ERA-20C [2]，20CRv3[3]） が作られています。

本書ではサイズが小さく入手が容易なNCEP/NCAR再解析[4]を用います。 NCEP/NCAR再解析は，京都大学生存圏研究所グローバル大気観測データや米国国立大気研究センター（NCAR）研究データアーカイブ等から入手可能です。

Note

より新しい再解析データセットとして3つあげておきます。

• JRA-55[5]: 気象庁が作成した長期再解析。より長期のデータセット作成中でJRA-3Qが近日完了予定。

• MERRA-2[6]: NASAによる衛星データを駆使した再解析

• ERA-5[7]: 精度に定評があるECMWFが作成したデータセット。積雪深が過大[8]

!curl -O http://database.rish.kyoto-u.ac.jp/arch/ncep/data/ncep.reanalysis.derived/surface/slp.mon.mean.nc

NCEP/NCAR再解析データはNetCDF形式で提供されています。 NetCDFはデータと共に座標などのメタデータが一つのファイルの中に保存されています。 xarrayはnetcdf4ライブラリを使ってNetCDFファイルを読み書きすることができます。 データを読むにはxarray.open_dataset()を用います。

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr


ds = xr.open_dataset("slp.mon.mean.nc")

dsを印字するとメタ情報が表示されます。

ds

<xarray.Dataset>
Dimensions:  (lat: 73, lon: 144, time: 904)
Coordinates:
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2023-04-01
Data variables:
    slp      (time, lat, lon) float32 ...
Attributes:
    description:    Data is from NMC initialized reanalysis\n(4x/day).  These...
    platform:       Model
    Conventions:    COARDS
    NCO:            20121012
    history:        Thu May  4 18:12:35 2000: ncrcat -d time,0,622 /Datasets/...
    title:          monthly mean slp from the NCEP Reanalysis
    dataset_title:  NCEP-NCAR Reanalysis 1
    References:     http://www.psl.noaa.gov/data/gridded/data.ncep.reanalysis...

xarray.Dataset

• Dimensions:
  • lat: 73
  • lon: 144
  • time: 904
• Coordinates: (3)
  • lat
    (lat)
    float32
    90.0 87.5 85.0 ... -87.5 -90.0

    units :

    degrees_north

    actual_range :

    [ 90. -90.]

    long_name :

    Latitude

    standard_name :

    latitude

    axis :

    Y

    array([ 90. ,  87.5,  85. ,  82.5,  80. ,  77.5,  75. ,  72.5,  70. ,  67.5,
            65. ,  62.5,  60. ,  57.5,  55. ,  52.5,  50. ,  47.5,  45. ,  42.5,
            40. ,  37.5,  35. ,  32.5,  30. ,  27.5,  25. ,  22.5,  20. ,  17.5,
            15. ,  12.5,  10. ,   7.5,   5. ,   2.5,   0. ,  -2.5,  -5. ,  -7.5,
           -10. , -12.5, -15. , -17.5, -20. , -22.5, -25. , -27.5, -30. , -32.5,
           -35. , -37.5, -40. , -42.5, -45. , -47.5, -50. , -52.5, -55. , -57.5,
           -60. , -62.5, -65. , -67.5, -70. , -72.5, -75. , -77.5, -80. , -82.5,
           -85. , -87.5, -90. ], dtype=float32)

  • lon
    (lon)
    float32
    0.0 2.5 5.0 ... 352.5 355.0 357.5

    units :

    degrees_east

    long_name :

    Longitude

    actual_range :

    [ 0. 357.5]

    standard_name :

    longitude

    axis :

    X

    array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
            25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
            50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
            75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
           100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
           125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
           150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
           175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
           200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
           275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
           300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
           325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
           350. , 352.5, 355. , 357.5], dtype=float32)

  • time
    (time)
    datetime64[ns]
    1948-01-01 ... 2023-04-01

    long_name :

    Time

    delta_t :

    0000-01-00 00:00:00

    prev_avg_period :

    0000-00-01 00:00:00

    standard_name :

    time

    axis :

    T

    actual_range :

    [1297320. 1956936.]

    array(['1948-01-01T00:00:00.000000000', '1948-02-01T00:00:00.000000000',
           '1948-03-01T00:00:00.000000000', ..., '2023-02-01T00:00:00.000000000',
           '2023-03-01T00:00:00.000000000', '2023-04-01T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • slp
    (time, lat, lon)
    float32
    ...

    long_name :

    Sea Level Pressure

    valid_range :

    [ 870. 1150.]

    units :

    millibars

    precision :

    1

    var_desc :

    Sea Level Pressure

    level_desc :

    Sea Level

    statistic :

    Mean

    parent_stat :

    Other

    dataset :

    NCEP Reanalysis Derived Products

    actual_range :

    [ 955.56085 1082.5582 ]

    [9502848 values with dtype=float32]

• Attributes: (8)

  description :

  Data is from NMC initialized reanalysis (4x/day). These are the 0.9950 sigma level values.

  platform :

  Model

  Conventions :

  COARDS

  NCO :

  20121012

  history :

  Thu May 4 18:12:35 2000: ncrcat -d time,0,622 /Datasets/ncep.reanalysis.derived/surface/slp.mon.mean.nc ./surface/slp.mon.mean.nc Mon Jul 5 23:22:35 1999: ncrcat slp.mon.mean.nc /Datasets/ncep.reanalysis.derived/surface/slp.mon.mean.nc /dm/dmwork/nmc.rean.ingest/combinedMMs/slp.mon.mean.nc /home/hoop/crdc/cpreanjuke2farm/cpreanjuke2farm Thu Oct 26 23:42:16 1995 from pre.sig995.85.nc created 95/02/06 by Hoop (netCDF2.3) Converted to chunked, deflated non-packed NetCDF4 2014/09

  title :

  monthly mean slp from the NCEP Reanalysis

  dataset_title :

  NCEP-NCAR Reanalysis 1

  References :

  http://www.psl.noaa.gov/data/gridded/data.ncep.reanalysis.derived.html

.sel()メソッドを用いると配列の添え字ではなく座標値を用いた選択ができます。 東経135度，北緯35度における8月の海面気圧の年々変動を調べてみましょう。 添え字は0から始まるので[7::12]となっていることに注意してください。

ds.slp.sel(lon=135, lat=35)[7::12].plot(figsize=[7,5])
plt.show()

Cartopy

海面気圧を正射図法で描いてみます。 正射図法を用いると宇宙から見た地球や地球儀のように描かれます。 地図投影にはCartopyを用います。 CartopyはPython以外で書かれた外部ライブラリに依存しています。 自分のパソコンを使っている場合は，インストール方法や付録を参考にしてください。 Google Colaboratoryの場合は!を使ってシェルを呼び出し，apt-getとpipを使って依存ライブラリとCartopyをインストールします。

!apt-get install -qq libgdal-dev libproj-dev
!pip install cartopy

Caution

Google Colaboratoryで以下のようなエラーが出る場合はshapelyのバージョンがcartopyと整合していません。一旦shapelyをアンインストールしてソースからコンパイルします。

Geometry must be a Point or LineString

!pip uninstall --yes shapely
!pip install shapely --no-binary shapely
!pip install cartopy

• 1行目: 地図投影が定義されているcartopy.crsをccrsという名前でインポートします。

• 5行目: 全球データを描画するとグリニッジ子午線や日付変更線で切れてしまうことがあります。cartopy.util.add_cyclic_pointを利用して経度方向が周期的になるように点を追加します。

• 8行目: 投影をprojectionに指定します。 ccrs.Orthographic()は正射図法で衛星から地球を眺めたような投影法です。

• 10行目: 周期的になるように点を追加した経度，緯度，データを与えて陰影付等値線をcontourf()で描きます。 projectionは用いる地図投影を指定するのに対し，描画関数に与えるtransformは，データが定義されている座標系を指定します。 地図投影は正射図法ですが，座標を経度緯度で与えているのでtransformはPlateCarree()とします。 ccrs.PlateCarree()は標準緯度を赤道とした正距円筒図法で，経度と緯度をそのまま使います。

• 11行目: 海岸線を描きます。

• 12行目: 描画範囲を全球とします。

• 13行目: カラーバーを描きます。

Caution

Cartopyは初めて必要となったときに海岸線等のデータをダウンロードします。サーバがダウンしているとデータが取得できないことがあります。 coastlines()はNatural Earthのshapefileを簡単に描くためのメソッドです。 CartopyはNatural Earth以外のデータを取得して描くことができます。 coastlines()がうまくいかないときは，その代わりにGSHHSのデータを使ってみましょう。 地図に描く部品をCartopyではfeatureと呼んでいます。 GSHHSFeature()をadd_feature()メソッドで追加することにより海岸線を描きます。

import cartopy.feature as cfeature

ax.add_feature(cfeature.GSHHSFeature())

import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point


wdata, wlon = add_cyclic_point(ds.slp.sel(time="2020-12")[0], ds.lon)

fig = plt.figure(figsize=[9, 7])
ax = fig.add_subplot(111, projection=ccrs.Orthographic(135,35))

p = ax.contourf(wlon, ds.lat, wdata, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_global()
fig.colorbar(p)
plt.show()

細かなパラメタを指定せずに簡単に使えように，NaturalEarthの陸LAND，海OCEAN，湖水LAKES等が定義されています。

• 5行目: ccrs.LambertConformal()はランベルト図法を投影に指定します。中心となる経度と緯度を与えます。

• 7～9行目: NaturalEarthの海，湖水，陸のデータを利用します。初回はデータをダウンロードするので，メッセージが表示されます。

• 12行目: 等値線を指定するために4 hPa間隔で960から1036 hPaまでの配列を作りlevelsに与えます。等値線の色は黒にします。

• 14行目: ax.clabel()メソッドを呼び，等値線ラベルを自動生成（manual=False）し，等値線上に配置（inline=True）します。

import cartopy.feature as cfeature


fig = plt.figure(figsize=[7,7])
ax = fig.add_subplot(111, projection=ccrs.LambertConformal(135,35))

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.LAND)
ax.coastlines()
p = ax.contour(wlon, ds.lat, wdata,
               levels=np.arange(960, 1040, 4), colors=["k"],
               transform=ccrs.PlateCarree())
ax.clabel(p, manual=False, inline=True)
ax.set_extent([90, 180, 0, 70], ccrs.PlateCarree())
ax.gridlines()
plt.show()

大陸上のシベリア高気圧とアリューシャン低気圧がコントラストが明瞭です。