遥感数据机器学习的准备工作:python将栅格数据提取至EXCEL
大部分我们处理的降水、气温等栅格数据的格式是nc形式,需要我们将他转换成栅格数据并导入至Arcgis中,进行下一步操作。转好之后的数据我们就可以拿来用了,下面我们通过python代码实现快速的栅格属性导出至EXCEL表中。之后我们根据自己的研究区进行裁剪【Spatial Analyst 工具-提取分析-按掩膜提取】下面我们将裁剪好的研究区进行栅格转点【转换工具-由栅格转出-栅格转点】
·
大部分我们处理的降水、气温等栅格数据的格式是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")
更多推荐
已为社区贡献3条内容
所有评论(0)