Python + GDAL处理数据(3):提取栅格数据到点和栅格重采样

  在进行Meta-Analysis时，经常需要整合文献报道的数据，但大多数时候我们是无法完全获取到这些信息的，比如在研究降水对生态系统生产力的影响时，可能就很少会报道土壤氢离子浓度指数(pH) 或者土壤容重(BD) 等信息，这时我们可能会需要从一些可信赖的数据源去获取这些信息。一般来说，最优的数据获取方式是直接联系论文的作者，但有时候这并不是最有效的方式；其次就是通过其他相同位点的研究报道来获取，但很多时候，相同位点进行研究报道的情况可能是一致的，也就是说，同样都是研究降雨的影响的话，可能都会忽略对土壤属性的报道；因而，我们需要提取报道的栅格数据来填补空缺数据。
  当然，ArcMap很容易就能完成这些工作，但由于ArcGis不是免费的，在一些情况下，可能是没法实现的，而Python开源的属性和强大的第三方包，学会Python实现这些功能很重要(哪怕你学会了后不用，至少可以避免一些利益纠纷)。

def read_img(filename): #读取影像
    dataset = gdal.Open(filename)  # 打开文件
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
    im_proj = dataset.GetProjection()  # 地图投影信息
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height).astype(np.float)  # 将数据写成数组，对应栅格矩阵
    return dataset,im_proj, im_geotrans, im_data, im_height, im_width #输出文件、坐标系、仿射矩阵、栅格数据、栅格的高和宽
def GetPointsData(ds,Colname): #获取点上的数据
    global df_values #设置df_values为全局变量，因而在每次使用时都需要确保df_values为所需要的数据
    #获取放射变换信息
    bands = ds.RasterCount #获取栅格的波段数
    transform = ds.GetGeoTransform() #获取栅格的仿射矩阵
    xOrigin = transform[0] #获取栅格左上角栅格的点位信息
    yOrigin = transform[3] #获取栅格左上角栅格的点位信息
    pixelWidth = transform[1] #获取栅格的宽度
    pixelHeight = transform[5] #获取栅格的高度
    for i in df_values.index: #对df_values的每个数据进行数据提取
        xOffset = int((df_values.loc[i,"Longitude"]-xOrigin)/pixelWidth) #获取点所在的位置
        yOffset = int((df_values.loc[i,"Latitude"]-yOrigin)/pixelHeight) #获取点所在的位置
        df_values.loc[i,"Tif_Xs"] = xOffset 
        df_values.loc[i,"Tif_Ys"] = yOffset
        for  k in range(bands): #分别求不同波段的数据
            band = ds.GetRasterBand(k+1) #获取k+1波段的数据，由于python是0开始，gdal是1开始，因而在python中使用gdal时需要加1
            data = band.ReadAsArray(xOffset, yOffset,1,1) #读取栅格位置的数据 #xsize,ysize,是栅格的列数和行数
            value = data[0,0] #data为二维的数据，需要提取出来
            if value != 1e+20: #CMIP6使用1e+20作为缺省值，因而对于非缺省值，应该保存他的数据
                df_values.loc[i,Colname] = value
            else:
                df_values.loc[i,Colname] = " " #如果是缺省值，直接填空格来占位
    del ds #删除内存占用栅格
    return df_values #返回df_values
if __name__ == "__main__":
	Ta_path = "./Ta_2015.tif"
	ds,proj, geotrans, values, row, column = read_img(Ta_path) #读取温度栅格，获取数据array以及栅格信息
	GetPointsData(ds,"Ta"+RCP+str(YearsNum[Ind])) #获取df_values中每个点每年的温度数据

  这样，我们只需要有需要数据的有效栅格就可以获取到我们需要的点的相关数据。但很多时候，我们下载到的数据精度是不一致，因而，将数据进行重采样就非常重要了。因而，这里也对数据重采样进行处理。

def ReprojectImages(referencePath,inputfilePath,outputfilePath): #影像重采样
	# 获取参考影像信息, 其实可以自定义这些信息，有参考的话就不用查这些参数了
	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
    # 获取输入影像信息
    inputrasfile = gdal.Open(inputfilePath, gdal.GA_ReadOnly) #打开输入影像
    inputProj = inputrasfile.GetProjection() #获取输入影像的坐标系
    # 创建重采样输出文件（设置投影及六参数）
    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,)
if __name__ == "__main__":
	inputfilePath = "Tiff Store/Ta_2015.tif"
	outputfilePath = './Tiff Store/Ta_2015_10KM.tif"
	referencePath = './Tiff Store/Ta_1990_10km_wgs1984.tif' # 重采样参考文件
	ReprojectImages(referencePath,inputfilePath,outputfilePath)

关于GDAL的函数参数，可以参考：
[ https://www.osgeo.cn/python_gdal_utah_tutorial/index.html ]