python使用cartopy画地图

最新推荐文章于 2025-03-02 16:15:41 发布

萝卜一下

最新推荐文章于 2025-03-02 16:15:41 发布

阅读量2.9k

点赞数 4

文章标签： python 开发语言

本文链接：https://blog.csdn.net/weixin_41796280/article/details/130759462

版权

1、环境准备

创建conda环境

conda create -n map_cartopy python=3.9

激活环境

conda activate map_cartopy

安装cartopy，使用conda install 进行安装，一步搞定，pip安装需要先下载依赖库，conda安装直接将依赖库包含在内

conda install cartopy

下载其他依赖库

xarray # 读取气象数据 .nc文件
matplotlib
pykrige # 调用kriging插值算法
frykit # 大神写的依赖cartopy包画中国的包
openpyxl
cnmaps
scipy
basemap #需要离线下载whl进行安装 目前已经不再更新 可以不安装

这里安装的frykit程序包是大佬自己写的一个画中国地图的程序包，我们在此包的基础上进行画图，修改其中一些函数进行封装。
basemap需要下载whl轮子进行离线安装，首先去官网下载对应的basemap和pyproj，官网地址为https://www.lfd.uci.edu/~gohlke/pythonlibs/，该库目前已经不更新，本教程不使用basemap进行画图，但会在后续代码中用到。

2、画中国地图

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import frykit.plot as fplt
import frykit.shp as fshp

# from scipy.interpolate import Rbf


#设置画图字体的大小
plt.rcParams.update({'font.size':20})
#解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
plt.rc("font", family='Microsoft YaHei')

# 读取shp文件记录.
names = fshp.get_cn_province_names()
provinces = fshp.get_cn_shp(level='省')

# 设置地图范围和刻度.
extents1 = [74, 136, 14, 56]
extents2 = [105, 122, 2, 25]
xticks = np.arange(-180, 181, 10)
yticks = np.arange(-90, 91, 10)

# 设置投影.
map_crs = ccrs.LambertConformal(
    central_longitude=105,
    standard_parallels=(25, 47)
)
data_crs = ccrs.PlateCarree()

# 设置刻度风格.
plt.rc('xtick.major', size=8, width=0.9)
plt.rc('ytick.major', size=8, width=0.9)
plt.rc('xtick', labelsize=8, top=False, labeltop=False)
plt.rc('ytick', labelsize=8, right=False, labelright=False)

# 准备主地图.
fig = plt.figure(figsize=(10, 6))
##################
fig.canvas.draw()# 运行程序必须加上这句，为了将画布进行缓存，否则运行出错，去掉这句在shell调试阶段仍能运行。因为调试阶段都是缓存的数据
###################
ax1 = fig.add_subplot(projection=map_crs)
fplt.set_extent_and_ticks(
    ax1, extents=extents1,
    xticks=xticks, yticks=yticks
)
ax1.gridlines(
    xlocs=xticks, ylocs=yticks,
    lw=0.5, ls='--', color='gray'
)

# 添加要素.
# 画海，陆地，河流，湖泊
# ax1.add_feature(cfeature.OCEAN.with_scale('50m'))
# ax1.add_feature(cfeature.LAND.with_scale('50m'))
# ax1.add_feature(cfeature.RIVERS.with_scale('50m'))
# ax1.add_feature(cfeature.LAKES.with_scale('50m'))
ax1.add_feature(cfeature.LAND.with_scale('50m'), fc='floralwhite')
ax1.add_feature(cfeature.OCEAN.with_scale('50m'), fc='skyblue')
fplt.add_cn_province(ax1, lw=0.3)
fplt.add_nine_line(ax1, lw=0.5)

# 添加指北针和比例尺.
fplt.add_compass(ax1, 0.92, 0.85, size=15, style='star')
scale = fplt.add_map_scale(ax1, 0.05, 0.1, length=1000)
scale.set_xticks([0, 500, 1000])

# 设置标题.
ax1.set_title(
    '中国地图', y=1.01,
    fontsize='large', weight='bold'
)

# 准备小地图.
ax2 = fig.add_subplot(projection=map_crs)
ax2.set_extent(extents2, crs=data_crs)
fplt.move_axes_to_corner(ax2, ax1, shrink=0.3, loc='bottom right')
ax2.gridlines(
    xlocs=xticks, ylocs=yticks,
    lw=0.5, ls='--', color='gray'
)

