znlgis 博客

GIS开发与技术分享

第20章：仿射变换与坐标操作

仿射变换是几何学中的基本操作，能够对空间要素进行平移、旋转、缩放和倾斜等变换，同时保持直线和平行性不变。本章将全面介绍 GeoPandas 中仿射变换的实现方法，以及底层坐标的提取、修改和精度控制等操作。

20.1 仿射变换概述

20.1.1 什么是仿射变换

仿射变换（Affine Transformation）是线性变换加平移的组合，包含以下基本类型：

  • 平移（Translation）：将几何体沿 X/Y 方向移动指定距离
  • 旋转（Rotation）：将几何体绕某一点旋转指定角度
  • 缩放（Scaling）：将几何体按比例放大或缩小
  • 倾斜（Skew）：将几何体沿某一轴进行剪切变换
性质 说明
保持直线性 直线变换后仍为直线
保持平行性 平行线变换后仍然平行
保持比例关系 同一直线上点的比值不变
不保持角度 角度可能变化（纯旋转除外）

20.1.2 仿射变换矩阵

二维仿射变换用 6 个参数 [a, b, d, e, xoff, yoff] 表示，变换公式为：

x' = a * x + b * y + xoff
y' = d * x + e * y + yoff
参数 含义 恒等变换值
a X 方向缩放因子 1
b Y 对 X 的倾斜因子 0
d X 对 Y 的倾斜因子 0
e Y 方向缩放因子 1
xoff X 方向平移量 0
yoff Y 方向平移量 0

注意： 恒等变换（不做任何改变）对应的矩阵参数为 [1, 0, 0, 1, 0, 0]。

20.2 GeoSeries.affine_transform() 方法

20.2.1 方法签名与参数

GeoSeries.affine_transform(matrix)
参数 类型 说明
matrix list/tuple 6 元素仿射矩阵 [a, b, d, e, xoff, yoff]

返回值为变换后的 GeoSeries。

20.2.2 基础用法示例

import geopandas as gpd
from shapely.geometry import Polygon

gdf = gpd.GeoDataFrame(
    {"name": ["矩形A", "矩形B"]},
    geometry=[
        Polygon([(0, 0), (2, 0), (2, 1), (0, 1)]),
        Polygon([(3, 0), (5, 0), (5, 1), (3, 1)])
    ]
)

# 应用恒等变换
result = gdf.geometry.affine_transform([1, 0, 0, 1, 0, 0])
print(result)

输出：

0    POLYGON ((0 0, 2 0, 2 1, 0 1, 0 0))
1    POLYGON ((3 0, 5 0, 5 1, 3 1, 3 0))
dtype: geometry

20.3 平移变换（translate）

20.3.1 使用 affine_transform 实现平移

平移矩阵为 [1, 0, 0, 1, xoff, yoff]：

import geopandas as gpd
from shapely.geometry import Polygon

gdf = gpd.GeoDataFrame(
    {"name": ["建筑物"]},
    geometry=[Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])]
)

# 向右平移 10，向上平移 5
gdf["translated"] = gdf.geometry.affine_transform([1, 0, 0, 1, 10, 5])
print(gdf["translated"].iloc[0])

输出：

POLYGON ((10 5, 14 5, 14 8, 10 8, 10 5))

20.3.2 使用 Shapely translate 函数

from shapely.affinity import translate

# 逐要素平移
gdf["translated"] = gdf.geometry.apply(
    lambda geom: translate(geom, xoff=10, yoff=5)
)

注意： affine_transform() 是向量化操作，性能优于逐要素 apply。数据量较大时应优先使用 affine_transform()。

20.4 旋转变换（rotate）

20.4.1 旋转矩阵推导

绕原点逆时针旋转角度 θ 对应参数 [cos(θ), -sin(θ), sin(θ), cos(θ), 0, 0]。

20.4.2 绕原点旋转

import math
import geopandas as gpd
from shapely.geometry import Polygon

square = Polygon([(1, 0), (2, 0), (2, 1), (1, 1)])
gdf = gpd.GeoDataFrame({"name": ["正方形"]}, geometry=[square])

# 绕原点逆时针旋转 90 度
theta = math.radians(90)
cos_t, sin_t = math.cos(theta), math.sin(theta)
gdf["rotated"] = gdf.geometry.affine_transform(
    [cos_t, -sin_t, sin_t, cos_t, 0, 0]
)

print("原始：", gdf.geometry.iloc[0])
print("旋转后：", gdf["rotated"].iloc[0])

输出：

原始： POLYGON ((1 0, 2 0, 2 1, 1 1, 1 0))
旋转后： POLYGON ((0 1, 0 2, -1 2, -1 1, 0 1))

20.4.3 绕任意点旋转

Shapely 的 rotate 函数可直接指定旋转中心：

