系列文章目录: ArcGIS自定义脚本编程
一、引言
栅格(.tif)文件是一种常见的数据存储格式,在空间分析中的过程中,我们常常需要将栅格文件中包含的数据提取出来,导出为类似于{(x1,y1,v1), (x2,y2,v2)…}。其中,x、y分别为栅格文件中某个像元中心处对应的横坐标和纵坐标,v为此像元对应的值。
针对这一过程,即提取栅格数据并添加x、y字段的过程,通常可以利用ArcGIS进行,即对于浮点型栅格,先用栅格转点工具转为点要素,再使用添加XY坐标工具,再打开点要素的属性表导出。
但是,在处理文件大的栅格文件时,ArcGIS的速度简直是慢的令人发指。比如说提取全中国1km的高程栅格数据,光是栅格转点这一步就用了8分钟左右。顺带一提,用QGIS执行相同的操作,栅格转点只需10s,这种方法推荐给对python不熟悉的同学,具体操作如下:
本文利用gdal库实现了这一过程,提取全中国1km的高程栅格数据的所需时间与ArcGIS、QGIS对比如下:gdal(约2s)<< QGIS(约30s)<< ArcGIS(约15min)。
二、脚本代码
import time
from osgeo import gdal
import numpy as np
import pandas as pd
import os
def rasterToPoints(rasterfile, nodata=None, v_name=None):
"""
:param rasterfile: 待执行栅格转点的栅格文件
:param nodata:栅格中的无数据值,默认取栅格的最小值
:param v_name:导出表格中栅格值所在列的名称,默认为栅格的文件名
:return:x、y、value
"""
# numpy禁用科学计数法,pandas中存储浮点型时只保留四位小数
np.set_printoptions(suppress=True)
pd.set_option('display.float_format', lambda x: '%.4f' % x)
rds = gdal.Open(rasterfile) # type:gdal.Dataset
if rds.RasterCount != 1:
print("Warning, RasterCount > 1")
cols = rds.RasterXSize
rows = rds.RasterYSize
band = rds.GetRasterBand(1) # type:gdal.Band
transform = rds.GetGeoTransform()
print(transform)
x_origin = transform[0]
y_origin = transform[3]
pixel_width = transform[1]
pixel_height = transform[5]
if (pixel_height + pixel_width) != 0:
print("Warning, pixelWidth != pixelHeight")
# 读取栅格
values = np.array(band.ReadAsArray())
x = np.arange(x_origin + pixel_width * 0.5, x_origin + (cols + 0.5) * pixel_width, pixel_width)
y = np.arange(y_origin + pixel_height*0.5, y_origin + (rows+0.5) * pixel_height, pixel_height)
px, py = np.meshgrid(x, y)
if v_name is None:
v_name = os.path.splitext(os.path.split(rasterfile)[1])[0]
dataset = {"x": px.ravel(),
"y": py.ravel(),
v_name: values.ravel()}
df_temp = pd.DataFrame(dataset, dtype="float32")
# 删除缺失值
if nodata is None:
nodata = df_temp[v_name].min()
df_temp = df_temp[df_temp[v_name] != nodata]
else:
df_temp = df_temp[df_temp[v_name] != nodata]
df_temp.index = range(len(df_temp))
return df_temp
if __name__ == "__main__":
# 禁用科学计数法
np.set_printoptions(suppress=True)
pd.set_option('display.float_format', lambda x: '%.4f' % x)
# 执行栅格转点,并计时
s = time.time()
in_tif = r"D:\ChinaGW\rawData\surface_variables\dem.tif"
outfile = rasterToPoints(in_tif)
outfile.to_pickle("timeused.pkl")
e = time.time()
print("time used {0}s".format(e-s))
三、运行结果
输入:
待提取的栅格文件dem.tif
输出
pkl文件,包含x、y字段以及栅格的值
四、讨论
(一)使用例
例一:如何保存为csv文件
if __name__ == "__main__":
np.set_printoptions(suppress=True)
pd.set_option('display.float_format', lambda x: '%.4f' % x)
in_tif = r"D:\ChinaGW\rawData\surface_variables\dem.tif"
outfile = rasterToPoints(in_tif)
outfile.to_csv("timeused.csv") # 修改这一行
例二:修改栅格值所在列的列名
if __name__ == "__main__":
np.set_printoptions(suppress=True)
pd.set_option('display.float_format', lambda x: '%.4f' % x)
in_tif = r"D:\ChinaGW\rawData\surface_variables\dem.tif"
outfile = rasterToPoints(in_tif,v_name="newName")
outfile.to_pickle("timeused.pkl")
例三:更改nodata值
这里nodata值,可以通过ArcGIS查看栅格的属性查看。之所以本人函数中nodata值选择为像元值的最小值,是因为本人习惯把栅格的NoData值设置为一个小到离谱的值,方便识别。所以,在使用时请自行修改nodata的值。
if __name__ == "__main__":
np.set_printoptions(suppress=True)
pd.set_option('display.float_format', lambda x: '%.4f' % x)
in_tif = r"D:\ChinaGW\rawData\surface_variables\dem.tif"
outfile = rasterToPoints(in_tif, v_name="newName", nodata=-999)
outfile.to_pickle("timeused.pkl")
(二)不足
不足一:没有进度条显示
不足二:泛用性不足
仅适用于:栅格波段数为1;仿射地理变换参数GeoTransform[2]、和GeoTransform[4]这两个参数的值为0。