地理绘图 | Python 绘制土地利用分类图——CLCD 数据

geodata/
├── __init__.py              # 包入口
├── loader.py                # 核心加载器类
└── data/                    # 数据仓库
    ├── DEM数据/
    │   ├── dem_1km.tif     # 1公里分辨率
    │   └── dem_250m.tif    # 250米分辨率
    ├── 90mDEM/
    │   └── 90mdem.tif      # 90米分辨率
    ├── DEM30m/
    │   └── COPDEM.tif      # 30米高精度
    ├── 土地利用/
    │   └── CLCD_v01_2024_albert.tif  # CLCD 2024
    └── 矢量数据/
        ├── 省_Areas.shp
        ├── 市_Areas.shp
        └── 县_Areas.shp

from geodata import GeoDataLoader

# 创建加载器实例
loader = GeoDataLoader()

# 方法一：加载 CLCD 土地利用数据
clcd = loader.load_clcd(2024)  # 加载 2024 年数据

# 方法二：加载矢量边界数据
provinces = loader.load_vector("省", "Areas")  # 省级边界
cities = loader.load_vector("市", "Areas")     # 市级边界
counties = loader.load_vector("县", "Areas")   # 县级边界

# 方法三：加载 DEM 数据（可选）
dem_250m = loader.load_dem("250m")  # 250米分辨率
dem_30m = loader.load_dem("30m")    # 30米高精度

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.patches as mpatches
import cartopy.crs as ccrs
import geopandas as gpd
import rasterio
from rasterio.mask import mask
import numpy as np
from geodata import GeoDataLoader

# 创建数据加载器
loader = GeoDataLoader()

# 目标区域配置
TARGET_LEVEL = "市"# 行政级别：省/市/县
TARGET_NAME = "毕节市"# 区域名称
CLCD_YEAR = 2024# 数据年份

# 图框配置
EXPAND_RATIO = 0.05# 边界外扩 5%
LON_INTERVAL = 30/60# 经度间隔 30 分钟
LAT_INTERVAL = 30/60# 纬度间隔 30 分钟

# 字体配置
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['axes.unicode_minus'] = False

# geodata/loader.py 的核心实现
classGeoDataLoader:
def__init__(self):
        self.data_dir = Path(__file__).parent / "data"
        self.clcd_dir = self.data_dir / "土地利用"
        self.vector_dir = self.data_dir / "矢量数据"

defload_clcd(self, year=2024):
        clcd_path = self.clcd_dir / f"CLCD_v01_{year}_albert.tif"
return rasterio.open(clcd_path)

# CLCD 官方分类配色方案
CLCD_CLASSES = {
1: {'name': '农田', 'color': '#fae3c8'},      # 浅黄色
2: {'name': '森林', 'color': '#446f33'},      # 深绿色
3: {'name': '灌木', 'color': '#33a02c'},      # 中绿色
4: {'name': '草地', 'color': '#abd37b'},      # 浅绿色
5: {'name': '水域', 'color': '#1e69b4'},      # 蓝色
6: {'name': '冰雪', 'color': '#a6cee3'},      # 浅蓝色
7: {'name': '裸地', 'color': '#cfbda3'},      # 灰褐色
8: {'name': '不透水面', 'color': '#e24c42'},  # 红色（城市）
9: {'name': '湿地', 'color': '#2898e8'}       # 青蓝色
}

defcreate_clcd_cmap():
"""创建 CLCD 专用配色方案"""
    colors = [CLCD_CLASSES[i]['color'] for i in sorted(CLCD_CLASSES.keys())]
    cmap = mcolors.ListedColormap(colors)
# 边界值: 0.5, 1.5, 2.5, ..., 9.5
    bounds = [i - 0.5for i in range(1, len(CLCD_CLASSES) + 2)]
    norm = mcolors.BoundaryNorm(bounds, cmap.N)
return cmap, norm

配色的科学性

这套配色方案经过精心设计：

• 植被梯度：森林（深绿）→ 灌木（中绿）→ 草地（浅绿），颜色深浅反映植被密度
• 水体系统：水域（蓝色）→ 湿地（青蓝色）→ 冰雪（浅蓝色），都在蓝色系中
• 人工地表：农田（浅黄色）和不透水面（红色）使用暖色调，与自然植被的冷色调形成对比
• 裸地：灰褐色，中性色调，表示无植被覆盖的自然地表