from shapely.affinity import rotate
from shapely.geometry import Polygon
import geopandas as gpd

rect = Polygon([(0, 0), (4, 0), (4, 2), (0, 2)])
gdf = gpd.GeoDataFrame({"name": ["矩形"]}, geometry=[rect])

# 绕质心旋转 45 度
gdf["rotated_centroid"] = gdf.geometry.apply(
    lambda geom: rotate(geom, angle=45, origin='centroid')
)

# 绕指定点旋转
gdf["rotated_point"] = gdf.geometry.apply(
    lambda geom: rotate(geom, angle=45, origin=(2, 1))
)

print("绕质心旋转 45°：")
print(gdf["rotated_centroid"].iloc[0])

输出：

绕质心旋转 45°：
POLYGON ((0.5857864376269051 -0.4142135623730949, 3.414213562373095 -0.4142135623730951, 3.4142135623730954 2.414213562373095, 0.5857864376269053 2.4142135623730947, 0.5857864376269051 -0.4142135623730949))

注意： rotate() 的 origin 参数支持 'center'（边界框中心）、'centroid'（质心）或 (x, y) 坐标元组。角度为逆时针方向，单位为度。

20.5 缩放变换（scale）

20.5.1 均匀缩放与非均匀缩放

缩放矩阵为 [sx, 0, 0, sy, 0, 0]：

import geopandas as gpd
from shapely.geometry import Polygon

rect = Polygon([(0, 0), (2, 0), (2, 1), (0, 1)])
gdf = gpd.GeoDataFrame({"name": ["矩形"]}, geometry=[rect])

# 均匀缩放 2 倍
gdf["uniform"] = gdf.geometry.affine_transform([2, 0, 0, 2, 0, 0])

# 非均匀缩放：X 放大 3 倍，Y 放大 1.5 倍
gdf["nonuniform"] = gdf.geometry.affine_transform([3, 0, 0, 1.5, 0, 0])

print("原始面积：", gdf.geometry.iloc[0].area)
print("均匀 2x 面积：", gdf["uniform"].iloc[0].area)
print("非均匀面积：", gdf["nonuniform"].iloc[0].area)

输出：

原始面积： 2.0
均匀 2x 面积： 8.0
非均匀面积： 9.0

20.5.2 使用 Shapely scale 函数

from shapely.affinity import scale

# 以质心为中心放大 2 倍
gdf["scaled"] = gdf.geometry.apply(
    lambda geom: scale(geom, xfact=2, yfact=2, origin='centroid')
)
print(gdf["scaled"].iloc[0])

输出：

POLYGON ((-1 -0.5, 3 -0.5, 3 1.5, -1 1.5, -1 -0.5))

注意： affine_transform 的缩放中心始终是原点。若要绕其他点缩放，可使用 Shapely 的 scale 函数。

20.6 组合变换

20.6.1 变换顺序的重要性

仿射变换不可交换，顺序不同结果通常也不同：

import math
import geopandas as gpd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame({"name": ["点A"]}, geometry=[Point(1, 0)])
theta = math.radians(90)
cos_t, sin_t = math.cos(theta), math.sin(theta)

# 先旋转 90° 再平移 (5, 0)
r = gdf.geometry.affine_transform([cos_t, -sin_t, sin_t, cos_t, 0, 0])
result1 = r.affine_transform([1, 0, 0, 1, 5, 0])

# 先平移 (5, 0) 再旋转 90°
t = gdf.geometry.affine_transform([1, 0, 0, 1, 5, 0])
result2 = t.affine_transform([cos_t, -sin_t, sin_t, cos_t, 0, 0])

print("先旋转再平移：", result1.iloc[0])
print("先平移再旋转：", result2.iloc[0])

输出：

先旋转再平移： POINT (5 1)
先平移再旋转： POINT (0 6)

20.6.2 矩阵乘法实现组合变换

使用 NumPy 矩阵乘法将多个变换合并为一个矩阵：

import numpy as np
import geopandas as gpd
from shapely.geometry import Polygon
import math

def make_matrix(a, b, d, e, xoff, yoff):
    """6 参数转 3x3 仿射矩阵"""
    return np.array([[a, b, xoff], [d, e, yoff], [0, 0, 1]])

def to_params(m):
    """3x3 矩阵转 6 参数列表"""
    return [m[0,0], m[0,1], m[1,0], m[1,1], m[0,2], m[1,2]]

rect = Polygon([(0, 0), (2, 0), (2, 1), (0, 1)])
gdf = gpd.GeoDataFrame({"name": ["矩形"]}, geometry=[rect])

# 缩放 2 倍 -> 旋转 45° -> 平移 (10, 5)
scale_m = make_matrix(2, 0, 0, 2, 0, 0)
theta = math.radians(45)
rotate_m = make_matrix(
    math.cos(theta), -math.sin(theta),
    math.sin(theta), math.cos(theta), 0, 0
)
translate_m = make_matrix(1, 0, 0, 1, 10, 5)

