【python+GIS】python实现多个栅格值提取到点，写入EXCEL表

一、应用背景

下载了FLUXNET的站点，也下载了地表温度数据，想要将所有的温度数据（tif）的值提取到点上。

二、参考文献

《python快速实现shp点提取栅格值 (两种方法对比)》
《Python+GDAL | 基于矢量（点）读取栅格数据》

三、代码

from osgeo import ogr
import os, sys
from osgeo import gdal
from osgeo.gdalconst import *
import csv
import xlwt 
inputSHP=r'D:\DATA\pointToraster\points.shp'  #点数据文件
InputRasterFolder=r'D:\DATA\pointToraster\tif'  #放栅格数据的文件夹
# 设置Excel编码
file = xlwt.Workbook('encoding = utf-8') 
# 创建sheet工作表
sheet1 = file.add_sheet('sheet1',cell_overwrite_ok=True) 
#改变工作空间
#############获取矢量点位的经纬度
#设置driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#打开矢量
ds = driver.Open(inputSHP, 0)
#获取图层
layer = ds.GetLayer()

#获取要素及要素地理位置
xValues = []
yValues = []
feature = layer.GetNextFeature()
while feature:
    geometry = feature.GetGeometryRef()
    x = geometry.GetX()
    y = geometry.GetY()
    xValues.append(x)
    yValues.append(y)  
    feature = layer.GetNextFeature()

#############获取点位所在像元的栅格值
#读取栅格
#获取注册类
#打开栅格数据

input_folder_list=os.listdir(InputRasterFolder) #读取文件夹里所有文件
tif_files=list() #创建一个只装tif格式的列表
for filename in input_folder_list:  #遍历
    if os.path.splitext(filename)[1] == '.tif':  #不管文件名里面有多少个tif，都只认最后一个tif
        tif_files.append(filename)  #将文件夹里的tif文件加入只有tif的列表
print(tif_files)
sheet1.write(0, 0, "Lon") #excel表的第1列为经度
sheet1.write(0, 1, "Lat") #excel表的第2列为纬度

for i in range(0,len(tif_files)): #遍历tif
    sheet1.write(0, i+2, filename) #在表格第一行设置列名
    ds = gdal.Open(InputRasterFolder+'\\'+tif_files[i], GA_ReadOnly)
    #获取行列、波段
    rows = ds.RasterYSize
    cols = ds.RasterXSize
    bands = ds.RasterCount
    #获取放射变换信息
    transform = ds.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]
    #
    values = []
    for j in range(len(xValues)):       #遍历所有点
        x = xValues[j]
        y = yValues[j]
        #获取点位所在栅格的位置
        xOffset = int((x-xOrigin)/pixelWidth)
        yOffset = int((y-yOrigin)/pixelHeight)    
        band = ds.GetRasterBand(1) #读取影像
        data = band.ReadAsArray(xOffset, yOffset,1,1) #读出从(xoffset,yoffset)开始，大小为(xsize,ysize)的矩阵
        value = data[0,0]*0.01 #乘以参数，这个根据自己的数据情况做出修改
        #将数据经纬度和对应栅格数值写入excel表
        sheet1.write(j + 1, 0, x)    #第j+1行，第1列 
        sheet1.write(j + 1, 1, y)    #第j+1行，第2列 
        sheet1.write(j + 1, i+2, value) 


file.save(r'D:\DATA\pointToraster\pointToraster.xls') #保存表格

print("OK")

四、结果

最终的结果呈现为一个EXCEL表：

为了验证是否正确，将它与ARCGIS执行【多值提取到点】后的结果进行比较，验证通过：