第3章：矢量与栅格数据模型

3.1 引言：为什么既需要“矢量”也需要“栅格”？

在实际的地理问题中，我们经常会遇到类似的困惑：用点线面描述河流与道路很直观，但对于植被覆盖、降雨强度、地形高程这些“连续变化”的信息，点线面却显得“颗粒太粗”。因此，GIS 形成了两种互补的数据模型：矢量（点/线/面）与栅格（规则像元网格）。本章通过问题引入与代码示例，讲清两者的分工、优势与局限，以及在同一工作流中如何相互转化与协同使用。

3.2 先修要求

• 基础 Python（with、函数、包管理）
• 基础 GIS 概念（坐标系、投影、空间参考）
• 简单的地理数据认识（例如 Shapefile、GeoJSON、TIFF）

3.3 核心概念与原理说明

• 矢量模型（Vector）：用点/线/面表达离散对象，适合道路网、行政区、河流要素等；属性结构清晰，拓扑可维护；空间精度高，但表达连续表面能力有限。
• 栅格模型（Raster）：规则网格像元表达连续现象（每个像元一个数值）。适合地形高程、遥感影像、降雨温度等；易做卷积、重采样、邻域分析；但边界“锯齿”且存储体量往往较大。
• 坐标参考与投影：两类数据必须共享一致的坐标参考与投影，才能正确叠加与分析。常见误区是“叠在一起但不重合”，多源于坐标系不一致或未正确重投影。
• 拓扑约束：矢量在叠置分析前应确保几何合法（无自交、无重复环等），否则会导致面积错误或分析异常。
• 分辨率与像元大小：栅格的像元大小决定空间细节与存储成本；过粗会丢失信息，过细则成本过高。分辨率选择要服务于问题尺度。
• 模型协同：常见工作流是“矢量区域筛选 → 栅格统计/卷积 → 输出区域统计表”。掌握两者互操作是工程落地的关键。

3.4 知识点与学习目标

• 认识矢量与栅格的适用场景与差异
• 掌握基本读写：geopandas 读写矢量，rasterio 读写栅格
• 学会基本空间分析：缓冲、叠置、栅格统计、重采样
• 掌握坐标系一致性与重投影的要点
• 能在一个任务中混合使用两种模型完成目标

3.5 内容讲解与示例（文章式展开）

选择合适的数据模型：从问题出发

当你的问题关注“对象的几何边界与关系”（如道路连通、行政区比较、河网结构），矢量更直观；当你的问题关注“空间连续分布的强度或状态”（如温度、降雨、植被），栅格更契合。一个常见错误是“强行用矢量表达连续现象”，导致分析复杂度陡增。

读写与坐标系一致性

在工程实践中，第一步往往是把数据正确读入与统一坐标参考：

# 代码片段：读取矢量与栅格并统一坐标系
import geopandas as gpd
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling# 读取矢量行政区
gdf = gpd.read_file("data/admin_boundary.geojson")
# 读取栅格
src = rasterio.open("data/ndvi.tif")# 将矢量重投影到栅格的坐标系（保持栅格原真值）
# 注意：矢量的 crs 属性可能为 None，需要手动指定或查阅数据文档
gdf = gdf.to_crs(src.crs)
print(gdf.crs)
print(src.crs)

经验要点：如果叠加后不重合，先检查两者 crs 是否一致，再检查是否存在坐标轴顺序问题（经纬度 vs 纬经度）。

矢量分析：缓冲、叠置与拓扑修复

矢量分析常见步骤是先修复几何，再做缓冲或叠置：

# 几何合法化，避免自交、多边形裂缝导致面积/叠置异常
gdf["geometry"] = gdf.buffer(0)# 对关键设施点位做缓冲区（例如 1000m 影响半径）
facilities = gpd.read_file("data/facilities.geojson").to_crs(src.crs)
buffer_1km = facilities.buffer(1000)# 叠置：求行政区内的设施缓冲覆盖区（intersection）
cover_in_admin = gpd.overlay(gpd.GeoDataFrame(geometry=buffer_1km, crs=src.crs),gdf,how="intersection"
)
# 输出用于后续统计
cover_in_admin.to_file("outputs/cover_in_admin.geojson", driver="GeoJSON")

工程提示：buffer(0) 并非万能，但能快速修复常见自交问题；复杂拓扑应使用专门工具（如 PostGIS 拓扑或 shapely 高级操作）。

栅格分析：裁剪、统计与重采样

