GeoPandas入门：从环境搭建到空间分析全流程实战

1. 为什么GeoPandas是地理空间分析绕不开的“第一块砖”

你手头有一份城市地铁站列表，带经纬度；一份行政区划边界图；一份过去三年的空气质量监测点数据；甚至是一张用手机随手拍的带GPS信息的照片——这些都不是普通表格，它们天然带着“位置”这个强属性。每天产生的数据里，有超过60%都隐含着地理坐标，但绝大多数人拿到后第一反应还是打开Excel，把经度纬度当两列普通数字处理。这就像拿着一把瑞士军刀，却只用它来拧螺丝——不是不行，但浪费了90%的功能。

GeoPandas就是那把被真正“解锁”的瑞士军刀。它不是从零造轮子，而是站在pandas这个数据科学界事实标准的肩膀上，给DataFrame装上了“空间感知能力”。你不需要重学一套语法，写 df.head() 和 gdf.head() 的体验几乎一模一样；但当你敲下 gdf.buffer(500) ，它就能瞬间给你画出每个地铁站500米半径的服务圈；用 gdf.sjoin(other_gdf, predicate='within') ，三秒内就能告诉你哪些监测点落在了某个生态保护区内部。这种平滑的学习曲线，正是它在科研、城市规划、物流调度、环境评估等领域迅速成为首选工具的核心原因。

我带过不少刚转行的数据分析师，他们最常问的问题不是“GeoPandas能做什么”，而是“我到底该不该现在就学？”。我的回答永远是：如果你的业务里哪怕出现一次“在哪”“多远”“是否重叠”“怎么画在地图上”这类问题，那就别犹豫。因为一旦你用pandas处理完数据清洗、特征工程，再用GeoPandas做空间关联，最后用matplotlib或contextily出图，整条链路完全在Python生态内闭环，没有格式转换、没有平台切换、没有数据丢失。我去年帮一家共享单车公司做热力图分析，从原始GPS轨迹点到最终呈现给运营团队的“高需求缺口区域”地图，整个流程脚本不到80行，而之前他们用GIS软件手动导出、配准、叠加、截图，一个图要花两小时。这不是炫技，是生产力的真实跃迁。

关键词“GeoPandas”“地理空间分析”“空间数据科学”背后，本质是一套思维范式的切换：从“记录是什么”转向“记录在哪里，以及它和周围的关系是什么”。这篇教程不讲虚的，接下来每一行代码、每一个参数、每一次报错，都是我在真实项目里踩过、修过、验证过的路径。你不需要是GIS专家，也不需要懂投影坐标系的数学推导，只要你会用pandas读CSV、算均值、做分组，你就能立刻上手解决实际问题。

2. 环境搭建：为什么conda是GeoPandas安装的“唯一推荐路径”

2.1 地理空间库的依赖地狱，不是传说

GeoPandas表面看只是一个Python包，但它背后是一个精密咬合的“地理空间技术栈”。它本身不直接处理坐标变换、不解析Shapefile二进制结构、不计算多边形面积——这些硬核活儿，全由四个底层C/C++库默默承担：

• Shapely ：负责所有几何对象（点、线、面）的创建、计算与关系判断。比如 point.within(polygon) 返回True/False，背后是Shapely调用GEOS库的复杂算法。
• Fiona ：地理空间数据的“读写引擎”。它让 geopandas.read_file() 能同时打开.shp、.geojson、.gpkg甚至.kml文件，靠的就是Fiona对GDAL/OGR库的封装。
• PyProj ：坐标系转换的“翻译官”。当你执行 gdf.to_crs(epsg=32633) ，PyProj在后台调用PROJ库，把WGS84经纬度精确换算成UTM坐标系下的米制单位。
• Rtree ：空间索引的“快递分拣中心”。没有它，判断10万个点是否落在某个省界内，得做10万次穷举计算；有了Rtree，它先快速筛出“可能相关”的几百个候选，再精准计算，速度提升百倍。

