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計算環境を作ってみる に示したように、Miniforge3 をインストールしているため、conda-forge / packages / cartopy にしたがって

conda install -c conda-forge cartopy 

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

参考Webページ

公式ページ

Using cartopy with matplotlib

地図

Cartopyで地理データを可視化する1
Cartopyで地理データを可視化する2
Cartopyで地理データを可視化する3

プロットと線

【Python】　Cartopyでプロット　－マーカー－

軸について

Tick Labels
cartopyで地図を描くときに軸ラベルとグリッドを整える

円を描く

【Science】水平線までどのくらいの距離？Part2ーマップで表示

作図事例

プログラム例

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

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()

以　上