使用 Python 绘制历年中国省市区县地图(小地图版本+长版)

一、背景与原理

之前我们使用 R 语言的 sf + ggplot2 体系绘制历年中国省市区县地图(小地图版本 + 长版),本讲义将同样的功能用 Python 实现。Python 版本使用 geopandas + matplotlib 体系,提供了完整的模块化绘图函数。

附件中提供了 1949~2021 年每年的省市区县地图数据,包含 mini 和 long 两套:

mini 版本(小地图版本) long 版本(长版)

R 与 Python 的包映射

R 包 Python 等价 功能
sf geopandas 矢量数据读写与操作
ggplot2 + geom_sf matplotlib 地图绘制
ggspatial(比例尺、指北针) matplotlib(手动实现) 地图装饰
raster rasterio 栅格数据读取
readxl pandas(read_excel) Excel 数据读取
dplyr pandas 数据操作
scico matplotlib colormaps 科学配色

二、使用 reticulate 创建与管理 Python 虚拟环境

2.1 安装 reticulate

options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
install.packages("reticulate")

2.2 虚拟环境初始化原理

reticulate 会在 R 会话中绑定一个 Python 解释器,绑定后不可更改。因此我们必须在任何 {python} chunk 运行之前,通过 Sys.setenv(RETICULATE_PYTHON = ...) 环境变量锁定 Python 路径。setup chunk 中的代码已完成这一操作。

2.3 安装 Python 包

py_pkgs <- c("numpy", "pandas", "geopandas", "matplotlib",
"rasterio", "shapely", "openpyxl", "affine")
installed <- py_list_packages(".venv")$package
need_install <- setdiff(py_pkgs, installed)
if (length(need_install) > 0) {
virtualenv_install(".venv", packages = need_install)
cat("已安装:", paste(need_install, collapse = ", "), "\n")
} else {
cat("所有包已安装完毕。\n")
}

2.4 验证激活状态

py_config()

2.5 查看已安装的包

py_list_packages(".venv")

2.6 虚拟环境管理常用命令

# 删除虚拟环境(如需重建)
virtualenv_remove(".venv")

# 重建虚拟环境
virtualenv_create(".venv")
virtualenv_install(".venv", packages = c("numpy", "pandas", "geopandas", ...))

三、环境配置与数据准备

3.1 导入 Python 包

import os
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.font_manager as fm
from matplotlib.patches import Rectangle
from shapely.geometry import box, Point
from shapely.ops import unary_union
import warnings
warnings.filterwarnings("ignore")

3.2 全局配置

# 中国地图常用的 Albers 等面积投影
MY_CRS = "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 小地图坐标范围(MY_CRS 投影下的坐标)
SMALL_BBOX = {"xmin": 120000, "xmax": 1766004.1,
"ymax": 2557786.0, "ymin": 320000}

# 小地图缩放和平移参数
SMALL_SCALE = 0.5
SMALL_OFFSET_X = 2100000
SMALL_OFFSET_Y = 1665139

# ---- 字体配置:使用附件中的 LXGWWenKai 字体 ----
_font_path = os.path.join(os.getcwd(), "LXGWWenKai-Regular.ttf")
fm.fontManager.addfont(_font_path)
_font_prop = fm.FontProperties(fname=_font_path)
CN_FONT = _font_prop.get_name()

plt.rcParams["font.family"] = CN_FONT
plt.rcParams["axes.unicode_minus"] = False
plt.rcParams["font.size"] = 10

# 配色方案
LINE_COLORS = {
"九段线": "#A29AC4",
"海岸线": "#0055AA",
"小地图框格": "black",
"省份": "black",
}
LINE_WIDTHS = {
"九段线": 0.6,
"海岸线": 0.3,
"小地图框格": 0.3,
"省份": 0.3,
}

3.3 比例尺与装饰工具函数

def add_scale_bar(ax, provmap, location="bl",
height_scale=0.016, width_scale=0.2):
"""在地图左下角添加比例尺。"""
bounds = provmap.total_bounds
map_width_m = bounds[2] - bounds[0]
target_km = map_width_m * width_scale / 1000
nice_values = [500, 1000, 2000, 3000, 5000, 10000, 20000, 50000, 100000]
bar_km = min(nice_values, key=lambda x: abs(x - target_km))
bar_m = bar_km * 1000

