Python实现绘制Matlab格式的地图边框的示例代码

顾星河

发布时间：2026-10-16 | 3人浏览过

原创

本文介绍如何使用Python的Cartopy库绘制具有Matlab风格的地图边框。内容涵盖色斑图绘制、经纬度插值、比例尺与指南针的添加方法，以及自定义边框样式的实现。通过具体代码示例，展示如何配置字体、处理地理数据并生成高清地图输出，适用于地理信息可视化需求。

首先，使用Python绘制色斑图。代码中引入了matplotlib.pyplot、numpy、cartopy相关模块（如ccrs、cfeature）以及pandas和maskout等库。配置了字体为Times New Roman，字号16，数学文本字体集为stix。读取河南地区的Excel数据，获取纬度、经度和PM2.5数据。利用scipy.interpolate中的Rbf函数进行经纬度插值，生成网格数据。插值方法可选linear、nearest、cubic等。


import matplotlib.pyplot as plt
import numpy as np
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import pandas as pd
import maskout
import matplotlib.patches as mpatches
from cartopy.io.shapereader import Reader
from matplotlib import rcParams
config = {"font.family":'Times New Roman',"font.size": 16,"mathtext.fontset":'stix'}
rcParams.update(config)
df1 = pd.read_excel(r"F:/Rpython/lp37/henanmap/henan.xlsx")
lat = df1['lat']
lon = df1['lon']
PM25 = df1['PM25']
# 设置经纬度并用scipy包中的函数插值
from scipy.interpolate import Rbf
olon = np.linspace(110, 117, 100)
olat = np.linspace(31, 37, 100)
olon, olat = np.meshgrid(olon, olat)
# 插值处理,‘linear'，‘nearest'，‘cubic'
# cubic, gaussian, inverse_multiquadric, linear, multiquadric, quintic, thin_plate
# rain_data_new = griddata((lon,lat), data, (olon,olat), method='linear')
func1=Rbf(lon,lat,PM25,function='linear')
PM25=func1(olon,olat)
print(olon)
print(olat)
print(PM25)
print(olat.shape)
print(PM25.shape)
# 建立画布
fig2 = plt.figure(figsize=(16, 12))
proj = ccrs.PlateCarree()
leftlon, rightlon, lowerlat, upperlat = (110, 117, 31, 37)  # 根据上下限确定范围，至少为10°
# 春 ------------------------------------------------------------------------------------
#左 底 宽 高
ax=fig2.add_subplot(1,1,1,projection=ccrs.PlateCarree())
# 在画布的绝对坐标建立子图
ax.set_extent([leftlon, rightlon, lowerlat, upperlat],crs=ccrs.PlateCarree())
# 设置地图属性:加载国界、海岸线，land 为灰色
ax.stock_img()
# 绘制河南省行政边界
ticks=np.arange(30,80,5)
cf=ax.contourf(olon,olat,PM25,levels=ticks,cmap='gist_rainbow',transform=ccrs.PlateCarree(),extend='both')
clip1=maskout.shp2clip(cf,ax,'F:/Rpython/lp37/henanmap/henan1')
cf=plt.colorbar(cf,ticks=ticks,shrink=0.96,orientation='vertical',extend='both',pad=0.01,aspect=35)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth=0.5,zorder=2,color='k')# 添加海岸线
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_geometries(Reader(r'F:/Rpython/lp37/henanmap/henan1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2,zorder=1)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/china1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='r',linewidth=1.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/china2.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='r',linewidth=0.8)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/1级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='RoyalBlue',linewidth=0.4)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/2级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DodgerBlue',linewidth=0.3)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/3级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DeepSkyBlue',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/4级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='SkyBlue',linewidth=0.15)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/5级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='LightSkyBlue',linewidth=0.05)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/主要湖泊.shp').geometries(),ccrs.PlateCarree(),edgecolor='none',linewidth=0,facecolor='#BEE8FF')
# 以下6条语句是定义地理坐标标签格式
ax.set_xticks(np.arange(leftlon, rightlon+0.1, 1),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat, upperlat+0.1, 1),crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.savefig('F:/Rpython/lp37/plot185.4.png',dpi=600,bbox_inches='tight',pad_inches=0)
plt.show()

接着，在地图上添加比例尺和指南针。定义add_north函数，用于绘制带'N'文字注释的指南针。参数包括坐标轴ax、文字大小labelsize、位置loc_x和loc_y、宽度width、高度height和间隙pad。函数计算坐标轴的范围，确定指南针三角形和箭头的顶点坐标，并绘制多边形和线条。同时，添加文字注释显示'N'。


