Met Post

Python、プログラミングのメモ、気象、海洋のこととか

Cartopyで地理データを可視化する1

http://scitools.org.uk/images/cartopy_small.png

図がきれいに書けて細かいカスタマイズも出きるのでデータの可視化はmatplotlibを使うことが多い。地理データ(主に気象、海洋データ)をプロットするときには、Basemapを主に使ってきたのだが、最近CartopyというPythonモジュールを見つけたのでドキュメントを参考にいろいろ触ってみた。

Cartopyとは

Cartopyはイギリス気象局(MetOffice)によってオープンソースライブラリScitoolsの一部として開発されている地図描画用のツールを提供するPythonライブラリ。オブジェクト指向で点や線、画像などを様々な地図座標系に変換したり、matplotlibを使ったプロットへのインタフェースを備え、Shapelyを使ってベクトルデータを扱うこともできるとのこと。

Scitoolsには格子データを扱うモジュールirisも含まれる。特にirisはgribやnetCDFといった気象で使う格子データフォーマットをPythonで扱うモジュールで、ドキュメントを見た感じ便利な機能がたくさんありそうだったのでデータ解析などで使ってみたいなーというのが触ってみた主な理由。

インストール

anacondaで簡単にインストールできる。

conda install -c scitools cartopy

地図を描いてみる

公式レファレンスにランベルト図法とかポーラーステレオ図法とか基本的な投影法が紹介されている。ほとんどmatplotlibと同じ感覚で使えるので扱いやすい。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
import matplotlib.path as mpath
import numpy as np

# 50m解像度用の陸地データ作成
land_50m = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])

正距円筒図法 (PlateCaree)

経緯線は垂直線と等間隔の水平線で表される。経度、緯度をそのままX座標、Y座標に読み変えた図法。赤道付近でのひずみは小さいが、赤道から離れるにつれ東西は拡大されひずみが大きくなる。地図から経緯度を読み取ることが目的であれば最も適した図法。

plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
plt.show()

f:id:ttmych:20151105165642p:plain

plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=140.))
ax.set_extent([110,170,10,60])
ax.add_feature(land_50m)
ax.coastlines(resolution='50m', lw=0.5)
ax.gridlines(xlocs=mticker.MultipleLocator(10), 
             ylocs=mticker.MultipleLocator(10), 
             linestyle='-', 
             color='gray')
plt.show()

f:id:ttmych:20151105171124p:plain

ランベルト正角円錐図法 (LambertConformal)

緯線は同心円、経線は中心からの放射状直線で表される。緯線間隔は正角となるように補正され、頂点(極)は点で表される。割線間隔を0度にすれば、割線は1本で中心緯度と一致する。正角円錐図法は正角円筒図法(メルカトル図法)に比べて面積のひずみを小さく押えることができる。

plt.figure()
ax = plt.axes(projection=ccrs.LambertConformal())
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
plt.show()

f:id:ttmych:20151105171126p:plain

plt.figure()
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=140.0,
                                               central_latitude=30,
                                               standard_parallels=(30,60)))
ax.set_extent([110,170,10,60],ccrs.PlateCarree())
ax.add_feature(land_50m)
ax.coastlines(resolution='50m', lw=0.5)
ax.gridlines(xlocs=mticker.MultipleLocator(10), 
             ylocs=mticker.MultipleLocator(10), 
             linestyle='-', 
             color='gray')
plt.show()

f:id:ttmych:20151105171128p:plain

メルカトル図法(Mercator)

経緯線は垂直線と水平線で表される。赤道付近でのひずみは小さくなり、経線方向のひずみと緯線方向のひずみが等しくなるように、緯線間隔を補正して正角にしている。

plt.figure()
ax = plt.axes(projection=ccrs.Mercator())
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
plt.show()

f:id:ttmych:20151105171130p:plain

plt.figure()
ax = plt.axes(projection=ccrs.Mercator(central_longitude=140.0))
ax.set_extent([110,170,10,60],ccrs.PlateCarree())
ax.add_feature(land_50m)
ax.coastlines(resolution='10m', lw=0.5)
ax.gridlines(xlocs=mticker.MultipleLocator(10), 
             ylocs=mticker.MultipleLocator(10), 
             linestyle='-', 
             color='gray')
plt.show()

f:id:ttmych:20151105171131p:plain

正射方位図法(Orthographic)

地球の無限遠方に光源を置いて平面に投影する図法。中心位置から見て半球を投影できる。 cartopyの場合は、(緯度, 経度) = (central_latitude, central_longitude) で表した点を中心とした半球を投影する。

plt.figure()
ax = plt.axes(projection=ccrs.Orthographic())
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
plt.show()

f:id:ttmych:20151105171132p:plain

plt.figure()
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=0, central_latitude=90))
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
plt.show()

f:id:ttmych:20151105171134p:plain

平射方位図法, ステレオ図法(Stereographic, NorthPolarStereo, SouthPolarStereo)

投影面と反対側の地表に光源を置いて平面に投影する図法。半球を投影できる。中心から離れると距離と面積が拡大されていくが、正角図法なので部分ごとに見れば形状は正しくなる。半球の投影に適している。北極や南極を中心にして使うことがほとんどなので、NorthPolarStereoかSouthPolarStereoを使う。

plt.figure()
ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=-180))
ax.set_extent([-180,180,20,90], ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
plt.show()

f:id:ttmych:20151105171136p:plain

plt.figure()
ax = plt.axes(projection=ccrs.SouthPolarStereo())
ax.set_extent([-180,180,-90,-30], ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)
ax.gridlines(linestyle='-', color='gray')
    
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)

plt.show()

f:id:ttmych:20151105171137p:plain

まとめ

  • 良いところ
    • matplotlibのaxesと同じ感覚で扱えるので使いやすい。
    • Basemapより地図パラメータが少ないので指定しやすい。
    • irisとの連携(たぶんこれが一番のメリット)。
  • 残念なところ
    • Python3に対応していない。Cartopyは対応しているけど、irisが対応していない。
    • (英語のサイトも含めて)サンプルが少ない。