読者です 読者をやめる 読者になる 読者になる

Met Post

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

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

cartopy matplotlib python

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

今回もNatural Earthで公開されているshapeファイルデータを使って地図にいろいろ描画してみる。今回はshapeファイルからデータを取り出して加工して都道府県境界を描画してみる。

shapereader.natural_earthはカテゴリを指定することで、Natural Earthのshapeファイルのパス(URI)を返してくれる。以前にダウンロードしたことがあれば、キャッシュを使ってくれるので時間短縮できる。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader

shpfilename = shapereader.natural_earth(resolution='10m',
                                        category='cultural',
                                        name='admin_1_states_provinces')

reader = shapereader.Reader(shpfilename)
provinces = reader.records()

10m解像度のstates-provincesデータには各国の国や行政区分地域のデータが含まれている。ここから日本のデータのみを取り出してみる。各レコードのキーを調べてみる。

province = next(provinces)
print sorted(province.attributes.keys())
# キーadminを表示
print province.attributes['admin']
['OBJECTID_1', 'abbrev', 'adm0_a3', 'adm0_label', 'adm0_sr', 'adm1_cod_1', 'adm1_code', 'admin', 'area_sqkm', 'check_me', 'code_hasc', 'code_local', 'datarank', 'diss_me', 'featurecla', 'fips', 'fips_alt', 'gadm_level', 'geonunit', 'gn_a1_code', 'gn_id', 'gn_level', 'gn_name', 'gn_region', 'gns_adm1', 'gns_id', 'gns_lang', 'gns_level', 'gns_name', 'gns_region', 'gu_a3', 'hasc_maybe', 'iso_3166_2', 'iso_a2', 'labelrank', 'latitude', 'longitude', 'mapcolor13', 'mapcolor9', 'name', 'name_alt', 'name_len', 'name_local', 'note', 'postal', 'provnum_ne', 'region', 'region_cod', 'region_sub', 'sameascity', 'scalerank', 'sov_a3', 'sub_code', 'type', 'type_en', 'wikipedia', 'woe_id', 'woe_label', 'woe_name']

Aruba

調べてみると'admin'が国名になっているようだ。Reader.records()はシェープファイルに含まれるデータのRecordを取り出すイテレータを返す。Recordにはshapely.geometryのインスタンスとそのメタデータが含まれる。メタデータの情報から日本の都道府県データのRecordを取り出してみる。

import itertools

reader = shapereader.Reader(shpfilename)
provinces = reader.records()

# 条件に合致するRecordへのイテレータ
provinces_of_japan = itertools.ifilter(lambda province: province.attributes['admin'] == 'Japan', provinces)

# 都道府県名を確認
count = 0
for province in provinces_of_japan:
    print province.attributes['name']
    if count >= 5: break
    count = count + 1
Hiroshima
Okayama
Shimane
Tottori
Yamaguchi

上手く取り出せているようだ。取り出した都道府県のレコードからgeometryを取り出して、GeoAxes.add_geometries()を使って追加する。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shapereader
import itertools

# acquire the provinces dataset from Natural Earth
shpfilename = shapereader.natural_earth(resolution='10m',
                                        category='cultural',
                                        name='admin_1_states_provinces')

# get shapefile records
reader = shapereader.Reader(shpfilename)
provinces = reader.records()

# filter provinces of Japan
provinces_of_japan = itertools.ifilter(lambda province: province.attributes['admin'] == 'Japan', provinces)

# plot
colors = itertools.cycle(['red', 'blue', 'green', 'lime', 'orange', 'cyan', 'purple'])

plt.figure(figsize=[8,8])
ax = plt.axes(projection=ccrs.PlateCarree())

for province in provinces_of_japan:
    geometry = province.geometry
    ax.add_geometries(geometry, ccrs.PlateCarree(), facecolor=next(colors))
ax.set_extent([135, 142, 33, 40])
plt.show()

f:id:ttmych:20160604160706p:plain

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shapereader
import itertools

