【Python&RS】基于Python栅格数据/遥感影像投影转换
有ENVI出现的地方，就一定会有Python的身影，为了解放双手方便批量投影转换，最近同步研究了一下如何利用Python实现遥感影像的投影转换，主打的就是一个懒。
之前还分享过矢量数据的投影转换，这样矢量+栅格的全家桶不就有了嘛。感兴趣的可以自己查看：
一、查看影像信息
为了水文章，这个其实不太重要，但有肯定更好，因为有些数据在ENVI中能看到地理坐标系，但GDAL却识别不到坐标系，这就会导致投影转换的失败。所以最好还是提前看一看，GDAL能不能读出投影信息吧。
这里需要注意ds_geo和ds_prj两个输出，仿射地理变换参数和投影坐标系一样重要！
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) # 以数组的形式读取整个数据集二、转换函数
这里用到了GDAL中的Warp函数，是不是有点熟悉。之前裁剪也是用的这个函数，不得不说这个函数是真的强大，之前有文章介绍了其中的函数，大家有兴趣可以去看下，同时给个赞吧！
Warp函数中的dstSRS参数就是目标的投影，这里采用EPSG编码，32651表示UTM/WGS84 51N投影坐标系。其他坐标系对应的代码可以查看
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三、回调函数（没啥用）
上一篇博文已经介绍过了，就两点：1.水字数，2.记录进度。省的看着代码发呆。
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)
无论是用ENVI、ArcGIS这类软件进行的投影转换还是用Python的三方库，其实都是用一些内置的参数就行转换的，所以只适用一些常见的坐标系之间的转换。最精确的肯定还是四参数、七参数进行转换，所以只能说能转成功就皆大欢喜，转不成功就算了，强扭瓜不甜。
如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！
