GDAL：利用一幅遥感影像裁剪另一幅遥感影像

1. 背景

最近在利用GEE下载多源遥感影像（一幅为Sentinel，10m分辨率；另一幅为GOES静止卫星影像，2000米分辨率）时，下载的中心点坐标与下载半径设置的完全一致，但下载的数据却存在几个像素的偏差，这阻碍了自己下一步工作的开展。最终，我决定利用GDAL写一个算法，将范围较大的遥感影像，利用另一幅较小的遥感影像进行裁剪，进而使两者的空间范围重合。代码如有不恰当的地方，欢迎大家讨论交流。说明：本文暂时只适用于“一幅遥感影像完全在另一幅遥感影像范围内的情况”

2. 代码实现

from osgeo import gdal
import datetime
'''
Description:Using one image to clip another image (only for one image is completely within the range of the other)
Author: Dawn
Date: 2023/03/22
'''
def get_time():
    return datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')

# KeyM: get spatial range of the Image
def getSpatialRage(image):
    geo_transform = image.GetGeoTransform()
    pixel_width = geo_transform[1]
    pixel_height = geo_transform[5]
    width = image.RasterXSize; height = image.RasterYSize
    min_x = geo_transform[0]
    max_y = geo_transform[3]
    max_x = min_x + width * pixel_width
    min_y = max_y + height * pixel_height
    return min_x, max_y, max_x, min_y

# keyM: set image baseline
def imageBaseline(ds1,ds2):
    # get the spatial region of image_1
    min_x1, max_y1, max_x1, min_y1 = getSpatialRage(ds1)
    # get the shape of image_2
    min_x2, max_y2, max_x2, min_y2 = getSpatialRage(ds2)
    # define clipImage and baseImage
    if (max_x1 - min_x1) > (max_x2 - min_x2):
        baseImage = ds1
        clipImage = ds2
    else:
        baseImage = ds2
        clipImage = ds1
    return baseImage, clipImage

# KeyM： clip method
def clipMethod(baseImage, clipImage):
    '''
    :param baseImage: large range image
    :param clipImage: small range image
    '''
    #Steps: transform spatial coordinate to pixle coordinate
    geo_transform = baseImage.GetGeoTransform()
    geo_proj = baseImage.GetProjection()
    min_x, max_y, max_x, min_y = getSpatialRage(clipImage) # spatial coordinates
    pixel_minx = int((min_x - geo_transform[0]) / geo_transform[1])  # col
    pixel_maxy = int((min_y - geo_transform[3]) / geo_transform[5])  # row
    pixel_maxx = int((max_x - geo_transform[0]) / geo_transform[1])
    pixel_miny = int((max_y - geo_transform[3]) / geo_transform[5])
    new_ClipImage_Height = pixel_maxy - pixel_miny
    new_ClipImage_Weight = pixel_maxx - pixel_minx

    # create output data (clip_data)
    output_data = gdal.GetDriverByName('GTiff').Create('../filePath/new_ClipImage.tif',new_ClipImage_Weight,new_ClipImage_Height,1,gdal.GDT_Int16)
    # define new geoTransform, six parameters are left-up x, horizontal resolution, rotation angle, left-up y, rotation angle, vertical resolution, respectively
    new_transform = (min_x,geo_transform[1],geo_transform[2],max_y,geo_transform[4],geo_transform[5])
    output_data.SetGeoTransform(new_transform)
    output_data.SetProjection(geo_proj)
    data_Inf = baseImage.ReadAsArray(pixel_minx, pixel_miny, new_ClipImage_Weight, new_ClipImage_Height)
    output_data.GetRasterBand(1).WriteArray(data_Inf)

# TODO:设置全局变量
# KeyV:open the first image
ds1 = gdal.Open('../filePath/2020LC.tif')
# KeyV:open the second image
ds2 = gdal.Open("../filePath/geoImage.tif")
if __name__ == '__main__':
    print('>= start clippling:', get_time())
    # 自动根据两幅影像的空间范围决定哪一幅为“被裁剪影像”，哪一幅为“裁剪影像”
    baseImage, clipImage = imageBaseline(ds1, ds2)
    clipMethod(baseImage, clipImage)
    print('>=finish Clipping:',get_time())

3. 裁剪前后对比

通过目视对比，以及查看图像的属性信息，发现裁剪结果完全符合预期。

4. 存在的问题

在查看输出影像的结果时，发现裁剪后导出的影像数据量增大了几十倍。但是通过查看图像的属性，发现无论是影像的长、宽，还是像素大小，和裁剪前的影像是完全一致的。关于这点并没有解决，欢迎大家在留言中指教，后续找到解决方案后，我也会将其放到留言中。