基于大区域栅格数字高程模型的地表表面面积计算方法

摘要

本发明公开了一种基于大区域栅格数字高程模型的地表表面面积计算方法，它首先将原始栅格数字高程模型分割为若干个栅格数据，然后计算每个栅格数据的地表表面面积，最后累加所有栅格数据的表面面积得到原始栅格数字高程模型的表面面积。由于采用了大数据分块的策略，本发明能够高效出计算大区域栅格数字高程模型的地表表面面积，并且在像素级别尺度构建三角网，以及精确处理边缘像元，比现有计算方法计算精度更高，不仅提高了表面面积算法的精度，而且解决了计算大数据量数字高程模型表面面积计算机内存不足的问题。

说明书

技术领域

本发明涉及一种基于大区域栅格数字高程模型的地表表面面积计算方法，属于地理信息空间分析技术领域。

背景技术

随着三维技术的飞速发展，三维地形数据大量出现在我们的工作和生活中，其中三维地形数字分析提供给我们丰富可靠的数据，其中地表表面面积是其中经常用到的一种数字分析。地表表面面积很难直接量测和计算，一般使用投影面积代替表面面积。在地形起伏较小的平原地区，地表真实面积和投影面积相差不大，但在地形起伏比较复杂的山区，两者差异明显。

目前有很多方法计算投影面积，运用投影面和坡度正割值计算表面面积、逐步分割多边形的方法计算表面面积、将一个像元分割为2个三角形统计DEM的面积等，这些方法存在一些缺点：计算精度不够或者没有考虑大数据量下的计算问题，并不能很好满足我们高效精确计算栅格数字高程模型的地表表面面积计的需求。

发明内容

本发明的目的是提供一种大区域栅格数字高程模型的地表表面面积计算方法，提高现有表面面积算法的精度，以及解决计算大数据量数字高程模型表面面积计算机内存不足的问题。

本发明解决其技术问题所采取的技术方案是：一种基于大区域栅格数字高程模型的地表表面面积计算方法，其特征是：它首先将原始栅格数字高程模型分割为若干个栅格数据，然后计算每个栅格数据的地表表面面积，最后得到原始栅格数字高程模型的表面面积。

进一步地，所述每个栅格数据的地表表面面积的计算具体过程为：

S2，查找每个栅格数据的边界像元并计算边界像元面积；

S3，对每个栅格数据构建三角网并计算三角网面积；

将步骤S2中计算得到的边界像元面积与步骤S3中计算得到的三角网面积进行相加，得到该栅格数据的地表表面面积。

进一步地，所述的原始栅格数字高程模型的表面面积通过累加所有栅格数据的表面面积得到。

进一步地，所述步骤S2的具体步骤如下：

S2.1，将每个栅格数据转化为高程数组Array[height][width]，with为栅格数据的水平方向像元个数，height为栅格数据的垂直方向像元个数，Array[x][y]表示栅格数据第x行第y列像元的高程值，如果像元[x,y]为空值，则取Array[x][y]为浮点数的最小值，记作noData；通过x,y和栅格数据左上角坐标计算出第x行第y列像元中心点的地理坐标；

记计算像元8邻域像元个数的函数为NeighbourCount，NeighbourCount(x,y)为计算位置[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]和[x-1,y-1]存在非空像元的个数，NeighbourCount(x,y)<8则代表像元(x,y)为边界像元；

声明边界像元面积字典Dictionary，键为边界点坐标，值为边界像元面积；Dictionary[x,y]＝S表示坐标为[x,y]的边界像元的面积为S；

获取一个像元的平面面积，记为S₀；

设当前行i＝0,当前列j＝0：

S2.2，比较i和height的值：

如果i＝height，则退出计算；

如果i<height，则进入下一步；

S2.3，比较j和width的值：

如果j＝width，i＝i+1，j＝0，返回步骤S2.2；

如果j<width，则进入下一步；

S2.4，判断像元(i,j)是否为边界像元：

如果NeighbourCount(i,j)＝8，j，则j＝j+1，返回步骤S2.3；

如果NeighbourCount(i,j)<8，则进入下一步；

