MODIS数据进行重投影

由于MODIS数据采用的是SIN正弦投影 ，我们平常一般都是采用地理坐标，一般我们都会对MODIS数据进行重投影。MODIS Reprojection Tools（MRT）是专门用来对MODIS数据进行处理的，但是总感觉这软件操作起来麻烦。 所以今天我们就介绍一下两种基于Python中的GDAL对MODIS进行重投影的方法。

gdal.Warp

gdal.Warp是一个很好用的函数们可以用来重投影、影像裁剪等。用它对MODIS数据进行重投影很简单。

from osgeo import gdal
import numpy as np
from osgeo import osr
#使用gdal.Warp对MODIS数据进行重投影。
modis_lai=gdal.Open(r'D:\DATA\MODIS\LAI\MCD15A3H.A2020181.h13v11.006.2020188204620.hdf')
subdataset_one = modis_lai.GetSubDatasets()[0][0] # 第一个子数据集合
gdal.Warp(r'D:\code\reprojection.tif', subdataset_one, dstSRS='EPSG:4326')

gdal.ReprojectImage

gdal.ReprojectImage需要我们自己计算转换参数。

from osgeo import gdal
import numpy as np
from osgeo import osr
#使用gdal.Warp对MODIS数据进行重投影。
modis_lai=gdal.Open(r'D:\DATA\MODIS\LAI\MCD15A3H.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