A method for generating a spatial grid based on geometric constraints

Abstract

The application discloses a kind of space grid generation methods based on geometric constraint, specifically includes: obtaining target area space range data, execution longitude and latitude analysis generates benchmark grid set;Construct Cartesian coordinate system, convert longitude and latitude into spatial coordinate set;Based on preset side length, execution recursive quadric dissection generates multistage equidistant grid cell set;With GeoSOT the 15th level grid establishes hierarchical mapping and forms fusion grid structure;Establishes encoding structure and generates grid encoding sequence;Establishes geographical entity mapping relationship and generates grid set;Based on grid set, execution distance calculation, area calculation and intersection operation generates topological relationship result.The present application realizes the equidistant dissection and fusion encoding calculation of space grid, significantly reduces the complexity of space distance and area calculation and improves the calculation accuracy and processing stability.

Description

Technical Field

[0001] This invention relates to the technical field of spatial grid coding and spatial computing, and in particular to a spatial grid generation method based on geometric constraints. Background Technology

[0002] With the development of spatial information technology and digital applications, spatial data is widely used in surveying and mapping, remote sensing, geographic information systems, and digital twins. To achieve efficient organization, management, and computation of spatial data, it is typically necessary to discretize continuous space into a grid structure to support operations such as spatial positioning, distance calculation, area calculation, and topological relationship analysis. Existing technologies widely employ global grid systems based on latitude and longitude, such as the GeoSOT grid. These methods achieve multi-level spatial partitioning through recursive quadtree partitioning, possessing a unified coding structure and multi-scale expression capabilities, and can meet the needs of large-scale spatial data management.

[0003] However, in practical applications, grid systems based on latitude and longitude coordinates are affected by the Earth's curvature and projection methods, resulting in inconsistent grid scales across different latitude regions. This leads to differences in the actual spatial size of grid cells. When performing distance and area calculations, complex projection transformations and curvature corrections are often required, increasing computational complexity and power consumption. In high-precision calculation scenarios for local areas, this type of method struggles to balance computational efficiency and accuracy, especially in large-scale spatial data processing, where computational performance is significantly limited.

[0004] Furthermore, existing spatial computing methods typically perform calculations based on latitude and longitude vector coordinates, lacking a unified computing mechanism based on discrete grids, which hinders efficient data indexing and batch computing operations. Differences in grid partitioning rules and encoding methods between different systems also increase the difficulty of spatial data fusion and sharing. Therefore, how to construct a spatial grid generation and computing method that balances local isometry with global grid system compatibility has become an urgent technical problem to be solved. Summary of the Invention

[0005] One objective of this invention is to propose a spatial grid generation method based on geometric constraints. This invention introduces a local equidistant grid partitioning mechanism and a grid coding analytical calculation model to perform grid positioning and coding-driven calculations on spatial data of the target region. It constructs a spatial distance calculation and surface accumulation calculation process based on grid coordinates, forming a unified spatial calculation structure with the characteristics of high computational efficiency, strong result consistency, and applicability to large-scale spatial data processing.

[0006] A spatial mesh generation method based on geometric constraints according to an embodiment of the present invention includes the following steps: Acquire spatial range data of the target area, perform latitude and longitude analysis calculations, extract boundary coordinate values, and generate a reference grid set; Using the coordinates of the lower left corner of the baseline grid as a reference point, a regional Cartesian coordinate system is constructed, and the latitude and longitude values ​​of each point in the target area are converted into Cartesian coordinate values ​​to form a set of spatial coordinates. In the Cartesian coordinate system, a recursive quadrilateral subdivision operation is performed on the spatial coordinate set using a preset side length as the 0th level grid scale. The row number and column number are recorded at each level of subdivision to generate a multi-level equidistant grid cell set. Using GeoSOT level 15 mesh as the reference unit, a local partitioning origin is established at the reference unit location. Mapping calculations are performed on the equidistant mesh level and the reference mesh level to form a fused mesh structure. Establish an encoding structure at the grid cell position, write hexadecimal encoding for the 0th level cell, write quaternary encoding for each level of subdivision, record column offset and row offset values ​​at the encoding position, and generate a grid encoding sequence. Establish geographic entity mapping relationships at the encoding location, perform grid coverage parsing on the geographic entities, and generate a grid set; Perform calculations at the locations of the grid set, calculate the distance using Euclidean distance to the grid center coordinates, calculate the area using surface accumulation to generate the area value, and perform intersection and union calculations on the grid set to generate topological relationship results.

[0007] Optionally, the latitude and longitude resolution calculation specifically includes: Establish a boundary point index sequence for the spatial range data of the target area, extract longitude and latitude values ​​at the index positions, write the longitude and latitude component record table at the value positions, perform minimum and maximum value calculations at the longitude value positions, and perform minimum and maximum value calculations at the latitude value positions to generate longitude range values ​​and latitude range values. Perform integer alignment operations on the longitude range values ​​to generate longitude boundary values, perform integer alignment operations on the latitude range values ​​to generate latitude boundary values, construct an outer rectangular range between the longitude boundary values ​​and the latitude boundary values, and generate a set of region boundary coordinates. Perform a grid matching operation at the location of the region boundary coordinate set to map the region boundary coordinate set to the GeoSOT level 15 grid. Generate a reference grid index sequence at the mapping location, perform a deduplication operation at the index location, retain a unique grid identifier, and form a reference grid set.

[0008] Optionally, the construction of the regional Cartesian coordinate system includes the following steps: Extract the longitude and latitude values ​​corresponding to the lower left corner grid cell at the reference grid set location, establish the reference point coordinates at this location, perform offset calculations to generate the origin of the region partitioning at the reference point location in the west and south directions, and write the origin identifier at the origin point location. The coordinate axis directions are defined by establishing the origin of the region partitioning, establishing the X-axis in the horizontal direction and the Y-axis in the vertical direction, and writing the direction labels at the coordinate axis positions to form a Cartesian coordinate system structure; The latitude and longitude component record table reads the longitude and latitude values, performs geodetic coordinate transformation at the value location, maps the latitude and longitude values ​​to plane coordinate values, and generates the corresponding X coordinate value and Y coordinate value at the mapping location. The coordinate values ​​are used to perform a range determination operation on the X and Y coordinate values. If the position meets the condition of the first quadrant, a valid identifier is written and the set of spatial coordinates of the region is output.

[0009] Optionally, the quadrilateral partitioning operation specifically includes: Establish a grid index sequence at the location of the regional spatial coordinate set, extract the X and Y coordinate values ​​corresponding to each point at the index location, set the initial grid side length value at the value location as the 0th level grid scale, and construct the 0th level grid coverage area at the scale location. Perform a binary partitioning calculation at the coverage area of ​​the level 0 grid, dividing the current grid cell into the lower left sub-cell, lower right sub-cell, upper left sub-cell, and upper right sub-cell. Write the column index value and row index value at the sub-cell position respectively, and form a sub-cell identifier set at the index position. Establish a hierarchical index identifier at the grid index sequence position, record the current subdivision hierarchical value at the hierarchical position, perform coordinate range calculation on each sub-unit at the hierarchical position, write the sub-unit boundary coordinate value at the range position, perform an inclusion determination operation on the sub-unit boundary coordinate value and the regional spatial coordinate set at the index position, and retain the sub-unit identifier at the position that meets the inclusion condition. Recursive partitioning is performed at the reserved sub-cell identifier position. At the recursive position, partitioning calculation and range determination calculation are repeatedly performed on each sub-cell. At each level of partitioning position, column index values ​​and row index values ​​are written into the hierarchical record table in hierarchical order to form a multi-level grid hierarchical index sequence. Perform edge length update calculation on the grid cell at the hierarchical record table position, perform binary search calculation on the current edge length value at the update position, write the next level edge length value at the calculation position, and establish the correspondence between grid cell and edge length value at the mapping position. Execute a termination decision operation at the multi-level grid level index sequence position. At the decision position, compare the current level value with the preset maximum level value. Stop the recursive subdivision at the position that meets the condition. Output the final set of grid cell identifiers at the termination position. Perform equidistant verification calculations at the final grid cell identifier set location, perform consistency judgment on the side length values ​​of adjacent grid cells at the verification location, and output multi-level equidistant grid cell sets at the location that meets the consistency conditions.

