CN104821013A - Method and system for specific surface area extraction based on geodetic coordinate system digital elevation model - Google Patents
Method and system for specific surface area extraction based on geodetic coordinate system digital elevation model Download PDFInfo
- Publication number
- CN104821013A CN104821013A CN201510235733.8A CN201510235733A CN104821013A CN 104821013 A CN104821013 A CN 104821013A CN 201510235733 A CN201510235733 A CN 201510235733A CN 104821013 A CN104821013 A CN 104821013A
- Authority
- CN
- China
- Prior art keywords
- raster
- cell
- surface area
- point
- grid
- 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.)
- Pending
Links
Landscapes
- Processing Or Creating Images (AREA)
Abstract
The invention provides a method and a system for specific surface area extraction based on a geodetic coordinate system digital elevation model (DEM). The method comprises: an original DEM using the geodetic coordinate system to store, creating a grid file; traversing grid pixels of the original DEM in sequence, selecting adjacent four points to establish a regular triangle, and calculating surface area of a new grid pixel region, defining the obtained grid file as a digital specific surface area model SAM; giving any vector polygon, performing spatial overlaying on the vector polygon and the SAM, and summarizing pixel values of the grid pixels which are in the vector polygon in the SAM, to obtain the specific surface area of the vector polygon. As a whole, the method and the system are characterized by simple realization and high efficiency. Compared with the prior art, the precision is greatly improved.
Description
Technical field
The invention belongs to geographical information space analysis technical field, particularly relate to a kind of earth's surface surface area extractive technique scheme based on earth coordinates digital elevation model.
Background technology
In geographical information space analysis technical field, as follows to the correlative study briefing of earth's surface surface area at present:
(1) earth's surface surface area
Earth's surface surface area refers to the area of the earth surface of Combining with terrain condition.In technical field of geographic information, usually the earth is reduced to the spheroid of a standard when modeling is carried out to the earth, and does not consider the fluctuating situation of earth surface reality.Therefore, when carrying out map making and area measuring, all kinds of atural object on earth's surface and key element are conceptualized as the geometric object of the different shapes such as the point, line, surface be distributed on spheroid usually.When carrying out the area statistics of atural object, statistics is the ellipsoid area of such atural object on standard spheroid body usually, but not the real surface area of atural object.Such as, in the Second National land investigation project of Ministry of Land and Resources, the area statistics of all land use patterns such as arable land, field, forest land are all ellipsoid areas, instead of land use pattern real surface area at the earth's surface.
Due to earth surface highland big rise and fall, cause earth's surface surface area and ellipsoid area can there is larger difference in mountain area and knob.Such as, can knob in the larger mountain area of topographic relief, earth's surface surface area differs can reach more than 20% with ellipsoid area.Ellipsoid area is the area of atural object on this theoretical model of earth ellipsoid, and surface area is only earth object true area at the earth's surface.Therefore, only represent that the area of atural object is the needs that can not meet the practical applications such as resource exploration (territory investigation, Forestry Investigation) with the ellipsoid area of atural object.As can be seen here, utilize Local topographic relief data, determine that all kinds of atural object and key element surface area is at the earth's surface very necessary.Therefore, presided over by State Council's first time geographical national conditions census office, the whole nation first time geographical national conditions of carrying out in 2013-2015 years generally investigate with Basic Statistics and Tables, the formal earth's surface surface area using all kinds of for earth's surface key element calculate as one of formal geographical national conditions statistical study index (list of references 1: State Council's first time office of national geographical national conditions generaI investigation leading group. geographical national conditions generaI investigation basic statistics technical stipulation (GDPJ02-2013)).
(2) based on the earth's surface surface area extracting method of digital elevation model
At area of geographic information, extracting method at present for earth's surface surface area mainly contains following several: 1) utilize zone digit elevation model (DEM) to build the TIN (TIN) in region, superposition needs the vector polygon of gauging surface area on this basis, obtains each polygonal surface area by TIN and polygonal superposition calculation, (list of references 2: Xie Chenglei, Zhao Rong, Liang Yong. based on geographical national conditions statistic unit surface area accurate Calculation [J] of geographic coordinate. sensor information, 2014, 04:47-51.) 2) directly regular triangular net is built based on Law of DEM Data, superposition needs the vector polygon of gauging surface area on this basis, by regular triangular net and polygonal spatial relation, the area gathering the spatial triangle fallen in polygon obtains each polygonal surface area (list of references 3: State Council's first time national geographical national conditions generaI investigation leading group office .GDPJ 02-2013 geographical national conditions generaI investigation basic statistics technical stipulation, 2013.9.), 3) based on adjacent 4 DEM grid cell, raster cells, integral formula is utilized to extrapolate the surface area in the region be made up of 4 adjacent picture elements, and then according to this curved surface and polygonal spatial relation, gather the surface area fallen in polygon, thus obtain polygonal surface area (list of references 4: Jiang Fan, Lv Xiaohua, Wang Zhonglan. based on the DEM surface area Algorithm Analysis [J] of Composed Calculus Algorithms. Institute of Surveying and Mapping's journal, 2005,04:263-265.).4) DEM upper 3 is utilized to take advantage of the pixel window of 3 to build 8 spatial triangles, calculate the surface area in the pixel region, center in 3*3 region, and then utilize grid cell, raster cell to superpose with polygon, gather and to fall in polygon grid cell, raster cell value as polygonal surface area.(list of references 5: Xie Chenglei, Zhao Rong, Liang Yong. based on geographical national conditions statistic unit surface area accurate Calculation [J] of geographic coordinate. sensor information, 2014,04:47-51.; List of references 6:Jenness, JS.Calculating landscape surface area from digital elevation models, WILDLIFE SOCIETYBULLETIN, 2004,03,829-839)
In above-mentioned four kinds of methods, method 1 is the method be most widely used at present, also be the method used in Geographic Information System (GIS) the software ArcGIS of international popular, its main thought is, build TIN by region DEM, approach the truth on earth's surface with multiple irregular spatial triangle.Regional feature is more complicated, then triangle is closeer, and landform is more smooth, then triangle is more sparse.When calculating polygon surface area, TIN is superposed the calculating that can complete surface area with polygon; Its principal feature is, inessential grid cell, raster cell in the DEM that utilized TIN to cast out, can go out the main landform relief feature in region, to a certain degree decrease the calculated amount of spatial triangle area with less Based on Feature Points.But built by DEM in the process of TIN, owing to having given up so-called inessential grid cell, raster cell, the loss of terrain model precision can be caused, and then it is inaccurate to cause surface area to calculate; On the other hand, the structure of TIN needs to calculate based on projected coordinate system, namely needs ellipsoid latitude and longitude coordinates represented to be projected as X in plane in units of rice, Y-coordinate.As can be seen here, the geometric configuration after projection of the geometric object on ellipsoid also can change, and then causes the result of calculation of earth's surface surface area also likely to there is error.
Method 2 is a kind of directly based on DEM structure regular triangular net, and then calculates the method for area of a polygon.Its principal feature is: directly build based on the DEM under earth coordinates, do not have loss of significance; But, when carrying out regular triangular net and polygon is vector superposed, there is the problems such as calculated amount is larger.
Method 3 is a kind of directly based on DEM and then utilize Composed Calculus Algorithms fitting surface, and then calculates the method for area of a polygon.The problems such as its principal feature is: precision is higher, calculation of complex, and calculated amount is large, the main determination be applicable to small-scale (polygon number is few or zoning is little) earth's surface surface area.
Method 4 is also a kind of directly based on DEM and then utilize Composed Calculus Algorithms fitting surface, and then calculates the method for area of a polygon.Its principal feature is: precision is higher, calculates the advantages such as simple.Wherein, list of references 4 and list of references 5 difference to some extent.List of references 4 can directly utilize earth coordinates to calculate, and list of references 5 mainly calculates for plane right-angle coordinate; List of references 4 has carried out further process for the error produced in vector raster overlay, but the method is comparatively complicated, adds calculated amount.Thus the high precision computation demand to small-scale (polygon number is few or zoning is little) earth's surface surface area is applicable to, and to (polygon quantity reaches more than 100,000 on a large scale, or region area reaches more than 10,000 square kilometres) there is the deficiencies such as computing time is long in solving of problem, thus limits the scope of its application.
In sum, the deficiency that the present invention will exist for existing method, the advantage of comprehensive various method, invents a kind of earth's surface surface area extractive technique scheme based on earth coordinates digital elevation model, meets extensive earth's surface surface area measurement demand.
For ease of understanding technical solution of the present invention, first explain the theoretical foundation that the present invention relates to:
(1) conventional in geography information several coordinate systems
There is a position for point to describe earth surface in three dimensions accurately, usually needing to be accomplished by certain specific coordinate system.Usually, conventional coordinate system mainly can be divided into 3 kinds: 1) earth coordinates, namely the position of certain point is represented with latitude and longitude coordinates, usually with (B, L, H) mode is expressed, and wherein, the longitude that B represents the latitude of this point, L represents this point, H represent the elevation (height above sea level) of this point; In earth coordinates, the coordinate unit of B, L is degree, and the unit of H is rice; 2) plane projection coordinate system, on the basis of earth coordinates, by using certain projecting method (as equiareal projection, equal angle projection's method or equidistant projection method), based on specific central meridian or standard parallel, earth ellipsoid face is projected as plane, and the position that represents a little is carried out with planimetric rectangular coordinates, as (X, Y, Z), wherein, X, Y represent this off-set value relative to origin of projection coordinates, and Z is the height value of this point.Wherein, Z is the particular community eigenwert of certain point in plane right-angle coordinate, does not form a dimension in coordinate system.The coordinate unit of X, Y, Z is rice.3) rectangular coordinate system in space take earth center as true origin, states certain some position in rectangular coordinate system in space in the mode of (X, Y, Z).In rectangular coordinate system in space, the unit of X, Y, Z is rice.Earth reference ellipsoid parameter (major axis, excentricity etc. of earth ellipsoid) conversion and inverse mutually can be passed through between three coordinate systems.
In above-mentioned coordinate system: 1) under earth coordinates, the position of each point directly represents based on latitude and longitude coordinates, and the position of point accurately, does not offset and distortion; 2) under projected coordinate system, the position of each point is a kind of planimetric position, and according to the difference of projecting method, the geometric configuration before and after projection may all produce certain distortion and skew in angle, length and area etc.; 3) rectangular coordinate system in space is a kind of three-dimensional system of coordinate, and the position of each point is described by real X, Y, the Z coordinate under three-dimensional system of coordinate, does not have because project the change brought.In three kinds of coordinate systems, the some position that earth coordinates and rectangular coordinate system in space represent and geometric configuration less than distortion, and have distortion or skew based on the some position of plane right-angle coordinate and geometric configuration.
(2) digital elevation model
Digital elevation model (Digital Elevation Model), vehicle economy M, utilizes grating image to represent the height fluctuating situation of earth surface elevation (also referred to as " height above sea level ") usually.Grating image stores based on specific geographical co-ordinate system usually, and each grid is according to the specific region of certain block of its coordinate number corresponding earth surface in geographical co-ordinate system.The height value in this region of numeric representation of grid.The basic thought of digital elevation model approaches earth's surface real landform height fluctuating situation by using grid that is discrete, rule.In DEM, the resolution (size of grid) of grid cell, raster cell determines the precision of DEM: raster resolution higher (grid is less) usually, then its precision is higher, more can approach earth's surface true form really; Raster resolution lower (grid is larger), then its precision is lower.Raster resolution is higher, then the data volume of DEM is larger, and consequent calculated amount is also larger; Otherwise then data volume is less, calculated amount is less.When carrying out the spatial analysis based on DEM, usually needing the accuracy requirement according to result of calculation, under the prerequisite ensureing computational solution precision, selecting the DEM of adequate resolution.
(3) vector data and raster data
Vector model and Raster Data Model are two kinds of different modes of expressing geospatial object in GIS.In vector model, geospatial object is conceptualized as the geometric object of the different shapes such as polygon, line and point, and with the information of accurate coordinate storage in conjunction with node each on object; In Raster Data Model, usually represent the spatial distribution characteristic of geographic object with discrete grid cell, raster cell and property value thereof.In GIS, two kinds of models respectively have quality: vector model can the shape of point-device description geographic object, border, scope and position, but inefficiency when doing spatial operation based on vector model, calculated amount is large; Raster Data Model can not describe the shape of geographic object, scope, border and position accurately, but when carrying out spatial operation, counting yield is high, calculated amount is little.
Summary of the invention
For the limitation of existing earth's surface surface area extractive technique in computational accuracy and efficiency, the invention provides a kind of high precision, high efficiency earth's surface surface area extracting method, improve the efficiency extracted towards extensive earth's surface surface area.
Technical scheme provided by the invention comprises a kind of earth's surface surface area extracting method based on earth coordinates digital elevation model, carries out following steps,
Step one, establish original DEM adopt earth coordinates store, DEM grid cell, raster cell line number is M, and columns is N, and raster resolution is C, and described DEM is digital elevation model; Create (M-1) × 2 row, the raster file of sky that (N-1) × 2 arrange, in raster file, the value of each grid cell, raster cell is defaulted as 0, and the coordinate system of raster file is identical with original DEM; If the starting point latitude and longitude coordinates of original DEM is (B1, L1), then the starting point latitude and longitude coordinates of raster file is (B1+C, L1+C);
Step 2, the grid cell, raster cell of original DEM to be traveled through successively, builds regular triangle and to calculate the surface area in new grid cell, raster cell region as follows, carry out following sub-step,
Step 2.1, if the row number that i represents the line number of current traversal, j represents current traversal, then the grid cell, raster cell ranks coordinate of current traversal is (i, j), and making i initial value be 2, j initial value is 2,
Step 2.2, determines current grid pixel point according to current (i, j), and obtains ranks number for (i-1, j-1), and (i, j-1) and (i-1, j) grid cell, raster cell point, forms 4 adjacent points;
Step 2.3,4 points adjacent according to step 2.2 gained obtain a quadrilateral in the plane, and interpolation goes out the mid point height value of quadrilateral four edges and the height value of diagonal line intersection point, obtains 5 points;
Step 2.4, in 4 points adjacent by step 2.2 gained and step 2.3 gained 4 corresponding quadrilaterals of point, the mid point of four edges and diagonal line intersection point, obtain 8 spatially adjacent triangles;
Step 2.5, to the latitude and longitude coordinates of step 2.2 and step 2.3 gained totally 9 points, is converted to rectangular space coordinate according to following formula,
X＝(N+H)×Cos(B)×Cos(L)
Y＝(N+H)×Cos(B)×Sin(L)
Z＝[N×(1-e
2)+H]×Sin(B)
In formula, B, L, H represent latitude, longitude and height value of certain point respectively, and X, Y, Z represent this coordinate figure in three dimensions rectangular coordinate system respectively; The cosine value that Cos (B) is latitude, Sin (B) is the sine value of latitude, and Sin (L) is the sine value of longitude, and Cos (L) is the cosine value of longitude; e
2for earth's spheroid the first excentricity square, the computing formula of Parameter N is as follows,
In formula, a is the major semi-axis of earth's spheroid;
Step 2.6, utilizes Heron's formula to calculate 8 leg-of-mutton areas respectively, comprises and obtains area according to the rectangular space coordinate on an Atria summit;
Step 2.7, obtains quadrilateral by adjacent triangle combination of two, obtains 4 new quadrilaterals altogether;
Step 2.8, according to the component relationship of 8 triangles and 4 new quadrilaterals, gathers the area that leg-of-mutton area obtains 4 new quadrilaterals respectively;
Step 2.9, is written in the corresponding ranks number (2i-3,2j-3) of new raster file, (2i-2,2j-3), (2i-2,2j-2), (2i-3,2j-2) by the area of 4 new quadrilaterals; The ranks coordinate getting the next one traversal of original DEM is new (i, j), and return step 2.2, until traveled through, the raster file obtained is defined as digitally surface area model, is labeled as SAM, enters step 3;
Step 3, given any vector polygon, carry out Spatial Overlap by vector polygon and SAM, gathers in SAM the pixel value dropping on grid cell, raster cell in vector polygon and obtain V, using the surface area of V as this vector polygon; Wherein, whether grid cell, raster cell drops in vector polygon judges according to following rule,
If the central point of grid cell, raster cell drops in vector polygon or on vector polygon limit, then this grid cell, raster cell is calculated in vector polygon, otherwise does not exist.
The present invention also provides a kind of earth's surface surface area extraction system based on earth coordinates digital elevation model, comprises with lower unit,
Raster file sets up unit, earth coordinates are adopted to store for establishing original DEM, DEM grid cell, raster cell line number is M, columns is N, raster resolution is C, and described DEM is digital elevation model, creates (M-1) × 2 row, the raster file of sky that (N-1) × 2 arrange, in raster file, the value of each grid cell, raster cell is defaulted as 0, and the coordinate system of raster file is identical with original DEM; If the starting point latitude and longitude coordinates of original DEM is (B1, L1), then the starting point latitude and longitude coordinates of raster file is (B1+C, L1+C);
Grid cell, raster cell region surface area extraction unit, for traveling through successively the grid cell, raster cell of original DEM, building regular triangle and calculating the surface area in new grid cell, raster cell region, comprising with lower module,
Initialization module, for establishing the row number that i represents the line number of current traversal, j represents current traversal, then the grid cell, raster cell ranks coordinate of current traversal is (i, j), and making i initial value be 2, j initial value is 2,
Consecutive point extraction module, for determining current grid pixel point according to current (i, j), and obtain ranks number for (i-1, j-1), (i, j-1) and (i-1, j) grid cell, raster cell point, forms 4 adjacent points;
Interpolating module, obtain a quadrilateral in the plane for 4 points adjacent according to consecutive point extraction module gained, interpolation goes out the mid point height value of quadrilateral four edges and the height value of diagonal line intersection point, obtains 5 points;
Triangular modeling block, for by the mid point of four edges in adjacent 4 points of consecutive point extraction module gained and interpolating module gained 4 corresponding quadrilaterals of point and diagonal line intersection point, obtains 8 spatially adjacent triangles;
Rectangular space coordinate modular converter, for the latitude and longitude coordinates to consecutive point extraction module and interpolating module gained 9 points, is converted to rectangular space coordinate according to following formula,
X＝(N+H)×Cos(B)×Cos(L)
Y＝(N+H)×Cos(B)×Sin(L)
Z＝[N×(1-e
2)+H]×Sin(B)
In formula, B, L, H represent latitude, longitude and height value of certain point respectively, and X, Y, Z represent this coordinate figure in three dimensions rectangular coordinate system respectively; The cosine value that Cos (B) is latitude, Sin (B) is the sine value of latitude, and Sin (L) is the sine value of longitude, and Cos (L) is the cosine value of longitude; e
2for earth's spheroid the first excentricity square, the computing formula of Parameter N is as follows,
In formula, a is the major semi-axis of earth's spheroid;
Triangle area extraction module, for utilizing Heron's formula to calculate 8 leg-of-mutton areas respectively, comprising and obtaining area according to the rectangular space coordinate on an Atria summit;
Quadrilateral builds module, for obtaining quadrilateral by adjacent triangle combination of two, obtains 4 new quadrilaterals altogether;
Quadrilateral area extraction module, for the component relationship according to 8 triangles and 4 new quadrilaterals, gathers the area that leg-of-mutton area obtains 4 new quadrilaterals respectively;
Digitally surface area model generation module, for the area of 4 new quadrilaterals being written to corresponding ranks number (2i-3,2j-3), (2i-2,2j-3), (2i-2 of new raster file, 2j-2), in (2i-3,2j-2); The ranks coordinate getting the next one traversal of original DEM is new (i, j), and order consecutive point extraction module works, until traveled through, the raster file obtained is defined as digitally surface area model, is labeled as SAM, and order surface area determination module works;
Surface area determination module, for given any vector polygon, carries out Spatial Overlap by vector polygon and SAM, gathers in SAM the pixel value dropping on grid cell, raster cell in vector polygon and obtains V, using the surface area of V as this vector polygon; Wherein, whether grid cell, raster cell drops in vector polygon judges according to following rule,
If the central point of grid cell, raster cell drops in vector polygon or on vector polygon limit, then this grid cell, raster cell is calculated in vector polygon, otherwise does not exist.
The present invention proposes a kind of earth's surface surface area extractive technique scheme based on earth coordinates digital elevation model, by terrestrial coordinate is converted to rectangular coordinate system in space, build rule space network of triangle on this basis, and then build surface area raster file by the adjacent space triangulation network; Space overlapping based on surface area raster file and vector polygon operates, and completes the extraction of the surface area to arbitrary polygon.The present invention has generally and realizes simple, efficiency high.Relative to the earth's surface surface area extracting method of conventional at present (TIN) based on TIN, precision of the present invention also improves a lot.The surface area extracting method that the present invention proposes no matter efficiency and precision can both meet extensive geometric object earth's surface surface area extraction demand in practical engineering application, automatically can determine multiple type of ground objects surface area fast, can extensively use in resource management, economic development decision support and national conditions investigation.
Accompanying drawing explanation
Fig. 1 is that in the embodiment of the present invention, grid cell, raster cell superposes schematic diagram with vector polygon;
Fig. 2 is adjacent 4 the DEM grid cell point schematic diagram of the embodiment of the present invention;
Fig. 3 is embodiment of the present invention DEM interpolation theory schematic diagram;
Fig. 4 is that the new quadrilateral of the embodiment of the present invention builds principle schematic.
Embodiment
Below in conjunction with drawings and Examples, technical solution of the present invention is specifically described.
The prerequisite calculating earth's surface surface area is carried out numerical simulation to surface, earth's surface and approaches.DEM simulates earth's surface with grid that is discrete, rule.But earth's surface is a continuous print surface, cannot directly obtain continuous print surface based on DEM.Therefore, the present invention proposes, when calculating surface, earth's surface, the interconnective triangulation network can be built (in rectangular coordinate system in space based on grid point discrete on DEM, any 3 points can ensure in same plane, and the point of more than 3 cannot ensure in same plane, also cannot reference area).The triangulation network formed with interconnective multiple triangle is to approach the real continuous surface in earth's surface.When calculating the earth's surface surface area of arbitrary polygon, by carrying out Spatial Overlap Analysis with the triangulation network of vector form, the statistics triangle area dropped in polygon can obtain polygonal surface area.Do computing to same area two VectorLayers and two raster map layers, the time needed for grid computing even only needs one of percentage of vector data computing.And when vector graphics is more, the odds for effectiveness of grid computing is more obvious.
The surface area extracting method based on earth coordinates digital elevation model that the embodiment of the present invention provides, can adopt computer software technology to realize automatically running.
Various coordinate system has different application in different application scenario.Such as when carrying out map making, usually need ellipsoid to be projected as plane right-angle coordinate to be plotted on two dimensional surface by the geographic element on ellipsoid; And when carrying out 3 D stereo and analyzing, then can directly adopt three dimensions rectangular coordinate system to avoid loss of significance.In the present invention, for ensureing not produce distortion, original DEM adopts earth coordinates to store, if the line number of grid cell, raster cell is M, columns is N, and raster resolution is C, specifically comprises the following steps:
One, the raster file of the sky of establishment (M-1) × 2 row, (N-1) × 2 row, in raster file, the value of each grid cell, raster cell is defaulted as 0.Grid cell, raster cell Value Types is double, and the coordinate system of raster file is identical with original DEM, is namely also earth coordinates; If the starting point latitude and longitude coordinates of original DEM is (B1, L1), then the starting point latitude and longitude coordinates of new raster file is (B1+C, L1+C).
Two, the grid cell, raster cell of original DEM is traveled through successively, builds regular triangle in accordance with the following methods and calculate the surface area in new grid cell, raster cell region:
(1) establish the row number that i represents the line number of current traversal, j represents current traversal, then the grid cell, raster cell ranks coordinate of current traversal is (i, j).Making i initial value be 2, j initial value is 2, namely travels through from the second row secondary series of DEM.
(2) see Fig. 2, current grid pixel point is determined according to current (i, j), and obtain ranks number for (i-1, j-1), (i, j-1) and (i-1, j) grid cell, raster cell point, form 4 adjacent points and be numbered respectively.Current grid pixel point is No. 3 points, and ranks number are (i-1, j-1), and the grid cell, raster cell point of (i, j-1) and (i-1, j) is respectively No. 1 point, No. 4 points, No. 2 points.
(3) 4 points having (2) gained adjacent can obtain a quadrilateral in the plane, and quadrilateral four edges divides other mid point and diagonal line intersection point totally 5 points, and interpolation goes out the mid point height value of quadrilateral four edges and the height value of diagonal line intersection point.See Fig. 3, No. 1 point and No. 2 points form limit mid point be No. 6 points, No. 2 points and No. 3 points form limit mid point be No. 7 points, No. 3 points and No. 4 points form limit mid point be No. 8 points, No. 1 point and No. 4 points form limit mid point be No. 9 points.The mid point on any limit, as 6,7,8, the height value of No. 9 points and latitude and longitude coordinates are obtained by the mean value calculation of height value of 2 on its limit.Such as, the height value of No. 6 points is the mean value of the height value of No. 1 point and No. 2 points; The like can obtain 7,8, the height value of No. 9 points.The height value of central point No. 5 points by 1,2,3, the mean value of the height value of No. 4 points obtains.
(4) can obtain 8 spatially adjacent spatial triangles by (2) and (3), number consecutively: I, II, III, IV, V, VI, VII, VIII, as accompanying drawing 3.
(5) by 9 points (No. 1 point, No. 2 points ... No. 9 points) latitude and longitude coordinates be converted to rectangular space coordinate according to following formula:
X＝(N+H)×Cos(B)×Cos(L)
Y＝(N+H)×Cos(B)×Sin(L)
Z＝[N×(1-e
2)+H]×Sin(B)
In formula: B, L, H represent certain latitude put, longitude and height value respectively; X, Y, Z represent this coordinate figure in three dimensions rectangular coordinate system respectively; The cosine value that Cos (B) is latitude, Sin (B) is the sine value of latitude, and Sin (L) is the sine value of longitude, and Cos (L) is the cosine value of longitude.Constant e
2for earth's spheroid the first excentricity square, the known constant as earth's spheroid inputs, such as, country 2000 the earth ellipsoids e
2for 6.6943800229008E-03.The computing formula of Parameter N is as follows:
In formula: a is the major semi-axis of earth's spheroid.A and e
2be constant, those skilled in the art can search and obtain in the definition of reference ellipsoid of earth.
(6) according to (5) acquired results, utilize Heron's formula to calculate the area of 8 spatial triangles respectively, namely obtain area according to the coordinate on spatial triangle three summits.Be implemented as follows:
For any spatial triangle, the A on given spatial triangle summit, B, C and rectangular space coordinate (X thereof
a, Y
a, Z
a), (X
b, Y
b, Z
b) and (X
c, Y
c, Z
c), its area S is drawn by following formula:
In formula: variable p=(a+b+c)/2.A, b, c represent the leg-of-mutton length of side respectively.The length of side of spatial triangle can adopt following formula to calculate:
In formula: l
ijrepresent any two points i, the three dimensions distance of j, (X
i, Y
i, Z
i) and (X
j, Y
j, Z
j) represent the rectangular space coordinate of i, j at 2 respectively.
(7) a new plane quadrilateral can be obtained by triangle I, II to be numbered (one), triangle III, IV can obtain a new plane quadrilateral and be numbered (two), triangle V, VI can obtain a new plane quadrilateral and be numbered (three), triangle VII, VIII can obtain a new plane quadrilateral and be numbered (four), as shown in Figure 4.
(8) according to the component relationship of 8 spatial triangles and 4 new quadrilaterals, gather respectively triangle (), (two), (three), (four) area obtain the area of 4 new quadrilaterals, be designated as S1 respectively, S2, S3, S4.
(9) area of 4 new quadrilaterals is written in the corresponding ranks number of new raster file, namely obtains the surface area in new grid cell, raster cell region.Wherein, S1 writes the ranks coordinate (2i-3,2j-3) of new grid, S2 writes the ranks coordinate (2i-3,2j-2) of new grid, and S3 writes the ranks coordinate (2i-2 of new grid, 2j-2), S4 writes the ranks coordinate (2i-2,2j-3) of new grid.The ranks coordinate getting the next one traversal of original DEM is new (i, j), returns (2), such repeated execution of steps (2)-step (8), travels through ranks coordinate (2,2), (2 successively, 3) ... (2, N), (3,2), (3,3) ... (3, N),, (M, 2) are to (M, N), until after having traveled through, step 3 has been entered.
From accompanying drawing 4, the spatial resolution of new raster file is original half, and spatial accuracy is higher.In the new raster file completed, the value of grid cell, raster cell illustrates this district's pixel and to take up space the earth's surface surface area in region, and this raster file is defined as digitally surface area model, is labeled as SAM.
Three, given any vector polygon, carries out Spatial Overlap by polygon and SAM, gathers in SAM the pixel value dropping on grid cell, raster cell in polygon and obtains V, using V as this polygonal surface area.
Wherein, whether grid cell, raster cell drops in polygon judges according to following rule: if the central point of grid cell, raster cell drops in polygon or on polygon limit, then this pixel is calculated in polygon, otherwise does not exist.As shown in Figure 1, a rectangle represents a grid cell, raster cell, and the round dot in rectangle represents the central point of grid cell, raster cell.When doing spatial operation based on grid, usually whole grid is replaced to improve counting yield to simplify to calculate with the central point of grid cell, raster cell.
Example is calculated as with the land-use map spot surface area in county of Hubei Province.The Land_use change vector figure data of existing this county, marked the land-use style of this figure spot in the attribute of each vector graphics, Land_use change polar plot spot number is about 14.4 ten thousand, and is covered with whole region.The ellipsoid area of known this county is 1980.85 square kilometres, and the ellipsoid area in arable land is 598 square kilometres, the ellipsoid area in field is 165 square kilometres, forest land area is 850 square kilometres.The spatial resolution of DEM is 10 meters.DEM and Land_use change polar plot spot all adopt national 2000 earth coordinates to store.
The hardware configuration of test computer is: CPU Intel Core i5-4590, dominant frequency 3.3GHZ; Internal memory 4GB, hard disk free space is greater than 100G.Computer software is: Windows 7, ArcGIS 10.2.1.
Utilize C# programming language, based on ArcGIS software platform, realize the method that the embodiment of the present invention provides.First, utilize embodiment method to build the surface area raster file in whole county, 67 seconds consuming time.On this basis, the surface area completing 14.4 ten thousand vector graphics for 25 seconds consuming time calculates.Gather the surface area of polar plot spot according to type, the surface area obtaining whole county is 2070 square kilometres, and arable land surface area is 609 square kilometres, field area 173 square kilometres, and forest land area is 912 square kilometres.In the case of this county, the total surface area of this county has had more about 90 square kilometres than ellipsoid area, and wherein arable land surface area has had more 11 square kilometres than arable land ellipsoid area, and forest land surface area has had more 62 square kilometres than forest land ellipsoid area.As can be seen here, extract various type of ground objects surface area based on high resolution DEM, be very important in resource management, economic development decision support and national conditions investigation.
Relative to the method 1 mentioned in background technology, main advantage of the present invention is: (1) directly calculates based on the DEM under earth coordinates, avoids the loss of significance needing to carry out plane projection and DEM grid point vacuate when building TIN; (2) operand that the operand utilizing vector polygon and grid surface file to carry out space overlapping carries out space overlapping than vector polygon and TIN is little.When needing the polygon quantity of gauging surface area more, the odds for effectiveness of the inventive method is more obvious.
Relative to the method 2 mentioned in background technology, main advantage of the present invention is: avoid the huge calculated amount produced when the vector polygon Sum fanction triangulation network (being also vector) superposes, economize on resources, odds for effectiveness is obvious.
Relative to the method 3 mentioned in background technology, main advantage of the present invention is: calculate simple, calculated amount is little, and without the need to the integral operation of complexity, required system resource is few, is applicable to surface area on a large scale and extracts solving of task.
Relative to the method 4 mentioned in background technology, main advantage of the present invention is: the spatial resolution generating the new raster file obtained is the half of original DEM, and spatial accuracy when vector and raster overlay computing is promoted further.For 10*10 rice resolution, according to method 4 carry out grid judge time, the numerical value of each area grid is at least 100 square metres; When judging according to the present invention, original grid is divided into four by one, and the numerical value of each grid is at least 25 square metres.Grid division more detailed rules and regulations more can describe the border of geometric object accurately.
During concrete enforcement, software modularity mode can also be adopted to provide corresponding system.The invention provides a kind of earth's surface surface area extraction system based on earth coordinates digital elevation model, comprise with lower unit,
Raster file sets up unit, earth coordinates are adopted to store for establishing original DEM, DEM grid cell, raster cell line number is M, columns is N, raster resolution is C, and described DEM is digital elevation model, creates (M-1) × 2 row, the raster file of sky that (N-1) × 2 arrange, in raster file, the value of each grid cell, raster cell is defaulted as 0, and the coordinate system of raster file is identical with original DEM; If the starting point latitude and longitude coordinates of original DEM is (B1, L1), then the starting point latitude and longitude coordinates of raster file is (B1+C, L1+C);
Grid cell, raster cell region surface area extraction unit, for traveling through successively the grid cell, raster cell of original DEM, building regular triangle and the surface area in computation grid pixel region, comprising with lower module,
Initialization module, for establishing the row number that i represents the line number of current traversal, j represents current traversal, then the grid cell, raster cell ranks coordinate of current traversal is (i, j), and making i initial value be 2, j initial value is 2,
Consecutive point extraction module, for determining current grid pixel point according to current (i, j), and obtain ranks number for (i-1, j-1), (i, j-1) and (i-1, j) grid cell, raster cell point, forms 4 adjacent points;
Interpolating module, obtain a quadrilateral in the plane for 4 points adjacent according to consecutive point extraction module gained, interpolation goes out the mid point height value of quadrilateral four edges and the height value of diagonal line intersection point, obtains 5 points;
Triangular modeling block, for by the mid point of four edges in adjacent 4 points of consecutive point extraction module gained and interpolating module gained 4 corresponding quadrilaterals of point and diagonal line intersection point, obtains 8 spatially adjacent triangles;
Rectangular space coordinate modular converter, for the latitude and longitude coordinates to consecutive point extraction module and interpolating module gained 9 points, is converted to rectangular space coordinate according to following formula,
X＝(N+H)×Cos(B)×Cos(L)
Y＝(N+H)×Cos(B)×Sin(L)
Z＝[N×(1-e
2)+H]×Sin(B)
In formula, B, L, H represent latitude, longitude and height value of certain point respectively, and X, Y, Z represent this coordinate figure in three dimensions rectangular coordinate system respectively; The cosine value that Cos (B) is latitude, Sin (B) is the sine value of latitude, and Sin (L) is the sine value of longitude, and Cos (L) is the cosine value of longitude; e
2for earth's spheroid the first excentricity square, the computing formula of Parameter N is as follows,
In formula, a is the major semi-axis of earth's spheroid;
Triangle area extraction module, for utilizing Heron's formula to calculate 8 leg-of-mutton areas respectively, comprising and obtaining area according to the rectangular space coordinate on an Atria summit;
Quadrilateral builds module, for obtaining quadrilateral by adjacent triangle combination of two, obtains 4 new quadrilaterals altogether;
Quadrilateral area extraction module, for the component relationship according to 8 triangles and 4 new quadrilaterals, gathers the area that leg-of-mutton area obtains 4 new quadrilaterals respectively;
Digitally surface area model generation module, for the area of 4 new quadrilaterals being written to corresponding ranks number (2i-3,2j-3), (2i-2,2j-3), (2i-2 of new raster file, 2j-2), in (2i-3,2j-2); The ranks coordinate getting the next one traversal of original DEM is new (i, j), and order consecutive point extraction module works, until traveled through, the raster file obtained is defined as digitally surface area model, is labeled as SAM, and order surface area determination module works;
Surface area determination module, for given any vector polygon, carries out Spatial Overlap by vector polygon and SAM, gathers in SAM the pixel value dropping on grid cell, raster cell in vector polygon and obtains V, using the surface area of V as this vector polygon; Wherein, whether grid cell, raster cell drops in vector polygon judges according to following rule,
If the central point of grid cell, raster cell drops in vector polygon or on vector polygon limit, then this grid cell, raster cell is calculated in vector polygon, otherwise does not exist.
Specific embodiment described herein is only to the explanation for example of the present invention's spirit.Those skilled in the art can make various amendment or supplement or adopt similar mode to substitute to described specific embodiment, but can't depart from spirit of the present invention or surmount the scope that appended claims defines.
Claims (2)
1., based on an earth's surface surface area extracting method for earth coordinates digital elevation model, it is characterized in that: carry out following steps,
Step one, establish original DEM adopt earth coordinates store, DEM grid cell, raster cell line number is M, and columns is N, and raster resolution is C, and described DEM is digital elevation model; Create (M-1) × 2 row, the raster file of sky that (N-1) × 2 arrange, in raster file, the value of each grid cell, raster cell is defaulted as 0, and the coordinate system of raster file is identical with original DEM; If the starting point latitude and longitude coordinates of original DEM is (B1, L1), then the starting point latitude and longitude coordinates of raster file is (B1+C, L1+C);
Step 2, the grid cell, raster cell of original DEM to be traveled through successively, builds regular triangle and to calculate the surface area in new grid cell, raster cell region as follows, carry out following sub-step,
Step 2.1, if the row number that i represents the line number of current traversal, j represents current traversal, then the grid cell, raster cell ranks coordinate of current traversal is (i, j), and making i initial value be 2, j initial value is 2,
Step 2.2, determines current grid pixel point according to current (i, j), and obtains ranks number for (i-1, j-1), and (i, j-1) and (i-1, j) grid cell, raster cell point, forms 4 adjacent points;
Step 2.3,4 points adjacent according to step 2.2 gained obtain a quadrilateral in the plane, and interpolation goes out the mid point height value of quadrilateral four edges and the height value of diagonal line intersection point, obtains 5 points;
Step 2.4, in 4 points adjacent by step 2.2 gained and step 2.3 gained 4 corresponding quadrilaterals of point, the mid point of four edges and diagonal line intersection point, obtain 8 spatially adjacent triangles;
Step 2.5, to the latitude and longitude coordinates of step 2.2 and step 2.3 gained totally 9 points, is converted to rectangular space coordinate according to following formula,
X＝(N+H)×Cos(B)×Cos(L)
Y＝(N+H)×Cos(B)×Sin(L)
Z＝[N×(1-e
2)+H]×Sin(B)
In formula, B, L, H represent latitude, longitude and height value of certain point respectively, and X, Y, Z represent this coordinate figure in three dimensions rectangular coordinate system respectively; The cosine value that Cos (B) is latitude, Sin (B) is the sine value of latitude, and Sin (L) is the sine value of longitude, and Cos (L) is the cosine value of longitude; e
2for earth's spheroid the first excentricity square, the computing formula of Parameter N is as follows,
In formula, a is the major semi-axis of earth's spheroid;
Step 2.6, utilizes Heron's formula to calculate 8 leg-of-mutton areas respectively, comprises and obtains area according to the rectangular space coordinate on an Atria summit;
Step 2.7, obtains quadrilateral by adjacent triangle combination of two, obtains 4 new quadrilaterals altogether;
Step 2.8, according to the component relationship of 8 triangles and 4 new quadrilaterals, gathers the area that leg-of-mutton area obtains 4 new quadrilaterals respectively;
Step 2.9, is written in the corresponding ranks number (2i-3,2j-3) of new raster file, (2i-2,2j-3), (2i-2,2j-2), (2i-3,2j-2) by the area of 4 new quadrilaterals; The ranks coordinate getting the next one traversal of original DEM is new (i, j), and return step 2.2, until traveled through, the raster file obtained is defined as digitally surface area model, is labeled as SAM, enters step 3;
Step 3, given any vector polygon, carry out Spatial Overlap by vector polygon and SAM, gathers in SAM the pixel value dropping on grid cell, raster cell in vector polygon and obtain V, using the surface area of V as this vector polygon; Wherein, whether grid cell, raster cell drops in vector polygon judges according to following rule,
If the central point of grid cell, raster cell drops in vector polygon or on vector polygon limit, then this grid cell, raster cell is calculated in vector polygon, otherwise does not exist.
2., based on an earth's surface surface area extraction system for earth coordinates digital elevation model, it is characterized in that: comprise with lower unit,
Raster file sets up unit, earth coordinates are adopted to store for establishing original DEM, DEM grid cell, raster cell line number is M, columns is N, raster resolution is C, and described DEM is digital elevation model, creates (M-1) × 2 row, the raster file of sky that (N-1) × 2 arrange, in raster file, the value of each grid cell, raster cell is defaulted as 0, and the coordinate system of raster file is identical with original DEM; If the starting point latitude and longitude coordinates of original DEM is (B1, L1), then the starting point latitude and longitude coordinates of raster file is (B1+C, L1+C);
Grid cell, raster cell region surface area extraction unit, for traveling through successively the grid cell, raster cell of original DEM, building regular triangle and the surface area in computation grid pixel region, comprising with lower module,
Initialization module, for establishing the row number that i represents the line number of current traversal, j represents current traversal, then the grid cell, raster cell ranks coordinate of current traversal is (i, j), and making i initial value be 2, j initial value is 2,
Consecutive point extraction module, for determining current grid pixel point according to current (i, j), and obtain ranks number for (i-1, j-1), (i, j-1) and (i-1, j) grid cell, raster cell point, forms 4 adjacent points;
Interpolating module, obtain a quadrilateral in the plane for 4 points adjacent according to consecutive point extraction module gained, interpolation goes out the mid point height value of quadrilateral four edges and the height value of diagonal line intersection point, obtains 5 points;
Triangular modeling block, for by the mid point of four edges in adjacent 4 points of consecutive point extraction module gained and interpolating module gained 4 corresponding quadrilaterals of point and diagonal line intersection point, obtains 8 spatially adjacent triangles;
Rectangular space coordinate modular converter, for the latitude and longitude coordinates to consecutive point extraction module and interpolating module gained 9 points, is converted to rectangular space coordinate according to following formula,
X＝(N+H)×Cos(B)×Cos(L)
Y＝(N+H)×Cos(B)×Sin(L)
Z＝[N×(1-e
2)+H]×Sin(B)
In formula, B, L, H represent latitude, longitude and height value of certain point respectively, and X, Y, Z represent this coordinate figure in three dimensions rectangular coordinate system respectively; The cosine value that Cos (B) is latitude, Sin (B) is the sine value of latitude, and Sin (L) is the sine value of longitude, and Cos (L) is the cosine value of longitude; e
2for earth's spheroid the first excentricity square, the computing formula of Parameter N is as follows,
In formula, a is the major semi-axis of earth's spheroid;
Triangle area extraction module, for utilizing Heron's formula to calculate 8 leg-of-mutton areas respectively, comprising and obtaining area according to the rectangular space coordinate on an Atria summit;
Quadrilateral builds module, for obtaining quadrilateral by adjacent triangle combination of two, obtains 4 new quadrilaterals altogether;
Quadrilateral area extraction module, for the component relationship according to 8 triangles and 4 new quadrilaterals, gathers the area that leg-of-mutton area obtains 4 new quadrilaterals respectively;
Digitally surface area model generation module, for the area of 4 new quadrilaterals being written to corresponding ranks number (2i-3,2j-3), (2i-2,2j-3), (2i-2 of new raster file, 2j-2), in (2i-3,2j-2); The ranks coordinate getting the next one traversal of original DEM is new (i, j), and order consecutive point extraction module works, until traveled through, the raster file obtained is defined as digitally surface area model, is labeled as SAM, and order surface area determination module works;
Surface area determination module, for given any vector polygon, carries out Spatial Overlap by vector polygon and SAM, gathers in SAM the pixel value dropping on grid cell, raster cell in vector polygon and obtains V, using the surface area of V as this vector polygon; Wherein, whether grid cell, raster cell drops in vector polygon judges according to following rule,
If the central point of grid cell, raster cell drops in vector polygon or on vector polygon limit, then this grid cell, raster cell is calculated in vector polygon, otherwise does not exist.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510235733.8A CN104821013A (en) | 2015-05-11 | 2015-05-11 | Method and system for specific surface area extraction based on geodetic coordinate system digital elevation model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510235733.8A CN104821013A (en) | 2015-05-11 | 2015-05-11 | Method and system for specific surface area extraction based on geodetic coordinate system digital elevation model |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN104821013A true CN104821013A (en) | 2015-08-05 |
Family
ID=53731296
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510235733.8A Pending CN104821013A (en) | 2015-05-11 | 2015-05-11 | Method and system for specific surface area extraction based on geodetic coordinate system digital elevation model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104821013A (en) |
Cited By (31)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105651209A (en) * | 2016-02-05 | 2016-06-08 | 中测新图(北京)遥感技术有限责任公司 | Emergency obtaining method and device for designated region area |
| CN105718751A (en) * | 2016-02-03 | 2016-06-29 | 西北农林科技大学 | Deep scarification operation area calculating method based on area grid point coverage |
| CN105956994A (en) * | 2016-05-13 | 2016-09-21 | 山东理工大学 | Graph processing method and device based on rasterized superposition analysis |
| CN106023311A (en) * | 2016-05-09 | 2016-10-12 | 黄河勘测规划设计有限公司 | Method for improving three-dimensional landform generation precision |
| CN106123812A (en) * | 2016-08-14 | 2016-11-16 | 覃泽林 | The method and device of relief surface sugarcane acreage is obtained based on remote sensing image |
| CN106777117A (en) * | 2016-12-15 | 2017-05-31 | 南京师范大学 | A kind of automatic identifying method of horizontal stratum tectonic landform |
| CN106952324A (en) * | 2017-04-07 | 2017-07-14 | 山东理工大学 | The parallel overlap-add procedure device and method of vector polygon rasterizing |
| CN107239664A (en) * | 2017-06-08 | 2017-10-10 | 国家基础地理信息中心 | A kind of ground mulching surface area rasterizing computational methods and device |
| CN107958209A (en) * | 2017-11-16 | 2018-04-24 | 深圳天眼激光科技有限公司 | Illegal construction identification method and system and electronic equipment |
| CN108595572A (en) * | 2018-04-16 | 2018-09-28 | 中国海洋大学 | One kind being suitable for urbanized area extracted with high accuracy water system and Watershed Unit method |
| CN108648271A (en) * | 2018-04-11 | 2018-10-12 | 华中科技大学 | A kind of interpolation method generating complicated landform grid model based on GIS data |
| CN108765560A (en) * | 2018-05-18 | 2018-11-06 | 云南电网有限责任公司电力科学研究院 | A kind of method and device converting dem data to grid model |
| CN109493423A (en) * | 2018-09-27 | 2019-03-19 | 北京市天元网络技术股份有限公司 | The calculation method and device of the midpoint of three-dimensional earth model surface two o'clock |
| CN109559377A (en) * | 2018-11-19 | 2019-04-02 | 成都智库二八六信息技术有限公司 | A method of it is cut using longitude and latitude polygon and generates three-dimensional map model |
| CN109656998A (en) * | 2018-12-04 | 2019-04-19 | 长江勘测规划设计研究有限责任公司 | The method for calculating turbine-generator units efficiency based on operational monitoring data |
| CN109813255A (en) * | 2019-01-16 | 2019-05-28 | 苏州科技大学 | The area computation method of big trapezoidal segment on a kind of earth ellipsoid face |
| CN110245202A (en) * | 2019-06-03 | 2019-09-17 | 哈尔滨工程大学 | A kind of polar coordinate system target grid set method for solving closing on search based on iteration |
| CN110322556A (en) * | 2019-04-29 | 2019-10-11 | 武汉大学 | A kind of high-speed, high precision arrow grid Overlap Analysis method based on border cuts |
| CN110490830A (en) * | 2019-08-22 | 2019-11-22 | 中国农业科学院农业信息研究所 | A kind of agricultural remote sensing method for correcting image and system |
| CN110610539A (en) * | 2019-08-13 | 2019-12-24 | 北京中色地科测绘有限公司 | Stratum curved surface construction method, device, equipment and storage medium |
| CN111127646A (en) * | 2019-12-26 | 2020-05-08 | 西南林业大学 | Construction method and system of rasterized elevation curved surface for measuring height difference of landform |
| CN111125893A (en) * | 2019-12-12 | 2020-05-08 | 河海大学 | Non-dispersive water flow path simulation method based on DEM and flow collection |
| CN111462322A (en) * | 2020-03-27 | 2020-07-28 | 昆明市测绘研究院 | DEM-based urban plane coordinate system establishment method |
| CN111738119A (en) * | 2020-06-12 | 2020-10-02 | 云南电网有限责任公司电力科学研究院 | Forest fire area estimation method and device based on MODIS satellite monitoring |
| CN112632209A (en) * | 2020-12-29 | 2021-04-09 | 智慧足迹数据科技有限公司 | Data association method and device, computer equipment and readable storage medium |
| CN112651937A (en) * | 2020-12-24 | 2021-04-13 | 中国林业科学研究院资源信息研究所 | Layered crown identification method |
| CN112686945A (en) * | 2021-01-11 | 2021-04-20 | 阿坝州自然资源与科技信息研究所 | Web-based three-dimensional terrain surface area measurement method |
| CN113884036A (en) * | 2021-09-22 | 2022-01-04 | 浙江省水利河口研究院(浙江省海洋规划设计研究院) | Coordinate conversion method and device, electronic equipment and storage medium |
| CN114155350A (en) * | 2022-02-07 | 2022-03-08 | 北京恒济引航科技股份有限公司 | Underwater terrain erosion and deposition evolution analysis method and device for channel |
| CN116680871A (en) * | 2023-05-11 | 2023-09-01 | 中国科学院空天信息创新研究院 | Global background radiation data and abundance acquisition method, device and equipment |
| CN117351066A (en) * | 2023-12-05 | 2024-01-05 | 深圳市法本信息技术股份有限公司 | Land area measurement method, apparatus, terminal device and storage medium |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102945570A (en) * | 2012-11-23 | 2013-02-27 | 华东师范大学 | Method for constructing full-space three-dimensional digital earth model |
-
2015
- 2015-05-11 CN CN201510235733.8A patent/CN104821013A/en active Pending
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102945570A (en) * | 2012-11-23 | 2013-02-27 | 华东师范大学 | Method for constructing full-space three-dimensional digital earth model |
Non-Patent Citations (5)
| Title |
|---|
| JEFF S.JENNESS: "Calculating landscape surface area from digital elevation models", 《WILDLIFE SOCIETY BULLETIN》 * |
| 国务院第一次全国地理国情普查领导小组办公室: "地理国情普查基本统计技术规定", 《第一次全国地理国情普查》 * |
| 曹振宇 等: "大地坐标系中多边形面积计算方法研究", 《测绘通报》 * |
| 江帆 等: "基于复化公式的DEM表面积算法分析", 《测绘学院学报》 * |
| 谢成磊 等: "基于地理坐标的地理国情统计单元表面面积精确计算", 《遥感信息》 * |
Cited By (49)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105718751A (en) * | 2016-02-03 | 2016-06-29 | 西北农林科技大学 | Deep scarification operation area calculating method based on area grid point coverage |
| CN105718751B (en) * | 2016-02-03 | 2017-04-12 | 西北农林科技大学 | Deep scarification operation area calculating method based on area grid point coverage |
| CN105651209B (en) * | 2016-02-05 | 2018-10-02 | 中测新图(北京)遥感技术有限责任公司 | The emergent acquisition methods and device of specified region area |
| CN105651209A (en) * | 2016-02-05 | 2016-06-08 | 中测新图(北京)遥感技术有限责任公司 | Emergency obtaining method and device for designated region area |
| CN106023311A (en) * | 2016-05-09 | 2016-10-12 | 黄河勘测规划设计有限公司 | Method for improving three-dimensional landform generation precision |
| CN106023311B (en) * | 2016-05-09 | 2018-12-28 | 黄河勘测规划设计有限公司 | Improve the method that dimensional topography generates precision |
| CN105956994A (en) * | 2016-05-13 | 2016-09-21 | 山东理工大学 | Graph processing method and device based on rasterized superposition analysis |
| CN106123812A (en) * | 2016-08-14 | 2016-11-16 | 覃泽林 | The method and device of relief surface sugarcane acreage is obtained based on remote sensing image |
| CN106777117A (en) * | 2016-12-15 | 2017-05-31 | 南京师范大学 | A kind of automatic identifying method of horizontal stratum tectonic landform |
| CN106777117B (en) * | 2016-12-15 | 2020-04-03 | 南京师范大学 | Automatic identification method for horizontal rock stratum structure landform |
| CN106952324A (en) * | 2017-04-07 | 2017-07-14 | 山东理工大学 | The parallel overlap-add procedure device and method of vector polygon rasterizing |
| CN107239664A (en) * | 2017-06-08 | 2017-10-10 | 国家基础地理信息中心 | A kind of ground mulching surface area rasterizing computational methods and device |
| CN107958209A (en) * | 2017-11-16 | 2018-04-24 | 深圳天眼激光科技有限公司 | Illegal construction identification method and system and electronic equipment |
| CN107958209B (en) * | 2017-11-16 | 2021-10-29 | 深圳天眼激光科技有限公司 | Illegal construction identification method and system and electronic equipment |
| CN108648271A (en) * | 2018-04-11 | 2018-10-12 | 华中科技大学 | A kind of interpolation method generating complicated landform grid model based on GIS data |
| CN108648271B (en) * | 2018-04-11 | 2020-07-14 | 华中科技大学 | Interpolation method for generating complex terrain grid model based on GIS data |
| CN108595572A (en) * | 2018-04-16 | 2018-09-28 | 中国海洋大学 | One kind being suitable for urbanized area extracted with high accuracy water system and Watershed Unit method |
| CN108765560A (en) * | 2018-05-18 | 2018-11-06 | 云南电网有限责任公司电力科学研究院 | A kind of method and device converting dem data to grid model |
| CN109493423A (en) * | 2018-09-27 | 2019-03-19 | 北京市天元网络技术股份有限公司 | The calculation method and device of the midpoint of three-dimensional earth model surface two o'clock |
| CN109493423B (en) * | 2018-09-27 | 2022-08-19 | 北京市天元网络技术股份有限公司 | Method and device for calculating midpoint positions of two points on surface of three-dimensional earth model |
| CN109559377A (en) * | 2018-11-19 | 2019-04-02 | 成都智库二八六信息技术有限公司 | A method of it is cut using longitude and latitude polygon and generates three-dimensional map model |
| CN109656998A (en) * | 2018-12-04 | 2019-04-19 | 长江勘测规划设计研究有限责任公司 | The method for calculating turbine-generator units efficiency based on operational monitoring data |
| CN109656998B (en) * | 2018-12-04 | 2023-04-07 | 长江勘测规划设计研究有限责任公司 | Method for calculating efficiency of hydroelectric generating set based on operation monitoring data |
| CN109813255B (en) * | 2019-01-16 | 2024-04-09 | 苏州科技大学 | Method for calculating area of large trapezoid block on earth ellipsoid |
| CN109813255A (en) * | 2019-01-16 | 2019-05-28 | 苏州科技大学 | The area computation method of big trapezoidal segment on a kind of earth ellipsoid face |
| CN110322556A (en) * | 2019-04-29 | 2019-10-11 | 武汉大学 | A kind of high-speed, high precision arrow grid Overlap Analysis method based on border cuts |
| CN110245202B (en) * | 2019-06-03 | 2023-07-21 | 哈尔滨工程大学 | Polar coordinate system target grid set solving method based on iterative proximity search |
| CN110245202A (en) * | 2019-06-03 | 2019-09-17 | 哈尔滨工程大学 | A kind of polar coordinate system target grid set method for solving closing on search based on iteration |
| CN110610539A (en) * | 2019-08-13 | 2019-12-24 | 北京中色地科测绘有限公司 | Stratum curved surface construction method, device, equipment and storage medium |
| CN110490830B (en) * | 2019-08-22 | 2021-09-24 | 中国农业科学院农业信息研究所 | Agricultural remote sensing image correction method and system |
| CN110490830A (en) * | 2019-08-22 | 2019-11-22 | 中国农业科学院农业信息研究所 | A kind of agricultural remote sensing method for correcting image and system |
| CN111125893A (en) * | 2019-12-12 | 2020-05-08 | 河海大学 | Non-dispersive water flow path simulation method based on DEM and flow collection |
| CN111127646A (en) * | 2019-12-26 | 2020-05-08 | 西南林业大学 | Construction method and system of rasterized elevation curved surface for measuring height difference of landform |
| CN111127646B (en) * | 2019-12-26 | 2023-03-14 | 西南林业大学 | Construction method and system of rasterized elevation curved surface for measuring height difference of landform |
| CN111462322A (en) * | 2020-03-27 | 2020-07-28 | 昆明市测绘研究院 | DEM-based urban plane coordinate system establishment method |
| CN111462322B (en) * | 2020-03-27 | 2021-04-02 | 昆明市测绘研究院 | DEM-based urban plane coordinate system establishment method |
| CN111738119B (en) * | 2020-06-12 | 2024-02-09 | 云南电网有限责任公司电力科学研究院 | Forest mountain fire area estimation method and device based on MODIS satellite monitoring |
| CN111738119A (en) * | 2020-06-12 | 2020-10-02 | 云南电网有限责任公司电力科学研究院 | Forest fire area estimation method and device based on MODIS satellite monitoring |
| CN112651937B (en) * | 2020-12-24 | 2024-02-20 | 中国林业科学研究院资源信息研究所 | Layered crown identification method |
| CN112651937A (en) * | 2020-12-24 | 2021-04-13 | 中国林业科学研究院资源信息研究所 | Layered crown identification method |
| CN112632209A (en) * | 2020-12-29 | 2021-04-09 | 智慧足迹数据科技有限公司 | Data association method and device, computer equipment and readable storage medium |
| CN112686945A (en) * | 2021-01-11 | 2021-04-20 | 阿坝州自然资源与科技信息研究所 | Web-based three-dimensional terrain surface area measurement method |
| CN113884036A (en) * | 2021-09-22 | 2022-01-04 | 浙江省水利河口研究院(浙江省海洋规划设计研究院) | Coordinate conversion method and device, electronic equipment and storage medium |
| CN113884036B (en) * | 2021-09-22 | 2023-09-29 | 浙江省水利河口研究院(浙江省海洋规划设计研究院) | Coordinate conversion method and device, electronic equipment and storage medium |
| CN114155350B (en) * | 2022-02-07 | 2022-04-22 | 北京恒济引航科技股份有限公司 | Underwater terrain erosion and deposition evolution analysis method and device for channel |
| CN114155350A (en) * | 2022-02-07 | 2022-03-08 | 北京恒济引航科技股份有限公司 | Underwater terrain erosion and deposition evolution analysis method and device for channel |
| CN116680871A (en) * | 2023-05-11 | 2023-09-01 | 中国科学院空天信息创新研究院 | Global background radiation data and abundance acquisition method, device and equipment |
| CN116680871B (en) * | 2023-05-11 | 2024-03-12 | 中国科学院空天信息创新研究院 | Global background radiation data and abundance acquisition method, device and equipment |
| CN117351066A (en) * | 2023-12-05 | 2024-01-05 | 深圳市法本信息技术股份有限公司 | Land area measurement method, apparatus, terminal device and storage medium |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104821013A (en) | Method and system for specific surface area extraction based on geodetic coordinate system digital elevation model | |
| US11226431B2 (en) | Method and device for filling invalid regions of terrain elevation model data | |
| CN107180450A (en) | A kind of algorithm of the river valley transverse shape based on DEM | |
| Li et al. | Using CA-Markov model to model the spatiotemporal change of land use/cover in Fuxian Lake for decision support | |
| CN108763825B (en) | Numerical simulation method for simulating wind field of complex terrain | |
| CN103544390A (en) | Cellular automata based rapid outburst flood routing simulation and analysis method | |
| CN107665242A (en) | A kind of multiple dimensioned lattice encoding method of regional space and device | |
| CN110544305B (en) | Terrain steeple ridge line information fusion method constructed by DEM (digital elevation model) oriented to regular grids | |
| CN115796712B (en) | Regional land ecosystem carbon reserve estimation method and device and electronic equipment | |
| Meng | Raster data projection transformation based-on Kriging interpolation approximate grid algorithm | |
| CN102289581A (en) | Method for simulating city expansion based on space function division | |
| Herzog et al. | Least-cost kernel density estimation and interpolation-based density analysis applied to survey data | |
| CN105761310A (en) | Simulated analysis and image display method of digital map of sky visible range | |
| CN103455709B (en) | A kind of super-resolution method for digital elevation model and system thereof | |
| CN112528508A (en) | Electromagnetic visualization method and device | |
| Cervilla et al. | Total 3D‐Viewshed Map: Quantifying the Visible Volume in Digital Elevation Models | |
| Miky et al. | A combined contour lines iteration algorithm and Delaunay triangulation for terrain modeling enhancement | |
| CN102880753B (en) | Based on the land utilization space characteristic dimension conversion method of fractal dimension | |
| CN117036644A (en) | Hexagonal grid construction method for remote sensing satellite region observation task | |
| CN106600691A (en) | Multipath 2D video image fusion correction method and system in 3D geographical space | |
| Bothale et al. | Evaluation and comparison of multi resolution DEM derived through Cartosat-1 stereo pair–A case study of Damanganga Basin | |
| CN114969944A (en) | High-precision road DEM construction method | |
| Yao | Application of GIS remote sensing information integration in eco-environmental quality monitoring | |
| Smith et al. | Map Comparison Methods for Three‐Dimensional Space and Time Voxel Data | |
| CN113743027A (en) | Method and device for drawing wind resource map based on CFD technology |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| EXSB | Decision made by sipo to initiate substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20150805 |
|
| WD01 | Invention patent application deemed withdrawn after publication |