【Python&RS】基于Python栅格数据/遥感影像投影转换

RS迷途小书童

已于 2023-09-12 11:41:34 修改

阅读量828

点赞数 3

分类专栏： Python&Remote Sensing 文章标签： python 开发语言数据分析经验分享

于 2023-06-28 10:17:14 首次发布

本文链接：https://blog.csdn.net/m0_56729804/article/details/131430570

版权

Python&Remote Sensing 专栏收录该内容

22 篇文章 21 订阅

订阅专栏

有ENVI出现的地方，就一定会有Python的身影，为了解放双手方便批量投影转换，最近同步研究了一下如何利用Python实现遥感影像的投影转换，主打的就是一个懒。

之前还分享过矢量数据的投影转换，这样矢量+栅格的全家桶不就有了嘛。感兴趣的可以自己查看：【Python&GIS】矢量数据投影转换（WGS84转地方坐标系）

一、查看影像信息

为了水文章，这个其实不太重要，但有肯定更好，因为有些数据在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函数，是不是有点熟悉。之前裁剪也是用的这个函数，不得不说这个函数是真的强大，之前有文章介绍了其中的函数，大家有兴趣可以去看下，同时给个赞吧！【Python&RS】GDAL批量裁剪遥感影像/栅格数据

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

三、回调函数（没啥用）

上一篇博文已经介绍过了，就两点：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)