ListedColormap 与 BoundaryNorm

ListedColormap 创建一个离散的配色方案，每个类别对应一个固定的颜色。BoundaryNorm 定义了类别之间的边界，确保值 1 映射到第一个颜色，值 2 映射到第二个颜色，以此类推。这种组合非常适合分类数据的可视化。

第三步：加载矢量边界与 CLCD 数据

接下来，需要加载两类数据：矢量边界（定义研究区域）和 CLCD 栅格数据（土地利用分类）。

# 加载矢量边界
gdf = loader.load_vector(TARGET_LEVEL, "Areas")
target_gdf = gdf[gdf['name'] == TARGET_NAME]

# 加载 CLCD 数据
src = loader.load_clcd(CLCD_YEAR)
print(f"CLCD 原始坐标系: {src.crs}")

# 将矢量转换到 CLCD 的坐标系进行裁剪
target_gdf_proj = target_gdf.to_crs(src.crs)
temp_nodata = int(src.nodata if src.nodata isnotNoneelse255)

# 裁剪 CLCD 数据到目标区域
out_image, out_transform = mask(
    src, 
    target_gdf_proj.geometry, 
    crop=True, 
    nodata=temp_nodata
)

坐标系转换的必要性

CLCD 数据使用 Albers 等积投影（一种适合中国的投影坐标系），而矢量边界通常使用 WGS84 地理坐标系。要用矢量边界裁剪 CLCD，必须先将矢量转换到 CLCD 的坐标系。

mask 函数的作用

mask() 函数是 rasterio 库的核心功能之一。它接收一个栅格数据集和一个矢量几何体，返回裁剪后的栅格数据。crop=True 参数让它裁剪到最小边界矩形，大大减少数据量，提高后续处理速度。

第四步：坐标系重投影

裁剪后的 CLCD 数据仍然在 Albers 投影坐标系中，但为了与底图和矢量边界正确对齐，需要将它重投影到 WGS84 地理坐标系。这是整个流程中最关键的技术步骤。

from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'EPSG:4326'# WGS84 地理坐标系

# 计算目标转换参数
transform_wgs84, width_wgs84, height_wgs84 = calculate_default_transform(
    src.crs,                    # 源坐标系
    dst_crs,                    # 目标坐标系
    out_image.shape[2],         # 源宽度
    out_image.shape[1],         # 源高度
    *rasterio.transform.array_bounds(out_image.shape[1], out_image.shape[2], out_transform)
)

# 创建目标数组
clcd_wgs84 = np.empty((height_wgs84, width_wgs84), dtype=out_image.dtype)

# 执行重投影
reproject(
    source=out_image[0],
    destination=clcd_wgs84,
    src_transform=out_transform,
    src_crs=src.crs,
    dst_transform=transform_wgs84,
    dst_crs=dst_crs,
    resampling=Resampling.nearest  # 最近邻插值，保持分类值不变
)

# 计算重投影后的地理范围
left, bottom, right, top = rasterio.transform.array_bounds(
    height_wgs84, width_wgs84, transform_wgs84
)
data_extent = [left, right, bottom, top]

print(f"重投影完成: {src.crs} -> {dst_crs}")

为什么要重投影？

如果不重投影，直接在 Albers 坐标系中绘制 CLCD 数据，会出现以下问题：

1. 与 WGS84 坐标系的底图无法对齐，导致 CLCD 数据显示在错误的位置
2. 与 WGS84 坐标系的矢量边界无法匹配，边界线和 CLCD 数据错位
3. 经纬度网格无法正确显示

Resampling.nearest 的重要性

对于分类数据，必须使用最近邻插值（nearest neighbor）。如果使用双线性或三次插值，会产生中间值（如 1.5、2.3），破坏分类的离散性。最近邻插值确保每个像素的值仍然是整数（1、2、3...），对应明确的土地类型。

calculate_default_transform 的作用

