遥感影像处理中，投影转换是一项常见步骤。ENVI能处理这类任务，但Python的功能更为出色，它能批量完成投影转换，大大减轻了我们的工作量。那么，Python是如何实现这一功能的？咱们来探究一下。
查看影像信息
查看影像资料看似简单，实则关键。在ENVI中，某些数据能明确显示地理坐标，但GDAL却无法识别，这可能导致投影转换直接失败。因此，最好事先确认GDAL能否读取投影信息。以前，我处理一批影像时，未先查看信息，导致转换时出现问题，浪费了大量时间。
在检查过程中，ds_geo和ds_prj这两个输出结果需特别关注。仿射地理变换的参数和投影坐标系的设定同样关键。这就像盖房子，如果基础参数设置错误，后续工程将难以顺利进行。投影转换亦是如此。
转换函数
这次操作使用了GDAL库中的Warp函数，这个函数大家应该不陌生，之前进行裁剪时也用它来过。这个函数功能强大，很多文章都有对其功能的介绍。我在之前的【Python&RS】GDAL批量裁剪遥感影像/栅格数据文章中，详细解释了这个函数的使用方法。有兴趣的朋友可以查阅一下，觉得有帮助的话，请给我点个赞。
def Get_data(filepath):
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
print("影像的宽度为：" + str(ds_width))
print("影像的高度为：" + str(ds_height))
print("仿射地理变换参数为：" + str(ds_geo))
print("投影坐标系为：" + str(ds_prj))
return ds_geo
# data = ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集在Warp函数中，dstSRS参数指的是目标投影坐标系统，我们使用EPSG编码来标识，比如32651代表的是UTM/WGS84 51N投影坐标系。若要了解其他坐标系的编码，可以查阅【Python&GIS】中关于如何通过经纬度生成矢量点文件的资料。熟练掌握这些编码及函数操作，投影转换将变得简单易行。
回调函数（没啥用）
前文已经提及了回调函数的相关内容，它主要承担两个职能。首先，它能够起到增加文字量的作用，当然这只是一个玩笑的说法；其次，它有助于记录操作进程。在进行代码投影转换时，若只是静静等待，确实会感到无聊。而引入回调函数后，我们便能够实时查看转换的进度，从而减轻了等待的枯燥感。以前没有使用回调函数时，代码毫无反应，让人误以为出现了错误。但自从添加了回调函数，我们便对整个过程有了更清晰的把握。
尽管对投影转换的结果影响不大，但它的存在让人感到安心，仿佛每一步都有明确的路径。这就像在长跑中，清楚自己已经跑了多远，以及距离终点还有多远，这样的认知让人更有毅力继续前进。
完整代码
def Projection_Transform(path_input1, path_output1):
"""
:param path_input1: 待转换的数据路径
:param path_output1: 转换后的数据路径
:return: None
"""
ds_input = gdal.Open(path_input1, gdal.GA_ReadOnly) # 打开数据集dataset
ds = gdal.Warp(path_output1, ds_input, dstSRS="EPSG:32651", resampleAlg='cubic', callback=Show_Progress)
# 输出路径、输入影像、目标投影、重采样方法(最邻近内插双线性内插三次卷积等)、回调函数
ds.BuildOverviews(overviewlist=[2, 4, 8, 16])
# 创建金字塔
ds = None无论是借助ENVI或ArcGIS等软件进行坐标投影的转换，还是采用Python的相关库，本质上都是通过预设的参数来实现的。这类方法通常适用于一些常见的坐标系转换。但要达到最高精度，还是得使用四参数或七参数转换。当然，有时这些方法可能无法成功，但也不必过分烦恼，毕竟过分追求未必有所得。
学习Python或RS时，若遇到难题，可直接留言询问。若对批量处理有兴趣，不妨告诉我，我会把相关代码拿出来供大家共同学习。大家多在实践中探索，时间久了，自然能熟练运用其中的方法。
过往经验分享
以前我提到过矢量数据的投影转换技术。现在有了矢量数据与栅格数据的投影转换方法，这不就是一套完整的工具了吗？记得之前帮朋友处理了一大批矢量数据，我用Python完成了投影转换，他夸我效率很高。掌握了这些技巧，无论是工作还是学习，都能节省很多时间和精力。
def Show_Progress(percent, msg, tag):
"""
:param percent: 进度，0~1
:param msg:
:param tag:
:return:
"""
print("开始进程......")
if 25 <= percent*100 <= 26:
print("进度：" + "%.2f" % (percent*100) + "%")
if 50 <= percent*100 <= 51:
print("进度：" + "%.2f" % (percent*100) + "%")
if 75 <= percent*100 <= 76:
print("进度：" + "%.2f" % (percent*100) + "%")将矢量与栅格数据相结合，在应用中效果显著。例如，在进行地理信息分析时，这两种数据类型都会用到。掌握这两种处理方式，操作起来会更加得心应手。希望各位能多加实践，体会其中的便捷。
# -*- coding: utf-8 -*-
"""
@Time ： 2023/6/26 17:03
@Auth ： RS迷途小书童
@File ：Raster Data Projection Transform.py
@IDE ：PyCharm
@Purpose ：栅格数据投影/坐标转换
"""
from osgeo import gdal
def Get_data(filepath):
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
print("影像的宽度为：" + str(ds_width))
print("影像的高度为：" + str(ds_height))
print("仿射地理变换参数为：" + str(ds_geo))
print("投影坐标系为：" + str(ds_prj))
return ds_geo
# data = ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集
def Show_Progress(percent, msg, tag):
"""
:param percent: 进度，0~1
:param msg:
:param tag:
:return:
"""
print("开始进程......")
if 25 <= percent*100 <= 26:
print("进度：" + "%.2f" % (percent*100) + "%")
if 50 <= percent*100 <= 51:
print("进度：" + "%.2f" % (percent*100) + "%")
if 75 <= percent*100 <= 76:
print("进度：" + "%.2f" % (percent*100) + "%")
def Projection_Transform(path_input1, path_output1):
"""
:param path_input1: 待转换的数据路径
:param path_output1: 转换后的数据路径
:return: None
"""
ds_input = gdal.Open(path_input1, gdal.GA_ReadOnly) # 打开数据集dataset
ds = gdal.Warp(path_output1, ds_input, dstSRS="EPSG:32651", resampleAlg='cubic', callback=Show_Progress)
# 输出路径、输入影像、目标投影、重采样方法(最邻近内插双线性内插三次卷积等)、回调函数
ds.BuildOverviews(overviewlist=[2, 4, 8, 16])
# 创建金字塔
ds = None
if __name__ == "__main__":
path_input = "B:/sharpening2_proj"
path_output = r"B:/sharpening1.tif"
Get_data(path_input)
Projection_Transform(path_input, path_output)
学习建议
学习Python处理遥感图像投影转换时，需多加练习和探索。可从基础的图像开始，逐步掌握查阅数据、运用转换工具等操作。遇到难题不要急躁，将错误信息记录下来，上网搜索或与同行交流，多数问题都能找到解决之道。
平日里应多浏览相关资料和文章，比如我之前撰写过的，通过学习别人的经验，我们也能更快地提升自己。通过实践中的不断总结，逐步积累，相信大家最终都能熟练地运用Python进行遥感影像的投影转换。大家是否有意尝试使用Python进行投影转换？若觉得这篇文章有所帮助，请不要忘记点赞和转发！