当前位置：首页 > news >正文

【图像处理基石】GIS图像处理入门：4个核心算法与Python实现（附完整代码）

news 2025/10/10 11:50:01

在这里插入图片描述

最近在项目里需要处理遥感影像、电子地图这类GIS数据，踩了不少“入门坑”后，整理了这篇偏向实操的入门博客——从工具搭建到核心算法实现，全程用Python代码落地，适合快速上手。

1. 前言：先搞懂“GIS图像处理”到底是什么？

和普通的图片（比如JPG、PNG）不同，GIS图像（遥感影像、DEM、电子地图等）自带地理坐标和投影信息——比如一张卫星图，每个像素都对应地球表面的具体经纬度。

我们处理GIS图像，本质是解决这些问题：

把“raw数据”（比如卫星传回的DN值）转成可解读的物理量（反射率、温度）；
消除图像畸变（比如卫星姿态、地形导致的几何偏差）；
提取有用信息（比如裁剪出城市范围、识别耕地面积）。

常用场景：遥感监测、城市规划、灾害评估、农业估产，比如最近很火的“高分卫星影像解译”，核心就是GIS图像处理。

2. 必备工具：3分钟搭好Python环境

GIS图像处理不像普通CV用OpenCV就行，需要专门处理“地理信息”的库。核心工具就3个，直接用pip安装（如果GDAL装不上，建议用conda conda install gdal）：

# 核心库：读写GIS影像、处理地理坐标
pip install gdal
# 数值计算：处理图像像素数组
pip install numpy
# 可视化：显示图像、画结果图
pip install matplotlib

简单介绍下这3个库的分工：

GDAL/OGR：GIS领域的“瑞士军刀”，支持TIFF、SHP（矢量文件）等90%以上的GIS格式，能读取图像的投影、地理坐标；
NumPy：把GIS图像转成数组（比如256*256的像素矩阵），做数值计算；
Matplotlib：替代OpenCV的imshow，方便显示带地理信息的图像。

3. 核心入门算法：原理+代码落地

这部分是重点，我选了4个“入门必会”的算法——从“读图像”到“预处理”，覆盖80%的基础场景，每个算法都附可运行代码。

3.1 算法1：GIS图像读取与显示（最基础但有坑）

普通图像用OpenCV的cv2.imread()就能读，但GIS图像不行——会丢失地理坐标！必须用GDAL读。

原理

用GDAL的Open()函数打开影像文件（比如TIFF）；
获取图像的“地理变换参数”（知道每个像素对应哪个经纬度）和“投影信息”（比如WGS84、UTM）；
把图像数据转成NumPy数组，用Matplotlib显示。

