[学习笔记]python绘制图中图（绘制站点分布图）

背景

在绘制站点分布图时，有时需要采用图中图的方式，以便于在一张图中尽可能多的表达信息。此处记录一下利用python matplotlib绘制图中图的脚本，方便然后查询。

包含数据

该绘图脚本中包含以下数据：

CMONOC站点分布（蓝色点）
CMONOC穿刺点分布（灰色点）
某研究中采用的位于湖北省附近的一些地面跟踪站分布（红色点）

绘制思路

首先将CMONOC站点以及穿刺点分布画上，再通过plt.axes在图上再加一块画布绘制小范围的地图，接着在小范围的地图上标点。为了方便在大范围地图上找到图中图的位置，也要在相应位置上用红框标记下，这里采用了plt.vlines和plt.hlines两者结合的方法。

代码

注：此处仅提供绘图脚本，其中用到的部分站点坐标文件因一些原因不能公开。

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib.patches as mpathes
import numpy as np
import os
from adjustText import adjust_text  # 导入文字调整的库函数

outputDir = './cmonoc_ipp/'
if not os.path.exists(outputDir):  # True/False
    os.mkdir(outputDir)

colors_lst = ['blue', 'magenta', 'darkolivegreen', 'mediumpurple', 'palevioletred', 'cadetblue']
site_list = r'xxxxxxxxxxxx/ofilelst.txt' # 测站名列表路径
ccl_file_dir = r'xxxxxxxxxxxxxxxxxxxxx/cclFileCmonoc/' # 穿刺点文件路径

# 读取测站列表文件（含经纬度及站点名称）
f = open(site_list, 'r')
ObsLines = f.readlines()

rec_name = []
for i in range(len(ObsLines)):
    fields = ObsLines[i].split()
    rec_name.append(fields[0])

plt.figure(figsize=(6, 4.5))  # 设置大小和分辨率
plt.rcParams['font.sans-serif'] = ['arial']

lat_range = range(-15, 60 + 10, 10)
lon_range = range(60, 150 + 10, 10)

m = Basemap(projection='cyl', lon_0=110, lat_0=20, resolution='h', llcrnrlon=60, urcrnrlon=150,
            llcrnrlat=-10, urcrnrlat=60)
m.drawcoastlines(color='black', linewidth=0.8)
# draw parallels and meridians.
m.drawmeridians(range(70, 150 + 20, 20), labels=[0, 0, 1, 1], color='gray', linewidth=0.8, font='arial', fontsize=10)
m.drawparallels(range(0, 60, 10), labels=[1, 1, 0, 0], color='gray', linewidth=0.8, font='arial', fontsize=10)

for i in range(len(rec_name)):

    site = rec_name[i][2:6]
    print('Porcessing site: ', site, '...')

    file_name = site + '_2023_091.ccl'
    arc_file = ccl_file_dir + file_name
    print('Processing the file: ', file_name)
    sec = []
    ipp = []
    elev = []
    site_lst = []

    with open(arc_file, "r") as f:
        for line in f.readlines():
            line = line.split()
            # if line[4] != 'W05':
            #     continue
            sec.append(int(line[2]))
            ipp.append([float(line[20]), float(line[21])])
            elev.append(float(line[5]))
            site_lst.append(site)

    lon = []
    lat = []
    for k in range(len(sec)):
        lon_tmp, lat_tmp = m(ipp[k][1], ipp[k][0])
        lon.append(lon_tmp)
        lat.append(lat_tmp)

    plt.scatter(lon, lat, s=0.0001, c='silver', zorder=100)
    # plt.plot(lon, lat, marker='o', color='grey', markersize=0.1, zorder=100)

# CMONOC坐标列表路径
listFile = r'E:/DoctoralStudy/2python_prog/plotSiteMap0615/sitelist/site_pos_cmonoc.txt'

# 读取测站列表文件（含经纬度及站点名称）
f = open(listFile, 'r')
ObsLines = f.readlines()

lon = []
lat = []
staname = []
for i in range(len(ObsLines)):
    fields = ObsLines[i].split()
    # print(fields[0], fields[1], fields[2])
    lon.append(float(fields[0]))
    lat.append(float(fields[1]))
    staname.append(fields[2])

lon, lat = m(lon, lat)
# ----------在地图上绘制坐标点------------#
for i in range(len(staname)):
    # m.scatter(lon[i],
    #           lat[i],
    #           s=20,
    #           c='red',
    #           marker='o')
    plt.plot(lon[i], lat[i], marker='s', color='blue', markersize=3, zorder=100) # CMONOC站点

# 标记图中图位置，红色线
plt.vlines(108, ymin=24, ymax=36, colors='r', zorder=100)
plt.vlines(124, ymin=24, ymax=36, colors='r', zorder=100)
plt.hlines(24, xmin=108, xmax=124, colors='r', zorder=100)
plt.hlines(36, xmin=108, xmax=124, colors='r', zorder=100)

# 绘制图中图
plt.axes([0.54, 0.15, 0.34, 0.28])
plt.rcParams['axes.facecolor'] = 'white'
m = Basemap(projection='cyl', lon_0=110, lat_0=20, resolution='h', llcrnrlon=108, urcrnrlon=124, llcrnrlat=24,
            urcrnrlat=36)
m.drawcoastlines(color='grey', linewidth=0.8)
# 图中图经纬度标注受大图影响，标注时打开bbox使其有白色背景，增强可读性
m.drawmeridians(np.arange(100, 124+4, 4), labels=[0, 0, 0, 1], color='gray', linewidth=0.8, font='arial',
                fontsize=10, bbox=dict(facecolor="white", edgecolor="white", pad=0.4))
m.drawparallels(np.arange(26, 36+4, 4), labels=[1, 0, 0, 0], color='gray', linewidth=0.8, font='arial',
                fontsize=10, bbox=dict(facecolor="white", edgecolor="white", pad=0.4))


# 读取测站列表文件（含经纬度及站点名称）
listFile = r'D:/TempDataBackup/LEOData/091DataProc/site_pos_leo.txt'
f = open(listFile, 'r')
ObsLines = f.readlines()

lon = []
lat = []
staname = []
for i in range(len(ObsLines)):
    fields = ObsLines[i].split()
    # print(fields[0], fields[1], fields[2])
    lon.append(float(fields[0]))
    lat.append(float(fields[1]))
    staname.append(fields[2])

lon, lat = m(lon, lat)
# ----------在地图上绘制坐标点，添加文字------------#
for i in range(len(staname)):
    # m.scatter(lon[i],
    #           lat[i],
    #           s=20,
    #           c='red',
    #           marker='o')
    plt.plot(lon[i], lat[i], marker='o', color='red', markersize=5, zorder=100)
# 此处可以直接使用 m.scatter(lon,lat)

# ---------在坐标点添加文字------#
texts = []
for i in range(len(staname)):
    texts.append(
        plt.text(lon[i],  # 坐标
                 lat[i],
                 staname[i],  # 文字字符串
                 font='arial',
                 fontsize=10,
                 style="italic",
                 weight="normal",
                 verticalalignment='center',
                 horizontalalignment='right',
                 rotation=0, zorder=100))  # 给散点加标签
adjust_text(texts, only_move={'text': 'xy'},)

fig_name = outputDir + 'cmonoc091IPPTrack.tif'
plt.savefig(fig_name, dpi=600, bbox_inches='tight', pad_inches=0.1)  # 输出地图，并设置边框空白紧密
plt.show()  # 将图像显示出来