Rectangular grid rainfall isoline generation and filling method and system based on digital analog simulation engine

Abstract

The invention relates to the technical field of scientific calculation visualization, discloses a rainfall isoline generation and filling method and system of a rectangular grid based on a digital analog simulation engine, and mainly realizes research on a generation algorithm and a filling algorithm of a rainfall isoline so as to achieve the purposes of improving the operation efficiency and improving the drawing accuracy. A drawing process of a contour line is divided into two major parts of generation and filling, generation of the contour line is realized step by step, and on the basis of previous research, a contour line tracking algorithm is improved, and a back-tracking function is added into a traditional tracking algorithm, so that secondary traversal of a rectangular grid is avoided, the operation efficiency is improved, and meanwhile, the calculation efficiency is improved. In the step of contour line smoothing, multiple Bezier curves are selected to fit a contour line, so that the problem of contradiction between line accuracy and smoothness is solved, the calculation amount is reduced, and the working intensity is relieved.

Description

Technical Field

[0001] The invention relates to the technical field of scientific computing visualization, and in particular to a method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine. Background Art

[0002] As an important input parameter of the hydrological model, rainfall is also regarded by researchers as an important reference factor for flood control, flood fighting, disaster reduction, irrigation and other issues. It is the main object of investigation and research in the field of hydrological, water conservancy and meteorological forecasting. The rainfall problem has always been a key issue for scholars in the field of hydrology and water resources. The most common method to obtain rainfall is to obtain the average surface rainfall in the region by setting up rainfall observation stations. Although this method has many disadvantages, it is still widely used in the field of water conservancy. In the work of flood control and drought relief, it is extremely important to publish and share rainfall data in real time so that users can obtain rainfall data in the area of ​​concern in real time. Rainfall contour lines are tools for visually observing rainfall conditions. Rainfall contour lines represent: the lines connecting points with equal rainfall values ​​within a given period of time, that is, the lines connecting equal points. Rainfall contour lines have also become an important tool for flood control, reservoir scheduling and other work due to their intuitive display of the geographical distribution of rainfall. They are also widely used in flood control command, hydrology and water conservancy departments.

[0003] Contour lines have been deeply and closely studied and applied in many fields such as meteorology, geology, and geological disaster monitoring. They effectively represent the size and distribution information of geometric field variables such as temperature, rainfall, groundwater level, and soil moisture. The distance between contour lines reflects the change in the size of field variables. It is widely used in the field of water conservancy information and plays a very important role in flood control decision-making. The earliest contour lines were calculated and drawn manually. This method of generating contour lines increased the workload of staff, repeated work, and was inefficient, which was not conducive to a wide range of professional applications and its commercial expansion. Some disadvantages of traditional manual drawing of contour lines: This low-end manual drawing method is slow and inefficient, and the quality of the final contour map is also unsatisfactory. Summary of the invention

[0004] 1. Technical issues to be resolved

[0005] In view of the deficiencies in the prior art, the present invention provides a method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine, which has the advantages of high efficiency and good quality in contour line drawing, and solves the problem of low efficiency in traditional manual contour line drawing.

[0006] (II) Technical solution

[0007] To achieve the above object, the present invention provides the following technical solutions: including grid modeling of discrete rainfall data, a generation algorithm of rainfall contour lines, and a filling algorithm of rainfall contour lines;

[0008] The grid modeling of discrete rainfall data includes a DEM model, a regular grid generated by discrete points, spatial interpolation of rainfall and Kriging interpolation, and the grid modeling of discrete rainfall data is built based on the DEM model, the regular grid generated by discrete points, spatial interpolation of rainfall and Kriging interpolation;

[0009] The generation algorithm of the rainfall contour line is based on the grid modeling of discrete rainfall data. The generation algorithm of the rainfall contour line includes contour lines, the generation and process of contour lines, rectangular grids and the process of realizing rectangular grids by the contour line algorithm. The generation algorithm of the rainfall contour line is based on the theoretical basis of contour lines, the generation and process of contour lines, rectangular grids and the process of realizing rectangular grids by the contour line algorithm.

[0010] The filling algorithm of the rainfall contour line is based on the generation algorithm of the rainfall contour line. The filling algorithm of the rainfall contour line includes the filling algorithm design and the annotation of the contour line. The filling algorithm of the rainfall contour line is based on the filling algorithm design and the annotation of the contour line as the practical basis to generate the filling algorithm.

[0011] Preferably, the rainfall of the DEM is expressed in a watershed distribution. The DEM model uses known three-dimensional coordinate points in an arbitrary coordinate system to perform relevant statistics on continuous ground. DEM is divided into three categories: one is a regular grid model, another is an irregular triangulated network model, and the other is a contour model. These three models can be converted and generated from each other.

[0012] The basin-distributed rainfall estimation requires the longitude and latitude information of each station, as well as the elevation information of each station;

[0013] DEM is a finite sequence of three-dimensional vectors representing area D:

[0014] {Vi=(xi,yi,zi),i=1,2,……,n}

[0015] Where Vi represents the value of the ground characteristic information corresponding to the pth grid; xi,yi represents the coordinate value of a data point on the xoy plane, and zi represents the elevation value corresponding to the data point (xi,yi).

[0016] Preferably, the data source of the digital elevation model is obtained by interpreting aerial photos, collecting GPS receivers, importing elevation data points and converting vector contour lines, and the spatial storage of the digital elevation model is completed in the form of regular grids, triangular patches and contour lines, each node of which has an x ​​and y coordinate value, the z value of the point determines the shape of the surface, and the broken line defines the change of the surface;

[0017] The basic content of the DEM is summarized by a set of ordered point pairs of a two-dimensional function:

[0018] Mp=fk(vp,wp)

[0019] k=1,2,3......m; p=1,2,3......n (2-1)

[0020] Wherein, Mp represents the value of the ground characteristic information corresponding to the pth grid point, vp, wp represent the two-dimensional plane coordinates of point p, and any plane coordinate system of map projection can be used, m is the number of ground characteristic information types, n is the number of ground points, Kp refers to a characteristic attribute of the ground, including soil moisture content, vegetation coverage and soil type. The digital rainfall model is a special model when m = 1 and the ground characteristic attribute represented by Mp is the rainfall value of each grid point;

[0021] The discrete point generated regular grid is a rectangular regular grid, the basic unit of the regular grid model is a square or rectangle, the spatial position of each grid unit is represented by the coordinate (x, y), and a value is assigned to each grid unit, and the discrete point generated regular grid is assigned a rainfall value;

[0022] The modeling steps of generating a regular grid from discrete points are as follows:

[0023] S1.1, first divide the data area into rectangular grids on a two-dimensional plane to form a plane grid covering the entire area;

[0024] S1.2. Using the sampling points distributed in the grid area, the rainfall value of each vertex of the grid unit is obtained through a specific mathematical algorithm, thereby obtaining the final rectangular grid model.

[0025] Preferably, the spatial interpolation of rainfall utilizes discrete monitoring data of a specific area to estimate non-monitoring data on regular grid points, and converts such discrete site observation data into regular grid data. The spatial interpolation of rainfall is divided into global interpolation and local interpolation. The local interpolation includes: Thiessen polygon method, inverse distance interpolation, spline function interpolation method and Kriging interpolation; the spatial interpolation of rainfall includes the interpolation principle of rectangular grids and the interpolation process of rectangular grids;

[0026] The interpolation principle of the rectangular grid is as follows: for the interpolation of the rectangular grid, a regular rectangular grid covering the point set Zi(xi,yi) is first established on the xoy plane, where Zi(xi,yi) is the projection point set of Zi(xi,yi,zi) on the xoy plane. For each grid node, an insertion point must be defined, and the defined insertion point satisfies the projection point on the xoy plane as the grid node. Find n points Zi(xi,yi) close to the grid node position, and realize the interpolation through the process of selecting the measured data points and calculating the influencing parameters w;

[0027] The interpolation process of the rectangular grid is:

[0028] (1) Select the measurement data points

[0029] When performing spatial data interpolation, measurement data points are selected from the vicinity of the insertion point to participate in the calculation. The calculation process has the following steps:

[0030] S2.1, taking the grid node as the center, find the projection point in the nearby area;

[0031] S2.2, determine whether the projection point found on the xoy plane is a valid point, and the valid point is the projection point of the measured data point on the xoy plane;

[0032] S2.3, find the corresponding measurement data point according to the found valid projection point;

[0033] (2) Calculation of influencing parameters

[0034] When calculating the rainfall value of an insertion point based on the found measurement data points, the farther the data point is from the insertion point, the smaller the influence. When it is far from the data point, the calculation of the influencing parameters is required;

[0035] Assume that the insertion point is Z(x,y,z), the data point sequence is Z1(x1,y1,z1)…Zn(xn,yn,zn), and the projection distance Di from the projection point of Z(x,y,z) to Zi(xi,yi,zi) is obtained by the following formula:

[0036] Di=(x-xi)*(x-xi)+(y-yi)*(y-yi) (2-2)

[0037] The impact factor wi of Zi is defined as

[0038] (3) Interpolation test standard.

[0039] Preferably, the interpolation test standard uses a cross-validation method to verify the effect of interpolation. The cross-validation method is to assume that the rainfall values ​​of some stations of known rainfall stations in the region are unknown variables, and then interpolate and estimate the rainfall values ​​of nearby stations, and then calculate the error between the actual monitoring value and the theoretical estimated value of the station to be estimated, and use this as a standard to judge the quality of the interpolation method;

[0040] The interpolation test standard uses average error, average absolute error, average relative error, mean square error and root mean square error as evaluation criteria for analysis, and the algorithm is as follows:

[0041]

[0042] Experimental comparison of several common interpolation methods: create a data table as the interpolation data source, the data table format is xls, convert xls into grid data format, and use this grid data to draw contour lines interpolated by different interpolation methods;