代码实现

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt# 解决GDAL中文路径问题
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")def read_gis_image(image_path):"""读取GIS图像（支持TIFF等格式）:param image_path: 图像路径（比如"test.tif"）:return: 像素数组、地理变换参数、投影信息"""# 1. 打开影像文件dataset = gdal.Open(image_path)if dataset is None:raise Exception(f"无法打开图像：{image_path}")# 2. 获取图像基本信息width = dataset.RasterXSize  # 图像宽度（像素数）height = dataset.RasterYSize  # 图像高度（像素数）band_count = dataset.RasterCount  # 波段数（比如RGB是3波段，遥感影像可能有10+波段）# 3. 获取地理信息（关键！普通图像没有这个）geo_transform = dataset.GetGeoTransform()  # 地理变换参数projection = dataset.GetProjection()  # 投影信息（WKT格式）# 4. 读取像素数据（这里以单波段为例，多波段可循环读取）band = dataset.GetRasterBand(1)  # 读取第1波段（索引从1开始，不是0！）pixel_array = band.ReadAsArray(0, 0, width, height)  # 转成NumPy数组# 5. 释放资源（重要，避免内存泄漏）dataset = Noneband = Nonereturn pixel_array, geo_transform, projectiondef show_gis_image(pixel_array, title="GIS图像"):"""显示GIS图像"""plt.figure(figsize=(8, 6))# 用灰度显示（多波段可调整为RGB组合）plt.imshow(pixel_array, cmap="gray")plt.title(title)plt.axis("off")  # 隐藏坐标轴（GIS图像看地理意义，不看像素坐标）plt.show()# ------------------- 测试代码 -------------------
if __name__ == "__main__":image_path = "你的GIS图像路径.tif"  # 替换成自己的TIFF文件pixel_array, geo_transform, projection = read_gis_image(image_path)# 打印关键信息print(f"图像尺寸：{pixel_array.shape}")  # (高度, 宽度)print(f"地理变换参数：{geo_transform}")  # (左上角x, 水平分辨率, 0, 左上角y, 0, 垂直分辨率)print(f"投影信息：{projection[:100]}...")  # 打印前100个字符（WKT格式很长）# 显示图像show_gis_image(pixel_array, title="原始GIS图像")

关键注意点

GDAL读取波段时，索引从1开始（不是Python的0！）；
geo_transform参数很重要：(x0, dx, 0, y0, 0, dy)，其中x0,y0是左上角像素的地理坐标，dx是水平分辨率（每个像素代表多少米/度），dy通常是负数（因为图像坐标系从上到下，地理坐标系从下到上）。

3.2 算法2：辐射定标（从“数字”到“物理量”）

卫星传回来的原始数据是“DN值”（0-255的数字），没有物理意义——我们需要把它转成“反射率”（衡量物体反射太阳光的能力），这个过程就是辐射定标。

原理

核心公式（线性定标，大部分卫星影像适用）：
反射率 = 增益（Gain）× DN值 + 偏移（Bias）
其中Gain和Bias从卫星影像的元数据文件（比如MTL.txt）中获取，不同波段的参数不同。

代码实现

def radiometric_calibration(pixel_array, gain, bias):"""辐射定标：DN值转反射率:param pixel_array: 原始DN值数组:param gain: 增益（从元数据获取）:param bias: 偏移（从元数据获取）:return: 反射率数组（值范围通常0-1）"""# 应用定标公式reflectance = gain * pixel_array + bias# 裁剪异常值（反射率不会小于0，通常不会大于1.5）reflectance = np.clip(reflectance, 0, 1.5)# 归一化到0-1（方便显示）reflectance = (reflectance - np.min(reflectance)) / (np.max(reflectance) - np.min(reflectance))return reflectance# ------------------- 测试代码 -------------------
if __name__ == "__main__":# 1. 先读取原始图像（DN值）image_path = "你的原始影像.tif"dn_array, geo_transform, projection = read_gis_image(image_path)# 2. 从元数据获取Gain和Bias（这里以Landsat 8为例，实际要换自己的参数）gain = 0.0001  # 示例值，实际从MTL.txt找"REFLECTANCE_MULT_BAND_x"bias = -0.1    # 示例值，实际从MTL.txt找"REFLECTANCE_ADD_BAND_x"# 3. 执行定标reflectance_array = radiometric_calibration(dn_array, gain, bias)# 4. 显示对比plt.figure(figsize=(12, 6))plt.subplot(1, 2, 1)plt.imshow(dn_array, cmap="gray")plt.title("原始DN值图像")plt.axis("off")plt.subplot(1, 2, 2)plt.imshow(reflectance_array, cmap="gray")plt.title("辐射定标后的反射率图像")plt.axis("off")plt.show()

关键注意点

不同卫星（Landsat、Sentinel、高分）的定标公式可能不同，必须查对应卫星的官方文档；
如果没有元数据，可从卫星官网（比如USGS）下载对应影像的MTL文件，里面有现成的Gain和Bias。

3.3 算法3：几何校正（消除图像“歪掉”的问题）

卫星拍摄时，会因为姿态倾斜、地形起伏导致图像“畸变”——比如地面上的正方形，在图像上变成梯形。几何校正就是把畸变的图像校正到正确的地理坐标上。

原理（多项式校正，入门常用）

选“控制点”：在畸变图像和“标准地图”（比如谷歌地图）上找相同的点（比如路口、山顶），记录它们的坐标；
用多项式拟合控制点，计算“校正矩阵”；
对原始图像的每个像素，用校正矩阵计算它在标准坐标系中的位置，再通过“重采样”填充像素值。

代码实现

from osgeo import gdal, osrdef geometric_correction(input_path, output_path, control_points):"""几何校正（多项式校正）:param input_path: 输入畸变图像路径:param output_path: 输出校正后图像路径:param control_points: 控制点列表，格式：[(畸变x1, 畸变y1, 标准x1, 标准y1), ...]"""# 1. 打开输入影像input_ds = gdal.Open(input_path)if input_ds is None:raise Exception(f"无法打开输入图像：{input_path}")# 2. 定义输出影像的投影（这里用WGS84投影，可根据需求修改）output_srs = osr.SpatialReference()output_srs.ImportFromEPSG(4326)  # 4326 = WGS84坐标系（经纬度）output_wkt = output_srs.ExportToWkt()# 3. 准备控制点（分离畸变坐标和标准坐标）src_points = np.array([(cp[0], cp[1]) for cp in control_points], dtype=np.float32)dst_points = np.array([(cp[2], cp[3]) for cp in control_points], dtype=np.float32)# 4. 计算多项式校正矩阵（2次多项式，需要至少6个控制点）if len(control_points) < 6:raise Exception("2次多项式校正需要至少6个控制点")poly_matrix = cv2.getPerspectiveTransform(src_points[:4], dst_points[:4])  # 简化用透视变换（4点），实际用多项式# 5. 创建输出影像（和输入影像同尺寸）driver = gdal.GetDriverByName("GTiff")output_ds = driver.Create(output_path,input_ds.RasterXSize,input_ds.RasterYSize,1,gdal.GDT_Float32  # 浮点型，避免精度丢失)# 6. 设置输出影像的地理变换和投影# 这里简化：假设标准坐标的分辨率是0.0001度（根据实际调整）output_geo = (min(dst_points[:, 0]), 0.0001, 0,max(dst_points[:, 1]), 0, -0.0001)output_ds.SetGeoTransform(output_geo)output_ds.SetProjection(output_wkt)# 7. 重采样（双线性插值，平衡速度和精度）gdal.ReprojectImage(input_ds,output_ds,input_ds.GetProjection(),output_wkt,gdal.GRA_Bilinear  # 重采样方法：GRA_NearestNeighbor（最快）、GRA_Bilinear（较优）)# 8. 释放资源input_ds = Noneoutput_ds = Noneprint(f"几何校正完成，输出路径：{output_path}")# ------------------- 测试代码 -------------------
if __name__ == "__main__":# 示例控制点（需要自己在图像上选！这里是模拟值）control_points = [(100, 200, 116.39, 39.90),  # (畸变图像x,y, 标准经纬度lon,lat)(200, 200, 116.40, 39.90),(200, 300, 116.40, 39.89),(100, 300, 116.39, 39.89),(150, 250, 116.395, 39.895),(250, 250, 116.405, 39.895)]# 执行校正geometric_correction(input_path="畸变影像.tif",output_path="校正后影像.tif",control_points=control_points)

关键注意点

控制点越多、越均匀，校正精度越高；
重采样方法选择：追求速度用“最近邻”，追求精度用“双线性插值”或“三次卷积”；
实际项目中，常用“地面控制点（GCP）”工具（比如QGIS）手动选点，比代码硬编码更灵活。

3.4 算法4：图像裁剪（提取你关心的区域）

拿到一张大的遥感影像，我们通常只需要其中一小块（比如某个城市、某个区县）——这时候就需要图像裁剪，支持按“矢量边界”（SHP文件）或“坐标范围”裁剪。

原理

按矢量边界：用SHP文件（比如城市边界）作为“裁剪框”，保留框内的像素；
按坐标范围：指定经纬度范围（比如lon:116.3-116.4, lat:39.8-39.9），裁剪出对应区域。

代码实现（按矢量边界裁剪）

def crop_by_shapefile(input_path, shapefile_path, output_path):"""按矢量边界（SHP文件）裁剪GIS图像:param input_path: 输入图像路径:param shapefile_path: SHP文件路径（裁剪边界）:param output_path: 输出裁剪后图像路径"""# 1. 打开输入影像input_ds = gdal.Open(input_path)if input_ds is None:raise Exception(f"无法打开输入图像：{input_path}")# 2. 调用GDAL的Warp函数进行裁剪（核心！）gdal.Warp(destNameOrDestDS=output_path,srcDSOrSrcDSTab=input_ds,cutlineDSName=shapefile_path,  # 裁剪边界的SHP文件cropToCutline=True,  # 严格按边界裁剪dstNodata=np.nan,  # 边界外的像素设为NaN（透明）format="GTiff"  # 输出格式)# 3. 释放资源input_ds = Noneprint(f"按SHP裁剪完成，输出路径：{output_path}")# 按坐标范围裁剪（补充）
def crop_by_coords(input_path, output_path, min_lon, max_lon, min_lat, max_lat):"""按经纬度范围裁剪:param min_lon/max_lon: 最小/最大经度:param min_lat/max_lat: 最小/最大纬度"""input_ds = gdal.Open(input_path)# 定义裁剪范围（格式：xmin, ymin, xmax, ymax）crop_extent = (min_lon, min_lat, max_lon, max_lat)gdal.Warp(destNameOrDestDS=output_path,srcDSOrSrcDSTab=input_ds,outputBounds=crop_extent,  # 坐标范围dstNodata=np.nan,format="GTiff")input_ds = Noneprint(f"按坐标裁剪完成，输出路径：{output_path}")# ------------------- 测试代码 -------------------
if __name__ == "__main__":# 1. 按SHP裁剪crop_by_shapefile(input_path="原始影像.tif",shapefile_path="北京市边界.shp",output_path="北京市裁剪影像.tif")# 2. 按坐标裁剪（示例：北京某区域）crop_by_coords(input_path="原始影像.tif",output_path="北京某区域裁剪影像.tif",min_lon=116.3, max_lon=116.4,min_lat=39.8, max_lat=39.9)

关键注意点

SHP文件和输入影像的投影必须一致（比如都用WGS84），否则会裁剪失败；
如果投影不一致，需要先对SHP文件做“投影转换”（用QGIS或GDAL的ogr2ogr工具）。

4. 进阶学习方向

入门后如果想深入，可以往这3个方向走：

图像增强：用直方图均衡化、高斯滤波等方法提升图像质量（推荐库：scikit-image）；
图像分类：把图像分成“耕地、建筑、水体”等类别（监督分类：SVM、随机森林；非监督分类：K-Means）；
深度学习：用CNN做遥感影像解译（比如用U-Net分割建筑区域，推荐框架：TensorFlow/PyTorch）。

推荐学习资源：

GDAL官方文档：https://gdal.org/ （最权威的GIS处理指南）；
USGS遥感教程：https://www.usgs.gov/landsat-missions/landsat-science （学习辐射定标、几何校正的原理）；
QGIS软件：免费开源的GIS工具，可可视化调试（比如手动选控制点、查看影像投影）。