Python地理处理如何应对DICOM影像？GIS坐标转换实战技巧（附：完整代码）

wowwwai GIS研习社 · 工具流程与项目排障

引言

《Python地理处理如何应对DICOM影像？GIS坐标转换实战技巧（附：完整代码）》要解决的是一个很常见但容易踩坑的问题：DICOM影像能不能像GeoTIFF一样直接放进QGIS、ArcGIS Pro或Python GIS流程里做坐标转换、叠加和分析？

答案是：多数DICOM影像不能直接当作带地理坐标的栅格使用。DICOM通常保存的是医学影像或工业影像中的像素矩阵、像素间距、切片方向和患者坐标系信息，而不是EPSG:4326、CGCS2000、高斯投影这类GIS坐标参考系统。

本文以Python地理处理为主线，演示如何读取DICOM影像、理解DICOM坐标与GIS坐标转换的关系，并给出一套可复用的完整代码：先把DICOM转换为带局部坐标的GeoTIFF，再根据控制点或已知配准参数转换到真实GIS坐标系。

背景：DICOM影像为什么不能直接当作GIS栅格？

DICOM，全称Digital Imaging and Communications in Medicine，是医学影像领域常用的数据标准。它不仅保存像素值，还保存设备、扫描参数、切片位置、像素间距、方向矩阵等大量元数据。

GIS读者最容易误解的一点是：DICOM里的“坐标”并不等于GIS里的“地理坐标”。

• DICOM坐标：常见为患者坐标系，通常使用LPS方向，即Left、Posterior、Superior。
• GIS坐标：通常为地理坐标系或投影坐标系，例如EPSG:4326、EPSG:3857、CGCS2000、UTM等。
• GeoTIFF坐标：通过仿射变换和CRS定义像素与真实空间位置之间的关系。

因此，Python地理处理DICOM影像时，第一步不是盲目“转换坐标系”，而是先判断DICOM影像是否有足够的信息可以建立像素坐标到空间坐标的映射。

原理：DICOM坐标转换到GIS坐标的核心逻辑

在GIS中，一幅栅格影像能够正确显示，至少需要两个关键信息：

• 像素到空间坐标的仿射变换：说明第几行第几列的像素对应哪个空间位置。
• 坐标参考系统CRS：说明这些空间坐标属于哪个坐标系。

DICOM影像中常见的关键标签包括：

这些信息可以帮助我们建立DICOM内部的局部空间关系，但它们通常仍然不是地球坐标。也就是说，DICOM坐标转换到GIS坐标一般分为两步：

1. 从DICOM像素坐标转换到局部物理坐标，例如以毫米或米为单位的局部平面坐标。
2. 通过控制点、已知仿射参数或外部测量信息，把局部坐标配准到真实GIS坐标系。

如果没有控制点、外部定位文件或设备提供的地理参考信息，就不能凭空把普通DICOM影像转换成真实地理坐标。此时最多只能生成局部坐标GeoTIFF，用于测量、切片浏览或后续配准。

步骤：Python读取DICOM并生成可用于GIS处理的GeoTIFF

步骤1：准备Python环境

建议使用独立虚拟环境，避免GDAL、Rasterio和NumPy版本冲突。

pip install pydicom numpy rasterio pyproj scikit-image

如果你使用Conda环境，Rasterio通常更稳定：

conda create -n dicom-gis python=3.11
conda activate dicom-gis
conda install -c conda-forge pydicom numpy rasterio pyproj scikit-image

步骤2：读取DICOM像素和关键元数据

下面代码读取单张DICOM影像，提取像素矩阵、像素间距、影像位置和影像方向。对于GIS读者来说，这一步相当于读取一幅“还没有GIS地理参考”的栅格。

import pydicom
import numpy as np

dicom_path = "input.dcm"

ds = pydicom.dcmread(dicom_path)
arr = ds.pixel_array.astype(np.float32)

print("Rows:", ds.Rows)
print("Columns:", ds.Columns)

pixel_spacing = getattr(ds, "PixelSpacing", None)
image_position = getattr(ds, "ImagePositionPatient", None)
image_orientation = getattr(ds, "ImageOrientationPatient", None)