[0010] Optionally, the step of performing mapping calculations on the equidistant grid levels and the reference grid level to form a fused grid structure includes the following steps: Extract the GeoSOT Level 15 grid cell identifier from the reference grid set location, read the latitude and longitude range values ​​of the corresponding grid cell at the identifier location, write the reference grid boundary coordinate values ​​at the range location, and establish the reference grid index sequence at the coordinate location. Read the boundary coordinate values ​​of each equidistant grid cell at the location of the regional spatial coordinate set, extract the coordinate values ​​of the lower left corner of the grid at the boundary location, perform coordinate inverse calculation at the coordinate location, convert the planar coordinate values ​​into latitude and longitude values, and generate the corresponding latitude and longitude coordinate values ​​at the conversion location. Perform grid affiliation determination operation at the latitude and longitude coordinate value location. At the determination location, perform interval matching calculation between the converted latitude and longitude coordinate values ​​and the reference grid boundary coordinate values. Write the reference grid identifier at the location that meets the interval inclusion condition. Form the reference grid affiliation sequence at the equidistant grid cell location. Extract the level values ​​corresponding to each level grid cell at the position of the equidistant grid level index sequence, read the grid side length value at the level position, read the spatial scale value corresponding to the reference grid at the position of the reference grid index sequence, perform the proportional relationship calculation at the value position, write the level mapping coefficient value at the calculation position, and establish the scale correspondence between the equidistant grid level and the reference grid level at the mapping position. At the mapping relationship location, perform an association write operation on the reference grid identifier and the equidistant grid cell identifier, establish a mapping index table at the write location, and record the reference grid identifier, equidistant grid level value, and equidistant grid cell identifier at the index table location to form a set of grid mapping relationships; Perform a consistency check operation at the location of the grid mapping relationship set, perform coverage detection on the equidistant grid cells corresponding to the same reference grid at the check location, determine whether each equidistant grid cell completely falls into the corresponding reference grid range at the detection location, and write a valid identifier at the location that meets the conditions. At the valid identification location, perform a structural reorganization operation on the mapping relationship set, perform aggregation calculation on the equidistant grid cells according to the reference grid identification at the reorganization location, and generate a fused grid structure data set at the aggregation location.

[0011] Optionally, the generation of the grid-coded sequence specifically includes: Establish an coded index sequence at the grid cell index position, extract the 0th level grid cell identifier at the index position, convert the grid cell number value to a hexadecimal character value at the identifier position, write the coded prefix value at the character position, and form the initial coded field at the coded position. Extract the column index value and row index value corresponding to each level grid cell at the position of the equidistant grid level index sequence. Perform the sub-cell position determination operation at the index position according to the reverse Z order rule. Generate 4-ary code value at the determination position based on the column index value and row index value. Write the code field in the order of the level at the code position to form a level code sequence at the field position. At the encoding field position, perform concatenation calculation on the hexadecimal encoded prefix value and the hierarchical encoded sequence, generate a complete grid encoded sequence at the concatenation position, perform a uniqueness determination operation at the encoding position, and output the grid encoded sequence at the position that satisfies the uniqueness condition.

[0012] Optionally, the step of performing Euclidean distance calculation on the grid center coordinates to generate distance values ​​includes the following steps: Extract the target grid code value at the grid code sequence position. Perform hexadecimal parsing operation on the code prefix at the code position to generate the column number and row number value of the level 0 grid. Perform bit-by-bit parsing operation on the 4-ary code sequence at the code position. Extract the column offset value and row offset value according to the reverse Z order rule at the parsing position. Perform recursive calculation on the column number value and row number value at the parsing position to generate the final column number value and final row number value corresponding to the target grid cell. Combine the grid side length value with the column number position and row number position to perform coordinate calculation to generate the coordinate value of the lower left corner of the grid. At the coordinate position, the center offset calculation is performed between the coordinate value of the lower left corner of the grid and the grid side length to generate the center coordinate value of the grid. At the center coordinate position, the center coordinate values ​​corresponding to the two grid cells are extracted. At the coordinate position, the difference calculation is performed on the corresponding X coordinate value and Y coordinate value respectively. At the difference position, the square calculation is performed on the difference value. At the square position, the sum calculation is performed. At the calculation position, the square root operation is performed on the sum result to generate the distance value.

[0013] Optionally, the step of performing surface accumulation and calculation on the mesh set to generate area values ​​specifically includes: Extract the grid code values ​​at the grid set index position, perform hexadecimal prefix parsing operation at the code position to generate initial column number and row number values, perform bit-by-bit parsing operation on the 4-ary code sequence at the code position, extract column offset and row offset values ​​according to the reverse Z-order rule at the parsing position, perform recursive calculation on the column number and row number values ​​at the parsing position to generate final column number and final row number values, perform coordinate calculation on the column number and row number positions combined with the grid side length value to generate the coordinate value of the lower left corner of the grid, and construct the coordinate values ​​of the four vertices of the grid at the coordinate position; At the grid set index position, a coverage determination operation is performed on each grid cell. At the determination position, a spatial inclusion detection is performed on the coordinate values ​​of the four vertices of the grid. A complete grid identifier is written at the position where all vertices meet the target area range condition. An edge grid identifier is written at the position where some vertices meet the condition. At the complete grid identifier position, an area calculation is performed according to the grid side length value to generate a complete area value. At the edge grid identifier position, a subdivision calculation is performed on the grid cell. At the subdivision position, the grid cell is divided into several sub-cells. The center coordinate value is extracted at the sub-cell position. A spatial inclusion detection is performed at the coordinate position. At the detection position, the number of sub-cells that meet the condition is counted. At the count position, the coverage ratio value is calculated. At the ratio position, the area is calculated by combining the grid side length value to generate an edge area value. At the area position, the complete area value and the edge area value are accumulated to generate a total area value.

[0014] The beneficial effects of this invention are: (1) The present invention constructs a local equidistant grid system based on the Cartesian coordinate system, transforming the spatial calculation process from latitude and longitude vector calculation to grid coordinate calculation. In the process of distance calculation and area calculation, the coordinate transformation and curvature correction steps are reduced, the calculation complexity is reduced, and the efficiency of spatial data processing is improved.

[0015] (2) By establishing a hierarchical mapping relationship between the reference grid and the local equidistant grid, this invention achieves a unified expression of the global grid system and the local computational grid, forming a consistent coding and computational structure in the process of spatial positioning, distance calculation and topological relationship analysis, thereby enhancing the stability and applicability of spatial data processing. Attached Figure Description