S2.5，在边界像元面积字典中加入像元[i,j]，并初始化面积为0，即Dictionary[i,j]＝0；

S2.6，j＝j+1,返回步骤S2.3。

进一步地，所述步骤S3的具体步骤如下：

S3.1，获取步骤S2.1高程数组Array[height][width]，像元空值noData和边界像元面积字典Dictionary；

构建三角网：提取每个像元的中心点作为顶点，每个顶点与相邻的8个像元中心点中的2个形成三角形；

计算边界像元面积：每个边界像元按照对角线和对边中位线分成8份，每一份面积为S₀/8，计算在构建三角网之后剩余的份数，然后计算剩余份数的平面面积，即得到该边界像元的面积；

设当前行i＝0，当前列j＝0，栅格数据的地表表面面积为S_all；

S3.2，比较i和height的值：

如果i＝height，则进入步骤S3.8；

如果i<height，则进入下一步；

S3.3，比较j和width的值：

如果j＝width，i＝i+1，j＝0，返回步骤S3.2；

如果j<width，则进入下一步；

S3.4，比较Array[i][j]和noData的值：

如果Array[i][j]＝noData，j＝j+1，返回步骤S3.3；

如果Array[i][j]>noData，则进入下一步；

S3.5，判断第i-1行数据是否存在：

如果i-1<0，则进入步骤S3.6；

如果i-1>＝0，则进入下一步；

S3.5.1，判断像元Array[i-1][j+1]、Array[i][j+1]和Array[i-1][j]：

如果Array[i-1][j+1]和Array[i][j+1]存在且不等于noData，而且Array[i-1][j]等于noData，则以索引为[i,j]、[i,j+1]和[i-1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.5.3；

S3.5.2判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i，j+1]，则Dictionary[i,j+1]＝Dictionary[i,j+1]+2*S₀/8；

如果Dictionary中存在像元[i-1，j+1]，则Dictionary[i-1，j+1]＝Dictionary[i-1，j+1]+S₀/8；

S3.5.3，判断像元Array[i-1][j]、Array[i-1][j+1]和Array[i][j+1]：

如果Array[i-1][j]和Array[i-1][j+1]存在且不等于noData，而且Array[i][j+1]不存在或者等于noData，则以索引为[i,j]、[i-1,j]和[i-1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.6

S3.5.4判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i-1,j]，则Dictionary[i-1,j]＝Dictionary[i-1,j]+2*S₀/8；

如果Dictionary中存在像元[i-1,j+1]，则Dictionary[i-1,j+1]＝Dictionary[i-1，j+1]+S₀/8；

S3.6，判断第i+1行数据是否存在：

如果i+1>＝height，则退出计算；

如果i+1<height，则进入下一步；

S3.6.1，判断像元Array[i][j+1]和Array[i+1][j+1]：

如果Array[i][j+1]和Array[i+1][j+1]存在且不等于noData，则以索引为[i,j]、[i,j+1]和[i+1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.6.3；

S3.6.2，判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i,j+1]，则Dictionary[i,j+1]＝Dictionary[i,j+1]+2*S₀/8；

如果Dictionary中存在像元[i+1,j+1]，则Dictionary[i+1,j+1]＝Dictionary[i+1,j+1]+S₀/8；

S3.6.3，判断像元Array[i+1][j]和Array[i+1][j+1]：

如果Array[i+1][j]和Array[i+1][j+1]存在且不等于noData，则以索引为[i,j]、[i+1,j]和[i+1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.7；

S3.6.4，判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i+1,j]，则Dictionary[i+1,j]＝Dictionary[i+1,j]+2*S₀/8；

如果Dictionary中存在像元[i+1,j+1]，则Dictionary[i+1,j+1]＝Dictionary[i+1,j+1]+S₀/8；

S3.7，j＝j+1,返回步骤S3.3；

S3.8，对边界像元面积字典Dictionary中所有项求和得到S_dic，并加入到S_all，即得到整个栅格数据的地表表面面积S_all＝S_all+S_dic。

本发明提供的另一种基于大区域栅格数字高程模型的地表表面面积计算方法，其特征是：包括以下步骤：

S1，对原始栅格数字高程模型进行分块处理，生成若干个栅格数据；

