第十届全国大学生GIS应用技能大赛试题及参考解题过程(A下午)
作者:互联网
本次操作使用的ArcGIS版本为 10.8 。
第Ⅰ部分:试题
一、 案例背景
\(\quad\quad\)太阳能是一种可再生资源,是指太阳的热辐射能。太阳能资源丰富,既可免费使用,又无需运输,对环境无任何污染。太阳能的利用还不是很普及,太阳能的使用受到昼夜、季节、地理纬度和海拔高度等自然条件的限制以及晴、阴、云、雨等随机因素的影响。
\(\quad\quad\)某住宅小区希望在屋顶安装太阳能,供给家庭日常用电。你将为该小区评估是否适合安装太阳能电池板。
二、 数据说明(见 “上午A” 文件夹中的 “数据” 文件夹)
- Building.shp:建筑物数据。
- DSM.tif:数字表面模型。
- DTM.tif:数字地面模型。
三、 分析要求(100分)
-
根据建筑物修正DTM(15分)
DTM 一般指数字地面模型,在本项目中,你可以理解为 DEM。
DSM 是数字地面模型,包括地形及地表的所有对象 (例如树和建筑物等) 。
本项目中的 DTM 和 DSM 数据都是通过雷达数据生成,固有一定的误差。正确的 DTM 应该在建筑物处是平整的,而不是倾斜或者高低起伏。
1) 请为 Building 添加属性字段【基本高度】,并计算数值。 (5分)
\(\quad\quad\)基本高度指建筑物底部高程值,本题请计算每栋建筑物范围内的 DTM 的平均值,保留 2 位小数位数。2) 修正 DTM 数据,将结果命名为 “DTM 修正” 。 (10分)
\(\quad\quad\)建筑物所在位置对应的 DTM 应是平整的,高程值为建筑物的基本高度。 -
计算每栋房屋的其它基础信息(20分)
1) 请为 Building 添加属性字段【最大高度】,并计算数值。 (5分)
\(\quad\quad\)最大高度指建筑物屋顶的最大高程值,保留 2 位小数位数。
2) 请为 Building 添加属性字段【建筑物高度】,并计算数值。 (5分)
\(\quad\quad\)建筑物高度指建筑物本身的高度,保留 2 位小数位数。
3) 填写下表中的建筑物屋顶形态。 (10)
\(\quad\quad\)屋顶形态分为平屋顶、双坡屋顶、四坡屋顶。屋顶朝向分为平面、东、南、西、北、东南、西南、东北和西北。如果建筑物是双坡屋顶,请填写 2 个屋顶朝向;如果建筑物是四坡屋顶,请填写 4 个屋顶朝向。序号 BuildingBM 屋顶形态 屋顶朝向 1 2116 2 2156 3 2161 4 2165 5 2171 -
创建房屋屋顶区域 2021 年每月预计获得太阳辐射量栅格数据,在环境设置中,将 Building 作为掩膜(15分)
1) 【太阳辐射区域】工具可以计算从栅格表面获得的入射太阳辐射。输出总辐射栅格用于表示为输入表面的每个位置所计算的全局辐射或全部日照入射量(直射 + 散射)的输出栅格。输出单位为瓦特每小时每平方米(WH/m2)。
\(\quad\quad\)a)本题工具主参数中天空大小为 200,间隔小时数 0.5,其他主参数请自行研究。
\(\quad\quad\)b)本题工具地形参数中地形方向为 16,其它保持不变。
2) 太阳辐射区域工具计算的太阳辐射量为假设该地区全部晴天可以获得的太阳辐射量。实际太阳辐射受晴、阴、云、雨等随机因素的影响。为了简化计算,我们把天气分为晴天和非晴天,晴天可以获得太阳辐射区域工具计算的全部太阳辐射量,非晴天则无法获得任何太阳辐射量。
3) 该地区 2021 年每月晴天数据如下:月份 1月 2月 3月 4月 5月 6月 7月 8月 9月 10月 11月 12月 晴天 25 20 24 23 20 15 22 26 27 20 25 26 4) 根据每月晴天数据,创建房屋屋顶范围内 2021 年每月太阳辐射栅格数据,像元值表示瓦特每小时每平方米(WH/m2)。请命名为“太阳辐射 1 月”、“太阳辐射 2 月”、“太阳辐射 3 月”,以此类推,完成下表。
月份 1月 2月 3月 4月 5月 6月 7月 8月 9月 10月 11月 12月 低值 高值 -
计算 8 月份可用房屋屋顶范围内的太阳辐射量,在环境设置中,将 Building 作为掩膜(30分)
考虑到太阳能电池板的能效,符合以下条件的屋顶才适合安装电池板。
1) 如果屋顶坡度小于等于 15 度,则任何屋顶方向都可以安装太阳能电池板。
2) 如果屋顶坡度大于 15 度,且小于等于 30 度,则屋顶朝向不可以朝北(北方为大于等于 337.5,或者小于等于 22.5 度)。
3) 如果屋顶坡度大于 30 度,且小于等于 45 度,则屋顶朝向不可以朝北、东北和西北(东北为大于 22.5 度且小于等于 67.5 度,西北为大于等于 292.5 度且小于 337.5 度)。
4) 创建 2021 年 8 月,符合以上条件要求的房屋屋顶范围内的太阳辐射量栅格数据,命名为“可用辐射 8 月”,像元值表示瓦特每小时每平方米(WH/m2)。 -
计算 8 月份每栋房屋可接收的太阳辐射量,在环境设置中,将 Building 作为掩膜
1) 为 Building 添加属性字段【可用面积】,计算每栋房屋屋顶可用太阳辐射区域的面积。 (10分)
\(\quad\quad\)面积单位为平方米,保留 2 位小数位数。1 平方英尺 = 0.093 平方米。
2) 为 Building 添加属性字段【可用辐射 8 月】,计算每栋房屋屋顶 8 月可接受太阳辐射量。 (10月)
\(\quad\quad\)a) 如果每栋房屋可用面积小于 25 平方米,则该栋房屋不适合安装太阳能电池板,无需计算可接受太阳辐射量。
\(\quad\quad\)b) 可用辐射量单位为千瓦时(KWH),保留 2 位小数位数。
第Ⅱ部分:参考解题过程
一、 导入原始数据并配置地理处理空间
\(\quad\quad\)打开【地理处理】|【环境】,设置工作空间(我的数据文件夹就放在桌面上),并设置输出栅格像元大小与 DTM 一致。
\(\quad\quad\)导入原始数据 “Building.shp” 、 “DTM.tif” 以及 “DSM.tif” 进入显示窗口中,操作前查看数据的坐标系是一个好习惯,发现三者的投影坐标系均一致,无需特别处理。
二、根据建筑物修正DTM
-
计算【基本高度】字段
\(\quad\quad\)我们首先需要计算建筑物底部的高程值,其数值为房屋区域 DTM 的平均值,可以想到需要使用分区统计之类的工具,又涉及属性表添加字段,所以选定【ArcToolBox】|【Spatial Analyst】|【区域分析】|【以表格显示分区统计】该工具。
\(\quad\quad\)首先打开 Building 的属性表并且添加双精度类型字段 “base_tall” ,我选的精度是 12,小数位数被规定为 2。添加字段之后在属性表中右键字段名称更改字段属性,将【别名】设置为中文名称 “基本高度”。\(\quad\quad\)然后打开【以表格显示分区统计】该工具,以 Building 为分区的依据对每栋房屋范围内的 DTM 进行均值统计,该工具的输出结果是一个表格,【内容列表】需要 “按源列出” 才可以看到并打开。工具详细设置界面如下:
\(\quad\quad\)上述操作的输出结果是一张表格,内部有 “ROWID”、“BUILDINGBM”、“COUNT”、“AREA”、“MEAN” 这 5 个字段。将 Building.shp 与表格 mean_dtm 根据 “BUILDINGBM” 字段建立连接。然后使用【字段计算器】让 “mean_dtm:MEAN” 对 “Building:base_tall”进行赋值。
【基本高度】字段已经赋值完毕,结果如下,解除 Building 与 mean_dtm 之间的连接。
-
修正 DTM 数据
\(\quad\quad\)使用【ArcToolBox】|【Spatial Analyst】|【区域分析】|【分区统计】工具,设置界面如下:\(\quad\quad\)输出结果如下,可以看出,每栋房屋的 mean_DTM 是均匀一致的。
\(\quad\quad\)然后将 mean_DTM.tif 与 DTM 【镶嵌至新栅格】,需要注意镶嵌的顺序:如果输入数据中 mean_DTM 在前,那么镶嵌方式选择为 “FIRST”,反之则选择镶嵌方式为 “LAST”,因为我们的目的是希望 mean_DTM 代替局部的 DTM,来产生 “DTM 修正”。除此之外,镶嵌至新栅格需要设置波段数、像元大小、像素类型等参数,这些参数均参照 DTM。
\(\quad\quad\)DTM 修正结果如下:
三、计算每栋房屋的其它基础信息
\(\quad\quad\)首先对 Building 添加字段【max_tall】和【build_tall】,保留 2 位小数,并且将别名分别设置为【最大高度】和【建筑物高度】。
-
计算房屋的最大高度
\(\quad\quad\)获取建筑物的最大高度与上一步骤获取建筑物底部高程值步骤一致,先使用【以表格显示分区统计】工具,对 DSM 的每一栋房屋屋顶范围内的最大值进行统计,输出表格后与 Building 建立连接,然后对字段【最大高度】进行赋值,最后解除连接。\(\quad\quad\)连接属性表以及使用【字段计算器】进行赋值的步骤已省略,计算【最大高度】与之前计算【基本高度】的步骤相似,可参照前文。
-
计算房屋的建筑物高度
\(\quad\quad\)建筑物高度等于最大高度减去基本高度,因此有:\(build\_tall = max\_tall - base\_tall\)。\(\quad\quad\)最终得到的【建筑物高度】连同上一步得到的【最大高度】如下图所示:
-
获取指定房屋的屋顶形态
\(\quad\quad\)由题目提示可以知道,判断房屋屋顶形态的依据是坡向。因此先计算建筑物屋顶范围内的坡向。输入数据为 DSM。\(\quad\quad\)先不忙点击【确定】执行工具,点击【环境】设置【掩膜】为 Building,这样可以只提取屋顶范围内的坡向。
\(\quad\quad\)完成设置并执行【坡向】工具后得到如下结果:
\(\quad\quad\)到这里可以使用 Building.shp 辅助定位特定编号的房屋。使用【按属性选择】工具,根据 Building.shp 的 BuildingBM 属性分别先后筛选编号为 “2116”、“2156”、“2161”、“2165”、“2171” 这 5 栋房屋。然后根据 Building.shp 的高亮提示,结合 “Aspect_DSM.tif” 可以判断该房屋的屋顶坡向分布情况,如上图所示,编号为 “2116” 的房屋屋顶为双坡屋顶,并且朝向为西北与东南。最终完成下表:
序号 BuildingBM 屋顶形态 屋顶朝向 1 2116 双坡屋顶 西北、东南 2 2156 四坡屋顶 东、南、西、北 3 2161 四坡屋顶 西北、东北、东南、西南 4 2165 四坡屋顶 西北、东北、东南、西南 5 2171 四坡屋顶 西北、东北、东南、西南
四、创建房屋屋顶区域 2021 年每月预计获得太阳辐射量栅格数据,在环境设置中,将 Building 作为掩膜
\(\quad\quad\)计算房屋屋顶区域 12 个月每一个月的太阳辐射量栅格数据主要使用【太阳辐射区域】工具,该工具位于【ArcToolBox】|【Spatial Analyst】|【太阳辐射】目录下,并且题目中对于参数已经做了限定,其余参数保持默认即可。尝试了几次该工具后我发现,其计算速度比较慢,如果手动计算 12 次要等很久,于是想到了使用 Arcpy 的 Python 模块来编写脚本。这样在 Python 代码运行的过程中可以继续使用 ArcMap 进行操作。
\(\quad\quad\)使用 ArcGIS 的 Python 脚本功能主要有 2 种方式:
- 使用 ArcMap 主界面的 Python 脚本编辑器;
- 使用外部的 Python 开发环境引入 Arcpy 模块。
\(\quad\quad\)这里我选择后者以方便调试。我打开了 Pycharm,新建 Python 项目,然后设置编译器的位置为 ArcGIS 安装位置处的 Python 2.7.18 环境,开发环境就这样搭建好啦。
\(\quad\quad\)那么接下来的问题是如何编写代码,由于 ArcMap 的【模型构建器】具有一键导出模型为 Python 脚本的功能,因此我先在 ArcToolBox 中拖动【太阳辐射区域】工具进入 【模型构建器】编辑界面,然后设置参数,然后一键生成代码,ArcMap 这个功能可以说是我的福音啊哈哈哈。
\(\quad\quad\)上图中可以看出我事先设置模型属性中的掩膜为 Building.shp。太阳辐射区域的设置界面如下所示(以计算屋顶区域 2021 年 1 月份太阳辐射为例):
\(\quad\quad\)设置完工具之后回到【模型构建器】的编辑界面,点击【导出至 Python 脚本】:
\(\quad\quad\)当然,导出的脚本仍然有一些不完善的地方,例如需要针对每个月晴天天数对结果进行乘法修正,需要得到每个月太阳辐射的低值和高值。这时可以参照 ArcGIS 的帮助文档修改代码,帮助文档对 Python 脚本也是有详细介绍的:
\(\quad\quad\)最后,这是我的脚本的最终版本:
# -*- coding:UTF-8 -*-
# Name: AreaSolarRadiation_example02.py
# Description: Derives incoming solar radiation from a raster surface.
# Outputs a global radiation raster and optional direct, diffuse and direct duration rasters
# for a specified time period. (April to July).
#
# Requirements: Spatial Analyst Extension
# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *
import numpy as np
# Set environment settings
env.workspace = "C:/Users/Fangh/Desktop/data" # 工作空间,需更改
# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
# Set Geoprocessing environments
env.cellSizeProjectionMethod = "CONVERT_UNITS"
env.cellSize = "DTM.tif" # 像元大小与DTM.tif一致
env.mask = "Building.shp" # 设置屋顶区域Building.shp为掩膜
# Set local variables
inRaster = "C:/Users/Fangh/Desktop/data/DSM.tif"
latitude = 41.77686382537 # 纬度
skySize = 200 # 天空大小
timeConfig = [TimeMultipleDays(2021, 1, 31), # Jenuary
TimeMultipleDays(2021, 32, 59), # February
TimeMultipleDays(2021, 60, 90), # March
TimeMultipleDays(2021, 91, 120), # April
TimeMultipleDays(2021, 121, 151), # May
TimeMultipleDays(2021, 152, 181), # June
TimeMultipleDays(2021, 182, 212), # July
TimeMultipleDays(2021, 213, 243), # August
TimeMultipleDays(2021, 244, 273), # September
TimeMultipleDays(2021, 274, 304), # October
TimeMultipleDays(2021, 305, 334), # November
TimeMultipleDays(2021, 335, 365)] # December
dayInterval = 14 # 间隔天数
hourInterval = 0.5 # 间隔小时数
zFactor = 1 # 地形参数/Z因子
calcDirections = 16 # 地形参数/计算方向
zenithDivisions = 8 # 辐射参数/天顶分割
azimuthDivisions = 8 # 辐射参数/方位角分割
diffuseProp = 0.3 # 辐射参数/散射比例
transmittivity = 0.5 # 辐射参数/透射率
outDirectRad = "" # 无需"输出直接辐射栅格数据"
outDiffuseRad = "" # 无需"输出散射辐射栅格数据"
outDirectDur = "" # 无需"直接辐射持续时间栅格数据"
# Sunny days in every month
SunnyDay = [25.0/31.0, 20.0/28.0, 24.0/31.0,
23.0/30.0, 20.0/31.0, 15.0/30.0,
22.0/31.0, 26.0/31.0, 27.0/30.0,
20.0/31.0, 25.0/30.0, 26.0/31.0]
# statistic the maximum and minimum of every raster
mni_max_value = np.zeros([2, 12])
outputcsv = r"C:/Users/Fangh/Desktop/data/statics.csv"
outputfile = open(outputcsv, 'w')
for i in range(12):
# Execute AreaSolarRadiation
outGlobalRad = AreaSolarRadiation(inRaster, latitude, skySize, timeConfig[i],
dayInterval, hourInterval, "NOINTERVAL", zFactor, "FROM_DEM",
calcDirections, zenithDivisions, azimuthDivisions, "UNIFORM_SKY",
diffuseProp, transmittivity, outDirectRad, outDiffuseRad, outDirectDur)
# Consider the sunny day
outGlobalRad *= SunnyDay[i]
# Save the output
outGlobalRad.save("C:/Users/Fangh/Desktop/data/Rad_"+str(i+1)+".tif")
# Get the minimum and maximum
minimumvalueInfo = arcpy.GetRasterProperties_management(outGlobalRad, "MINIMUM")
mni_max_value[0, i] = minimumvalueInfo.getOutput(0)
maximumvalueInfo = arcpy.GetRasterProperties_management(outGlobalRad, "MAXIMUM")
mni_max_value[1, i] = maximumvalueInfo.getOutput(0)
# Write the outputfile
outputfile.write('Month'+','+'January'+','+'February'+','+'March'+','+'April'+','+'May'+','+'June'+','+'July'+','+'August'+','+'September'+','+'October'+','+'November'+','+'December'+'\n')
outputfile.write('Minimum')
for i in range(12):
outputfile.write(','+str(mni_max_value[0, i]))
outputfile.write('\n')
outputfile.write('Maximum')
for i in range(12):
outputfile.write(','+str(mni_max_value[1, i]))
outputfile.write('\n')
outputfile.close()
print("All Job Were Done!")
\(\quad\quad\)这份代码在计算出 12 个月份的屋顶区域太阳辐射量栅格数据之后,生成了一份 statistic.csv 逗号分隔值文件(用 Excel 打开),与数据存放在相同文件夹之下。于是得到房屋屋顶区域 2021 年每月的太阳辐射量栅格数据以及每一个月的低值和高值如下表所示:
月份 | 1月 | 2月 | 3月 | 4月 | 5月 | 6月 | 7月 | 8月 | 9月 | 10月 | 11月 | 12月 |
---|---|---|---|---|---|---|---|---|---|---|---|---|
低值 | 4.23221254349 | 5.31080770493 | 9.62105178833 | 12.2721033096 | 12.555768013 | 10.0049877167 | 14.2893257141 | 14.9097099304 | 12.1840105057 | 6.18915987015 | 4.79556846619 | 3.67597270012 |
高值 | 39828.0976563 | 48683.125 | 85742.09375 | 106786.398438 | 108502.523438 | 85074.40625 | 123380.328125 | 129046.210938 | 107483.84375 | 56247.484375 | 44828.5390625 | 28945.9609375 |
\(\quad\quad\)最后,手动拖入 12 份栅格数据进入 ArcMap 查看有无错误,并修改中文名称为 “太阳辐射 * 月”。结果如下:
五、计算 8 月份可用房屋屋顶范围内的太阳辐射量,在环境设置中,将 Building 作为掩膜
\(\quad\quad\)这一要求可以使用 “与或非” 逻辑解决,借助【ArcToolBox】|【Spatial Analyst】|【地图代数】|【栅格计算器】工具,当坡度位于于 \([0, 15]\) 度时,可以在任意坡向方向安装太阳能电池板,或者坡度位于 \((15, 30]\) 度之间时,且坡向位于 \((22.5, 337.5)\) 度之间,或者坡度位于 \((30, 45]\) 度之间时,且坡向位于 \((67.5, 292.5)\) 度之间。
-
计算坡度
\(\quad\quad\)打开【ArcToolBox】|【Spatial Analyst】|【表面分析】|【坡度】,输入 DSM,设置 Building 为掩膜: -
计算 8 月份可用房屋屋顶范围内的太阳辐射量
\(\quad\quad\)表达式的的全文如下(其中“太阳辐射量\”表示图层组):
( ("Slope_DSM.tif" <= 15) + ("Slope_DSM.tif" > 15)*("Slope_DSM.tif" <= 30)*("Aspect_DSM.tif" > 22.5)*("Aspect_DSM.tif" < 337.5) + ("Slope_DSM.tif" > 30)*("Slope_DSM.tif" <= 45)*("Aspect_DSM.tif" > 67.5)*("Aspect_DSM.tif" < 292.5) ) * "太阳辐射量\太阳辐射8月.tif"
-
将 0 值转化为 NoData
\(\quad\quad\)因为上述逻辑表达式的结果是 “0” 或者 “1” 与 “太阳辐射量\太阳辐射8月.tif” 进行乘法操作,所以结果中出现了 “0” ,现在需要将 8 月份可用房屋屋顶范围内的零太阳辐射值变为 NoData。仍然使用【栅格计算器】,使用 “SetNull” 函数,使用之前查看帮助文档中的语法:\(\quad\quad\)因此,表达式为:
SetNull("batter_rad_8.tif", "batter_rad_8.tif", "Value = 0")
\(\quad\quad\)其中,“batter_rad_8.tif” 为上一步的计算结果,此步骤的输出结果为 “valid_rad_8.tif”,修改中文名称为 “可用辐射 8 月”。
六、计算 8 月份每栋房屋可接收的太阳辐射量,在环境设置中,将 Building 作为掩膜
\(\quad\quad\)首先添加字段【可用面积】与【可用辐射量 8 月】如下所示,小数位数为 2 位,数据类型为双精度,精度为 12 位:
-
计算【可用面积】
\(\quad\quad\)上一步骤已经得到可用房屋屋顶范围内的太阳辐射量,此步骤需要计算其面积,而栅格数据计算面积通常是像元总数乘以单个像元大小。因此先查看单个像元大小,发现为 “1, 1”,且线性单位为英尺,因此一个像元大小为 1 平方英尺。\(\quad\quad\)对 “可用辐射 8 月” 进行区域求和操作,就可以得到每一栋房屋可用屋顶范围内的像元总数,在栅格数据中,它通常是由 “Count” 字段来表示。设置界面如下:
\(\quad\quad\)打开 “area_rad” 表格可以发现,除了 “COUNT” 字段,分区统计工具自动生成了一个 “AREA” 属性字段,并且两者相等没有差别,单位是平方英尺,后面我们使用 “AREA” 继续计算。此外,“SUM” 字段表示可用房屋屋顶范围内的太阳辐射量之和(此时 SUM 没有乘以面积,因此没有具体意义,但在后续过程会用到)。
\(\quad\quad\)连接 Building 与 area_rad 表格:
\(\quad\quad\)使用【字段计算器】将面积的单位从平方英尺转换为平方米,表达式如下所示:
-
计算【可用辐射量 8 月】
\(\quad\quad\)首先【按属性选择】可用面积大于 25 平方米的房屋屋顶,然后在属性表中点击【显示所选记录】按钮,如下图标注:\(\quad\quad\)然后右键字段名称【可用辐射量 8 月】,打开【字段计算器】,我理解是这里使用微积分的思想计算每栋房屋的太阳辐射量,即:
\[\]\[ \]\(\quad\quad\)输出结果后取消【按属性选择】,回到正常的属性表界面。可以看到如下结果:
\(\quad\quad\)最后取消 Building 与 area_rad 之间的连接。
结尾:参赛收获及数据下载地址
\(\quad\quad\)这道题目做下来觉得操作量挺大的,耗时也比较长,例如使用 Arcpy 那一部分在我的电脑上运行了半个小时左右。在规定时间内完成全部题目很有难度。当然,如果操作者对 ArcGIS 分析流程更熟悉的话就容易很多,我还要继续加油。
\(\quad\quad\)数据下载地址(提取密码:wx7d)
标签:Building,GIS,试题,DTM,辐射量,解题,屋顶,quad,tif 来源: https://www.cnblogs.com/XiongHaiyang/p/16512709.html