这四个库，每一个都依赖特定版本的系统级C库（如GEOS、GDAL、PROJ），而这些C库在Windows/macOS/Linux上的编译环境、动态链接库路径、甚至编译器版本都千差万别。用pip硬装，相当于让你自己去组装一台发动机——理论上可行，但现实中90%的人会在第一步 pip install shapely 时卡死在“Microsoft Visual C++ 14.0 is required”或者“gdal-config not found”。

提示：我见过最典型的失败案例，是一位同事在Windows上用pip装完所有包，代码跑起来报错 OSError: Could not find libspatialite.dll 。查了三天才发现，他装的spatialite版本和GDAL版本不兼容，而spatialite根本不是GeoPandas的直接依赖，是Fiona的间接依赖……这种嵌套式依赖冲突，手动解决成本极高。

2.2 conda：预编译二进制的“开箱即用”方案

conda的优势，在于它不是一个单纯的Python包管理器，而是一个 跨语言、跨平台的环境与依赖管理器 。Anaconda/Miniconda发行版早已为GeoPandas生态打包好了所有组合：特定版本的Python + 特定版本的GDAL + 特定版本的GEOS + 特定版本的PROJ，全部经过严格测试，确保二进制兼容。它提供的不是源码，而是编译好的 .so （Linux）、 .dylib （macOS）或 .dll （Windows）文件，直接扔进你的环境就能跑。

安装步骤极其简单，且 必须按顺序执行 ：

# 第一步：安装Miniconda（轻量版，比Anaconda快得多）
# Windows用户：下载 https://repo.anaconda.com/miniconda/Miniconda3-latest-Windows-x86_64.exe
# macOS用户：下载 https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-arm64.pkg (Apple Silicon) 或 Miniconda3-latest-MacOSX-x86_64.pkg (Intel)
# Linux用户：wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh && bash Miniconda3-latest-Linux-x86_64.sh

# 第二步：创建一个干净、专用的地理空间环境（强烈建议！）
conda create -n geo_env python=3.10
conda activate geo_env

# 第三步：用conda-forge频道安装GeoPandas（官方推荐，更新最及时）
conda install -c conda-forge geopandas

# 第四步：验证安装（关键！）
python -c "import geopandas as gpd; print(gpd.__version__)"

实测下来，这套流程在Windows 10/11、macOS Monterey/Ventura/Sonoma、Ubuntu 20.04/22.04上，成功率接近100%。我自己的工作流里，永远会为每个新项目新建一个conda环境，比如 conda create -n traffic_analysis python=3.11 ，这样不同项目的GDAL版本互不干扰，出了问题也能快速回滚。

注意：不要用 pip install geopandas 在已有的conda环境中覆盖安装！这会导致conda精心维护的二进制依赖被pip的源码编译版本破坏，后续大概率出现 ImportError: DLL load failed 。如果真要用pip（比如公司内网限制conda），请务必使用 pip install --only-binary=all geopandas 强制只装二进制轮子，并提前 conda install -c conda-forge gdal fiona pyproj rtree 打好底层基础。

2.3 验证环境：三行代码，排除90%的后续问题

安装完成后，别急着写分析代码，先运行这三行“体检”命令：

import geopandas as gpd
import matplotlib.pyplot as plt

# 1. 检查核心功能是否正常
print("GeoPandas版本:", gpd.__version__)
print("Shapely是否可用:", gpd.datasets.get_path('naturalearth_lowres') is not None)

# 2. 加载一个内置小数据集，测试读取和基础绘图
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
print("世界国家数据加载成功，共", len(world), "个国家")

# 3. 尝试最简单的绘图（不显示，只生成对象）
ax = world.plot()
print("绘图对象创建成功，类型:", type(ax))

如果这三行都顺利输出，恭喜，你的GeoPandas环境已经稳如磐石。如果卡在第二步，说明Fiona/GDAL读取失败，大概率是conda环境没激活或路径问题；如果卡在第三步绘图，可能是matplotlib后端配置问题（常见于服务器无GUI环境），此时加一句 plt.switch_backend('Agg') 即可。这些看似琐碎的验证，能帮你把问题消灭在萌芽，避免在后续复杂分析中陷入“不知道是代码错了还是环境坏了”的绝望循环。

3. GeoDataFrame：理解它，就理解了GeoPandas的全部灵魂

