Optimization of Numerical Methods for Transforming UTM Plane Coordinates to Lambert Plane Coordinates

: The rapid transformation from UTM (Universal Transverse Mecator) projection to Lambert projection helps to realize timely merging, inversion, and analysis of high-frequency partitioned remote sensing images. In this study, the transformation error and the efficiency of the linear rule approximation method, the improved linear rule approximation method, the hyperbolic transformation method, and the conformal transformation method were compared in transforming the coordinates of sample points on WGS84 (The World Geodetic System 1984)-UTM zonal projections to WGS84-Lambert projection coordinates. The effect of the grid aspect ratio on the coordinate transformation error of the conformal transformation method was examined. In addition, the conformal transformation method-based error spatial pattern of the sample points was analyzed. The results show that the conformal transformation method can better balance error and efficiency than other numerical methods. The error of the conformal transformation method is less affected by grid size. The maximum x-error is less than 0.36 m and the maximum y-error is less than 1.22 m when the grid size reaches 300 km × 300 km. The x- and y-error values decrease when square grids are used; namely, setting the grid aspect ratio close to 1 helps to weaken the effect of increasing grid area on the error. The dispersion of the error distribution and the maximum error of sample points both decrease relative to their minimum distance to the grid edge and stabilize at a minimum distance equal to 70 km. This study can support the rapid integration of massive remote sensing data over large areas.


Introduction
In recent years, the rapid growth in the volume of remote sensing data has raised high demands for the performance of spatial data coordinate transformations [1][2][3][4][5][6]. Multisource remote sensing images are generally organized on the WGS1984 (The World Geodetic System 1984)-UTM (Universal Transverse Mecator) projection coordinate system with nonuniform segmentation rules and diverse spatial and temporal resolutions. The rapid transformation of spatial data from the UTM zonal projection coordinate system (applied to the local area) to a conical or cylindrical projection coordinate system (e.g., Lambert projection coordinate system) is important for realizing timely merging, retrieval, and analysis of high-frequency segmented remote sensing images.
Coordinate transformations can be performed by either analytical or numerical methods. Analytical methods have a rigorous foundation of mathematics and require known reference ellipsoid parameters and analytical equations; they are generally more accurate algorithm was promoted for transformations between different conformal projections and for forward and inverse transformation of conformal projection coordinates, showing low parameter requirements, low error, and high stability [9]. In addition, the implementation of the conformal polynomial method is more flexible and easier than the finite element method and the difference quotient method [20,21].
In addition to the type of polynomial and its method of solving the parameters, the trade-off between coordinate transformation error and efficiency is comprehensively influenced by the local transformation area and the number and distribution of control points [13]. The transformation error generally increases as the area of the transformation region increases. In practice, the study area is split into several subregions by regular rectangular or triangular grids to mitigate error propagation. For each subregion, independent numerical polynomial parameters are fitted by using the grid vertices as control points. The transformation error can be reduced to some extent by decreasing the area of the subregion. When the number of subregions increases to a certain level, the efficiency of the coordinate transformation decreases as the time spent on parameter retrieval increases, and the error stabilizes.
Research on numerical methods and implementation pathways has provided important theoretical and empirical support for this study. However, the applicability of different numerical methods varies with the coordinate transformation scenario, the area, and the shape of the subregion. Which numerical method is more suitable for transforming raster data from UTM zonal projection to Lambert projection? What is the appropriate grid layout scheme? Fewer studies have addressed these questions.
In this study, the linear rule approximation method, the improved linear rule approximation method, the hyperbolic transformation method, and the conformal transformation method were applied to transform the coordinates of the points in the three WGS84-UTM projection zones (i.e., UTM Zone 48N; UTM Zone 50N; UTM Zone 52N) into WGS84-Lambert projection coordinates. For each study area, we mimicked the raster data organization by deploying 33 million uniformly distributed sample points at 200 m intervals. Fifteen square grid segmentation schemes (i.e., 1 × 1; 5 × 5; 10 × 10; 25 × 25; 50 × 50; 75 × 75; 100 × 100; 125 × 125; 150 × 150; 175 × 175; 200 × 200; 225 × 225; 275 × 275; and 300 × 300; unit: km 2 ) were applied to the study area, and then the coordinate transformation errors and efficiency of each numerical method were tested when applying each scheme. The effect of the grid aspect ratio on the error of the conformal transformation method was examined by calculating the coordinate transformation error by setting the grid length in the x-direction to 100, 200, 300, and 400 km and y-direction to 300, 350, 400, 450, 500, 550, and 600 km. In addition, the error distribution pattern of sample points based on the conformal transformation method was analyzed. This study can support the rapid integration of massive remote sensing data in large areas.
The rest of the paper is organized as follows: Section 2 describes in detail the four numerical methods and experimental scenarios involved in this study. Section 3 presents the experimental results of transforming WGS1984 UTM plane coordinates to WGS1984 Lambert plane coordinates using the four numerical methods. Finally, the comparison and application scenarios of the four methods are discussed and summarized in Sections 4 and 5.