[0016] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is a flowchart of a spatial mesh generation method based on geometric constraints proposed in this invention; Figure 2 This is a schematic diagram of local equidistant mesh generation for a spatial mesh generation method based on geometric constraints proposed in this invention; Figure 3 This is a schematic diagram of the GeoSOT fusion mechanism of a spatial mesh generation method based on geometric constraints proposed in this invention.

[0017] Figure 4 This is a schematic diagram of the level 0 computational grid encoding rule for a spatial grid generation method based on geometric constraints proposed in this invention. Detailed Implementation

[0018] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0019] refer to Figures 1-3 A spatial mesh generation method based on geometric constraints includes the following steps: Acquire spatial range data of the target area, perform latitude and longitude analysis calculations, extract boundary coordinate values, and generate a reference grid set; Using the coordinates of the lower left corner of the baseline grid as a reference point, a regional Cartesian coordinate system is constructed, and the latitude and longitude values ​​of each point in the target area are converted into Cartesian coordinate values ​​to form a set of spatial coordinates. In the Cartesian coordinate system, a recursive quadrilateral subdivision operation is performed on the spatial coordinate set using a preset side length as the 0th level grid scale. The row number and column number are recorded at each level of subdivision to generate a multi-level equidistant grid cell set. Using GeoSOT level 15 mesh as the reference unit, a local partitioning origin is established at the reference unit location. Mapping calculations are performed on the equidistant mesh level and the reference mesh level to form a fused mesh structure. Establish an encoding structure at the grid cell position, write hexadecimal encoding for the 0th level cell, write quaternary encoding for each level of subdivision, record column offset and row offset values ​​at the encoding position, and generate a grid encoding sequence. Establish geographic entity mapping relationships at the encoding location, perform grid coverage parsing on the geographic entities, and generate a grid set; Perform calculations at the locations of the grid set, calculate the distance using Euclidean distance to the grid center coordinates, calculate the area using surface accumulation to generate the area value, and perform intersection and union calculations on the grid set to generate topological relationship results.

[0020] In this embodiment, performing latitude and longitude resolution calculations specifically includes: Establish a boundary point index sequence for the spatial range data of the target area, extract longitude and latitude values ​​at the index positions, write the longitude and latitude component record table at the value positions, perform minimum and maximum value calculations at the longitude value positions, and perform minimum and maximum value calculations at the latitude value positions to generate longitude range values ​​and latitude range values. Perform integer alignment operations on the longitude range values ​​to generate longitude boundary values, perform integer alignment operations on the latitude range values ​​to generate latitude boundary values, construct an outer rectangular range between the longitude boundary values ​​and the latitude boundary values, and generate a set of region boundary coordinates. Perform a grid matching operation at the location of the region boundary coordinate set to map the region boundary coordinate set to the GeoSOT level 15 grid. Generate a reference grid index sequence at the mapping location, perform a deduplication operation at the index location, retain a unique grid identifier, and form a reference grid set.

[0021] Specifically, performing integer alignment operations to generate longitude boundary values ​​within a longitude range includes: Read the minimum and maximum longitude values ​​at the longitude range numerical positions, perform a degree splitting operation on the longitude values ​​at the numerical positions to decompose the longitude values ​​into degree, minute and second structures, and extract the component values ​​at the component positions. At the minimum longitude value position, the component value is rounded down. At the rounding position, the component value is adjusted to an integer minute value. At the minute position, the second value is set to zero. At the calculation position, the degree value, minute value and second value are recombined to generate the lower longitude boundary value. At the position of the maximum longitude value, the component values ​​are rounded up. At the rounding position, if the minute value has a decimal part or the second value is greater than zero, the minute value is incremented by one. At the minute position, the second value is set to zero. When the minute value meets the carry condition, the degree value is incremented by one. At the calculation position, the degree value, minute value and second value are recombined to generate the upper boundary value of longitude. A range write operation is performed at the positions of the lower and upper longitude boundaries to form whole-point aligned longitude boundary values.

[0022] Specifically, performing mesh matching operations at the coordinate set location of the region boundary includes: Read the longitude and latitude values ​​corresponding to each boundary point at the region boundary coordinate set location, establish an index sequence at the boundary point location, read the GeoSOT level 15 grid division parameters at the index location, extract the longitude and latitude division interval values ​​at the parameter location, perform difference calculation and division operation on the longitude values ​​at the boundary point location, perform rounding operation to generate longitude direction column number values ​​at the result location, perform the same operation to generate latitude direction row number values ​​at the latitude value location, construct grid identifiers at the column number location and row number location, form a grid index sequence at the identifier location, perform traversal calculation between the minimum and maximum column number values ​​at the index sequence location, perform traversal calculation between the minimum and maximum row number values, generate a grid identifier set for the covered area at the traversal location, perform deduplication judgment operation at the set location, and generate a baseline grid set at the output location.

[0023] In this embodiment, constructing the regional Cartesian coordinate system includes the following steps: Extract the longitude and latitude values ​​corresponding to the lower left corner grid cell at the reference grid set location, establish the reference point coordinates at this location, perform offset calculations in the west and south directions at the reference point location to generate the origin of the region partitioning, and write the origin identifier at the origin point location. The offset calculation is to perform a fixed distance translation operation at the reference point location along the direction of decreasing longitude and decreasing latitude to generate the origin of the region partitioning. The coordinate axis directions are defined by establishing the origin of the region partitioning, establishing the X-axis in the horizontal direction and the Y-axis in the vertical direction, and writing the direction labels at the coordinate axis positions to form a Cartesian coordinate system structure; The latitude and longitude component record table reads the longitude and latitude values ​​at the location, performs geodetic coordinate transformation operation at the value location, maps the latitude and longitude values ​​to plane coordinate values, and generates the corresponding X coordinate value and Y coordinate value at the mapping location. The geodetic coordinate transformation operation is the process of converting the latitude and longitude coordinate values ​​to plane rectangular coordinate values ​​at the value location. The coordinate values ​​are used to perform a range determination operation on the X and Y coordinate values. If the position meets the condition of the first quadrant, a valid identifier is written and the set of spatial coordinates of the region is output.

[0024] The geodetic coordinate transformation involves reading longitude and latitude values ​​at the specified locations, extracting the longitude and latitude reference values ​​corresponding to the origin of the regional subdivision at the specified locations, performing a difference calculation between the longitude values ​​and the longitude reference values ​​at the longitude value locations to generate longitude difference values, performing a difference calculation between the latitude values ​​and the latitude reference values ​​at the latitude value locations to generate latitude difference values, multiplying the difference values ​​with preset distance conversion coefficients at the difference locations to generate corresponding planar distance values, writing the longitude distance value into the X-coordinate value at the calculation location, writing the latitude distance value into the Y-coordinate value at the calculation location, and forming a planar rectangular coordinate representation at the coordinate location.

