【Python】地理空間データを扱えるCartopyのインストール方法【Cartopy】

Pythonで地理空間データ (geography information system data, GIS data または geospatial data)（位置が緯度，経度で表された地球上のある値に関するデータ）を扱いたいと思ったことはありませんか？たとえば，気象データ（各地点のアメダスデータ，数値予報が出力する数値モデル上の格子点データ（GPVデータ））はGISデータの代表的なものの一つです。

別の記事でも folium というライブラリを使用してGISデータをプロットする方法を紹介しています。
今回はCartopyというGISデータを扱うためのパッケージとその基本的な使い方を紹介します。

Cartopyを使うと，例えば異なる図法（メルカトル図法，正角方位図法，ランベルト正角円錐図法…）でGISデータを非常に簡単にプロットすることができます。

動作検証済み環境

macOS BigSur (11.7.1), python3.9.15, Cartopy0.21.0

Pythonでやってみる気象データ解析の技術書を販売中

¥2,500 → ¥1,500 今なら40%OFF!!

気象の勉強をされている方必読！
気象予報士でもある著者がわかりやすく解説
プログラムを公開しているのでお手元で試せます！

技術書ページへ

Cartopyとは？

Cartopyとは，GISデータをPythonで扱うためのパッケージです。もともとはイギリスの気象庁 (UK MetOffice) で開発が始まったそうです。
地球は球体（正確には回転楕円体というべきか）ですから，通常私達が用いるデカルト座標系は局所的にしか満たされません。
地球上のデータを二次元で表すためには，様々な図法が用いられおり，データを然るべき方法で変換がする必要があります（例えば，メルカトル図法で極地方で大陸が引き伸ばされていたことを思い出してください）。

CartopyはGISデータの面倒な変換を行ってくれます。それでは，早速インストールしてみましょう！

Cartopyのインストール方法

Cartopyはpipやcondaでインストールできます。今回は，公式サイトで説明されている方法のうち，homebrewとpipを使用する方法を説明します。

Cartopy本体をインストールする前に，Cartopyが依存している（内部で呼び出している）GEOS，Shapely，pyshpをインストールする必要があります。これらの間の依存関係やリンクはややこしいため，説明されている順序で行う必要があります。

なにはともあれ，まずは，homebrewをアップデートしておきましょう（ $マークは「プロンプト」です。入力する必要はありません。）。

$ brew update

しばらくアップデートしていないと，時間がかかる場合があります。アップデートが終わったら，GEOSをインストールします。

$ brew install geos

インストールが終わったら，pipでpyshpをインストールします。

$ pip3 install --upgrade pyshp

次に，Shapely ですが，ソースコードからビルドする必要があります。これは，上でインストールしたGEOSとリンクさせるためです。したがって，すでにShapelyをインストールしている方は，アンインストールする必要があります。インストールされているかいないかは，pip3 list とすると出てくるインストール済みパッケージリストで確認することができます。インストールされている場合は，コメントアウトしているpip3 uninstall shapelyを実行してアンインストールしてください。

$ pip3 list
$ #pip3 uninstall shapely
$ pip3 install "shapely<2" --no-binary shapely

Cartopyのビルドの際に必要になるpkg-configをインストールします。インストール後，PKG_CONFIG_PATHを設定しておきます。

$ brew install pkg-config
$ export PKG_CONFIG_PATH=/usr/local/bin/pkgconfig

ここまでできれば，pipでCartopyが入れられます。

$ pip3 install cartopy

以上でインストールは完了です！！
現在 (2022/11/30) この方法でインストールすると，Cartopy0.21.0が入ります。

日本周辺の地図のプロット

正距円筒図法でプロットするプログラム

公式サイトのギャラリーを参考に日本周辺の地図をプロットし，通天閣があるところに点を打ってみましょう！