Numerical Methods
As shown in Figure 1, the four corner points of the rectangular grid are set as A, B, C, and D, and point G is an arbitrary point within the grid. Their coordinates in the WGS84-based UTM projection coordinate system are ( Their corresponding coordinates in the Lambert projection coordinate system are ( On this basis, the application of the linear rule approximation method, the improved linear rule approximation method, the hyperbolic transformation method, and the conformal trans- approximation method, the improved linear rule approximation method, the hyperbolic transformation method, and the conformal transformation method are illustrated by calculating the coordinates of point G (i.e., ( , )).
The ILRA method adds a correction for deformation in the x-and y-directions based on the LRA method. Its formula is shown in Equations (7) and (8).

Hyperbolic Transformation (HT) Method
The equation for the HT method is shown in Equations (9) and (10). show the distribution of control points based on the UTM projection coordinate system and Lambert projection coordinate system. The corner points A, B, C, and D of the grid and corner points B and D of the adjacent grids are the control points, and G is the sample point whose coordinate is transformed within the grid.

Linear Rule Approximation (LRA) Method
The linear rule approximation method equations are listed in Equations (1) and (2).
where the parameters m, n, p, q are calculated by Equations (3)-(6), ∆x = x B − x A , and The ILRA method adds a correction for deformation in the x-and y-directions based on the LRA method. Its formula is shown in Equations (7) and (8).

Hyperbolic Transformation (HT) Method
The equation for the HT method is shown in Equations (9) and (10).
The eight parameters a 1 -a 4 and b 1 -b 4 can be solved by bringing points A, B, C, and D as control points into the two systems of linear equations, as shown in Equations (11) and (12).
where Q i and P i are calculated as shown in Equations (15)- (18).
The six parameters a 1 -a 3 and b 1 -b 3 can be solved directly by bringing the values of the coordinates of points A, B, C, and D into the system of linear equations in Equation (9) as follows, where the subscript j denotes the number of the three control points and takes values 1, 2, and 3. The three control points are selected as points B, C, and D, in turn. P i1 and Q i1 , P i2 and Q i2 , and P i3 and Q i3 can be calculated in Equations (19) and (20) by replacing x G and y G with x B and y B , x C and y C , and x D and y D , respectively.

Test Scenarios
In this study, the reference datum for both the UTM projection coordinate system and the Lambert projection coordinate system is the WGS-84 ellipsoid. The projection parameters (Table 1) applicable to China are applied to the Lambert projection coordinate system, whose central meridian longitude is the same as that of UTM Zone 48N. Since the Lambert projection deformation is symmetrical on both sides of the central meridian and is influenced by the distance from the central meridian, the experiments were performed in three projection zones: UTM Zone 48N (90 • N-96 • N), UTM Zone 50N (102 • N-108 • N), and UTM Zone 52N (114 • N-120 • N). For each projection zone, the spatial extent of the study area was designed to be 300-700 km (X-coordinate) and 2450-5750 km (Y-coordinate), as shown in Figure 2. For each study area, 33 million uniformly distributed sample points were deployed at 200 m intervals, mimicking the raster data organization. Fifteen groups of square grids based on row-column ordering (i.e., 1 × 1, 5 × 5, 10 × 10, 25 × 25, 50 × 50, 75 × 75, 100 × 100, 125 × 125, 150 × 150, 175 × 175, 200 × 200, 225 × 225, 275 × 275, and 300 × 300, unit: km 2 ) were applied to divide the entire study area into numerous subareas. On this basis, the linear rule approximation method, improved linear rule approximation method, hyperbolic transformation method, and conformal transformation method were applied to transform the coordinates of the sample points from the UTM projection coordinate system into the Lambert coordinate system. In addition, the coordinate transformation of sample points was implemented as follows: 1.
In the first step, for a particular numerical method, the four vertices of each subzone are used as control points to calculate their corresponding polynomial parameters. The polynomial parameters are stored in a memory array.

2.
In the second step, all sample points are grouped according to the spatial extent of the subzone. Sample points within the same subzone are assigned to the same group.

3.
In the third step, for each group, the polynomial parameters are called from the memory array, and the coordinate transformation is performed in bulk for the sample points belonging to that group. Another experiment was carried out to test the effect of the grid aspect ratio on the error of the conformal transformation method. A total of 4 × 7 rectangular grid divisions of the study area were designed with lengths of 100, 200, 300, and 400 km in the x-direction and 300, 350, 400, 450, 500, 550, and 600 km in the y-direction. Remote Sens. 2022, 14, x FOR PEER REVIEW 7 of 18

Hardware and Software Environment
The hardware and software environment for this study is shown in Table 2. The results of the projection coordinate transformation based on the map projection library PROJ are considered 'true values' and used to calculate the polynomial parameters and coordinate transformation errors for each numerical method.

Hardware and Software Environment
The hardware and software environment for this study is shown in Table 2. The results of the projection coordinate transformation based on the map projection library PROJ are considered 'true values' and used to calculate the polynomial parameters and coordinate transformation errors for each numerical method.

Error Statistics of Multiple Numerical Methods
The maximum, minimum, mean, and median errors of the linear rule approximation (LRA) method, improved linear rule approximation (ILRA) method, hyperbolic transformation (HT) method, and conformal transformation (CT) method for multiple grid sizes and multiple conversion areas were calculated and expressed as box plots, as shown in Figure 3. First, the x-error and y-error of the CT method are significantly lower than those of the other three methods for the same grid size. However, this does not mean that the other methods are not applicable because the errors can be controlled within the threshold value by controlling the grid size when they are used. Second, the x-error and y-error increase as the grid size increases, and the magnitude of the error increase gradually increases. The CT method has the smallest error increase, the maximum x-error is less than 0.36 m and the maximum y-error is less than 1.22 m when the grid size reaches 300 km × 300 km. The variation characteristics of the x-error of the LRA method and the ILRA method are almost the same, but the y-error of the LRA method is significantly higher than that of the ILRA method. The error is significantly smaller and more stable than that of the ILRA method. For the LRA and ILRA methods, the x-error is larger than the y-error; for the HT and CT methods, the x-error is smaller than the y-error. Third, as the distance between the conversion area and the central meridian of the Lambert projection increases, the maximum value and the dispersion of x-and y-errors of the four methods increase or decrease to different degrees. Among them, the effect of the conversion area change on the HT and CT methods is mainly focused on the x-error, while for the LRA and ILRA methods, the y-error is more affected.

Error Statistics of Multiple Numerical Methods
The maximum, minimum, mean, and median errors of the linear rule approximation (LRA) method, improved linear rule approximation (ILRA) method, hyperbolic transformation (HT) method, and conformal transformation (CT) method for multiple grid sizes and multiple conversion areas were calculated and expressed as box plots, as shown in Figure 3. First, the x-error and y-error of the CT method are significantly lower than those of the other three methods for the same grid size. However, this does not mean that the other methods are not applicable because the errors can be controlled within the threshold value by controlling the grid size when they are used. Second, the x-error and y-error increase as the grid size increases, and the magnitude of the error increase gradually increases. The CT method has the smallest error increase, the maximum x-error is less than 0.36 m and the maximum y-error is less than 1.22 m when the grid size reaches 300 km × 300 km. The variation characteristics of the x-error of the LRA method and the ILRA method are almost the same, but the y-error of the LRA method is significantly higher than that of the ILRA method. The error is significantly smaller and more stable than that of the ILRA method. For the LRA and ILRA methods, the x-error is larger than the y-error; for the HT and CT methods, the x-error is smaller than the y-error. Third, as the distance between the conversion area and the central meridian of the Lambert projection increases, the maximum value and the dispersion of x-and y-errors of the four methods increase or decrease to different degrees. Among them, the effect of the conversion area change on the HT and CT methods is mainly focused on the x-error, while for the LRA and ILRA methods, the y-error is more affected.  Where the x error = |x0 − x|, y error = |y0 − y|. x0, y0 are calculated by PROJ, and x, y are calculated by the above numerical method. To constrain the maximum error of coordinate transformation to be less than 5 m, we set the grid size of the LRA method, ILRA method, HT method, and CT method to 10 km × 10 km, 10 km × 10 km, 25 km × 25 km, and 275 km × 275 km, respectively. The cumulative error frequency of 33,000,000 sample points in CA2 was calculated, as shown in Figure 4. In terms of x-errors, the CT method has the fastest convergence rate, and the x-error of all sample points is less than or equal to 0.5 m. The percentages of sample points with x-errors less than 0.5 m for the HT method, LRA method, and ILRA method are 98.93%, 63.69%, and 63.42%, respectively. In terms of y-errors, the CT method and the LRA method converged to 1.25 m almost simultaneously. A total of 96.59% and 87.32% of the sample points are less than 1.25 m for the ILRA method and the HT method, respectively. The CT method has the highest comprehensive stability, even though its grid area is 121 times larger than that of the HT method and 756.25 times larger than that of the LRA and ILRA methods. To constrain the maximum error of coordinate transformation to be less than 5 m, we set the grid size of the LRA method, ILRA method, HT method, and CT method to 10 km × 10 km, 10 km × 10 km, 25 km × 25 km, and 275 km × 275 km, respectively. The cumulative error frequency of 33,000,000 sample points in CA2 was calculated, as shown in Figure 4. In terms of x-errors, the CT method has the fastest convergence rate, and the x-error of all sample points is less than or equal to 0.5 m. The percentages of sample points with xerrors less than 0.5 m for the HT method, LRA method, and ILRA method are 98.93%, 63.69%, and 63.42%, respectively. In terms of y-errors, the CT method and the LRA method converged to 1.25 m almost simultaneously. A total of 96.59% and 87.32% of the sample points are less than 1.25 m for the ILRA method and the HT method, respectively. The CT method has the highest comprehensive stability, even though its grid area is 121 times larger than that of the HT method and 756.25 times larger than that of the LRA and ILRA methods.

Comparison of the Transformation Efficiency of Multiple Numerical Methods
The time elapsed of coordinate transformation for the four numerical methods at grid sizes of 10 km × 10 km, 10 km × 10 km, 25 km × 25 km, and 275 km × 275 km was calculated, as shown in Figure 5. The number of sample points is 1.32 × 10 6 , 5.28 × 10 7 , 1.056 × 10 8 , 1.65 × 10 8 , 2.112 × 10 8 , 2.64 × 10 8 , and 3.3 × 10 8 . The coordinate transformation was implemented by a single thread, and the time elapsed was taken as the average of the calculated results for different conversion areas (i.e., CA1; CA2; CA3). The computational time elapsed of the four methods increases linearly with the increase of the number of sample points and is steadily expressed as HT < LRA < ILRA < CT. The time elapsed of the CT method is longer than that of the other methods but with a small difference. The number of grid points of a Sentinel 2 multispectral remote sensing image (spatial resolution: 10 m) is nearly 100 million, and the difference in computation time elapsed is less than 1.5 s when the four numerical methods are executed on a single thread. sizes of 10 km× 10 km, 10 km × 10 km, 25 km × 25 km, and 275 km × 275 km was calculated, as shown in Figure 5. The number of sample points is 1.32 × 10 6 , 5.28 × 10 7 , 1.056 × 10 8 , 1.65 × 10 8 , 2.112 × 10 8 , 2.64 × 10 8 , and 3.3 × 10 8 . The coordinate transformation was implemented by a single thread, and the time elapsed was taken as the average of the calculated results for different conversion areas (i.e., CA1; CA2; CA3). The computational time elapsed of the four methods increases linearly with the increase of the number of sample points and is steadily expressed as HT < LRA < ILRA < CT. The time elapsed of the CT method is longer than that of the other methods but with a small difference. The number of grid points of a Sentinel 2 multispectral remote sensing image (spatial resolution: 10 m) is nearly 100 million, and the difference in computation time elapsed is less than 1.5 s when the four numerical methods are executed on a single thread.

Effect of Grid Shape Change on the Maximum Coordinate Transformation Error of the CT Method
We used CA1 to calculate the x maximum error and y maximum error of the CT method under different grid shapes and areas. The grid shape and area were adjusted by setting its length in the x-direction to 100 km, 200 km, 300 km, and 400 km, and its length in the y-direction to 300 km, 350 km, 400 km, 450 km, 500 km, and 600 km, as shown in Figure 6a,b.
First, the maximum x-error and maximum y-error both increase with increasing length in the x-direction of the grid, and the increase becomes more drastic when the length in the y-direction of the grid is greater than 450 km. Second, the x maximum error and y maximum error increase with increasing length in the x-direction of the grid when the length in the y-direction of the grid is the same. This phenomenon is mainly caused by the change in the grid area. Third, for the CT method, the grid area is not the only dominant factor affecting the error; another is the grid shape. Setting the grid aspect ratio close to 1 helps to weaken the effect of increasing grid area on the error. For example, in Figure 6a,b, the maximum x-error and maximum y-error of the 300 × 300 (unit: km 2 ) grid

Effect of Grid Shape Change on the Maximum Coordinate Transformation Error of the CT Method
We used CA1 to calculate the x maximum error and y maximum error of the CT method under different grid shapes and areas. The grid shape and area were adjusted by setting its length in the x-direction to 100 km, 200 km, 300 km, and 400 km, and its length in the y-direction to 300 km, 350 km, 400 km, 450 km, 500 km, and 600 km, as shown in Figure 6a are smaller than those of the 200 × 400 (unit: km 2 ) grid and the 100 × w (w = 450, 500, 550, 600) (unit: km 2 ) grid, although the former has a larger grid area than the latter. Similarly, the maximum x-error and maximum y-error of the 400 × 400 (unit: km 2 ) grid are almost the same as those of the 100 × 600 (unit: km 2 ) grid, although the area of the former grid is nearly 2.7 times larger than that of the latter. Fourth, setting the grid size to 300 × 300 (unit: km 2 ) is more suitable for applying the CT method to perform the transformation of the sample point set from the UTM planar coordinate system to the Lambert planar coordinate system because it considers both the grid size and the errors.
In another experiment, the coordinate transformation error of the CT method was calculated by setting the ratio of the length in the y-direction to the length in the x-direction of the grid (RYX) to 16, 6.25, 4, 1.5625, 1, 0.64, 0.25, 0.16, and 0.0625, while keeping the grid area equal to 10,000 km 2 , as shown in Figure 6c,d. The maximum x-error and maximum y-error are minimized when RYX is 1. The error increases as the ratio deviates from 1. In addition, keeping the grid area constant and stretching the grid along the x-direction and the y-direction in equal proportions, the variation of errors is not consistent. The grid with a ratio equal to 1/16 can be regarded as the grid with a ratio equal to 16 placed horizontally, and their stretching ratios are equal, but the maximum x-error of the former is smaller than that of the latter, and the maximum y-error of the former is larger than that of the latter. Figure 6. (a,b) show the variation of max error of CT method when fixing the length in the x-direction (100, 200, 300, and 400 km) and the length in the y-direction of the grid; (c,d) show the variation of max error of CT method for different RYX with constant grid area (10,000 km 2 ). For a specific grid system, RYX represents the ratio of the grid length in the y-direction to that in the x-direction.

Spatial Distribution of Coordinate Transformation Error of the CT Method
The spatial distribution of the x-error and y-error of the sample point set was expressed by the raster data model. For each conversion area (i.e., CA1, CA2, and CA3), 33,000,000 sample points were arranged in a matrix of 2000 × 16,500 with a 200 m interval. Figure 6. (a,b) show the variation of max error of CT method when fixing the length in the x-direction (100, 200, 300, and 400 km) and the length in the y-direction of the grid; (c,d) show the variation of max error of CT method for different RYX with constant grid area (10,000 km 2 ). For a specific grid system, RYX represents the ratio of the grid length in the y-direction to that in the x-direction.
First, the maximum x-error and maximum y-error both increase with increasing length in the x-direction of the grid, and the increase becomes more drastic when the length in the y-direction of the grid is greater than 450 km. Second, the x maximum error and y maximum error increase with increasing length in the x-direction of the grid when the length in the y-direction of the grid is the same. This phenomenon is mainly caused by the change in the grid area. Third, for the CT method, the grid area is not the only dominant factor affecting the error; another is the grid shape. Setting the grid aspect ratio close to 1 helps to weaken the effect of increasing grid area on the error. For example, in Figure 6a,b, the maximum x-error and maximum y-error of the 300 × 300 (unit: km 2 ) grid are smaller than those of the 200 × 400 (unit: km 2 ) grid and the 100 × w (w = 450, 500, 550, 600) (unit: km 2 ) grid, although the former has a larger grid area than the latter. Similarly, the maximum x-error and maximum y-error of the 400 × 400 (unit: km 2 ) grid are almost the same as those of the 100 × 600 (unit: km 2 ) grid, although the area of the former grid is nearly 2.7 times larger than that of the latter. Fourth, setting the grid size to 300 × 300 (unit: km 2 ) is more suitable for applying the CT method to perform the transformation of the sample point set from the UTM planar coordinate system to the Lambert planar coordinate system because it considers both the grid size and the errors.
In another experiment, the coordinate transformation error of the CT method was calculated by setting the ratio of the length in the y-direction to the length in the x-direction of the grid (RYX) to 16, 6.25, 4, 1.5625, 1, 0.64, 0.25, 0.16, and 0.0625, while keeping the grid area equal to 10,000 km 2 , as shown in Figure 6c,d. The maximum x-error and maximum y-error are minimized when RYX is 1. The error increases as the ratio deviates from 1. In addition, keeping the grid area constant and stretching the grid along the x-direction and the y-direction in equal proportions, the variation of errors is not consistent. The grid with a ratio equal to 1/16 can be regarded as the grid with a ratio equal to 16 placed horizontally, and their stretching ratios are equal, but the maximum x-error of the former is smaller than that of the latter, and the maximum y-error of the former is larger than that of the latter.

Spatial Distribution of Coordinate Transformation Error of the CT Method
The spatial distribution of the x-error and y-error of the sample point set was expressed by the raster data model. For each conversion area (i.e., CA1, CA2, and CA3), 33,000,000 sample points were arranged in a matrix of 2000 × 16,500 with a 200 m interval. The applied numerical method was the CT method, and the grid size was 300 × 300 (unit: km 2 ). Considering that the x-error and y-error of the sample point follow the heavytailed distribution feature, the head/tail breaks model [22] was applied to the classification rendering of the sample point error.
As shown in Figure 7, first, the x-error and y-error of the sample points all exhibit an overall increasing trend with increasing y-coordinates. For the three transformation regions, the sample points with the largest errors are all distributed in the region with y-coordinates > 5450 km. Second, the x-errors of the sample point exhibit an overall increasing trend as the distance from the sample point to the central meridian of the Lambert projection increases. However, the effect of the distance from the sample point to the central meridian of the Lambert projection on the y-error of the sample point is reversed and concentrated at high latitudes. In the region with y-coordinates > 5450 km, the y-error of some sample points decreases with their deviation from the central meridian. Third, the coordinate transformation error of sample points is also influenced by their positions inside the grid. For each grid, the sample points with the largest and smallest x-error levels are both distributed in a "bow" shape at the edges of the grid. The sample points with the smallest y-error level are concentrated near the four vertices of the grid, and the sample points with the largest y-error level are concentrated in a "bow" shape in the middle of the four edges of the grid. In addition, the distribution of errors within the grid shows symmetric characteristics, with the x-error of the sample point being centrosymmetric and the y-error being axisymmetric. side the grid. For each grid, the sample points with the largest and smallest x-error levels are both distributed in a "bow" shape at the edges of the grid. The sample points with the smallest y-error level are concentrated near the four vertices of the grid, and the sample points with the largest y-error level are concentrated in a "bow" shape in the middle of the four edges of the grid. In addition, the distribution of errors within the grid shows symmetric characteristics, with the x-error of the sample point being centrosymmetric and the y-error being axisymmetric. Figure 7. (a,b) show the spatial distribution of the x-error and the y-error of CT method in CA1, CA2, and CA3, with a grid size of 300 km × 300 km. For each conversion area, sample points were organized as 400 km × 3300 km raster data with 200 m as the spatial resolution.
To further investigate the error distribution within the grid, the two groups of 90,000 sample points were arranged in a 300 × 300 matrix within the specific grid of 300 km × 300 km in CA1. The sample spacing is 1 km, and the coordinates of the vertices of the first grid are (300, 5450), (600, 5450), (600, 5750), and (300, 5750) (unit: km); the coordinates of the vertices of the second grid are (300, 2450), (600, 2450), (600, 2750), and (300, 2750) (unit: km). The x-coordinate intervals of these two grids are the same, but they are 3300 km apart in the y-axis direction. For each sample point in the first grid, the relationship between the coordinate transformation error and its minimum distance to the grid edge (MinDGE, which can be interpreted as the minimum value of the distance from the sample point to the four edges of the grid) is expressed in Figure 8a,b. The dispersion of the error distribution and the maximum error both decrease with increasing MinDGE and stabilize at MinDGE equal to 70 km. The maximum x-error is 25% of the global x-error, and the maximum y-error is 80% of the global maximum y-error for the sample point with MinDGE equal to 70 km. As MinDGE increases, the minimum x-and y-errors of sample points increase, the mean x-error decreases, and the mean y-error increases. As shown in Figure 8c,d, the coordinate transformation error of the sample point near the grid center converges stably and does not change with the y-coordinate of the conversion area. In the two grids, the transformation errors of the sample point stabilize at MinDGE equal to 70 km. coordinate transformation error and its minimum distance to the grid edge (MinDGE, which can be interpreted as the minimum value of the distance from the sample point to the four edges of the grid) is expressed in Figure 8a,b. The dispersion of the error distribution and the maximum error both decrease with increasing MinDGE and stabilize at MinDGE equal to 70 km. The maximum x-error is 25% of the global x-error, and the maximum y-error is 80% of the global maximum y-error for the sample point with MinDGE equal to 70 km. As MinDGE increases, the minimum x-and y-errors of sample points increase, the mean x-error decreases, and the mean y-error increases. As shown in Figure  8c,d, the coordinate transformation error of the sample point near the grid center converges stably and does not change with the y-coordinate of the conversion area. In the two grids, the transformation errors of the sample point stabilize at MinDGE equal to 70 km.

Comparison of Multiple Numerical Methods
In this study, the error and efficiency of four numerical methods (i.e., the linear rule approximation method, improved linear rule approximation method, hyperbolic transformation method, and conformal transformation method) were compared in transforming sample point coordinates from the UTM coordinate system to the Lambert coordinate system. The CT method and HT method have their own advantages in applications. The CT method is more suitable for applications with discrete sample point distributions and higher requirements for coordinate transformation accuracy. The HT method can meet the requirements for applications with a relatively concentrated sample point distribution and high requirements for coordinate transformation efficiency. The LRA method and ILRA method have shown poor suitability in our experiments because they are worse than the HT method in terms of coordinate transformation speed and worse than the CT method in terms of transformation accuracy and error convergence speed. The ILRA method is designed to weaken the transformation error of LRA by increasing the number of items in the Taylor series. In this study, however, the error and time elapsed of the ILRA method are higher than those of the LRA method.
The coordinate transformation error of the CT method is influenced by both the xand y-coordinates of the sample points. The effect of the y-coordinate can be weakened by reducing the grid size as the y-coordinate increases. Another important factor affecting the coordinate transformation error of the CT method is the distance from the sample points to the center of the grid because the sample points with the largest errors are close to the edges of the grid. This error can be reduced by constructing an overlapping grid system. In this strategy, sample points close to the grid edges (i.e., overlapping areas) are inside multiple grids and assigned to suitable grids on their distances to different grid edges. Table 3 lists the application comparison of the CT method, HT method, and LRA method. All applications met their needs. Yang et al. [9] applied the CT method to the projection transformation of topographic maps, and the transformation in a larger range still maintained low error. Its application value is underestimated, and it can be applied to remote sensing data processing. The HT method and LRA method have been applied to map projections in map digitization and map projections or inverse projections of large amounts of data [12,13].

Analysis of Potential Application Scenario
The Gauss-Krueger 3 degree/6 degree zonal projection coordinate system is widely used in the organization of regional thematic remote sensing data products in China. The merger and analysis of large-scale thematic raster data require transforming sample point coordinates from the Gauss-Krueger projection coordinate system into the Lambert projection coordinate system. This application requirement can benefit from our results because the Gauss-Krueger projection differs from the UTM projection only in the scale factor.
The Raster Dataset Clean and Reconstitution Multigrid (RDCRMG) was designed for efficient remote sensing data processing and management and has provided support for multiple large data volume studies (i.e., monitoring of vegetation dryness; crop type identification; and arable land quality analysis) [2,[23][24][25][26][27][28][29][30]. According to the architecture of RDCRMG, different types of raster images were subdivided into several independent blocks and distributed for storage in different data nodes by using the multigrid as a consistent partition unit. The World Geodetic System 1984 (WGS 84)-based Universal Transverse Mercator (UTM) 6-degree strip division projection coordinate system has been designated as an RDCRMG spatial reference. Partitioned calculation result data should be transformed into the Lambert projection coordinate system for seamless data integration. On that basis, the CT method was most suitable for integration into RDCRMG. In addition, the CT method can be used to improve the coordinate transformation efficiency between the projection coordinate system and global discrete grid system.

Conclusions
In this study, four numerical methods (i.e., the LRA method, ILRA method, HT method, and CT method) have been used to convert the raster simulated data from the UTM projection coordinate system into the Lambert projection coordinate system. The transformation error and efficiency of these methods were compared. The relationship between multiple grid layout schemes and the error was explored. The x-error and y-error of the CT method are significantly lower than those of the other three methods at the same grid size. The maximum x-error is less than 0.36 m and the maximum y-error is less than 1.22 m when the grid size reaches 300 km × 300 km. Constraining the maximum error of coordinate transformation to less than 5 m, the CT method has the highest comprehensive stability, even though its grid size is the largest. The computational time elapsed of the four methods increases linearly as the number of sample points increases and can be expressed as HT < LRA < ILRA < CT. The time consumption for coordinate transformations of approximately 100 million points by the four methods are all less than 10 s, showing high efficiency. The difference in time consumption among these methods is less than 1.5 s. Considering efficiency and error, the CT method is more suitable than other methods for the rapid transformation of spatial data from UTM projection into Lambert projection. For the CT method, in addition to the grid area, the grid shape also significantly affects the error. Setting the grid aspect ratio close to 1 helps to weaken the effect of increasing grid area on the error. The x-error and y-error of the sample points both exhibit an overall increasing trend as their y-coordinates increase. The x-error of the sample points shows an overall increasing trend as the distance from the sample point to the central meridian of the Lambert projection increases, while the y-error is the opposite. The coordinate transformation error of the sample points near the grid center converges stably. The dispersion of the error distribution and the maximum error of the sample points both decrease as their minimum distance to the grid edge decreases. Within the 300 km × 300 km grid at different latitude zones, the transformation errors of the sample points stabilize when the minimum distance equals to 70 km. This research can provide a reference for realizing the timely merging, inversion, and analysis of high-frequency remote sensing images.