[0025] The specific calculations for determining the range of X and Y coordinate values ​​include: Read the minimum boundary values ​​in the X and Y directions corresponding to the origin of the region partitioning at the coordinate position, and read the maximum boundary values ​​in the X and Y directions corresponding to the spatial range of the region at the boundary position. At the coordinate value position, perform difference calculation between the X coordinate value and the minimum boundary value in the X direction, and determine whether the calculation result is greater than or equal to zero at the difference position. At the coordinate value position, perform difference calculation between the X coordinate value and the maximum boundary value in the X direction, and determine whether the calculation result is less than or equal to zero at the difference position. At the coordinate value position, perform difference calculation between the Y coordinate value and the minimum boundary value in the Y direction, and determine whether the calculation result is greater than or equal to zero at the difference position. At the coordinate value position, perform difference calculation between the Y coordinate value and the maximum boundary value in the Y direction, and determine whether the calculation result is less than or equal to zero at the difference position. Perform a combined judgment operation on each judgment result at the judgment position, write a valid flag at the position that meets the condition, write an invalid flag at the position that does not meet the condition, and form a range judgment result at the output position.

[0026] In this embodiment, the quadrilateral partitioning operation specifically includes: Establish a grid index sequence at the location of the regional spatial coordinate set, extract the X and Y coordinate values ​​corresponding to each point at the index location, set the initial grid side length value at the value location as the 0th level grid scale, and construct the 0th level grid coverage area at the scale location. Perform a binary partitioning calculation at the coverage area of ​​the level 0 grid, dividing the current grid cell into the lower left sub-cell, lower right sub-cell, upper left sub-cell, and upper right sub-cell. Write the column index value and row index value at the sub-cell position respectively, and form a sub-cell identifier set at the index position. Establish a hierarchical index identifier at the grid index sequence position, record the current subdivision hierarchical value at the hierarchical position, perform coordinate range calculation on each sub-unit at the hierarchical position, write the sub-unit boundary coordinate value at the range position, perform an inclusion determination operation on the sub-unit boundary coordinate value and the regional spatial coordinate set at the index position, and retain the sub-unit identifier at the position that meets the inclusion condition. Recursive partitioning is performed at the reserved sub-cell identifier position. At the recursive position, partitioning calculation and range determination calculation are repeatedly performed on each sub-cell. At each level of partitioning position, column index values ​​and row index values ​​are written into the hierarchical record table in hierarchical order to form a multi-level grid hierarchical index sequence. The hierarchical record table performs edge length update calculations on the grid cells at the update position, performs binary search calculations on the current edge length value at the update position, writes the next level edge length value at the calculation position, and establishes the correspondence between grid cells and edge length values ​​at the mapping position. The multi-level grid hierarchy index sequence position performs a termination judgment operation. At the judgment position, the current level value is compared with the preset maximum level value. The recursive subdivision stops at the position that meets the condition, and the final set of grid cell identifiers is output at the termination position. Perform equidistant verification calculations at the final grid cell identifier set location, perform consistency judgment on the side length values ​​of adjacent grid cells at the verification location, and output multi-level equidistant grid cell sets at the location that meets the consistency conditions.

[0027] Specifically, the execution of the decision operation includes: Read the minimum X coordinate value, maximum X coordinate value, minimum Y coordinate value and maximum Y coordinate value of the sub-unit at the sub-unit boundary coordinate position, and form the sub-unit boundary range data at the coordinate position; Establish a coordinate point index sequence at the location of the regional spatial coordinate set, extract the X coordinate value and Y coordinate value corresponding to each coordinate point one by one at the index location, perform difference calculation on the X coordinate value and the minimum X coordinate value at the value location, determine whether the calculation result is greater than or equal to zero at the difference location, perform difference calculation on the X coordinate value and the maximum X coordinate value at the value location, and determine whether the calculation result is less than or equal to zero at the difference location. At the numerical position, perform a difference calculation between the Y coordinate value and the minimum Y coordinate value, and at the difference position, determine whether the calculation result is greater than or equal to zero; at the numerical position, perform a difference calculation between the Y coordinate value and the maximum Y coordinate value, and at the difference position, determine whether the calculation result is less than or equal to zero. At the decision position, perform a combined decision operation on each decision result, write an inclusion flag at the position that meets the condition, write a non-inclusion flag at the position that does not meet the condition, and form a sub-unit inclusion decision result at the output position.

[0028] Specifically, performing recursive subdivision operations at the reserved sub-cell identifier positions includes: Establish a sub-unit index sequence at the reserved sub-unit identifier position, extract the boundary coordinate values ​​corresponding to each sub-unit at the index position, read the current sub-unit side length value at the coordinate position, and record the current section level value at the value position. The sub-cell index position performs a binary partitioning calculation on the current sub-cell. At the partitioning position, the coordinates of the lower left corner of the sub-cell are used as a reference point to divide the current sub-cell into the lower left sub-cell, the lower right sub-cell, the upper left sub-cell, and the upper right sub-cell. The column index value and the row index value are written at the sub-cell position respectively, and a new layer of sub-cell identifier set is formed at the index position. The new layer of sub-unit identifier sets the location for each sub-unit to perform boundary coordinate calculations, update the minimum X coordinate value, maximum X coordinate value, minimum Y coordinate value and maximum Y coordinate value at the coordinate location, and generate the sub-unit boundary coordinate value at the coordinate location. Perform an inclusion decision operation at the sub-unit identifier set location, filter sub-unit identifiers that meet the regional spatial range conditions at the decision location, and form the next level of recursive sub-unit set at the retention location; At the level position, perform incremental calculation on the current level value. At the decision position, compare the incremented level value with the preset maximum subdivision level value. At the position where the termination condition is not met, repeat the partitioning and filtering calculation on the next level recursive subunit set. At the position where the termination condition is met, stop the recursive operation. At the output position, generate the final subunit set.

[0029] Specifically, the calculation of edge length update for grid cells based on the location of the hierarchical record table includes: Read the level value corresponding to the current grid cell in the level record table, read the side length value of the level 0 grid in the value position, and establish the side length value sequence in the value position. At each level position, the side length is calculated in hierarchical order. At the calculation position, the side length value of the previous level is divided by 2. At the result position, the side length value corresponding to the current level is written. At each level position, the corresponding side length value is generated level by level. Extract the grid cell identifier at the grid cell index position, read the corresponding level value at the identifier position, read the side length value corresponding to the level at the value position, and write the side length value into the grid cell record table at the corresponding position. In the grid cell record table, the side length values ​​of the corresponding grid cells at the same level are compared and calculated. A valid identifier is written at the position where the values ​​are consistent, and the updated grid side length value is generated at the output position.

[0030] Specifically, performing consistency checks includes: Establish a grid cell index sequence at the final grid cell identifier set location, extract the level value corresponding to each grid cell at the index location, read the grid side length value corresponding to the level at the value location, and form a grid side length value sequence at the record location. At the grid cell index sequence position, perform comparison calculation on adjacent grid cells, read the side length values ​​corresponding to the two grid cells at the comparison position, perform difference calculation at the value position, determine whether the difference value is equal to zero or within the preset error range at the difference position, write a consistency flag at the position that meets the condition, and write an inconsistency flag at the position that does not meet the condition. The grid cell index sequence position performs item-by-item comparison calculation on the side length values ​​corresponding to all grid cells in the same level. At the comparison position, the difference calculation is performed on each side length value. At the difference position, it is determined whether the difference value meets the consistency condition. At the position where all conditions are met, a level consistency flag is written. At the position where the conditions are not met, a level inconsistency flag is written. The determination position performs a combined determination calculation on the consistency identifier and the hierarchical consistency identifier. If the condition is met, an equidistant identifier is written; if the condition is not met, a non-equidistant identifier is written. An equidistant verification result is generated at the output position.