print("PixelSpacing:", pixel_spacing)
print("ImagePositionPatient:", image_position)
print("ImageOrientationPatient:", image_orientation)

如果PixelSpacing、ImagePositionPatient或ImageOrientationPatient不存在，说明这份DICOM影像缺少建立局部空间坐标所需的关键标签，需要通过其他资料补充，例如设备导出的定位文件、影像说明、人工控制点等。

步骤3：处理DICOM像素值缩放

许多医学DICOM影像会使用RescaleSlope和RescaleIntercept对像素值进行线性转换。例如CT影像常见的HU值就依赖这两个参数。即使本文目标是GIS坐标转换，也建议先正确处理像素值，否则后续分类、拉伸、阈值分析都会出错。

def get_scaled_dicom_array(ds):
    arr = ds.pixel_array.astype(np.float32)

    slope = float(getattr(ds, "RescaleSlope", 1.0))
    intercept = float(getattr(ds, "RescaleIntercept", 0.0))

    arr = arr * slope + intercept
    return arr

arr = get_scaled_dicom_array(ds)

步骤4：用DICOM标签构建局部仿射变换

Rasterio使用仿射变换描述像素坐标到空间坐标的关系。DICOM的ImageOrientationPatient提供了行方向和列方向的方向余弦，PixelSpacing提供像素间距。

下面代码适合处理单张二维DICOM影像，生成以DICOM患者坐标为基础的局部坐标GeoTIFF。注意这里的单位通常是毫米，不是米，也不是经纬度。

from affine import Affine

def build_affine_from_dicom(ds):
    if not hasattr(ds, "PixelSpacing"):
        raise ValueError("DICOM缺少PixelSpacing，无法计算像素尺寸。")

    if not hasattr(ds, "ImagePositionPatient"):
        raise ValueError("DICOM缺少ImagePositionPatient，无法确定影像原点。")

    if not hasattr(ds, "ImageOrientationPatient"):
        raise ValueError("DICOM缺少ImageOrientationPatient，无法确定影像方向。")

    row_spacing = float(ds.PixelSpacing[0])
    col_spacing = float(ds.PixelSpacing[1])

    origin = np.array(ds.ImagePositionPatient, dtype=np.float64)
    orientation = np.array(ds.ImageOrientationPatient, dtype=np.float64)

    row_cosines = orientation[0:3]
    col_cosines = orientation[3:6]

    x_origin = origin[0]
    y_origin = origin[1]

    x_col_step = col_cosines[0] * col_spacing
    y_col_step = col_cosines[1] * col_spacing

    x_row_step = row_cosines[0] * row_spacing
    y_row_step = row_cosines[1] * row_spacing

    transform = Affine(
        x_col_step, x_row_step, x_origin,
        y_col_step, y_row_step, y_origin
    )

    return transform

local_transform = build_affine_from_dicom(ds)
print(local_transform)

这里容易出现一个理解问题：Rasterio中的行列方向与DICOM中的行方向、列方向不是一句“左上角加像素大小”就能完全描述。对于旋转影像，必须保留方向余弦，否则输出影像在GIS软件中会发生方向错误。

步骤5：输出局部坐标GeoTIFF

如果没有真实GIS坐标系，建议不要随便写EPSG:4326或EPSG:3857。可以先输出不带CRS的GeoTIFF，或者使用工程内部自定义的局部坐标系说明。

import rasterio

output_tif = "dicom_local.tif"

with rasterio.open(
    output_tif,
    "w",
    driver="GTiff",
    height=arr.shape[0],
    width=arr.shape[1],
    count=1,
    dtype=arr.dtype,
    transform=local_transform,
    crs=None
) as dst:
    dst.write(arr, 1)

print("已输出:", output_tif)

此时生成的dicom_local.tif可以在QGIS中打开，但它还不一定能和真实地图、遥感影像或矢量边界正确叠加。它只是保留了DICOM内部的局部空间关系。

步骤6：使用控制点把DICOM影像配准到真实GIS坐标系

