栅格数据投影转换

本是古典何须时尚 2022-05-23 04:35 598阅读 0赞

# 栅格数据投影转换 #

作者：阿振

邮箱：tanzhenyugis@163.com

博客：[https://blog.csdn.net/theonegis/article/details/80089375][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][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反射变换的六参数模型：

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

Xgeo=GT(0)\+Xpixel∗GT(1)\+Yline∗GT(2)Ygeo=GT(3)\+Xpixel∗GT(4)\+Yline∗GT(5)  X g e o = G T ( 0 ) + X p i x e l ∗ G T ( 1 ) + Y l i n e ∗ G T ( 2 ) Y g e o = G T ( 3 ) + X p i x e l ∗ G T ( 4 ) + Y l i n e ∗ G T ( 5 )

（  Xge0  X g e 0  ,  Yge0  Y g e 0  ）表示对应于图上坐标（  Xpixel  X p i x e l  ,  Yline  Y l i n e  ）的实际地理坐标。对一个上北下南的图像，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)

[https_blog.csdn.net_theonegis_article_details_80089375]: https://blog.csdn.net/theonegis/article/details/80089375
[http_spatialreference.org]: http://spatialreference.org