# acquire the provinces dataset from Natural Earth
shpfilename = shapereader.natural_earth(resolution='10m',
                                        category='cultural',
                                        name='admin_1_states_provinces_lines')

# get shapefile records
reader = shapereader.Reader(shpfilename)
provinces = reader.records()

# filter provinces of Japan
provinces_of_japan = itertools.ifilter(lambda province: province.attributes['adm0_name'] == 'Japan', provinces)

# plot
plt.figure(figsize=[8,8])
ax = plt.axes(projection=ccrs.PlateCarree())

for province in provinces_of_japan:
    geometry = province.geometry
    ax.add_geometries(geometry, ccrs.PlateCarree(), facecolor='none', linestyle=':')
ax.coastlines(resolution='10m')
ax.set_extent([138, 141, 34, 37])
plt.show()

f:id:ttmych:20160604160711p:plain

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

cartopy matplotlib python

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

cartopyにはNatural Earthで公開されているshapeファイルデータを使うためのモジュールが実装されている。Natural Earthのデータを使って地図を描いてみた。

海岸線を描画する

GeoAxes.coastlines()を使えば、作成したGeoAxesクラスの地図にNaturalEarthの海岸線データを描画することができる。NaturalEarthの海岸線ベクトルデータを使用している。解像度は3種類ある。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

plt.figure(figsize=[15,5])

## 110m resolution (default)
ax1 = plt.subplot(1, 3, 1, projection=ccrs.PlateCarree())
ax1.coastlines(resolution='110m')
ax1.set_extent([120,150,20,50], ccrs.PlateCarree())
ax1.set_title('110m coastline')

## 50m resolution
ax1 = plt.subplot(1, 3, 2, projection=ccrs.PlateCarree())
ax1.coastlines(resolution='50m')
ax1.set_extent([120,150,20,50], ccrs.PlateCarree())
ax1.set_title('50m coastline')

## 10m resolution
ax1 = plt.subplot(1, 3, 3, projection=ccrs.PlateCarree())
ax1.coastlines(resolution='10m')
ax1.set_extent([120,150,20,50], ccrs.PlateCarree())
ax1.set_title('10m coastline')

f:id:ttmych:20160604155435p:plain

様々なベクトルデータを描画する

keyword description
edgecolor/ec 境界線の色
facecolor/fc 塗りつぶしの色

cartopy.feature.NaturalEarthFeatureを使うと、NaturalEarthのデータをダウンロードしてポリゴンを作ることができる。作成したポリゴンをGeoAxesに追加すればいい。

使い方は、

NaturalEarthFeature(category, name, scale, **kwargs)

で作成する。引数のcategory, name, scaleは、NaturalEarthのカテゴリ、データセット名、解像度に対応している。詳細はNaturalEarthのダウンロードのページを参照。

例えば、

NaturalEarthFeature('physical', 'land', '50m')

の場合は、「download/50m/cultural/ne_50m_land.zip」がダウンロードされて描画される。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

plt.figure(figsize=(10, 10))

# physical category
land_50m  = cfeature.NaturalEarthFeature('physical', 'land', '50m', 
                                         edgecolor='face',                   # same color with facecolor
                                         facecolor=cfeature.COLORS['land'])  # use predefiend color of cartopy
ocean_50m = cfeature.NaturalEarthFeature('physical', 'ocean', '50m', 
                                         edgecolor='face',                   # same color with facecolor
                                         facecolor=cfeature.COLORS['water']) # use predefiend color of cartopy
lakes_50m = cfeature.NaturalEarthFeature('physical', 'lakes', '50m', 
                                         edgecolor='face',                   # same color with facecolor
                                         facecolor=cfeature.COLORS['water']) # use predefiend color of cartopy
river_50m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '50m',
                                         edgecolor=cfeature.COLORS['water'], # use predefiend color of cartopy
                                         facecolor='none')                   # no filled color

# cultural category
countries_50m  = cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', '50m', 
                                              edgecolor='gray',
                                              facecolor='none')  # no filled color

states_50m  = cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '50m', 
                                           edgecolor='gray',
                                           facecolor='none')  # no filled color

