栅格数据投影转换
使用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)+Xpixel∗GT(1)+Yline∗GT(2)Ygeo=GT(3)+Xpixel∗GT(4)+Yline∗GT(5)Xgeo=GT(0)+Xpixel∗GT(1)+Yline∗GT(2)Ygeo=GT(3)+Xpixel∗GT(4)+Yline∗GT(5)
\begin{matrix} X_{geo} = GT(0) + X_{pixel} * GT(1) + Y_{line} * GT(2) \\ Y_{geo} = GT(3) + X_{pixel} * GT(4) + Y_{line} * GT(5) \\ \end{matrix} （Xge0Xge0X_{ge0}, Yge0Yge0Y_{ge0}）表示对应于图上坐标（XpixelXpixelX_{pixel}, YlineYlineY_{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)
腾讯云开发者