x0_frac = 0.08 if location == "bl" else 0.62
xlim = ax.get_xlim()
ylim = ax.get_ylim()
ax_width = xlim[1] - xlim[0]
x0_data = xlim[0] + x0_frac * ax_width
y0_data = ylim[0] + 0.04 * (ylim[1] - ylim[0])
bar_height_m = map_width_m * height_scale
n_seg = 4
seg_w = bar_m / n_seg
for i in range(n_seg):
color = "black" if i % 2 == 0 else "white"
rect = Rectangle(
(x0_data + i * seg_w, y0_data), seg_w, bar_height_m,
facecolor=color, edgecolor="black", linewidth=0.5,
clip_on=False, zorder=10,
)
ax.add_patch(rect)
ax.text(x0_data + bar_m / 2, y0_data - bar_height_m * 0.3,
f"{bar_km:,} km", fontsize=7, ha="center", va="top",
color="black", zorder=10)

def add_north_arrow(ax):
"""在地图右上角添加指北针。"""
ax.annotate("N", xy=(0.98, 0.95), xycoords="axes fraction",
fontsize=10, fontweight="bold", ha="center", va="bottom")
ax.annotate("", xy=(0.98, 0.95), xycoords="axes fraction",
xytext=(0.98, 0.88), textcoords="axes fraction",
arrowprops=dict(arrowstyle="->", color="black", lw=1.5))

def draw_lines(ax, provlinemap):
"""绘制线条元素。"""
for _, row in provlinemap.iterrows():
cls = row["class"]
color = LINE_COLORS.get(cls, "gray")
lw = LINE_WIDTHS.get(cls, 0.3)
geom = row["geometry"]
if geom.geom_type == "MultiLineString":
for line in geom.geoms:
xs, ys = line.xy
ax.plot(xs, ys, color=color, linewidth=lw, zorder=3)
elif geom.geom_type == "LineString":
xs, ys = geom.xy
ax.plot(xs, ys, color=color, linewidth=lw, zorder=3)

def draw_province_labels(ax, provmap, centroids):
"""绘制省名标注。"""
for i, row in provmap.iterrows():
cx, cy = centroids.iloc[i].x, centroids.iloc[i].y
ax.text(cx, cy, row["省"], fontsize=5, color="gray",
ha="center", va="center", zorder=4)

def add_vertical_colorbar(ax, fig, cmap, norm, label=None,
ticks=None, ticklabels=None,
n_segments=4):
"""在地图左下角、比例尺上方绘制竖向分段色条。"""
xlim = ax.get_xlim()
ylim = ax.get_ylim()
ax_width = xlim[1] - xlim[0]
ax_height = ylim[1] - ylim[0]

cb_width_m = ax_width * 0.0075
cb_height_m = ax_height * 0.15
seg_height_m = cb_height_m / n_segments

x0_frac = 0.08
x0_data = xlim[0] + x0_frac * ax_width
scale_y0 = ylim[0] + 0.04 * ax_height
bar_height_m = ax_width * 0.016
cb_bottom = scale_y0 + bar_height_m + ax_height * 0.025

vmin, vmax = norm.vmin, norm.vmax
for i in range(n_segments):
frac_lo = i / n_segments
frac_hi = (i + 1) / n_segments
color = cmap((frac_lo + frac_hi) / 2)
y_bottom = cb_bottom + i * seg_height_m
rect = Rectangle(
(x0_data, y_bottom), cb_width_m, seg_height_m,
facecolor=color, edgecolor="black", linewidth=0.5,
clip_on=False, zorder=10,
)
ax.add_patch(rect)

if ticks is not None and ticklabels is not None:
for tick_val, tick_lbl in zip(ticks, ticklabels):
frac = (tick_val - vmin) / (vmax - vmin) if vmax > vmin else 0
frac = max(0, min(1, frac))
y_pos = cb_bottom + frac * cb_height_m
ax.plot(
[x0_data + cb_width_m, x0_data + cb_width_m + cb_width_m * 0.3],
[y_pos, y_pos],
color="black", linewidth=0.5, clip_on=False, zorder=10,
)
ax.text(
x0_data + cb_width_m + cb_width_m * 0.4, y_pos,
tick_lbl, fontsize=6, ha="left", va="center",
color="black", zorder=10,
)