# 添加要素.
# ax2.add_feature(cfeature.OCEAN.with_scale('50m'))
# ax2.add_feature(cfeature.LAND.with_scale('50m'))
# ax2.add_feature(cfeature.RIVERS.with_scale('50m'))
# ax2.add_feature(cfeature.LAKES.with_scale('50m'))
ax2.add_feature(cfeature.LAND.with_scale('50m'), fc='floralwhite')
ax2.add_feature(cfeature.OCEAN.with_scale('50m'), fc='skyblue')
fplt.add_nine_line(ax2, lw=0.5)
fplt.add_cn_province(ax2, lw=0.3)
fplt.add_map_scale(ax2, 0.4, 0.15, length=500)

# 简化名称.
for i, name in enumerate(names):
    name = fshp.simplify_province_name(name)
    if name == '香港' or name == '澳门':
        name = ''
    names[i] = name

# 添加省名.
for name, province in zip(names, provinces):
    point = province.representative_point()
    ax1.text(
        point.x, point.y, name,
        ha='center', va='center',
        fontsize=5,
        transform=data_crs
    )
# plt.show()
# 保存图片.
fig.savefig('./image/cn_map.png', dpi=300, bbox_inches='tight')
plt.close(fig)

结果展示
在这里插入图片描述

3、利用克里金插值算法画全国AQI分布图

利用338城市小时数据进行画图，先画一个散点图看一下大概分布情况。

import pandas as pd
import matplotlib.pyplot as plt

# 散点图测试
df = pd.read_excel(r'./example/test_pm10.xlsx')

x = df['经度'].values
y = df['纬度'].values
z = df['AQI'].values
fig = plt.figure()
plt.scatter(x, y, c=z, vmax=160, vmin=15,cmap='turbo')
plt.show()
fig.savefig('./image/contourf_AQI_scatter.png', dpi=300, bbox_inches='tight')
plt.close(fig)

结果如下：
在这里插入图片描述
接下来导入数据，进行克里金插值，再将插值数据映射到中国地图上。

import numpy as np
import xarray as xr
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import frykit.plot as fplt
import frykit.shp as fshp
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import pandas as pd
# from scipy.interpolate import Rbf
from pykrige.ok import OrdinaryKriging
import copy

#设置画图字体的大小
plt.rcParams.update({'font.size':20})
#解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
plt.rc("font", family='Microsoft YaHei')

# 读取shp文件记录.
names = fshp.get_cn_province_names()
provinces = fshp.get_cn_shp(level='省')

#############################################################
# 读取数据.
df = pd.read_excel(r'./example/test_pm10.xlsx')

Krin = OrdinaryKriging(df['经度'], df['纬度'], df['AQI'], variogram_model='linear')

# lon = df['经度'].values
# lat = df['纬度'].values
# temp = df['temp'].values

# #对经纬度进行插值，这里作者讲经纬度插值成了0.05°*0.05°的网格
# olon = np.linspace(70,140,100)
# olat = np.linspace(15,60,100)
olon = np.arange(74,136,1.0)
olat = np.arange(14,56,1.0)
data, ssl = Krin.execute('grid', olon, olat)
olon,olat = np.meshgrid(olon,olat)

# func = Rbf(lon,lat,temp,function='cubic')
# #对var进行插值，得到格点var_new1
# temp_new1 = func(olon,olat)
# # #筛选数据清洗
# temp_new1[temp_new1 < 15] = 15
# temp_new1[temp_new1 > 160] = 160
# X = olon
# Y = olat
# Z = temp_new1 
X = olon
Y = olat
Z = data

###########################################################

# 设置地图范围和刻度.
extents1 = [74, 136, 14, 56]
extents2 = [105, 122, 2, 25]
xticks = np.arange(-180, 181, 10)
yticks = np.arange(-90, 91, 10)

# 设置投影.
map_crs = ccrs.LambertConformal(
    central_longitude=105,
    standard_parallels=(25, 47)
)
data_crs = ccrs.PlateCarree()

# 设置刻度风格.
plt.rc('xtick.major', size=8, width=0.9)
plt.rc('ytick.major', size=8, width=0.9)
plt.rc('xtick', labelsize=8, top=False, labeltop=False)
plt.rc('ytick', labelsize=8, right=False, labelright=False)

