A directional hierarchical mesh element interpolation method

Abstract

The present invention discloses a directional hierarchical mesh element interpolation method, which reads river data and known point data; classifies river data; defines category weights, flow direction weights and upstream and downstream weights; grids rivers to obtain grid rivers and grid river categories; constructs a directed network; associates known point IDs with river names; calculates spatialized category weights, spatialized flow direction weights, and spatialized upstream and downstream weights corresponding to known points; and performs interpolation calculations using known point values ​​and corresponding spatialized category weights, spatialized flow direction weights, and spatialized upstream and downstream weights. The present invention overcomes the shortcomings of general interpolation methods that do not consider element boundaries, different levels, and directions when interpolating on networkable elements.

Description

Technical Field

[0001] The present invention relates to the field of spatial analysis technology, and more specifically to a directional hierarchical mesh element interpolation method. Background Art

[0002] Interpolation is one of the commonly used methods in spatial analysis. It is the process of inferring unknown points associated with known points through known point data. It is mainly used to deal with problems such as missing data and data spatialization. It is widely used in environmental assessment, climate forecasting, geography, digital image processing, and mechanical manufacturing. In spatial analysis, common usage methods include filling in missing sequence data, interpolating spatially distributed point data into continuous surface data, and removing noise from image or video data.

[0003] Common interpolation methods, such as the inverse distance weight method and the kernel function method, calculate the weight by substituting the distance into the corresponding power exponential function or kernel function, and calculate the attribute value of the unknown point by combining the attribute value of the known point. These methods can convert the data of meteorological monitoring points and soil sampling points into surface data. However, common interpolation methods often consider the restrictions of spatial elements less. For example, when the element is a mesh structure, the change of attributes is usually along the network path, rather than along the radius direction of the known point. In some studies, the path direction is approximated by connecting the starting point and the end point for calculation, but this will have a relatively large error when the element is more curved. At the same time, the network often has different hierarchical structures. For example, the capacity of a certain element of the first-level network and the capacity of a certain element of the second-level network have an order of magnitude difference. If the interpolation calculation is performed without distinction, the value of the first-level network area will far exceed the limit. Moreover, networks such as rivers have directionality. In some processing, the weight difference between upstream and downstream or upstream and downstream directions must also be considered. Summary of the invention

[0004] The technical problem to be solved by the present invention is that, in view of the shortcomings of existing spatial interpolation methods in which the boundaries of elements, the hierarchical structure and direction of the network are less considered when interpolating on networkable elements, a directional hierarchical mesh element interpolation method is provided.

[0005] The above-mentioned object of the present invention is achieved by the following technical solutions:

[0006] A directional hierarchical mesh element interpolation method includes the following steps:

[0007] Step S1, reading river data and known point data;

[0008] Step S2, classifying the river data, and recording the categories obtained by classification into the categories of the river data;

[0009] Step S3, defining the category weights of river data of different categories, defining the flow direction weights of known points, and defining the upstream and downstream weights of known points;

[0010] Step S4, rasterizing the river with the unique river name and the category of the river data to obtain raster rivers and raster river categories;

[0011] Step S5: construct a directed network according to the river structure;

[0012] Step S6, associating the known point ID with the river name, and recording the river name of the river where the known point is located into the river name attribute of the known point data;

[0013] Step S7, calculating the corresponding path distance through the known point and the grid river network, and calculating the spatialized category weight corresponding to the known point, the spatialized flow direction weight corresponding to the known point, and the spatialized upstream and downstream weight corresponding to the known point;

[0014] Step S8: perform interpolation calculation using known point values ​​and corresponding spatialized category weights, spatialized flow direction weights, and spatialized upstream and downstream weights.

[0015] The step S5 as described above comprises the following steps:

[0016] Step S51, determining the river section ID in the downstream direction of the river according to the river structure, and recording the river section ID in the downstream direction of the river into the downstream direction of the river section name attribute of the river data;

[0017] Step S52: Establish a directed network according to the river name and the downstream river section name attributes.

[0018] The step S7 as described above comprises the following steps:

[0019] S71, determining the row and column numbers of the known point on the grid where the grid river is located by using the coordinates of the known point;

[0020] S72, assuming that the coordinates of the known point p1 are (x1, y1), calculate the row and column number (h1, c1) of the known point p1 on the grid where the grid river is located, where the distance value between the known point p1 and the grid is set to 0*hindrance1, hindrance1 is the resistance value at the row and column number (h1, c1), and add the grid with the row and column number (h1, c1) as the starting grid to the completed record sequence;

