一、簡介

　　影象重取樣就是從高解析度遙感影像中提取出低解析度影像，或者從低解析度影像中提取高解析度影像的過程。常用的方法有最鄰近內插法、雙線性內插法、三次卷積法等

二、重取樣方法

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,)