当我们只关心行政区内的栅格信息时，应先裁剪再统计，以减少计算量并保证结果“在区域之内”：

import numpy as np
import rasterio.mask# 使用行政区边界对 NDVI 栅格裁剪
shapes = [geom for geom in gdf.geometry]
clipped_image, clipped_transform = rasterio.mask.mask(src, shapes, crop=True)
ndvi = clipped_image[0]# 统计：忽略 nodata，计算均值与百分位
nodata = src.nodata if src.nodata is not None else np.nan
ndvi = np.where(ndvi == nodata, np.nan, ndvi)
mean_ndvi = np.nanmean(ndvi)
p90_ndvi = np.nanpercentile(ndvi, 90)
print("mean=", mean_ndvi, "p90=", p90_ndvi)# 重采样：将 NDVI 重采样到更粗分辨率以加速下游分析
dst_crs = src.crs
transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds
)
dst_profile = src.profile.copy()
dst_profile.update({"transform": transform, "width": width, "height": height})dst_data = np.empty((src.count, height, width), dtype=src.dtypes[0])
for i in range(src.count):reproject(source=rasterio.band(src, i + 1),destination=dst_data[i],src_transform=src.transform,src_crs=src.crs,dst_transform=transform,dst_crs=dst_crs,resampling=Resampling.bilinear,)
# 保存重采样结果
with rasterio.open("outputs/ndvi_resampled.tif", "w", **dst_profile) as dst:dst.write(dst_data)

实践经验：统计前一定要正确处理 nodata，否则均值与百分位会被无效像元污染。重采样方法需根据数据类型选择：分类栅格常用最近邻，连续栅格常用双线性或立方卷积。

模型协同：区域-栅格汇总表

很多业务需要“在每个行政区计算一个栅格指标”，这就需要把栅格统计结果回写到矢量属性表：

import pandas as pdstats = []
for idx, row in gdf.iterrows():geo = [row.geometry]img, _ = rasterio.mask.mask(src, geo, crop=True)arr = img[0].astype(float)nodata = src.nodata if src.nodata is not None else np.nanarr = np.where(arr == nodata, np.nan, arr)stats.append({"admin_id": row.get("id", idx),"mean_ndvi": np.nanmean(arr),"p90_ndvi": np.nanpercentile(arr, 90)})df = pd.DataFrame(stats)
# 合并回矢量属性
gdf = gdf.merge(df, left_on=gdf.index, right_on=df.index, how="left")
gdf.to_file("outputs/admin_ndvi.geojson", driver="GeoJSON")

这样，一个既含有边界几何又含有统计指标的矢量数据就准备好了，下游可以直接用于制图或指标分类。

常见误区与排障思路（对应 3.8）

• 叠加不重合：检查 crs、是否经纬度顺序错误；必要时显式重投影。
• 面积/结果异常：排查矢量几何是否非法（自交、洞不闭合），使用 buffer(0) 尝试修复。
• 统计偏差：确保 nodata 被正确处理；分类栅格不要用双线性重采样。
• 性能问题：先裁剪再统计；按行政区分块处理；合理选择分辨率。

工程建议：数据组织与版本化

• 保持 data/（原始）与 outputs/（中间/结果）分离，避免覆盖原始数据。
• 在文件名中编码坐标系与分辨率信息（如 ndvi_wgs84_10m.tif）。
• 记录处理流水线（脚本名、参数、时间），便于复现与追踪。

3.6 实践任务（可运行代码）

• 任务一：将行政区与 NDVI 栅格对齐并输出区域均值与 P90 指标
• 任务二：为设施点生成 1km 缓冲并与行政区叠置，计算覆盖面积比例
• 任务三：将 NDVI 重采样到更粗分辨率并比较统计差异

参考脚本：gis_examples/modules/01_basics/vector_raster_demo.py

代码要点摘录：

# 读取矢量与栅格 → 坐标系对齐 → 栅格裁剪 → 统计 → 合并属性
# 脚本中包含必要的异常检查与注释，便于直接运行与修改

3.6.1 道路矢量转栅格热度与区域评分（端到端）

为方便在城市尺度表达道路通行状况，常将道路矢量转成“格网热度栅格”，随后按行政区统计并生成评分表：

