栅格数据投影转换

栅格数据投影转换

作者：阿振

邮箱：tanzhenyugis@163.com

博客：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反射变换的六参数模型：

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

Xgeo=GT(0)+Xpix