CN104574449B - 基于dem的投影面积计算方法 - Google Patents
基于dem的投影面积计算方法 Download PDFInfo
- Publication number
- CN104574449B CN104574449B CN201510043819.0A CN201510043819A CN104574449B CN 104574449 B CN104574449 B CN 104574449B CN 201510043819 A CN201510043819 A CN 201510043819A CN 104574449 B CN104574449 B CN 104574449B
- Authority
- CN
- China
- Prior art keywords
- intersection point
- area
- grid
- dem
- point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Geometry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Analysis (AREA)
Abstract
本发明提供了一种基于DEM的投影面积计算方法,步骤为:提出DEM数据及封闭区域边界,计算交点坐标;交点数据预处理:在得到交点坐标基础上,需要对交点数据预处理,判断交点的排列顺序,计算交点所在的行列数;边界区域的有效面积计算,按区间依次计算,中间整格网面积计算按照交点在每一行面积,与传统的技术相比,计算简单,且通过计算任意多边形的投影面反推该图斑的实际地表面积,该方法大大缩短计算时间,提高计算效率,且计算精度与构TIN方法相当。
Description
技术领域
本发明涉及土地测量技术领域,特别涉及一种基于DEM的投影面积计算方法。
背景技术
随着计算机的普及与GIS的迅速发展,图斑面积成为GIS中的一项重要统计内容。2013年国务院颁布了第一次地理国情普查的通知,其中一项重要内容就是土地面积的统计计算。在面积计算过程中,必须考虑计算机的内存问题。一幅DEM数据文件包括大量格网数据,若将数据全部读入将造成计算机空间的浪费;常用的面积计算方法构造三角网,该计算方法模型严密,计算精度高,但计算复杂、计算工作量大。而基于DEM的投影面积计算方法根据投影面上图斑面积,利用实际面积与投影图斑面积的比例因子,反算图斑时间面积,从而大大降低计算量,提高计算效率。
发明内容
针对现有技术中的上述缺陷和问题,本发明实施例的目的是提供一种更好的基于DEM的投影面积计算方法,大大缩短计算时间,且计算精度与构TIN方法相当,可以大大提高计算效率。
为了达到上述目的,本发明实施例提供如下技术方案:
一种基于DEM的投影面积计算方法,其特征在于,步骤如下:
Q1:提取计算区域的DEM数据及封闭区域边界,并计算交点坐标;
Q2:交点数据预处理:在得到交点数据基础上,需要对交点数据预处理,判断交点的排列顺序,计算交点所在的行列数;
Q21:交点的排列顺序,计算区域中最北、最东、最南、最西点坐标,统计这四个点在交点文件中的排列位置,然后根据排列位置判断交点的排列顺序,顺序为顺时针或逆时针。
Q22:交点的行列数计算:采用逐行计算区域面积,需要计算交点在整幅图中的行列号,统计每一个格网内的交点个数;每个交点统计的信息包括:交点在交点文件所属于的区域;交点在DEM中的行号;交点在DEM所在的列号;交点在交点文件中的行号;交点在本行交点序列中的位置;交点的属性,0为与格网的交点,1为格网内部点;
Q3:边界区域的有效面积计算,按区间依次计算,
Q31:当多条折线有公共相交部分时,需要判断格网的几个顶点属于多边形,并将格网顶点与交点按一定重新顺序排列,交点的排列顺序为顺时针时,符合左手法则,四指从线段起点指向终点,大拇指垂直于其它四指,大拇指指向为区域内方向);反之,符合右手法则;计算公式为:
(n为封闭区域交点总个数,i为n个交点中第i个交点)
Q32:当多条折线没有公共相交部分时,将当前格网分割成面积更小的格网,仍采用步骤Q31中的计算公式来计算格网中的不规则多边形的面积;
Q4:中间整格网面积计算:
按照Q22步骤中得到的交点的计算在每一行面积,
Q41:先判断属于区域有效整格网的位置;
Q42:对于每一个整格网,用中心点为顶点的四个三角形的面积代替四边形面积,三角形面积计算采用海伦公式
其中,a、b、c为三角形的边长;SΔ为三角形面积;S四边形为四边形面积。
作为上述技术方案的优选,所述Q21步骤中对排列位置判断交点的排列顺序,顺序为顺时针或逆时针,当若最北、最东、最南、最西点分别为1,20,50,88,交点排列顺序为顺时针;当最北、最东、最南、最西点分别为20,1,88,50,交点排列顺序为逆时针。
作为上述技术方案的优选,所述Q22步骤中还需要根据交点所在格网的连续性将交点分成若干区间,当第i行包含N个连续区间,则将该行交点分为N个区间,N为大于1的自然数。
作为上述技术方案的优选,所述Q3步骤应用于地表面积计算时,当对于高分辨的格网,且地形变化幅度不大时,用多边形的投影面积推算其地表面积,计算公式与步骤Q31中的计算公式相同,其投影面积与地表面积之间的函数关系为:其中,S地表为地表面积,S为使用步骤Q31中的计算公式得出的多边形的投影面积,SABCD,SA1B1C1D1为在高分辨率的格网内,投影面积与地表面积分别所在的四边形的面积。
作为上述技术方案的优选,所述Q41步骤判断属于区域有效整格网的位置的具体方法为,利用点与多边形的所属关系即点在多边形内或外判断方法,然后再利用计算射线与多边形的交点个数,判断点与多边形的所属关系。
本发明实施例提供的一种基于DEM的投影面积计算方法,与传统的技术相比,通过计算任意多边形的投影面反推该图斑的实际地表面积,该方法大大缩短计算时间,提高计算效率,且计算精度与构TIN方法相当。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例1的基于DEM的投影面积计算方法的流程步骤示意图。
图2为本发明实施例1的基于DEM的投影面积计算方法的被拆分成网格的移动窗口的结构示意图。
图3为本发明实施例1的基于DEM的投影面积计算方法的Q21步骤的交点逆时针排列的结构示意图。
图4为本发明实施例1的基于DEM的投影面积计算方法的Q21步骤的交点顺时针排列的结构示意图。
图5为本发明实施例1的基于DEM的投影面积计算方法的Q31步骤的多边形顶点逆时针排列的结构示意图。
图6为本发明实施例1的基于DEM的投影面积计算方法的Q31步骤的多边形顶点顺时针排列的结构示意图。
图7为本发明实施例1的基于DEM的投影面积计算方法的Q32步骤的当多条折线没有公共相交部分时,多条折线组成的多边形的结构示意图。
图8为本发明实施例1的基于DEM的投影面积计算方法的Q3步骤应用于地表面积计算时,多边形投影面积与地表面积的结构示意图。
具体实施方式
下面将结合本发明的附图,对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例1
如图1-8所示,本发明实施例提供的一种基于DEM的投影面积计算方法,如下:
交点数据预处理
在得到根据DEM文件内插得到的交点数据基础上,需要对交点数据预处理,判断交点的排列顺序,计算交点所在的行列数。
一)交点的排列顺序
计算区域中最北、最东、最南、最西点坐标,统计这四个点在交点文件中的排列位置,然后根据排列位置判断交点的排列顺序,顺时针或逆时针,如图3~4所示,,如若最北、最东、最南、最西点分别为1,20,50,88,交点排列顺序为顺时针;若最北、最东、最南、最西点分别为20,1,88,50,交点排列顺序为逆时针。
二)交点的行列数计算
若采用逐行计算区域面积,需要计算交点在整幅图中的行列号,统计每一个格网内的交点个数。交点数据预处理文件格式为
其中,第1列为某交点在交点文件所属于的区域;第2列为该交点在DEM中的行号;第3列为该交点在DEM所在的列号;第4列为该交点在交点文件中的行号;第5列为该交点在本行交点序列中的位置;第6列为该交点的属性,0为与格网的交点,1为格网内部点。
另外,需要根据交点所在格网的连续性将交点分成若干区间。如图2所示,第i行包含两个连续区间,因此,将该行交点分为两个区间。
边界区域的有效面积计算
边界区域的有效面积计算是面积计算的核心,同时也是面积计算的难点。在计算边界区域有效面积时,按区间依次计算。以第i行第j交点区间第k格网为例,说明边界区域有效面积计算。
一)多条折线有公共相交部分
若第k格网中有n个交点,每一个交点首尾相连,如图5~6所示。在计算不规则多边形面积时,采用
不规则多边形计算公式要求顶点坐标按顺序排列,因此,需要判断格网的几个顶点属于多边形,并将格网顶点与交点按一定重新顺序排列。
为了提高排序的可靠性,在第k格网中的交点为已排好序列的基础上,选取一个与格网相交的交点i和属于多边形的格网顶点进行排序,这样大大降低了多边形顶点重新排序造成的误判率。如图5、6中,将折线由4个交点构成的折线看作为端点为i与i+3线段,按一定法则判断格网4个顶点与线段的关系,最终判定一个顶点属于多边形。线段与格网顶点的关系满足一定的规律,但这种判断需先判断交点的排列顺序,即交点的排列顺序为顺时针时,符合左手法则(四指从线段起点指向终点,大拇指垂直于其它四指,大拇指指向为区域内方向);反之,符合右手法则。
二)多条折线没有公共相交部分
对于由3条以上折线没有公共相交部分时,采用以上方法较确定有效区域。如图7所示,三条折线没有公共相交部分,格网内有效面积由两个多变形构成,这种图形结构采用以上方法容易造成格网顶点的误判,导致面积计算错误。因此,提出了基于更高分辨率的计算方法,将格网分成更小分辨率的格网,仍采用公式(1)计算格网中的不规则多边形的面积,这样就避免了一个格网有两个以上多边形的情况。
三)地表面积计算
对于高分辨的格网,若地形变化幅度不大,可以用多边形的投影面积推算其地表面积,但如何确定二者可靠合理的函数关系是比较复杂的。在多次比较分析的基础上,提出了一种基于比例因子的新面积计算方法。
如图8所示,在高分辨率的格网内,四边形A1B1C1D1与四边形ABCD的面积比与不规则多边形S地表与S的面积比一致,即
若采用公式(1)得到不规则多边形S,则S地表为
中间整格网面积计算
在交点数据预处理中,已经交点进行区间划分。在每一行面积计算时,要判断属于区域有效整格网的位置。此处采用根据一种组合方法判断整格网的在该行位置,该方法利用点与线段的关系判断方法(线段与格网顶点判断原理一致),然后再利用陈瑞卿在西安交通大学学报的41卷第1期提出的一种快速算法进一步判断。另外,对于每一个整格网,用中心点为顶点的四个三角形的面积代替四边形面积。三角形面积计算采用海伦公式
实际应用
以陕西省吴堡县数据为例,分析该面积量算方法的计算效果。该数据包括26704个数据文件,约1300万交点个数,最终构成了17710个图斑,共50多km2。
采用新的面积计算方法与构TIN计算方法两种方案,26704个数据文件(17710个图斑)所用计算时间相差较大,其中新的面积计算方法用了3077.329秒,构TIN算法共用了12600秒;总面积互差如下表所示。
两方案总面积互差(m2) | 万分比 |
31392.5045 | 0.62 |
新的面积计算方法根据图斑在投影面上的面积,利用自适应因子推算图斑的地表面积,其精度与严密构TIN方法得到的面积精度相当,计算效率大大提高。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应所述以权利要求的保护范围为准。
Claims (5)
1.一种基于DEM的投影面积计算方法,其特征在于,步骤如下:
Q1:提取计算区域的DEM数据及封闭区域边界,并计算交点坐标;
Q2:交点数据预处理:在得到交点数据基础上,需要对交点数据预处理,判断交点的排列顺序,计算交点所在的行列数;
Q21:交点的排列顺序,计算区域中最北、最东、最南、最西点坐标,统计这四个点在交点文件中的排列位置,然后根据排列位置判断交点的排列顺序,顺序为顺时针或逆时针;
Q22:交点的行列数计算:采用逐行计算区域面积,需要计算交点在整幅图中的行列号,统计每一个格网内的交点个数;每个交点统计的信息包括:交点在交点文件所属于的区域;交点在DEM中的行号;交点在DEM所在的列号;交点在交点文件中的行号;交点在本行交点序列中的位置;交点的属性,0为与格网的交点,1为格网内部点;
Q3:边界区域的有效面积计算,按区间依次计算,
Q31:当多条折线有公共相交部分时,需要判断格网的几个顶点属于多边形,并将格网顶点与交点按一定重新顺序排列,交点的排列顺序为顺时针时,符合左手法则,四指从线段起点指向终点,大拇指垂直于其它四指,大拇指指向为区域内方向-;反之,符合右手法则;计算公式为:
n为封闭区域交点总个数,i为n个交点中第i个交点,S为投影面积;
Q32:当多条折线没有公共相交部分时,将当前格网分割成面积更小的格网,仍采用步骤Q31中的计算公式来计算格网中的不规则多边形的面积;
Q4:中间整格网面积计算:
按照Q22步骤中得到的交点的计算在每一行面积,
Q41:先判断属于区域有效整格网的位置;
Q42:对于每一个整格网,用中心点为顶点的四个三角形的面积代替四边形面积,三角形面积计算采用海伦公式
其中,a、b、c为三角形的边长;SΔ为三角形面积;S四边形为四边形面积。
2.根据权利要求1所述的一种基于DEM的投影面积计算方法,其特征在于,所述Q21步骤中对排列位置判断交点的排列顺序,顺序为顺时针或逆时针,当若最北、最东、最南、最西点分别为1,20,50,88,交点排列顺序为顺时针;当最北、最东、最南、最西点分别为20,1,88,50,交点排列顺序为逆时针。
3.根据权利要求1所述的一种基于DEM的投影面积计算方法,其特征在于,所述Q22步骤中还需要根据交点所在格网的连续性将交点分成若干区间,当第i行包含N个连续区间,则将该行交点分为N个区间,N为大于1的自然数。
4.根据权利要求1所述的一种基于DEM的投影面积计算方法,其特征在于,所述Q3步骤应用于地表面积计算时,当对于高分辨的格网,且地形变化幅度不大时,用多边形的投影面积推算其地表面积,计算公式与步骤Q31中的计算公式相同,其投影面积与地表面积之间的函数关系为: 其中,S地表为地表面积,S为使用步骤Q31中的计算公式得出的多边形的投影面积,SABCD,SA1B1C1D1为在高分辨率的格网内,投影面积与地表面积分别所在的四边形的面积。
5.根据权利要求1所述的一种基于DEM的投影面积计算方法,其特征在于,所 述Q41步骤判断属于区域有效整格网的位置的具体方法为,利用点与多边形的所属关系即点在多边形内或外判断方法,然后再利用计算射线与多边形的交点个数,判断点与多边形的所属关系。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510043819.0A CN104574449B (zh) | 2015-01-27 | 2015-01-27 | 基于dem的投影面积计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510043819.0A CN104574449B (zh) | 2015-01-27 | 2015-01-27 | 基于dem的投影面积计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104574449A CN104574449A (zh) | 2015-04-29 |
CN104574449B true CN104574449B (zh) | 2017-11-21 |
Family
ID=53090416
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510043819.0A Active CN104574449B (zh) | 2015-01-27 | 2015-01-27 | 基于dem的投影面积计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104574449B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106127820A (zh) * | 2016-06-30 | 2016-11-16 | 北京合众思壮科技股份有限公司 | 一种不规则地形的面积确定方法及装置 |
CN107239664B (zh) * | 2017-06-08 | 2021-03-26 | 国家基础地理信息中心 | 一种地表覆盖表面积栅格化计算方法及装置 |
CN108648275A (zh) * | 2018-05-10 | 2018-10-12 | 四川省冶地工程勘察设计有限公司 | 城市变化立体遥感影像匹配可疑区域自动化检验方法 |
CN110378951A (zh) * | 2018-09-10 | 2019-10-25 | 天津京东深拓机器人科技有限公司 | 确定堆垛物品支撑面积的方法和装置 |
CN110162650B (zh) * | 2019-04-18 | 2020-10-27 | 中国测绘科学研究院 | 一种兼顾局部最优与整体面积平衡的小图斑融解方法 |
CN114155350B (zh) * | 2022-02-07 | 2022-04-22 | 北京恒济引航科技股份有限公司 | 一种航道的水下地形冲淤演变分析方法及装置 |
CN116612249B (zh) * | 2023-05-22 | 2024-04-19 | 浪潮智慧科技有限公司 | 一种线段与dem相交点的分析方法、设备及介质 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102074050A (zh) * | 2011-03-01 | 2011-05-25 | 哈尔滨工程大学 | 大规模地形绘制的分形多分辨率简化方法 |
CN103077400A (zh) * | 2012-12-26 | 2013-05-01 | 中国土地勘测规划院 | 土地利用数据库支持的地类信息遥感自动识别方法 |
CN103871102A (zh) * | 2014-03-28 | 2014-06-18 | 南京大学 | 一种基于高程点和道路轮廓面的道路三维精细建模方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7310606B2 (en) * | 2006-05-12 | 2007-12-18 | Harris Corporation | Method and system for generating an image-textured digital surface model (DSM) for a geographical area of interest |
-
2015
- 2015-01-27 CN CN201510043819.0A patent/CN104574449B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102074050A (zh) * | 2011-03-01 | 2011-05-25 | 哈尔滨工程大学 | 大规模地形绘制的分形多分辨率简化方法 |
CN103077400A (zh) * | 2012-12-26 | 2013-05-01 | 中国土地勘测规划院 | 土地利用数据库支持的地类信息遥感自动识别方法 |
CN103871102A (zh) * | 2014-03-28 | 2014-06-18 | 南京大学 | 一种基于高程点和道路轮廓面的道路三维精细建模方法 |
Non-Patent Citations (3)
Title |
---|
DEM在土地勘测定界中的应用;张述清 等;《地矿测绘》;20101231;第26卷(第1期);全文 * |
基于SRTM DEM数据的全国地表面积计算研究;张伟 等;《地理与地理信息科学》;20140531;第30卷(第3期);全文 * |
基于地表粗糙度地学意义的地表真实面积计算;曾珍 等;《遥感技术与应用》;20141031;第29卷(第5期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN104574449A (zh) | 2015-04-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104574449B (zh) | 基于dem的投影面积计算方法 | |
CN103871102B (zh) | 一种基于高程点和道路轮廓面的道路三维精细建模方法 | |
CN104298734B (zh) | 土地利用现状变更调查中变更图斑移动更新方法 | |
CN108573279A (zh) | 图像标注方法及终端设备 | |
CN103336783B (zh) | 联合泰森多边形与反距离加权的密度图制图方法 | |
CN105844709A (zh) | 复杂河道地形流域洪水演进虚拟仿真的淹没线追踪方法 | |
CN102853803B (zh) | 一种文物破损面积的测定方法 | |
CN104899329B (zh) | 一种采用最小外矩形框进行宗地四至查找的方法 | |
CN102902844A (zh) | 基于大数据量dem数据的子流域划分方法 | |
CN103679807A (zh) | 一种带边界约束的散乱点云重构方法 | |
CN109636870A (zh) | 一种狭长图斑分块融解方法及装置 | |
CN108460832B (zh) | 一种基于建筑信息模型的外壳提取方法 | |
CN108917769A (zh) | 一种基于九叉树的机器人自适应栅格地图创建方法 | |
CN107918957A (zh) | 一种保持结构和纹理特征的三维建筑模型化简方法 | |
CN109949692A (zh) | 路网匹配方法、装置、计算机设备及存储介质 | |
CN110135180A (zh) | 满足节点差分隐私的度分布直方图发布方法 | |
CN103268342A (zh) | 基于cuda的dem动态可视化加速***和方法 | |
CN104731885B (zh) | 一种基于层次‑语义的多尺度空间数据拓扑关系维护方法 | |
CN107452059A (zh) | 一种三维插值建模方法、设备及其存储设备 | |
CN109858402A (zh) | 一种图像检测方法、装置、终端以及存储介质 | |
CN106373192B (zh) | 一种非拓扑一致性三维网格块体追踪算法 | |
CN107506485A (zh) | 一种基于gis的离散站点分级显示的方法及装置 | |
CN104331389A (zh) | 基于八点法的等值线追踪算法 | |
CN112116709A (zh) | 一种提高地形表达精度的地形特征线处理方法 | |
CN106097446A (zh) | 一种基于大规模地震数据的三维地质曲面重建方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |