遙感影像重采樣

一、簡介

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

二、重采樣方法

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