[0043] The Kriging interpolation method first needs to analyze the spatial field, propose a variogram model, and then carry out Kriging calculation based on the variogram model;

[0044] The key factors of the Kriging interpolation include regionalized variables and variograms. The spatial information with dual attributes of structure and uncertainty is called regionalized variables. If Z(x) is a regionalized variable and satisfies the second-order stationary assumption, the following formula is obtained:

[0045]

[0046] Among them, the determination of the weight coefficient λ is the core of the Kriging method. When solving the weight coefficient λ, it is necessary to fully ensure that the linear unbiasedness of the estimator Ze(u) and the minimum estimation variance are met.

[0047] Preferably, the contour lines:

[0048] (1) An isoline is a smooth, continuous curve;

[0049] (2) For a given attribute value, the corresponding number of contour lines is the corresponding number;

[0050] (3) Contour lines can be closed or open;

[0051] (4) The contour lines will not intersect each other;

[0052] The generation principle and process of the contour lines are as follows:

[0053] The data type for contour drawing is discrete two-dimensional data. The domain of discrete data is three-dimensional domain, and its data itself is a three-dimensional coordinate vector (x, y, z), where (x, y) is a plane vector and z is an attribute value used to represent the rainfall value. There are two main forms of discrete data distribution: one is a rectangular grid distribution data field, and the other is a scattered irregular distribution data field. The process of generating contour lines using the rectangular grid method is: first, a set of three-dimensional discrete data points (x, y, z) is interpolated by an algorithm to obtain the attribute values ​​of each rectangular grid point, and its equivalent points are obtained by linear interpolation on the boundary of the grid cell; then, the contour lines are tracked according to the tracking principle of the rectangular grid method and the contour line map is drawn. The generation process of the contour line is divided into three steps:

[0054] S3.1, Grid processing of discrete data;

[0055] S3.2, search for contour lines;

[0056] S3.3, display of contour map;

[0057] The rectangular grid method uses the collected discrete data points to establish a surface model on a regular area close to the geographic surface. The discrete points in the obtained model are interpolated by an interpolation algorithm to obtain the attribute value of each grid point. Then, according to the attribute values ​​of the adjacent points on the edge of each cell in the rectangular grid, the specific positions of all the equal-value points on the edge of each cell are obtained by linear interpolation. Finally, according to the tracking principle of the rectangular grid method, all the equal-value points of all the isolines are searched out and connected in sequence into smooth curves, namely: isolines;

[0058] The implementation of the contour line regular rectangular grid method includes the calculation of equal value points, the establishment of the connection order of equal value points and the tracking of contour lines, and its process is as follows:

[0059] Step 1: The generation of contour lines first requires obtaining discrete data through corresponding algorithms based on actual rainfall data;

[0060] Step 2: Then, interpolate the discrete data to obtain gridded data;

[0061] Step 3: Calculate the gridded data to obtain equivalent points;

[0062] Step 4: Track and smooth the points with equal value, and finally generate the equal value lines;

[0063] Calculation of the equivalent points:

[0064] The drawing area obtained by rectangular grid modeling is composed of m×n grid units. The xoy plane coordinate system is established, the coordinate of the lower left corner of the rectangular area is (0,0), and the basic unit lengths of the directional grid units are dx and dy respectively; according to formula (3-5-1) and formula (3-5-2), the horizontal and vertical coordinate values ​​of any grid point can be obtained as follows:

[0065] X=i×dx (i=0,1,2……m) (3-5-1)

[0066] Y=j×dy (j=0,1,2……n) (3-5-2)

[0067] Among them, dx and dy represent the length of the grid unit in the X direction and the Y direction respectively, i and j represent the number of grid units in the X direction and the Y direction respectively, and the formula is calculated based on the horizontal side pentagon and the vertical side pentagon obtained by connecting the two-sided triangles, and the formula is:

[0068] Horizontal side pentagon:

[0069]

[0070] Get the intersection coordinates

[0071] y=y A =y B (3-2)

[0072] Vertical pentagon

[0073]

[0074] Get the intersection coordinates x=x A =x B (3-3)

[0075]

[0076] According to the above steps, the equivalent points on all grid edges can be calculated;

[0077] The establishment of the connection order of the equivalent points: the attribute values ​​on the grid unit edges are linearly distributed. At this time, each grid unit edge has a maximum of 4 intersections with the equivalent line, and can have 2 intersections. At the same time, the proof by contradiction is used to prove that there is no one or three intersections between the equivalent line and the grid unit edge.

[0078] Preferably, the tracing of the contour lines includes contour line smoothing, cubic spline function, tension spline function and Bezier function;

[0079] The method for smoothing the contour line includes: the contour line smoothing algorithm includes tension spline function, linear iteration method, piecewise cubic polynomial method, Bezier curve function, B-spline method and cubic Bezier method, and the contour line smoothing algorithm satisfies: (1) forming a piecewise linear function; (2) forming a smooth curve through the fitted curve;

[0080] The cubic spline function:

[0081] Given a set of type value points Pi = (xi, yi), i = 0, 1, ..., n, construct two interpolating cubic spline functions with respect to the parameters:

[0082] x=x(u),y=y(u)

[0083] Interpolate at the point sets (u,x) and (u,y), i=0,1,2...n respectively, and then merge the three to form a cubic parametric spline curve:

[0084] P(u)=(x(u),y(y))

[0085] Determine contours by constructing parametric equations:

[0086]

[0087] The tension spline function is a piecewise interpolation function which is a linear combination of a straight line interpolation and two hyperbolic functions shσx and chσx:

[0088] f(x)=c1+c2x+c3shσx+c4chσx (3-6)

[0089] The whole has a continuous second-order derivative. The role of σ is to control the position of the inflection point and the shape of the curve. Through mathematical derivation, the function expression of the tension spline function is obtained:

[0090]

[0091] The Bezier function is composed of line segments and nodes, where nodes are draggable curves. The Bezier curve is defined as: given n+1 control vertices, its expression is:

[0092] P(t)=∑B i ,n(t)P i u∈[0,1] (3-8)

[0093] Where: Bi,n(t) is called basis function;

[0094] The harmonic function expression between the vertex position vectors of the Bezier polygon is:

[0095]

[0096] where: Bi,n(k) is the Bernstein basis function;

[0097] For the smooth condition of the splicing of the multi-segment Bezier curves, when using the Bezier curve fitting in segments, the following conditions need to be satisfied:

[0098] (1) The positions of two adjacent Bezier curves need to be continuous, and the first and last points of the front and back two curves should be connected;

[0099] (2) The slopes of two adjacent Bezier curves need to be continuous, that is, the starting points where the front and back two curves coincide should be collinear with two adjacent points;

[0100] (3) The curvatures of two adjacent Bezier curves need to be continuous, that is, the two curves should ensure that five points are coplanar. Preferably, for the design of the multi-segment Bezier curve fitting algorithm, the algorithm process of using the multi-segment Bezier curve to fit the contour line is as follows:

[0101] S4.1. Assume n + 1 discrete data points P0, P1, ……, Pn on the contour line. First, obtain two points P0 and P3 of the first segment of the Bezier curve. The initial point k for each segment to start fitting is 0 for the first segment. An initial step length is used to fit the number of discrete data points of one segment of the Bezier curve, and an error limit value E.

[0102] An initial value of M is equal to one step;

[0103] S4.3. Take the points from the k-th to the k + m-th in the discrete data point sequence. When k + m > n, then m = n - k points, that is, Pk, Pk+1 …… Pk+m, and use the least squares method to approximate the contour line segment with a Bezier curve;

[0104] S4.4. Calculate the error e between the points on the Bezier curve and the points on the corresponding contour line segment;

[0105] When e < E, the fitting of this segment is successful, and k = k + m; when k ≥ n, it ends, otherwise continue to execute step 2; when e ≥ E, during the fitting of this segment, by reducing the number of discrete data points for fitting m = 2*m / 3, continue to execute step 2 for fitting. The amount of discrete point data is reduced by 1 / 3. First, obtain two points P0 and P3 of the first segment of the Bezier curve.

[0106] Assume that the Bezier curve passes through n points on the broken line at ti, and the n points are S0, S1, …, Sn respectively. Then there is the formula:

[0107]

[0108] In the formula, let S0 = P0, S3 = P3, that is, find the two points P0 and P3 of the first section of the curve. After finding points P0 and P3, find points P1 and P2, and the formula is:

[0109]

[0110] Use the multi-segment idea to smooth the contour lines. According to the smooth condition 1, we can get T0 = P3. Use the least square method to find three control points, which are: T1, T2, T3. Let: T1 = T0 + G*R, where T0 = P3;

[0111] There is a formula:

[0112]

[0113] Take a vertex of the broken line as a double vertex, plus a point on the two preceding and following line segments to form the control polygon of the cubic Bezier curve, determine the start and end points of each cubic Bezier segment, and from the properties of the Bezier curve, it can be concluded that the start and end points of each cubic Bezier curve are the first and last points of the control polygon respectively.

[0114] Preferably, the filling algorithm is designed to divide the filling work into three steps:

[0115] S5.1, construct connected regions;

[0116] S5.2, sort the regions in order of their area size;

[0117] S5.3, filling each area with color according to the determined filling parameters;

[0118] The filling algorithm divides the contour lines into 11 categories according to the relationship between the contour lines and the boundaries:

[0119] (a) The intersection of the contour lines as boundaries is 0;

[0120] (b) The contour line only intersects the lower boundary;

[0121] (c) The contour line only intersects the right boundary;

[0122] (d) The contour line only intersects the left boundary;

[0123] (e) The contour lines intersect only the upper boundary;

[0124] (f) The contour lines intersect the right and lower boundaries;

[0125] (g) The contour lines intersect the left and lower boundaries;

[0126] (h) The contour lines intersect the left and upper boundaries;

