gdal重采样深度讲解（前后不改变数据总量算法）（比arcgis好用）

从本科学习arcgis开始，就一直没有离开过重采样，arcgis里重采样方法大概有一下4种：

• NEAREST— The fastest resampling method; it minimizes changes to pixel values. Suitable for discrete data, such as land cover.
• BILINEAR— Calculates the value of each pixel by averaging (weighted for distance) the values of the surrounding 4 pixels. Suitable for continuous data.
• CUBIC— Calculates the value of each pixel by fitting a smooth curve based on the surrounding 16 pixels. Produces the smoothest image but can create values outside of the range found in the source data. Suitable for continuous data.
• MAJORITY— Determines the value of each pixel based on the most popular value within a 3 by 3 window. Suitable for discrete data

但是可以很容易发现这4种方法其实都对改变了原始数据计算的结果：比如本来有这样一个4*4网格，如果网格的值代表了NPP（g/m2），每个网格1m2，这块土地的NPP总量为15+10=25g，当采用Nearest的方法重采样到5*5的网格，得到的结果如下：

此时NPP的总量为24*0.64+10*0.64=21.76g，和之前有3.24g的区别，这个误差还是值得考虑的。所以真正的无误差算法应该是计算新网格内旧网格值对应面积的加权值，arcgis里的4种方法都不是这样的，但是最近发现gdal函数里有一个Average方法，这个函数可以完全实现总量不变：

计算后的NPP总量：0.64*（10+3.25*2+1.5625+21）=25，每个网格的值相当于把5*5的网格放在4*4的网格上，再计算出的平均值，十分适合这种需要总量不变的情况。

当然啦，重采样的方法还是要看你处理数据的需要了，如果处理土地利用数据，当然要用Nearest方法啦。另外还有Bilinear、Cubic等各种重采样方法，我这里就不详细展开了。

最后附上我的测试代码（只要更改路径就可以直接运行）：

import numpy as np
import gdal
import gdalconst
#---------------------更改输入路径----------------------------------
pathin = ''
#------------------------------------------------------------------
proj_ori = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
gt_ori = (0.0, 0.25, 0.0, 0.0, 0.0, -0.25)
#-----------------创建实验mask 4*4-----------------------------------------
out_ds = gdal.GetDriverByName('GTiff').Create(pathin+'test.tif' ,4,4,1,gdalconst.GDT_Float32,
                                                 ['COMPRESS=LZW'])#path, im_width, im_height, im_bands, datatype
out_ds.SetProjection(proj_ori)
GetGeoTransform = [0,0.25,0,0,0,-0.25]
out_ds.SetGeoTransform(GetGeoTransform)
ra = np.ones((4,4))
ra[0,0]=10#给位于（1,1）位置的像元赋值10，方便后面解释算法
out_band = out_ds.GetRasterBand(1).WriteArray(ra)
del out_ds
#-----------------读取实验mask的参数以传递给重采样后的mask---------------------
#Nearest
ras_ori  = gdal.Open(pathin + 'test.tif')
proj_ori = ras_ori.GetProjection()
gt_ori   = ras_ori.GetGeoTransform()
gt_dst = list(gt_ori)
gt_dst[1] = 0.2#更改分辨率
gt_dst[5] = -0.2
gt_dst = tuple(gt_dst)                        #得到geotransform_dst
ras_dst = gdal.GetDriverByName('GTiff').Create(pathin+'Nearest.tif',5,5,
                                               1,gdalconst.GDT_Float32,
                                               ['COMPRESS=LZW'])
ras_dst.SetGeoTransform(gt_dst)
ras_dst.SetProjection(proj_ori)
gdal.ReprojectImage(ras_ori,ras_dst,proj_ori,proj_ori,gdalconst.GRA_NearestNeighbour)
del ras_ori
del ras_dst



#Average
ras_ori  = gdal.Open(pathin + 'test.tif')
proj_ori = ras_ori.GetProjection()
gt_ori   = ras_ori.GetGeoTransform()
gt_dst = list(gt_ori)
gt_dst[1] = 0.2#更改分辨率
gt_dst[5] = -0.2
gt_dst = tuple(gt_dst)                        #得到geotransform_dst
ras_dst = gdal.GetDriverByName('GTiff').Create(pathin+'average.tif',5,5,
                                               1,gdalconst.GDT_Float32,
                                               ['COMPRESS=LZW'])
ras_dst.SetGeoTransform(gt_dst)
ras_dst.SetProjection(proj_ori)
gdal.ReprojectImage(ras_ori,ras_dst,proj_ori,proj_ori,gdalconst.GRA_Average)
del ras_ori
del ras_dst