Python地学分析 — GDAL对遥感影像重采样

欢迎关注博主的微信公众号：“智能遥感”。

该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。

本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

“ 对遥感影像重采样是遥感数据预处理常用的手段之一，本节重点讲解高、低分辨率图像重采样原理与方法。”

思来想去，如果一味地给大家直接上代码，会显得有些枯燥无聊，可读性会有所降低。为了提高大家阅读的积极性，从本节开始，给大家展现一些高清遥感图像实例（封面图）。这样既可以学习到遥感科学在实际生活中的用途，还可以体会遥感之美。怎么突然觉得自己好有情怀？？？下面请先看封面简介：

图像来自LandSat官网

厄勒海峡大桥：1991年，丹麦和瑞典政府同意建立一座连接厄勒海峡两岸国家间的桥梁。瑞典马尔默（右）和丹麦哥本哈根（左）之间16公里长的厄勒海峡干线于2000年建成并通车。随着冰河时代末期海平面上升，切断了两者之间的陆地连接，丹麦和瑞典7000年后再次建立起两国的联系。

厄勒海峡链接有三个主要部分。在丹麦方面，链接从一条3,510米的水下隧道开始。隧道从水下延伸到一个长4055米的人工岛Peberholm上的一条道路上，该岛在图像中的天然岛屿南部呈现出明亮的白色形状。电缆支撑的厄勒海峡大桥横跨海峡东部向瑞典延伸7,845米，在整个图像上划出一条细长的白线（文字描述译自LandSat官网）。

遥感影像重采样

针对不同的遥感业务场景或者科研需求，对应的不同分辨率的遥感影像有着不同的用处。什么叫遥感影像的重采样呢？其定义是：根据各相邻的原采样点内插出新采样点的过程。内插的方法有双线性插值法、双三次卷积法和最邻近像元法等。重采样的定义是不是很专业？咳咳，好像有点抽象唉，没事没事，看看下面的示意图，你就能明白了。

上图过程表示从低分辨率图像往高分辨率图像进行重采样。在编写程序时，需要一个更大的数组来存储重采样后的新图像。左边图像将每个像元内插为对应的四个小像元，从6个大像元重采样成24个小像元。看到这里，相信大家就基本明白什么是重采样了吧。下面来看看具体的Python程序实战，其实重采样就这么简单。

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/10/23 15:09'import os
from osgeo import gdalos.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 * 2gtiff_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

这个例子有一些值得注意的重要事项。首先，在创建新数据集时，你将行数和列数加倍，并将这些相同的数字作为参数传递给ReadAsArray。这样就可以确保输入数据维度与输出数据维度相匹配，并且还可以将数据重新采样到较大的维度。你可以使用buf_obj参数的现有数组并获得相同的结果，而不是使用buf_xsize和buf_ysize参数。你还可以提供win_xsize和win_ysize参数，但默认为行和列的原始大小。

这说明在不改变地理转换大小的情况下，对比重新采样前后的结果。左上角的较小图像是正确的。较大的一个是通过使用输入图像中未编辑的地理转换创建的，与矢量图不匹配。

上面讲述了如何重低分辨率重采样成高分辨率，而重高分辨率重采样成低分辨率，则是实际遥感应用中比较多的情况。最近邻插值常用于将图像重采样到较小尺寸时，这种情况下，在输出中使用每四个像元块的右下像元值。

# _*_ coding: utf-8 _*_
# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/10/23 15:50'import osimport numpy as np
from osgeo import gdalos.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.RasterCountgtiff_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

该代码与第一个代码类似，不同之处在于输出行和列的数量减半而不是加倍，并且像元大小加倍而不是减半。请注意，在这种情况下，你要确保行数和列数是整数，因为除法的结果可能是浮点数，如果不是整型数据，程序很可能报错，可以自己尝试一下。

重采样前、后对比图，下图变得模糊，其实是分辨率变低了