[0021] S73, find the neighboring grid of the starting grid (h1, c1), calculate the grid distance between the starting grid (h1, c1) and the neighboring grid that is not empty and not in the completed record sequence, the row and column numbers of the neighboring grid are recorded as (hn, cn), and the grid distance of the neighboring grid in the horizontal and vertical positions is calculated as np.sqrt(cellsize 2 +(dm1-dmn) 2 )*hindranceN, np.sqrt is the square root function. For the grid distance of the diagonal domain grid, the calculation formula is np.sqrt(2cellsize 2 +(dm1-dmn) 2 )*hindranceN, where dm1 is the elevation value of the starting grid (h1, c1), dmn is the elevation value of the neighborhood grid (hn, cn), hindranceN is the resistance value of the neighborhood grid, select the grid closest to the known point p1 in the neighborhood grid, record it as the nearest grid (h2, c2), add the nearest grid (h2, c2) to the completed record sequence completed, and record the grid distance from the starting grid (h1, c1) to all neighborhood grids to the dictionary D;

[0022] S74. Take the nearest grid (h2, c2) as the new starting grid, calculate the grid distance between the new starting grid and the neighboring grid that is not empty and not in the completed record sequence, and record the neighboring grid as (hn, cn). For the grid distance of the horizontal and vertical position grids, the calculation formula is d2+np.sqrt(cellsize 2 +(dm1-dmn) 2 )*hindranceN, for the diagonal grid distance, the calculation formula is d2+np.sqrt(2cellsize 2 +(dm1-dmn) 2 )*hindranceN, where d2 is the grid distance from the new starting grid to the original starting grid. For the neighborhood grids that have already recorded grid distances, if the newly calculated grid distance is less than the recorded minimum grid distance corresponding to the neighborhood grid, then the corresponding minimum grid distance recorded by the neighborhood grid is updated to the newly calculated grid distance. Otherwise, no processing is done. For the neighborhood grids that are not recorded in the dictionary D, the grid distance of the neighborhood grid is taken as the minimum grid distance and recorded in the dictionary D. The neighborhood grid with the smallest grid distance in the dictionary D is selected as the new starting grid.

[0023] S75, looping through step S74 until there is no new starting point grid, and the path distance calculation is completed;

[0024] S76, querying the grid river name corresponding to the grid river through the known point ID, and querying the category corresponding to the river section to which the known point belongs using the grid river name;

[0025] S77, query the category obtained in step S76, query the corresponding category weight, and calculate the spatialized category weight of the known point in combination with the raster river category mask, recorded as Wh;

[0026] S78, querying the downstream river section of the river section in the directed network by the river name, and calculating and obtaining the spatialized flow direction weight of the known point by combining the grid river and the flow direction weight, which is recorded as Wd;

[0027] S79. Calculate the reverse graph of the directed network through the directed network, obtain the upstream river section of the river section through the reverse graph of the directed network, and combine the grid river and the upstream and downstream weights to calculate and obtain the spatialized upstream and downstream weights of the known point, which is recorded as Wud.

[0028] The step S8 described above comprises the following steps:

[0029] Step S81, calculate the comprehensive weight W(i) of the i-th known point, W(i) = 1 / D(i)*Wh(i)*Wd(i)*Wud(i), Wh(i) is the spatialized category weight corresponding to the i-th known point, Wd(i) is the spatialized flow direction weight corresponding to the i-th known point, Wud(i) is the spatialized upstream and downstream weight corresponding to the i-th known point, D(i) is the shortest distance from the grid where the known point i is located to the domain grid,

[0030] S82. Calculate the interpolation result result, result = sum(W(i)*value(i) / sum(W(i))), value(i) is the value of the i-th known point, and sum is the summation function.

[0031] Compared with the prior art, the present invention has the following advantages:

[0032] The invention implements a directional hierarchical mesh element interpolation method, which realizes interpolation that takes into account the boundaries, levels and directions of networkable elements. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 is a flow chart of the present invention;

[0034] Figure 2 It is a schematic diagram of the river and known points;

[0035] Figure 3 It is a schematic diagram of a directed network;

[0036] Figure 4A simple example diagram for path distance calculation;

[0037] Figure 5 Schematic diagram of the results. DETAILED DESCRIPTION

[0038] In order to facilitate those skilled in the art to understand and implement the present invention, the present invention is further described in detail below in conjunction with embodiments. It should be understood that the implementation examples described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0039] like Figure 1 As shown, a directional hierarchical mesh element interpolation method includes the following steps:

[0040] S1. Read the river data CJriver and the known point data CJpoint.

