cgiからmatplotlibで画像を表示させる
簡単な部内サイトなどを作る時、matplotlibで書いた画像を動的に表示させたい場合cgiスクリプトで動かしたい場合があると思う。 matplotlibの標準出力のサンプルを参考に、レスポンスヘッダを付加してすれば良いと思って、
import sys import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt plt.plot([1, 2, 3]) print('Content-type: image/png') print() plt.savefig(sys.stdout.buffer, format='png')
のように書いたところ、エラーで表示されない。Apacheのログをみると
malformed header from script 'image.cgi': Bad header: \x89PNG
となっていた。ヘッダが悪い？と思っていろいろ調べてみると、stackoverflowにそれっぽい回答が。print関数の出力が標準出力にうまく返されていない？のが原因らしい。
回答の通り以下のように書き換えたらうまく表示された。
import sys import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt plt.plot([1, 2, 3]) print('Content-type: image/png') print() sys.stdout.flush() plt.savefig(sys.stdout.buffer, format='png'))
もしくは、
import sys import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt plt.plot([1, 2, 3]) sys.stdout.buffer.write(b'Content-type: image/png\n\n') plt.savefig(sys.stdout.buffer, format='png'))
少し前に試した時は、こんなテクニック必要なかったような気がするけどバージョン間の違いだろうか。
開発環境はpyenv & minicondaで仮想環境を構築。バージョンは、python=3.6.2, matplotlib=2.1.0, Apache=2.4
GeoPandasで国土数値情報を加工してみた 2
引き続き、地理情報データをPandasで扱うための拡張ライブラリgeopandasを利用して国土数値情報データを加工してみた。
引き続き、国土数値情報データのダウンロードページの土地利用区分の3次メッシュデータを利用する。今回は、2次メッシュに変換したときの土地利用区分を求めてみる。
まずは、データの読み込みを行う。
# 処理したいshapeファイルのリスト shpfiles = [shpfile for shpfile in glob.glob('./*/*jgd*/*.shp')] # shapeファイルを読み込んで結合させる df = pd.concat([gpd.read_file(shpfile, encoding='shift-jis') for shpfile in shpfiles]) # 後の処理のため3次メッシュを明示的に整数にしておく df['メッシュ'] = df['メッシュ'].apply(int)
前回は、3次メッシュのまま表示させていたので気にならなかったが、このままだと3次メッシュを結合させる際にPolygonの境界が小数点以下でずれていることが問題となる。元データのせいなのか、Pythonでの読み込みの際の型変換が問題なのかはわからないが、隣接格子で微妙に座標がずれてしまいPolygonを連結させることができなかった。
これを修正するために、Polygonの座標を四捨五入して小数点以下を整えることにする。以下のエントリを参考にした。
from shapely.geometry import shape, mapping import numpy as np def set_precision_shape(shape_in, precision): try: geojson = mapping(shape_in) coords = geojson['coordinates'] geojson['coordinates'] = np.round(np.array(geojson['coordinates']), precision) except: return shape_in return shape(geojson) df['geometry'] = df['geometry'].apply(lambda g: set_precision_shape(g, 6))
次に2次メッシュでの最大値を求める。geometryの集約については、geopandas.dissolveで処理できる。
# 2次メッシュに変換 df['メッシュ'] = df['メッシュ'].apply(lambda x: x//100) # geometryを集約 df_geo = df[['メッシュ','geometry']].dissolve(by='メッシュ') # geometry以外の要素の合計を求める df_attr = df.groupby('メッシュ').sum() # メッシュコードをキーにして結合させる df_mesh2 = pd.concat([df_geo, df_attr], axis=1).reset_index()
あとは前回と同じく、メッシュ毎の最大の利用区分を求めてshapeファイルに出力する。
cols = [e for e in df_mesh2.columns.tolist() if e not in ['geometry', 'メッシュ']] df_mesh2['landtype'] = df_mesh2[cols].apply(lambda x: x.argmax(), axis=1) df_mesh2 = df_mesh2[['メッシュ','landtype','geometry']]
QGISで表示すると以下のようになった。2次メッシュへの変換はうまく計算できているようだ。
GeoPandasで国土数値情報を加工してみた1
地理情報データをPandasで扱うための拡張ライブラリgeopandasを利用して国土数値情報データを加工してみた。
国土数値情報データのダウンロードページの土地利用区分の3次メッシュデータを利用する。このデータには各メッシュにおける土地利用区分毎の面積がデータとして格納されている。 今回は、このデータを加工して各メッシュで最も利用面積の大きな利用区分を求めたデータを作成してみる。また、国土数値情報のファイルは1次メッシュ毎にファイル分割されていて県別になっていないため、神奈川県に含まれるメッシュのみ抽出してみた。
まずは、データの読み込みを行う。
# 処理したいshapeファイルのリスト shpfiles = [shpfile for shpfile in glob.glob('./*/*jgd*/*.shp')] # shapeファイルを読み込んで結合させる df = pd.concat([gpd.read_file(shpfile, encoding='shift-jis') for shpfile in shpfiles]) # 後の処理のため3次メッシュを明示的に整数にしておく df['メッシュ'] = df['メッシュ'].apply(int)
読み込んだDataFrameには、shapeファイルのshapeデータ（polygon、pointなど）と対応する属性が格納されている。
総務省統計局のホームページ で公開されている、都道府県毎の3次メッシュコードのcsvファイルを利用して、神奈川県に該当するメッシュのみ取り出してみる。
# 神奈川県のメッシュコード一覧を読み込み df_city = pd.read_csv('csv/city-mesh_14.csv', encoding='shift-jis', dtype={'基準メッシュコード': np.int64}) # 3次メッシュコードが一覧に含まれるレコードのみ抽出 df = df[df['メッシュ'].isin(df_city['基準メッシュコード'].tolist())]
各メッシュ毎に最大面積をもつ利用区分を求める。利用区分名が列名になっているので行ごとに処理する。 その他の列は不要なので、3次メッシュコード、土地利用区の列を取り出す。
cols = [e for e in df.columns.tolist() if e not in ['geometry', 'メッシュ']] df['landtype'] = df[cols].apply(lambda x: x.argmax(), axis=1) df = df[['メッシュ','landtype','geometry']]
最後にDataFrameをshapeファイル形式で出力する。dbfファイル、shxファイルも同時に出力される。
df.to_file('土地利用区分.shp', driver='ESRI Shapefile', encoding='shift-jis')
出力したshapeファイルをQISで表示すると以下のようになる。
QGISでguiから編集・加工できるが、pythonで加工できるようにしておくとまとめて加工するときはスクリプトで自動化できて便利だと思う。Pandasベースなので統計処理とかも応用ができるしいろいろ試していきたい。
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()
Cartopyで地理データを可視化する2
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')
様々なベクトルデータを描画する
| 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')
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())
雨量予想をプロットしてみた
xrayとcartopyで雨量予想をプロットする。ピークの値を表示してそれっぽくプロットする方法を調べてみた。
ピークの検出
以下のFAQが参考になった。
隣接する格子でのピークを参照するために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) # 各格子を隣接５格子のうち最大の値でフィルターして置き換える 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()
プロットしてみる
メソモデルの１時間雨量の予想をプロットしてみた。
import cartopy.crs as ccrs import matplotlib.pyplot as plt import matplotlib.patheffects as PathEffects import xray import numpy as np ## select data # １時間雨量から東日本の範囲を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()
メモ
地図座標への変換は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でヒートマップ的なプロットを書いてみた
xrayとcartopyでヒートマップ的なプロットをしてみた。雨量などの格子データのプロットに使いたい。
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()
