python使用gdal实现影像重采样

gdal实现影像的重采样

对于tif格式的遥感影像而言，通常使用左上角那个像元真实的x坐标值，像素宽度，x像素旋转角度，原点y坐标，y像素旋转角度，像素高度（负值）共六个参数来表示整幅影像像素点的坐标。

当影像分辨率改变后，x坐标值和y坐标值度都不会改变，但像素高度和宽度会改变，而像素高度和宽度改变的话，影像的行列数也会改变。

1.gdal实现影像重采样至高分辨率

这里，将一幅影像的像元高度和宽度都变为原来的0.5倍，实现代码如下：

'''
Created on 2020年2月10日

@author: Sun Strong
'''
from osgeo import gdal
import os
path=r"C:\Users\Sun Strong\Desktop\gdal" 
os.chdir(path)
in_ds=gdal.Open("p047r027_7t20000730_z10_nn10.tif")
in_band=in_ds.GetRasterBand(1)
xsize=in_band.XSize
ysize=in_band.YSize
geotrans=list(in_ds.GetGeoTransform())
geotrans[1]/=2#像元宽度变为原来的0.5倍
geotrans[5]/=2#像元高度变为原来的0.5倍
#重采样后的影像
if os.path.exists('resampled.tif'):#如果已存在同名影像
    os.remove('resampled.tif')#则删除之
out_ds=in_ds.GetDriver().Create('resampled.tif',xsize*2,ysize*2,1,in_band.DataType)#创建一幅重采样后的影像的句柄，行列数都变成原来的2倍
out_ds.SetProjection(in_ds.GetProjection())#设置投影坐标
out_ds.SetGeoTransform(geotrans)#设置地理变换参数
data=in_band.ReadAsArray(buf_xsize=xsize*2, buf_ysize=ysize*2)#使用更大的缓冲读取影像，与重采样后影像行列对应
out_band=out_ds.GetRasterBand(1)
out_band.WriteArray(data)#写入数据到新影像中
out_band.FlushCache()
out_band.ComputeBandStats(False)#计算统计信息
out_ds.BuildOverviews('average',[1,2,4,8,16,32])#构建金字塔
del out_ds#删除句柄
del in_ds
print("This process has succeeded!")

2.gdal实现重采样至更低分辨率

这里，将一幅影像的像元高度和宽度都变为原来的2倍，实现代码如下：

'''
Created on 2020年2月10日

@author: Sun Strong
'''
from osgeo import gdal
import os
import numpy as np
path=r"C:\Users\Sun Strong\Desktop\gdal"
os.chdir(path)
image_name='p047r027_7t20000730_z10_nn30.tif'
in_ds=gdal.Open(image_name)
geotrans=list(in_ds.GetGeoTransform())
geotrans[1]*=2#像元宽度变为原来的两倍
geotrans[5]*=2#像元高度也变为原来的两倍
in_band=in_ds.GetRasterBand(1)
xsize=in_band.XSize
ysize=in_band.YSize
x_resolution=int(xsize/2)#影像的行列都变为原来的一半
y_resolution=int(ysize/2)
if os.path.exists('LowerResolutionImage.tif'):#如果存在重采样后的影像，则删除之
    os.remove('LowerResolutionImage.tif')
out_ds=in_ds.GetDriver().Create('LowerResolutionImage.tif',x_resolution,y_resolution,1,in_band.DataType)#创建一个构建重采样影像的句柄
out_ds.SetProjection(in_ds.GetProjection())#设置投影信息
out_ds.SetGeoTransform(geotrans)#设置地理变换信息
data=np.empty((y_resolution,x_resolution),np.int)#设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值
in_band.ReadAsArray(buf_obj=data)
out_band=out_ds.GetRasterBand(1)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False)#计算统计信息
out_ds.BuildOverviews('average',[1,2,4,8,16,32])#构建金字塔
del in_ds#删除句柄
del out_ds
print("This process has succeeded!")

3.总结

影像的重采样，需要注意重采样前后六个地理变换参数的变化，还有重采样后影像像元的宽度、高度以及行列数目。上述例子提供了两种读取像元值的方法（对应重采样影像的行列数目）：
一、指定重采样后影像的行列数目来读取像元值
data=in_band.ReadAsArray(buf_xsize=xsize2, buf_ysize=ysize2)
二、准备一个numpy数组来接收像元值
data=np.empty((y_resolution,x_resolution),np.int)
in_band.ReadAsArray(buf_obj=data)