states_10m  = cfeature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', '10m', 
                                           edgecolor='gray',
                                           facecolor='none')  # no filled color

ax1 = plt.subplot(2, 2, 1, projection=ccrs.PlateCarree())
ax1.set_extent([-20, 60,-40, 40], ccrs.PlateCarree())
ax1.add_feature(land_50m)
ax1.add_feature(ocean_50m)
ax1.add_feature(lakes_50m)
ax1.add_feature(river_50m)
ax1.set_title('physical')

ax2 = plt.subplot(2, 2, 2, projection=ccrs.PlateCarree())
ax2.set_extent([-20, 60,-40, 40], ccrs.PlateCarree())
ax2.coastlines(resolution='50m')
ax2.add_feature(countries_50m)
ax2.set_title('countries_border')

ax3 = plt.subplot(2, 2, 3, projection=ccrs.PlateCarree())
ax3.set_extent([-130, -70, 20, 80], ccrs.PlateCarree())
ax3.coastlines(resolution='50m')
ax3.add_feature(countries_50m, edgecolor='k', linewidth=0.5)
ax3.add_feature(states_50m)
ax3.set_title('states_border')

ax4 = plt.subplot(2, 2, 4, projection=ccrs.PlateCarree())
ax4.set_extent([130, 140, 30, 40], ccrs.PlateCarree())
ax4.coastlines(resolution='10m')
ax4.add_feature(states_10m)
ax4.set_title('Japan')

f:id:ttmych:20160604155432p:plain

110m解像度のデータの一部は簡単に読み込めるように実装されている。

Name Description
cartopy.feature.BORDERS 国境線。
cartopy.feature.COASTLINE 主たる島を含む海岸線。
cartopy.feature.LAKES 湖。
cartopy.feature.LAND 大陸と主たる島。
cartopy.feature.OCEAN 海洋。
cartopy.feature.RIVERS 河川。
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

plt.figure(figsize=(8,8))

ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.BORDERS, edgecolor='gray')
ax.add_feature(cfeature.COASTLINE)

ax.set_extent([-20, 60,-40, 40], ccrs.PlateCarree())

f:id:ttmych:20160604155438p:plain

雨量予想をプロットしてみた

cartopy matplotlib python xray

xraycartopyで雨量予想をプロットする。ピークの値を表示してそれっぽくプロットする方法を調べてみた。

ピークの検出

以下のFAQが参考になった。

stackoverflow.com

隣接する格子でのピークを参照するためにScipyのndimage.filters.maximum_filterを使った。

maximum_filterで隣接格子の最大値でフィルターしたデータを用意して、元データと比較すればピークの格子を特定できる。

import numpy as np
import numpy.random as random
from scipy.ndimage.filters import maximum_filter

# ランダムデータ
data = random.rand(25,25)
# 各格子を隣接5格子のうち最大の値でフィルターして置き換える
local_max = maximum_filter(data, footprint=np.ones((5,5)), mode='constant')
# 元のデータとフィルターしたデータで値が同じ格子(ピーク格子)以外をマスクする
peak = np.ma.array(data, mask=~(data==local_max))

import matplotlib.pyplot as plt
plt.subplot(2,2,1)
plt.pcolormesh(data, vmin=0, vmax=1)
plt.title('Original')
plt.colorbar()

plt.subplot(2,2,2)
plt.pcolormesh(local_max, vmin=0, vmax=1)
plt.title('Filtered')
plt.colorbar()

plt.subplot(2,2,3)
plt.pcolormesh(peak, vmin=0, vmax=1)
plt.title('Peak')
plt.colorbar()

plt.show()

f:id:ttmych:20160117212342p:plain

プロットしてみる

メソモデルの1時間雨量の予想をプロットしてみた。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patheffects as PathEffects

import xray
import numpy as np

## select data
# 1時間雨量から東日本の範囲をxray.dataArrayで取得
ncfile = 'http://opendap.rish.kyoto-u.ac.jp/opendap/hyrax/jmadata/gpv/latest/MSM-S.nc'
ds = xray.open_dataset(ncfile)
darray = ds['r1h'].sel(lat=slice(38,32),lon=slice(133,143))[18,:,:]