if label:
ax.text(
x0_data + cb_width_m / 2, cb_bottom + cb_height_m + ax_height * 0.01,
label, fontsize=7, ha="center", va="bottom",
color="black", zorder=10,
)

四、小地图版本

4.1 数据读取

读取小地图版本的中国省级地图数据:

# 读取省级多边形
provmap = gpd.read_file(
"provmapdata/minishp/chinaprov2019mini/chinaprov2019mini.shp",
encoding="utf-8"
)
provmap = provmap.dropna(subset=["省代码"]).reset_index(drop=True)
print(f"省级多边形: {len(provmap)} 个区域")
print(provmap.head())

读取线条数据(九段线、海岸线、小地图框格):

# 读取线条
provlinemap = gpd.read_file(
"provmapdata/minishp/chinaprov2019mini/chinaprov2019mini_line.shp",
encoding="utf-8"
)
keep_classes = ["九段线", "海岸线", "小地图框格"]
provlinemap = provlinemap[
provlinemap["class"].isin(keep_classes)
][["class", "geometry"]].reset_index(drop=True)
print(f"线条: {len(provlinemap)} 条")
print(provlinemap["class"].value_counts())

4.2 在地图上添加散点图

读取经纬度数据并转换为 GeoDataFrame:

# 读取散点数据
pointdf = pd.read_excel("瞪羚、独角兽、创新型企业经纬度数据.xlsx")
pointdf = pointdf.dropna(subset=["经度"]).reset_index(drop=True)
print(f"散点数据: {len(pointdf)} 条")
print(pointdf.head())
# 转换为 GeoDataFrame(WGS84 → MY_CRS)
pointdfsf = gpd.GeoDataFrame(
pointdf,
geometry=gpd.points_from_xy(pointdf["经度"], pointdf["纬度"]),
crs="EPSG:4326",
)
pointdfsf = pointdfsf.to_crs(MY_CRS)
print(pointdfsf.crs)

4.3 小地图的散点处理原理

小地图框格显示的是南海诸岛放大图。对于散点数据,需要将落在小地图范围内的点平移到框格位置。处理流程如下:

  1. 用 st_bbox 定义小地图在主图中的范围
  2. 用 st_intersection 提取落在该范围内的散点
  3. 对这些散点执行:先以范围左下角为原点缩小到 50%,再平移到小地图框格位置
  4. 将平移后的散点与原始散点合并
# 步骤 1:定义小地图范围
small_bbox = box(
SMALL_BBOX["xmin"], SMALL_BBOX["ymin"],
SMALL_BBOX["xmax"], SMALL_BBOX["ymax"]
)

# 步骤 2:提取小地图范围内的散点
pointdfsf_small = pointdfsf[
pointdfsf.geometry.intersects(small_bbox)
].copy()
print(f"小地图范围内的散点: {len(pointdfsf_small)} 个")

# 步骤 3:平移到小地图位置
# 先将坐标原点移到 (xmin, ymin),再缩放 0.5,最后平移到目标位置
gdf_small = pointdfsf_small.copy()
gdf_small["geometry"] = gdf_small.geometry.translate(
xoff=-SMALL_BBOX["xmin"],
yoff=-SMALL_BBOX["ymin"],
)
gdf_small["geometry"] = gdf_small.geometry.scale(
xfact=SMALL_SCALE, yfact=SMALL_SCALE, origin=(0, 0),
)
gdf_small["geometry"] = gdf_small.geometry.translate(
xoff=SMALL_OFFSET_X + SMALL_BBOX["xmin"] * SMALL_SCALE,
yoff=SMALL_OFFSET_Y + SMALL_BBOX["ymin"] * SMALL_SCALE,
)

# 步骤 4:合并
pointdfsfall = pd.concat([pointdfsf, gdf_small], ignore_index=True)
print(f"合并后散点总数: {len(pointdfsfall)} 个")

4.4 绘制散点地图

# 获取省名标注坐标(质心)
centroids = provmap.geometry.representative_point()