# scripts/ch3/vector_to_raster_road_heat.py
import os
import numpy as np
import geopandas as gpd
import rasterio
from rasterio.transform import from_origin
from rasterio.features import rasterize
import rasterio.maskdef roads_to_heat(roads_path, cell_m=100, out_tif='outputs/ch3_road_heat.tif'):roads = gpd.read_file(roads_path).to_crs('EPSG:3857')minx, miny, maxx, maxy = roads.total_boundswidth = int((maxx - minx) // cell_m) + 1height = int((maxy - miny) // cell_m) + 1transform = from_origin(minx, maxy, cell_m, cell_m)# 线要素栅格化：占用叠加（简单近似，可替换为长度加权）shapes = [(geom, 1) for geom in roads.geometry if geom is not None]heat = rasterize(shapes, out_shape=(height, width), transform=transform, all_touched=True,merge_alg=rasterio.enums.MergeAlg.add).astype(np.float32)# 平滑近邻：作为“道路密度”直觉的简化实现from scipy.signal import convolve2dkernel = np.array([[1,1,1],[1,2,1],[1,1,1]], dtype=np.float32)heat_smooth = convolve2d(heat, kernel, mode='same', boundary='symm')profile = {'driver':'GTiff','dtype':'float32','count':1,'width':width,'height':height,'crs':'EPSG:3857','transform':transform}os.makedirs(os.path.dirname(out_tif), exist_ok=True)with rasterio.open(out_tif, 'w', **profile) as dst:dst.write(heat_smooth, 1)return out_tifdef score_accessibility(admin_path, heat_tif, out_csv='outputs/ch3_admin_accessibility.csv'):import pandas as pdgdf = gpd.read_file(admin_path).to_crs('EPSG:3857')with rasterio.open(heat_tif) as src:rows = []for idx, row in gdf.iterrows():img, _ = rasterio.mask.mask(src, [row.geometry], crop=True)arr = img[0]mean_heat = float(np.nanmean(arr))p90_heat = float(np.nanpercentile(arr, 90))score = 0.6*mean_heat + 0.4*p90_heatrows.append({'admin_id': row.get('id', idx), 'mean_heat': mean_heat, 'p90_heat': p90_heat, 'score': score})import pandas as pddf = pd.DataFrame(rows)os.makedirs(os.path.dirname(out_csv), exist_ok=True)df.to_csv(out_csv, index=False)return out_csvif __name__ == '__main__':roads = 'data/roads.geojson'admin = 'data/admin_boundary.geojson'tif = roads_to_heat(roads, cell_m=100)csv = score_accessibility(admin, tif)print('Saved', tif, csv)

运行：

• python scripts/ch3/vector_to_raster_road_heat.py
• 输出：outputs/ch3_road_heat.tif, outputs/ch3_admin_accessibility.csv

参数与调优建议：

• cell_m 设定网格大小；城市常见 50–200m，兼顾细节与计算量。
• 若需长度/等级加权，可对线要素按属性赋值后栅格化；或分段按长度权重。
• 平滑核仅作演示，工程中可替换为核密度估计/距离衰减以更符合交通直觉。

排障：

• 热度为零/道路未烧录：检查 all_touched、坐标系一致性与道路几何合法性。
• 行政区裁剪异常：检查几何并尝试 buffer(0) 修复；必要时简化几何提高鲁棒性。

3.7 输出与评估标准

• 统计结果完整：每个行政区都有 mean_ndvi 与 p90_ndvi
• 可复现：脚本能在新环境中按 README 步骤跑通
• 质量控制：坐标系一致、无明显几何错误、nodata 处理正确
• 表达清晰：输出文件命名规范、指标含义明确

3.8 常见错误与排障

• 错误：CRS mismatch → 解决：统一 to_crs 或对栅格重投影
• 错误：TopologyException → 解决：buffer(0) 或修复几何
• 错误：统计含 nodata → 解决：在统计前替换或屏蔽
• 错误：重采样方法不当 → 解决：按数据类型选择最近邻/双线性等

3.9 延伸阅读与资源

• 文档：geopandas、rasterio、shapely 官方指南
• 实践：PostGIS 栅格扩展、GDAL 工具链
• 案例：土地覆盖分类、区域生态指数计算

3.10 本章总结

本章从问题引入出发，系统阐释了矢量与栅格两种互补的数据模型及其协同方式。关键在于：让几何表达离散对象，让像元表达连续现象；通过坐标系一致、拓扑合法与合适的重采样策略，把两者在同一工作流中无缝拼接。实践部分给出裁剪-统计-合并的标准范式，便于在行政区尺度生成栅格摘要指标，为后续制图与决策提供坚实基础。