S2，查找每个栅格数据的边界像元并计算边界像元面积；

S3，对每个栅格数据构建三角网并计算三角网面积；

S4，计算每个栅格数据的地表表面面积；

S5，累加所有栅格数据的表面面积即得到原始栅格数字高程模型的表面面积。

进一步地，在步骤S1中，对原始栅格数字高程模型进行格网化分割，生成若干个大小相同的栅格数据。

进一步地，在步骤S4中，每个栅格数据的地表表面面积为每块栅格数据的三角网面积和边界像元面积之和。

进一步地，所述步骤S2的具体步骤如下：

S2.1，将每个栅格数据转化为高程数组Array[height][width]，with为栅格数据的水平方向像元个数，height为栅格数据的垂直方向像元个数，Array[x][y]表示栅格数据第x行第y列像元的高程值，如果像元[x,y]为空值，则取Array[x][y]为浮点数的最小值，记作noData；通过x,y和栅格数据左上角坐标计算出第x行第y列像元中心点的地理坐标；

记计算像元8邻域像元个数的函数为NeighbourCount，NeighbourCount(x,y)为计算位置[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]和[x-1,y-1]存在非空像元的个数，NeighbourCount(x,y)<8则代表像元(x,y)为边界像元；

声明边界像元面积字典Dictionary，键为边界点坐标，值为边界像元面积；Dictionary[x,y]＝S表示坐标为[x,y]的边界像元的面积为S；

获取一个像元的平面面积，记为S₀；

设当前行i＝0,当前列j＝0：

S2.2，比较i和height的值：

如果i＝height，则退出计算；

如果i<height，则进入下一步；

S2.3，比较j和width的值：

如果j＝width，i＝i+1，j＝0，返回步骤S2.2；

如果j<width，则进入下一步；

S2.4，判断像元(i,j)是否为边界像元：

如果NeighbourCount(i,j)＝8，j，则j＝j+1，返回步骤S2.3；

如果NeighbourCount(i,j)<8，则进入下一步；

S2.5，在边界像元面积字典中加入像元[i,j]，并初始化面积为0，即Dictionary[i,j]＝0；

S2.6，j＝j+1,返回步骤S2.3。

进一步地，所述步骤S3的具体步骤如下：

S3.1，获取步骤S2.1高程数组Array[height][width]，像元空值noData和边界像元面积字典Dictionary；

构建三角网：提取每个像元的中心点作为顶点，每个顶点与相邻的8个像元中心点中的2个形成三角形；

计算边界像元面积：每个边界像元按照对角线和对边中位线分成8份，每一份面积为S₀/8，计算在构建三角网之后剩余的份数，然后计算剩余份数的平面面积，即得到该边界像元的面积；

设当前行i＝0，当前列j＝0，栅格数据的地表表面面积为S_all；

S3.2，比较i和height的值：

如果i＝height，则进入步骤S3.8；

如果i<height，则进入下一步；

S3.3，比较j和width的值：

如果j＝width，i＝i+1，j＝0，返回步骤S3.2；

如果j<width，则进入下一步；

S3.4，比较Array[i][j]和noData的值：

如果Array[i][j]＝noData，j＝j+1，返回步骤S3.3；

如果Array[i][j]>noData，则进入下一步；

S3.5，判断第i-1行数据是否存在：

如果i-1<0，则进入步骤S3.6；

如果i-1>＝0，则进入下一步；

S3.5.1，判断像元Array[i-1][j+1]、Array[i][j+1]和Array[i-1][j]：

如果Array[i-1][j+1]和Array[i][j+1]存在且不等于noData，而且Array[i-1][j]等于noData，则以索引为[i,j]、[i,j+1]和[i-1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.5.3；

S3.5.2判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i，j+1]，则Dictionary[i,j+1]＝Dictionary[i,j+1]+2*S₀/8；

如果Dictionary中存在像元[i-1，j+1]，则Dictionary[i-1，j+1]＝Dictionary[i-1，j+1]+S₀/8；

S3.5.3，判断像元Array[i-1][j]、Array[i-1][j+1]和Array[i][j+1]：

如果Array[i-1][j]和Array[i-1][j+1]存在且不等于noData，而且Array[i][j+1]不存在或者等于noData，则以索引为[i,j]、[i-1,j]和[i-1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.6

S3.5.4判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i-1,j]，则Dictionary[i-1,j]＝Dictionary[i-1,j]+2*S₀/8；

如果Dictionary中存在像元[i-1,j+1]，则Dictionary[i-1,j+1]＝Dictionary[i-1，j+1]+S₀/8；

S3.6，判断第i+1行数据是否存在：

如果i+1>＝height，则退出计算；

如果i+1<height，则进入下一步；

S3.6.1，判断像元Array[i][j+1]和Array[i+1][j+1]：

如果Array[i][j+1]和Array[i+1][j+1]存在且不等于noData，则以索引为[i,j]、[i,j+1]和[i+1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.6.3；

S3.6.2，判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i,j+1]，则Dictionary[i,j+1]＝Dictionary[i,j+1]+2*S₀/8；

如果Dictionary中存在像元[i+1,j+1]，则Dictionary[i+1,j+1]＝Dictionary[i+1,j+1]+S₀/8；

S3.6.3，判断像元Array[i+1][j]和Array[i+1][j+1]：

如果Array[i+1][j]和Array[i+1][j+1]存在且不等于noData，则以索引为[i,j]、[i+1,j]和[i+1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.7；

S3.6.4，判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i+1,j]，则Dictionary[i+1,j]＝Dictionary[i+1,j]+2*S₀/8；

如果Dictionary中存在像元[i+1,j+1]，则Dictionary[i+1,j+1]＝Dictionary[i+1,j+1]+S₀/8；

S3.7，j＝j+1,返回步骤S3.3；

S3.8，对边界像元面积字典Dictionary中所有项求和得到S_dic，并加入到S_all，即得到整个栅格数据的地表表面面积S_all＝S_all+S_dic。

本发明的有益效果是：本发明首先对栅格数字高程模型分块；然后通过查找每块数据边界像元以及对每块数据构建三角网来计算每块数据的地表表面面积；最后对所有分块数据面积求和，从而得到原始栅格数字高程模型的表面面积。由于采用了大数据分块的策略，本发明能够高效出计算大区域栅格数字高程模型的地表表面面积，并且在像素级别尺度构建三角网，以及精确处理边缘像元，比现有计算方法计算精度更高，不仅提高了表面面积算法的精度，而且解决了计算大数据量数字高程模型表面面积计算机内存不足的问题。

附图说明

图1是本发明的一种方法流程图；

图2是本发明的另一种方法流程图；

图3是本发明对原始栅格数字高程模型进行格网化分割示意图；

图4是本发明所述网格数据的示意图；

图5(a)为本发明所述网格数据中像元与其正左方和正下方像元构建三角形的示例图；

图5(b)为本发明所述网格数据中像元与其与正下方和右下方像元构建三角形示例图；

图5(c)为本发明所述网格数据中像元与其与正下方和左下方像元构建三角形示例图；

图5(d)为本发明所述网格数据中像元与其与正右方和正下方像元构建三角形的示例图；

图6为本发明的栅格数字高程模型与生成的三角网的二维示意图；

图7为本发明的栅格数字高程模型的三维示意图；

图8为本发明的生成三角网的三维示意图。

具体实施方式

为能清楚说明本方案的技术特点，下面通过具体实施方式，并结合其附图，对本发明进行详细阐述。下文的公开提供了许多不同的实施例或例子用来实现本发明的不同结构。为了简化本发明的公开，下文中对特定例子的部件和设置进行描述。此外，本发明可以在不同例子中重复参考数字和/或字母。这种重复是为了简化和清楚的目的，其本身不指示所讨论各种实施例和/或设置之间的关系。应当注意，在附图中所图示的部件不一定按比例绘制。本发明省略了对公知组件和处理技术及工艺的描述以避免不必要地限制本发明。

针对现有的评价公共设施均衡性方法的不足，为了实现城市公共设施的公平配置和分布均衡的目标，本发明合理地考虑了人与公共设施之间作用的跨区域性问题，既适用于区域之间也适用于区域内部，直观易懂，能够准确地评价一个城市的公共设施的服务辐射范围、服务能力和均衡性。