如果你已经知道DICOM影像上几个像素点对应的真实GIS坐标，可以通过控制点建立仿射变换。控制点至少需要3个，建议使用4个以上，并尽量分布在影像四周。

下面示例假设你已经有像素坐标与目标投影坐标的对应关系，目标坐标系为EPSG:4547。请根据你的项目改成实际EPSG代码。

import numpy as np
from affine import Affine
import rasterio

def affine_from_gcps(pixel_points, map_points):
    """
    pixel_points: [(col, row), ...]
    map_points: [(x, y), ...]
    返回从像素坐标到地图坐标的Affine
    """
    if len(pixel_points) < 3:
        raise ValueError("仿射配准至少需要3个控制点。")

    A = []
    B = []

    for (col, row), (x, y) in zip(pixel_points, map_points):
        A.append([col, row, 1, 0, 0, 0])
        A.append([0, 0, 0, col, row, 1])
        B.append(x)
        B.append(y)

    A = np.array(A, dtype=np.float64)
    B = np.array(B, dtype=np.float64)

    params, residuals, rank, s = np.linalg.lstsq(A, B, rcond=None)

    a, b, c, d, e, f = params
    return Affine(a, b, c, d, e, f), residuals

pixel_points = [
    (120, 80),
    (820, 95),
    (135, 690),
    (805, 675)
]

map_points = [
    (384512.35, 3435120.80),
    (384862.10, 3435114.25),
    (384518.90, 3434816.40),
    (384855.55, 3434822.15)
]

gis_transform, residuals = affine_from_gcps(pixel_points, map_points)

print("GIS仿射变换:", gis_transform)
print("残差:", residuals)

target_crs = "EPSG:4547"

with rasterio.open(
    "dicom_gis_registered.tif",
    "w",
    driver="GTiff",
    height=arr.shape[0],
    width=arr.shape[1],
    count=1,
    dtype=arr.dtype,
    transform=gis_transform,
    crs=target_crs
) as dst:
    dst.write(arr, 1)

print("已输出: dicom_gis_registered.tif")

这一步才是真正意义上的DICOM影像GIS坐标转换。它不是简单修改坐标系定义，而是通过控制点建立像素位置和真实地图坐标之间的关系。

步骤7：在QGIS或ArcGIS Pro中验证结果

输出GeoTIFF后，不要只看文件能否打开，还要验证空间位置是否正确。

1. 在QGIS或ArcGIS Pro中加载dicom_gis_registered.tif。
2. 加载同一坐标系下的参考数据，例如正射影像、边界线、测量点或已有栅格。
3. 检查影像四角和关键特征点是否与参考数据重合。
4. 使用识别工具查看坐标值，确认单位和范围是否符合目标投影。
5. 如果出现整体偏移，优先检查控制点输入顺序、行列号和坐标系EPSG代码。

常见坑：Python处理DICOM影像坐标转换时最容易错在哪里？

坑1：把DICOM患者坐标当成经纬度

DICOM中的ImagePositionPatient看起来像三维坐标，但它通常不是经纬度。直接把这些值写成EPSG:4326，会导致影像出现在错误位置，甚至完全不可见。

坑2：忽略LPS与GIS常用坐标轴方向差异

DICOM常见LPS坐标方向与GIS中的东、北、高并不是同一套空间定义。处理三维医学影像时，如果需要与其他空间系统结合，必须明确轴方向转换关系。

坑3：只用PixelSpacing，不用ImageOrientationPatient

如果影像存在旋转，只使用像素大小会导致方向错误。正确做法是同时使用像素间距、原点和方向余弦构建仿射变换。

坑4：以为修改CRS就完成了配准

在GIS中，“定义投影”和“坐标转换”不是一回事。给没有真实地理参考的DICOM影像强行指定EPSG，只是贴了一个标签，并不会让影像自动跑到正确位置。

坑5：控制点分布太集中

控制点如果都集中在影像一角，仿射配准结果可能在远处产生明显误差。建议控制点覆盖影像四周，并选择清晰、稳定、可重复识别的点。

坑6：没有处理多帧或多切片DICOM