[0127] (i) The contour line intersects the right and upper boundaries;

[0128] (j) The contour line intersects the left and right boundaries;

[0129] (k) The contour lines intersect the upper and lower boundaries;

[0130] The filling algorithm process is as follows:

[0131] S6.1, initialize the iso-value area, and classify the iso-value lines according to the relationship between the iso-value lines and the boundaries;

[0132] S6.2. Sort the contour lines from large to small according to their classification. If the classification is the same, in order to ensure that the colors can be filled in order, they must be sorted according to certain rules. Assume that the coordinates of the first point of the contour line are x, y, and the maximum values ​​of the points above, below, left and right of the contour line are maxy, miny, maxx, maxy;

[0133] S6.3. Determine the parent-child relationship between contour lines of the same type. Judgment of the parent-child relationship: Among them, in the 6th to 11th contour lines, in the same type, the nth contour line is the parent contour line of the n+1th contour line, that is, the n+1th contour line is the child contour line of the nth contour line. In the 1st to 5th contour lines, start from the last contour line of the same type and search from the back to the front. Let the current contour line be a, and search forward in sequence. When encountering the contour line b, it can simultaneously meet the following conditions:

[0134] minx(a)<=minx(b)miny(a)<=miny(b),maxx(a)>=maxx(b),maxy(a)>=maxy(b), then b is the parent contour of a;

[0135] S6.4. Construct a closed area with contour lines. Contour lines and boundaries form a closed area. Add corresponding boundary points according to different contour line types.

[0136] The construction of closed area is divided into a non-closed area algorithm and a closed area algorithm;

[0137] The algorithm steps of the non-closed area are as follows:

[0138] S7.1. In the boundary point sequence, let point pi(xi,yi) be an unused point and save its sequence number and coordinate value as the current point, with sequence number i.

[0139] S7.2, add the point pi(xi,yi) to the region and mark it as used, find the vertex pi+1(xi+1,yi+1) of the next region boundary line, record the point number i+1, and add it to the region;

[0140] S7.3, continue to search for the next point... When the starting point of the contour line is found, add the contour line to the region; when the end point of the contour line is found, add the contour line to the region in the reverse direction, and take the other end point qj (xj, yj) as the current point, and save the sequence number j;

[0141] S7.4, loop search until returning to the starting point pi(xi,yi), and finally achieving the purpose of forming a closed connected area, that is, a non-closed area;

[0142] S7.5, sequentially taking out the points pi(xi,yi) marked as unused in steps 2 and 3, repeating the above steps, and finally forming a number of non-closed regions;

[0143] S7.6, when all the points marked as unused have been used, the construction process is terminated and all the non-closed regions are obtained;

[0144] The steps of forming a closed area algorithm are as follows:

[0145] S8.1. Find the point with the smallest horizontal coordinate value in the closed contour line, record it as minpi(xi,yi), and save the sequence number i of this point;

[0146] S8.2. Sort the closed contour lines in the order of the size of the points minp1(x1,y1)…minpn(xn,yn) and save the sequence.

[0147] Preferably, the area coloring is divided into the following steps:

[0148] S9.1. Suppose a closed area has N nodes with coordinates (xi, yi), (i = 1, 2, 3 ... n), then the area calculation formula is:

[0149]

[0150] Calculate the area of ​​the region according to this formula and sort it by area size:

[0151] S9.2. When a region contains a sub-region, its fill parameter is the smaller value of the contour line and the contour line of the sub-region; when it does not contain a sub-region, it is the value of its own contour line;

[0152] S9.3. Select different colors according to the different values ​​of all regional boundaries. The same color can only fill the areas with the same rainfall value, and different rainfall values ​​are filled with different colors;

[0153] The labeling of the contour lines has the following steps:

[0154] S10.1, determine whether the contour line is a closed contour line, if yes, proceed to step 2, and determine whether the contour line has been marked, if no, proceed to the next step;

[0155] S10.2. Find three consecutive equal-value points on the contour line, set as Q1, Q2, Q3, and calculate the angle a between each two line segments between these three points. When the angle value is greater than a given value α, that is, a>α, where α is set to 130°, it is considered that the rainfall value can be marked at this location;

[0156] S10.3. Determine whether the point is a rainfall station. If it is, continue to search for the next point Q4 in the original direction. If it is not a rainfall station, continue to mark the point.

[0157] S10.4. When the isohyetal line is closed, the closed isohyetal line segment is disconnected and the marked value is added at the break. When the isohyetal line is not closed, the mark is added at both ends of the line.

[0158] Compared with the prior art, the present invention provides a method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine, which has the following beneficial effects:

[0159] 1. The present invention conducts a detailed study and implementation of the generation and filling algorithm of rainfall contour lines based on rectangular grids. Specifically, a two-person backtracking function is added to the traditional contour line tracking algorithm to make improvements, thereby solving the problem of waste of secondary traversal caused by the traditional method. The idea of ​​reverse tracking with a two-way queue is used to improve the resource consumption phenomenon in the tracking process. In the contour line smoothing process, a smoothing algorithm of multi-segment Bezier curve fitting is used to smooth the contour lines, making the smoothing process more flexible and the final smoothing effect more ideal. BRIEF DESCRIPTION OF THE DRAWINGS

[0160] Figure 1 A flow chart for generating contour lines using the rectangular grid method of the present invention;

[0161] Figure 2-1 shows the distribution of rainfall observation sites in some areas of Jiangsu Province;

[0162] Figure 2-2 is a diagram of a regular grid model;

[0163] Figure 2-3 Schematic diagram of interpolation of discrete data;

[0164] Figure 2-4 Projection point set and grid node diagram of the original point on plane xoy;

[0165] Figure 2-5 Kriging interpolation effect diagram;

[0166] Figure 2-6 is a, b pentagonal graph;

[0167] Figure 2-7 is a grid edge contour map;

[0168] Figure 3-1 is a diagram connecting equal value points a and b;

[0169] Figure 3-2 is a smoothing effect diagram of contour lines;

[0170] Figure 3-3 is a diagram showing the connection of two Bezier curves; Figure 3.3(a) is a diagram of isovalue points on the edge of a grid; Figure 3.3(b) shows a diagram of three equal-value points on the grid edge;

[0171] Figure 3-4 is a schematic diagram of the smoothed contour lines;

[0172] Figure 4-1 is a flow chart of the filling algorithm;

[0173] Figure 4-2 is a flow chart of the filling algorithm;

[0174] Figure 4-3 Schematic diagram of color ratio and filled contour lines;

[0175] Figure 4-4 Annotated diagram of contour lines;

[0176] Figure 5-1 System flow chart;

[0177] Figure 5-2 Contour map;

[0178] Figure 5-3 Contour map of rainfall in Jiangsu Province in February;

[0179] Figure 5-4 Rainfall isosurface map of Jiangsu Province in February. DETAILED DESCRIPTION

[0180] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0181] See also Figure 1 -5, including grid modeling of discrete rainfall data, generation algorithm of rainfall contour lines and filling algorithm of rainfall contour lines;

[0182] The grid modeling of discrete rainfall data includes a DEM model, a regular grid generated by discrete points, spatial interpolation of rainfall and Kriging interpolation, and the grid modeling of discrete rainfall data is built based on the DEM model, the regular grid generated by discrete points, spatial interpolation of rainfall and Kriging interpolation;

[0183] The generation algorithm of rainfall contour lines is based on the grid modeling of discrete rainfall data. The generation algorithm of rainfall contour lines includes contour lines, the generation and process of contour lines, rectangular grids and the process of realizing rectangular grids by contour line algorithms. The generation algorithm of rainfall contour lines is based on the theoretical basis of contour lines, the generation and process of contour lines, rectangular grids and the process of realizing rectangular grids by contour line algorithms.

[0184] The filling algorithm of rainfall contour lines is based on the generation algorithm of rainfall contour lines. The filling algorithm of rainfall contour lines includes filling algorithm design and contour line annotation. The filling algorithm of rainfall contour lines is based on the filling algorithm design and contour line annotation as the practical basis to generate the filling algorithm.

[0185] The rainfall in DEM is expressed in a basin distribution. The DEM model uses known three-dimensional coordinate points in an arbitrary coordinate system to perform relevant statistics on continuous ground. DEM is divided into three categories: one is a regular grid model, another is an irregular triangulated network model, and the other is a contour model. These three models can be converted and generated from each other.

[0186] The estimation of the distributed rainfall in the watershed requires the longitude and latitude information of each station, as well as the elevation information of each station;

[0187] DEM is a finite sequence of three-dimensional vectors representing area D:

[0188] {Vi=(xi,yi,zi),i=1,2,……,n}

[0189] Where Vi represents the value of the ground characteristic information corresponding to the pth grid; xi,yi represents the coordinate value of a data point on the xoy plane, and zi represents the elevation value corresponding to the data point (xi,yi).

[0190] The data source of the digital elevation model is obtained by interpreting aerial photos, collecting GPS receivers, importing elevation data into points and converting vector contour lines. The spatial storage of the digital elevation model is completed in the form of regular grids, triangular patches and contour lines. Each node has an x ​​and y coordinate value. The z value of the point determines the shape of the surface, and the break line defines the change of the surface.

[0191] The basic content of DEM is summarized by a set of ordered point pairs of a two-dimensional function:

[0192] Mp=fk(vp,wp)

[0193] k=1,2,3......m; p=1,2,3......n (2-1)

[0194] Wherein, Mp represents the value of the ground characteristic information corresponding to the pth grid point, vp, wp represent the two-dimensional plane coordinates of point p, and any plane coordinate system of map projection can be used, m is the number of ground characteristic information types, n is the number of ground points, Kp refers to a certain characteristic attribute of the ground, and the attributes include soil moisture content, vegetation coverage and soil type. The digital rainfall model is a special model when m=1 and the ground characteristic attribute represented by Mp is the rainfall value of each grid point. The present invention takes the rainfall distribution of some cities and counties in Jiangsu Province as an example, as shown in Figure 2-1;