[0041] Specifically, in this embodiment, the river data CJriver is the river data of the main water system network in the Yangtze River Basin, and the known point data CJpoint is the corresponding virtual point data (refer to Figure 2 ) to perform interpolation calculation.

[0042] S2. Classify the river data CJriver according to actual needs, and record the classified categories into the category class of the river data CJriver.

[0043] Specifically, when implementing step S2, those skilled in the art know that the category determination can be automatically determined by combining multiple attribute parameters of the river data CJriver, such as clustering or classification of multiple attribute parameters, or the actual category can be directly assigned. For example, when the category is represented by water quality grade, the water quality grade value is directly used as the category.

[0044] S3. Define the category weights of different categories of river data CJriver, define the flow direction weights of known points, and define the upstream and downstream weights of known points.

[0045] Specifically, when implementing step S3, class weights are set separately for different classes.

[0046] For example, the category class of the river data CJriver is [1,2,3], and the category weight of the river data CJriver with category class 1 is [w1,w2,w3], that is, when the known point falls on the river section with category class 1, the calculation weight of the river sections of different categories during interpolation is [w1,w2,w3]; the category weight of the river data CJriver with category class 2 is [w4,w5,w6]; the category weight parameter of the river data CJriver with category class 3 is [w7,w8,w9];

[0047] The flow direction weights of known points include Wal and Wan. Wal represents the weight of the river section that can be reached in the downstream direction of the river section where the known point is located, and Wan represents the weight of the river section that cannot be reached in the downstream direction of the river section where the known point is located.

[0048] The upstream and downstream weights of a known point include Wup and Wdown. Wup represents the weight of the river section located upstream of the known point, and Wdown represents the weight of the river section located downstream of the known point.

[0049] S4, rasterizing the river using the unique river name riverID and the category class of the river data, and obtaining the raster river rasterRiver and the raster river category rasterClass;

[0050] Specifically, when implementing step S4, rasterization can be performed through software (ArcGIS) or programming (GDAL), and the size of the raster is determined by actual needs.

[0051] S5. According to the river structure, construct a directed network riverG (refer to Figure 3 ):

[0052] S51. Determine the river section ID in the downstream direction of the river according to the river structure, and record the river section ID in the downstream direction of the river into the downstream direction river section name attribute toRiverID of the river data CJriver. If there is no downstream record, the value is -9999.

[0053] S52. A directed network is established according to the river name riverID and the downstream direction river section name attribute toRiverID.

[0054] S6. Associate the known point ID with the river name riverID, and record the river name riverID of the river where the known point is located into the river name attribute onRiverID of the known point data CJpoint.

[0055] S7. Obtain weights by calculating the path distance between the known point and the grid river network and the directed network riverG:

[0056] S71. Determine the row and column number (h, c) of the known point on the grid where the rasterRiver is located by using the coordinates (x, y) of the known point, where h is equal to ceil(abs(y-yc) / cellsize), ceil is a rounding function, abs() is an absolute value function, c is equal to ceil(abs(x-xc) / cellsize), (xc, yc) is the coordinate of the upper left corner of the grid where the rasterRiver is located, and cellsize is the resolution;

[0057] S72. Calculate the path distance. Take a single-point calculation as an example. Assume that the known point p1 is any known point. The coordinates of the known point p1 are (x1, y1). Calculate the row and column number (h1, c1) of the known point p1 on the grid where the rasterRiver is located. Here, the distance value between the known point p1 and the grid is set to 0*hindrance1. Hindrance1 is the resistance value at the row and column number (h1, c1). That is, it is assumed that the distance between the known point and the grid it falls into is 0. The grid with the row and column number (h1, c1) is added as the starting grid to the completed record sequence completed.

[0058] S73, find the neighboring grid of the starting grid (h1, c1), calculate the grid distance between the starting grid (h1, c1) and the neighboring grid that is not empty and not in the completed record sequence completed, the row and column numbers of the neighboring grid are recorded as (hn, cn), and the grid distance of the neighboring grid in the horizontal and vertical positions is calculated as np.sqrt(cellsize 2 +(dm1-dmn) 2 )*hindranceN, np.sqrt is the square root function. For the grid distance of the diagonal domain grid, the calculation formula is np.sqrt(2cellsize 2 +(dm1-dmn) 2 )*hindranceN, where dm1 is the elevation value of the starting grid (h1, c1), dmn is the elevation value of the neighborhood grid (hn, cn), hindranceN is the resistance value of the neighborhood grid, select the grid closest to the known point p1 in the neighborhood grid, record it as the nearest grid (h2, c2), add the nearest grid (h2, c2) to the completed record sequence completed, and record the grid distances d1, d2, d3, d4... from the starting grid (h1, c1) to all neighborhood grids to the dictionary D;

