遥感影像图像处理：
重采样：
图像重采样就是从高分辨率遥感影像中提取出低分辨率影像，或者从低分辨率影像中提取高分辨率影像的过程。常用的方法有最邻近内插法、双线性内插法、三次卷积法等
二、重采样方法
1 使用ReadAsArray函数
def ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_obj=None,
buf_xsize = None, buf_ysize = None, buf_type = None,
resample_alg = GRIORA_NearestNeighbour,
callback = None,
callback_data = None)•xoff=0, yoff=0，指定从原图像波段数据中的哪个位置开始读取。
•win_xsize=None, win_ysize=None，指定从原图像波段中读取的行数和列数。
•buf_xsize=None, buf_ysize=None，指定暂存在内存中的新图像的行数和列数。
•buf_type=None，指定新图像的像素值的类型。
•buf_obj=None，指定新图像像素值数组的变量，因为整个方法也会返回一个新图像像素值的数组，用这两种方式获取重采样后的数组都可以。
•resample_alg=GRIORA_NearestNeighbour，重采样方法，默认为最近邻方法。
•callback=None,callback_data=None，回调函数和数据。
该函数的作用在于将一部分数据读取到已定义的一个数组中。从其参数 resample_alg来看，该函数可以完成重采样功能。但是需要对重采样后的地理变换进行重新设置。地理变换中包含像素大小等信息，重采样后，像素大小发生变化，地理变换也要随之更新
低分辨率重采样成高分辨率
# _*_ coding: utf-8 _*_
import os
from osgeo import gdal
os.chdir(r'D:\osgeopy-data\Landsat\Washington')
in_ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
in_band = in_ds.GetRasterBand(1)
out_rows = in_band.YSize * 2
out_columns = in_band.XSize * 2
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('band1_resampled.tif',
out_columns, out_rows)
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] /= 2
geotransform[5] /= 2
out_ds.SetGeoTransform(geotransform)
data = in_band.ReadAsArray(
buf_xsize=out_columns, buf_ysize=out_rows)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32, 64])
del out_ds高分辨率重采样成低分辨率
# _*_ coding: utf-8 _*_
import os
import numpy as np
from osgeo import gdal
os.chdir(r'D:\osgeopy-data\Landsat\Washington')
in_ds = gdal.Open('nat_color.tif')
out_rows = int(in_ds.RasterYSize / 2)
out_columns = int(in_ds.RasterXSize / 2)
num_bands = in_ds.RasterCount
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('nat_color_resampled.tif',
out_columns, out_rows, num_bands)
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] *= 2
geotransform[5] *= 2
out_ds.SetGeoTransform(geotransform)
data = in_ds.ReadRaster(
buf_xsize=out_columns, buf_ysize=out_rows)
out_ds.WriteRaster(0, 0, out_columns, out_rows, data)
out_ds.FlushCache()
for i in range(num_bands):
out_ds.GetRasterBand(i + 1).ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16])
del out_ds注意，在这种情况下，要确保行数和列数是整数，因为除法的结果可能是浮点数，如果不是整型数据，程序很可能报错。
2 使用warp函数
Gdal的Warp函数，该函数的作用是“图像重投影和变形”，函数中也有一个resampleAlg参数，可以用来指定重采样的方法，如果不指定resampleAlg，则默认使用最近邻方法，
#重采样方法为双线性重采样
gdal.Warp("resampletif.tif",dataset,width=newCols, height=newRows, resampleAlg=gdalconst.GRIORA_Bilinear)参数详解（未列完）
srcSRS 源坐标系统
dstSRS 目标坐标系统
resampleAllg 重采样方法
multeThread 多线程
cutLineDSname 裁剪mask矢量数据集名字
format 输出格式 eg GTIFF
cutLineLayername 裁剪mask图层名
cutLinewhere 裁剪where语句
例如下面的代码实现了使用warp函数进行重采样的功能(常用作处理时序影像)。
def ReprojectImages2():
# 若采用gdal.Warp()方法进行重采样
# 获取输出影像信息
inputrasfile = gdal.Open(inputfilePath, gdal.GA_ReadOnly)
inputProj = inputrasfile.GetProjection()
# 获取参考影像信息
referencefile = gdal.Open(referencefilefilePath, gdal.GA_ReadOnly)
referencefileProj = referencefile.GetProjection()
referencefileTrans = referencefile.GetGeoTransform()
bandreferencefile = referencefile.GetRasterBand(1)
x = referencefile.RasterXSize
y = referencefile.RasterYSize
nbands = referencefile.RasterCount
# 创建重采样输出文件（设置投影及六参数）
driver = gdal.GetDriverByName('GTiff')
output = driver.Create(outputfilePath, x, y, nbands, bandreferencefile.DataType)
output.SetGeoTransform(referencefileTrans)
output.SetProjection(referencefileProj)
options = gdal.WarpOptions(srcSRS=inputProj, dstSRS=referencefileProj, resampleAlg=gdalconst.GRA_Bilinear)
gdal.Warp(output, inputfilePath, options=options)3 使用gdal.ReprojectImage()进行重采样
参数说明（未列完）：
Dataset src_ds 输入数据集
Dataset dst_ds 输出文件
GDALResampleAlg eResampleAlg 重采样方法(最邻近内插\双线性内插\三次卷积等)
GDALProgressFunc 回调函数
char const * src_wkt=None 输入投影（原始投影）
char const * dst_wkt=None 参考投影（目标投影）
代码实现：
outputfilePath = 'G:/studyprojects/gdal/GdalStudy/Files/images/ReprojectImage.tif'
inputfilePath='G:/studyprojects/gdal/GdalStudy/Files/images/2016CHA.tif'
referencefilefilePath='G:/studyprojects/gdal/GdalStudy/Files/images/2018CHA.tif'
def ReprojectImages():
# 获取输出影像信息
inputrasfile = gdal.Open(inputfilePath, gdal.GA_ReadOnly)
inputProj = inputrasfile.GetProjection()
# 获取参考影像信息
referencefile = gdal.Open(referencefilefilePath, gdal.GA_ReadOnly)
referencefileProj = referencefile.GetProjection()
referencefileTrans = referencefile.GetGeoTransform()
bandreferencefile = referencefile.GetRasterBand(1)
Width= referencefile.RasterXSize
Height = referencefile.RasterYSize
nbands = referencefile.RasterCount
# 创建重采样输出文件（设置投影及六参数）
driver = gdal.GetDriverByName('GTiff')
output = driver.Create(outputfilePath, Width,Height, nbands, bandreferencefile.DataType)
output.SetGeoTransform(referencefileTrans)
output.SetProjection(referencefileProj)
# 参数说明 输入数据集、输出文件、输入投影、参考投影、重采样方法(最邻近内插\双线性内插\三次卷积等)、回调函数
gdal.ReprojectImage(inputrasfile, output, inputProj, referencefileProj, gdalconst.GRA_Bilinear,0.0,0.0,)
