一次比较完整的区域放大图

笔者因为一些原因进行了一次区域位置放大图的绘制,由于笔者能力有限,在进行绘制的时候进行了大量参考和借鉴,文末会附上借鉴链接。

"""
作者:久相处
日期:2022年03月03日
"""
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import numpy as np
import shapefile
from mpl_toolkits.axes_grid.inset_locator import mark_inset
from mpl_toolkits.axes_grid.inset_locator import inset_axes
from mpl_toolkits.basemap import Basemap
import matplotlib.patches as mpatches



#  设置画图字体的大小
plt.rcParams.update({'font.size': 20})
#  解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#  解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
#  设置画布和绘画区
fig = plt.figure(figsize=[8, 12])
# 画大图
ax1 = fig.add_subplot(111)
#  读取shp格式的地图,有三个文件分别为.dbf,.shp,.shx缺一不可
#  这里的shp文件其实就是一堆点,把各个国家圈出来
#  sf得到的是一个类,存放了全球的地理信息
sf = shapefile.Reader(r'D:\Shp\country1')
#  sf.shapeRecords()里面放了各个国家地区的信息
#  循环的目的是单个拿出来,找到中国  shape_rec是单个
for shape_rec in sf.shapeRecords():
    # shape_rec.record内部存放了单一国家的很多信息,比如名字,货币等等
    # 其中shape_rec.record[2]放的是国家名
    # 可以print(shape_rec.record)看看
    if shape_rec.record[2] == "China" and "Vietnam":
        codes = []  # 用于存放移动路径(画图动作)
        #  shape_rec.shape是一个类,图形类
        #  里面三个属性shapeType:图形类型 points图形边界坐标 parts图形起始索引
        #  解释一下parts属性, 一个国家的边界条件可能不是全连在一起的, 会分成一块一块,那么就相当于多个图形,存在shp文件不连续,parts里面就放了每个区域的起始索引(下标)
        pts = shape_rec.shape.points
        #  上文已经说过parts的意义,设想一下,两块区域的起始坐标中间夹的不就是一块区域的所有坐标嘛,但是最后一块区域没有结束值,所有加了[len(pts)],这就是最后一个点的索引。
        prt = list(shape_rec.shape.parts)+[len(pts)]
        #  下面的循环主要作用是:建立一个绘图路径,利用区块起始点的索引生成一个画图动作
        for i in range(len(prt)-1):
            codes += [Path.MOVETO]  # 点移动
            codes += [Path.LINETO]*(prt[i+1]-prt[i]-2)  # 画线
            codes += [Path.CLOSEPOLY]  # 这块画布,循环结束,下一块
        clip = Path(pts, codes)  # 利用数据和路径生成一个画图动作
        clip = PathPatch(clip, transform=ax1.transData)  # 再加入ax的变换

#  生成一个Basemap类,用来 变换坐标,画出合适尺寸的图
m = Basemap(llcrnrlon=105.5,  # 地图左边经度
            llcrnrlat=17.5,  # 地图下方纬度
            urcrnrlon=111.5,  # 地图上方经度
            urcrnrlat=23.0,  # 地图上方纬度
            resolution=None,  # 分辨率,这里设置为无
            lat_0=20,
            projection='cyl')  # 投影方式:默认,圆柱投影
#  读取图形文件,画出中国行政边界
m.readshapefile("bou2_4l", "China Map", color="k", linewidth=1.2)
m.readshapefile(r'D:\Shp\vn_shp\VNM_adm0', "Vietnam Map", color='k', linewidth=1.2)#添加了越南区域,若想继续添加其他国家区域,可一直使用m.readshapefile函数
#  再进行一些修饰,让图片更好看
#  画纬度
parallels = np.arange(17.5, 23, 1)
m.drawparallels(parallels, labels=[True, False, True, True], color='dimgrey', dashes=[1, 3])
#  画经度
meridians = np.arange(105.5, 111.5, 2)
m.drawmeridians(meridians, labels=[True, True, False, True], color='dimgrey', dashes=[1, 3])#这里的参数,大家可以适当改动,自然知道这些参数对应的意思
# y移除原来的坐标轴标签
plt.ylabel(" ")
plt.xlabel(" ")
#  设置标题
plt.rcParams.update({"font.size": 30})
ax1.set_title(u'中国西南部地区', color='black', fontsize=25)
#  画南海
with open(r"D:\Shp\CN-border-La.dat") as src:
    context = src.read()
    blocks = [cnt for cnt in context.split(">") if len(cnt) > 0]
    borders = [np.fromstring(block,  dtype=float, sep=' ')for block in blocks]
sub_ax = fig.add_axes([0.628, 0.11, 0.15, 0.15], projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115))
for line in borders:
    sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k', transform=ccrs.Geodetic())