## for annotation
# 数値プロットのために座標を2dのnumpy.ndarrayにしておく
ylab, xlab = darray.dims
xval = darray[xlab].values
yval = darray[ylab].values
xval, yval = np.meshgrid(xval, yval)

## detect peak
# 隣接15格子のピークかつR1=1mm以上の座標を算出
from scipy.ndimage.filters import maximum_filter

zval = darray.to_masked_array(copy=False)
local_max = maximum_filter(zval, footprint=np.ones((15,15)), mode='constant')
local_max[0,:]  = np.ma.masked
local_max[-1,:] = np.ma.masked
local_max[:,0]  = np.ma.masked
local_max[:,-1] = np.ma.masked
mask = (zval == local_max) & (zval >=1)

## for colormap
colors = ['#FFFFFF', '#00FFFF', '#00B0FF', '#0070FF', '#228B22', '#00FF00', '#FFFF00', '#FF8000', '#FF0000', '#FF00FF']
levels = [1,3,5,10,20,30,40,50,80]

## plot
plt.figure(figsize=(8,5))
ax = plt.axes(projection=ccrs.Mercator())
transform = ccrs.PlateCarree()._as_mpl_transform(ax)
darray.plot(ax=ax, levels=levels, colors=colors, extend='both', transform=transform)
for x, y, val in zip(xval[mask].flat, yval[mask].flat, zval[mask].flat):
    val = '{}'.format(int(val))
    ax.plot(x, y, '+', color='k', mew=1, ms=6, transform=transform)
    ax.annotate(val, xy=(x,y), xytext=(1, -1), xycoords=transform,
                textcoords='offset points',  ha='left', va='top', fontsize=10,
                path_effects=[PathEffects.withStroke(linewidth=2, foreground='w')])
ax.coastlines('10m')
plt.show()

f:id:ttmych:20160117212345p:plain

メモ

地図座標への変換はCartopyを使っている。matplotlibのannotateを使って変換する場合は次のようにすればいいらしい。

ax = plt.axes(projection=ccrs.Mercator())
transform = ccrs.PlateCarree()._as_mpl_transform(ax)
ax.annotate(val, xy=(x,y), xycoords=transform)

こちらのFAQを参考にした。 stackoverflow.com

xrayとcartopyでヒートマップ的なプロットを書いてみた

cartopy xray python matplotlib

github.com

github.com

xraycartopyでヒートマップ的なプロットをしてみた。雨量などの格子データのプロットに使いたい。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

import xray
import numpy as np
import numpy.random as random
from datetime import datetime

# toy weather data
lon = 136.5 + np.arange(20)*0.25
lat = 38. - np.arange(20)*0.25
time = datetime(2016,1,1,0)
data = random.rand(len(lat), len(lon))*100
rain = xray.DataArray(data, coords={'time': time, 'lat': lat, 'lon': lon}, dims=['lat', 'lon'])

# prepare data for annotation
darray = rain
ylab, xlab = rain.dims
xval = darray[xlab].values
yval = darray[ylab].values
zval = darray.to_masked_array(copy=False)
xval, yval = np.meshgrid(xval, yval)

# for colormap
colorlist = ['#FFFFFF', '#00FFFF', '#00B0FF', '#0070FF', '#228B22', '#00FF00', '#FFFF00', '#FF8000', '#FF0000', '#FF00FF']
levels = [1,3,5,10,20,30,40,50,80]

# plot
plt.figure(figsize=(8,6))
ax = plt.axes(projection=ccrs.PlateCarree())
rain.plot(ax=ax, levels=levels, colors=colorlist, extend='both')
for x, y, val in zip(xval.flat, yval.flat, zval.flat):
    val = '{}'.format(int(val))
    ax.text(x, y, val, ha='center', va='center')
ax.coastlines(resolution='10m')
plt.show()

f:id:ttmych:20160110150633p:plain

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が対応していない。
    • (英語のサイトも含めて)サンプルが少ない。

この記事のソースコード