3.1 它不是DataFrame的“加强版”，而是“空间化”的继承体

很多初学者看到 GeoDataFrame ，下意识觉得它是pandas.DataFrame的“升级版”，仿佛装了更多功能的豪华版。这是个危险的误解。准确地说， GeoDataFrame 是 pandas.DataFrame 的一个 子类（subclass） ，它完全继承了DataFrame的所有方法（ head() , groupby() , merge() , query() ），但额外增加了一个核心属性： geometry列 。

你可以把它想象成一辆普通轿车（DataFrame）和一辆加装了GPS导航与激光雷达的智能车（GeoDataFrame）。普通车能载人、能加速、能转弯（对应DataFrame的各种数据操作）；智能车除了这些，还能实时知道“我在哪”（geometry列）、“前方500米有路口”（空间关系计算）、“我的行驶轨迹是一条线”（LineString几何类型）。关键在于，智能车的“导航系统”（geometry列）必须被明确指定，否则它就退化成普通车。

import geopandas as gpd
import pandas as pd

# 创建一个纯属性的DataFrame（没有空间信息）
df = pd.DataFrame({
    'name': ['A', 'B', 'C'],
    'population': [1000, 2000, 1500]
})

# 创建一个带坐标的GeoDataFrame（必须指定geometry列）
gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy([2.17, 2.18, 2.19], [41.40, 41.41, 41.42]), # 经度、纬度
    crs="EPSG:4326" # 必须声明坐标系！
)

print("普通DataFrame类型:", type(df))
print("GeoDataFrame类型:", type(gdf))
print("GeoDataFrame的geometry列类型:", type(gdf.geometry))
print("GeoDataFrame的crs信息:", gdf.crs)

这段代码清晰展示了核心差异： gdf 比 df 多了一个 geometry 属性，其类型是 GeoSeries ，并且 gdf.crs 存储了坐标系信息。没有 crs ，GeoPandas就不知道这些数字代表的是经纬度还是平面直角坐标，所有空间计算都会出错。

3.2 Geometry列：空间数据的“心脏”，必须且只能有一个

一个GeoDataFrame可以有多个列包含几何对象（比如 geometry 列存行政区边界， centroid 列存中心点， buffer 列存服务范围），但 有且仅有一个列被标记为“active geometry” ，也就是当前默认用于所有空间操作的列。这个列的名字默认是 geometry ，但你可以通过 gdf.set_geometry('centroid', inplace=True) 随时切换。

为什么必须唯一？因为所有空间操作都需要一个明确的“主体”。比如 gdf.buffer(1000) ，它是在问：“以哪个几何对象为中心，向外缓冲1000米？” 如果有多个几何列，GeoPandas无法猜测你的意图。这就像你不能对一个人同时执行“测量身高”和“称体重”两个动作而不指定目标——虽然都是身体指标，但操作对象必须明确。

# 基于Barcelona数据集，演示geometry列的设置与切换
districts = gpd.read_file("districts.geojson")

# 查看当前active geometry（通常是'geometry'）
print("当前active geometry列名:", districts._geometry_column_name)
print("当前geometry列前3个值:", districts.geometry.head(3))

# 计算每个区的中心点，并存入新列'centroid'
districts['centroid'] = districts.centroid

# 将'centroid'设为新的active geometry
districts = districts.set_geometry('centroid')

# 现在buffer操作将作用于中心点，而非原始多边形边界
districts_buffered = districts.buffer(500) # 得到500米半径的圆
print("缓冲后几何类型:", type(districts_buffered)) # <class 'geopandas.geoseries.GeoSeries'>

# 切换回原始geometry列
districts = districts.set_geometry('geometry')

这个切换机制非常实用。例如，在分析城市设施可达性时，你可能先用 districts.centroid 计算到最近地铁站的距离，再用 districts.geometry 计算该区被绿地覆盖的面积比例。灵活切换active geometry，让你无需反复复制数据，就能针对不同空间实体进行计算。

3.3 CRS：坐标系不是可选项，而是空间计算的“空气”

