栅格数据投影转换

Posted 2022-03-08 卡尔曼和玻尔兹曼谁曼

栅格数据投影转换

作者：阿振

邮箱：[e-mail]

博客：https://blog.csdn.net/theonegis/article/details/80089375

修改时间：2018-06-01

声明：本文为博主原创文章，转载请注明原文出处

---

使用GDAL提供的命令行工具进行转换

GDAL提供了gdalwarp命令可以方便地让我们进行影像拼接，重投影，裁剪，格式转换等功能

比如，我们需要将MODIS数据的Sinusoidal投影转为UTM投影，我们可以这样操作。

我需要转换的地区位于UTM的49度带内，我查看了一下其EPSG的编码为：EPSG:32649（WGS 84 / UTM zone 49N）

注：推荐大家一个网站，可以查阅各种投影的定义：http://spatialreference.org

然后，终端中执行如下命令：

gdalinfo MOD09A1.A2017361.h28v06.006.2018005034659.hdf （用于查看MODIS数据中的波段名称与地址，这里我们只转换第一波段）

gdalwarp -t_srs EPSG:32649 HDF4_EOS:EOS_GRID:"MOD09A1.A2017361.h28v06.006.2018005034659.hdf":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01 MODSI_WARP_32649.tif（-t_srs参数用于指定输出投影信息，可以是EPSG，或者OGC WKT，或者PROJ4格式，后面分别是输入数据和输出数据文件名）

使用代码进行转换

使用命令行转换，当然有两种方法啦：

第一，直接在代码中调用命令行工具的接口（比较懒的人，像我，当然直接用第一种啦，有现成的工具为什么不用）；

第二，自己做投影转换之后的坐标计算，主要是计算重投影之后的GeoTransform参数，有了GeoTransform参数以及投影的定义，我们就可以通过SetGeoTransform()和SetProjection()投影转换了.

下面我给出具体的实现代码：

第一种方法直接调用gdal.Warp()方法，该方法其实就是对gdalwarp命令的封装，第一个参数是输出文件，第二个参数是输入文件或者输入的Dataset，后面的都是可选参数（dstSRS参数指定输出投影）

from osgeo import gdal

root_ds = gdal.Open('MOD09A1.A2017361.h28v06.006.2018005034659.hdf')
# 返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息
# tuple中的第一个元素描述的是数据子集的全路径
ds_list = root_ds.GetSubDatasets()

# 取出第1个数据子集（MODIS反射率产品的第一个波段）进行转换
# 第一个参数是输出数据，第二个参数是输入数据，后面可以跟多个可选项
gdal.Warp('reprojection.tif', ds_list[0][0], dstSRS='EPSG:32649')

# 关闭数据集
root_ds = None

在介绍第二种方法之前，我们有必要回忆一下之前说过的GDAL反射变换的六参数模型：

放射变换使用如下的公式表示栅格图上坐标和地理坐标的关系：

$$\\beginmatrix X_geo = GT(0) + X_pixel * GT(1) + Y_line * GT(2) \\\\ Y_geo = GT(3) + X_pixel * GT(4) + Y_line * GT(5) \\\\ \\endmatrix$$
（ $X_ge0$, $Y_ge0$）表示对应于图上坐标（ $X_pixel$, $Y_line$）的实际地理坐标。对一个上北下南的图像，GT(2)和GT(4)等于0， GT(1)是像元的宽度, GT(5)是像元的高度的相反数。（GT(0),GT(3)）坐标对表示左上角像元的左上角坐标。

通过这个放射变换，我们可以得到图上所有像元对应的地理坐标。

好了，所以我们需要计算对于上面的六参数，我们主要需要计算重投影以后图像左上角的坐标（最小的X坐标值和最大的Y坐标值），这个转换我们可以通过osr.CoordinateTransformation类进行，下面给出实现代码：

from osgeo import gdal
from osgeo import osr

# root_ds = gdal.Open('/Users/tanzhenyu/Resources/DataWare/MODIS/MOD09A1.A2017361.h28v06.006.2018005034659.hdf')
# # 返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息
# # tuple中的第一个元素描述的是数据子集的全路径
# ds_list = root_ds.GetSubDatasets()
#
# # 取出第1个数据子集（MODIS反射率产品的第一个波段）进行转换
# # 第一个参数是输出数据，第二个参数是输入数据，后面可以跟多个可选项
# gdal.Warp('reprojection.tif', ds_list[0][0], dstSRS='EPSG:32649')
#
# # 关闭数据集
# root_ds = None


def reproject(src_file, dst_file, p_width, p_height, epsg_to):
    """
    :param src_file: 输入文件
    :param dst_file: 输出文件
    :param p_width: 输出图像像素宽度
    :param p_height: 输出图像像素高度
    :param epsg_to: 输出图像EPSG坐标代码
    :return:
    """
    # 首先，读取输入数据，然后获得输入数据的投影，放射变换参数，以及图像宽高等信息
    src_ds = gdal.Open(src_file)
    src_srs = osr.SpatialReference()
    src_srs.ImportFromWkt(src_ds.GetProjection())

    srs_trans = src_ds.GetGeoTransform()
    x_size = src_ds.RasterXSize
    y_size = src_ds.RasterYSize
    d_type = src_ds.GetRasterBand(1).DataType

    # 获得输出数据的投影，建立两个投影直接的转换关系
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(epsg_to)
    tx = osr.CoordinateTransformation(src_srs, dst_srs)

    # 计算输出图像四个角点的坐标
    (ulx, uly, _) = tx.TransformPoint(srs_trans[0], srs_trans[3])
    (urx, ury, _) = tx.TransformPoint(srs_trans[0] + srs_trans[1] * x_size, srs_trans[3])
    (llx, lly, _) = tx.TransformPoint(srs_trans[0], srs_trans[3] + srs_trans[5] * y_size)
    (lrx, lry, _) = tx.TransformPoint(srs_trans[0] + srs_trans[1] * x_size + srs_trans[2] * y_size,
                                      srs_trans[3] + srs_trans[4] * x_size + srs_trans[5] * y_size)

    min_x = min(ulx, urx, llx, lrx)
    max_x = max(ulx, urx, llx, lrx)
    min_y = min(uly, ury, lly, lry)
    max_y = max(uly, ury, lly, lry)

    # 创建输出图像，需要计算输出图像的尺寸（重投影以后图像的尺寸会发生变化）
    driver = gdal.GetDriverByName('GTiff')
    dst_ds = driver.Create(dst_file,
                           int((max_x - min_x) / p_width),
                           int((max_y - min_y) / p_height),
                           1, d_type)
    dst_trans = (min_x, p_width, srs_trans[2],
                 max_y, srs_trans[4], -p_height)

    # 设置GeoTransform和Projection信息
    dst_ds.SetGeoTransform(dst_trans)
    dst_ds.SetProjection(dst_srs.ExportToWkt())
    # 进行投影转换
    gdal.ReprojectImage(src_ds, dst_ds,
                        src_srs.ExportToWkt(), dst_srs.ExportToWkt(),
                        gdal.GRA_Bilinear)
    dst_ds.GetRasterBand(1).SetNoDataValue(0)  # 设置NoData值
    dst_ds.FlushCache()

    del src_ds
    del dst_ds


if __name__ == '__main__':
    src_file = 'HDF4_EOS:EOS_GRID:"MOD09A1.A2017361.h28v06.006.2018005034659.hdf":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01'
    dst_file = 'reprojection.tif'
    reproject(src_file, dst_file, 450, 450, 32649)