はじめに
本記事は、地理情報を扱うために学習した内容をまとめた記事になります。
学習教材として、ヘルシンキ大学の Automating GIS-processes 2018 という講座を利用しました。
Google 翻訳を駆使すれば、問題なく演習が可能です。
この講座を受講すれば、『Python』 を使って『無料』で以下のようなことができるようになります。
人口増減率が大きいエリアの導出
営業マンが効率よく回れるルートの導出
新規出店に適したエリアの導出
※上記画像は esri ジャパン様 の画像を引用しています。
※参考のためイメージを載せています。
※本記事に誤りがあれば、ご指摘を頂けると幸いです。
本記事の対象者
以下の条件を満たしていること
- Dockerfile を作成してイメージ作成が可能なこと
- Jupyter Notebook の利用経験があること
- Python での GISデータハンドリングに興味があること
環境
- OS: Windows10 Pro(64ビット)
- CPU: Intel(R) Core(TM) i7-2600@3.40GHz 3.40GHz
- メモリ: 16.0GB
- ツール
- Docker for Windows: 2.0.0.3
GIS とは
GIS（ジー アイ エス）とは、Geographic Information System の略称です。
日本語では地理情報システムと訳されます。
地球上に存在する地物や事象はすべて地理情報と言えますが、これらを
コンピューターの地図上に可視化して、情報の関係性、パターン、傾向を
わかりやすいかたちで導き出すのが GIS の大きな役割です。
GIS の詳細は ESRI ジャパン GIS 基礎解説 を参照してください。
「GIS(地理情報システム)とは」「GIS 基礎解説」を読めば、GIS の概要を把握できます。
GIS データフォーマット
1.シェープファイル
シェープファイルは GIS データ フォーマットの 1 つです。
病院などの目標物や道路や建物などの位置や形状、属性情報を持つ
ベクターデータ（ポイント、ライン、ポリゴン）を格納することができます。
ベクターデータについては後述します。
画像はESRI ジャパンより引用
シェープファイルは複数のファイルから構成されています。
必須のファイルは ３ つあり、このうち 1 つでも欠けると
GIS アプリケーション上でシェープファイルと認識することができません。
また、シェープファイルを構成するファイルのサイズは
それぞれ 2GB のサイズ制限があります。
| 拡張子 | 概要 |
|---|---|
| .shp | エリアや地点を表す図形の情報を格納するファイル（必須） |
| .dbf | 図形の属性情報(国名・県名など)を格納するファイル（必須） |
| .shx | .shp と .dbf の対応関係を示すファイル（必須） |
| .prj | 図形の持つ座標系の定義情報を格納するファイル |
データ処理をする上で、シェープファイルを「GeoJSON」「TopoJSON」という
データ形式に変換することが多々あります。
GIS データを Web アプリで活用する場合は、特にこれらの形式で扱われることが多いです。
本記事で使用するfoliumライブラリ(詳細は後述)では「GeoJSON」に変換することで、
容易にインタラクティブマップを作成可能となります。
TopoJSON は利用したことがないのですが、以下に説明を記載しておきます。
2.GeoJSON
JavaScript Object Notation (JSON) を用いて空間データをエンコードし
非空間属性を関連付けるファイルフォーマットです。
属性にはポイント（住所や座標）、ライン（各種道路や境界線）、 ポリゴン（国や地域）などが含まれます。
ファイルフォーマット
{ "type": "FeatureCollection",
"features": [
{ "type": "Feature",
"geometry": {"type": "Point", "coordinates": [102.0, 0.5]},
"properties": {"prop0": "value0"}
},
{ "type": "Feature",
"geometry": {
"type": "LineString",
"coordinates": [
[102.0, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]
]
},
"properties": {
"prop0": "value0",
"prop1": 0.0
}
},
]
}
3.TopoJSON
地理空間情報をエンコードするためのフォーマットである GeoJSON の拡張形式です。
詳細はこちらを参照してください。
TopoJSON ファイルは GeoJSON よりずっと軽量になります。
上のアメリカ合衆国のシェイプファイルは GeoJSON の場合 2.2M ですが、TopoJSON では、
単純化していない境界線メッシュ形式でも 436K、すなわち 80.4% サイズが削減されます。
TopoJSON では描写効率も改善されます。共有点を一度しか描かなくて済むからです。
さらなるファイル容量削減のために、TopoJSON は整数座標のデルタエンコーディングのために、
浮動小数点ではなく固定小数点を用います。これにより、正確性を犠牲にすることなく
座標値の精度の丸め処理を省略できます（例：LilJSON )。また、GeoJSON と同様、
ファイルは簡単にテキストエディタで編集したり、gzip で圧縮することができます。
GISデータモデル
空間的な情報についてはさまざまな表現形態があり、それらに応じた GIS データモデルがあります。
GISデータモデルには、ベクターデータ、ラスターデータ、3Dデータ、といった種類があります。
ここでは、本記事で使用するベクターデータについて説明します。
ベクターデータとは
現実世界に存在する地物（目に見えないものも含む）をの3つの要素で表現したものです。
格納される座標値は、緯度・経度、あるいは地図投影法により定義された値です。
詳細は esri ジャパンのページ を参照して下さい。
| 種類 | 説明 | 代表 |
|---|---|---|
| ポイント（点） | 図形情報を（x, y）の座標値で格納 | 電柱、信号、ATM |
| ライン（線） | （x, y）の座標値を結んだ線分として格納 | 鉄道路線、排水管 |
| ポリゴン（面） | （x, y）の座標値を結んだ閉じた線分として格納 | 建物、行政界 |
画像は Automating GIS-processes 2018 より引用
環境構築手順
以下の2ファイルを利用して、GIS データハンドリング用の Python 環境を構築します。
| ファイル名 | 用途 |
|---|---|
| Dockerfile | コンテナのビルド |
| env_name.yml | 環境設定ファイル |
1. Dockerfile の準備
FROM continuumio/miniconda3:4.5.12
COPY env_name.yml /tmp/
RUN conda env create -f=/tmp/env_name.yml
RUN echo "source activate gis" > ~/.bashrc
ENV PATH /opt/conda/envs/gis/bin/:$PATH
WORKDIR /opt/notebooks
EXPOSE 8888
CMD jupyter notebook --notebook-dir=/opt/notebooks --ip=0.0.0.0 --port=8888 --no-browser --allow-root --NotebookApp.token='gis'
2. 環境設定ファイルの準備
Anacondaでは、環境設定をYAML形式で書き出すことで、そこから簡単に環境を再構築することができます。
使用するライブラリは このページ を参考にしています。
conda または pip で、versionを指定せずに1つずつインストールすると、ライブラリ間の依存関係が壊れます。
依存関係を解決するために、ここを参考にしました。(他に良い方法があったら教えてください)
name: gis
channels:
- conda-forge
- defaults
dependencies:
- geopandas==0.4.0
- geoplot==0.2.3
- geoviews==1.6.2
- haversine==2.0.0
- osmnx==0.9
- pysal==1.14.4
- contextily==1.0rc1
- rasterio==1.0.18
- rasterstats==0.13.0
- shapely==1.6.4
- scikit-learn==0.20.2
- lightgbm==2.2.2
- pip:
- dash==0.19.0
- dash-core-components==0.14.0
- dash-html-components==0.8.0
- dash-renderer==0.11.1
- pycrs==1.0.0
- mplleaflet==0.0.5
- jupyter==1.0.0
- bokeh==1.0.4
- vega_datasets==0.7.0
3. Jupyter Notebook の起動
- PowerShell を起動。 Dockerfile と env_name.yml が格納されているディレクトリで以下のコマンドを実行します。
docker build -t gis .
docker run -d --rm --name gis -p 8888:8888 -v C:/path//:/opt/notebooks gis
proxy 環境の場合は以下のコマンドを実行
docker build --build-arg http_proxy=http://proxysrv:port --build-arg https_proxy=http://proxysrv:port --no-cache -t gis .
docker run -d --rm --name gis -p 8888:8888 -v C:/path/:/opt/notebooks gis
※http://proxysrv:port の proxysrv と port は適宜設定すること。
※C:/path/に関してもボリュームしたいディレクトリを指定すること。
- ブラウザを起動して以下のリンクを開くと Jupyter Notebook が起動
GISデータハンドリングで利用するライブラリ
特に geopandas, folium, osmnx は GIS データのハンドリングで利用されるので、
利用方法を抑えておきたいです。
| ライブラリ名 | 説明 |
|---|---|
| geopandas | 地理空間データを操作するのが簡単になり、Pandasの機能と形の組み合わせが可能 |
| folium | データを地図上に可視化するライブラリ |
| osmnx | OpenStreetMapから街路ネットワークを取得、構築、分析、および視覚化する |
| shapely | 平面幾何学的オブジェクトの操作と分析のためのPythonパッケージ |
| fiona | 空間データの読み書き（geopandasの代替） |
| pyproj | 地図作成変換と測地計算を実行（PROJ.4に基づく） |
| geopy | Geocoding library：住所への座標 <->座標への住所 |
講座内容
全部で7講座あります。
私は Lesson1~6 を翻訳しながら写経しました。
Lesson7 は業務上不要だったため実施していません。
また、不明点があったときは、都度ググりながら調べて進めました。
| Lesson | Learning goals |
|---|---|
| 1 | ・Python で GIS を実行するためにどのようなツールが利用可能か分かること ・geometry オブジェクトが利用可能で、それらの用途を理解できること ・shapely を使用してさまざまな種類の形状を作成できること ・ファイルから座標を読み取り、それらに基づきポイントを作成できること |
| 2 | ・シェープファイルの Read/Write ができること ・GeoDataFrame に geometry を作成できること ・データの座標系を変更できること |
| 3 | ・住所 ⇔ ポイントの変換ができること ・Point in Polygon（PIP）クエリを実行できること ・KMLファイルを読み込めること ・レイヤー間で空間結合ができること ・特定の点からの最近傍点を探せること |
| 4 | ・空間データの分類ができること ・オーバーレイ解析ができること ・ジオメトリオブジェクトの集約ができること ・ジオメトリの単純化ができること |
| 5 | ・Geopandasを使用して背景ベースマップで静的マップを作成できること ・Foliumを使用して簡単なインタラクティブマップを作成できること ・GitHubページを使ってインターネット上で地図（静的/インタラクティブ）を共有できること |
| 6 | ・Pythonを使ってOpenStreetMapからデータを取得して保存できること ・単純な街路ネットワークの特性と統計を抽出できること ・osmnx / networkxで最短経路アルゴリズムを使用して単純な経路最適化ができること |
| 7 | ・QGIS Processing Toolboxでスクリプト用のシンプルなユーザーインターフェースを作成できること ・QGISアルゴリズムを実行できること ・PyQGISとQGIS Processingのドキュメントをオンラインで読む方法を |
GIS データハンドリング
ヘルシンキ大学の講座から、私の仕事で使いそうな部分をピックアップしました。
また、自分で調べた内容も一部記載しています。
| # | 見出し |
|---|---|
| 1 | シェープファイルのread/write |
| 2 | 座標参照系の設定 (Map projections) |
| 3 | Point in Polygon(PIP) |
| 4 | 点間の距離計算 |
| 5 | ポイント ⇔ 住所(名称) の変換 |
| 6 | 最近傍点の探索 |
| 7 | オーバーレイ解析 |
| 8 | インタラクティブマップの作成 |
1. シェープファイルのread/write
1.1. データ準備
L2_data.zip をダウンロードします。
zip を展開すると以下のような構成になっています。
L2_data
├ DAMSELFISH_distributions.cpg
├ DAMSELFISH_distributions.dbf
├ DAMSELFISH_distributions.prj
├ DAMSELFISH_distributions.shp
├ DAMSELFISH_distributions.shx
├ Europe_borders.cpg
├ Europe_borders.dbf
├ Europe_borders.prj
├ Europe_borders.sbn
├ Europe_borders.sbx
├ Europe_borders.shp
└ Europe_borders.shx
1.2. read
シェープファイルを読み込みます。
geometry 列にはベクターデータが格納されています。
import geopandas as gpd
fp = 'L2_data/DAMSELFISH_distributions.shp'
data = gpd.read_file(fp)
display(data.head(2))
| ID_NO | BINOMIAL | ORIGIN | COMPILER | YEAR | CITATION | SOURCE | DIST_COMM | ISLAND | SUBSPECIES | ... | RL_UPDATE | KINGDOM_NA | PHYLUM_NAM | CLASS_NAME | ORDER_NAME | FAMILY_NAM | GENUS_NAME | SPECIES_NA | CATEGORY | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 183963.0 | Stegastes leucorus | 1 | IUCN | 2010 | International Union for Conservation of Nature... | None | None | None | None | ... | 2012.1 | ANIMALIA | CHORDATA | ACTINOPTERYGII | PERCIFORMES | POMACENTRIDAE | Stegastes | leucorus | VU | POLYGON ((-115.6437454219999 29.71392059300007... |
| 1 | 183963.0 | Stegastes leucorus | 1 | IUCN | 2010 | International Union for Conservation of Nature... | None | None | None | None | ... | 2012.1 | ANIMALIA | CHORDATA | ACTINOPTERYGII | PERCIFORMES | POMACENTRIDAE | Stegastes | leucorus | VU | POLYGON ((-105.589950704 21.89339825500002, -1... |
1.3. write(読み込んだファイルから書き込み)
さきほど読み込んだファイルから50行抽出したものをファイル出力します。
outfp = "L2_data/DAMSELFISH_distributions_SELECTION.shp"
selection = data[0:50]
selection.to_file(outfp)
ファイルが出力されていることを確認できます。
L2_data
├ ・・・
├ DAMSELFISH_distributions_SELECTION.cpg
├ DAMSELFISH_distributions_SELECTION.dbf
├ DAMSELFISH_distributions_SELECTION.prj
├ DAMSELFISH_distributions_SELECTION.shp
├ DAMSELFISH_distributions_SELECTION.shx
├ ・・・
└ Europe_borders.shx
1.4. write(ゼロベース)
GeoDataFrame を作ってシェープファイルの書き込みをします。
import geopandas as gpd
from shapely.geometry import Point, Polygon
# 列「geometry」を追加
newdata = gpd.GeoDataFrame()
newdata['geometry'] = None
# 列「geometry」に Polygon オブジェクトを設定
coordinates = [(24.950899, 60.169158), (24.953492, 60.169158), (24.953510, 60.170104), (24.950958, 60.169990)]
poly = Polygon(coordinates)
newdata.loc[0, 'geometry'] = poly
# シェープファイルの出力
outfp = "L2_data/Senaatintori.shp"
newdata.to_file(outfp)
display(newdata.head())
| geometry | |
|---|---|
| 0 | POLYGON ((24.950899 60.169158, 24.953492 60.16... |
2. 座標参照系の設定 (Map projections)
2.1. CRSとは
画像は Automating GIS-processes 2018 より引用
GeoDataFrame内のベクターデータは、任意の空間内の座標の集合であるため、
座標参照系（CRS）の設定が必要になります。
CRS とは 「Coordinate reference system」の略称です。
CRSは、座標が地球上の場所とどのように関連しているかをPythonに伝えます。
地図投影法（または投影座標系）は、緯度と経度を平面単位に体系的に変換したもので、
単位は通常（メートル単位ではなく）メートルで表されます.
この変換は、平面の2次元マップ上に3次元の地球を表示するために使用されます。
データセットの CRS は異なることが多いです。
CRS が同一になるように再定義するのが一般的な手法になります。
レイヤー間の空間的関係を分析するためには、レイヤーが同じ座標参照系を持つことが重要です。
Shapefileには、.prj ファイルに座標参照系に関する情報が格納されています。
Geopandasを使用してデータをGeoDataFrameに読み込むと、
自動的にGeoDataFrameの.crs属性に座標参照系が格納されます。
長々と書きましたが、Geopandasでは投影法の定義は簡単です。
Europe_borders.shpファイルからデータを読み取り、crsを確認してみましょう。
import geopandas as gpd
fp = 'L2_data/Europe_borders.shp'
data = gpd.read_file(fp)
print(data.crs)
{'init': 'epsg:4326'}
CRSの値を確認すると、epsg4326 という値が設定されています。
epsgは「EPSGコード」を指しており、これに関しては後述します。
2.2 EPSG コードとは
EPSGコードとは、GISで使用される様々な要素に必要なパラメータを1つにまとめ、
そのパラメータの集合体どうしを区別するためにコードを割り振ったコード体系のことです。
様々な要素とは、座標参照系、測地系、本初子午線、地図投影法等があります。
EPSG（European Petroleum Survey Group）という団体によって作成されましたが、
現在では様々なGISで使用されています。
詳細は 「EPSGコードとは 」を参照してください。
座標参照系を変更して違いを確認してみましょう。
まずは、変更前の座標参照系で geometry 列を確認します。
print(data['geometry'].head(2))
0 POLYGON ((8.457777976989746 54.56236267089844,...
1 POLYGON ((8.71992015838623 47.69664382934571, ...
Name: geometry, dtype: object
次に、座標参照系を変更して geometry 列を確認します。
# 座標参照系_変更後
orig = data.copy() # 退避
data = data.to_crs(epsg=3035)
print(data['geometry'].head(2))
0 POLYGON ((4221214.558088431 3496203.404338956,...
1 POLYGON ((4224860.478308966 2732279.319617757,...
Name: geometry, dtype: object
座標の値が変更されたことが確認できました。
2.3. 座標参照系の変更前後をそれぞれ投影
座標参照系の変更前後をプロットして、違いを確認してみましょう。
import matplotlib.pyplot as plt
%matplotlib inline
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 8))
orig.plot(ax=ax1, facecolor='gray');
ax1.set_title("WGS84");
data.plot(ax=ax2, facecolor='blue');
ax2.set_title("ETRS Lambert Azimuthal Equal Area projection");
plt.tight_layout()
きちんと変わっていますね。
3. Point in Polygon(PIP)
基本的な地理空間操作として以下2点があります。
- ある点が領域の内側にあるか外側にあるかを調べること
- 線が別の線または多角形と交差するかどうかを調べること
これらを調べるために、Python の shapley というライブラリを使用します。
shapely で PIP を実行するためには、以下2点の方法があります。
- 点が多角形内にあるかどうかを調べる
within（）という関数を使用 - 多角形に点が含まれているかどうかを調べる
contains（）という関数を使用
説明を割愛しますが、LineString または Polygon が別の Polygon の内側にあるかどうかの確認も可能です。
まずは、座標タプルのリストと Point オブジェクトを使って Polygon を作成しましょう。
from shapely.geometry import Point, Polygon
# Point オブジェクトの作成
p1 = Point(24.952242, 60.1696017)
p2 = Point(24.976567, 60.1612500)
# Polygon オブジェクトの作成
coords = [(24.950899, 60.169158), (24.953492, 60.169158), (24.953510, 60.170104), (24.950958, 60.169990)]
poly = Polygon(coords)
# 確認
print(p1)
print(p2)
print(poly)
POINT (24.952242 60.1696017)
POINT (24.976567 60.16125)
POLYGON ((24.950899 60.169158, 24.953492 60.169158, 24.95351 60.170104, 24.950958 60.16999, 24.950899 60.169158))
これらの点がポリゴン内にあるかどうかを確認しましょう
print(p1.within(poly))
print(p2.within(poly))
True
False
最初の点はその多角形の内側にあり、他の点はそうではないことがわかります。
最初の点は、点の位置を多角形の重心と比較するとわかるように、多角形の中心に近いです。
print(p1)
print(poly.centroid)
POINT (24.952242 60.1696017)
POINT (24.95224242849236 60.16960179038188)
ポリゴンに点が含まれているかどうかを確認してみましょう。
print(poly.contains(p1))
print(poly.contains(p2))
True
False
4. 点間の距離計算
Europe_borders.shp を使って、すべてのヨーロッパ諸国の重心（中点）から
フィンランドのヘルシンキまでのユークリッド距離を計算します。
ヘルシンキと他のヨーロッパ諸国との間の距離は,メートル単位の距離を表す
正距方位図法（Azimuthal Equidistant-projection）を使用して計算します。
4.1. GeoDataFrame の作成
まず、フィンランドのヘルシンキの場所を表す単一のPointを含むGeoDataFrameを作成しましょう。
from shapely.geometry import Point
import pycrs
hki_lon = 24.9417
hki_lat = 60.1666
helsinki = gpd.GeoDataFrame([[Point(hki_lon, hki_lat)]], geometry='geometry', crs={'init': 'epsg:4326'}, columns=['geometry'])
4.2. 正距方位図法への変換
この GeoDataFrameを正距方位図法に変換します。
変換を実行するために、pyproj と呼ばれる座標参照系のライブラリを使用します。
proj = 'aeqd'は正距方位図法を指します。
ellps = 'WGS84'は地球楕円体を指します。
World Geodetic System(世界測地系）1984の略語で、
米国が構築・維持している世界測地系です。GPSの軌道情報で使われているほか、
GPSによるナビゲ－ションの位置表示の基準として使われています。
datum = 'WGS84'は測地基準系を指します。
地球上の位置を経度・緯度で表わすための基準を「測地基準系（測地系）」といい、
地球の形に最も近い回転楕円体で定義されています。
lat_0は中心点の緯度座標です。
lon_0は中心点の経度座標です。
import pyproj
aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84', lat_0=hki_lat, lon_0=hki_lon).srs
helsinki = helsinki.to_crs(crs=aeqd)
print(helsinki)
print('\nCRS:\n', helsinki.crs)
geometry
0 POINT (0 0)
CRS:
+units=m +proj=aeqd +ellps=WGS84 +datum=WGS84 +lat_0=60.1666 +lon_0=24.9417
x と y の座標を確認すると、ともに 0 になりました。
ヘルシンキを中心に、ポイントが投影されていることが分かります。
次に、Europe_borders.shp を正距方位図法に変換します。
fp = 'L2_data/Europe_borders.shp'
europe_borders_aeqd = gpd.read_file(fp)
europe_borders_aeqd = europe_borders_aeqd.to_crs(crs=aeqd)
print(europe_borders_aeqd.head(2))
TZID geometry
0 Europe/Berlin POLYGON ((-1057542.597130521 -493724.801682817...
1 Europe/Berlin POLYGON ((-1216418.435359255 -1243831.63520634...
列「geometry」の座標は、ヘルシンキからさまざまな方向へのメートル単位の距離を表しているので、
かなり大きな数であることがわかります。
4.3. ヨーロッパの国境とヘルシンキの位置をプロット
次に、ヨーロッパの国境とヘルシンキの位置をプロットして、座標の再投影が機能したか確認してみましょう。
fig, ax = plt.subplots(figsize=(10,10))
europe_borders_aeqd.plot(ax=ax)
helsinki.plot(ax=ax, color='black', markersize=10)
地図からわかるように、x軸とy軸の0位置はヘルシンキが配置されている場所にあるため、
投影は実際にヘルシンキの中心にあります。
座標値は、ヘルシンキ（黒点）からさまざまな方向（南,北,東,西）までの距離をメートルで示しています。
4.4. すべての国からヘルシンキへの距離を計算
次に、すべての国からヘルシンキへの距離を計算してみましょう。
計算するために、ヨーロッパ諸国の境界を表すすべてのPolygonの重心を計算します。
geopandas では、centroid を使用することで簡単に計算できます。
europe_borders_aeqd['centroid'] = europe_borders_aeqd.centroid
print(europe_borders_aeqd.head(2))
TZID geometry \
0 Europe/Berlin POLYGON ((-1057542.597130521 -493724.801682817...
1 Europe/Berlin POLYGON ((-1216418.435359255 -1243831.63520634...
centroid
0 POINT (-1057718.135423443 -492420.5658204998)
1 POINT (-1218235.216971495 -1242668.589667923)
各Polygonの重心を表すPointジオメトリを持つcentroidという新しい列を作成しました。
次に、国の重心とヘルシンキの間の距離を計算します。
まず,calculate_distance（）という距離を計算するための専用関数を作成します.
def calculate_distance(row, dest_geom, src_col='geometry', target_col='distance'):
"""
Calculates the distance between Point geometries.
Parameters
----------
dest_geom : shapely.Point
A single Shapely Point geometry to which the distances will be calculated to.
src_col : str
A name of the column that has the Shapely Point objects from where the distances will be calculated from.
target_col : str
A name of the target column where the result will be stored.
Returns
-------
Distance in kilometers that will be stored in 'target_col'.
"""
# Calculate the distances
dist = row[src_col].distance(dest_geom)
# Convert into kilometers
dist_km = dist / 1000
# Assign the distance to the original data
row[target_col] = dist_km
return row
ヘルシンキの中心点(dest_geomパラメータに)の座標を取得します。
helsinki_geom = helsinki.loc[0, 'geometry']
print(helsinki_geom)
POINT (0 0)
すべての準備が整ったので、各国の重心とヘルシンキの間の距離を計算してみましょう。
europe_borders_aeqd = europe_borders_aeqd.apply(calculate_distance,
dest_geom=helsinki_geom,
src_col='centroid',
target_col='dist_to_Hki',
axis=1)
print(europe_borders_aeqd.head(10))
TZID geometry \
0 Europe/Berlin POLYGON ((-1057542.597130521 -493724.801682817...
1 Europe/Berlin POLYGON ((-1216418.435359255 -1243831.63520634...
2 Europe/Berlin POLYGON ((-1194521.638643211 -571726.459328881...
3 Europe/Berlin POLYGON ((-1185933.276237431 -571780.052772927...
4 Europe/Berlin POLYGON ((-1182416.219837206 -569097.571220090...
5 Europe/Berlin POLYGON ((-1172799.400650826 -565749.438722840...
6 Europe/Berlin POLYGON ((-1162805.427517967 -563558.434197177...
7 Europe/Berlin POLYGON ((-1129053.540904054 -568388.469960322...
8 Europe/Berlin POLYGON ((-1109126.532709598 -570899.989413469...
9 Europe/Berlin POLYGON ((-703490.1465879158 -664009.791857453...
centroid dist_to_Hki
0 POINT (-1057718.135423443 -492420.5658204998) 1166.724332
1 POINT (-1218235.216971495 -1242668.589667923) 1740.207536
2 POINT (-1194210.78929945 -568987.1532380251) 1322.832487
3 POINT (-1185320.605845876 -571340.3134827728) 1315.832319
4 POINT (-1182191.163363772 -567293.763983083) 1311.258236
5 POINT (-1175758.089721244 -564846.4759828376) 1304.399719
6 POINT (-1157868.191291224 -565162.3607219103) 1288.435968
7 POINT (-1124883.959682355 -567850.4028534386) 1260.086506
8 POINT (-1113376.309723986 -569037.768865205) 1250.364263
9 POINT (-703970.559180466 -663710.5756044324) 967.515517
これで、ヨーロッパ諸国の重心とヘルシンキの間の距離を計算できました。
他のヨーロッパ諸国の重心からヘルシンキまでの最長距離および平均距離を確認してみましょう。
max_dist = europe_borders_aeqd['dist_to_Hki'].max()
mean_dist = europe_borders_aeqd['dist_to_Hki'].mean()
print("Maximum distance to Helsinki is %.0f km, and the mean distance is %.0f km." % (max_dist, mean_dist))
Maximum distance to Helsinki is 3470 km, and the mean distance is 1177 km.
4.5. 補足
今回はヨーロッパ諸国間の距離を計算しましたが、
世界中の複数の場所間の距離を計算したい場合は、Python の harversine パッケージを使いましょう。
環境構築の際にこのパッケージをインストールしているので、興味のある方は試してください。
5. ポイント ⇔ 住所(名称) の変換
OpenStreetMap データに基づくジオコーダーである Nominatim を使って、ポイント ⇔ 住所 の変換をします。
Nominatim は小規模のジオコーディングジョブに使用される場合、サービスの使用にAPIキーを必要としません。
理由は、1秒あたり1リクエスト（3600 /時）に制限されるためです。
大規模のジオコーディングジョブを実行する場合は、Google API を利用するのが良いでしょう。
5.1. 住所(名称) ⇒ ポイント
import geopy.geocoders
from geopy.geocoders import Nominatim
geolocator = Nominatim(format_string="%s, 日本")
location = geolocator.geocode("東京タワー")
print('location.address: ', location.address)
print('location.altitude: ', location.altitude)
print('location.latitude: ', location.latitude)
print('location.longitude: ', location.longitude)
print('location.point: ', location.point)
location.address: 東京タワー, 東京タワー通り, 南青山6, 港区, 東京都, 関東地方, 105-0011, 日本
location.altitude: 0.0
location.latitude: 35.65858645
location.longitude: 139.745440057962
location.point: 35 39m 30.9112s N, 139 44m 43.5842s E
プロキシ環境の場合は、プロキシを設定してから実行しましょう。
proxies= {'http' : 'http://proxysrv:port',
'https': 'http://proxysrv:port'
}
geopy.geocoders.options.default_user_agent = 'mamurata0924'
geopy.geocoders.options.default_timeout = 7
geopy.geocoders.options.default_proxies = proxies
5.2. ポイント ⇒ 住所(名称)
location = geolocator.reverse("35.6806303 139.768196161")
print('location.address: ', location.address)
print('location.altitude: ', location.altitude)
print('location.latitude: ', location.latitude)
print('location.longitude: ', location.longitude)
print('location.point: ', location.point)
location.address: 東京駅一番街, 1, 外堀通り, 中央区, 東京都, 関東地方, 100-6701, 日本
location.altitude: 0.0
location.latitude: 35.68190185
location.longitude: 139.768425758104
location.point: 35 40m 54.8467s N, 139 46m 6.33273s E
5.3. Geocoding in Geopandas
geopandasでは、統合された geopy 機能を使用してジオコーディングを行うことができます。
geopandasにはgeocode()という名前の関数があり、住所のリスト（文字列）をジオコーディングして、
列「geometry」に結果のポイントオブジェクトを含むGeoDataFrameを返すことができます。
import pandas as pd
import geopandas as gpd
from geopandas.tools import geocode
data = [
[1000, 'Itamerenkatu 14, 00101 Helsinki, Finland'],
[1001, 'Kampinkuja 1, 00100 Helsinki, Finland'],
[1002, 'Kaivokatu 8, 00101 Helsinki, Finland'],
[1003, 'Hermannin rantatie 1, 00580 Helsinki, Finland'],
]
columns = ['id', 'addr']
df = pd.DataFrame(data=data, columns=columns)
geo = geocode(data['addr'], provider='nominatim', user_agent='autogis_student_xx')
# プロキシ版
# geo = geocode(data['addr'], provider='nominatim', user_agent='autogis_student_xx', proxies=proxies)
geo.head()
| geometry | address | |
|---|---|---|
| 0 | POINT (24.9155624 60.1632015) | Ruoholahti, 14, Itämerenkatu, Ruoholahti, Läns... |
| 1 | POINT (24.9316914 60.1690222) | Kamppi, 1, Kampinkuja, Kamppi, Eteläinen suurp... |
| 2 | POINT (24.9416849 60.1699637) | Bangkok9, 8, Kaivokatu, Keskusta, Kluuvi, Etel... |
| 3 | POINT (24.9787146 60.1909075) | Hermannin rantatie, Kalasatama, Sörnäinen, Hel... |
プロキシ環境の場合は、プロキシを設定してから実行しましょう。
import urllib.request
proxies= {'http' : 'http://proxysrv:port',
'https': 'http://proxysrv:port'
}
proxy_support = urllib.request.ProxyHandler(proxies)
opener = urllib.request.build_opener(proxy_support)
urllib.request.install_opener(opener)
# プロキシ版
geo = geocode(data['addr'], provider='nominatim', user_agent='autogis_student_xx', proxies=proxies)
geo.head()
6. 最近傍点の探索
GISタスクの1つに、最近傍点の探索が挙げられます。
shapely と geopandas を使って最近傍点の探索を見ていきましょう。
6.1. shapely 編
始点から最も近い終点を見つけることができるようにするには、
終点からMultiPointオブジェクトを作成する必要があります.
from shapely.geometry import Point, MultiPoint
from shapely.ops import nearest_points
orig = Point(1, 1.67)
dest1, dest2, dest3 = Point(0, 1.45), Point(2, 2), Point(0, 2.5)
destinations = MultiPoint([dest1, dest2, dest3])
nearest_geoms = nearest_points(orig, destinations) # 最近傍点
print(nearest_geoms[0])
print(nearest_geoms[1])
POINT (1 1.67)
POINT (0 1.45)
nearest_points()は,最初の項目（インデックス0）が原点のジオメトリで、
2番目の項目（インデックス1）が目的点から実際に最も近いジオメトリのタプルを返します。
つまり、座標(1, 1.67) から最も近い座標は (0, 1.45)となります。
6.2 geopandas 編
GeoDataFrameを使用して、一連の原点から一連の終点に
最も近い点を見つける方法を見ていきましょう。
検索元の GeoDataFrame をdf1、検索先の GeoDataFrame をdf2とします。
data = [
['Suur-Espoonlahti', Point(24.76754037242762, 60.0440879200116)],
['Suur-Kauklahti', Point(24.57415010885406, 60.19764302289445)],
['Pohjoinen', Point(24.94156264995636, 60.24654213027523)],
]
columns = ['Name', 'centroid']
df1 = gpd.GeoDataFrame(data=data, columns=columns)
display(df1)
| Name | centroid | |
|---|---|---|
| 0 | Suur-Espoonlahti | POINT(24.76754037242762 60.0440879200116) |
| 1 | Suur-Kauklahti | POINT(24.57415010885406 60.19764302289445) |
| 2 | Pohjoinen | POINT(24.94156264995636 60.24654213027523) |
data = [
[1000, Point(24.9155624, 60.1632015)],
[1001, Point(24.9316914, 60.1690222)],
[1002, Point(24.9416849, 60.1699637)],
[1003, Point(24.9787146, 60.1909075)],
]
columns = ['id', 'geometry']
df2 = gpd.GeoDataFrame(data=data, columns=columns)
display(df2)
| id | geometry | |
|---|---|---|
| 0 | 1000 | POINT (24.9155624 60.1632015) |
| 1 | 1001 | POINT (24.9316914 60.1690222) |
| 2 | 1002 | POINT (24.9416849 60.1699637) |
| 3 | 1003 | POINT (24.9787146 60.1909075) |
以下の関数を使用して、df2からdf1の重心に最も近いポイント（id列から値を取得）を見つけましょう。
from shapely.geometry import Point, MultiPoint
from shapely.ops import nearest_points
def nearest(row, geom_union, df1, df2, geom1_col='geometry', geom2_col='geometry', src_column=None):
# Find the nearest point and return the corresponding value from specified column.
# Find the geometry that is closest
nearest = df2[geom2_col] == nearest_points(row[geom1_col], geom_union)[1]
# Get the corresponding value from df2 (matching is based on the geometry)
value = df2[nearest][src_column].get_values()[0]
return value
unary_union = df2.unary_union # MultiPoint オブジェクトの作成
df1['nearest_id'] = df1.apply(nearest, geom_union=unary_union, df1=df1, df2=df2, geom1_col='centroid', src_column='id', axis=1)
display(df1.head())
| Name | centroid | nearest_id | |
|---|---|---|---|
| 0 | Suur-Espoonlahti | POINT(24.76754037242762 60.0440879200116) | 1000 |
| 1 | Suur-Kauklahti | POINT(24.57415010885406 60.19764302289445) | 1000 |
| 2 | Pohjoinen | POINT(24.94156264995636 60.24654213027523) | 1003 |
7. オーバーレイ解析
オーバーレイ解析とは、地理的な要素を重ね合わせることにより、
情報同士の関係性や新たな発見を導き出す解析を指します。
詳細は ESRI ジャパンのオーバーレイ解析のページを参照してください。
基本的なオーバーレイ操作
画像は Automating GIS-processes 2018 より引用
7.1. データ準備
L4_data.zip をダウンロードします。
zip を展開すると以下のような構成になっています。
L4_data
├ Amazon_river.cpg
├ Amazon_river.dbf
├ Amazon_river.prj
├ Amazon_river.shp
├ Amazon_river.shx
├ Helsinki_borders.cpg
├ Helsinki_borders.dbf
├ Helsinki_borders.prj
├ Helsinki_borders.shp
├ Helsinki_borders.shx
├ TravelTimes_to_5975375_RailwayStation.cpg
├ TravelTimes_to_5975375_RailwayStation.dbf
├ TravelTimes_to_5975375_RailwayStation.prj
├ TravelTimes_to_5975375_RailwayStation.shp
└ TravelTimes_to_5975375_RailwayStation.shx
7.2. データ読込
以下のファイルを読み込んで、可視化してみましょう。
TravelTimes_to_5975375_RailwayStation_Helsinki.shpHelsinki_borders.shp
import geopandas as gpd
import matplotlib.pyplot as plt
import shapely.speedups
%matplotlib inline
# shape ファイルの読込が速くなるおまじない
shapely.speedups.enable()
border_fp = './lesson4/Helsinki_borders.shp'
grid_fp = './lesson4/TravelTimes_to_5975375_RailwayStation.shp'
grid = gpd.read_file(grid_fp)
hel = gpd.read_file(border_fp)
ax = grid.plot(facecolor='gray')
hel.plot(ax=ax, facecolor='None', edgecolor='blue')
オーバーレイ解析をして、ヘルシンキレイヤ(青い領域)と交差する
グリッドのポリゴンから新しいレイヤを作成しましょう。
7.3. CRS の一致確認
解析の前に、CRSの確認をしましょう。
複数のレイヤーを含む GIS 操作を実行する場合は、レイヤーの CRS が一致している必要があります。
hel.crs == grid.crs
True
一致していますね。
7.4. オーバーレイ解析
geopandas を使えば、簡単にオーバーレイ解析が可能です。
ここでは、intersection の操作をしてみましょう。
intersection = gpd.overlay(grid, hel, how='intersection')
intersection.plot(color="b")
ヘルシンキの境界線と交差するグリッドセルのみが表示されるようになりました。
グリッドセルは境界に基づいて切り取られています。
データセットには両方の入力レイヤからの属性が含まれています。
display(intersection.head())
| car_m_d | car_m_t | car_r_d | car_r_t | from_id | pt_m_d | pt_m_t | pt_m_tt | pt_r_d | pt_r_t | pt_r_tt | to_id | walk_d | walk_t | GML_ID | NAMEFIN | NAMESWE | NATCODE | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 94 | 29476 | 41 | 29483 | 46 | 5876274 | 29990 | 76 | 95 | 24984 | 77 | 99 | 5975375 | 25532 | 365 | 27517366 | Helsinki | Helsingfors | 091 | POLYGON ((402250.0001322526 6685750.000039577,... |
| 95 | 29456 | 41 | 29462 | 46 | 5876275 | 29866 | 74 | 95 | 24860 | 75 | 93 | 5975375 | 25408 | 363 | 27517366 | Helsinki | Helsingfors | 091 | POLYGON ((402367.8898583531 6685750.000039573,... |
| 96 | 36772 | 50 | 36778 | 56 | 5876278 | 33541 | 116 | 137 | 44265 | 130 | 146 | 5975375 | 31110 | 444 | 27517366 | Helsinki | Helsingfors | 091 | POLYGON ((403250.000132058 6685750.000039549, ... |
| 97 | 36898 | 49 | 36904 | 56 | 5876279 | 33720 | 119 | 141 | 44444 | 132 | 155 | 5975375 | 31289 | 447 | 27517366 | Helsinki | Helsingfors | 091 | POLYGON ((403456.4840652301 6685750.000039542,... |
| 193 | 29411 | 40 | 29418 | 44 | 5878128 | 29944 | 75 | 95 | 24938 | 76 | 99 | 5975375 | 25486 | 364 | 27517366 | Helsinki | Helsingfors | 091 | POLYGON ((402000.0001323065 6685500.000039622,... |
7.5. 解析結果を GeoJSON に出力
shape ファイルに出力でも良いのですが、
今回は試しに、GeoJSON に出力します。
※注意
GeoPandas の仕様上、geometry の 列タイプが 'MultiPolygon' の必要あり
2019/2/18時点で、issue 867 で対応中
from shapely.geometry.multipolygon import MultiPolygon
intersection['geometry'] = intersection['geometry'].apply(lambda x: MultiPolygon([x]))
outfp = "./L4_data/TravelTimes_to_5975375_RailwayStation_Helsinki.geojson"
intersection.to_file(outfp, driver="GeoJSON")
8. インタラクティブマップの作成
インタラクティブマップの作成には、Python ライブラリの folium を使用します。
マップには、国土交通省が作成した「地形分類図」「土壌図」を表示していきます。
また、マップ作成には folium の Gallery を参考にしました。
Jupyter Notebook 形式で非常に分かりやすいので、これから folium に触れる方は必見です。
8.1. データ準備
国土調査 土地分類基本調査 50万分の１土地分類基本調査から、
「地形分類図」「土壌図」をダウンロードします。
landform
├ landform.avl
├ landform.dbf
├ landform.shp
└ landform.shx
soil
├ soil.avl
├ soil.dbf
├ soil.shp
└ soil.shx
8.2. データ前処理
ダウンロードしたデータを読み込み、ラベルエンコーディングします。
今回は東京都を表示対象にします。
from sklearn.preprocessing import LabelEncoder
import pandas as pd
import geopandas as gpd
import folium
import branca
def label_encode(df, column_name):
le = LabelEncoder()
le = le.fit(df[column_name])
return le.transform(df[column_name]).astype(int)
# 地形分類図の読み込み
gdf_land = gpd.read_file('./landform/landform.shp', encoding='shift-jis')
gdf_land = gdf_land[gdf_land['PREF']=='東京都']
gdf_land['terrain'] = label_encode(gdf_land, '地形大区分')
display(gdf_land.head(3))
| ID | NO | PREF | 地形区分名 | 地形大区分 | geometry | terrain | |
|---|---|---|---|---|---|---|---|
| 10605 | 0 | c | 東京都 | 埋立地 | 埋立地 | POLYGON ((139.8134010579003 35.59388260180398,... | 3 |
| 10606 | 0 | c | 東京都 | 埋立地 | 埋立地 | POLYGON ((139.8051275277766 35.59496175790707,... | 3 |
| 10607 | 0 | c | 東京都 | 埋立地 | 埋立地 | POLYGON ((139.7799472187043 35.61276783360819,... | 3 |
# 土壌図の読み込み
gdf_soil = gpd.read_file('./soil/soil.shp', encoding='shift-jis')
gdf_soil = gdf_soil[gdf_soil['PREF']=='東京都']
gdf_soil['soil_type'] = label_encode(gdf_soil, '土壌区分名')
display(gdf_soil.head(3))
| ID | NO | PREF | 土壌区分名 | SOILNAME11 | 土壌大区分 | SOILNAME12 | SOILNAME22 | SOILNAME21 | geometry | soil_type | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 5429 | 0 | c | 東京都 | 埋立地 | c | 埋立地 | c | c | c | POLYGON ((139.8134010579003 35.59388260180398,... | 3 |
| 5430 | 0 | c | 東京都 | 埋立地 | c | 埋立地 | c | c | c | POLYGON ((139.8051275277766 35.59496175790707,... | 3 |
| 5431 | 0 | c | 東京都 | 埋立地 | c | 埋立地 | c | c | c | POLYGON ((139.7799472187043 35.61276783360819,... | 3 |
8.3. カラーマップ作成
branca というライブラリを使ってカラーマップを作成します。
このライブラリは folium から派生したもので、マップ固有ではない機能を提供します。
詳しくはこちらのページを参照してください。
まずは、地形分類図における地形区分のカラーマップを作成します。
import branca.colormap as cm
linear_land_max = gdf_land['terrain'].max()
linear_land_min = gdf_land['terrain'].min()
linear_land_length = len(gdf_land['terrain'].unique())
linear_land = cm.linear.YlGnBu_09.scale(linear_land_min, linear_land_max).to_step(linear_land_length)
linear_land.caption = "地形区分"
display(linear_land)
次に、土壌図における土壌区分のカラーマップを作成します。
linear_soil_max = gdf_soil['soil_type'].max()
linear_soil_min = gdf_soil['soil_type'].min()
linear_soil_length = len(gdf_soil['soil_type'].unique())
linear_soil = cm.linear.YlOrBr_07.scale(linear_soil_min, linear_soil_max).to_step(linear_soil_length)
linear_soil.caption = "土壌区分"
display(linear_soil)
8.4. EPSGコードの設定
デフォルトだと crs が None になっているため、EPSG コードを設定します。
設定値は ここ を参考にしました。
from fiona.crs import from_epsg
gdf_land.crs = from_epsg(2450)
gdf_soil.crs = from_epsg(2450)
8.5. コロプレス図(地形分類図編)の作成
コロプレス図とは、統計数値を地図に表現する場合の一方法として、数値をいくつかの階級に区分し、
区域単位、たとえば市区町村別に、各階級に応じた色彩または明暗によって表す地図を指します。
ここでは、地形分類図のデータをコロプレス図として表現します。
folium には コロプレス図を作成する choropleth というメソッドがありますが、
このメソッドはマップ作成の自由度が低いです。
そのため、GeoJSONというメソッドを使用して、コロプレス図を作成します。
# ベースマップの作成
tokyo23_location = [35.658593, 139.745441]
m = folium.Map(location=tokyo23_location,tiles='cartodbpositron',zoom_start=10)
# コロプレス図の追加
folium.GeoJson(
data=gdf_land.to_json(),
style_function=lambda feature: {
'fillColor': linear_land(feature['properties']['terrain']),
'color': 'gray',
'weight': 2,
'fill_opacity': 0.6,
'line_opacity': 0.2,
},
highlight_function=lambda feature: {
'fillColor': linear_land(feature['properties']['terrain']),
'color': 'gray',
'weight': 4,
'dashArray': '1, 1',
'fill_opacity': 0.6,
'line_opacity': 0.2,
},
tooltip=folium.features.GeoJsonTooltip(
fields=['terrain', '地形区分名'],
aliases=['地形区分コード', '地形区分名'],
localize=True
),
name="【図】地形分類"
).add_to(m)
# 凡例の追加
linear_land.add_to(m)
# マップ表示
display(m)
コードの解説をします。
最初にベースマップについて説明します。
tiles には cartodbpositron を指定していますが、他の値でも設定可能です。
以下の記事を見ながら、タイルを切り替えてマップを作成するとよいでしょう。
Python folium 指定できる 地図の タイル について
次に、folium.GeoJson 内の引数について説明します。
data=gdf_land.to_json() の gdf_land.to_json() は、data に GeoJSON を設定するという意味です。
geopandas の to_json()関数は、GeoDataFrame を GeoJSON に変換します。
{
"type": "FeatureCollection",
"features": [{
"id": "10605",
"type": "Feature",
"properties": {
"ID": 0,
"NO": "c",
"PREF": "東京都",
"地形区分名": "埋立地",
"地形大区分": "埋立地",
"terrain": 3
},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[138.5015508233247, 36.00032080614344],
[138.50155081335006, 36.00032082883919],
[138.50155088118132, 36.000320862881736],
[138.50155093704217, 36.00032088395551],
[138.5015509729516, 36.000320838563816],
[138.50155099090588, 36.00032079479348],
[138.50155098691576, 36.00032078993018],
[138.5015509769406, 36.00032078506696],
[138.50155092107954, 36.00032075588761],
[138.50155089115398, 36.000320739676845],
[138.50155087918387, 36.000320739676994],
[138.50155086920887, 36.00032074291936],
[138.50155085324923, 36.00032076885741],
[138.5015508233247, 36.00032080614344]
]
]
},
"bbox": [138.50155081335006,
36.000320739676845,
138.50155099090588,
36.00032088395551
]
}]
}
style_function と highlight_function の feature['properties']['terrain'] の部分は
GeoJSON上の feature['properties']['terrain']と紐づいています。
この部分を理解すれば、コロプレス図でのハイライト・スタイル設定は容易になります。
上記以外の箇所については、公式ドキュメントを確認してください。
次に表示されたマップの内容を見ていきましょう。
地形区分に対応した色でエリアが表現されていますね。
右上にはカラーマップの凡例があり、フォーカスを合わせると、
ハイライトされた上でツールチップが表示されます。
8.6. コロプレス図(土壌図編)の作成
同様に土壌図のデータもコロプレス図にしてみましょう。
# コロプレス図の追加
folium.GeoJson(
data=gdf_soil.to_json(),
style_function=lambda feature: {
'fillColor': linear_soil(feature['properties']['soil_type']),
'color': 'gray',
'weight': 2,
'fill_opacity': 0.6,
'line_opacity': 0.2,
},
highlight_function=lambda feature: {
'fillColor': linear_land(feature['properties']['soil_type']),
'color': 'gray',
'weight': 4,
'fill_opacity': 0.6,
'line_opacity': 0.2,
},
tooltip=folium.features.GeoJsonTooltip(
fields=['soil_type', '土壌区分名'],
aliases=['土壌区分コード', '土壌区分名'],
localize=True
),
name="【図】土壌"
).add_to(m)
# 凡例の追加
linear_soil.add_to(m)
# マップ表示
display(m)
土壌図のデータをコロプレス図にできました。
地形分類図のレイヤーと重なっているため、非常に見にくいという問題があります。
この問題を解消するために、次にレイヤコントロールの説明をします。
8.7. レイヤーコントロールの追加
レイヤコントロールの追加は非常に簡単です。
たった1行で済みます。
folium.LayerControl(collapsed=False).add_to(m)
display(m) # マップ表示
マップにチェックボックスが追加されました。
これで図の種類を切り替え可能です。
8.8. マーカーの追加
地図上にマーカーを追加します。
マーカーのアイコン設定にはFont Awesome Iconsを参照してください。
import folium.plugins
# ベースマップ
tokyo23_location = [35.658593, 139.745441]
m = folium.Map(location=tokyo23_location,tiles='cartodbpositron',zoom_start=10)
# マーカークラスター
mcg = folium.plugins.MarkerCluster(control=False)
m.add_child(mcg)
# グループ分け
g1 = folium.plugins.FeatureGroupSubGroup(mcg, '【敷地】訓練用')
m.add_child(g1)
g2 = folium.plugins.FeatureGroupSubGroup(mcg, '【敷地】テスト用')
m.add_child(g2)
# マーカー作成
folium.Marker(
location=[35.7, 139.80],
tooltip='訓練用_1',
icon=folium.Icon(color='green', icon='heart')
).add_to(g1)
folium.Marker(
location=[35.75, 139.85],
tooltip='訓練用_2',
icon=folium.Icon(color='green', icon='heart')
).add_to(g1)
folium.Marker(
location=[35.80, 139.90],
tooltip='テスト用_1',
icon=folium.Icon(color='blue', icon='info-sign')
).add_to(g2)
folium.Marker(
location=[35.85, 139.95],
tooltip='テスト用_2',
icon=folium.Icon(color='blue', icon='info-sign')
).add_to(g2)
# レイヤーコントロールの追加
folium.LayerControl(collapsed=False).add_to(m)
# マップ表示
display(m)
マーカークラスターにグループを2つ登録したので、
地図上でもそのようになっていますね。
地図を拡大すると、マーカーが4つ打たれていることが分かります。
画像では表示していませんが、マーカーにフォーカスを合わせるとツールチップが表示されます。
また、右上にあるグループ切り替えのチェックボックスで、マーカーの表示を切り替えられます。
8.9. タイルの追加
今までのタイルは cartodbpositron でした。
google maps と google street view を追加してみましょう。
tokyo23_location = [35.658593, 139.745441]
m = folium.Map(location=tokyo23_location,tiles='cartodbpositron',zoom_start=10)
folium.raster_layers.TileLayer(
tiles='http://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
attr='google',
name='google maps',
max_zoom=20,
subdomains=['mt0', 'mt1', 'mt2', 'mt3'],
overlay=False,
control=True,
).add_to(m)
folium.raster_layers.TileLayer(
tiles='http://{s}.google.com/vt/lyrs=m&x={x}&y={y}&z={z}',
attr='google',
name='google street view',
max_zoom=20,
subdomains=['mt0', 'mt1', 'mt2', 'mt3'],
overlay=False,
control=True,
).add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
display(m)
タイルをラジオボタンで切り替えられるようになりました。
8.10 ミニマップの追加
「東京を見ていたけど、次は北海道を見たい」というときに、
マップ上でスライドして移動するのは非常に面倒です。
そんな場合はミニマップを活用しましょう。
tokyo23_location = [35.658593, 139.745441]
m = folium.Map(location=tokyo23_location,tiles='cartodbpositron',zoom_start=10)
minimap = folium.plugins.MiniMap(toggle_display=True, width=200, height=200)
m.add_child(minimap)
display(m)
このようになります。右下の矢印をクリックするとミニマップの表示切替が可能です。
8.11. マウスポジションの追加
マップ上でマウスオーバーしているときに、その箇所の座標を知りたいときがあります。(きっとあるはずです・・・)
そういったときは、マウスポジションを追加しましょう。
tokyo23_location = [35.658593, 139.745441]
m = folium.Map(location=tokyo23_location,tiles='cartodbpositron',zoom_start=10)
formatter = "function(num) {return L.Util.formatNum(num, 3) + ' º ';};"
folium.plugins.MousePosition(
position='topright',
separator=' | ',
empty_string='NaN',
lng_first=True,
num_digits=20,
prefix='Coordinates:',
lat_formatter=formatter,
lng_formatter=formatter
).add_to(m)
display(m)
マップの右上に座標が表示されるようになりました。
8.12. HTMLに出力
マップに可視化した結果を共有したい場合があります。
folium では HTML に出力して保存が可能です。
m.save('./test.html')
まとめ
ヘルシンキ大学講座には、この記事に記載していない、役立つ内容が盛りだくさんです。
Python で GIS データを活用したい人は、ぜひ受講して下さい。