[0195] The regular grid generated by the discrete points is a rectangular regular grid, the basic unit of the regular grid model is a square or rectangle, the spatial position of each grid unit is represented by the coordinates (x, y), and a value is assigned to each grid unit, and the value assigned to each grid unit of the regular grid generated by the discrete points is a rainfall value;

[0196] The modeling steps for generating regular grids from discrete points are as follows:

[0197] S1.1, first divide the data area into rectangular grids on a two-dimensional plane to form a plane grid covering the entire area;

[0198] S1.2. Using the sampling points distributed in the grid area, the rainfall values ​​of each vertex of the grid unit are obtained through a specific mathematical algorithm to obtain the final rectangular grid model, as shown in Figure 2-2.

[0199] Spatial interpolation of rainfall: The accurate estimation and inversion of spatial data are of great significance to the research and development of hydrology and water conservancy. Meteorological observation data is a kind of spatial data. The discrete monitoring data of a specific area is used to estimate the non-monitoring data on the regular grid points, so that the discrete station observation data can be converted into regular grid data. This is "spatial interpolation";

[0200] In order to draw rainfall contour lines using observation data from discrete rainfall stations, after completing the modeling of rainfall values, it is necessary to perform spatial interpolation on rainfall data, because the rainfall values ​​obtained from the stations are irregular and irregular data. The two-dimensional spatial coordinates of discrete data are generally determined by the plane coordinates of discrete sampling points, and the third dimension is the rainfall value. The distribution of stations in rainfall observation areas is usually uneven, so the data points obtained have uneven distribution characteristics. It is for this reason that the interpolation algorithm can be used to estimate the data in areas without stations. The spatial interpolation theory has a premise assumption, namely: in terms of spatial position, the closer the distance between points, the closer the characteristic values ​​and the higher the correlation; on the contrary, the farther the distance between points, the weaker the similarity and correlation. Based on this premise and assumption, rainfall data is related to ground elevation and the distance from known rainfall stations to a certain extent, as shown in Figure 2-3;

[0201] The spatial interpolation of rainfall uses the discrete monitoring data of a specific area to estimate the non-monitoring data on regular grid points, and converts the discrete station observation data into regular grid data. The spatial interpolation of rainfall is divided into global interpolation and local interpolation. Local interpolation includes: Thiessen polygon method, inverse distance interpolation, spline function interpolation method and Kriging interpolation; the spatial interpolation of rainfall includes the interpolation principle and process of rectangular grids;

[0202] The interpolation principle of the rectangular grid is as follows: For the interpolation of the rectangular grid, first establish a regular rectangular grid covering the point set Zi(xi,yi) on the xoy plane, where Zi(xi,yi) is the projection point set of Zi(xi,yi,zi) on the xoy plane. For each grid node, an insertion point must be defined. The defined insertion point satisfies the projection point on the xoy plane as the grid node. Find n points Zi(xi,yi) close to the grid node position. To realize interpolation, the process of selecting the measured data points and calculating the influencing parameters w is required;

[0203] The interpolation process of the rectangular grid is:

[0204] (1) Select the measurement data points

[0205] When performing spatial data interpolation, measurement data points are selected from the vicinity of the insertion point to participate in the calculation. The calculation process has the following steps:

[0206] S2.1, taking the grid node as the center, find the projection point in the nearby area;

[0207] S2.2, determine whether the projection point found on the xoy plane is a valid point. The valid point is the projection point of the measured data point on the xoy plane, as shown in Figure 2-4;

[0208] S2.3, find the corresponding measurement data point according to the found valid projection point;

[0209] (2) Calculation of influencing parameters

[0210] When calculating the rainfall value of an insertion point based on the found measurement data points, the farther the data point is from the insertion point, the smaller the influence. When it is far from the data point, the calculation of the influencing parameters is required;

[0211] Assume that the insertion point is Z(x,y,z), the data point sequence is Z1(x1,y1,z1)…Zn(xn,yn,zn), and the projection distance Di from the projection point of Z(x,y,z) to Zi(xi,yi,zi) is obtained by the following formula:

[0212] Di=(x-xi)*(x-xi)+(y-yi)*(y-yi) (2-2)

[0213] The impact factor wi of Zi is defined as

[0214] (3) Interpolation test standard.

[0215] The interpolation test standard uses the cross-validation method to verify the effect of interpolation. The cross-validation method assumes that the rainfall values ​​of some stations with known rainfall in the area are unknown variables, and then interpolates the rainfall values ​​of nearby stations. The error between the actual monitoring value and the theoretical estimated value of the station to be estimated is calculated, and this is used as a standard to judge the quality of the interpolation method.

[0216] The interpolation test standard uses average error, mean absolute error, mean relative error, mean square error and root mean square error as evaluation criteria for analysis. The algorithm is as follows:

[0217]

[0218] Experimental comparison of several common interpolation methods: Create a data table as the interpolation data source in xls format, convert xls to grid data format, and use this grid data to draw contour lines interpolated using different interpolation methods, as shown in Table 1:

[0219] Table 1 Annual rainfall at rain gauge stations in some areas of Jiangsu

[0220]

[0221]

[0222] The Kriging interpolation method first needs to analyze the spatial field, propose a variogram model, and then carry out Kriging calculations based on this variogram model. This method makes full use of the properties of the data space field and the correlation between data points. Therefore, in the interpolation process, it can not only reflect the anisotropy of the spatial field, but also automatically identify the spatial distribution of sampling points, eliminating the errors caused by uneven sampling. At the same time, the Kriging interpolation method can produce maximum and minimum values ​​beyond the estimated area, so the contour lines drawn are more accurate and in line with reality. The contour line map drawn by this algorithm is shown in Figure 2-5:

[0223] Kriging interpolation algorithm description:

[0224] for(int i = 0; i <Num;i++){myData[i].distance=(myData[i].x-importX)*(myData[i].x-importX)+(myData[i].y-importY)*(myData[i].y-importY);}