如图1、图3至图8所示，本发明的一种基于大区域栅格数字高程模型的地表表面面积计算方法，它首先将原始栅格数字高程模型分割为若干个栅格数据(栅格数据的大小可以根据计算机性能调整，一般取512*512)，然后计算每个栅格数据的地表表面面积，最后累加所有栅格数据的表面面积得到原始栅格数字高程模型的表面面积。

进一步地，所述每个栅格数据的地表表面面积的计算具体过程为：

S2，查找每个栅格数据的边界像元并计算边界像元面积；

S3，对每个栅格数据构建三角网并计算三角网面积；

将步骤S2中计算得到的边界像元面积与步骤S3中计算得到的三角网面积进行相加，得到该栅格数据的地表表面面积。

如图2至图8所示，本发明提供的另一种基于大区域栅格数字高程模型的地表表面面积计算方法，它包括以下步骤：

S1，对原始栅格数字高程模型进行分块处理，生成若干个栅格数据(栅格数据的大小可以根据计算机性能调整，一般取512*512)；

S2，查找每个栅格数据的边界像元并计算边界像元面积；

S3，对每个栅格数据构建三角网并计算三角网面积；

S4，计算每个栅格数据的地表表面面积；

S5，重复步骤S2至步骤S4直至遍历完所述生成的栅格数据为止，并累加所有栅格数据的表面面积即得到原始栅格数字高程模型的表面面积。

进一步地，在步骤S1中，对原始栅格数字高程模型进行格网化分割，生成若干个大小相同的栅格数据。

进一步地，在步骤S4中，每个栅格数据的地表表面面积为每块栅格数据的三角网面积和边界像元面积之和。

进一步地，在上述两种计算方法中，所述步骤S2的具体步骤如下：

S2.1，将每个栅格数据转化为高程数组Array[height][width]，with为栅格数据的水平方向像元个数，height为栅格数据的垂直方向像元个数，Array[x][y]表示栅格数据第x行第y列像元的高程值，如果像元[x,y]为空值，则取Array[x][y]为浮点数的最小值，记作noData；通过x,y和栅格数据左上角坐标计算出第x行第y列像元中心点的地理坐标；

记计算像元8邻域像元个数的函数为NeighbourCount，NeighbourCount(x,y)为计算位置[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]、[x-1,y-1]和[x-1,y-1]存在非空像元的个数，NeighbourCount(x,y)<8则代表像元(x,y)为边界像元；

声明边界像元面积字典Dictionary，键为边界点坐标，值为边界像元面积；Dictionary[x,y]＝S表示坐标为[x,y]的边界像元的面积为S；

获取一个像元的平面面积，记为S₀；

设当前行i＝0,当前列j＝0：

S2.2，比较i和height的值：

如果i＝height，则退出计算；

如果i<height，则进入下一步；

S2.3，比较j和width的值：

如果j＝width，i＝i+1，j＝0，返回步骤S2.2；

如果j<width，则进入下一步；

S2.4，判断像元(i,j)是否为边界像元：

如果NeighbourCount(i,j)＝8，j，则j＝j+1，返回步骤S2.3；

如果NeighbourCount(i,j)<8，则进入下一步；

S2.5，在边界像元面积字典中加入像元[i,j]，并初始化面积为0，即Dictionary[i,j]＝0；

S2.6，j＝j+1,返回步骤S2.3。

进一步地，在上述两种计算方法中，所述步骤S3的具体步骤如下：

S3.1，获取步骤S2.1高程数组Array[height][width]，像元空值noData和边界像元面积字典Dictionary；

构建三角网：提取每个像元的中心点作为顶点，每个顶点与相邻的8个像元中心点中的2个形成三角形；

计算边界像元面积：每个边界像元按照对角线和对边中位线分成8份，每一份面积为S₀/8，计算在构建三角网之后剩余的份数，然后计算剩余份数的平面面积，即得到该边界像元的面积；

设当前行i＝0，当前列j＝0，栅格数据的地表表面面积为S_all；

S3.2，比较i和height的值：

如果i＝height，则进入步骤S3.8；

