1.简介
最近在做实验,想用geopandas和matplotlib绘制研究区的地图,效果如下:
弄了一天,只绘制了比较简洁的地图要素,有些字体也没有调整过来,希望有大佬可以指教一下!!
难点是绘制指北针和研究区格网划分,需要预先定义函数。
1.1格网划分:
def make_grid(boundary,radius):# 构建正方形格网
# 以边界的中心点按照radius半径进行外扩
minx, miny, maxx, maxy = boundary.total_bounds # 边界范围
# 求边界的中心点
deltax = ((maxx-minx)%(radius*2))/2
deltay = ((maxy-miny)%(radius*2))/2
# 以边界的中心点按照radius半径进行外扩
minx -= deltax
maxx += deltax
miny -= deltay
maxy += deltay
cells=[] # 格网列表
for x in np.arange(minx,maxx,radius*2): # 以radius*2为步长
for y in np.arange(miny,maxy,radius*2): # 以radius*2为步长
cells.append(geo.Polygon([(x,y),(x,y+radius*2),(x+radius*2,y+radius*2),(x+radius*2,y),(x,y)])) # 构建格网
return gpd.GeoDataFrame(cells,columns=['geometry'],crs=boundary.crs).sjoin(boundary,how='inner') # 与边界进行空间连接
1.2指北针
def add_compass(ax, labelsize=10, loc_x=0.08, loc_y=0.92, width=0.05, height=0.05):
# labelsize控制指北针标签的字体大小,loc_x和loc_y控制指北针的位置,width和height控制指北针的大小
minx, maxx = ax.get_xlim()
miny, maxy = ax.get_ylim()
ylen = maxy - miny
xlen = maxx - minx
# 外接矩形的四个点
LTop = [minx + xlen*(loc_x - width), miny + ylen*(loc_y + height)]
RTop = [minx + xlen*(loc_x + width), miny + ylen*(loc_y + height)]
LBottom = [minx + xlen*(loc_x - width), miny + ylen*(loc_y - height)]
RBottom = [minx + xlen*(loc_x + width), miny + ylen*(loc_y - height)]
# 内接矩形的四个点
LTop1 = [minx + xlen*(loc_x - width*.18), miny + ylen*(loc_y + height*.18)]
RTop1 = [minx + xlen*(loc_x + width*.18), miny + ylen*(loc_y + height*.18)]
LBottom1 = [minx + xlen*(loc_x - width*.18), miny + ylen*(loc_y - height*.18)]
RBottom1 = [minx + xlen*(loc_x + width*.18), miny + ylen*(loc_y - height*.18)]
# 十字型指北针的五个点
Top = [(LTop[0] + RTop[0])/2, LTop[1]]
Bottom = [(LBottom[0] + RBottom[0])/2, LBottom[1]]
Left = [LTop[0], (LTop[1] + LBottom[1])/2]
Right = [RTop[0], (RTop[1] + RBottom[1])/2]
Center = [(LTop[0] + RTop[0])/2, (LBottom[1] + LTop[1])/2]
# 绘制指北针
triangle1 = mpatches.Polygon([Left,Center,LTop1], edgecolor='black', fill=False, alpha=0.8)
triangle2 = mpatches.Polygon([Left,Center,LBottom1], color='black', alpha=0.7)
triangle3 = mpatches.Polygon([Top,Center,LTop1], color='black', alpha=0.7)
triangle4 = mpatches.Polygon([Top,Center,RTop1], edgecolor='black', fill=False, alpha=0.8)
triangle5 = mpatches.Polygon([Right,Center,RTop1], color='black', alpha=0.7)
triangle6 = mpatches.Polygon([Right,Center,RBottom1], edgecolor='black', fill=False, alpha=0.8)
triangle7 = mpatches.Polygon([Bottom,Center,RBottom1], color='black', alpha=0.7)
triangle8 = mpatches.Polygon([Bottom,Center,LBottom1], edgecolor='black', fill=False, alpha=0.8)
ax.add_patch(triangle1)
ax.add_patch(triangle2)
ax.add_patch(triangle3)
ax.add_patch(triangle4)
ax.add_patch(triangle5)
ax.add_patch(triangle6)
ax.add_patch(triangle7)
ax.add_patch(triangle8)
# 绘制指北针标签
ax.text(s='N',
x=minx + xlen*loc_x,
y=miny + ylen*(loc_y + height),
fontsize=labelsize,
horizontalalignment='center', # 水平对齐方式
verticalalignment='bottom') # 垂直对齐方式
ax.text(s='S',
x=minx + xlen*loc_x,
y=miny + ylen*(loc_y - height),
fontsize=labelsize,
horizontalalignment='center', # 水平对齐方式
verticalalignment='top') # 垂直对齐方式
ax.text(s='W',
x=minx + xlen*(loc_x - width),
y=miny + ylen*(loc_y),
fontsize=labelsize,
horizontalalignment='right', # 水平对齐方式
verticalalignment='center') # 垂直对齐方式
ax.text(s='E',
x=minx + xlen*(loc_x + width),
y=miny + ylen*(loc_y),
fontsize=labelsize,
horizontalalignment='left', # 水平对齐方式
verticalalignment='center') # 垂直对齐方式
2.完整代码
'''导入包'''
import numpy as np
import geopandas as gpd
from xpinyin import Pinyin
from shapely import geometry as geo
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
'''读取数据'''
gdf = gpd.read_file('data/上海市内三环.shp') # 上海市内三环——线
gdf2 = gpd.read_file('data/310000_full.shp',encoding='utf-8') # 上海市——面
'''预处理'''
def make_grid(boundary,radius):# 构建正方形格网
# 以边界的中心点按照radius半径进行外扩
minx, miny, maxx, maxy = boundary.total_bounds # 边界范围
# 求边界的中心点
deltax = ((maxx-minx)%(radius*2))/2
deltay = ((maxy-miny)%(radius*2))/2
# 以边界的中心点按照radius半径进行外扩
minx -= deltax
maxx += deltax
miny -= deltay
maxy += deltay
cells=[] # 格网列表
for x in np.arange(minx,maxx,radius*2): # 以radius*2为步长
for y in np.arange(miny,maxy,radius*2): # 以radius*2为步长
cells.append(geo.Polygon([(x,y),(x,y+radius*2),(x+radius*2,y+radius*2),(x+radius*2,y),(x,y)])) # 构建格网
return gpd.GeoDataFrame(cells,columns=['geometry'],crs=boundary.crs).sjoin(boundary,how='inner') # 与边界进行空间连接
# 根据闭合的线获取面然后进行格网划分
gdf3 = gdf[gdf['Name']=='外环线'].copy()
gdf4 = gdf3['geometry'].apply(lambda x: geo.Polygon(x)) # 外环面
gdf4 = gpd.GeoDataFrame(gdf4, geometry='geometry', crs='EPSG:4326')
gdf5 = make_grid(gdf4.to_crs(32651), 500) # 1km格网
# 区域名称转拼音
p = Pinyin()
gdf2['name1'] = gdf2['name'].apply(lambda x: p.get_pinyin(x, '')) # 拼音
gdf2['name2'] = gdf2['name1'].apply(lambda x: x[:-2]) # 去掉后缀
gdf2['name3'] = gdf2['name2'].apply(lambda x: x.capitalize()) # 首字母大写
# 密集区域换成简写
gdf2.loc[gdf2['name3'] == 'Xuhui', 'name3'] = 'XH'
gdf2.loc[gdf2['name3'] == 'Yangpu', 'name3'] = 'YP'
gdf2.loc[gdf2['name3'] == 'Hongkou', 'name3'] = 'HK'
gdf2.loc[gdf2['name3'] == 'Huangpu', 'name3'] = 'HP'
gdf2.loc[gdf2['name3'] == 'Changning', 'name3'] = 'CN'
gdf2.loc[gdf2['name3'] == 'Putuo', 'name3'] = 'PT'
gdf2.loc[gdf2['name3'] == 'Jingan', 'name3'] = 'JA'
'''绘制地图'''
# 画布
fig = plt.figure(figsize=(6, 4))
plt.rc('font',family='Times New Roman') # 设置全局字体
'''主图'''
# 创建填满画布的坐标轴
ax = fig.add_axes((0, 0, 1.5, 1.5)) # (左, 下, 右, 上)
ax = gdf2.boundary.plot(ax=ax,
edgecolor='black', # 边界颜色
linewidth=0.5, # 边界宽度
alpha=0.8) # 透明度
ax = gdf4.plot(ax=ax,
facecolor='grey', # 填充颜色
edgecolor='white', # 边界颜色
linewidth=0.5, # 边界宽度
alpha=0.6) # 透明度
# 注释
ax.text(x=120.71,#文本x轴坐标
y=31.95, #文本y轴坐标
s='Shanghai City', #文本内容
ha='left', # 文字水平对齐方式
va='baseline', # 文字垂直对齐方式
fontdict=dict(fontsize=15,
family='serif',
weight='heavy') # 设置字体属性
)
ax.text(x=120.78,#文本x轴坐标
y=30.6, #文本y轴坐标
s='* YP:Yangpu HK:Hongkou HP:Huangpu CN:Changning PT:Putuo XU:Xuhui JA:Jingan', #文本内容
ha='left', # 文字水平对齐方式
va='baseline', # 文字垂直对齐方式
fontdict=dict(fontsize=8,
family='serif',
weight='heavy') # 设置字体属性
)
ax.text(x=120.763,#文本x轴坐标
y=30.56, #文本y轴坐标
s='** The study area in Shanghai is denoted in grey', #文本内容
ha='left', # 文字水平对齐方式
va='baseline', # 文字垂直对齐方式
fontdict=dict(fontsize=8,
color='black',
family='serif',
weight='heavy') # 设置字体属性
)
# 显示每个区的名称
for idx, _ in enumerate(gdf2.representative_point()):
# 提取名称
region = gdf2.loc[idx, 'name3'] # 从gdf2中提取区域名称
ax.text(_.x, _.y, region, ha="center", va="center", size=10) # 绘制区域名称
# 设置图布范围
ax.set_xlim(120.7, 122.5)
ax.set_ylim(30.5, 32)
# 移除坐标轴
ax.axis('off')
'''子图'''
# 绘制子图
ax_child = fig.add_axes((1, 0, 1.5, 1.5)) # (左, 下, 右, 上)
ax_child = gdf5.to_crs(4326).geometry.plot(ax=ax_child,
facecolor='grey', # 填充颜色
edgecolor='white', # 边界颜色
label='Urban Districts', # 图例名称
alpha=0.8) # 透明度
ax_child = gdf[gdf['Name']!='外环线'].plot(ax=ax_child,
label='Ring Road', # 图例名称
color='black', # 边界颜色
linewidth=0.8, # 边界宽度
alpha=0.8) # 透明度
# 设置图布范围
ax_child.set_ylim(31.07, 31.41)
# 比例尺
scalebar = ScaleBar(dx=100000, # 1格网 = 0.2m
dimension="si-length", # 比例尺单位
units="m", # 比例尺单位
location="lower left", # 比例尺位置
frameon=False) # 是否显示边框
ax_child.add_artist(scalebar) # 添加比例尺
# 指北针
def add_compass(ax, labelsize=7, loc_x=0.07, loc_y=0.93, width=0.04, height=0.04):
# labelsize控制指北针标签的字体大小,loc_x和loc_y控制指北针的位置,width和height控制指北针的大小
minx, maxx = ax.get_xlim()
miny, maxy = ax.get_ylim()
ylen = maxy - miny
xlen = maxx - minx
# 外界矩形的四个点
LTop = [minx + xlen*(loc_x - width), miny + ylen*(loc_y + height)]
RTop = [minx + xlen*(loc_x + width), miny + ylen*(loc_y + height)]
LBottom = [minx + xlen*(loc_x - width), miny + ylen*(loc_y - height)]
RBottom = [minx + xlen*(loc_x + width), miny + ylen*(loc_y - height)]
# 内接矩形的四个点
LTop1 = [minx + xlen*(loc_x - width*.18), miny + ylen*(loc_y + height*.18)]
RTop1 = [minx + xlen*(loc_x + width*.18), miny + ylen*(loc_y + height*.18)]
LBottom1 = [minx + xlen*(loc_x - width*.18), miny + ylen*(loc_y - height*.18)]
RBottom1 = [minx + xlen*(loc_x + width*.18), miny + ylen*(loc_y - height*.18)]
# 十字型指北针的五个点
Top = [(LTop[0] + RTop[0])/2, LTop[1]]
Bottom = [(LBottom[0] + RBottom[0])/2, LBottom[1]]
Left = [LTop[0], (LTop[1] + LBottom[1])/2]
Right = [RTop[0], (RTop[1] + RBottom[1])/2]
Center = [(LTop[0] + RTop[0])/2, (LBottom[1] + LTop[1])/2]
# 绘制指北针
triangle1 = mpatches.Polygon([Left,Center,LTop1], edgecolor='black', fill=False, alpha=0.8)
triangle2 = mpatches.Polygon([Left,Center,LBottom1], color='black', alpha=0.7)
triangle3 = mpatches.Polygon([Top,Center,LTop1], color='black', alpha=0.7)
triangle4 = mpatches.Polygon([Top,Center,RTop1], edgecolor='black', fill=False, alpha=0.8)
triangle5 = mpatches.Polygon([Right,Center,RTop1], color='black', alpha=0.7)
triangle6 = mpatches.Polygon([Right,Center,RBottom1], edgecolor='black', fill=False, alpha=0.8)
triangle7 = mpatches.Polygon([Bottom,Center,RBottom1], color='black', alpha=0.7)
triangle8 = mpatches.Polygon([Bottom,Center,LBottom1], edgecolor='black', fill=False, alpha=0.8)
ax.add_patch(triangle1)
ax.add_patch(triangle2)
ax.add_patch(triangle3)
ax.add_patch(triangle4)
ax.add_patch(triangle5)
ax.add_patch(triangle6)
ax.add_patch(triangle7)
ax.add_patch(triangle8)
# 绘制指北针标签
ax.text(s='N',
x=minx + xlen*loc_x,
y=miny + ylen*(loc_y + height),
fontsize=labelsize,
horizontalalignment='center', # 水平对齐方式
verticalalignment='bottom') # 垂直对齐方式
ax.text(s='S',
x=minx + xlen*loc_x,
y=miny + ylen*(loc_y - height),
fontsize=labelsize,
horizontalalignment='center', # 水平对齐方式
verticalalignment='top') # 垂直对齐方式
ax.text(s='W',
x=minx + xlen*(loc_x - width),
y=miny + ylen*(loc_y),
fontsize=labelsize,
horizontalalignment='right', # 水平对齐方式
verticalalignment='center') # 垂直对齐方式
ax.text(s='E',
x=minx + xlen*(loc_x + width),
y=miny + ylen*(loc_y),
fontsize=labelsize,
horizontalalignment='left', # 水平对齐方式
verticalalignment='center') # 垂直对齐方式
add_compass(ax_child) # 添加指北针
# 设置图例
# 创建自定义图例对象
line_patch = mlines.Line2D([], [], color='black', linewidth=0.8, label='Ring Road') # 线图例对象
polygon_patch = mpatches.Patch(facecolor='grey', edgecolor='white', label='Urban Districts') # 面图例对象
legend = ax_child.legend(handles=[line_patch,polygon_patch], # 图例对象列表
loc="lower right", # 图例位置
ncol=1, # 图例列数
shadow=False) # 是否带有阴影
legend.set_title('Legend', {'size': 10, 'family': 'serif', 'weight': 'heavy'}) # 设置图例标题和大小
# 注释
ax_child.text(x=121.605,#文本x轴坐标
y=31.395, #文本y轴坐标
s='Study Area', #文本内容
ha='left', # 文字水平对齐方式
va='baseline', # 文字垂直对齐方式
fontdict=dict(fontsize=10,
color='black',
family='serif',
weight='heavy') # 设置字体属性
)
# 转度分秒的函数
def dms_formatter(x, pos, direction): # 方向: 'latitude' or 'longitude'
degrees = int(x)
minutes = int((x - degrees) * 60)
seconds = int(((x - degrees) * 60 - minutes) * 60)
if direction == "latitude":
return f"{degrees}°{minutes}'{seconds}\"N"
elif direction == "longitude":
return f"{degrees}°{minutes}'{seconds}\"E"
# 获取地图范围
x_min, x_max = ax_child.get_xlim()
y_min, y_max = ax_child.get_ylim()
# 刻度
x = (x_min + x_max) / 2
y = (y_min + y_max) / 2
# 计算刻度位置
x_ticks = [x - 0.1, x, x + 0.1]
y_ticks = [y - 0.15, y - 0.075, y, y + 0.075, y + 0.15]
# 设置刻度标签
ax_child.set_xticks(x_ticks) # 设置x轴刻度
ax_child.set_yticks(y_ticks) # 设置y轴刻度
# 设置刻度标签文本
ax_child.set_xticklabels([dms_formatter(x, None, "longitude") for x in x_ticks]) # 设置x轴刻度标签文本
ax_child.set_yticklabels([dms_formatter(y, None, "latitude") for y in y_ticks]) # 设置y轴刻度标签文本
# 设置纬度刻度标签为竖向,居中对齐
ax_child.set_yticklabels(ax_child.get_yticklabels(),
rotation=90, # 设置刻度标签旋转角度
va="center") # 设置刻度标签对齐方式
# 显示刻度标签和刻度线
ax_child.tick_params(axis="both", # 设置刻度线和刻度标签
which="both", # 设置刻度线和刻度标签
labeltop=True, # 设置上刻度标签
labelbottom=True, # 设置下刻度标签
labelleft=True, # 设置左刻度标签
labelright=True, # 设置右刻度标签
top=True, # 设置上刻度线
bottom=True, # 设置下刻度线
left=True, # 设置左刻度线
right=True) # 设置右刻度线
# 保存地图
plt.savefig('data/研究区.png', dpi=300, bbox_inches='tight')
3.结果
实验数据:https://pan.baidu.com/s/1JsTOB5ngC1USDRmJj0y1_w?pwd=1111