假装有目录
所用软件:
GMS
QGIS
实现功能:
将GMS中计算的网格结果如 Head
ET
Recharge
等转化成GIS软件中可操作的 .SHP 格式文件,然后按照 GIS 的土地利用类型(栅格数据)进行分类汇总,生成每个土地利用类型的 mean
sum
max
等统计值
首先将3D网格转化成2D网格,此时2D网格只有默认的数据
然后在 Grid 中将 3D 数据转化成 2D 的
此时会出现相应选项,可以导出最大值、平均值、最高活跃网格值等
此时的 2D Grid 数据是无法导出成 GIS 可用的 .SHP 格式,因此继续将 2D Grid 数据转化成 2D Scatter 数据,2D Scatter 数据导出成 .SHP 文件后拖入到 QGIS 中,采用Sample raster value功能将土地利用类型数据提取到散点中,这样提取的好处在于土地利用类型的代码是直接被提取,不受插值的影响
提取完成以后得到新的散点数据,里面包含了一个名为 r_value1 新数据列,采用QGIS的 Statistics by categories 工具进行自动分类汇总,注意分类数据和目标数据的选择
最后生成的汇总表格如图,可导出成CSV用EXCEL进行后续操作
出现问题:
数据表中包含了空值-999,直接进行 Statistics by categories 操作会把空值计算到各个类型中,而由于数据表太大,在QGIS的属性表中删除这些空值会让程序运行崩溃。需要利用 GDAL/OGR 进行直接读取属性表操作。
安装 gdal
,在cmd中输入命令
conda install gdal
此处参考python shp文件操作
导出所需字段和要素的数据表文件,生成CSV
import pandas as pd
#导入库
try:
from osgeo import ogr
except:
import ogr
#加载相应数据类型的驱动,相当于初始化一个对象
driver = ogr.GetDriverByName('ESRI Shapefile')
def read_dbf_to_csv(name):
#打开数据
fileName="samplepoint%s.shp"%name
dataSource = driver.Open(fileName,0) #0是只读,1可写
if dataSource is None:
print('could not open')
sys.exit(1)
#获取图层
layer = dataSource.GetLayer(0)
#查看图层的信息
print("图层描述 :{0}".format(layer.GetDescription()))
print("图层范围 :{0}".format(layer.GetExtent()))
print("要素数量 :{0}".format(layer.GetFeatureCount()))
print("元数据描述 :{0}".format(layer.GetMetadata()))
print("空间参考 :{0}".format(layer.GetSpatialRef()))
#获取图层的字段信息
layerDefinition = layer.GetLayerDefn()
print("字段名 \t字段类型 \t字段长度 \t字段精度 \t字段类型代码")
for i in range(layerDefinition.GetFieldCount()):
fieldName = layerDefinition.GetFieldDefn(i).GetName()
fieldTypeCode = layerDefinition.GetFieldDefn(i).GetType()
fieldType = layerDefinition.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
fieldWidth = layerDefinition.GetFieldDefn(i).GetWidth()
GetPrecision = layerDefinition.GetFieldDefn(i).GetPrecision()
print("{0} \t{1} \t{2} \t{3} \t{4}".format(fieldName,fieldType,fieldWidth,GetPrecision,fieldTypeCode,))
#获取要素,并获取要素相应的字段
f_count=layer.GetFeatureCount()
dbf=pd.DataFrame(columns=['LUCC','ET','RCH'])
#获取属性表
for i in range(f_count):
feat = layer.GetFeature(i) # 获取要素
RCH=feat.GetField(4)
if not (RCH==-999):
LUCC=feat.GetField(5) # 获取要素的字段
ET=feat.GetField(3)
#geom = feat.GetGeometryRef() # 获取要素坐标
temp=pd.DataFrame([[LUCC, ET, RCH]],columns=['LUCC','ET','RCH'])
dbf=dbf.append(temp,ignore_index=True)
print('{:.4}%'.format(i/f_count*100)) #打印属性表获取进度
dbf.to_csv('%s.csv'%name)
dataSource.Destroy() #关闭数据源
for name in ['2000','2010','2015']:
read_dbf_to_csv(name)
导出成为 CSV 文件后进行分类汇总
此处参考资料分类汇总功能实现Group
import pandas as pd
df=pd.read_table('2015.csv',sep=',',index_col=0)
categories=df.groupby(by='LUCC')
print(categories.count())
print(categories.sum())
print(categories.mean())
最后就得到根据土地利用类型分类汇总的各项数据了