CRS（Coordinate Reference System，坐标参考系统）是GeoPandas里最容易被忽视、也最致命的概念。你可以把它理解为空间计算的“空气”——平时感觉不到，但一旦缺失，所有操作都会窒息。

• 地理坐标系（Geographic CRS） ：如 EPSG:4326 （WGS84），用经纬度表示位置。它的单位是“度”，但1度经度在赤道约111公里，在北极点约0公里。这意味着，直接用 df.geometry.distance(other_point) 计算两个经纬度点的距离，结果是“度”，毫无地理意义。
• 投影坐标系（Projected CRS） ：如 EPSG:32633 （UTM Zone 33N），把球面“压平”到二维平面，单位是“米”。在这里，1米就是1米，距离、面积、长度计算才真正准确。

我踩过最大的坑，就是在分析北京通勤数据时，直接用 EPSG:4326 计算了所有地铁站到居民区的距离，得到的结果是0.005、0.012这样的数字，我以为单位是公里，结果汇报时被领导当场指出：“0.005度？那是550米？你确定？” —— 这就是没转投影的典型后果。

正确的做法永远是：

1. 读取时明确声明CRS ： gpd.read_file(..., crs="EPSG:4326")
2. 分析前统一转投影 ： gdf = gdf.to_crs(epsg=32633)
3. 计算后再转回地理坐标系出图 （如果需要）： gdf_plot = gdf.to_crs(epsg=4326)

# 正确的空间距离计算流程
from shapely.geometry import Point

# 1. 创建一个点（Sagrada Familia），明确指定其CRS
sagrada_point = gpd.GeoSeries([Point(2.174368, 41.403657)], crs="EPSG:4326")

# 2. 将districts和sagrada_point都转到同一投影坐标系（西班牙适用）
districts_projected = districts.to_crs(epsg=2062) # EPSG:2062 是西班牙常用投影
sagrada_projected = sagrada_point.to_crs(epsg=2062)

# 3. 计算距离（单位：米）
districts_projected['dist_to_sagrada_m'] = districts_projected.centroid.distance(sagrada_projected.iloc[0])

# 4. 转回地理坐标系用于绘图（保持地图形状正确）
districts_for_plot = districts_projected.to_crs(epsg=4326)

记住这个铁律： 任何涉及距离、面积、长度、缓冲区（buffer）的操作，都必须在投影坐标系下进行。 否则，你的分析结果再漂亮，也是空中楼阁。

4. 核心空间操作：从读取到可视化，一条完整链路实录

4.1 数据读取： read_file() 为何能“自动识别”上百种格式？

geopandas.read_file() 的魔力，源于它对Fiona库的深度封装。Fiona本身支持超过100种矢量数据格式，包括：

• Shapefile (.shp) ：行业标准，但实际是三个文件（.shp, .shx, .dbf）的组合，必须放在一起。
• GeoJSON (.geojson) ：纯文本，单文件，Web友好，开发调试首选。
• GeoPackage (.gpkg) ：SQLite数据库，可存储多层、栅格、属性表，现代GIS推荐格式。
• KML (.kml) ：Google Earth格式，适合展示。
• FileGDB (.gdb) ：Esri私有格式，需额外驱动。

read_file() 的“自动识别”逻辑很简单：它先检查文件扩展名，再尝试用对应驱动打开；如果失败，就遍历所有支持的驱动，直到找到能成功读取的那个。所以，即使你传入一个没有扩展名的文件，它也可能猜对。

# 实战技巧：如何优雅地处理各种来源的数据？
import os

# 情况1：本地Shapefile（注意：传入.shp文件路径，read_file会自动找同目录的.shx和.dbf）
# districts_shp = gpd.read_file("data/barcelona_districts.shp")

# 情况2：GitHub上的GeoJSON（直接传URL，无需下载）
url_geojson = "https://raw.githubusercontent.com/jcanalesluna/bcn-geodata/master/districtes/districtes.geojson"
districts = gpd.read_file(url_geojson)

# 情况3：压缩包里的Shapefile（zip://协议）
# districts_zip = gpd.read_file("zip://data/districts.zip!districts.shp")

# 情况4：内存中的GeoJSON字符串（适用于API返回）
geojson_str = '{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[2.17,41.40]},"properties":{"name":"Sagrada"}}]}'
districts_mem = gpd.read_file(geojson_str)