[0031] In this embodiment, performing mapping calculations between the equidistant grid level and the reference grid level to form a fused grid structure includes the following steps: Extract the GeoSOT Level 15 grid cell identifier from the reference grid set location, read the latitude and longitude range values ​​of the corresponding grid cell at the identifier location, write the reference grid boundary coordinate values ​​at the range location, and establish the reference grid index sequence at the coordinate location. Read the boundary coordinate values ​​of each equidistant grid cell at the location of the regional spatial coordinate set, extract the coordinate values ​​of the lower left corner of the grid at the boundary location, perform coordinate inverse calculation at the coordinate location, convert the planar coordinate values ​​into latitude and longitude values, and generate the corresponding latitude and longitude coordinate values ​​at the conversion location. Perform grid affiliation determination operation at the latitude and longitude coordinate value location. At the determination location, perform interval matching calculation between the converted latitude and longitude coordinate values ​​and the reference grid boundary coordinate values. Write the reference grid identifier at the location that meets the interval inclusion condition. Form the reference grid affiliation sequence at the equidistant grid cell location. Extract the level values ​​corresponding to each level grid cell at the position of the equidistant grid level index sequence, read the grid side length value at the level position, read the spatial scale value corresponding to the reference grid at the position of the reference grid index sequence, perform the proportional relationship calculation at the value position, write the level mapping coefficient value at the calculation position, and establish the scale correspondence between the equidistant grid level and the reference grid level at the mapping position. At the mapping relationship location, perform an association write operation on the reference grid identifier and the equidistant grid cell identifier, establish a mapping index table at the write location, and record the reference grid identifier, equidistant grid level value, and equidistant grid cell identifier at the index table location to form a set of grid mapping relationships; Perform a consistency check operation at the location of the grid mapping relationship set, perform coverage detection on the equidistant grid cells corresponding to the same reference grid at the check location, determine whether each equidistant grid cell completely falls into the corresponding reference grid range at the detection location, and write a valid identifier at the location that meets the conditions. At the valid identification location, perform a structural reorganization operation on the mapping relationship set, perform aggregation calculation on the equidistant grid cells according to the reference grid identification at the reorganization location, and generate a fused grid structure data set at the aggregation location.

[0032] Specifically, performing coordinate inverse calculations at coordinate positions includes: Read the X and Y coordinate values ​​corresponding to the grid cell at the coordinate position, read the X and Y reference values ​​corresponding to the origin of the region partition at the value position, perform addition calculation on the X coordinate value and X reference value at the coordinate position, generate the longitude distance value at the calculation position, perform addition calculation on the Y coordinate value and Y reference value at the coordinate position, and generate the latitude distance value at the calculation position. Read the preset distance conversion coefficient at the distance value location, perform division calculation on the distance value in the longitude direction at the value location to generate the longitude difference value, and perform division calculation on the distance value in the latitude direction at the value location to generate the latitude difference value. The system reads the longitude and latitude reference values ​​corresponding to the origin of the region subdivision at the difference value position. It then performs an addition calculation on the longitude difference value and the longitude reference value at the value position to generate the longitude value. Finally, it performs an addition calculation on the latitude difference value and the latitude reference value at the value position to generate the latitude value. The system then outputs the latitude and longitude coordinate results.

[0033] Specifically, performing proportional calculations at numerical positions includes: Read the level value corresponding to the equidistant grid at the numerical location, read the grid side length value corresponding to the level at the numerical location, read the spatial scale value corresponding to the reference grid at the numerical location, and form a scale value sequence at the recording location. At the numerical location, a division calculation is performed between the grid side length value and the reference grid spatial scale value. A proportional value is generated at the calculation location, and the proportional calculation result is recorded at the numerical location. Perform a range determination operation on the proportional value at the proportional value position, determine whether the proportional value meets the preset range conditions at the determination position, write a valid identifier at the position where the conditions are met, and form the proportional relationship calculation result at the output position.

[0034] Specifically, performing a structure reorganization operation on the mapping relationship set at the valid identifier location includes: Read the base grid identifier value, equidistant grid level value and equidistant grid cell identifier value corresponding to each record at the mapping relationship set location, and establish a record index sequence at the value location; At the record index sequence position, perform classification calculation on the base grid identifier value, divide the records corresponding to the same base grid identifier into the same group at the classification position, and form several record subsets at the grouping position; Perform sequential calculations on the equidistant grid cell identifier values ​​at each record subset location, sort them according to the level value, column number value, and row number value at the sorting location, and generate an ordered grid cell sequence at the sorting location. At the position of the ordered grid cell sequence, a continuity detection calculation is performed on the adjacent grid cells. At the detection position, the adjacency relationship of the grid cells in spatial coordinates is determined. At the position that meets the adjacency condition, a merging calculation is performed. At the calculation position, a combined grid cell set is generated. At the location of the combined grid cell set, a summary calculation is performed on the results of each group. At the summary location, a fused grid structure data set indexed by the reference grid identifier is generated, and at the output location, the structure reorganization result is formed.

[0035] Specifically, performing aggregation calculations on equidistant grid cells at the reorganization location according to the baseline grid identifier includes: Read the reference grid identifier value and the equidistant grid cell identifier value corresponding to each record at the mapping relationship set location, and establish a record index sequence at the value location; At the record index sequence position, perform classification calculation on the reference grid identifier value, divide the equidistant grid cell identifiers corresponding to the same reference grid identifier into the same group at the classification position, and form multiple equidistant grid cell sets at the grouping position; Extract the row and column numbers of the corresponding grid cells at each equidistant grid cell set location, perform sorting calculations on the row and column numbers at the value locations, and generate an ordered grid cell sequence at the sorting locations. At the ordered grid cell sequence position, the adjacent relationship detection calculation is performed on the adjacent grid cells. At the detection position, it is determined whether the difference between the row number value and the difference between the column number value meet the adjacent condition. At the position where the condition is met, the merging calculation is performed. At the calculation position, the adjacent grid cells are identified and combined into a continuous grid block. Perform summary calculations on each grid block at the continuous grid block location, generate an aggregated grid set under the corresponding reference grid identifier at the summary location, and form the aggregated calculation result at the output location.

[0036] In this embodiment, generating the grid-coded sequence specifically includes: Establish an coded index sequence at the grid cell index position, extract the 0th level grid cell identifier at the index position, convert the grid cell number value to a hexadecimal character value at the identifier position, write the coded prefix value at the character position, and form the initial coded field at the coded position. Extract the column index value and row index value corresponding to each level grid cell at the position of the equidistant grid level index sequence. Perform the sub-cell position determination operation at the index position according to the reverse Z order rule. Generate 4-ary code value at the determination position based on the column index value and row index value. Write the code field in the order of the level at the code position to form a level code sequence at the field position. At the encoding field position, perform concatenation calculation on the hexadecimal encoded prefix value and the hierarchical encoded sequence, generate a complete grid encoded sequence at the concatenation position, perform a uniqueness determination operation at the encoding position, and output the grid encoded sequence at the position that satisfies the uniqueness condition.

[0037] Specifically, performing a uniqueness determination operation at the encoding position includes: performing a step-by-step comparison calculation on the grid encoding values ​​at the encoding position, deleting codes with the same values, and forming a unique encoding sequence at the output position.

