笔者因为一些原因进行了一次区域位置放大图的绘制,由于笔者能力有限,在进行绘制的时候进行了大量参考和借鉴,文末会附上借鉴链接。
"""
作者:久相处
日期: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进行添加,奈何笔者水平有限,未能成功。
在这里特别感谢博主:野生的气象小流星对我的各种帮助,给我提供资料,并且耐心教育我,其文章也让我收益匪浅。
最后附上参考文章链接。本文只为记录和给别人行个方便,若要学习,还望点击下方链接处,一探究竟。