# 从右到左应用：先缩放，再旋转，最后平移
combined = translate_m @ rotate_m @ scale_m
params = to_params(combined)
print("组合参数：", [round(p, 4) for p in params])

gdf["combined"] = gdf.geometry.affine_transform(params)
print("结果：", gdf["combined"].iloc[0])

输出：

组合参数： [1.4142, -1.4142, 1.4142, 1.4142, 10.0, 5.0]
结果： POLYGON ((10 5, 12.82842712474619 6.414213562373095, 11.414213562373096 7.828427124746191, 8.585786437626904 6.414213562373096, 10 5))

注意： 矩阵乘法顺序从右到左：translate @ rotate @ scale 表示先缩放、再旋转、最后平移。

20.7 坐标操作（Coordinate Manipulation）

20.7.1 提取坐标

import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon

gdf = gpd.GeoDataFrame(
    {"name": ["点", "线", "面"]},
    geometry=[
        Point(1, 2),
        LineString([(0, 0), (1, 1), (2, 0)]),
        Polygon([(0, 0), (3, 0), (3, 2), (0, 2)])
    ]
)

# 提取所有坐标为 DataFrame
coords = gdf.get_coordinates()
print(coords)

输出：

     x    y
0  1.0  2.0
1  0.0  0.0
1  1.0  1.0
1  2.0  0.0
2  0.0  0.0
2  3.0  0.0
2  3.0  2.0
2  0.0  2.0
2  0.0  0.0

20.7.2 修改坐标

使用 Shapely 的 set_coordinates() 直接替换坐标：

import geopandas as gpd
from shapely.geometry import Polygon
from shapely import get_coordinates, set_coordinates

rect = Polygon([(0, 0), (4, 0), (4, 3), (0, 3)])
gdf = gpd.GeoDataFrame({"name": ["矩形"]}, geometry=[rect])

# 提取并修改坐标
coords = get_coordinates(gdf.geometry.values)
new_coords = coords.copy()
new_coords[:, 0] += 100  # X 偏移
new_coords[:, 1] += 200  # Y 偏移

# 设置新坐标
new_geoms = set_coordinates(gdf.geometry.values, new_coords)
gdf["shifted"] = gpd.GeoSeries(new_geoms)
print(gdf["shifted"].iloc[0])

输出：

POLYGON ((100 200, 104 200, 104 203, 100 203, 100 200))

20.7.3 坐标精度控制

import geopandas as gpd
import numpy as np
from shapely.geometry import Point
from shapely import get_coordinates, set_coordinates

points = [
    Point(116.39128456789, 39.90735123456),
    Point(121.47365234567, 31.23041567890)
]
gdf = gpd.GeoDataFrame({"city": ["北京", "上海"]}, geometry=points)

# 四舍五入到 6 位小数（约 0.1 米精度）
coords = get_coordinates(gdf.geometry.values)
rounded = np.round(coords, decimals=6)
new_geoms = set_coordinates(gdf.geometry.values.copy(), rounded)
gdf["rounded"] = gpd.GeoSeries(new_geoms)

print("原始：", get_coordinates(gdf.geometry.values))
print("精度控制后：", get_coordinates(gdf["rounded"].values))

输出：

原始： [[116.39128457  39.90735123]
 [121.47365235  31.23041568]]
精度控制后： [[116.391285  39.907351]
 [121.473652  31.230416]]

注意： 经纬度保留 6 位小数约对应 0.1 米精度，通常满足大部分 GIS 应用需求。

20.8 坐标转换与 CRS 变换的区别

20.8.1 仿射变换 vs CRS 投影变换

对比维度 仿射变换 CRS 投影变换
数学基础 线性代数（矩阵运算） 大地测量学（椭球体投影）
变换性质 保持直线和平行性 可能引入变形
参数来源 自定义矩阵参数 标准 CRS 定义（EPSG 等）
适用场景 本地坐标调整、图形变换 地理坐标与投影坐标转换
使用方法 affine_transform() to_crs()

20.8.2 何时使用仿射变换

场景 推荐方法 原因
WGS84 转 Web 墨卡托 to_crs() 标准投影变换
CAD 本地坐标对齐到 GIS 仿射变换 需要平移、旋转、缩放对齐
建筑平面图定位 仿射变换 局部坐标到地理坐标的映射
不同投影间转换 to_crs() 涉及椭球体参数和投影公式
图形镜像翻转 仿射变换 纯几何操作
批量坐标偏移修正 仿射变换 简单的平移操作

20.9 实际应用案例

20.9.1 CAD 数据对齐到地理坐标

import math
import geopandas as gpd
from shapely.geometry import Polygon