# 准备主地图.
fig = plt.figure(figsize=(10, 6))
##################
fig.canvas.draw()# 运行程序必须加上这句，为了将画布进行缓存，否则运行出错，去掉这句在shell调试阶段仍能运行。因为调试阶段都是缓存的数据
###################
ax1 = fig.add_subplot(projection=map_crs)
fplt.set_extent_and_ticks(
    ax1, extents=extents1,
    xticks=xticks, yticks=yticks
)
ax1.gridlines(
    xlocs=xticks, ylocs=yticks,
    lw=0.5, ls='--', color='gray'
)

# 添加要素.
ax1.add_feature(cfeature.LAND.with_scale('50m'), fc='floralwhite')
ax1.add_feature(cfeature.OCEAN.with_scale('50m'), fc='skyblue')
fplt.add_cn_province(ax1, lw=0.3)
fplt.add_nine_line(ax1, lw=0.5)

# 绘制填色图.
levels = np.linspace(10, 160, 50)#(0,300,50)

cf = ax1.contourf(
    X, Y, Z, levels, cmap='turbo', extend='both',
    transform=data_crs, transform_first=True
)
fplt.clip_by_cn_border(cf)


# 绘制colorbar.
# cbar = fig.colorbar(
#     cf, ax=ax1, orientation='vertical',
#     shrink=0.6, pad=0.1, aspect=20,
#     ticks=cticks, extendfrac=0
# )
# cbar.ax.tick_params(length=4, labelsize=8)

# 添加一个新绘图，标注在主图的哪个位置 这里修改colorbar的刻度
# cticks = np.linspace(0, 300, 15, dtype=int)
# ax3 = fig.add_axes([ax1.get_position().x1-0.03,ax1.get_position().y0,0.02,ax1.get_position().height*0.25])
# cmap4 = cm.get_cmap('turbo', len(cticks)-1)
# norm4 = mcolors.BoundaryNorm(cticks, len(cticks)-1)
# im4 = cm.ScalarMappable(norm=norm4, cmap=cmap4)
# # 使用BoundaryNorm时,colorbar会自动按bins标出刻度.
# cbar4 = fig.colorbar(im4, cax=ax3, aspect=10)
# fplt.move_axes_to_corner(ax3, ax1, shrink=0.3, loc='bottom right')

# 另一种colorbar方式
ax3 = fig.add_axes([ax1.get_position().x1-0.03,ax1.get_position().y0,0.02,ax1.get_position().height*0.25])
cmap1 = copy.copy(cm.turbo)
norm1 = mcolors.Normalize(vmin=10, vmax=160)
im1 = cm.ScalarMappable(norm=norm1, cmap=cmap1)
cbar1 = fig.colorbar(im1, cax=ax3, ticks=np.linspace(10, 160, 6))#(0,300,7)
fplt.move_axes_to_corner(ax3, ax1, shrink=0.3, loc='bottom right')


# 添加指北针和比例尺.
fplt.add_compass(ax1, 0.92, 0.85, size=15, style='star')
# scale = fplt.add_map_scale(ax1, 0.95, 1, length=1000)
# scale.set_xticks([0, 500, 1000])

# 设置标题.
ax1.set_title(
    '338城市5月2日10点AQI分布图', y=1.01,
    fontsize=10, weight='bold'
)

# 准备小地图.
ax2 = fig.add_subplot(projection=map_crs)
ax2.set_extent(extents2, crs=data_crs)
fplt.move_axes_to_corner(ax2, ax1, shrink=0.3, loc='bottom left')
ax2.gridlines(
    xlocs=xticks, ylocs=yticks,
    lw=0.5, ls='--', color='gray'
)

# 添加要素.
ax2.add_feature(cfeature.LAND.with_scale('50m'), fc='floralwhite')
ax2.add_feature(cfeature.OCEAN.with_scale('50m'), fc='skyblue')
fplt.add_nine_line(ax2, lw=0.5)
fplt.add_cn_province(ax2, lw=0.3)
# fplt.add_map_scale(ax2, 0.4, 0.15, length=500)

# 绘制填色图.
cf = ax2.contourf(
    X, Y, Z, levels, cmap='turbo', extend='both',
    transform=data_crs, transform_first=True
)
fplt.clip_by_cn_border(cf)

