Python: Cartopy で地図を書いてみる

# はじめに
これまで地図関係の作図は、GMT (generic mapping Tools) を用いて行ってきたが、そろそろ Python で試してみようと思い、Cartopy をインストールし、使用してみることにした。結構いい感じの出力ができる模様。

#### 環境
+ M1 MacBook Air
+ macOS Big Sur (Version 11.5.2).
+ python version 3.9.1
+ conda version 4.10.3

# インストール
私の場合は、以前の投稿、[M1搭載MacBook Airが来た！手っ取り早くPython計算環境を作ってみる](https://qiita.com/damyarou/items/f0e0b4026bdda94df46f) に示したように、Miniforge3 をインストールしているため、[conda-forge  / packages / cartopy](https://anaconda.org/conda-forge/cartopy) にしたがって

```
conda install -c conda-forge cartopy 
```
により cartopy をインストール。特段問題なくインストール完了。
なお、もちろん、numpy, matplotlib は既にインストールしてある。


# 参考Webページ
#### 公式ページ
[Using cartopy with matplotlib](https://scitools.org.uk/cartopy/docs/v0.15/matplotlib/intro.html)

#### 地図
[Cartopyで地理データを可視化する1](https://metpost.hatenablog.com/entry/2015/11/05/180006)
[Cartopyで地理データを可視化する2](https://metpost.hatenablog.com/entry/2016/06/04/155751)
[Cartopyで地理データを可視化する3](https://metpost.hatenablog.com/entry/2016/06/04/161222)

#### プロットと線
[【Python】　Cartopyでプロット　－マーカー－](https://www.mitsumatado.com/zen/cartopy2/)

#### 軸について
[Tick Labels](https://scitools.org.uk/cartopy/docs/latest/gallery/gridlines_and_labels/tick_labels.html#sphx-glr-gallery-gridlines-and-labels-tick-labels-py)
[cartopyで地図を描くときに軸ラベルとグリッドを整える](https://qiita.com/earth06/items/f5958a89a546dce00c36)


#### 円を描く
[【Science】水平線までどのくらいの距離？Part2ーマップで表示](https://runtascience.hatenablog.com/entry/2020/06/23/%E6%B0%B4%E5%B9%B3%E7%B7%9A%E3%81%BE%E3%81%A7%E3%81%A9%E3%81%AE%E3%81%8F%E3%82%89%E3%81%84%E3%81%AE%E8%B7%9D%E9%9B%A2%EF%BC%9FPart2%E3%83%BC%E3%83%9E%E3%83%83%E3%83%97%E3%81%A7%E8%A1%A8%E7%A4%BA)


# 作図事例

![fig_map.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/129300/c5a3b789-8f91-fb96-31fe-1496e5b8bc67.jpeg)


# プログラム例
上記作図事例を出力するプログラム例。クアラルンプールを中心とした円は、大体あっているはずですが、厳密ではないかもしれません。

```python
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np


def circ(c_lon,c_lat,rr):
    ra=6378.140 # equatorial radius (km)
    rb=6356.755 # polar radius (km)
    da=2*np.pi*ra/360
    db=2*np.pi*rb/360
    theta=np.linspace(0,2*np.pi,180)
    x1, y1=c_lon*da, c_lat*db
    x=rr*np.cos(theta)+x1
    y=rr*np.sin(theta)+y1
    r_lon, r_lat=x/da, y/db
    return r_lon,r_lat


def main():
    scl='10m' # 1:10,000,000
    #scl='50m' # 1:50,000,000
    #scl='110m' # 1:110,000,000
    # physical category
    land  = cfeature.NaturalEarthFeature('physical', 'land', scl, 
                                             edgecolor='face',                   # same color with facecolor
                                             facecolor=cfeature.COLORS['land'])  # use predefiend color of cartopy
    ocean = cfeature.NaturalEarthFeature('physical', 'ocean', scl, 
                                             edgecolor='face',                   # same color with facecolor
                                             facecolor=cfeature.COLORS['water']) # use predefiend color of cartopy
    lakes = cfeature.NaturalEarthFeature('physical', 'lakes', scl, 
                                             edgecolor='face',                   # same color with facecolor
                                             facecolor=cfeature.COLORS['water']) # use predefiend color of cartopy
    river = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', scl,
                                             edgecolor=cfeature.COLORS['water'], # use predefiend color of cartopy
                                             facecolor='none')                   # no filled color
    # cultural category
    countries  = cfeature.NaturalEarthFeature('cultural', 'admin_0_countries', scl, 
                                                  edgecolor='gray',
                                                  facecolor='none')  # no filled color
    fsz=12 # fontsize
    xmin, xmax, dx = 90, 120, 5 # longitude
    ymin, ymax, dy = -10, 20, 5  # latitude
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    iw=10
    ih=iw/(xmax-xmin)*(ymax-ymin)
    plt.figure(figsize=(iw, ih))
    ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent([xmin, xmax, ymin, ymax], ccrs.PlateCarree())
    ax.set_xticks(list(range(xmin,xmax+dx,dx)), crs=ccrs.PlateCarree())
    ax.set_yticks(list(range(ymin,ymax+dy,dy)), crs=ccrs.PlateCarree())
    ax.add_feature(land)
    ax.add_feature(ocean)
    ax.add_feature(lakes)
    ax.add_feature(river)
    ax.add_feature(countries)
    #ax1.coastlines(resolution=scl)
    ax.set_title('Map')
    ax.gridlines(linestyle='-', color='#777777')
    
    c_lon,c_lat=101.688, 3.1357 # Kuala Lumpur
    rr1=500  # radius of circle (1)
    rr2=1000 # radius of circle (2)
    r1_lon,r1_lat=circ(c_lon,c_lat,rr1)
    r2_lon,r2_lat=circ(c_lon,c_lat,rr2)
    ax.plot([c_lon],[c_lat],'+',color='#0000ff')
    ax.plot(r1_lon,r1_lat,'-',color='#0000ff',lw=1)
    ax.plot(r2_lon,r2_lat,'--',color='#0000ff',lw=1)
    ax.text(c_lon,c_lat,'Kuala Lumpur',ha='center',va='bottom')
    ax.text(c_lon,np.max(r1_lat),'{0:.0f}km'.format(rr1),ha='center',va='bottom')
    ax.text(c_lon,np.max(r2_lat),'{0:.0f}km'.format(rr2),ha='center',va='bottom')

    #plt.tight_layout()
    fnameF='fig_map.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


#---------------
# Execute
#---------------
if __name__ == '__main__': main()
```

以　上