[0225] for(int i = 0; i<Num;i++){If((myData[i].x-importX> 0)&&(myData[i].y-importY>0)){myData[i].ID=1;}

[0226] If((myData[i].x-importX<0)&&(myData[i].y-importY>0)){myData[i].ID=2;twoQuadrant++;}

[0227] If((myData[i].x-importX<0)&&(myData[i].y-importY<0)){myData[i].ID=3; threeQuadrant++;}

[0228] If((myData[i].x-importX<0)&&(myData[i].y

[0229] -importY<0)){myData[i].ID=3; threeQuadrant++;}

[0230] If((myData[i].x-importX>0)&&(myData[i].y-importY<0)){myData[i].ID=4;}

[0231] The Kriging interpolation method first needs to analyze the spatial field, propose a variogram model, and then carry out Kriging calculation based on this variogram model;

[0232] The key factors of Kriging interpolation include regionalized variables and variograms. Spatial information with dual attributes of structure and uncertainty is called regionalized variables. Suppose Z(x) is a regionalized variable and satisfies the second-order stationary assumption, then the following formula is obtained:

[0233]

[0234] Among them, the determination of the weight coefficient λ is the core of the Kriging method. When solving the weight coefficient λ, it is necessary to fully ensure that the linear unbiasedness of the estimator Ze(u) and the minimum estimation variance are met.

[0235] Contour lines:

[0236] (1) An isoline is a smooth, continuous curve;

[0237] (2) For a given attribute value, the corresponding number of contour lines is the corresponding number;

[0238] (3) Contour lines can be closed or open;

[0239] (4) The contour lines will not intersect each other;

[0240] The principle and process of generating contour lines:

[0241] The data type for contour drawing is discrete two-dimensional data. The domain of discrete data is three-dimensional domain, and its data itself is a three-dimensional coordinate vector (x, y, z), where (x, y) is a plane vector and z is an attribute value used to represent the rainfall value. There are two main forms of discrete data distribution: one is a rectangular grid distribution data field, and the other is a scattered irregular distribution data field. The process of generating contour lines using the rectangular grid method is: first, a set of three-dimensional discrete data points (x, y, z) is interpolated by an algorithm to obtain the attribute values ​​of each rectangular grid point, and its equivalent points are obtained by linear interpolation on the boundary of the grid cell; then, the contour lines are tracked according to the tracking principle of the rectangular grid method and the contour line map is drawn. The generation process of the contour line is divided into three steps:

[0242] S3.1, Grid processing of discrete data;

[0243] S3.2, search for contour lines;

[0244] S3.3, display of contour map;

[0245] The rectangular grid method uses the collected discrete data points to establish a surface model on a regular area close to the geographic surface. The model uses an interpolation algorithm to interpolate the discrete points to obtain the attribute values ​​of each grid point. Then, according to the attribute values ​​of the adjacent points on the edge of each cell in the rectangular grid, the specific positions of all the equal-value points on the edge of each cell are obtained through linear interpolation. Finally, according to the tracking principle of the rectangular grid method, all the equal-value points of all the isolines are searched out and connected in sequence into smooth curves, namely: isolines;

[0246] The implementation of the contour line regular rectangular grid method includes the calculation of equal value points, the establishment of the connection order of equal value points and the tracking of contour lines. The process is as follows:

[0247] Step 1: The generation of contour lines first requires obtaining discrete data through corresponding algorithms based on actual rainfall data;

[0248] Step 2: Then, interpolate the discrete data to obtain gridded data;

[0249] Step 3: Calculate the gridded data to obtain equivalent points;

[0250] Step 4: Track and smooth the points with equal value, and finally generate the equal value lines;

[0251] For a flowchart of the rectangular grid method for generating contour lines, see Figure 1 ;

[0252] Calculation of equivalent points:

[0253] The drawing area obtained by rectangular grid modeling is composed of m×n grid units. The xoy plane coordinate system is established, the coordinate of the lower left corner of the rectangular area is (0,0), and the basic unit lengths of the directional grid units are dx and dy respectively; according to formula (3-5-1) and formula (3-5-2), the horizontal and vertical coordinate values ​​of any grid point can be obtained as follows:

[0254] X=i×dx (i=0,1,2……m) (3-5-1)

[0255] Y=j×dy(j=0,1,2...n)(3-5-2)

[0256] Among them, dx and dy represent the length of the grid unit in the X direction and the Y direction respectively, i and j represent the number of grid units in the X direction and the Y direction respectively, and the formula is deduced based on the horizontal side pentagon a (left side of Figure 2-6) and the vertical side pentagon b (right side of Figure 2-6) obtained by connecting the two-sided triangles. The a, b pentagon is shown in Figure 2-6, and the formula is;

[0257] Horizontal side pentagon a:

[0258]

[0259] Get the intersection coordinates

[0260] y=y A =y B (3-2)

[0261] Vertical side pentagon b:

[0262]

[0263] Get the intersection coordinates x=x A =x B (3-3)

[0264]

[0265] According to the above steps, the equivalent points on all grid edges can be calculated;

[0266] Establishment of the connection order of contour points: The attribute values on the edges of grid cells are linearly distributed. At this time, each edge of the grid cell has at most 4 intersections with the contour line and can have 2 intersections. At the same time, it is proved by the method of contradiction that there is no case of one or three intersections between the contour line and the edge of the grid cell, as shown in Figure 2-7:

[0267] Proof: Assume that the given attribute value is w. Let the attribute value of point A be wA, the attribute value of point B be wB, the attribute value of point C be wC, and the attribute value of point D be wD. Assume that wA>w.

[0268] ∵ wA>w, and there are contour points on the grid cell edge AD, ∴ wD<w

[0269] At this time, there may be the following 4 cases for points B and C:

[0270] (1) wC>w and wB>w, then there should be another contour point between the CD edges;

[0271] (2) wC>w and wB<w, then there is another contour point on each of the AB, CD, and BC edges;

[0272] (3) wC<w and wB<w, then there is another contour point between the AB edges;

[0273] (4) wC<w and wB>w, then there is another contour point between BC.

[0274] To sum up, it is contradictory to the fact that there is only one contour point in Figure 3.3(a), indicating that the situation of only one intersection does not exist.

[0275] Assume that wA>w. ∵ wA>w, and there are contour points on the grid cell edge AD ∴ wD<w

[0276] Also, ∵ there is one contour point between each of the edges AB, BC, and DA, ∴ wB<w, wD<w, wC>w;

[0277] ∵ wC>w and wD<w

[0278] ∴ There is another contour point on the CD edge, which is contradictory to the assumption. Therefore, the situation in Figure 3.3(b) is also impossible. Therefore, the number of contour points on a grid cell can only be 0, 2, or 4.

[0279] The two situations in Figures 3-1(a) and 3-1(b) will cause the contour lines of different rainfall values to intersect.

[0280] The tracing of the contour line is improved on the traditional method;

[0281] The tracing of the contour line includes the smoothing of the contour line, cubic spline function, tension spline function, and Bezier function;

[0282] The methods for smoothing contour lines include: contour line smoothing algorithms include tension spline function, linear iteration method, piecewise cubic polynomial method, Bezier curve function, B-spline method and cubic Bezier method. The contour line smoothing algorithms meet the following requirements: (1) forming a piecewise linear function; (2) forming a smooth curve through the fitted curve;

[0283] Cubic spline function:

[0284] Given a set of type value points Pi = (xi, yi), i = 0, 1, ..., n, construct two interpolating cubic spline functions with respect to the parameters:

[0285] x=x(u),y=y(u)

[0286] Interpolate at the point sets (u,x) and (u,y), i=0,1,2...n respectively, and then merge the three to form a cubic parametric spline curve:

[0287] P(u)=(x(u),y(y))

[0288] Determine contours by constructing parametric equations:

[0289]

[0290] The tension spline function is a linear combination of a piecewise interpolation function and two hyperbolic functions shσx and chσx:

[0291] f(x)=c1+c2x+c3shσx+c4chσx (3-6)

[0292] The whole has a continuous second-order derivative. The role of σ is to control the position of the inflection point and the shape of the curve. Through mathematical derivation, the function expression of the tension spline function is obtained:

[0293]

[0294] The Bezier function is composed of line segments and nodes. The nodes are draggable curves. The Bezier curve is defined as: given n+1 control vertices, its expression is:

[0295] P(t)=∑B i ,n(t)P i u∈[0,1] (3-8)

[0296] Where: Bi,n(t) is called basis function;

[0297] The harmonic function expression between the vertex position vectors of the Bezier polygon is:

[0298]

[0299] where: Bi,n(k) is the Bernstein basis function;

[0300] For the smooth condition of the multi-segment Bezier curve splicing, when using the Bezier curve fitting in segments, the following conditions need to be satisfied:

[0301] (1) The positions of two adjacent Bezier curves need to be continuous, and the start and end points of the front and back two curves should be connected;

[0302] (2) The slopes of two adjacent Bezier curves need to be continuous, that is, the starting point where the front and back two curves coincide should be collinear with two adjacent points;

[0303] (3) The curvatures of two adjacent Bezier curves need to be continuous, that is, the two curves should ensure that five points are coplanar.

[0304] Design of the multi-segment Bezier curve fitting algorithm. The algorithm process of using the multi-segment Bezier curve to fit the isocurve is as follows:

[0305] S4.1. Assume n + 1 discrete data points P0, P1,..., Pn on the isocurve. First, obtain two points P0 and P3 of the first segment of the Bezier curve. The initial point k for each segment to start fitting is 0 for the first segment. An initial step size is used to fit the number of discrete data points of a segment of the Bezier curve, and an error limit value E;

[0306] S4.2. An initial value of M is equal to one step;

[0307] S4.3. Take the points from the k-th to the k + m-th in the discrete data point sequence. When k + m > n, then m = n - k points, that is, Pk, Pk+1,..., Pk+m, and use the least squares method to approximate the isosegment with a Bezier curve;

[0308] S4.4. Calculate the error e between the points on the Bezier curve and the points on the corresponding isosegment;

[0309] When e < E, this segment of fitting is successful, and k = k + m; when k ≥ n, then end, otherwise continue to execute step 2; when e ≥ E, during this segment of fitting, by reducing the number of discrete data points for fitting m = 2*m / 3, continue to execute step 2 for fitting. The amount of discrete point data is reduced by 1 / 3. First, obtain two points P0 and P3 of the first segment of the Bezier curve.

[0310] Assume that the Bezier curve passes through n points on the polyline at ti, and the n points are S0, S1,..., Sn respectively. Then there is the formula:

[0311]

[0312] In the formula, let S0 = P0, S3 = P3, that is, find the two points P0 and P3 of the first section of the curve. After finding points P0 and P3, find points P1 and P2, and the formula is:

[0313]

[0314] Use the multi-segment idea to smooth the contour lines, as shown in Figure 3-3. From the smooth condition 1, we can get T0 = P3. Use the least squares method to find three control points, which are: T1, T2, T3. Let: T1 = T0 + G*R, where T0 = P3;

[0315] There is a formula:

[0316]

[0317] Take a vertex of the broken line as a double vertex, plus a point on the two preceding and following line segments to form the control polygon of the cubic Bezier curve, determine the start and end points of each cubic Bezier segment, and from the properties of the Bezier curve, it can be concluded that the start and end points of each cubic Bezier curve are the first and last points of the control polygon respectively.

[0318] Use Bezier curve fitting contour code to achieve: int i; float step = 1.0f / (float)(pnum-1), t, f1, f2, f3; float bx, by; for (i = 0; i <pnum;i++){t=i*step;f1=(float)0.5*(t-1)*(t-1);f2=(float)0.5*((-2)*t*t+2*t+1);f3=(float)0.5*t*t;bx=f1*ctlx[0]+f2*ctlx[1]

[0319] +f3*ctlx[2]; by=f1*ctly[0]+f2*ctly[1]+f3*ctly[2]; Bctlx.Add(bx); Bctly.Add(by);}……

[0320] for(i=0;i <src_vec_len-2;i++){tmctlx[0]=ctlx[i];tmctly[0]=ctly[i];tmctlx[1]=ctlx[i+1];

[0321] tmctly[1]=ctly[i+1]; tmctlx[2]=ctlx[i+2]; tmctly[2]=ctly[i+2]; tmBctlx.Clear(); tmBctly.Clear(); B_points(tmctlx,tmctly,pnum,tmBctlx,tmBctly); for(j=0;j <tmBctlx.Count;j++){Bctlx.Add(tmBctlx[j]);Bctly.Add(tmBctly[j]);}

[0322] The filling algorithm is designed to divide the filling work into three steps:

[0323] S5.1, construct connected regions;

[0324] S5.2, sort the regions in order of their area size;

[0325] S5.3, filling each area with color according to the determined filling parameters;

[0326] The filling algorithm divides the contour lines into 11 categories according to the relationship between the contour lines and the boundaries:

[0327] (a) The intersection of the contour lines as boundaries is 0;

[0328] (b) The contour line only intersects the lower boundary;

[0329] (c) The contour line only intersects the right boundary;

[0330] (d) The contour line only intersects the left boundary;

[0331] (e) The contour lines intersect only the upper boundary;

[0332] (f) The contour lines intersect the right and lower boundaries;

[0333] (g) The contour lines intersect the left and lower boundaries;

[0334] (h) The contour lines intersect the left and upper boundaries;

[0335] (i) The contour line intersects the right and upper boundaries;

[0336] (j) The contour line intersects the left and right boundaries;

[0337] (k) The contour lines intersect the upper and lower boundaries;

[0338] The filling algorithm flow is:

[0339] S6.1, initialize the iso-value area, and classify the iso-value lines according to the relationship between the iso-value lines and the boundaries;

[0340] S6.2. Sort the contour lines from large to small according to their classification. If the classification is the same, in order to ensure that the colors can be filled in order, they must be sorted according to certain rules. Assume that the coordinates of the first point of the contour line are x, y, and the maximum values ​​of the points above, below, left and right of the contour line are maxy, miny, maxx, maxy;

[0341] S6.3. Determine the parent-child relationship between contour lines of the same type. Judgment of the parent-child relationship: Among them, in the 6th to 11th contour lines, in the same type, the nth contour line is the parent contour line of the n+1th contour line, that is, the n+1th contour line is the child contour line of the nth contour line. In the 1st to 5th contour lines, start from the last contour line of the same type and search from the back to the front. Let the current contour line be a, and search forward in sequence. When encountering the contour line b, it can simultaneously meet the following conditions:

[0342] minx(a)<=minx(b)miny(a)<=miny(b),maxx(a)>=maxx(b),maxy(a)>=maxy(b), then b is the parent contour of a;

[0343] S6.4. Construct a closed area with contour lines. Contour lines and boundaries form a closed area. Add corresponding boundary points according to different contour line types.

[0344] Table 4.2 Sorting method for the same contour line type

[0345] type sort by Category 1 The ones with larger maxx-minx values ​​are placed first. Category 2 The ones with the largest maxy value are placed first Category 3 Miny values ​​with larger values ​​are placed in front Category 4 Minx values ​​with smaller values ​​are placed first Category 5 The ones with larger maxx values ​​are placed first Category 6 The ones with larger y values ​​are placed in front. Category 7 The ones with larger y values ​​are placed in front. Category 8 The smaller the y value, the higher the Category 9 The smaller the y value, the higher the Category 10 The ones with larger y values ​​are placed in front. Category 11 The x value with the largest value is placed first.

[0346] The construction of closed areas is divided into non-closed area algorithms and closed area algorithms;

[0347] The algorithm steps for non-closed areas are as follows:

[0348] S7.1. In the boundary point sequence, let point pi(xi,yi) be an unused point and save its sequence number and coordinate value as the current point, with sequence number i.

[0349] S7.2, add the point pi(xi,yi) to the region and mark it as used, find the vertex pi+1(xi+1,yi+1) of the next region boundary line, record the point number i+1, and add it to the region;

[0350] S7.3, continue to search for the next point... When the starting point of the contour line is found, add the contour line to the region; when the end point of the contour line is found, add the contour line to the region in the reverse direction, and take the other end point qj (xj, yj) as the current point, and save the sequence number j;

[0351] S7.4, loop search until returning to the starting point pi(xi,yi), and finally achieving the purpose of forming a closed connected area, that is, a non-closed area;

[0352] S7.5, sequentially taking out the points pi(xi,yi) marked as unused in steps 2 and 3, repeating the above steps, and finally forming a number of non-closed regions;

[0353] S7.6, when all the points marked as unused have been used, the construction process is terminated and all the non-closed regions are obtained;

[0354] The algorithm steps to form a closed area are as follows:

[0355] S8.1. Find the point with the smallest horizontal coordinate value in the closed contour line, record it as minpi(xi,yi), and save the sequence number i of this point;

[0356] S8.2. Arrange the closed contour lines in the order of the points minp1(x1,y1)…minpn(xn,yn);

[0357] The following steps are used to fill an area:

[0358] S9.1. Suppose a closed area has N nodes with coordinates (xi, yi), (i = 1, 2, 3 ... n), then the area calculation formula is:

[0359]

[0360] Calculate the area of ​​the region according to this formula and sort it by area size:

[0361] S9.2. When a region contains a sub-region, its fill parameter is the smaller value of the contour line and the contour line of the sub-region; when it does not contain a sub-region, it is the value of its own contour line;

[0362] S9.3. Select different colors according to the different values ​​of all regional boundaries. The same color can only fill the areas with the same rainfall value, and different rainfall values ​​are filled with different colors;

[0363] intcolor;PointF[]points;intfillingNum=0;for(inti=0;i<_closedContourLists.Count;i++){if(_cl osedContourLists[i].contourType==ContourType.ClosedContour){closedList.Add(_closedCont ourLists[i]); fillingNum++;}}SolidBrushbrush=newSolidBrush(Color.FromArgb(color));g.FillPolyg on(brush,points);}ScanLineClip_Calcu_Abs_Area_of 1 Polygon(AfterDotExpanddest_fxy z DzxCalcu,scal_w,scal_h,

[0364] out_all_city_bord_List[i].scrnx,out_all_city_bord_List[i].scrny,out_all_city_bord_List[i].scrnx.Count,areaeachgrid.areas_grid_from_top_to_bottum,

[0365] The following steps are involved in labeling contour lines:

[0366] S10.1, determine whether the contour line is a closed contour line, if yes, proceed to step 2, and determine whether the contour line has been marked, if no, proceed to the next step;

[0367] S10.2. Find three consecutive equal-value points on the contour line, set as Q1, Q2, Q3, and calculate the angle a between each two line segments between these three points. When the angle value is greater than a given value α, that is, a>α, where α is set to 130°, it is considered that the rainfall value can be marked at this location;

[0368] S10.3. Determine whether the point is a rainfall station. If it is, continue to search for the next point Q4 in the original direction. If it is not a rainfall station, continue to mark the point.

[0369] S10.4. When the isohyetal line is closed, the closed isohyetal line segment is disconnected and the marked value is added at the break. When the isohyetal line is not closed, the mark is added at both ends of the line.

[0370] Example 1

[0371] Based on the above algorithm principle, the algorithm is written in C# language, and the generation of contour lines is realized in the contour line subsystem of the Jiangsu Water Conservancy Information Service Platform, taking the rainfall in some areas of Jiangsu Province as an example. As one of the 15 subsystems of the Jiangsu Water Conservancy Geographic Information Service Platform software, the rainfall contour line production system calculates and generates contour lines, isosurfaces, and heat maps of water data such as rainfall, soil moisture, and groundwater according to the data preprocessing results, and provides a platform service interface to support other users to call.

[0372] 1. The main functions of the isobaric subsystem in the Jiangsu Provincial Water Conservancy Information Service Platform include:

[0373] (1) Data selection: Users request services to query the spatial location of sites by selecting by name, by fuzzy query, by spatial location, or by combined condition. The system then calls the data selection service of the platform rainfall contour service to implement the system's regional selection function and return the user's requested results. At the same time, users continue to request hydrological thematic attribute queries from the system to implement the time period selection function and extract the hydrological attribute data of the template site. After the system accepts the request, it calls the data query service and returns the result data such as time period precipitation, daily precipitation, and historical precipitation.

[0374] (2) Contour line drawing: The user requests the system to draw contour lines. The system calls the contour line drawing service provided by the platform rainfall contour line production, and calls the contour line generation service and the isosurface generation service according to the area selected by the user and the database selected for connection. After the contour line drawing is completed, the user continues to request symbolization of the generated results. The system calls the platform contour line drawing service to set the symbols, selects the color, width and line style of the generated contour lines, and returns the symbolized results. When the sites are densely populated, the user requests the system to provide a comprehensive map service. The system calls the multi-site stratification service in the platform contour line drawing service to use multi-site dynamic stratification to automatically display the data in layers according to the site level and return the results.

[0375] (3) Graphic output: The user requests the system to export result data. For the contour lines, isosurfaces, heat maps and other results generated by the system, the system calls the platform’s graphic output service to output the results in the form of vectors or XML or JSON.

[0376] 2. The module division of the contour line production subsystem includes two modules: data selection module and contour line drawing module;

[0378] Debugging phase of contour drawing module;

[0379] The contour drawing module calls the contour drawing algorithm through the service interface and defines the contour button as the contour drawing function:

[0380] clhjf1.Draw_Float_JWD_On_Bmp_Line(bord_jwdx,bord_jwdy,CDZXBord.Display_wid,CDZ XBord.Display_high,CDZXBord.Display_edge_shift); clhjf1.Draw_Float_JWD_On_Bmp_Line(bo rd_jwdx,bord_jwdy,"a1","a2",1200,1200,0);

[0381] The least squares control is used to test the fit:

[0382] public static double[]MultiLine(double[]arrX,double[]arrY,intlength,int dimension;){int n=dimension+1;double[,]Guass=new double[n,n+1];for(int i=0;i <n;i++){intj;

[0383] for(j=0; j <n;j++){Guass[i,j]=SumArr(arrX,j+i,length);}Guass[i,j]=SumArr(arrX,i,arrY,1,length);}return ComputGauss(Guass,n);}

[0384] The monitoring station names, longitude and latitude of each city and county, and daily rainfall are extracted. After the discrete data is rectangularly gridded, grid data is obtained. After equivalent point calculation and tracking, the contour map obtained after debugging and operation is shown in Figure 5-2:

[0385] 3. Operating Environment of Contour Line Drawing Subsystem

[0386] (1) Hardware equipment

[0387] The network equipment of the platform is configured, and the hardware equipment requirements of the platform are shown in Table 3 below:

[0388] Table 3 Hardware equipment list

[0389]

[0390] (2) Software Environment

[0391] The software environment is shown in Table 4 below:

[0392] Table 4 Software environment requirements

[0393] Software Type Software version requirements operating system Windows Server 2008 R2 Enterprise Edition x64 Database Systems Oracle10g

[0394] 4. The operating effect diagrams of the contour drawing subsystem are shown in Figures 5-3 and 5-4.

[0395] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine, characterized in that: Including grid modeling of discrete rainfall data, generation algorithm of rainfall contour lines and filling algorithm of rainfall contour lines; The grid modeling of discrete rainfall data includes a DEM model, a regular grid generated by discrete points, spatial interpolation of rainfall and Kriging interpolation, and the grid modeling of discrete rainfall data is built based on the DEM model, the regular grid generated by discrete points, spatial interpolation of rainfall and Kriging interpolation; The generation algorithm of the rainfall contour line is based on the grid modeling of discrete rainfall data. The generation algorithm of the rainfall contour line includes contour lines, the generation and process of contour lines, rectangular grids and the process of realizing rectangular grids by the contour line algorithm. The generation algorithm of the rainfall contour line is based on the theoretical basis of contour lines, the generation and process of contour lines, rectangular grids and the process of realizing rectangular grids by the contour line algorithm. The filling algorithm of the rainfall contour line is based on the generation algorithm of the rainfall contour line. The filling algorithm of the rainfall contour line includes the filling algorithm design and the annotation of the contour line. The filling algorithm of the rainfall contour line is based on the filling algorithm design and the annotation of the contour line as the practical basis to generate the filling algorithm.

2. The method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine according to claim 1, characterized in that: The rainfall of the DEM is expressed in a basin distribution. The DEM model uses known three-dimensional coordinate points in an arbitrary coordinate system to perform relevant statistics on continuous ground. DEM is divided into three categories: one is a regular grid model, another is an irregular triangulated network model, and the other is a contour line model. These three models can be converted and generated from each other. The basin-distributed rainfall estimation requires the longitude and latitude information of each station, as well as the elevation information of each station; DEM is a finite sequence of three-dimensional vectors representing area D: {Vi=(xi,yi,zi),i=1,2,……,n} Where Vi represents the value of the ground characteristic information corresponding to the pth grid; xi,yi represents the coordinate value of a data point on the xoy plane, and zi represents the elevation value corresponding to the data point (xi,yi).

3. The method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine according to claim 1, characterized in that: The data source of the digital elevation model is obtained by interpreting aerial photos, collecting GPS receivers, importing elevation data points and converting vector contour lines. The spatial storage of the digital elevation model is completed in the form of regular grids, triangular patches and contour lines, each node of which has an x ​​and y coordinate value. The z value of the point determines the shape of the surface, and the broken line defines the change of the surface. The basic content of the DEM is summarized by a set of ordered point pairs of a two-dimensional function: Mp=fk(vp,wp) k=1,2,3......m; p=1,2,3......n (2-1) Wherein, Mp represents the value of the ground characteristic information corresponding to the pth grid point, vp, wp represent the two-dimensional plane coordinates of point p, and any plane coordinate system of map projection can be used, m is the number of ground characteristic information types, n is the number of ground points, Kp refers to a characteristic attribute of the ground, including soil moisture content, vegetation coverage and soil type. The digital rainfall model is a special model when m = 1 and the ground characteristic attribute represented by Mp is the rainfall value of each grid point; The discrete point generated regular grid is a rectangular regular grid, the basic unit of the regular grid model is a square or rectangle, the spatial position of each grid unit is represented by the coordinate (x, y), and a value is assigned to each grid unit, and the discrete point generated regular grid is assigned a rainfall value; The modeling steps of generating a regular grid from discrete points are as follows: S1.1, first divide the data area into rectangular grids on a two-dimensional plane to form a plane grid covering the entire area; S1.

2. Using the sampling points distributed in the grid area, the rainfall value of each vertex of the grid unit is obtained through a specific mathematical algorithm, thereby obtaining the final rectangular grid model.

4. The method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine according to claim 1, characterized in that: The spatial interpolation of rainfall utilizes the discrete monitoring data of a specific area to estimate the non-monitoring data on a regular grid point, and converts the discrete site observation data into regular grid data. The spatial interpolation of rainfall is divided into global interpolation and local interpolation. The local interpolation includes: Thiessen polygon method, inverse distance interpolation, spline function interpolation method and Kriging interpolation; the spatial interpolation of rainfall includes the interpolation principle of a rectangular grid and the interpolation process of a rectangular grid; The interpolation principle of the rectangular grid is as follows: for the interpolation of the rectangular grid, a regular rectangular grid covering the point set Zi(xi,yi) is first established on the xoy plane, where Zi(xi,yi) is the projection point set of Zi(xi,yi,zi) on the xoy plane. An insertion point is defined for each grid node, and the defined insertion point satisfies the projection point on the xoy plane as the grid node. Find n points Zi(xi,yi) close to the grid node position, and realize the interpolation through the process of selecting the measured data points and calculating the influencing parameters w; The interpolation process of the rectangular grid is: (1) Select the measurement data points When performing spatial data interpolation, measurement data points are selected from the vicinity of the insertion point to participate in the calculation. The calculation process has the following steps: S2.1, taking the grid node as the center, find the projection point in the nearby area; S2.2, determine whether the projection point found on the xoy plane is a valid point, and the valid point is the projection point of the measured data point on the xoy plane; S2.3, find the corresponding measurement data point according to the found valid projection point; (2) Calculation of influencing parameters When calculating the rainfall value of an insertion point based on the found measurement data points, the farther the data point is from the insertion point, the smaller the influence. When it is far from the data point, the calculation of the influencing parameters is required; Assume that the insertion point is Z(x,y,z), the data point sequence is Z1(x1,y1,z1)…Zn(xn,yn,zn), and the projection distance Di from the projection point of Z(x,y,z) to Zi(xi,yi,zi) is obtained by the following formula: Di=(x-xi)*(x-xi)+(y-yi)*(y-yi) (2-2) The impact factor wi of Zi is defined as (3) Interpolation test standard.

5. The method and system for generating and filling rectangular grid rainfall contour lines based on a digital simulation engine according to claim 1, characterized in that: The interpolation test standard adopts the cross-validation method to verify the effect of interpolation. The cross-validation method is to assume that the rainfall values ​​of some stations of known rainfall stations in the area are unknown variables, and then interpolate the rainfall values ​​of nearby stations to make an estimate, and then calculate the error between the actual monitoring value and the theoretical estimated value of the station to be estimated, and use this as a standard to judge the quality of the interpolation method; The interpolation test standard uses average error, average absolute error, average relative error, mean square error and root mean square error as evaluation criteria for analysis, and the algorithm is as follows: Experimental comparison of several common interpolation methods: create a data table as the interpolation data source, the data table format is xls, convert xls into grid data format, and use this grid data to draw contour lines interpolated by different interpolation methods; The Kriging interpolation method first needs to analyze the spatial field, propose a variogram model, and then carry out Kriging calculation based on the variogram model; The key factors of the Kriging interpolation include regionalized variables and variograms. The spatial information with dual attributes of structure and uncertainty is called regionalized variables. If Z(x) is a regionalized variable and satisfies the second-order stationary assumption, the following formula is obtained: Among them, the determination of the weight coefficient λ is the core of the Kriging method. When solving the weight coefficient λ, it is necessary to fully ensure that the linear unbiasedness of the estimator Ze(u) and the minimum estimation variance are met.

6. The method and system for generating and filling rainfall contour lines of rectangular grids based on a digital simulation engine according to claim 1, characterized in that: The contour lines: (1) An isoline is a smooth, continuous curve; (2) For a given attribute value, the corresponding number of contour lines is the corresponding number; (3) Contour lines can be closed or open; (4) The contour lines will not intersect each other; The generation principle and process of the contour lines are as follows: The data type for contour drawing is discrete two-dimensional data. The domain of discrete data is three-dimensional domain, and its data itself is a three-dimensional coordinate vector (x, y, z), where (x, y) is a plane vector and z is an attribute value used to represent the rainfall value. There are two main forms of discrete data distribution: one is a rectangular grid distribution data field, and the other is a scattered irregular distribution data field. The process of generating contour lines using the rectangular grid method is: first, a set of three-dimensional discrete data points (x, y, z) is interpolated by an algorithm to obtain the attribute values ​​of each rectangular grid point, and its equivalent points are obtained by linear interpolation on the boundary of the grid cell; then, the contour lines are tracked according to the tracking principle of the rectangular grid method and the contour line map is drawn. The generation process of the contour line is divided into three steps: S3.1, Grid processing of discrete data; S3.2, search for contour lines; S3.3, display of contour map; The rectangular grid method uses the collected discrete data points to establish a surface model on a regular area close to the geographic surface. The discrete points in the obtained model are interpolated by an interpolation algorithm to obtain the attribute value of each grid point. Then, according to the attribute values ​​of the adjacent points on the edge of each cell in the rectangular grid, the specific positions of all the equal-value points on the edge of each cell are obtained by linear interpolation. Finally, according to the tracking principle of the rectangular grid method, all the equal-value points of all the isolines are searched out and connected in sequence into smooth curves, namely: isolines; The implementation of the contour line regular rectangular grid method includes the calculation of equal value points, the establishment of the connection order of equal value points and the tracking of contour lines, and its process is as follows: Step 1: The generation of contour lines first requires obtaining discrete data through corresponding algorithms based on actual rainfall data; Step 2: Then, interpolate the discrete data to obtain gridded data; Step 3: Calculate the gridded data to obtain equivalent points; Step 4: Track and smooth the points with equal value, and finally generate the equal value lines; Calculation of the equivalent points: The drawing area obtained by rectangular grid modeling is composed of m×n grid units. The xoy plane coordinate system is established, the coordinate of the lower left corner of the rectangular area is (0,0), and the basic unit lengths of the directional grid units are dx and dy respectively; according to formula (3-5-1) and formula (3-5-2), the horizontal and vertical coordinate values ​​of any grid point can be obtained as follows: X=i×dx(i=0,1,2...m)(3-5-1) Y=j×dy(j=0,1,2...n)(3-5-2) Among them, dx and dy represent the length of the grid unit in the X direction and the Y direction respectively, i and j represent the number of grid units in the X direction and the Y direction respectively, and the formula is calculated based on the horizontal side pentagon and the vertical side pentagon obtained by connecting the two-sided triangles, and the formula is: Horizontal side pentagon: Get the intersection coordinates y=y A =and B (3-2) Vertical pentagon Get the intersection coordinates x=x A =x B (3-3) According to the above steps, the equivalent points on all grid edges can be calculated; The establishment of the connection order of the equivalent points: the attribute values ​​on the grid unit edges are linearly distributed. At this time, each grid unit edge has a maximum of 4 intersections with the equivalent line, and can have 2 intersections. At the same time, the proof by contradiction is used to prove that there is no one or three intersections between the equivalent line and the grid unit edge.

7. The method and system for generating and filling rectangular grid rainfall contour lines based on a digital simulation engine according to claim 1, characterized in that: The contour line tracking includes contour line smoothing, cubic spline function, tension spline function and Bezier function; The method for smoothing the contour line includes: the contour line smoothing algorithm includes tension spline function, linear iteration method, piecewise cubic polynomial method, Bezier curve function, B-spline method and cubic Bezier method, and the contour line smoothing algorithm satisfies: (1) forming a piecewise linear function; (2) forming a smooth curve through the fitted curve; The cubic spline function: Given a set of type value points Pi = (xi, yi), i = 0, 1, ..., n, construct two interpolating cubic spline functions with respect to the parameters: x=x(u),y=y(u) Interpolate at the point sets (u,x) and (u,y), i=0,1,2...n respectively, and then merge the three to form a cubic parametric spline curve: P(u) = (x(u), y(y)) Determine the isocurve by constructing parametric equations: The tension spline function is a piecewise interpolation function, which is a linear combination of linear interpolation and two hyperbolic functions shσx and chσx: f(x) = c1 + c2x + c3shσx + c4chσx (3-6) It has a continuous second derivative as a whole. The role of σ is to control the position of the inflection point and the shape of the curve. Through mathematical derivation, the function expression of the tension spline function is obtained: The Bezier function is composed of line segments and nodes. The nodes are curves that can be dragged. The Bezier curve is defined as: Given n + 1 control vertices, its expression is: P(t)=∑B i ,n(t)P i u∈[0,1] (3-8) Where: Bi,n(t) is called the basis function; The harmonic function expression between the position vectors of the vertices of the Bezier polygon is: Where: Bi,n(k) is the Bernstein basis function; For the smooth condition of splicing multiple Bezier curves, when using Bezier curves for piecewise fitting, the following conditions need to be met: (1) The positions of two adjacent Bezier curves need to be continuous, and the first and last points of the front and back two curves need to be connected; (2) The slopes of two adjacent Bezier curves need to be continuous, that is, the starting points where the front and back two curves coincide need to be collinear with two adjacent points; (3) The curvatures of two adjacent Bezier curves need to be continuous, that is, the two curves need to ensure that five points are coplanar.

8. The method and system for generating and filling rectangular grid rainfall contour lines based on a digital simulation engine according to claim 1, characterized in that: For the design of the algorithm for fitting multiple Bezier curves, the algorithm process of using multiple Bezier curves to fit the isocurve is as follows: S4.

1. Assume n + 1 discrete data points P0, P1, ……, Pn on the isocurve. First, obtain two points P0 and P3 of the first Bezier curve. The initial point k for each segment to start fitting is 0 for the first segment. An initial step size is used to fit the number of discrete data points of a Bezier curve segment, and an error limit value E. S4.

2. An initial value of M is equal to one step; S4.

3. Take the k to k + m points in the discrete data point sequence. When k + m > n, then m = n - k points, that is, Pk, Pk+1 …… Pk+m are used to approximate the isosegment with a Bezier curve by the least squares method; S4.

4. Calculate the error e between the points on the Bezier curve and the corresponding points on the isosegment; When e < E, this segment is successfully fitted, and k = k + m; when k ≥ n, it ends, otherwise continue to execute step 2; when e ≥ E, during this segment fitting, by reducing the number of discrete data points for fitting m = 2*m / 3, continue to execute step 2 for fitting. The amount of discrete point data is reduced by 1 / 3. First, obtain two points P0 and P3 of the first Bezier curve. Assume that the Bezier curve passes through n points on the broken line at ti. The n points are S0, S1, …, Sn respectively, then there is a formula: Where, let S0 = P0, S3 = P3, that is, obtain two points P0 and P3 of the first segment curve. After obtaining points P0 and P3, to obtain points P1 and P2, there is a formula: Use the multi-segment idea to smooth the contour lines. According to the smooth condition 1, we can get T0 = P3. Use the least square method to find three control points, which are: T1, T2, T3. Let: T1 = T0 + G*R, where T0 = P3; There is a formula: Take a vertex of the broken line as a double vertex, plus a point on the two preceding and following line segments to form the control polygon of the cubic Bezier curve, determine the start and end points of each cubic Bezier segment, and from the properties of the Bezier curve, it can be concluded that the start and end points of each cubic Bezier curve are the first and last points of the control polygon respectively.

9. The method and system for generating and filling rectangular grid rainfall contour lines based on a digital simulation engine according to claim 1, characterized in that: The filling algorithm design divides the filling work into three steps: S5.1, construct connected regions; S5.

2. Sort by area size; S5.3, filling each area with color according to the determined filling parameters; The filling algorithm divides the contour lines into 11 categories according to the relationship between the contour lines and the boundaries: (a) The intersection of the contour lines as boundaries is 0; (b) The contour line only intersects the lower boundary; (c) The contour line only intersects the right boundary; (d) The contour line only intersects the left boundary; (e) The contour lines intersect only the upper boundary; (f) The contour lines intersect the right and lower boundaries; (g) The contour lines intersect the left and lower boundaries; (h) The contour lines intersect the left and upper boundaries; (i) The contour line intersects the right and upper boundaries; (j) The contour line intersects the left and right boundaries; (k) The contour lines intersect the upper and lower boundaries; The filling algorithm process is as follows: S6.

1. Initialize the iso-value area and classify the iso-value lines according to the relationship between the iso-value lines and the boundaries; S6.

2. Sort the contour lines from large to small according to their classification. If the classification is the same, in order to ensure that the colors can be filled in order, they must be sorted according to certain rules. Assume that the coordinates of the first point of the contour line are x, y, and the maximum values ​​of the points above, below, left and right of the contour line are maxy, miny, maxx, maxy; S6.

3. Determine the parent-child relationship between contour lines of the same type. Judgment of the parent-child relationship: Among them, in the 6th to 11th contour lines, in the same type, the nth contour line is the parent contour line of the n+1th contour line, that is, the n+1th contour line is the child contour line of the nth contour line. In the 1st to 5th contour lines, start from the last contour line of the same type and search from the back to the front. Let the current contour line be a, and search forward in sequence. When encountering the contour line b, it can simultaneously meet the following conditions: minx(a)<=minx(b)miny(a)<=miny(b),maxx(a)>=maxx(b),maxy(a)>=maxy(b), then b is the parent contour of a; S6.

4. Construct a closed area with contour lines. Contour lines and boundaries form a closed area. Add corresponding boundary points according to different contour line types. The construction of closed area is divided into a non-closed area algorithm and a closed area algorithm; The algorithm steps of the non-closed area are as follows: S7.

1. In the boundary point sequence, let point pi(xi,yi) be an unused point and save its sequence number and coordinate value as the current point, with sequence number i. S7.2, add the point pi(xi,yi) to the region and mark it as used, find the vertex pi+1(xi+1,yi+1) of the next region boundary line, record the point number i+1, and add it to the region; S7.3, continue to search for the next point... When the starting point of the contour line is found, add the contour line to the region; when the end point of the contour line is found, add the contour line to the region in the reverse direction, and take the other end point qj (xj, yj) as the current point, and save the sequence number j; S7.4, loop search until returning to the starting point pi(xi,yi), and finally achieving the purpose of forming a closed connected area, that is, a non-closed area; S7.5, sequentially taking out the points pi(xi,yi) marked as unused in steps 2 and 3, repeating the above steps, and finally forming a number of non-closed regions; S7.6, when all the points marked as unused have been used, the construction process is terminated and all the non-closed regions are obtained; The steps of forming a closed area algorithm are as follows: S8.

1. Find the point with the smallest horizontal coordinate value in the closed contour line, record it as minpi(xi,yi), and save the sequence number i of this point; S8.

2. Sort the closed contour lines in the order of the size of the points minp1(x1,y1)…minpn(xn,yn) and save the sequence.

10. The method and system for generating and filling rectangular grid rainfall contour lines based on a digital simulation engine according to claim 1, characterized in that: The area filling is divided into the following steps: S9.

1. Suppose a closed area has N nodes with coordinates (xi, yi), (i = 1, 2, 3 ... n), then the area calculation formula is: Calculate the area of ​​the region according to this formula and sort it by area size: S9.

2. When a region contains a sub-region, its fill parameter is the smaller value of the contour line and the contour line of the sub-region; when it does not contain a sub-region, it is the value of its own contour line; S9.

3. Select different colors according to the different values ​​of all regional boundaries. The same color can only fill the areas with the same rainfall value, and different rainfall values ​​are filled with different colors; The following steps are involved in labeling the contour lines: S10.1, determine whether the contour line is a closed contour line, if yes, proceed to step 2, and determine whether the contour line has been marked, if no, proceed to the next step; S10.

2. Find three consecutive equal-value points on the contour line, set as Q1, Q2, Q3, and calculate the angle a between each two line segments between these three points. When the angle value is greater than a given value α, that is, a>α, where α is set to 130°, it is considered that the rainfall value can be marked at this location; S10.

3. Determine whether the point is a rainfall station. If it is, continue to search for the next point Q4 in the original direction. If it is not a rainfall station, continue to mark the point. S10.

4. When the isohyetal line is closed, the closed isohyetal line segment is disconnected and the marked value is added at the break. When the isohyetal line is not closed, the mark is added at both ends of the line.

Images