print("数据源类型:", type(districts))
print("数据行数:", len(districts))
print("列名:", list(districts.columns))
print("几何类型分布:\n", districts.geometry.type.value_counts())

注意：读取大文件（如全国1:100万行政区划）时， read_file() 可能很慢。此时可以用 gpd.read_file(..., rows=1000) 先读前1000行快速预览，或用 fiona.open() 手动控制读取粒度。

4.2 空间属性计算： area , centroid , boundary , distance 的底层逻辑

这些看似简单的属性，背后是Shapely调用GEOS库的精密算法。理解它们的计算逻辑，能帮你避开无数陷阱。

• area ：对于 Polygon ，计算的是其在当前坐标系下的平面面积。如果CRS是 EPSG:4326 ，结果单位是“平方度”，毫无意义；必须先 to_crs() 到投影坐标系（如 EPSG:32633 ），结果才是平方米。 MultiPolygon 的面积是所有组成部分面积之和。

• centroid ：计算的是几何对象的“质心”（center of mass），不是“中心点”。对于一个细长的L型多边形，质心很可能落在图形外部！如果需要“视觉中心”，应该用 representative_point() ，它保证点一定在几何内部。

• boundary ：返回几何对象的“边界”。 Polygon 的边界是 LinearRing （闭合线）， LineString 的边界是其两个端点（ MultiPoint ）。这个操作常用于提取道路网络的节点。

• distance ：计算两个几何对象之间的 最小欧氏距离 。 Point.distance(Polygon) 返回点到多边形边界的最短距离； Polygon.distance(Polygon) 返回两个多边形边界的最短距离。注意，它不计算“路径距离”（如沿道路走的距离），那是网络分析的范畴。

# 深度实践：计算巴塞罗那各区的“紧凑度”（Compactness）
# 紧凑度 = 4 * π * area / perimeter²，值越接近1，形状越接近圆形
districts_proj = districts.to_crs(epsg=2062)

# 计算面积（平方米）和周长（米）
districts_proj['area_m2'] = districts_proj.area
districts_proj['length_m'] = districts_proj.length # 对Polygon，length是周长

# 计算紧凑度
districts_proj['compactness'] = (4 * 3.1415926 * districts_proj['area_m2']) / (districts_proj['length_m'] ** 2)

# 找出最“松散”和最“紧凑”的区
print("最松散的区（compactness最小）:", districts_proj.loc[districts_proj['compactness'].idxmin(), 'district_name'])
print("最紧凑的区（compactness最大）:", districts_proj.loc[districts_proj['compactness'].idxmax(), 'district_name'])

# 可视化紧凑度
ax = districts_proj.plot(column='compactness', legend=True, cmap='viridis', figsize=(10, 6))
plt.title("Barcelona Districts Compactness")
plt.show()

4.3 空间连接（Spatial Join）： sjoin() 的三种核心模式详解

gpd.sjoin() 是GeoPandas最强大的功能之一，它让“基于位置的关联”变得像SQL JOIN一样简单。但 predicate 参数的选择，直接决定了分析结果的业务含义。

# 场景：分析自行车道（bike_lane）与行政区（districts）的空间关系
# bike_lane 是 LineString，districts 是 Polygon

# 方式1：intersects（推荐用于线面关系）
lanes_in_districts = gpd.sjoin(bike_lane, districts, how="inner", predicate="intersects")
# 结果：每条自行车道对应它穿过的所有区，一条车道可能出现在多行

# 方式2：within（用于点面关系）
# 假设我们有自行车停放点（points）
# parking_points = gpd.read_file("parking.geojson")
# points_in_districts = gpd.sjoin(parking_points, districts, how="left", predicate="within")
# 结果：每个停车点对应一个区（或None），一行一记录

# 方式3：contains（用于面点关系）
# 找出哪些区“包含”了至少一个地铁站
# stations = gpd.read_file("stations.geojson")
# districts_with_stations = gpd.sjoin(districts, stations, how="inner", predicate="contains")
# 结果：每个有地铁站的区出现一次（即使有多个站）

