用 Python 把矢量等高线转成栅格网格

ytkz2026-01-052026-01-05

用 Python 把矢量等高线转成栅格网格——顺便聊聊乐高地形模型

空间数据主要有矢量和栅格两种形式。矢量擅长精确描述点、线、面，比如城市边界、道路或等高线；栅格则用规则的网格像素来表示，每个像素带一个值，常用于高程、温度之类的连续数据。实际分析时，经常需要把两者结合起来，比如把矢量等高线转成高程栅格。

今天拿杭州为例，演示怎么用 Python 把等高线（矢量）转成规则网格（栅格），最后再离散化成几个高度级别——这样就能想象用不同颜色的乐高砖来搭一个迷你杭州地形模型了。整个过程主要靠 GeoPandas 和 OSMnx，代码简单，几分钟就能跑通。

数据准备

杭州公开的等高线数据不太好找，我这里用一个常见的替代方案：从公开 DEM（数字高程模型）生成等高线，或者直接用开源高程数据。这里为了演示方便，我们先用 OSMnx 下载杭州市行政边界，然后假设有一份等高线数据（实际项目可以从自然资源部门或开源平台获取 DEM 后生成）。

先看边界：

import osmnx as ox
import matplotlib.pyplot as plt

# 下载杭州市行政边界
city = ox.geocode_to_gdf('杭州市, 中国')
city = city.to_crs(epsg=4326)

# 快速看一眼
city.plot(figsize=(10, 10), edgecolor='black')
plt.title('杭州市行政边界')
plt.axis('off')
plt.show()

结果就是熟悉的杭州轮廓，西湖、钱塘江一览无余。

生成规则网格

接下来，把这个不规则的行政区切成规则方格网格。这里写个小函数，根据边界框分成指定行列数的格子：

from shapely.geometry import box
import geopandas as gpd

def create_grid(bounds, rows, cols):
    minx, miny, maxx, maxy = bounds
    width = (maxx - minx) / cols
    height = (maxy - miny) / rows
    grid = []
    for i in range(cols):
        for j in range(rows):
            x1 = minx + i * width
            y1 = miny + j * height
            x2 = x1 + width
            y2 = y1 + height
            grid.append(box(x1, y1, x2, y2))
    return gpd.GeoDataFrame(grid, columns=['geometry'], crs='EPSG:4326')

bounds = city.total_bounds  # [minx, miny, maxx, maxy]
gdf_grid = create_grid(bounds, 30, 30)  # 30x30 的网格，可自行调整

# 叠加显示
ax = city.plot(figsize=(12, 10), color='none', edgecolor='black', linewidth=2)
gdf_grid.plot(ax=ax, color='none', edgecolor='gray', alpha=0.6)
plt.title('杭州边界 + 30x30 网格')
plt.axis('off')
plt.show()

网格均匀铺在杭州边界上，后续我们会把每个格子填充平均高程值。

把等高线值映射到网格

假设我们有一份等高线 GeoDataFrame（列名里有 ‘elevation’ 表示高度值），叫 gdf_contours。

实际操作中可以用空间连接（sjoin）找出每个网格包含哪些等高线，然后取平均值：

# 空间连接：网格内包含的等高线部分
joined = gpd.sjoin(gdf_contours, gdf_grid, how='inner', predicate='intersects')

# 按网格索引聚合，取平均高程
aggregated = joined.groupby('index_right')['elevation'].mean().reset_index()

# 合并回网格
gdf_grid = gdf_grid.reset_index().rename(columns={'index': 'grid_id'})
gdf_grid = gdf_grid.merge(aggregated, left_on='grid_id', right_on='index_right', how='left')

gdf_grid.head()

现在每个格子都有一个平均高程值了。如果某些格子没交到等高线，会是 NaN，可以填 0 或附近值。

离散化高度级别（为乐高做准备）

乐高砖高度有限制，我们把连续高程转成几个离散级别，比如 0~6 级，方便用不同层数的砖表示：

import numpy as np

# 先填补空值
gdf_grid['elevation'] = gdf_grid['elevation'].fillna(0)

# 简单线性离散化到 0-8 级
min_val = gdf_grid['elevation'].min()
max_val = gdf_grid['elevation'].max()
gdf_grid['level'] = np.floor(8 * (gdf_grid['elevation'] - min_val) / (max_val - min_val + 1e-6))
gdf_grid['level'] = gdf_grid['level'].clip(0, 7)  # 限制在 0-7

print(gdf_grid['level'].min(), gdf_grid['level'].max())

可视化结果

最后画出来看看效果，用 terrain 色带更直观：

ax = gdf_grid.plot(column='level', figsize=(12, 10), cmap='terrain', 
                   edgecolor='gray', linewidth=0.3, legend=True)

city.plot(ax=ax, color='none', edgecolor='black', linewidth=1.5)

plt.title('杭州地形高度离散网格（乐高风格预览）')
plt.axis('off')
plt.show()