# 为每个省分配颜色
all_colors = (
plt.colormaps.get_cmap("tab20").colors +
plt.colormaps.get_cmap("tab20b").colors +
plt.colormaps.get_cmap("Set3").colors
)
unique_provs = list(dict.fromkeys(pointdfsfall["省"].values))
color_map = {prov: all_colors[i % len(all_colors)] for i, prov in enumerate(unique_provs)}

# 绘图
此处代码需下载讲义材料查看~

五、绘制填充地图

5.1 读取市级地图并统计

虽然绘制市级地图,但我们仍然使用省级线条(突出省份区域):

# 读取市级地图
citymap = gpd.read_file(
"citymapdata/minishp/chinacity2021mini/chinacity2021mini.shp",
encoding="utf-8"
)
citymap = citymap.dropna(subset=["省代码"]).reset_index(drop=True)
print(f"市级多边形: {len(citymap)} 个区域")
# 省级线条(包含省份边界)
provlinemap2 = gpd.read_file(
"provmapdata/minishp/chinaprov2021mini/chinaprov2021mini_line.shp",
encoding="utf-8"
)
keep_classes2 = ["九段线", "海岸线", "小地图框格", "省份"]
provlinemap2 = provlinemap2[
provlinemap2["class"].isin(keep_classes2)
][["class", "geometry"]].reset_index(drop=True)
# 统计每个城市的公司数量
citydf = pointdf.groupby(["市", "市代码"]).size().reset_index(name="n")
print(citydf.head(10))
# 合并到地图数据
citymap2 = citymap.merge(citydf, on="市代码", how="left")
citymap2["n"] = citymap2["n"].fillna(0)
citymap2["plot_val"] = citymap2["n"] + 1 # log(0) 无定义
print(f"有数据的城市: {(citymap2['n'] > 0).sum()} 个")

5.2 连续变量填色地图

fig, ax = plt.subplots(1, 1, figsize=(10, 8.5), dpi=400)

cmap = plt.colormaps.get_cmap("YlOrBr")
norm = mcolors.LogNorm(vmin=citymap2["plot_val"].min(),
vmax=citymap2["plot_val"].max())

# 填充地图
citymap2.plot(ax=ax, column="plot_val", cmap=cmap, norm=norm,
edgecolor="gray", linewidth=0.01, zorder=1)

# 设置等比例,避免地图变形
ax.set_aspect("equal")
_xl = ax.get_xlim(); _yl = ax.get_ylim()
_dw = _xl[1] - _xl[0]; _dh = _yl[1] - _yl[0]
fig.set_size_inches(10, 10 * _dh / _dw)
ax.set_xlim(_xl); ax.set_ylim(_yl)

# 线条
draw_lines(ax, provlinemap2)

# 散点(叠加)
for prov, color in color_map.items():
mask = pointdfsfall["省"] == prov
pointdfsfall[mask].plot(ax=ax, color=color, markersize=0.5, zorder=2)

# 省名
draw_province_labels(ax, provmap, centroids)

# 装饰
此处代码需下载讲义材料查看~

5.3 连续变量填色地图(竖直图例)

也可以使用竖直图例展示连续变量:

fig, ax = plt.subplots(1, 1, figsize=(10, 8.5), dpi=400)

# 使用与 pic2 相同的 YlOrBr + LogNorm 配色
cmap3 = plt.colormaps.get_cmap("YlOrBr")
norm3 = mcolors.LogNorm(vmin=citymap2["plot_val"].min(),
vmax=citymap2["plot_val"].max())
citymap2.plot(ax=ax, column="plot_val", cmap=cmap3, norm=norm3,
edgecolor="gray", linewidth=0.01, legend=False, zorder=1)

# 设置等比例
ax.set_aspect("equal")
_xl3 = ax.get_xlim(); _yl3 = ax.get_ylim()
_dw3 = _xl3[1] - _xl3[0]; _dh3 = _yl3[1] - _yl3[0]
fig.set_size_inches(10, 10 * _dh3 / _dw3)
ax.set_xlim(_xl3); ax.set_ylim(_yl3)

# 线条
draw_lines(ax, provlinemap2)

# 散点
for prov, color in color_map.items():
mask = pointdfsfall["省"] == prov
pointdfsfall[mask].plot(ax=ax, color=color, markersize=0.5, zorder=2)