# 简化名称.
for i, name in enumerate(names):
    name = fshp.simplify_province_name(name)
    if name == '香港' or name == '澳门':
        name = ''
    names[i] = name

# 添加省名.
for name, province in zip(names, provinces):
    point = province.representative_point()
    ax1.text(
        point.x, point.y, name,
        ha='center', va='center',
        fontsize=5,
        transform=data_crs
    )

# 保存图片.
fig.savefig('./image/contourf_4.png', dpi=300, bbox_inches='tight')
plt.close(fig)

看一下结果：
在这里插入图片描述

3、利用克里金插值算法画京津冀AQI分布图

其中设计多省份切割、添加，城市切割、添加。

import numpy as np
import xarray as xr
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import frykit.plot as fplt
import frykit.shp as fshp
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import pandas as pd
# from scipy.interpolate import Rbf
from pykrige.ok import OrdinaryKriging
import copy
from cartopy.mpl.geoaxes import GeoAxes
from weakref import WeakValueDictionary, WeakKeyDictionary
from cartopy.mpl.feature_artist import _GeomKey
from matplotlib.collections import PathCollection
import shapely.geometry as sgeom



############################功能函数####################################
# 当polygon的引用计数为零时, 弱引用会自动清理缓存.
_key_to_polygon_cache = WeakValueDictionary()
_key_to_transformed_cache = WeakKeyDictionary()

_transform_geometry = None

def enable_fast_transform():
    '''启用快速坐标变换. 可能在地图边界产生错误的连线.'''
    global _transform_geometry
    _transform_geometry = fshp.transform_geometry
    _key_to_transformed_cache.clear()

def disable_fast_transform():
    '''关闭快速坐标变换. 变换结果更正确, 但速度可能较慢.'''
    global _transform_geometry
    def _transform_geometry(geom, crs_from, crs_to):
        return crs_to.project_geometry(geom, crs_from)
    _key_to_transformed_cache.clear()

enable_fast_transform()

def _cached_transform(polygon, crs_from, crs_to):
    '''调用_transform_geometry并缓存结果.'''
    if crs_from == crs_to:
        return polygon

    key = _GeomKey(polygon)
    _key_to_polygon_cache.setdefault(key, polygon)
    mapping = _key_to_transformed_cache.setdefault(key, {})
    value = mapping.get(crs_to)
    if value is None:
        value = _transform_geometry(polygon, crs_from, crs_to)
        mapping[crs_to] = value

    return value

def add_polygons(ax, polygons, crs=None, **kwargs):
    '''
    将一组多边形添加到ax上.

    Parameters
    ----------
    ax : Axes or GeoAxes
        目标Axes.

    polygons : list of Polygon or list of MultiPolygon
        多边形构成的列表.

    crs : CRS, optional
        多边形所处的坐标系. 默认为ccrs.PlateCarree().
        当ax是Axes时该参数无效, 不会对多边形进行坐标变换.

    **kwargs
        创建PathCollection时的关键字参数.
        例如facecolors, edgecolors, cmap, norm和array等.

    Returns
    -------
    pc : PathCollection
        代表一组多边形的集合对象.
    '''
    kwargs.setdefault('zorder', 1.5)
    array = kwargs.get('array', None)
    if array is not None and len(array) != len(polygons):
        raise ValueError('array的长度与polygons不匹配')

    if crs is None:
        crs = ccrs.PlateCarree()

    # GeoAxes会对多边形做坐标变换.
    trans = ax.transData
    to = fshp.polygon_to_path
    if isinstance(ax, GeoAxes):
        trans = ax.projection._as_mpl_transform(ax)
        to = lambda x: fshp.polygon_to_path(
            _cached_transform(x, crs, ax.projection)
        )

    # PathCollection比PathPatch更快
    paths = [to(polygon) for polygon in polygons]
    pc = PathCollection(paths, transform=trans, **kwargs)
    ax.add_collection(pc)

    return pc

def add_polygon(ax, polygon, crs=None, **kwargs):
    '''
    将一个多边形添加到ax上.

    Parameters
    ----------
    ax : Axes or GeoAxes
        目标Axes.

    polygons : list of Polygon or list of MultiPolygon
        多边形构成的列表.

    crs : CRS, optional
        多边形所处的坐标系. 默认为ccrs.PlateCarree().
        当ax是Axes时该参数无效, 不会对多边形进行坐标变换.

    **kwargs
        创建PathCollection时的关键字参数.
        例如facecolors, edgecolors, cmap, norm和array等.

    Returns
    -------
    pc : PathCollection
        只含一个多边形的集合对象.
    '''
    return add_polygons(ax, [polygon], crs, **kwargs)