# CAD 本地坐标系中的建筑轮廓
cad_buildings = gpd.GeoDataFrame(
    {"name": ["教学楼", "图书馆", "食堂"]},
    geometry=[
        Polygon([(100, 200), (150, 200), (150, 230), (100, 230)]),
        Polygon([(200, 300), (260, 300), (260, 340), (200, 340)]),
        Polygon([(50, 350), (90, 350), (90, 380), (50, 380)])
    ]
)

# 变换参数：缩放 0.00001，旋转 15°，平移到 (116.3, 39.9) 附近
theta = math.radians(15)
s = 0.00001
a = s * math.cos(theta)
b = -s * math.sin(theta)
d = s * math.sin(theta)
e = s * math.cos(theta)

transform_matrix = [a, b, d, e, 116.30, 39.90]
cad_buildings["geo_geom"] = cad_buildings.geometry.affine_transform(transform_matrix)

print("变换后地理坐标：")
for _, row in cad_buildings.iterrows():
    c = row["geo_geom"].centroid
    print(f"  {row['name']}: ({c.x:.6f}, {c.y:.6f})")

输出：

变换后地理坐标：
  教学楼: (116.300652, 39.902400)
  图书馆: (116.301392, 39.903651)
  食堂: (116.299721, 39.903706)

20.9.2 建筑平面图旋转与定位

import math
import numpy as np
import geopandas as gpd
from shapely.geometry import Polygon

# 建筑平面图（本地坐标，单位：米）
rooms = gpd.GeoDataFrame({
    "room": ["大厅", "办公室A", "办公室B", "会议室"],
    "geometry": [
        Polygon([(0, 0), (10, 0), (10, 8), (0, 8)]),
        Polygon([(10, 0), (16, 0), (16, 4), (10, 4)]),
        Polygon([(10, 4), (16, 4), (16, 8), (10, 8)]),
        Polygon([(0, 8), (16, 8), (16, 12), (0, 12)])
    ]
})

def make_matrix(a, b, d, e, xoff, yoff):
    return np.array([[a, b, xoff], [d, e, yoff], [0, 0, 1]])

def to_params(m):
    return [m[0,0], m[0,1], m[1,0], m[1,1], m[0,2], m[1,2]]

# 旋转 30° + 平移到投影坐标 (500000, 4400000)
theta = math.radians(30)
rotate_m = make_matrix(
    math.cos(theta), -math.sin(theta),
    math.sin(theta), math.cos(theta), 0, 0
)
translate_m = make_matrix(1, 0, 0, 1, 500000, 4400000)

combined = translate_m @ rotate_m
rooms["map_geom"] = rooms.geometry.affine_transform(to_params(combined))

print("平面图已定位到地图坐标系：")
for _, row in rooms.iterrows():
    c = row["map_geom"].centroid
    print(f"  {row['room']}: ({c.x:.2f}, {c.y:.2f})")
print(f"\n总面积: {rooms.geometry.unary_union.area:.1f} 平方米")

输出：

平面图已定位到地图坐标系：
  大厅: (500004.34, 4400006.83)
  办公室A: (500012.25, 4400003.60)
  办公室B: (500010.25, 4400007.06)
  会议室: (500003.07, 4400014.93)

总面积: 192.0 平方米

20.10 本章小结

主题 方法/函数 关键要点
仿射变换基础 [a, b, d, e, xoff, yoff] 线性变换 + 平移的组合
向量化变换 GeoSeries.affine_transform() 高性能批量变换
平移变换 [1, 0, 0, 1, xoff, yoff] 最简单的仿射变换
旋转变换 [cosθ, -sinθ, sinθ, cosθ, 0, 0] 注意旋转中心选择
缩放变换 [sx, 0, 0, sy, 0, 0] 均匀与非均匀缩放
组合变换 NumPy 矩阵乘法 顺序不可交换
坐标提取 get_coordinates() 返回 DataFrame/ndarray
坐标修改 set_coordinates() 直接替换坐标
精度控制 np.round() 6 位小数 ≈ 0.1 米
CRS 对比 to_crs() vs affine_transform() 投影变换与几何变换的区别

核心要诀：

  • 仿射变换保持直线和平行性，适用于局部坐标系调整和图形变换
  • affine_transform() 是向量化操作，大数据量时性能远优于逐要素 apply
  • 组合变换时矩阵乘法从右到左应用：最右边的变换最先执行
  • 绕任意点变换需组合平移操作，或直接使用 Shapely 的高级函数
  • 仿射变换与 CRS 投影变换有本质区别，不可混淆使用

掌握仿射变换与坐标操作后，可以灵活处理各种坐标系对齐和几何变换需求。下一章我们将进入数据可视化领域，学习 GeoPandas 的静态可视化方法，使用 Matplotlib 创建专业的地图图表。

↑