MODIS数据进行重投影
由于MODIS数据采用的是SIN正弦投影 ，我们平常一般都是采用地理坐标，一般我们都会对MODIS数据进行重投影。MODIS Reprojection Tools（MRT）是专门用来对MODIS数据进行处理的，但是总感觉这软件操作起来麻烦。 所以今天我们就介绍一下两种基于Python中的GDAL对MODIS进行重投影的方法。
gdal.Warp
gdal.Warp是一个很好用的函数们可以用来重投影、影像裁剪等。用它对MODIS数据进行重投影很简单。
代码语言：javascript复制from osgeo import gdal
import numpy as np
from osgeo import osr
#使用gdal.Warp对MODIS数据进行重投影。
modis_lai=gdal.Open(r'D:DATAMODISLAIMCD15A3H.A2020181.h13v11.006.2020188204620.hdf')
subdataset_one = modis_lai.GetSubDatasets()[0][0] # 第一个子数据集合
gdal.Warp(r'D:codereprojection.tif', subdataset_one, dstSRS='EPSG:4326')
gdal.ReprojectImage
gdal.ReprojectImage需要我们自己计算转换参数。
代码语言：javascript复制from osgeo import gdal
import numpy as np
from osgeo import osr
#使用gdal.Warp对MODIS数据进行重投影。
modis_lai=gdal.Open(r'D:DATAMODISLAIMCD15A3H.A2020181.h13v11.006.2020188204620.hdf')
subdataset_one = modis_lai.GetSubDatasets()[0][0] # 第一个子数据集合
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
# 获得输出数据的投影，建立两个投影直接的转换关系
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(epsg_to)
tx = osr.CoordinateTransformation(src_srs, dst_srs)
# 计算输出图像四个角点的坐标
(uly, ulx, _) = tx.TransformPoint(srs_trans[0], srs_trans[3])
(ury, urx, _) = tx.TransformPoint(srs_trans[0] srs_trans[1] * x_size, srs_trans[3])
(lly, llx, _) = tx.TransformPoint(srs_trans[0], srs_trans[3] srs_trans[5] * y_size)
(lry, lrx, _) = 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, gdal.GDT_Float32)
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
reproject(subdataset_one,
r'D:codereprojection2.tif', 0.005, 0.005, 4326)
这个reproject函数中，我们可以自己设置数据数值类型、空间分辨率，投影。上述两种方法都可以实现MODIS投影的功能，感觉会比MRT更加灵活一些吧。
