遥感影像重采样

一、简介

　　图像重采样就是从高分辨率遥感影像中提取出低分辨率影像，或者从低分辨率影像中提取高分辨率影像的过程。常用的方法有最邻近内插法、双线性内插法、三次卷积法等

二、重采样方法

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