【Python】geopandasで描画したシェープデータをflaskでhtml上に画像出力する

シェープファイルをhtmlに表示させる方法を調べていくうちにPythonが持っているgeopandasという地図出力ライブラリだとshpの出力が可能だとわかり、シェープファイル（shp形式）を表示させてみました。

ところが、html上に表示させるためflaskを活用しようと思ったのですが、geopandas flaskで検索してもなかなか見つからず難航していたところ、shpファイルに対し、matplotlib内のcanvas機能を使って画像化する方法があることを知りました。そこで、flaskでその画像を表示できないか試行錯誤していくうちに、ひとまず方法を実現できたので、備忘録として記しておきます。

なお、このページではシェープデータは全部ポリゴン（図形）を使用しています（他にドット、ラインがあります）。

# 1:html上に画像化したshpを表示させる
画像を表示させる手順は以下の通りです
1. matplotlibライブラリを使って、シェープデータを落とし込む表を作成する。
1. シェープデータを読み込み、plotメソッドで描画し、表に落とし込む。
1. 落とし込んだ表をcanvasで画像化する。
1. 画像化した表をflaskでhtml上に展開させる。

それぞれの方法はいろんなページで解説があるのですが、geopandasで描画したシェープデータをFlaskで表示させるという方法がどこにもなく、試行錯誤を繰り返しました。ですが、

- **geopandasで描画したシェープファイルは画像化が可能（参考記事より）**
- **flaskで画像データの転送は可能**

それらの情報を頼って、作成しひとまず完成したソースが以下となります。

```py.view_shp.py
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import geopandas as gpd
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from io import BytesIO
#Flaskを呼び出す
from flask import Flask, render_template,make_response
app = Flask(__name__)

fig, ax = plt.subplots(figsize=(10, 5))
fig.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 図の周りの余白削除
ax.set_axis_off()  # xy の目盛り軸削除
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot(ax=ax) #表にシェープファイルを落とし込む

#ルーティング
#htmlのルーティング
@app.route('/')
def shapelist():
	png_output = BytesIO()
	canvas = FigureCanvasAgg(fig) #canvasで画像化
	canvas.print_png(png_output)
	data = png_output.getvalue() #バイト列の取得
	response = make_response(data)
	response.headers['Content-Type'] = 'image/png'
	response.headers['Content-Length'] = len(data)
	return response

if __name__ == "__main__":
 app.run(debug=False, host='0.0.0.0', port=5000)

```

こんな感じでhtml上にシェープデータを表示できました。

![SnapCrab_No-0000.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/be741889-0563-fd07-f0b6-80628a05ab81.jpeg)

# 2A:画像化したshpをhtml上のimgタグに表示する
画像表示はできたので、次はよりシステマティックに操作していくため、imgタグの中に埋め込んで表示できるようにします。それにあたっては、このページを参考にしました。