def _get_boundary(ax):
    '''将GeoAxes.patch转为data坐标系下的多边形.'''
    patch = ax.patch
    ax.draw_artist(patch)  # 决定patch的形状.
    trans = patch.get_transform() - ax.transData
    # get_path比get_verts更可靠.
    path = patch.get_path().transformed(trans)
    boundary = sgeom.Polygon(path.vertices)

    return boundary

# 缓存常用数据.
_data_cache = {}

def _set_add_cn_kwargs(kwargs):
    '''初始化add_cn_xxx函数的参数.'''
    if not any(kw in kwargs for kw in ['facecolor', 'facecolors', 'fc']):
        kwargs['facecolors'] = 'none'
    if not any(kw in kwargs for kw in ['edgecolor', 'edgecolors', 'ec']):
        kwargs['edgecolors'] = 'black'

def _get_cn_city():
    '''获取中国市界并缓存结果.'''
    mapping = _data_cache.get('city')
    if mapping is None:
        names = fshp.get_cn_city_names()
        citys = fshp.get_cn_shp(level='市')
        mapping = dict(zip(names, citys))
        _data_cache['city'] = mapping
    return mapping

def add_cn_city(ax, name=None, **kwargs):
    '''
    将中国市界添加到ax上.

    Parameters
    ----------
    ax : Axes or GeoAxes
        目标Axes.

    name : str or list of str, optional
        市名, 可以是字符串或字符串构成的列表.
        默认为None, 表示添加所有市.

    **kwargs
        创建PathCollection时的关键字参数.
        例如facecolors, edgecolors, cmap, norm和array等.

    Returns
    -------
    pc : PathCollection
        代表市界的集合对象.
    '''
    mapping = _get_cn_city()
    if name is None:
        city = mapping.values()
    elif isinstance(name, str):
        city = [mapping[name]]
    else:
        city = [mapping[n] for n in name]
    _set_add_cn_kwargs(kwargs)
    pc = add_polygons(ax, city, **kwargs)

    return pc

def _get_province_city(name):
    '''获取中国市界并缓存结果.'''
    mapping = _data_cache.get('city')
    if mapping is None:
        # names = fshp.get_cn_city_names()
        citys = fshp.get_cn_shp(level='市',province = name)
        mapping = citys
        _data_cache['city'] = mapping
    return mapping

def add_province_city(ax, name=None, **kwargs):
    '''
    将中国某省的市界添加到ax上.

    Parameters
    ----------
    ax : Axes or GeoAxes
        目标Axes.

    name : str or list of str, optional
        省名, 字符串

    **kwargs
        创建PathCollection时的关键字参数.
        例如facecolors, edgecolors, cmap, norm和array等.

    Returns
    -------
    pc : PathCollection
        代表市界的集合对象.
    '''
    city = _get_province_city(name)
    _set_add_cn_kwargs(kwargs)
    pc = add_polygons(ax, city, **kwargs)

    return pc

################################################################


#设置画图字体的大小
plt.rcParams.update({'font.size':20})
#解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
#解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
plt.rc("font", family='Microsoft YaHei')

#############################################################
# 读取数据.
df = pd.read_excel(r'./example/河北省.xlsx')

Krin = OrdinaryKriging(df['经度'], df['纬度'], df['AQI'], variogram_model='gaussian')

# lon = df['经度'].values
# lat = df['纬度'].values
# temp = df['AQI'].values

#对经纬度进行插值，这里作者讲经纬度插值成了0.05°*0.05°的网格
# olon = np.linspace(70,140,100)
# olat = np.linspace(15,60,100)
olon = np.arange(113,120,0.2)
olat = np.arange(36,43,0.2)
data, ssl = Krin.execute('grid', olon, olat)
olon,olat = np.meshgrid(olon,olat)

# func = Rbf(lon,lat,temp,function='cubic')
# #对var进行插值，得到格点var_new1
# temp_new1 = func(olon,olat)
# #筛选数据清洗
# temp_new1[temp_new1 < 15] = 15
# temp_new1[temp_new1 > 160] = 160
# X = olon
# Y = olat
# Z = temp_new1 
X = olon
Y = olat
Z = data
##################################################