[0059] S74. Take the nearest grid (h2, c2) as the new starting grid, calculate the grid distance between the new starting grid and the neighboring grid that is not empty and not in the completed record sequence completed, and record the neighboring grid as (hn, cn). For the grid distance of the horizontal and vertical position grids, the calculation formula is d2+np.sqrt(cellsize 2 +(dm1-dmn) 2 )*hindranceN, for the diagonal grid distance, the calculation formula is d2+np.sqrt(2cellsize 2 +(dm1-dmn) 2 )*hindranceN, where d2 is the grid distance from the new starting grid to the original starting grid. For the neighborhood grids that have already recorded grid distances, if the newly calculated grid distance is less than the recorded minimum grid distance corresponding to the neighborhood grid, then the corresponding minimum grid distance recorded by the neighborhood grid is updated to the newly calculated grid distance. Otherwise, no processing is done. For the neighborhood grids that are not recorded in the dictionary D, the grid distance of the neighborhood grid is taken as the minimum grid distance and recorded in the dictionary D. The neighborhood grid with the smallest grid distance in the dictionary D is selected as the new starting grid, recorded as (h3,c3);

[0060] S75, loop through step S74 until there is no new starting grid, and the path distance calculation is completed ( Figure 4 Provide a simple example);

[0061] S76, querying the raster river name riverID of the corresponding raster river rasterRiver through the known point ID, and querying the class class corresponding to the river section to which the known point belongs using the raster river name riverID;

[0062] S77, query the corresponding category weight through the category obtained in step S76, and calculate the spatialized category weight of the known point in combination with the rasterClass mask of the river category, which is recorded as Wh;

[0063] S78, query the downstream river section of the river section in the directed network riverG through the river name riverID, and calculate and obtain the spatialized flow direction weight of the known point by combining the raster river rasterRiver and the flow direction weight, which is recorded as Wd;

[0064] S79. Calculate the reverse graph riverGre of riverG through riverG, obtain the upstream section of the river section through riverGre, and combine rasterRiver and the upstream and downstream weights to calculate and obtain the spatialized upstream and downstream weights of the known point, which is recorded as Wud.

[0065] S8. Interpolation calculation is performed using known point values ​​and corresponding spatialized category weights, spatialized flow direction weights, and spatialized upstream and downstream weights:

[0066] S81. Calculate the comprehensive weight W(i) of the ith known point, W(i) = 1 / D(i)*Wh(i)*Wd(i)*Wud(i). Here, you can choose to use a fixed distance or the nth-order neighborhood of the river section for masking. Wh(i) is the spatialized category weight corresponding to the ith known point, Wd(i) is the spatialized flow direction weight corresponding to the ith known point, Wud(i) is the spatialized upstream and downstream weight corresponding to the ith known point, and D(i) is the shortest distance from the grid where the known point i is located to the domain grid.

[0067] S82, calculate the interpolation result by weight (refer to Figure 5 ), result = sum(W(i)*value(i) / sum(W(i))), value(i) is the value of the i-th known point, and sum is the summation function.

[0068] The invention discloses a directional hierarchical mesh element interpolation method, which comprehensively considers element boundaries, network hierarchical structure and direction, etc., and realizes the interpolation of networkable elements.

[0069] It should be noted that the specific embodiments described in the present invention are merely examples of the spirit of the present invention. Those skilled in the art may make various modifications or additions to the specific embodiments described or replace them in similar ways, but they will not deviate from the spirit of the present invention or exceed the scope defined by the attached claims.

Claims