# 实用技巧：统计每个区的自行车道总长度
lanes_in_districts['length_m'] = lanes_in_districts.length
districts_lane_length = lanes_in_districts.groupby('district_name')['length_m'].sum().reset_index()
districts_lane_length.columns = ['district_name', 'total_bike_lane_m']
print(districts_lane_length.sort_values('total_bike_lane_m', ascending=False).head())

注意： sjoin() 的 how 参数和pandas的 merge() 类似： inner 只保留匹配项， left 保留左表所有项（右表无匹配则为NaN）， right 反之。对于“每个区有多少条车道”这种聚合需求， inner + groupby 是最自然的组合。

4.4 空间叠加（Overlay）： overlay() 的四种操作与业务映射

gpd.overlay() 是空间分析的“瑞士军刀”，它能对两个多边形图层进行布尔运算，生成全新的空间关系。它的 how 参数定义了运算类型：

# 场景：模拟“在每个区中心建500米半径公园”的规划方案
# 1. 创建公园缓冲区（基于区中心点）
parks_buffer = districts_proj.centroid.buffer(500) # 单位：米
parks_gdf = gpd.GeoDataFrame(geometry=parks_buffer, crs=districts_proj.crs)

# 2. overlay操作：intersection（公园实际落地的区域）
parks_on_land = gpd.overlay(districts_proj, parks_gdf, how='intersection')
print("intersection结果行数:", len(parks_on_land)) # 应等于districts_proj行数

# 3. overlay操作：difference（剩余未被公园占用的土地）
land_without_park = gpd.overlay(districts_proj, parks_gdf, how='difference')
print("difference结果行数:", len(land_without_park)) # 可能大于districts_proj行数（因切割）

# 4. 可视化对比
fig, axes = plt.subplots(1, 2, figsize=(15, 6))

parks_on_land.plot(ax=axes[0], column='district_name', alpha=0.7, legend=True)
axes[0].set_title("Park Footprint (Intersection)")

land_without_park.plot(ax=axes[1], column='district_name', alpha=0.7, legend=True)
axes[1].set_title("Land Without Park (Difference)")

plt.tight_layout()
plt.show()

# 5. 计算每个区被公园占用的比例
# 先计算原始区面积
districts_proj['orig_area'] = districts_proj.area

# 再计算intersection后的面积（即被占用面积）
parks_on_land['park_area'] = parks_on_land.area

# 合并回原表
districts_with_park = districts_proj.merge(
    parks_on_land[['district_name', 'park_area']], 
    on='district_name', 
    how='left'
).fillna({'park_area': 0})

districts_with_park['park_ratio'] = districts_with_park['park_area'] / districts_with_park['orig_area']
print(districts_with_park[['district_name', 'park_ratio']].sort_values('park_ratio', ascending=False))

这个例子完美展示了 overlay() 如何将抽象的规划设想，转化为可量化、可比较、可可视化的空间结果。它不是魔法，而是严谨的几何运算，每一步都经得起推敲。

5. 可视化进阶：从静态地图到专业级空间表达

5.1 plot() 方法的隐藏参数：超越默认的10个关键选项

gdf.plot() 看似简单，但它的参数足以满足90%的静态出图需求。掌握这些参数，能让你的图从“能看”变成“专业”。

# 一张图，展示所有关键参数的实战效果
fig, ax = plt.subplots(1, 1, figsize=(12, 8))

# 核心参数详解
districts_proj.plot(
    ax=ax,
    
    # 1. column: 按哪一列的数值着色（必须是数值列）
    column='population_22', 
    
    # 2. cmap: 颜色映射（推荐viridis, plasma, coolwarm）
    cmap='viridis',
    
    # 3. legend: 是否显示图例
    legend=True,
    
    # 4. legend_kwds: 图例的详细设置
    legend_kwds={
        'label': "Population in 2022",
        'orientation': "horizontal", # 横向图例
        'shrink': 0.6, # 缩小到60%
        'aspect': 30 # 宽高比
    },
    
    # 5. scheme: 分类方式（'quantiles', 'equal_interval', 'fisher_jenks'）
    scheme='fisher_jenks', # 自动寻找最优断点
    
    # 6. k: 分类数量（配合scheme使用）
    k=5,
    
    # 7. alpha: 透明度（0-1），用于叠加图层
    alpha=0.8,
    
    # 8. edgecolor & linewidth: 边界颜色和粗细
    edgecolor='white',
    linewidth=0.5,
    
    # 9. categorical: 如果column是字符串，设为True进行分类着色
    # categorical=True,
    
    # 10. missing_kwds: 如何处理缺失值（NaN）
    missing_kwds={'color': 'lightgrey', 'label': 'No Data'}
)