这个函数自动计算重投影所需的参数：新的仿射变换矩阵、新的宽度和高度。它确保重投影后的栅格覆盖相同的地理区域，且像素大小合理。

第五步：数据预处理与统计

重投影完成后，需要对数据进行预处理，并统计各土地类型的分布情况。

# 数据预处理
clcd_data = clcd_wgs84.astype('float32')
clcd_data[clcd_data == temp_nodata] = np.nan  # 无效值转为 NaN

# 统计各类别
unique_values, counts = np.unique(clcd_data[~np.isnan(clcd_data)], return_counts=True)
print(f"CLCD 数据形状: {clcd_data.shape}")
print(f"包含的土地类型: {[int(v) for v in unique_values]}")

# 计算地图范围（区域边界外扩 5%）
target_bounds = target_gdf.total_bounds
width = target_bounds[2] - target_bounds[0]
height = target_bounds[3] - target_bounds[1]

map_extent = [
    target_bounds[0] - width * EXPAND_RATIO,
    target_bounds[2] + width * EXPAND_RATIO,
    target_bounds[1] - height * EXPAND_RATIO,
    target_bounds[3] + height * EXPAND_RATIO
]

为什么要统计类别？

统计各土地类型有两个作用：

1. 验证数据：确认数据加载和处理正确，包含预期的土地类型
2. 优化图例：图例只显示数据中实际存在的类别，而不是全部 9 类。如果某个区域没有冰雪，图例就不显示冰雪，让图例更简洁

地图范围的外扩

外扩 5% 是一个经过实践检验的最佳值。它让地图不会显得过于拥挤，相邻区域的标注有足够的空间显示，同时又不会让目标区域显得太小。

第六步：绘制地图

将 CLCD 数据渲染成彩色的土地利用分类图。地图绘制是一个图层叠加的过程，每个图层都有其特定的作用和顺序。

# 创建图形和坐标轴
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent(map_extent, crs=ccrs.PlateCarree())

# 第一层：添加底图
ax.add_image(tiler, 8)

# 第二层：白色背景（确保 CLCD 区域可见）
ax.add_geometries(
    target_gdf.geometry,
    crs=ccrs.PlateCarree(),
    facecolor='white',
    edgecolor='none',
    zorder=4
)

# 第三层：绘制 CLCD 分类数据
cmap, norm = create_clcd_cmap()
clcd_masked = np.ma.masked_invalid(clcd_data)

im = ax.imshow(
    clcd_masked, 
    origin='upper', 
    extent=data_extent, 
    transform=ccrs.PlateCarree(),
    cmap=cmap,
    norm=norm,
    alpha=0.85,
    zorder=5,
    interpolation='nearest'
)

# 第四层：绘制目标区域边界
ax.add_geometries(
    target_gdf.geometry, 
    crs=ccrs.PlateCarree(), 
    facecolor='none', 
    edgecolor='black', 
    linewidth=1.5, 
    zorder=10
)

图层顺序的设计

这个图层顺序经过精心设计：

1. 底图（zorder=0）：提供地理背景，显示海洋、湖泊等水体
2. 白色背景（zorder=4）：在目标区域铺一层白色，确保 CLCD 数据在任何底图上都清晰可见
3. CLCD 数据（zorder=5）：土地利用分类的核心图层
4. 边界线（zorder=10）：清晰标识研究区域的范围

alpha=0.85 的微妙之处

将 CLCD 图层设置为 85% 不透明度，让底图的水体部分微微透出来。这样，CLCD 数据中的水域（蓝色）与底图的水体（也是蓝色）相互呼应，增强视觉效果。

interpolation='nearest' 的必要性

对于分类数据，必须使用最近邻插值。如果使用双线性或三次插值，matplotlib 会在不同类别之间创建平滑过渡，产生不存在的中间颜色，破坏分类的清晰性。

第七步：添加地图要素

一张专业的地图必须包含完整的地图要素：边界、标注、经纬度网格、指北针、比例尺和图例。这些元素让地图不仅美观，更具有科学价值。