以下のようなコードを書いて，map-japan.py という名前で，プログラムと生成される図を保存するディレクトリに（ここでは~/Desktop/LabCode/python/weather とします）に保存しましょう。

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(294/25.4, 210/25.4))
ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.PlateCarree())
# ---
# 別の図法で描きたい場合は例えば以下のようにする。
# たくさんの図法が用意されています。
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.Mercator())
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.Robinson())
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.LambertCylindrical())
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.LambertConformal())

ax.set_extent([120, 160, 20, 50], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.RIVERS)

# --- 通天閣の位置
ax.scatter(135+30/60+22.77/60/60, 34.0+39/60+9.13/60/60, transform=ccrs.PlateCarree())

# --- 緯線と経線
ax.gridlines(linestyle='dashed', color='gray', linewidth=0.5,  draw_labels=True)

plt.savefig('map-japan.png')

プログラムを実行する

ターミナルを開き，

$ cd Desktop/LabCode/python/weather/

と入力し，先程プログラムを保存したディレクトリを移動します。あとは以下のコマンドでmap-japan.pyを実行するだけです。（ $マークは「プロンプト」です。入力する必要はありません。）

$ python3 map-japan.py

実行結果

実行すると~/Desktop/LabCode/python/weather にmap-japan.png という画像ファイルが生成され，次のようになっていれば成功です。

初めてCartopyを使って地図をプロットするとき（初めてその解像度の地図を使うときも同様），ダウンロードが発生するため，時間がかかる場合があります。また，ターミナル上にwarningがでますが，心配いりません。

また，通天閣の位置に，青い点がプロットされている事がわかります！

コードの解説

上に書いたソースコードの解説をしていきます。

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

まず，必要なモジュールをインポートします。

fig = plt.figure(figsize=(294/25.4, 210/25.4))
ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.PlateCarree())
# ---
# 別の図法で描きたい場合は例えば以下のようにする。
# たくさんの図法が用意されています。
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.Mercator())
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.Robinson())
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.LambertCylindrical())
# ax = fig.add_axes([0.10, 0.10, 0.80, 0.80], projection=ccrs.LambertConformal(central_longitude=140))

用紙設定とプロット領域の設定をします。ここでは，A4サイズ横置き (縦210 mm，横294 mm）を想定します。

プロット領域は，用紙に対して左から10%，下から10%をプロット枠の左下隅として，プロット枠大きさを用紙の横80%，縦80%で指定します。正距円筒図法を使用するので，projection=ccrs.PlateCarree()をオプションにつけます。

別の図法で描きたい場合は，引数のprojectionにわたすものを変えてください。例えば，コメントアウトしたものが使用できます。ただし，ここに書いたものはあくまで一例です。他にもたくさん用意されています。

また，一番下のccrs.LambertConformal(central_longitude=140) のように，引数を入力できるようになっていて，デフォルト値を使うと思った図が出ない場合があります。

ax.set_extent([120, 160, 20, 50], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.RIVERS)

描画領域を東経120度から東経160度，北緯20度から北緯50度に設定します。crsの引数は，どの図法を使うときもccrs.PlateCarree()を指定してください。

地球上の要素である陸地 (land)，海洋 (ocean)，海岸線 (coastline)，国境線 (borders)，湖沼 (lakes)，河川 (rivers) を追加していきます。ちなみに最初に実行するときは，これらの要素をダウンロードしてくるため，ターミナル上にwarningが表示されます。

# --- 通天閣の位置
ax.scatter(135+30/60+22.77/60/60, 34.0+39/60+9.13/60/60, transform=ccrs.PlateCarree())

# --- 緯線と経線
ax.gridlines(linestyle='dashed', color='gray', linewidth=0.5, draw_labels=True)

plt.savefig('map-japan.png')

最後に，通天閣の位置（北緯34度39分9.13秒，東経135度30分22.77秒）をプロットします。通常のscatterプロットとは違って，transform=ccrs.PlateCarree() をつけることを忘れないでください。どの図法を利用しても，transform=ccrs.PlateCarree()でいいようです。

次の行のax.gridlines()では，緯線と経線を描くようにしています。

最後の行は，描いた図を map-japan.png という名前（とpngという形式）で保存するようにしています。