# 添加比例尺，指南针
def add_north(ax, labelsize=18, loc_x=0.95, loc_y=0.99, width=0.06, height=0.09, 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 = ax.get_xlim()
    miny, maxy = ax.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')
    ax.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom')
    ax.add_patch(triangle)
#-----------函数：添加比例尺--------------
def add_scalebar(ax,lon0,lat0,length,size=0.01):
    '''
    ax: 坐标轴
    lon0: 经度
    lat0: 纬度
    length: 长度
    size: 控制粗细和距离的
    '''
    # style 3
    ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1,linewidth=3, label='%d km' % (length))
    ax.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1,linewidth=3)
    ax.vlines(x = lon0+length/2/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1,linewidth=3)
    ax.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1,linewidth=3)
    ax.text(lon0+length/111,lat0+size+0.05,'%d' % (length),horizontalalignment = 'center')
    ax.text(lon0+length/2/111,lat0+size+0.05,'%d' % (length/2),horizontalalignment = 'center')
    ax.text(lon0,lat0+size+0.05,'0',horizontalalignment = 'center')
    ax.text(112.28,31.45,'km',horizontalalignment = 'center')

然后，绘制Matlab格式的地图边框。定义drow_the_scale函数，用于绘制比例尺。参数包括纬度y、起始经度x、文本text、长度length和线宽lw。函数计算步长，绘制黑白相间的水平线和垂直刻度线。根据经纬度范围，循环绘制网格线，区分黑色和白色线条。最后调用drow_the_scale函数绘制比例尺，并保存图像为PNG格式，分辨率为600dpi。


def drow_the_scale(y,x,text,length = 1.5,lw = 5):
    #画比例尺函数
    # y代表比例尺所在纬度
    # x代表比例尺开始的经度
    # text代表比例尺最后刻度值
    # length代表比例尺的长度，单位为多少个经度
    # lw代表比例尺的宽度
    step = length/5#计算步长，画五格
    #画黑白线五条
    plt.hlines(y=y,xmin=x,xmax=x + step,colors="black", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step,xmax=x + step*2,colors="white", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step*2,xmax=x + step*3,colors="black", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step*3,xmax=x + step*4,colors="white", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step*4,xmax=x + step*5,colors="black", ls="-", lw=lw)
    #画长刻度两个
    plt.vlines(x = x, ymin = y - (lw/100) *3, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + length, ymin = y - (lw/100) *3, ymax = y + lw/100, colors="black", ls="-", lw=1)
    #画段刻度四个
    plt.vlines(x = x + step, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + step*2, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + step*3, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + step*4, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    #写字，0，500，km
    plt.text(x,y - (lw/100) *7,'0',horizontalalignment = 'center')
    plt.text(x + length,y - (lw/100) *7,text,horizontalalignment = 'center')
    plt.text(x + length/2,y + (lw/100)*2,'km',horizontalalignment = 'center')
def drowscale(extent,scale_y,scale_x,scale_text,step = 5,lw = 10,scale_length = 1.5,scale_lw = 5):
    # 画地图黑白边框和比例尺
    # extent：表示四周经纬度[west, east, south, north]
    # scale_y,scale_x,scale_text：代表比例尺的位置，纬度，经度，刻度值
    # step：表示步长，一格代表几个经纬度
    # lw：代表边框宽度
    # scale_length：代表比例尺长度（单位为经度例如1.5个经度）
    # scale_lw：代表比例尺宽度
    for y in [extent[2],extent[3]] :#画上下两边框
        xmin = extent[0]
        while (xmin < extent[1]):
            plt.hlines(y=y,xmin=xmin,xmax=xmin+step,colors="white", ls="-", lw=lw)
            xmin = xmin+step*2
        xmin = extent[0]+step
        while (xmin < extent[1]):
            plt.hlines(y=y,xmin=xmin,xmax=xmin+step,colors="black", ls="-", lw=lw)
            xmin = xmin+step*2
    for x in [extent[0],extent[1]] :#画左右两边狂
        ymin = extent[2]
        while (ymin < extent[3]):
            plt.vlines(x = x, ymin = ymin, ymax = ymin+step, colors="black", ls="-", lw=lw)
            ymin = ymin+step*2
        ymin = extent[2]+step
        while (ymin < extent[3]):
            plt.vlines(x = x, ymin = ymin, ymax = ymin+step, colors="white", ls="-", lw=lw)
            ymin = ymin+step*2
    drow_the_scale(scale_y,scale_x,scale_text)#画比例尺