- [Matplotlib描画画像をFlaskによるWebアプリにてBytesIO を介して直接表示する方法](https://qiita.com/SatoshiTerasaki/items/8d79ec9d463bf0ce3595)

大事なポイントはbase64形式に変換が必要なことと、それをimgタグのsrcコード内に埋め込むことです。

```py.view_shp2.py
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import geopandas as gpd
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from io import BytesIO
import urllib
#Flaskを呼び出す
from flask import Flask, render_template
app = Flask(__name__)

fig, ax = plt.subplots(figsize=(10, 5))
fig.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 図の周りの余白削除
ax.set_axis_off()  # xy の目盛り軸削除
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot(ax=ax)

#ルーティング
#htmlのルーティング
@app.route('/')
def graph():
	png_output = BytesIO()
	canvas = FigureCanvasAgg(fig)
	canvas.print_png(png_output)
	img_data = urllib.parse.quote(png_output.getvalue()) #base64形式へのパース化
	return render_template("graph.html",title="shape描画テスト",img_data=img_data )

if __name__ == "__main__":
	app.run(debug=True, host='0.0.0.0', port=5000)
```

こうすれば、テンプレートに設定したimgタグのsrcプロパティに値が代入され、imgタグで表示できるようになります。

```html:graph.html
{% extends "layout.html" %}
{% block content %}
	<h1>shp画像の読込</h1>
	<img id="shape" src="data:image/png:base64,{{ img_data }}"></img>
{% endblock %}
```

![shp.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/39fa0e15-eeb2-7d93-d4bf-9a447793e3d4.jpeg)

# 2B:イベントから画像を呼び出してみる
これを更にAjaxで制御してみます。画像呼出ボタンを押下することで、shpファイルを表示させるようにします。FlaskでAjaxを使用する方法は、自分が以前作成した記事をそのまま活用しています。

- [【Python】Flask・DjangoでpartialなAjax（htmlを部品化してAjaxで制御）](https://qiita.com/BRSF/items/0fb6431cfb6acf81b028)

```py.view_shp3.py
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import geopandas as gpd
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from io import BytesIO
import urllib
#Flaskを呼び出す
from flask import Flask, render_template
app = Flask(__name__)

fig, ax = plt.subplots(figsize=(10, 5))
fig.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 図の周りの余白削除
ax.set_axis_off()  # xy の目盛り軸削除
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot(ax=ax)

#ルーティング
#index
@app.route('/')
def index():
	return render_template('index.html')

#Ajaxで部品をレスポンスする
@app.route('/_graph',methods=["POST"])
def graph():
	png_output = BytesIO()
	canvas = FigureCanvasAgg(fig)
	canvas.print_png(png_output)
	img_data = urllib.parse.quote(png_output.getvalue())
	return render_template("graph.html",title="shape描画テスト",img_data=img_data )

if __name__ == "__main__":
	app.run(debug=True, host='0.0.0.0', port=5000)
```

```html:index.html
{% extends "layout.html" %}
{% block content %}
<!-- ここにメインコンテンツを書く -->
<content class="box">
	<button id="bt">画像呼出</button>
	<article id="graph"></article>
</content>
{% endblock %}
```

```html:graph.html
	<h1>shp画像の読込</h1>
	<img id="shape" src="data:image/png:base64,{{ img_data }}"></img>
```

このように、ステータスがxhrになっていることがわかると思います。

![shp3.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/72ecf016-d0bf-3b7b-2195-918451770053.jpeg)

# 3:任意のシェープファイルを個別に描画し、画像出力する
任意のシェープファイルを表示させることができたので、今度はそれを個別に表示させてみます。shpファイルを分割して表示させる方法は以下のサイトを参考にしました。

- [How to split a single polyline shapefile into many shapefiles based on unique value using GeoPandas?](https://gis.stackexchange.com/questions/380724/how-to-split-a-single-polyline-shapefile-into-many-shapefiles-based-on-unique-va)

ここで、注意しなければいけないのは、逐一、図を初期化させることです。一見、回りくどいことを行っているようですが、これをやっておかないと、figにシェープが追記されてしまいます。

```py:shapes_split.py
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import geopandas as gpd
from matplotlib.backends.backend_agg import FigureCanvasAgg
from io import BytesIO
import urllib
#Flaskを呼び出す
from flask import Flask, render_template
app = Flask(__name__)

#表の描画クラス
class Shape:
	def set_ax(self):
		fig, ax = plt.subplots(figsize=(2, 2))
		fig.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 図の周りの余白削除
		ax.set_axis_off()  # xy の目盛り軸削除
		self.fig = fig
		return ax

	def draw(self):
		fig = self.fig
		png_output = BytesIO()
		canvas = FigureCanvasAgg(fig)
		canvas.print_png(png_output)
		img_data = urllib.parse.quote(png_output.getvalue())
		return img_data

path = 'tmp/hogehoge.shp'
df_shps = gpd.read_file(path,encoding='SHIFT-JIS')
df_shps_utm = df_shps.copy()
df_shps_utm.to_crs(epsg=3857)
shp = Shape()
img_datas = [] #ここの個別のシェイプを格納していく

for i in df_shps_utm.index:
	ax = shp.set_ax() #表を初期化する。この工程を省略するとaxは追記されてしまう。
	df_shps_utm.iloc[[i]].plot(ax=ax,color='limegreen',edgecolor='k',linewidth=0.6)
	data = shp.draw()
	img_datas.append(data)


#ルーティング
@app.route('/')
def graph():
	return render_template("graph.html",title="shape描画テスト",img_datas=img_datas )

if __name__ == "__main__":
	app.run(debug=True, host='0.0.0.0', port=5000)

```

```html:graph.html
	<h1>shp画像の読込</h1>
	{{path}}
	{% for img_data in img_datas %}
	<img class="shape" src="data:image/png:base64,{{ img_data }}"></img>
	{% endfor %}
```

こんな感じで、個別のシェープファイルが表示されるようになりました。

![split.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/a642d52b-6fc5-6a6a-0973-cdc52412a2cc.jpeg)

## 任意のシェープファイルを描画したい場合
任意のシェープファイルを描画して画像化する場合はシェープ描画の部分を以下のように記述すると大丈夫のようです。

```py:shapes.py
fig, ax = plt.subplots(figsize=(10, 5)) #表の作成
fig.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 表の周りの余白削除
ax.set_axis_off()  # xy の目盛り軸削除
path = 'hogehoge.shp'
df_shps = gpd.read_file(path,encoding='SHIFT-JIS')
df_shps_utm = df_shps.copy()
df_shps_utm.to_crs(epsg=3857) #測地系の調整
df_shps_utm.plot(ax=ax,color='none',edgecolor='k',linewidth=0.6) #シェイプを表に落とし込む
```

# 4A:シェープデータを地図に合わせる
任意のシェープファイルが表示できるようになったので、今度はそれをエリアの地図に合わせて表示していこうと思います。…その前に、大事なポイントがあり、地理学やGISでの知識となりますが、**測地系**をよく把握しておく必要があります。これを把握しておかないと、併合作業はうまくいきません。

## 測地系について
測地系とはなにか、その定義ですが、某測量大手の説明文を借りると

> 「緯度経度の座標軸を使って、地図画面上の特定の位置を示す際の基準となる前提条件」

のことです。

もっと平たくいえば地球は球面であるため、平面表示させるようにすると、どうしても面積や方角などで大小の歪みが生じます。そこでその歪みに対し、距離を極力合わせる、方角を極力合わせるなど用途、目的に合わせた測量方法によって、様々な平面化の方法が提唱され、それらは測地系として体系化されるようになりました。

また、地球全体でなくとも、自分の国だけあればいい、あるいは国内のその地域だけあればいい（エリアを狭めたほうが歪みは少なくなるため）場合も存在するため、日本の場合は更にエリアを限定したJGD2000、JGD2011といった日本測地系というものも用いられます（測量会社のCADデータなどは、この日本測地系を用いることが多いです）。

しかし、注意しなければいけないのは、複数のデータに対し、測地系が異なるということは、**座標情報が異なるために同じ場所であっても、想定通りの位置を示さなくなる**ということです。もっと噛み砕いていえば、測地系の対象となるエリアのサイズは大小様々、また基準点も測地系によって異なります。そして、右上の座標（Y=0,X=0）を基準点として位置を測定します。そのため、基準点が違えば、そこから同じ距離だけ移動したとしても、場所が全然違うのは一目瞭然です（例え話ですが、マンハッタンを基準点としてそこから南西に100キロ移動した場合の地点と、東京タワーを基準点として南西に100キロ移動した場合の地点を、同じ座標情報（y=100,x=100）と比較するようなものです）。

今回は測地系がEPSG:4612（測地座標系）の日本地図に、測地系がEPSG:3857（球面メルカトル）のローカルシェープデータを併合する作業となりますが、**前者が日本だけを基準として平面化**しているものに対し、**後者は世界全体を基準として平面化**しているため、そうなれば当然位置と距離に大幅なズレが生じます。そこで、その座標のズレを調整し再計算する処理が、以下の記述になります。

```py
# ポリゴンの日本地図を呼び出す
pref_poly = [Polygon(points) for points in pref_points(get_data())]
gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly) #測地系はEPSG:4612 
# ローカルのshp情報
shp_path = 'hogehoge.shp'
geos = gpd.read_file(shp_path,encoding='SHIFT-JIS')
geos.crs = f'epsg:{3857}' #元に指定された測地系（球面メルカトル）。
geos = geos.to_crs(epsg=4612) #変換された測地系（測地座標系）
```

## 併合したシェープデータを画像出力
それを踏まえた上で、日本地図にローカルのシェープ情報を併合させたプログラムが以下になります。

```py:shapes_bind.py
from flask import Flask, render_template,make_response
app = Flask(__name__)

from japanmap import pref_names, get_data, pref_code,pref_points
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
import geopandas as gpd
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from io import BytesIO
import urllib
# 表示用のfigure作成
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
fig.subplots_adjust(left=0,right=1,top=1,bottom=0)
ax.set_axis_off()
# 日本地図のポリゴンデータ作成しGeoDataFrameに格納
pref_poly = [Polygon(points) for points in pref_points(get_data())]
gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly)
gdf_pref['prefecture'] = pref_names[1:]  # 県名を格納
gdf_pref['pref_code'] = gdf_pref['prefecture'].apply(lambda x: pref_code(x))  # 県コードを格納
# 県コードを指定
gdf_pref = gdf_pref[gdf_pref['pref_code'] == xx] #xxは任意の都道府県番号（例：京都府なら26）
# 日本地図をプロット
gdf_pref.plot(ax = ax,color = 'gray')  # 塗りつぶし色を指定
#ローカルのシェープ情報
shp_path = 'hogehoge.shp'
geos = gpd.read_file(shp_path,encoding='SHIFT-JIS')
geos.crs = f'epsg:{3857}' #出力シェイプの測地系
geos = geos.to_crs(epsg=4612) #変換先の測地系（日本地図に合わせる）
geos.plot(ax=ax,color='b')
png_output = BytesIO()
canvas = FigureCanvasAgg(fig)
canvas.print_png(png_output)
img_datas = urllib.parse.quote(png_output.getvalue()) #base64形式へのパース化

#ルーティング
@app.route('/')
def graph():
	return render_template("shape.html",title="shape描画テスト",shapes=img_datas )

if __name__ == "__main__":
	app.run(debug=True, host='0.0.0.0', port=5000)
```

```html:shape.html
{% extends "layout.html" %}
{% block content %}
	<img class="shape" src="data:image/png:base64,{{ shapes }}"></img>
{% endblock %}
```

こんな感じに、ローカルのshpデータが日本地図に埋め込まれました(エリアは部分的に切り取っています)。

![test.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/cda168f8-0da8-979f-526b-429f77bb714c.jpeg)

# 4B:個別のローカルシェープデータをエリアの地図に併合して表示させる
では、上記のシェープデータを個別に表示させる方法とシェープデータを地図に併合させる方法を合わせてみます。

```py:shp_bind_split.py
from japanmap import pref_names, get_data, pref_code,pref_points
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
import geopandas as gpd
from matplotlib.backends.backend_agg import FigureCanvasAgg
from io import BytesIO
import urllib
#Flaskを呼び出す
from flask import Flask, render_template
app = Flask(__name__)

class Shape:
	def read_shape(self):
		path = 'hogehoge.shp'
		geos = gpd.read_file(path,encoding='SHIFT-JIS')
		geos.crs = f'epsg:{3857}'
		geos = geos.to_crs(epsg=4612)
		return geos
	
	def set_ax(self):
		fig, ax = plt.subplots(1, 1, figsize=(4, 4))
		fig.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 図の周りの余白削除
		ax.set_axis_off()  # xy の目盛り軸削除
		self.fig = fig
		self.ax = ax
	
	def area(self):
		# 日本地図のポリゴンデータ作成しGeoDataFrameに格納
		ax = self.ax
		pref_poly = [Polygon(points) for points in pref_points(get_data())]
		gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly)
		gdf_pref['prefecture'] = pref_names[1:]  # 県名を格納
		gdf_pref['pref_code'] = gdf_pref['prefecture'].apply(lambda x: pref_code(x))  # 県コードを格納
		gdf_pref = gdf_pref[gdf_pref['pref_code'] == xx] #xxは任意の都道府県コード（例：新潟県なら15）
		gdf_pref.plot(ax=ax,color='gray')
	
	def local(self,loc):
		ax = self.ax
		loc.plot(ax=ax,color='limegreen',edgecolor='k',linewidth=0.6)
		
	def draw(self):
		fig = self.fig
		png_output = BytesIO()
		canvas = FigureCanvasAgg(fig)
		canvas.print_png(png_output)
		img_data = urllib.parse.quote(png_output.getvalue())
		return img_data

Shp = Shape()
geos = Shp.read_shape() #シェープ情報の設定
img_datas = []
for i in geos.index:
	Shp.set_ax() #表を逐一作成させる
	Shp.area() #エリア指定
	Shp.local(geos.iloc[[i]]) #ローカルシェープ情報を落とし込み
	data = Shp.draw() #描画
	img_datas.append(data)

#ルーティング
@app.route('/')
def graph():
	return render_template("graph.html",title="shape描画テスト",img_datas=img_datas )

if __name__ == "__main__":
	app.run(debug=True, host='0.0.0.0', port=5000)
```

このように、エリアごとにポリゴンが埋め込まれていることがわかると思います。

![graph.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/7ed61193-8c22-9904-1432-8575535f7c25.jpeg)

ところが、これだと地図に合わせて表示させているので、ポリゴンが芥子粒みたいになってしまっています。今度はローカルシェープのサイズに合わせて地図を切り取って表示させたいと思います。

# 5：画像化したシェープデータを一定のサイズに切り取る
matplotlibで作成した表をトリミングする方法を調べてみたところ、それを一度画像化する必要があるみたいです。ですが、今までのようにFigureCanvasAggだとトリミングできないようなので、方法を調べたところ、pillow(PIL)で画像化してから再編集する方法に行き着きました。

それで作成してみたプログラムが以下の部分です。なお、画像切り取り、リサイズの方法は以下のページを参考にしています。

- [Python, Pillowで画像の一部をトリミング（切り出し/切り抜き）](https://note.nkmk.me/python-pillow-image-crop-trimming/)

```py:shape_trime.py
from flask import Flask, render_template,make_response
app = Flask(__name__)
from japanmap import pref_names, get_data, pref_code,pref_points
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
import geopandas as gpd
from io import BytesIO
import base64
from PIL import Image
import urllib
# 表示用のfigure作成
fig, ax = plt.subplots(1, 1,figsize=(20, 20))
ax.set_axis_off()
# 日本地図のポリゴンデータ作成しGeoDataFrameに格納
pref_poly = [Polygon(points) for points in pref_points(get_data())]
gdf_pref = gpd.GeoDataFrame(crs = 'epsg:4612', geometry=pref_poly)
gdf_pref['prefecture'] = pref_names[1:]  # 県名を格納
gdf_pref['pref_code'] = gdf_pref['prefecture'].apply(lambda x: pref_code(x))  # 県コードを格納
gdf_pref = gdf_pref[gdf_pref['pref_code'] == xx] #都道府県コードを指定。秋田県なら5。
# 日本地図をプロット
gdf_pref.plot(ax = ax,color = 'gray')  # 塗りつぶし色を指定
shp_path = 'hogehoge.shp'
geos = gpd.read_file(shp_path,encoding='SHIFT-JIS')
geos.crs = f'epsg:{3857}'
geos = geos.to_crs(epsg=4612)
geos.plot(ax=ax,color='green')
bytes = BytesIO() #原図の保存用
fig.savefig(bytes,format='png')
im = Image.open(bytes) #原図をPILで開く
bytes2 = BytesIO() #画像編集用の設定
im = im.crop((680,1200,900,1260)) #部分的に画像を切り取る
im = im.resize((im.width * 4,im.height * 4)) #切り取った画像を拡大
im.save(bytes2,'png')
data = urllib.parse.quote(bytes2.getvalue())
```

Pythonは全く詳しくないので原理はわからないのですが、matplotlibとPILは別々に使用するので、それに伴って**BytesIOを二回に分けて呼び出す必要がある**みたいです。そして、1回目はmatplotlibで作成した表に対し、画像化処理。次はPILで開いた画像に対する編集処理で、PILで呼び出した画像ファイルを編集していけば、適宜画像データを切り取れるみたいです。

![trim.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/d811f040-237f-6436-2e09-449ba9796328.jpeg)

これをループ化させれば、同じエリア内ならば、このように個別に作成することもできます。

![SnapCrab_No-0002.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272191/0d481356-9aaf-f54e-aa54-4b9da2a39ca9.jpeg)

ただ、これだと切り抜き位置を決め打ちで指定しているのと、極力シェープデータが崩れないように元の表（Figureサイズが20）をかなり大きく取っているので、読み込みに時間が掛かってしまっています。なので、これをピンポイントで切り取れないか方法を試行錯誤中です。

# 引き続き、やりたいこと
目的としてはまだ途中の段階ですが、今後はこのような実装を予定しており、成功した際には記事を加筆していく予定です。

- 任意のshpファイルに格納された個別のシェイプ情報を基準として、エリアの地図に埋め込んで表示


# 参考記事
- [geopandas で 世界地図を表示](https://qiita.com/XPT60/items/27911fc3190f3743e9d5)
- [shapeファイルを画像に変換して出力したい](https://teratail.com/questions/257556)
- [MatplotlibをFlaskで拡張して、誰でもPythonで作成したグラフを見られるようにする](https://qiita.com/kinuta_masaya/items/9bea70ade73cc48cbdf0)
- [Flask matplotlibを使ってグラフを表示する方法（JavaScriptのChart.jsを使ってグラフを描画）](https://start-python.hateblo.jp/entry/2020/01/27/090000)

ほか