# 省名
draw_province_labels(ax, provmap, centroids)

# 装饰
此处代码需下载讲义材料查看~

六、使用栅格数据绘制地图

6.1 高分辨率栅格 → 散点

高分辨率的栅格数据可以先聚合(降采样),再转换成密集的散点数据绘制。相当于 R 中 raster::aggregate() + rasterToPoints() 的组合:

import rasterio
from rasterio.warp import Resampling
from shapely.geometry import Point as ShpPoint

# 读取并聚合栅格(每隔 10 个像素取均值)
tif_path = "cn2022.tif"
aggregate_factor = 10

with rasterio.open(tif_path) as src:
new_width = src.width // aggregate_factor
new_height = src.height // aggregate_factor
data = src.read(
out_shape=(src.count, new_height, new_width),
resampling=Resampling.average,
)[0]
from affine import Affine
transform = Affine(
src.transform.a * aggregate_factor,
src.transform.b,
src.transform.c,
src.transform.d,
src.transform.e * aggregate_factor,
src.transform.f,
)
src_crs = src.crs
nodata = src.nodata

print(f"原始尺寸: {src.width} x {src.height}")
print(f"聚合后: {new_width} x {new_height}")

# 转为 GeoDataFrame
data_float = data.astype(np.float64)
if nodata is not None:
data_float[data == nodata] = np.nan

rows, cols = np.where(~np.isnan(data_float) & (data_float > 0))
xs = [transform * (c + 0.5, r + 0.5) for r, c in zip(rows, cols)]
values = data_float[rows, cols]

cn2022points = gpd.GeoDataFrame(
{"cn2022": values},
geometry=[ShpPoint(x[0], x[1]) for x in xs],
crs=src_crs,
)
cn2022points = cn2022points[cn2022points["cn2022"] > 0].reset_index(drop=True)
cn2022points = cn2022points.to_crs(MY_CRS)

print(f"有效散点: {len(cn2022points)} 个")
print(cn2022points.describe())
# 拆分小地图部分
def split_for_small_map(gdf):
small_bbox = box(
SMALL_BBOX["xmin"], SMALL_BBOX["ymin"],
SMALL_BBOX["xmax"], SMALL_BBOX["ymax"]
)
gdf_small = gdf[gdf.geometry.intersects(small_bbox)].copy()
if len(gdf_small) > 0:
gdf_small["geometry"] = gdf_small.geometry.translate(
xoff=-SMALL_BBOX["xmin"], yoff=-SMALL_BBOX["ymin"])
gdf_small["geometry"] = gdf_small.geometry.scale(
xfact=SMALL_SCALE, yfact=SMALL_SCALE, origin=(0, 0))
gdf_small["geometry"] = gdf_small.geometry.translate(
xoff=SMALL_OFFSET_X + SMALL_BBOX["xmin"] * SMALL_SCALE,
yoff=SMALL_OFFSET_Y + SMALL_BBOX["ymin"] * SMALL_SCALE)
return pd.concat([gdf, gdf_small], ignore_index=True)
return gdf

cn2022points_all = split_for_small_map(cn2022points)
print(f"合并后: {len(cn2022points_all)} 个散点")
# 绘图
fig, ax = plt.subplots(1, 1, figsize=(10, 8.5), dpi=400)

cmap = plt.colormaps.get_cmap("RdYlBu_r")
norm = mcolors.LogNorm(vmin=cn2022points_all["cn2022"].min(),
vmax=cn2022points_all["cn2022"].max())

# 散点
cn2022points_all.plot(ax=ax, column="cn2022", cmap=cmap, norm=norm,
markersize=0.1, zorder=1)

# 省级边界
provmap.plot(ax=ax, facecolor="none", edgecolor="gray",
linewidth=0.01, zorder=2)

# 设置等比例
ax.set_aspect("equal")
_xl4 = ax.get_xlim(); _yl4 = ax.get_ylim()
_dw4 = _xl4[1] - _xl4[0]; _dh4 = _yl4[1] - _yl4[0]
fig.set_size_inches(10, 10 * _dh4 / _dw4)
ax.set_xlim(_xl4); ax.set_ylim(_yl4)

# 线条
draw_lines(ax, provlinemap)

