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

今回も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()

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()