# 设置地图范围和刻度.
extents1 = [113, 120, 36, 43]

xticks = np.arange(-180, 181, 2)
yticks = np.arange(-90, 91, 2)

# 设置投影.
map_crs = ccrs.LambertConformal(
    central_longitude=105,
    standard_parallels=(25, 47)
)
data_crs = ccrs.PlateCarree()

# 设置刻度风格.
plt.rc('xtick.major', size=5, width=0.6)
plt.rc('ytick.major', size=5, width=0.6)
plt.rc('xtick', labelsize=8, top=False, labeltop=False)
plt.rc('ytick', labelsize=8, right=False, labelright=False)

# 准备主地图.
fig = plt.figure(figsize=(10, 6))
##################
fig.canvas.draw()
###################
ax1 = fig.add_subplot(projection=map_crs)
fplt.set_extent_and_ticks(
    ax1, extents=extents1,
    xticks=xticks, yticks=yticks
)
ax1.gridlines(
    xlocs=xticks, ylocs=yticks,
    lw=0.5, ls='--', color='gray'
)

# 添加要素.
ax1.add_feature(cfeature.LAND.with_scale('50m'), fc='floralwhite')
ax1.add_feature(cfeature.OCEAN.with_scale('50m'), fc='skyblue')
# hebei_names = ['石家庄市','唐山市','保定市','秦皇岛市','邢台市','衡水市','廊坊市','张家口市','沧州市','邯郸市','承德市']
# add_cn_city(ax1,lw=0.3)
add_province_city(ax1,name='河北省',lw=0.3)
name_BT = ['北京市','天津市','内蒙古自治区','山西省','河南省','山东省','辽宁省']
fplt.add_cn_province(ax1,name = name_BT,lw=0.3)

# 绘制填色图.
levels = np.linspace(15, 160, 50)

cf = ax1.contourf(
    X, Y, Z, levels, cmap='turbo', extend='both',
    transform=data_crs, transform_first=True
)
hebei = fshp.get_cn_shp(level='省',province='河北省')
beijing = fshp.get_cn_shp(level='省',province='北京市')
tianjin = fshp.get_cn_shp(level='省',province='天津市')
fplt.clip_by_polygons(cf,[hebei,beijing,tianjin])

#################################
# colorbar
# 另一种colorbar方式
ax3 = fig.add_axes([ax1.get_position().x1-0.03,ax1.get_position().y0,0.02,ax1.get_position().height*0.25])
cmap1 = copy.copy(cm.turbo)
norm1 = mcolors.Normalize(vmin=0, vmax=300)
im1 = cm.ScalarMappable(norm=norm1, cmap=cmap1)
cbar1 = fig.colorbar(im1, cax=ax3, ticks=np.linspace(0, 300, 7))
fplt.move_axes_to_corner(ax3, ax1, shrink=0.3, loc='bottom right')
#######################################################
# 添加标题
# 设置标题.
ax1.set_title(
    '河北省5月2日10点AQI分布图', y=1.01,
    fontsize=10, weight='bold'
)
###########################################
hebei_names_values = {
    '北京':[116.724762,39.904896],
    '天津':[117.201324,39.085068],
    '石家庄':[114.514778,38.04232],
    '唐山':[118.180222,39.630451],
    '保定':[115.464592,38.873959],
    '秦皇岛':[119.517807,39.88932],
    '邯郸':[114.539192,36.625561],
    '邢台':[114.504913,37.070572],
    '张家口':[114.885963,40.769062],
    '承德':[117.96286,40.952869],
    '沧州':[116.83886,38.304459],
    '衡水':[115.66922,37.739033],
    '廊坊':[116.68351,39.53791]
}
# 写出城市名称信息 
for key,value in hebei_names_values.items():
        # ax1.scatter(value[0] , value[1] , marker='.' , s=65 , color = "k" , zorder = 3,transform=data_crs)
        # ax1.text(value[0]-0.09 , value[1]+0.05 , key , fontsize = 9 , color = "k",transform=data_crs)
        ax1.text(value[0] , value[1] , key , ha='center', va='center', fontsize = 6 , color = "k",transform=data_crs)

# plt.show()
# 保存图片.
fig.savefig('./image/hebei.png', dpi=300, bbox_inches='tight')
plt.close(fig)