背景
最近在利用GEE下载多源遥感影像(一幅为Sentinel,10m分辨率;另一幅为GOES静止卫星影像,2000米分辨率)时,下载的中心点坐标与下载半径设置的完全一致,但下载的数据却存在几个像素的偏差,这阻碍了自己下一步工作的开展。最终,我决定利用GDAL写一个算法,将范围较大的遥感影像,利用另一幅较小的遥感影像进行裁剪,进而使两者的空间范围重合。代码如有不恰当的地方,欢迎大家讨论交流。说明:本文暂时只适用于“一幅遥感影像完全在另一幅遥感影像范围内的情况”
代码实现
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())
裁剪前后对比
通过目视对比,以及查看图像的属性信息,发现裁剪结果完全符合预期。
存在的问题
在查看输出影像的结果时,发现裁剪后导出的影像数据量增大了几十倍。但是通过查看图像的属性,发现无论是影像的长、宽,还是像素大小,和裁剪前的影像是完全一致的。关于这点并没有解决,欢迎大家在留言中指教,后续找到解决方案后,我也会将其放到留言中。