# 省名
for i, row in provmap.iterrows():
cx, cy = centroids.iloc[i].x, centroids.iloc[i].y
ax.text(cx, cy, row["省"], fontsize=5, color="gray",
ha="center", va="center", zorder=4)

ax.set_axis_off()
ax.set_title("2022 年中国夜间灯光地图", fontsize=14,
fontweight="bold", pad=30)
ax.text(0.5, 1.03, "绘图:微信公众号 RStata",
transform=ax.transAxes, fontsize=9, ha="center", va="bottom",
color="gray")
add_scale_bar(ax, provmap)
add_north_arrow(ax)

# 水平色条(左下角比例尺上方)
此处代码需下载讲义材料查看~

6.2 低分辨率栅格 → 多边形

低分辨率的栅格数据适合转换成 polygon 绘制:

from rasterio.features import shapes as rasterio_shapes
from shapely.geometry import shape as shapely_shape

tif_path2 = "NH3_em_anthro_2015_sector_ENE.tif"

with rasterio.open(tif_path2) as src:
data = src.read(1)
transform = src.transform
src_crs = src.crs
nodata = src.nodata

data_float = data.astype(np.float64)
if nodata is not None:
data_float[data == nodata] = np.nan

# 栅格转多边形
mask = ~np.isnan(data_float) & (data_float > 0)
results = [
{"properties": {"NH3": v}, "geometry": shapely_shape(s)}
for s, v in rasterio_shapes(data_float, transform=transform, mask=mask)
if v > 0
]

cn2022polygons = gpd.GeoDataFrame.from_features(results, crs=src_crs)
cn2022polygons = cn2022polygons.to_crs(MY_CRS)
print(f"多边形: {len(cn2022polygons)} 个")
print(cn2022polygons["NH3"].describe())
# 拆分小地图
cn2022polygons_all = split_for_small_map(cn2022polygons)

# 裁剪到中国范围
china_union = unary_union(provmap.geometry)
cn2022polygons_all = cn2022polygons_all[
cn2022polygons_all.geometry.intersects(china_union)
].copy()
cn2022polygons_all["geometry"] = cn2022polygons_all.geometry.intersection(china_union)

# 计算色标范围
values = cn2022polygons_all["NH3"].dropna().values
values = values[values > 0]
log_vals = np.log10(values)
rmin, rmax = log_vals.min(), log_vals.max()
vmin = 10 ** (rmin + 0.02 * (rmax - rmin))
vmax = 10 ** (rmax - 0.02 * (rmax - rmin))
vmedian = 10 ** np.mean([rmin, rmax])
print(f"色标范围: {vmin:.4f} ~ {vmax:.4f}, 中位数: {vmedian:.4f}")
fig, ax = plt.subplots(1, 1, figsize=(10, 8.5), dpi=400)

cmap = plt.colormaps.get_cmap("YlOrRd_r")
norm = mcolors.LogNorm(vmin=vmin, vmax=vmax)

# 填充地图
cn2022polygons_all.plot(ax=ax, column="NH3", cmap=cmap, norm=norm,
edgecolor="none", linewidth=0, zorder=1)

# 设置等比例
ax.set_aspect("equal")
_xl5 = ax.get_xlim(); _yl5 = ax.get_ylim()
_dw5 = _xl5[1] - _xl5[0]; _dh5 = _yl5[1] - _yl5[0]
fig.set_size_inches(10, 10 * _dh5 / _dw5)
ax.set_xlim(_xl5); ax.set_ylim(_yl5)

# 省级边界
provmap.plot(ax=ax, facecolor="none", edgecolor="gray",
linewidth=0.01, zorder=2)

# 线条
draw_lines(ax, provlinemap)

# 省名
for i, row in provmap.iterrows():
cx, cy = centroids.iloc[i].x, centroids.iloc[i].y
ax.text(cx, cy, row["省"], fontsize=5, color="gray",
ha="center", va="center", zorder=4)

ax.set_axis_off()
ax.set_title("各地区 NH₃ 排放分布(单位:kg/m²/yr)", fontsize=14,
fontweight="bold", pad=30)
ax.text(0.5, 1.03, "绘图:微信公众号 RStata",
transform=ax.transAxes, fontsize=9, ha="center", va="bottom",
color="gray")
add_scale_bar(ax, provmap)
add_north_arrow(ax)

