遥感数据机器学习的准备工作：python将栅格数据提取至EXCEL

大部分我们处理的降水、气温等栅格数据的格式是nc形式，需要我们将他转换成栅格数据并导入至Arcgis中，进行下一步操作。

之后我们根据自己的研究区进行裁剪【Spatial Analyst 工具-提取分析-按掩膜提取】

下面我们将裁剪好的研究区进行栅格转点【转换工具-由栅格转出-栅格转点】

转好之后的数据我们就可以拿来用了，下面我们通过python代码实现快速的栅格属性导出至EXCEL表中。

代码如下：

from osgeo import ogr
import os, sys
from osgeo import gdal
from osgeo.gdalconst import *
import csv
import xlwt

inputSHP = r'E:\CN\RasterT_tif1.shp'  #点数据文件
InputRasterFolder = r'E:\CN\train'  #放栅格数据的文件夹
# 设置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'E:\CN\RESULT\pointToraster.xls')  #保存表格

print("OK")