1. A directional hierarchical mesh element interpolation method, characterized in that: The following steps are involved: Step S1, reading river data and known point data; Step S2, classifying the river data, and recording the categories obtained by classification into the categories of the river data; Step S3, defining the category weights of river data of different categories, defining the flow direction weights of known points, and defining the upstream and downstream weights of known points; Step S4, rasterizing the river with the unique river name and the category of the river data to obtain raster rivers and raster river categories; Step S5: construct a directed network according to the river structure; Step S6, associating the known point ID with the river name, and recording the river name of the river where the known point is located into the river name attribute of the known point data; Step S7, calculating the corresponding path distance through the known point and the grid river network, and calculating the spatialized category weight corresponding to the known point, the spatialized flow direction weight corresponding to the known point, and the spatialized upstream and downstream weight corresponding to the known point; Step S8: interpolate the known point values ​​and the corresponding spatialized category weights, spatialized flow direction weights, and spatialized upstream and downstream weights. The step S7 comprises the following steps: S71, determining the row and column numbers of the known point on the grid where the grid river is located by using the coordinates of the known point; S72, assuming that the coordinates of the known point p1 are (x1, y1), calculate the row and column numbers (h1, c1) of the known point p1 on the grid where the river is located, where the distance value between the known point p1 and the grid is set to 0 hindrance1, hindrance1 is the resistance value at the row and column number (h1, c1), and the grid with the row and column number (h1, c1) is added as the starting grid to the completed record sequence; S73, find the neighboring grid of the starting grid (h1, c1), calculate the grid distance between the starting grid (h1, c1) and the neighboring grid that is not empty and not in the completed record sequence, the row and column numbers of the neighboring grid are recorded as (hn, cn), and the grid distance of the neighboring grid in the horizontal and vertical positions is calculated as np.sqrt(cellsize 2 +(dm1-dmn) 2 ) hindranceN, np.sqrt is the square root function. For the grid distance of the domain grid in the diagonal direction, the calculation formula is np.sqrt(2cellsize 2 +(dm1-dmn) 2 ) hindranceN, where dm1 is the elevation value of the starting grid (h1, c1), dmn is the elevation value of the neighborhood grid (hn, cn), hindranceN is the resistance value of the neighborhood grid, select the grid closest to the known point p1 in the neighborhood grid, record it as the nearest grid (h2, c2), add the nearest grid (h2, c2) to the completed record sequence completed, and record the grid distance from the starting grid (h1, c1) to all neighborhood grids to the dictionary D; S74. Take the nearest grid (h2, c2) as the new starting grid, calculate the grid distance between the new starting grid and the neighboring grid that is not empty and not in the completed record sequence, and record the neighboring grid as (hn, cn). For the grid distance of the horizontal and vertical position grids, the calculation formula is d2+np.sqrt(cellsize 2 +(dm1-dmn) 2 ) hindranceN, for the diagonal grid distance, the calculation formula is d2+np.sqrt(2cellsize 2 +(dm1-dmn) 2 ) hindranceN, where d2 is the grid distance from the new starting grid to the original starting grid. For the neighborhood grids that have already recorded grid distances, if the newly calculated grid distance is less than the recorded minimum grid distance corresponding to the neighborhood grid, the corresponding minimum grid distance recorded in the neighborhood grid is updated to the newly calculated grid distance. Otherwise, no processing is done. For the neighborhood grids that are not recorded in the dictionary D, the grid distance of the neighborhood grid is taken as the minimum grid distance and recorded in the dictionary D. The neighborhood grid with the smallest grid distance in the dictionary D is selected as the new starting grid. S75, looping through step S74 until there is no new starting point grid, and the path distance calculation is completed; S76, querying the grid river name corresponding to the grid river through the known point ID, and querying the category corresponding to the river section to which the known point belongs using the grid river name; S77, query the category obtained in step S76, query the corresponding category weight, and calculate the spatialized category weight of the known point in combination with the raster river category mask, recorded as Wh; S78, querying the downstream river section of the river section in the directed network by the river name, and calculating and obtaining the spatialized flow direction weight of the known point by combining the grid river and the flow direction weight, which is recorded as Wd; S79. Calculate the reverse graph of the directed network through the directed network, obtain the upstream river section of the river section through the reverse graph of the directed network, and combine the grid river and the upstream and downstream weights to calculate and obtain the spatialized upstream and downstream weights of the known point, which is recorded as Wud.

2. A directional hierarchical mesh element interpolation method according to claim 1, characterized in that: The step S5 comprises the following steps: Step S51, determining the river section ID in the downstream direction of the river according to the river structure, and recording the river section ID in the downstream direction of the river into the downstream direction of the river section name attribute of the river data; Step S52: Establish a directed network according to the river name and the downstream river section name attributes.

3. A directional hierarchical mesh element interpolation method according to claim 1, characterized in that: The step S8 comprises the following steps: Step S81: Calculate the comprehensive weight W(i) of the i-th known point, W(i) = 1 / D(i) Wh(i) Wd(i) Wud(i), Wh(i) is the spatialized category weight corresponding to the i-th known point, Wd(i) is the spatialized flow direction weight corresponding to the i-th known point, Wud(i) is the spatialized upstream and downstream weight corresponding to the i-th known point, and D(i) is the shortest distance from the grid where the known point i is located to the domain grid. S82, calculate the interpolation result result, result = sum(W(i) value(i) / sum(W(i))), value(i) is the value of the i-th known point, and sum is the summation function.

Images