# 添加区域标注
target_center = target_gdf.geometry.representative_point().iloc[0]
ax.text(
    target_center.x, target_center.y, TARGET_NAME, 
    transform=ccrs.PlateCarree(),
    fontfamily='SimSun', 
    fontsize=14, 
    fontweight='bold', 
    ha='center', va='center',
    zorder=15,
    path_effects=[pe.withStroke(linewidth=3, foreground='white', alpha=0.7)]
)

# 设置经纬度网格
fplt.set_map_ticks(ax, map_extent, dx=LON_INTERVAL, dy=LAT_INTERVAL)
ax.tick_params(axis='both', direction='out', length=5, width=1)

# 添加指北针和比例尺
fplt.add_compass(ax, 0.95, 0.88, size=20, style='arrow')

# 自动计算比例尺长度
map_width_km = (map_extent[1] - map_extent[0]) * 111
target_length = map_width_km / 4
standard_lengths = [5, 10, 20, 25, 50, 75, 100, 150, 200, 250, 500, 1000]
scale_length = min(standard_lengths, key=lambda x: abs(x - target_length))
fplt.add_scale_bar(ax, 0.05, 0.05, length=scale_length)

# 创建图例（只显示实际存在的类别）
legend_elements = []
for value in sorted(unique_values):
if int(value) in CLCD_CLASSES:
        class_info = CLCD_CLASSES[int(value)]
        legend_elements.append(
            mpatches.Patch(
                facecolor=class_info['color'], 
                edgecolor='black',
                linewidth=0.5,
                label=class_info['name']
            )
        )

legend = ax.legend(
    handles=legend_elements,
    loc='lower right',
    title=f'土地利用类型({CLCD_YEAR}年)',
    title_fontproperties={'family': 'SimSun', 'size': 9, 'weight': 'bold'},
    prop={'family': 'SimSun', 'size': 8},
    framealpha=0.95,
    edgecolor='black'
)
legend.set_zorder(20)  # 确保图例在所有元素之上

文字描边的重要性

path_effects=[pe.withStroke(linewidth=3, foreground='white', alpha=0.7)] 给文字添加白色描边。这个技巧让文字在任何背景上都清晰可读——无论背景是深色的森林还是浅色的农田。

比例尺的智能计算

比例尺长度根据地图宽度自动计算（约为地图宽度的 1/4），然后从标准长度列表中选择最接近的值。这确保比例尺长度总是"整数"（如 50km、100km），而不是奇怪的数字（如 73km）。

图例的优化

图例只显示数据中实际存在的土地类型。如果某个区域没有冰雪或湿地，图例就不显示这些类别，让图例更简洁、更有针对性。

第八步：导出与应用

最后一步是将地图导出为高分辨率图片，用于报告、论文或演示。

# 导出图片
filename = f"{TARGET_NAME}_CLCD_{CLCD_YEAR}.jpg"
plt.savefig(filename, dpi=1200, bbox_inches="tight", facecolor='white')
print(f"已保存: {filename}")
plt.close()

1200 DPI 的意义

DPI（Dots Per Inch）表示每英寸的像素数。1200 DPI 是高质量印刷的标准，远超屏幕显示所需的 72-150 DPI。这样的图片可以直接用于学术论文、书籍出版，甚至大幅海报。

bbox_inches='tight' 的作用

这个参数让 matplotlib 自动裁剪掉图片周围的空白区域，确保导出的图片紧凑而不浪费空间。

总结

通过这个内容，完成了从原始 CLCD 数据到专业土地利用分类图的完整流程。这个过程涉及多个关键技术：

数据处理技术：

• 矢量数据的加载和空间查询
• 栅格数据的裁剪和重投影
• 坐标系转换和对齐
• 无效值处理和数据统计

可视化技术：

• 离散配色方案的创建
• 多图层的叠加和顺序管理
• 地图要素的添加和布局
• 高分辨率图片的导出

设计思想：

• 数据与代码的解耦（geodata 包）
• 参数化配置，提高灵活性
• 自动化计算，减少手动调整
• 用户友好的输出和提示

基于PyGIS绘制研究区土地利用类型概况图，完全脱离GIS相关软件，相比传统的ArcGIS绘制研究区概况图更加省时省力，而且支持批量导出，批量渲染，更加方便快捷。

后续小编会更新不同尺度的研究区地形概况图。