前言

氣象、海洋、地球科學相關的分析必然少不了地圖的可視化。
Python中我常用的繪制地圖和空間信息分析的庫是Cartopy。
Cartopy有一個非常嚴重的問題，那就是自帶的中國國界數據有問題，這也是很多國外開源庫的普遍問題。
在做中國區域的分析時，由于九段線的位置很偏南，因此最標準的做法是同時繪制南海區域的子圖。
在做一些站點展示的時候，如果只單獨畫上幾個站點總覺得很丑，可以加上一些地形背景。

綜上，今天想要用一個小例子解決這3個問題：

1. 正確的中國國界線及九段線繪制
2. 南海小地圖繪制
3. 全球地形圖添加

準備工作

1. 獲取正確的中國矢量文件：公眾號后臺留言“中國行政區劃”    
   (這個矢量文件來自資源環境平臺，并和權威機構的標準地圖做了比對，吻合一致。)
2. 獲取全球地形圖像：公眾號后臺留言“全球地形”    
   (提供的是全球50m分辨率的tif圖，如果對分辨率要求不高可以直接使用stock_img())

代碼

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.io.shapereader import Reader, natural_earth
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from matplotlib.image import imread


def create_map():
    shp_path = './cn_shp/Province_9/'
    # --創建畫圖空間
    proj = ccrs.PlateCarree()  # 創建坐標系
    fig = plt.figure(figsize=(6, 8), dpi=400)  # 創建頁面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  

    # --設置地圖屬性
    provinces = cfeat.ShapelyFeature(
        Reader(shp_path + 'Province_9.shp').geometries(),
        proj, edgecolor='k',
        facecolor='none'
    )
    # 加載省界線
    ax.add_feature(provinces, linewidth=0.6, zorder=2)
    # 加載分辨率為50的海岸線
    ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=10)
    # 加載分辨率為50的河流~
    ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=10)
     # 加載分辨率為50的湖泊
    ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=10)
    ax.set_extent([105, 133, 15, 45])
    # ax.stock_img()
    ax.imshow(
        imread('./NE1_50M_SR_W.tif'),
        origin='upper',
        transform=proj,
        extent=[-180, 180, -90, 90]
    )
    # --設置網格點屬性
    gl = ax.gridlines(
        crs=ccrs.PlateCarree(),
        draw_labels=True,
        linewidth=1.2,
        color='k',
        alpha=0.5,
        linestyle='--'
    )
    gl.xlabels_top = False  # 關閉頂端的經緯度標簽
    gl.ylabels_right = False  # 關閉右側的經緯度標簽
    gl.xformatter = LONGITUDE_FORMATTER  # x軸設為經度的格式
    gl.yformatter = LATITUDE_FORMATTER  # y軸設為緯度的格式
    gl.xlocator = mticker.FixedLocator(np.arange(95, 145 + 5, 5))
    gl.ylocator = mticker.FixedLocator(np.arange(-5, 45 + 5, 5))

    # --設置小地圖
    left, bottom, width, height = 0.67, 0.15, 0.23, 0.27
    ax2 = fig.add_axes(
        [left, bottom, width, height], 
        projection=proj
    )
    ax2.add_feature(provinces, linewidth=0.6, zorder=2)
    ax2.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=10)
    ax2.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=10)
    ax2.add_feature(cfeat.LAKES.with_scale('50m'), zorder=10)
    ax2.set_extent([105, 125, 0, 25])
    # ax2.stock_img()
    ax2.imshow(
        imread('./NE1_50M_SR_W.tif'),
        origin='upper',
        transform=proj,
        extent=[-180, 180, -90, 90]
    )
    return ax

def main():
    ax = create_map()
    title = f'distribution of station around China'
    ax.set_title(title, fontsize=18)

    df = pd.read_csv('buyo_position.csv')
    df['lon'] = df['lon'].astype(np.float64)
    df['lat'] = df['lat'].astype(np.float64)
    ax.scatter(
        df['lon'].values,
        df['lat'].values,
        marker='o',
        s=10 ,
        color ="blue"
    )
    for i, j, k in list(zip(df['lon'].values, df['lat'].values, df['name'].values)):
        ax.text(i - 0.8, j + 0.2, k, fontsize=6)
    plt.savefig('station_distribute_map.png')

if __name__ == '__main__':
    main()