sub_ax.set_extent([106, 127, 0, 25], crs=ccrs.PlateCarree())  # 右下角的小图
# 画指北针
def add_north(ax1, labelsize=14, loc_x=0.92, loc_y=0.9, width=0.03, height=0.1, pad=0.14):
    """
    画一个比例尺带'N'文字注释
    主要参数如下
    :param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可
    :param labelsize: 显示'N'文字的大小
    :param loc_x: 以文字下部为中心的占整个ax横向比例
    :param loc_y: 以文字下部为中心的占整个ax纵向比例
    :param width: 指南针占ax比例宽度
    :param height: 指南针占ax比例高度
    :param pad: 文字符号占ax比例间隙
    :return: None
    """
    minx, maxx = ax1.get_xlim()
    miny, maxy = ax1.get_ylim()
    ylen = maxy - miny
    xlen = maxx - minx
    left = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]
    right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]
    top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]
    center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]
    triangle = mpatches.Polygon([left, top, right, center], color='k')
    ax1.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom')
    ax1.add_patch(triangle)


add_north(ax1)

#  
axins = inset_axes(ax1, width="55%", height="65%", loc='lower left', bbox_to_anchor=(0.96, 0.05, 1, 1), bbox_transform=ax1.transAxes)  # 设置了小图的框
#  -----------------------图像设置完毕-----------------------------
# 1(右上)2(左上)3(左下)4(右下)
mark_inset(ax1, axins, loc1=3, loc2=1, fc='none', ec='k', lw=1)
sf = shapefile.Reader(r'D:\Shp\广西壮族自治区\Guangxi_province')#读取shp地图文件
shapes = sf.shapes()  # 获取shapefile的类
codes = []  # 用来存放移动路径(画图动作)
pts = shapes[0].points  # 边界点,这里只有一个类所以用了shapes[0],没做循环
prt = list(shapes[0].parts) + [len(pts)]  # 区块起始索引
for i in range(len(prt) - 1):
    codes += [Path.MOVETO]  # 点移动
    codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)  # 画线
    codes += [Path.CLOSEPOLY]  # 这块画完,循环结束,下一块
clip = Path(pts, codes)  # 利用数据和路径生成一个画图动作
clip = PathPatch(clip, transform=axins.transData)  # 再加入ax的变换
m = Basemap(llcrnrlon=108,  # 地图左边经度
    llcrnrlat=20.8,  # 地图下方纬度
    urcrnrlon=109.7,  # 地图右边经度
    urcrnrlat =21.9,  # 地图上方纬度
    resolution=None,  # 分辨率,这里设置为无
    projection='cyl')  # 投影方式:默认,圆柱投影
#  ------Basemap类设置完毕--------
# 读取图形文件,画浙江行政边界
# 地级市界
m.readshapefile(r'D:\Shp\ditu\广西壮族自治区\Guangxi_city', 'Guangxi_city Map', color='k', linewidth=1.2)
# 省界
m.readshapefile(r'D:\Shp\ditu\广西壮族自治区\Guangxi_province', 'Guangxi Map', color='k', linewidth=1.2)
# ------Basemap类画地图完毕,下面是修饰工作--------

# 画纬度
parallels = np.arange(21, 22.3, 0.3)
m.drawparallels(parallels, labels=[False, False, True, True], color='dimgrey', dashes=[1, 3], fontsize=15)
# 画经度
meridians = np.arange(107.5, 110, 1.5)
m.drawmeridians(meridians, labels=[True, True, False, False], color='dimgrey', dashes=[1, 3], fontsize=15)

# 移除原来的坐标轴标签
plt.ylabel('')
plt.xlabel('')
# 设置标题
plt.rcParams.update({'font.size': 30})
axins.set_title(u' ', color='black', fontsize=15)
beihai0 = 108.2  # 北海经度
beihai1 = 21.78  # 北海纬度
plt.scatter(beihai0, beihai1, marker='.', s=20, color="green")  # 画个点


fangchenggang0 = 109.35  # 防城港经度
fangchenggang1 = 21.65  # 防城港纬度
plt.scatter(fangchenggang0, fangchenggang1, marker=".", s=20, color="green")  # 画个点

qinzhou0 = 108.75  # 钦州经度
qinzhou1 = 21.7  # 钦州纬度
plt.scatter(qinzhou0, qinzhou1, marker='.', s=20, color="green")  # 画个点

plt.rcParams.update({'font.size': 30})
plt.text(beihai0-0.15,beihai1+0.05, u"北海市", fontsize=10)  # 在斜上方写个字
plt.text(fangchenggang0-0.1, fangchenggang1+0.05, u"防城港市", fontsize=10)  # 在斜上方写个字
plt.text(qinzhou0-0.05, qinzhou1+0.1, u"钦州市", fontsize=10)  # 在斜上方写个字
plt.show()

写此文章,也是为了日后给他人行个方便,也是对于自己的一个记录,虽然最后还是被老师骂了.

最后形成图片如下:

 文章有几点遗憾:对于Basemap画图普遍采用cyl圆柱型投影,而笔者想要用m.drawmapscale添加比例尺但是提示不能用。而网上有基于墨卡托投影merc进行添加,奈何笔者水平有限,未能成功。

 

 在这里特别感谢博主:野生的气象小流星对我的各种帮助,给我提供资料,并且耐心教育我,其文章也让我收益匪浅。

最后附上参考文章链接。本文只为记录和给别人行个方便,若要学习,还望点击下方链接处,一探究竟。

小白学习Basemap气象画地图的第三天(中国温度分布图,mask外部)使Basemap添加指南针,显示经纬度刻度

 【Matplotlib】 局部放大图

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值