【Cartpy入門】Pythonで地図上にグラフを描画する

# 概要
気象や地震などのデータ解析をしていると，地図上にプロットしたりコンターを描いたりしたくなります．

そんなときに使えるのが，Pythonのライブラリ**Cartpy**です．

本記事では，Cartpyのドキュメント↓を参考にしながら基本的な使い方を解説します．

https://scitools.org.uk/cartopy/docs/latest/index.html

# とりあえず地図だけ描画してみる

## 基本形
地図（海岸線）だけ表示させるには，下のように書けばOKです．

```python
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

ax = plt.axes(projection=ccrs.PlateCarree())  # PlateCarreeは正距円筒図法
ax.coastlines()  # 海岸線を描画する
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/e7603314-40e4-48dc-e073-63e736138410.png)

他の地図投影法を使いたい場合は，`projection=`の部分を色々変えてみましょう．

- 正距円筒図法：`ccrs.PlateCarree()`
- メルカトル図法：`ccrs.Mercator()`
- 極投影図法：`ccrs.NorthPolarStereo()`
- ランベルト正角円錐図法：`ccrs.LambertConformal()`
- 正射投影図法：`ccrs.Orthographic()`
- ロビンソン図法：`ccrs.Robinson()`
- モルワイデ図法：`ccrs.Mollweide()`
- ランベルト正積円筒図法：`ccrs.LambertCylindrical()`
- ミラー図法：`ccrs.Miller()`

## 範囲指定＆経度線・緯度線の描画

全世界ではなく，特定の範囲だけを見たい場合は`ax.set_extent()`で範囲を指定します．

の引数の形式は

`[経度1, 経度2, 緯度1, 緯度2]`

です．

また，経度線・緯度線を描画したい場合は`ax.gridlines()`を使います．引数を指定すればいろいろカスタマイズできます．

`gridlines()`の詳細はこちら↓

https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html#cartopy.mpl.geoaxes.GeoAxes.gridlines

```python
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([120, 150, 20, 50], crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(
    draw_labels=True,  # 軸ラベルの設定
    xlocs=np.arange(120, 150, 5),  # x軸の目盛りの設定
    ylocs=np.arange(20, 50, 2.5),  # y軸の目盛りの設定
    dms=True,  # 度分秒表示の設定(True:度分秒表示, False:10進表示)
    )
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/8bedade9-8b56-557a-a445-a84c50ea0bc2.png)

## 国境線や河川を描く

`cartopy.feature`をインポートすると，国境線や河川を描けます．

```python
import cartopy.feature as cfeature

ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([120, 150, 20, 50], crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines()

ax.add_feature(cfeature.BORDERS, linestyle=':') # 国境線を描く
ax.add_feature(cfeature.RIVERS, edgecolor='red') # 河川を描く

plt.tight_layout()
plt.savefig('coastlines_gridlines_feature.png')
plt.show()
```

# 地図の上に点をプロット

基本的には，`ax.plot()`とか`ax.scatter()`を使えばOKです．

```python
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([120, 150, 20, 50], crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines()

# 緯度35度、経度135度の位置に赤い点を描く
ax.plot(135, 35, marker='o', color='red')

plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/f077459e-0348-3ede-589d-307df684614a.png)

正距円筒図法（`projection=ccrs.PlateCarree()`）であればこれでちゃんと描画できますが，他の図法を使うとおかしくなります．

例えば`projection`を`Metrcator()`に設定します．

```python
ax = plt.axes(projection=ccrs.Mercator()) # メルカトル図法
ax.set_extent([120, 150, 20, 50], crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines()

# 緯度35度、経度135度の位置に赤い点を描く
ax.plot(135, 35, marker='o', color='red')  # transformを指定しないと、projectionで指定した図法での座標になる

plt.show()
```

すると...

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/e266b958-b735-0f0a-23ba-7afa79f1e55b.png)

赤点がどこかに行ってしまいました．

これは，`135, 5`という数値がどの座標系での値なのかを指定してあげれば解決します．

つまり，グラフ描画部分の引数に`transform`を追加します．

```python
ax.plot(135, 35, marker='o', color='red', transform=ccrs.PlateCarree())
```

これはちょっと注意が必要ですね．
ドキュメントでもわかりやすく解説されています．

https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html

# ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．

本記事のコードは`05-cartpy`の中に入っています．

https://github.com/matsuda-tkm/qiita-calender-2023