# 水平色条(左下角比例尺上方,宽度为原来的 1/3)
此处代码需下载讲义材料查看~

七、长版地图

长版地图数据的使用更简单,这里仅演示散点图。长版不包含小地图框格,没有南海小地图:

# 读取长版省级地图
provmap_long = gpd.read_file(
"provmapdata/longshp/chinaprov2019long/chinaprov2019long.shp",
encoding="utf-8"
)
provmap_long = provmap_long.dropna(subset=["省代码"]).reset_index(drop=True)
print(f"长版省级多边形: {len(provmap_long)} 个")
# 长版线条(不含小地图框格)
provlinemap_long = gpd.read_file(
"provmapdata/longshp/chinaprov2019long/chinaprov2019long_line.shp",
encoding="utf-8"
)
keep_long = ["九段线", "海岸线"]
provlinemap_long = provlinemap_long[
provlinemap_long["class"].isin(keep_long)
][["class", "geometry"]].reset_index(drop=True)

# 查看范围
bounds = provmap_long.total_bounds
print(f"X 范围: {bounds[0]:.0f} ~ {bounds[2]:.0f}")
print(f"Y 范围: {bounds[1]:.0f} ~ {bounds[3]:.0f}")
# 绘图
fig, ax = plt.subplots(1, 1, figsize=(8, 9), dpi=400)

# 底图
provmap_long.plot(ax=ax, facecolor="white", edgecolor="gray",
linewidth=0.01, zorder=1)

# 设置等比例
ax.set_aspect("equal")
_xl6 = ax.get_xlim(); _yl6 = ax.get_ylim()
_dw6 = _xl6[1] - _xl6[0]; _dh6 = _yl6[1] - _yl6[0]
fig.set_size_inches(8, 8 * _dh6 / _dw6)
ax.set_xlim(_xl6); ax.set_ylim(_yl6)

# 散点
for prov, color in color_map.items():
mask = pointdfsf["省"] == prov
pointdfsf[mask].plot(ax=ax, color=color, markersize=2, zorder=2)

# 线条
draw_lines(ax, provlinemap_long)

# 省名
centroids_long = provmap_long.geometry.representative_point()
for i, row in provmap_long.iterrows():
cx, cy = centroids_long.iloc[i].x, centroids_long.iloc[i].y
ax.text(cx, cy, row["省"], fontsize=5, color="gray",
ha="center", va="center", zorder=4)

# 装饰
ax.set_axis_off()
ax.set_title("瞪羚、独角兽、创新型企业的地理分布",
fontsize=14, fontweight="bold", pad=30)
ax.text(0.5, 1.03, "绘图:微信公众号 RStata",
transform=ax.transAxes, fontsize=9, ha="center", va="bottom",
color="gray")
add_scale_bar(ax, provmap_long, height_scale=0.016, width_scale=0.3)
add_north_arrow(ax)
fig.text(0.5, 0.01, "数据来源:瞪羚云网站", ha="center",
fontsize=7, color="gray")

fig.savefig("pic6.png", bbox_inches="tight", dpi=400)
plt.close(fig)
print("图片已保存: pic6.png")

八、模块化绘图函数

上面的代码逐步演示了完整的绘图流程。为了方便复用,所有绘图逻辑已被封装到 draw_china_map.py 模块中,可以直接调用:

from draw_china_map import (
read_provmap, read_provline, read_citymap, read_point_data,
plot_mini_scatter_map, plot_mini_choropleth,
plot_mini_choropleth_discrete, plot_mini_raster_points,
plot_mini_raster_polygons, plot_long_scatter_map,
)

# 示例:一键绘制所有 6 张地图
provmap = read_provmap(2019, version="mini")
provline = read_provline(2019, version="mini")
pointdfsf = read_point_data("瞪羚、独角兽、创新型企业经纬度数据.xlsx")

plot_mini_scatter_map(provmap, pointdfsf, provline, save_path="pic1.png")

更多关于地图绘制的内容可以学习平台上的系列课程「使用 R 语言进行地理计算」。

点击这里跳转到 RStata 短书平台获取附件:使用 Python 绘制历年中国省市区县地图(小地图版本+长版)

评论