如果i<height，则进入下一步；

S3.3，比较j和width的值：

如果j＝width，i＝i+1，j＝0，返回步骤S3.2；

如果j<width，则进入下一步；

S3.4，比较Array[i][j]和noData的值：

如果Array[i][j]＝noData，j＝j+1，返回步骤S3.3；

如果Array[i][j]>noData，则进入下一步；

S3.5，判断第i-1行数据是否存在：

如果i-1<0，则进入步骤S3.6；

如果i-1>＝0，则进入下一步；

S3.5.1，判断像元Array[i-1][j+1]、Array[i][j+1]和Array[i-1][j]：

如果Array[i-1][j+1]和Array[i][j+1]存在且不等于noData，而且Array[i-1][j]等于noData，则以索引为[i,j]、[i,j+1]和[i-1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.5.3；

S3.5.2判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i，j+1]，则Dictionary[i,j+1]＝Dictionary[i,j+1]+2*S₀/8；

如果Dictionary中存在像元[i-1，j+1]，则Dictionary[i-1，j+1]＝Dictionary[i-1，j+1]+S₀/8；

S3.5.3，判断像元Array[i-1][j]、Array[i-1][j+1]和Array[i][j+1]：

如果Array[i-1][j]和Array[i-1][j+1]存在且不等于noData，而且Array[i][j+1]不存在或者等于noData，则以索引为[i,j]、[i-1,j]和[i-1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.6

S3.5.4判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i-1,j]，则Dictionary[i-1,j]＝Dictionary[i-1,j]+2*S₀/8；

如果Dictionary中存在像元[i-1,j+1]，则Dictionary[i-1,j+1]＝Dictionary[i-1，j+1]+S₀/8；

S3.6，判断第i+1行数据是否存在：

如果i+1>＝height，则退出计算；

如果i+1<height，则进入下一步；

S3.6.1，判断像元Array[i][j+1]和Array[i+1][j+1]：

如果Array[i][j+1]和Array[i+1][j+1]存在且不等于noData，则以索引为[i,j]、[i,j+1]和[i+1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.6.3；

S3.6.2，判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i,j+1]，则Dictionary[i,j+1]＝Dictionary[i,j+1]+2*S₀/8；

如果Dictionary中存在像元[i+1,j+1]，则Dictionary[i+1,j+1]＝Dictionary[i+1,j+1]+S₀/8；

S3.6.3，判断像元Array[i+1][j]和Array[i+1][j+1]：

如果Array[i+1][j]和Array[i+1][j+1]存在且不等于noData，则以索引为[i,j]、[i+1,j]和[i+1,j+1]的像元的中心点构建三角网，每个顶点的高程即对应像元的高程值，计算三角形面积为S_tra，S_all＝S_all+S_tra，则进入下一步；

否则，进入步骤S3.7；

S3.6.4，判断该三角形三个顶点是否为边缘像元，如果是，则计算其面积：

如果Dictionary中存在像元[i，j]，则Dictionary[i,j]＝Dictionary[i,j]+S₀/8；

如果Dictionary中存在像元[i+1,j]，则Dictionary[i+1,j]＝Dictionary[i+1,j]+2*S₀/8；

如果Dictionary中存在像元[i+1,j+1]，则Dictionary[i+1,j+1]＝Dictionary[i+1,j+1]+S₀/8；

S3.7，j＝j+1,返回步骤S3.3；

S3.8，对边界像元面积字典Dictionary中所有项求和得到S_dic，并加入到S_all，即得到整个栅格数据的地表表面面积S_all＝S_all+S_dic。

此外，本发明的应用范围不局限于说明书中描述的特定实施例的工艺、机构、制造、物质组成、手段、方法及步骤。从本发明的公开内容，作为本领域的普通技术人员将容易地理解，对于目前已存在或者以后即将开发出的工艺、机构、制造、物质组成、手段、方法或步骤，其中它们执行与本发明描述的对应实施例大体相同的功能或者获得大体相同的结果，依照本发明可以对它们进行应用。因此，本发明所附权利要求旨在将这些工艺、机构、制造、物质组成、手段、方法或步骤包含在其保护范围内。