PythonとCartopyで地図データを可視化する

PythonでGIS用シェープファイルを可視化するでは、日本地図とフェリー航路のシェープファイルをPythonライブラリのGeoPandasで可視化した。日本地図のシェープファイルは国土地理院のホームページからダウンロードしたが、地理空間データ生成用のPythonパッケージであるCartopyを使うと、自分で地図のシェープファイルを用意する必要がない。CartopyにはNatural Earthのシェープファイルをロードするメソッドが用意されている。今回はCartopyを使って地図データを可視化してみる。

環境

WSL2(Ubuntu20.04)とJupyter。

$ lsb_release -dr
Description:    Ubuntu 20.04.4 LTS
Release:        20.04
$ python3 -V
Python 3.8.10

Pythonライブラリのインストール

必要なPythonライブラリをインストールする。Cartopyのインストールは少し手間がかかるので後述。

$ pip3 install matplotlib geopandas

インストールされたバージョン。

$ pip show matplotlib
Name: matplotlib
Version: 3.5.1
...
$ pip show geopandas
Name: geopandas
Version: 0.10.2
...

Cartopyのインストール

Cartopyはpipでインストールというわけにはいかず、少し手間がかかる。Cartopyをインストールする前に、proj、geos、pyshp、shapelyをインストールしておく必要がある。

まずはprojのインストール。projはaptコマンドでもインストールできるが、バージョンが古いのでソースからインストールする。ビルドの前に必要なパッケージをインストールする。

$ sudo apt install cmake sqlite3 libsqlite3-dev 

続いてソースファイルをダウンロードしてビルドする。今回はバージョン8.2.1を使う。

$ curl -O https://download.osgeo.org/proj/proj-8.2.1.tar.gz
$ tar xzvf proj-8.2.1.tar.gz
$ cd proj-8.2.1
$ mkdir build
$ cd build
$ cmake ..
$ cmake --build .
$ sudo cmake --build . --target install

次にgoesのインストール。

$ sudo apt install libgeos-dev
$ pip3 install geos

pyshpとshapelyのインストール。

$ pip3 install pyshp shapely

インストールした各ライブラリのバージョン。

$ pip3 show geos
Name: geos
Version: 0.2.3
...
$ pip3 show pyshp
Name: pyshp
Version: 2.2.0
...
$ pip3 show shapely
Name: Shapely
Version: 1.8.1.post1
...

最後にCartopyのインストール。

$ pip3 install cartopy

インストールされたバージョン

$ pip3 show cartopy
Name: Cartopy
Version: 0.20.2
...

CartopyとGeoPandasで地図とフェリー航路を可視化する

地図の他に、PythonでGIS用シェープファイルを可視化すると同様にフェリー航路も可視化する。フェリー航路のシェープファイルは、国土地理院のホームページにある地球地図日本の第2.2版ベクタの全レイヤを使う。ファイルをダウンロードして解凍しておく。フェリー航路のシェープファイルはferryl_jpn.shp。

地図とフェリー航路を可視化するコードは以下の通り。

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

# 解凍したフェリー航路のシェープファイルのパス
path_ferry =  "gm-jpn-all_u_2_2/ferryl_jpn.shp"
df_ferry = gpd.read_file(path_ferry)

plt.figure(figsize=(20, 16))

# 地図投影法を指定
ax = plt.axes(projection=ccrs.PlateCarree())

# Natural Earthの地図シェープファイル（解像度10m）を読み込む
ax.coastlines(resolution='10m')

# フェリー航路をプロット
df_ferry.plot(ax=ax, color='blue', linewidth=1.0, label='Ferry routes')

以下のように地図とフェリー航路が表示される。クリッピングされる地図の位置は、読みこんだシェープファイル（今回はフェリー航路）の位置により自動調整される。

北海道が途中で切れているので北海道全体がクリッピング範囲に入るようにしてみる。Cartopyにはset_extentというメソッドが用意されていて、任意の場所でクリッピングできるらしいのだが、実行しても処理が途中で落ちてしまう。ここでは回避策として、透明なポリゴンの境界を配置して、クリッピング位置が自動調整されたときに北海道が切り取られないようにする。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import geopandas as gpd
from shapely.geometry import Polygon

# 解凍したフェリー航路のシェープファイルのパス
path_ferry =  "gm-jpn-all_u_2_2/ferryl_jpn.shp"
df_ferry = gpd.read_file(path_ferry)

plt.figure(figsize=(20, 16))

# 地図投影法を指定
ax = plt.axes(projection=ccrs.PlateCarree())

# Natural Earthの地図シェープファイル（解像度10m）を読み込む
ax.coastlines(resolution='10m')

# フェリー航路をプロット
df_ferry.plot(ax=ax, color='blue', linewidth=1.0, label='Ferry routes')

# 指定範囲で地図を切り取る
# set_extentがエラーになるので透明なポリゴンの境界で対処
# ax.set_extent([130, 150, 30, 40], crs=ccrs.PlateCarree())
polygon = Polygon([(130, 40), (150, 40), (150, 30), (130, 30), (130, 40)])
poly_gdf = gpd.GeoDataFrame([1], geometry=[polygon], crs=df_ferry.crs)
poly_gdf.boundary.plot(ax=ax, alpha=0)

以下のように北海道が途中で切れることなくプロットされた。