[0038] In this embodiment, the process of calculating the Euclidean distance to the grid center coordinates to generate a distance value includes the following steps: Step 1: Grid Location and Encoding Analysis By matching the latitude and longitude range of a geographic entity to its corresponding baseline grid (15th level GeoSOT grid), and then parsing the corresponding computational grid code (1 hexadecimal prefix + 16 quaternary main body), all grid cells covered by the entity are determined.

[0039] Step 2: Extracting grid coordinates (1) Let the side length of the level 0 computational grid (the initial, undivided grid) be... =1024m (default value in this invention); (2) Parse the hexadecimal prefix to determine the initial column and row numbers of the 0th level sub-unit. Convert the hexadecimal characters (0-F) into 4-bit binary numbers (padded with leading zeros if less than 4 bits). The first two bits of the binary number correspond to the "column number". (Range 0-3, from left to right column 0 to column 3), the last 2 bits of the binary number correspond to the "row number". (Range 0-3, from bottom to top, row 0 → row 3), the initial bottom left corner coordinates of the 0th level sub-unit are... ; (3) Parse the 4-ary main body, calculating the column and row offsets of each level of subdivision bit by bit. Each bit of the 4-ary main body (16 bits in total) corresponds to one reverse Z-order quadrature. It needs to be parsed bit by bit from left to right. The code value of each bit (0-3) corresponds to the "column offset + row offset" of the sub-unit (following the reverse Z-order rule: 0=bottom left, 1=bottom right, 2=top left, 3=top right). For the i-th bit of the 4-ary main body (i=1 to 16), let its code value be vi, then the column offset of this level of subdivision is (Corresponding to the column offsets in the table below), the row offset for this level of partitioning is... (Corresponding to the row offset in the table below), for each bit parsed, the column number and row number of the current grid need to be multiplied by 2 (because quad partitioning divides the current grid into 2 columns × 2 rows), plus the offset of that bit;

[0040] (4) Calculate the final column and row numbers, and derive the coordinates of the lower left corner of the grid. Combine the initial column and row numbers of the level 0 sub-cell with the offsets of each level of the 4-ary main body to obtain the final grid column number C and row number R. ; ; Finally, the coordinates of the lower left corner of the grid cell

[0041] ; ; ( This is the side length corresponding to the current grid level. If only the first m digits of the 4-ary main body are parsed, then... ;

[0042] Step 3: Distance Calculation Logic The distance between two points is taken as the center coordinates of the grids containing the two points. +L / 2,y n +L / 2) (where L is the side length of the grid cell) is calculated using the Euclidean distance formula, and the distance is... ; To calculate the distance between polygon features, it is necessary to calculate the centroid coordinates of the mesh sets covering the two polygon features. ; ; Sᵢ represents the area of ​​a single grid cell, and the centroid spacing is calculated using the Euclidean distance formula.