# 添加标题和去除坐标轴
ax.set_title("Barcelona Districts Population (2022)", fontsize=16, fontweight='bold')
ax.set_axis_off()

plt.show()

实操心得： scheme='fisher_jenks' 是我最常用的分类方式，它通过算法将数据分成k组，使得组内差异最小、组间差异最大，比等间距（ equal_interval ）更能反映数据的真实分布。对于人口、收入这类偏态分布数据，效果尤其好。

5.2 contextily ：为你的地图注入真实的“底图灵魂”

gdf.plot() 生成的是纯矢量图，没有街道、没有建筑、没有卫星影像。 contextily 的作用，就是为这张矢量图“铺上”一张来自OpenStreetMap、Stamen或CartoDB的在线瓦片底图，让分析结果瞬间拥有真实世界的上下文。

import contextily as ctx

# 关键：底图的CRS必须与你的GeoDataFrame一致！
# 如果你的gdf是EPSG:4326，ctx.add_basemap会自动适配
# 如果你的gdf是投影坐标系（如EPSG:2062），必须显式指定
ax = districts_proj.plot(figsize=(12, 6), column='population_22', cmap='plasma', legend=True)

# 添加底图（自动适配CRS）
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)

# 或者，使用更简洁的写法（推荐）
# ctx.add_basemap(ax)

ax.set_title("Barcelona Population with OSM Basemap")
ax.set_axis_off()
plt.show()

ctx.providers 提供了丰富的底图选择：

• OpenStreetMap.Mapnik ：标准OSM地图，清晰易读。
• Stamen.Terrain ：地形图，突出海拔和地貌。
• CartoDB.Positron ：深色背景浅色要素，适合数据密集型图。
• ESRI.WorldImagery ：卫星影像，适合验证空间精度。

注意：首次使用某个底图源时， contextily 会自动下载并缓存瓦片，可能稍慢。后续使用会极快。如果内网环境无法访问外网，可以提前下载离线瓦片包，或使用 ctx.bounds2img() 生成静态图片。

5.3 多图层叠加：构建专业级空间分析图

真正的空间分析图，从来不是单一层。它往往是“底图 + 分析结果 + 关键点位 + 比例尺 + 指北针”的组合。 matplotlib 的 ax 对象是这一切的枢纽。

import matplotlib.patches as mpatches

# 构建一张综合分析图
fig, ax = plt.subplots(1, 1, figsize=(14, 10))

# 1. 绘制底图（OSM）
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, crs=districts_proj.crs.to_string())

# 2. 绘制行政区（半透明填充，白色边框）
districts_proj.plot(ax=ax, facecolor='none', edgecolor='white', linewidth=1.2, alpha=0.7)

# 3. 绘制自行车道（绿色线条）
bike_lane_proj = bike_lane.to_crs(epsg=2062)
bike_lane_proj.plot(ax=ax, color='green', linewidth=2, label='Bike Lanes')

# 4. 绘制Sagrada Familia（红色星标）
sagrada_proj = sagrada_point.to_crs(epsg=2062)
sagrada_proj.plot(ax=ax, color='red', marker='*', markersize=150, label='Sagrada Familia')

# 5. 绘制各区分界线的中心点（蓝色圆点）
districts_proj.centroid.plot(ax=ax, color='blue', markersize=30, alpha=0.6, label='District Centroids')

# 6. 添加图例（自定义图例元素）
legend_elements = [
    mpatches.Patch(facecolor='none', edgecolor='white', label='District Boundaries'),
    mpatches.Patch(color='green', label='Bike Lanes'),
    m