有些DICOM文件包含多帧影像，或者一个序列由多张切片组成。本文示例主要面向单张二维DICOM。处理三维体数据时，还需要考虑SliceThickness、SpacingBetweenSlices、切片排序和体素坐标。

方法比较：DICOM影像进入GIS流程的几种方案

对于大多数GIS项目，如果DICOM影像本身没有地理坐标，推荐的实用路线是：先用Python读取DICOM并输出局部坐标GeoTIFF，再通过QGIS地理配准器或Python控制点脚本完成真实GIS坐标转换。

检查清单：开始转换前先确认这些信息

• 是否确认DICOM影像用途？是展示、测量、二维配准，还是三维分析？
• 是否存在PixelSpacing？没有它就难以确认像素实际尺寸。
• 是否存在ImagePositionPatient和ImageOrientationPatient？没有它们就难以构建可靠局部坐标。
• 是否知道目标GIS坐标系？例如EPSG:4326、EPSG:4547或当地工程坐标系。
• 是否有控制点、测量点、参考影像或外部定位文件？没有这些信息就不能完成真实地理配准。
• 控制点是否分布均匀？至少3个，实际项目建议4个以上。
• 是否检查了行列号顺序？像素点通常写作col、row，而数组索引通常是row、col。
• 是否在QGIS或ArcGIS Pro中与参考数据叠加验证？不要只依赖代码运行成功。
• 是否注意隐私与合规？医学DICOM可能包含患者信息，分享或发布前应清理敏感元数据。

FAQ

DICOM影像可以直接转换成GeoTIFF吗？

可以转换成GeoTIFF格式，但不代表一定有真实GIS坐标。如果DICOM只包含患者坐标或局部坐标，转换后的GeoTIFF也只能表示局部空间关系。要进入真实GIS坐标系，还需要控制点或外部地理参考信息。

Python地理处理DICOM影像时，应该使用哪个库？

读取DICOM建议使用pydicom，处理数组使用NumPy，输出GeoTIFF建议使用Rasterio。如果涉及坐标系转换，可以使用pyproj；如果涉及更复杂的影像重采样和配准，可结合GDAL、scikit-image或SimpleITK。

DICOM坐标转换到GIS坐标，为什么不能只写一个EPSG代码？

因为EPSG代码只是说明坐标属于哪个坐标参考系统，并不会自动改变影像每个像素的位置。如果像素到地图坐标的仿射关系错误，写入正确EPSG也无法让影像正确叠加。

ImagePositionPatient和ImageOrientationPatient一定存在吗？

不一定。不同设备、不同DICOM类型、不同导出方式可能保留的标签不同。如果缺少这些标签，需要检查原始序列、导出设置，或通过人工控制点重新建立空间关系。

控制点配准和DICOM内部坐标哪个更可靠？

这取决于数据来源。如果DICOM内部空间标签完整且与项目空间系统存在明确转换关系，可以使用内部坐标。如果目标是叠加到地图或工程坐标，通常仍需要控制点或外部定位信息来验证和校正。

多张DICOM切片如何处理成GIS数据？

多切片DICOM需要先按切片位置排序，结合PixelSpacing、SliceThickness或SpacingBetweenSlices建立三维体素坐标。若只需二维GIS叠加，可以选择目标切片或投影切片输出GeoTIFF；若要三维分析，则应使用体数据处理流程，再考虑与GIS或三维场景配准。

结论

Python地理处理DICOM影像的关键，不是把DICOM文件简单另存为GeoTIFF，而是先分清DICOM坐标、局部物理坐标和真实GIS坐标之间的关系。

实战中可以按三步走：第一，用pydicom读取像素和关键元数据；第二，根据PixelSpacing、ImagePositionPatient和ImageOrientationPatient生成局部坐标GeoTIFF；第三，在具备控制点或外部定位信息时，将影像配准到目标GIS坐标系。

只要记住一点，就能避开大多数错误：没有真实地理参考信息时，不要强行指定EPSG；有控制点或测量数据时，才可以把DICOM影像可靠地纳入QGIS、ArcGIS Pro和Python GIS分析流程。