[0043] In this embodiment, the specific steps for generating area values ​​by performing surface accumulation and calculation on the mesh set include: Step 1: Regional Gridding (1) Using the origin (0,0) of the region as the reference, the latitude and longitude (lon,lat) of all vertices of the surface feature are converted into regional Cartesian coordinates (X,Y) (unit: meters) through the geodetic coordinate transformation formula. ; ; in:( , (where ) represents the latitude and longitude of the origin of the region subdivision, and k is the conversion coefficient between latitude / longitude and distance (near the equator, k≈111319.9m / °). (2) Extract the bounding rectangle of the surface feature based on the converted Cartesian coordinates. ; (3) By extracting grid coordinates (inverse Z-order encoding parsing), determine the encoding of all computational grids within the bounding rectangle and the coordinates of the lower left corner. , forming a grid set (k is the total number of grids covered by the circumscribed rectangle).

[0044] Step 2: Determining and counting complete meshes (1) For each grid Gi in the grid set G, the Cartesian coordinates of its four vertices are: If all four vertices are inside the surface feature or on its boundary, then it is considered a complete mesh. (2) Use the "ray method" to verify whether the vertex is in the plane: emit a ray from the vertex in the positive X-axis direction and count the number of intersections between the ray and the boundary of the surface feature; if the number of intersections is odd, the vertex is in the plane; if it is even, it is outside the plane; if the ray coincides with the boundary, it is directly determined to be on the boundary and included in the complete mesh. (3) Count the number of all complete grids that meet the conditions. The area of ​​a single complete grid is ; Therefore, the total area of ​​the complete grid is ; Step 3: Determining the edge mesh and calculating the coverage ratio (1) In the mesh set G, after excluding complete meshes, if the bounding rectangle of a mesh intersects with the bounding rectangle of a surface feature, and does not meet the "complete mesh condition", it is determined to be an edge mesh, and the number is counted. ; (2) Subdivide the edge grid Gi into n×n tiny subgrids (e.g., n=100, subgrid side length l= / 100 (the higher the subdivision density, the higher the interpolation accuracy), for each tiny subgrid, take its center point. : The "ray-mapping method" is used to determine whether each center point is within the surface feature. The number of tiny sub-mesh units *ti* within the surface is counted, and the coverage ratio of the edge mesh *Gi* is determined. ; (If the edge grid is only tangent to the boundary of the surface feature, then) =0, not included in the total area); (3) The sum of the coverage areas of all edge grids is ; Step 4: Total Area Calculation and Error Correction (1) The preliminary calculated area of ​​a surface feature is the sum of the area of ​​the complete grid and the area of ​​the edge grid. ; (2) For higher precision, the mesh subdivision level m can be increased (e.g., from level 8 4m mesh to level 12 25cm mesh), or the subdivision density n of the edge mesh can be increased (e.g., from 100×100 to 200×200). The corrected area is... ; (in For correction factors, based on actual calibration experiments, generally | |≤0.001, i.e., error ≤0.1%.

[0045] Example 1: To verify the feasibility of this invention in practice, it was applied to a complex regional spatial data processing scenario. In this scenario, there is a large amount of spatial entity data from different sources with varying accuracy. Furthermore, traditional latitude-longitude-based methods require repeated coordinate transformations and curvature corrections during distance and area calculations, resulting in complex and time-consuming calculations. Additionally, significant calculation errors exist in localized areas, making it difficult to meet the needs of refined spatial analysis. In this scenario, the spatial range data of the target region is first acquired. Latitude-longitude analytical operations are then performed on the region boundaries to extract boundary coordinate values, forming a baseline grid set. Subsequently, using the lower left corner coordinates of the baseline grid as a reference point, a regional Cartesian coordinate system is constructed, converting the latitude and longitude values ​​of each point within the region into planar coordinate values, ensuring all spatial calculations are performed under a unified coordinate system. Based on this, a recursive quadrilateral partitioning operation is performed on the spatial coordinate set using a preset side length as the 0th-level grid scale, recording row and column numbers at each level to form a multi-level equidistant grid structure. After the mesh structure is established, GeoSOT level 15 mesh is used as the reference unit. A local partitioning origin is established at the reference unit location. Mapping calculations are performed on the equidistant mesh level and the reference mesh level to ensure that the local computational mesh is consistent with the standard mesh system. During the encoding stage, a combination of hexadecimal and quaternary codes is written to the mesh unit, and the column offset and row offset values ​​are recorded at the encoding position to form a unique mesh encoding sequence. Subsequently, spatial entities are mapped to the corresponding mesh sets, and spatial objects are expressed through mesh sets. During the calculation process, the mesh row and column numbers are obtained by parsing the mesh codes, and the mesh center coordinates are further calculated. Euclidean distance is used to calculate spatial distance, and the area value is calculated through full mesh counting and edge mesh subdivision. Throughout the calculation process, the complex latitude and longitude surface calculations are no longer relied upon, thus significantly reducing computational complexity. In this scenario, the computational effects of traditional methods and the method of this invention are compared, and the following data results are obtained.

[0046] Table 1: Comparison of Spatial Grid Calculation Results

[0047] Further detailed analysis can be conducted using the data in Table 1. Regarding the average distance error, the error is 2.35m for traditional method A, 2.68m for traditional method B, 2.12m for traditional method C, and 0.48m for the method of this invention. Numerical comparison shows that the method of this invention reduces the error by approximately 79.6% compared to traditional method A, approximately 82.1% compared to traditional method B, and approximately 77.4% compared to traditional method C. This change demonstrates that, under the same data conditions, mapping spatial locations to equidistant grids and performing calculations based on the grid center coordinates effectively avoids the accumulation of errors caused by surface factors in traditional latitude and longitude calculations, resulting in more stable and consistent distance calculation results.

[0048] Regarding the area calculation error index, traditional method A has an error of 3.21%, traditional method B has an error of 3.75%, traditional method C has an error of 2.98%, while the method of this invention has an error of 0.52%. The comparison shows that the method of this invention reduces the error by approximately 83.8% compared to traditional method A, approximately 86.1% compared to traditional method B, and approximately 82.5% compared to traditional method C. This result indicates that by representing the spatial region in a grid and combining full grid counting with edge grid subdivision calculation, area results closer to the true value can be obtained under complex boundary conditions, reducing the error impact caused by projection distortion and boundary irregularities in traditional methods.

[0049] In terms of computation time, traditional method A takes 185ms, traditional method B takes 210ms, traditional method C takes 172ms, while the method of this invention takes only 46ms. It can be seen that the method of this invention reduces computation time by approximately 75.1% compared to traditional method A, approximately 78.1% compared to traditional method B, and approximately 73.3% compared to traditional method C. This change demonstrates that replacing the traditional coordinate calculation process with grid encoding parsing and row / column number calculation can reduce complex calculation steps and significantly reduce computation time, especially in large-scale data processing scenarios where the advantages are even more pronounced.

[0050] Regarding the stability index, traditional method A achieves 0.72, traditional method B 0.68, and traditional method C 0.75, while the method of this invention reaches 0.91. This represents an improvement of approximately 33.8% compared to the lowest value of 0.68 and approximately 26.4% compared to the average value of approximately 0.72. This result indicates that the method of this invention exhibits minimal fluctuations in results during continuous computation, maintaining high consistency and repeatability, making it suitable for spatial computing tasks requiring long-term operation and batch processing.

[0051] Based on the above four indicators, it can be concluded that the method of the present invention is superior to the traditional method in terms of distance accuracy, area accuracy, computational efficiency and stability. All indicators show a significant improvement trend, which verifies the effectiveness and reliability of the local equidistant grid and coding calculation method in spatial data processing.

[0052] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A spatial mesh generation method based on geometric constraints, characterized in that, Includes the following steps: Acquire spatial range data of the target area, perform latitude and longitude analysis calculations, extract boundary coordinate values, and generate a reference grid set; Using the coordinates of the lower left corner of the baseline grid as a reference point, a regional Cartesian coordinate system is constructed, and the latitude and longitude values ​​of each point in the target area are converted into Cartesian coordinate values ​​to form a set of spatial coordinates. In the Cartesian coordinate system, a recursive quadrilateral subdivision operation is performed on the spatial coordinate set using a preset side length as the 0th level grid scale. The row number and column number are recorded at each level of subdivision to generate a multi-level equidistant grid cell set. Using GeoSOT level 15 mesh as the reference unit, a local partitioning origin is established at the reference unit location. Mapping calculations are performed on the equidistant mesh level and the reference mesh level to form a fused mesh structure. Establish an encoding structure at the grid cell position, write hexadecimal encoding for the 0th level cell, write quaternary encoding for each level of subdivision, record column offset and row offset values ​​at the encoding position, and generate a grid encoding sequence. Establish geographic entity mapping relationships at the encoding location, perform grid coverage parsing on the geographic entities, and generate a grid set; Perform calculations at the location of the grid set, calculate the Euclidean distance to the grid center coordinates to generate distance values, perform surface accumulation calculations to generate area values, and perform intersection and union calculations to generate topological relationship results. The process of performing mapping calculations between the equidistant grid levels and the reference grid levels to form a fused grid structure includes the following steps: Extract the GeoSOT Level 15 grid cell identifier from the reference grid set location, read the latitude and longitude range values ​​of the corresponding grid cell at the identifier location, write the reference grid boundary coordinate values ​​at the range location, and establish the reference grid index sequence at the coordinate location. Read the boundary coordinate values ​​of each equidistant grid cell at the location of the regional spatial coordinate set, extract the coordinate values ​​of the lower left corner of the grid at the boundary location, perform coordinate inverse calculation at the coordinate location, convert the planar coordinate values ​​into latitude and longitude values, and generate the corresponding latitude and longitude coordinate values ​​at the conversion location. Perform grid affiliation determination operation at the latitude and longitude coordinate value location. At the determination location, perform interval matching calculation between the converted latitude and longitude coordinate values ​​and the reference grid boundary coordinate values. Write the reference grid identifier at the location that meets the interval inclusion condition. Form the reference grid affiliation sequence at the equidistant grid cell location. Extract the level values ​​corresponding to each level grid cell at the position of the equidistant grid level index sequence, read the grid side length value at the level position, read the spatial scale value corresponding to the reference grid at the position of the reference grid index sequence, perform the proportional relationship calculation at the value position, write the level mapping coefficient value at the calculation position, and establish the scale correspondence between the equidistant grid level and the reference grid level at the mapping position. At the mapping relationship location, perform an association write operation on the reference grid identifier and the equidistant grid cell identifier, establish a mapping index table at the write location, and record the reference grid identifier, equidistant grid level value, and equidistant grid cell identifier at the index table location to form a set of grid mapping relationships; Perform a consistency check operation at the location of the grid mapping relationship set, perform coverage detection on the equidistant grid cells corresponding to the same reference grid at the check location, determine whether each equidistant grid cell completely falls into the corresponding reference grid range at the detection location, and write a valid identifier at the location that meets the conditions. At the valid identification location, perform a structural reorganization operation on the mapping relationship set, perform aggregation calculation on the equidistant grid cells according to the reference grid identification at the reorganization location, and generate a fused grid structure data set at the aggregation location.

2. The spatial mesh generation method based on geometric constraints according to claim 1, characterized in that, The latitude and longitude resolution calculation specifically includes: Establish a boundary point index sequence for the spatial range data of the target area, extract longitude and latitude values ​​at the index positions, write the longitude and latitude component record table at the value positions, perform minimum and maximum value calculations at the longitude value positions, and perform minimum and maximum value calculations at the latitude value positions to generate longitude range values ​​and latitude range values. Perform integer alignment operations on the longitude range values ​​to generate longitude boundary values, perform integer alignment operations on the latitude range values ​​to generate latitude boundary values, construct an outer rectangular range between the longitude boundary values ​​and the latitude boundary values, and generate a set of region boundary coordinates. Perform a grid matching operation at the location of the region boundary coordinate set to map the region boundary coordinate set to the GeoSOT level 15 grid. Generate a reference grid index sequence at the mapping location, perform a deduplication operation at the index location, retain a unique grid identifier, and form a reference grid set.

3. The spatial mesh generation method based on geometric constraints according to claim 2, characterized in that, The construction of the regional Cartesian coordinate system includes the following steps: Extract the longitude and latitude values ​​corresponding to the lower left corner grid cell at the reference grid set location, establish the reference point coordinates at this location, perform offset calculations to generate the origin of the region partitioning at the reference point location in the west and south directions, and write the origin identifier at the origin point location. The coordinate axis directions are defined by establishing the origin of the region partitioning, establishing the X-axis in the horizontal direction and the Y-axis in the vertical direction, and writing the direction labels at the coordinate axis positions to form a Cartesian coordinate system structure; The latitude and longitude component record table reads the longitude and latitude values, performs geodetic coordinate transformation at the value location, maps the latitude and longitude values ​​to plane coordinate values, and generates the corresponding X coordinate value and Y coordinate value at the mapping location. The coordinate values ​​are used to perform a range determination operation on the X and Y coordinate values. If the position meets the condition of the first quadrant, a valid identifier is written and the set of spatial coordinates of the region is output.

4. The spatial mesh generation method based on geometric constraints according to claim 3, characterized in that, The quadrilateral partitioning operation specifically includes: Establish a grid index sequence at the location of the regional spatial coordinate set, extract the X and Y coordinate values ​​corresponding to each point at the index location, set the initial grid side length value at the value location as the 0th level grid scale, and construct the 0th level grid coverage area at the scale location. Perform a binary partitioning calculation at the coverage area of ​​the level 0 grid, dividing the current grid cell into the lower left sub-cell, lower right sub-cell, upper left sub-cell, and upper right sub-cell. Write the column index value and row index value at the sub-cell position respectively, and form a sub-cell identifier set at the index position. Establish a hierarchical index identifier at the grid index sequence position, record the current subdivision hierarchical value at the hierarchical position, perform coordinate range calculation on each sub-unit at the hierarchical position, write the sub-unit boundary coordinate value at the range position, perform an inclusion determination operation on the sub-unit boundary coordinate value and the regional spatial coordinate set at the index position, and retain the sub-unit identifier at the position that meets the inclusion condition. Recursive partitioning is performed at the reserved sub-cell identifier position. At the recursive position, partitioning calculation and range determination calculation are repeatedly performed on each sub-cell. At each level of partitioning position, column index values ​​and row index values ​​are written into the hierarchical record table in hierarchical order to form a multi-level grid hierarchical index sequence. Perform edge length update calculation on the grid cell at the hierarchical record table position, perform binary search calculation on the current edge length value at the update position, write the next level edge length value at the calculation position, and establish the correspondence between grid cell and edge length value at the mapping position. Execute a termination decision operation at the multi-level grid level index sequence position. At the decision position, compare the current level value with the preset maximum level value. Stop the recursive subdivision at the position that meets the condition. Output the final set of grid cell identifiers at the termination position. Perform equidistant verification calculations at the final grid cell identifier set location, perform consistency judgment on the side length values ​​of adjacent grid cells at the verification location, and output multi-level equidistant grid cell sets at the location that meets the consistency conditions.

5. The spatial mesh generation method based on geometric constraints according to claim 1, characterized in that, The generated grid-coded sequence specifically includes: Establish an coded index sequence at the grid cell index position, extract the 0th level grid cell identifier at the index position, convert the grid cell number value to a hexadecimal character value at the identifier position, write the coded prefix value at the character position, and form the initial coded field at the coded position. Extract the column index value and row index value corresponding to each level grid cell at the position of the equidistant grid level index sequence. Perform the sub-cell position determination operation at the index position according to the reverse Z order rule. Generate 4-ary code value at the determination position based on the column index value and row index value. Write the code field in the order of the level at the code position to form a level code sequence at the field position. At the encoding field position, perform concatenation calculation on the hexadecimal encoded prefix value and the hierarchical encoded sequence, generate a complete grid encoded sequence at the concatenation position, perform a uniqueness determination operation at the encoding position, and output the grid encoded sequence at the position that satisfies the uniqueness condition.

6. The spatial mesh generation method based on geometric constraints according to claim 5, characterized in that, The process of calculating the Euclidean distance to the grid center coordinates to generate a distance value includes the following steps: Extract the target grid code value at the grid code sequence position. Perform hexadecimal parsing operation on the code prefix at the code position to generate the column number and row number value of the level 0 grid. Perform bit-by-bit parsing operation on the 4-ary code sequence at the code position. Extract the column offset value and row offset value according to the reverse Z order rule at the parsing position. Perform recursive calculation on the column number value and row number value at the parsing position to generate the final column number value and final row number value corresponding to the target grid cell. Combine the grid side length value with the column number position and row number position to perform coordinate calculation to generate the coordinate value of the lower left corner of the grid. At the coordinate position, the center offset calculation is performed between the coordinate value of the lower left corner of the grid and the grid side length to generate the center coordinate value of the grid. At the center coordinate position, the center coordinate values ​​corresponding to the two grid cells are extracted. At the coordinate position, the difference calculation is performed on the corresponding X coordinate value and Y coordinate value respectively. At the difference position, the square calculation is performed on the difference value. At the square position, the sum calculation is performed. At the calculation position, the square root operation is performed on the sum result to generate the distance value.

7. The spatial mesh generation method based on geometric constraints according to claim 6, characterized in that, The specific steps of performing surface accumulation and calculation on the mesh set to generate area values ​​include: Extract the grid code values ​​at the grid set index position, perform hexadecimal prefix parsing operation at the code position to generate initial column number and row number values, perform bit-by-bit parsing operation on the 4-ary code sequence at the code position, extract column offset and row offset values ​​according to the reverse Z-order rule at the parsing position, perform recursive calculation on the column number and row number values ​​at the parsing position to generate final column number and final row number values, perform coordinate calculation on the column number and row number positions combined with the grid side length value to generate the coordinate value of the lower left corner of the grid, and construct the coordinate values ​​of the four vertices of the grid at the coordinate position; At the grid set index position, a coverage determination operation is performed on each grid cell. At the determination position, a spatial inclusion detection is performed on the coordinate values ​​of the four vertices of the grid. A complete grid identifier is written at the position where all vertices meet the target area range condition. An edge grid identifier is written at the position where some vertices meet the condition. At the complete grid identifier position, an area calculation is performed according to the grid side length value to generate a complete area value. At the edge grid identifier position, a subdivision calculation is performed on the grid cell. At the subdivision position, the grid cell is divided into several sub-cells. The center coordinate value is extracted at the sub-cell position. A spatial inclusion detection is performed at the coordinate position. At the detection position, the number of sub-cells that meet the condition is counted. At the count position, the coverage ratio value is calculated. At the ratio position, the area is calculated by combining the grid side length value to generate an edge area value. At the area position, the complete area value and the edge area value are accumulated to generate a total area value.

Images