Determining the area corrections affecting the map areas in GIS applications

An auxiliary sphere in which the azimuth values are preserved is determined on the ellipsoidal reference surface (Figure 1). The ellipsoidal area is calculated using the Gauss–Bonnet theorem.

Figure 1

The F area in Figure 1 is calculated with the MATLAB code from Equations (2) and (3). Also, a, b, e and e′ are ellipsoid parameters. The solution of I₁ integral in Equation (5) can be calculated using math computer software or the serial coefficients produced by Karney method. A₀, A₁ and A₂ azimuths and the spherical geodesic line length are calculated by using spherical trigonometry in Equations (3–14): (2) F=R2(A2−A1)+(F(σ2)−F(σ1)) F = {R^2}\left({{A_2} - {A_1}} \right) + \left({{F_{\left({{\sigma _2}} \right)}} - {F_{\left({{\sigma _1}} \right)}}} \right) (3) R2=a22+b22tan h−1ee {R^2} = {{{a^2}} \over 2} + {{{b^2}} \over 2}{{\tan \,{h^{- 1}}e} \over e} (4) F(σi)=e2a2cos A0sin A0I1(σi) i=1,2 {F_{\left({{\sigma _i}} \right)}} = {e^2}{a^2}\cos \,{A_0}\sin \,{A_0}{I_1}\left({{\sigma _i}} \right)\,i = 1,2 (5) I1(σi)=−∫π/2σit(e′2)−t(k2 sin2 σi)e′2−k2sin2 σisin σi2dσ {I_1}\left({{\sigma _i}} \right) = - \int_{\pi /2}^{{\sigma _i}} {{{t\left({e{'^2}} \right) - t\left({{k^2}\,{{\sin}^2}\,{\sigma _i}} \right)} \over {e{'^2} - {k^2}{{\sin}^2}\,{\sigma _i}}}{{\sin \,{\sigma _i}} \over 2}d\sigma} (6) t(x)=x+x−1+1sin h−1x t\left(x \right) = x + \sqrt {{x^{- 1}} + 1} \sin \,{h^{- 1}}\sqrt x (7) k2=e′2cos2 A0 {k^2} = e{'^2}{\cos ^2}\,{A_0} (8) A1=ph(cos β1sin β2−sin β1 cos β2 cos Δλ+i cos β2 sin Δλ) {A_1} = ph\left({\cos \,{\beta _1}\sin \,{\beta _2} - \sin \,{\beta _1}\,\cos \,{\beta _2}\,\cos \,\Delta \lambda + i\,\cos \,{\beta _2}\,\sin \,\Delta \lambda} \right) (9) A2=ph(−sin β1cos β2+cos β1 sin β2 cos Δλ+i cos β1 sin Δλ) {A_2} = ph\left({- \sin \,{\beta _1}\cos \,{\beta _2} + \cos \,{\beta _1}\,\sin \,{\beta _2}\,\cos \,\Delta \lambda + i\,\cos \,{\beta _1}\,\sin \,\Delta \lambda} \right) (10) Δλ=(L2−L1)/1−e2(cos β1+cos β22)2 \Delta \lambda = \left({{L_2} - {L_1}} \right)/\sqrt {1 - {e^2}{{\left({{{\cos \,{\beta _1} + \cos \,{\beta _2}} \over 2}} \right)}^2}} (11) A0=ph(cos A1+i sin A1 sin β1+i sin A1 cos β1) {A_0} = ph\left({\cos \,{A_1} + i\,\sin \,{A_1}\,\sin \,{\beta _1} + i\,\sin \,{A_1}\,\cos \,{\beta _1}} \right) (12) σj=ph(cos Aj cos βj+i sin βj) j−1,2 {\sigma _j} = ph\left({\cos \,{A_j}\,\cos \,{\beta _j} + i\,\sin \,{\beta _j}} \right)\,j - 1,2 (13) |x+iy|=x2+y2 \left| {x + iy} \right| = \sqrt {{x^2} + {y^2}} (14) ph (x+iy)=arctan(yx) ph\,\left({x + iy} \right) = \arctan \left({{y \over x}} \right)

Maps can be produced using equal area projections to minimise area deformation, instead of Karney solution for area calculation on the ellipsoid reference surface. Large and medium-scaled topographic maps (1000–250,000) are produced using UTM or Lambert Conformal Conic (LCC) projection in the world. Area corrections (F – f) can be minimised with the selection of equal area projections. For calculating the area deformation, area correction formulas (F – f) calculated with the help of coordinates produced from the maps are presented to the users by using UTM and LCC systems. Area correction formulas (F–f) and area calculation on the ellipsoidal reference surface are not available in CAD and GIS applications. Therefore, instead of conformal projections, one of the projections that equal area projections should be selected for high accuracy calculation of a parcel area.

In this study, first, 15 forest parcels that varied between 300 and 3000 ha were selected within the boundaries of Kütahya Forest Regional Directorate (Figure 2). The area values of the 15 forest parcels selected in the study area were calculated with the equal area projections presented in Table 2 using GIS software. Thus, the most suitable equal area projection was determined by calculating the area deformations of the all projections (ITRF96 datum, 3° of longitude in width).

Figure 2

Areas of forest parcels were calculated by using the projections given in Table 2. At the same time, the areas of these forest parcels were calculated by the Karney method. The areas of the forest parcel calculated by the projections (Table 2) were compared with the areas calculated by the Karney method (Table 3). The distortion values are randomly distributed depending on the shape and size of the parcel and the distance of the parcel from the projection centre.

In Table 3, area corrections are calculated for the all equal area projections. According to Table 3 and Figure 3, two types of projections (Albers Equal Area Conic [AEAC] and Lambert Azimuthal Equal Area [LAEA]) with the least area corrections were determined. One of these projections can be preferred. AEAC projection was selected in this study in terms of the ease of calculation and processing time.

Figure 3

The area distortions have been also examined by employing Karney method in UTM. Thus, it is investigated whether the area reduction formula (F – f) is sufficient in the UTM system. The area on the ellipsoidal reference surface can be calculated by adding the reduction formula (F – f) to the area calculated from the map. Area distortion formulas (F–f) are defined from the surface of the sphere and the non-concave shape in UTM. In conformal transformation, the deformation value is directly related to the size of the area and its distance to the y-axis (Grossmann, 1976): (15) F−f=−fR2ym2 F - f = - {f \over {{R^2}}}y_m^2 where F, f, R, y_m represent the area on the ellipsoid reference surface, the area calculated from the map, the Gaussian mean radius and the average distance from the y-axis (prime meridian) of the parcel, respectively. Area reduction value (F – f) of application parcels in UTM was calculated using Equation (15) (y_m = 13 km, R = 6370 km). Area reduction values (F – f) of the parcels are given in Table 4.

The area reduction values calculated from Equation (4) in UTM were compared with the results obtained from the Karney method. In the light of the results, it has been observed that the results obtained from the Karney method are more reliable (Figure 4).

Figure 4

Therefore, in order to determine the area distortion in large areas, it is necessary to calculate the area with geographical coordinates on the ellipsoid reference surface or to transform to equal area projections.

The area difference that was calculated by using the Karney method and AEAC projection varies from 1 to 50 m². When comparing the parcel areas on the reference surface calculated using Equation (15) in UTM and AEAC projection, it has been found that there are differences from approximately 10 to 200 m² for each parcel (Figure 4).

If AEAC projection was chosen instead of UTM, a total of 800 m² distortion was observed in the total of all area distortion values in the parcels.

One of the important criteria in the area calculation is how the corner coordinates of the parcel are produced. In this context, the scale of the map where the corner coordinates of a parcel are produced is directly related to the area calculation of the parcel. Area calculations of the parcels, which cannot be measured in the terrain, are calculated on the scaled map. For this, the corner coordinates of the parcel are produced on a scaled map by digitisation technique. In this context, the larger the scale of the map produced in coordinate data, the higher the sensitivity values of the data produced on the map.

The area calculation performed using the parcel corner coordinates produced from the scaled map is formulated below: (16) 2f=∑i=1nxi(yi+1−yi−1) 2f = \sum\limits_{i = 1}^n {{x_i}\left({{y_{i + 1}} - {y_{i - 1}}} \right)} (17) 2f=∑i=1nyi(xi+1−xi−1) 2f = \sum\limits_{i = 1}^n {{y_i}\left({{x_{i + 1}} - {x_{i - 1}}} \right)} where n, f and (x,y) denote the total number of corners of the parcel, the area of the parcel on the map and the corner coordinates of the parcel. The parcel area is calculated using Equations (16) and (17). Then the root mean square error (RMSE) of the parcel area is calculated using following equations: (18) MF=mp22∑i=1nSi+1, i−12 {M_F} = {{{m_p}} \over {2\sqrt 2}}\sqrt {\sum\limits_{i = 1}^n {S_{i + 1,\,i - 1}^2}} (19) Si+1,i−12=(xi+1−xi−1)2+(yi−1−yi+1)2 S_{i + 1,i - 1}^2 = {\left({{x_{i + 1}} - {x_{i - 1}}} \right)^2} + {\left({{y_{i - 1}} - {y_{i + 1}}} \right)^2} (20) fscale=f±MF {f_{{\rm{scale}}}} = f \pm {M_F} where m_p and S represent the RMSE of the corner points and the square of the diagonal lengths of the parcel (Bogaert et al., 2005; Navratil and Feucht, 2008; Hejmanowska and Woźniak, 2009). Equations (18–20) are not available in GIS software. But users can add Equations (18–20) to the program with the help of the interface.

It is seen in Table 5 that as the scale of the map gets smaller and the area gets bigger, the value of the area correction due to the scale increases as a result of digitisation technique. For example, an area of 2500 ha, which is calculated using the corner coordinates produced from a 25,000-scale map, has been calculated as an area of 2.5 ha. When the parcels with large areas (forests, pastures, etc.) are taken into consideration, it is clearly seen that the area differences resulting from the use of medium-sized maps have high values. The areas calculated from the data obtained from the scaled maps using the digitisation technique are given in Table 5.

On the other hand, smooth square parcels are used for the application in Table 5. When using 15 forest parcels selected within the study area instead of these parcels, the results obtained are shown in Table 6.

When these forest parcels are examined, it is seen that the effect of digitisation technique used in scaled maps increases as the scale of the maps gets smaller. Therefore, large-scale map must be used for the area calculation with the digitisation technique in the parcels with large areas (Figure 5).

Figure 5

Today, area calculations are made on the surface of the ellipsoid or projection plane with a conform transformation from the ellipsoid. The ellipsoid selected is the reference ellipsoid determined in the geodetic datum. GRS80 is used in ITRF datum and Hayford ellipsoid in ED50 datum. The ellipsoid selected as reference is determined by the geoid. Ellipsoidal height (h) is measured in Global Navigation Satellite System (GNSS) measurements. The height type shown on the maps is the geoidal height. The transformation between these two types of heights is enabled by geoid undulation (N). The effect of the height factor in the calculation of a parcel area is determined using Equation (21): (21) Fh=F+F(1+2havgR) {F_h} = F + F\left({1 + {{2{h_{{\rm{avg}}}}} \over R}} \right) where F_h, F, h_avg denote the area of the parcel with a height from the ellipsoid, the area of the parcel on the surface of the ellipsoid and the average of the heights of the parcel corner points, respectively (Koçak, 1985). The slope factor is not taken into account in this calculation. Area correction of sample parcels with different heights and area values are shown in Table 7.

As seen in Table 7 and Figure 6, area deformations increase as the area and height values increase. When using 15 forest parcels selected within the study area instead of these parcels (for regular polygon parcels), the results obtained are shown in Table 8.

Figure 6

In the light of the results obtained in Table 8, it is observed that the area correction values increase as the average ellipsoidal height values of the parcels increase. Equation (18) that is formulated above is not available in GIS software. In applications that require high precision, GIS users should consider the height factor when using large areas of parcels.

In the calculation of the parcel areas, besides the height factor of the terrain, the slope of the terrain is an important factor to be taken into consideration. The effect of the slope factor on the area of the parcel with a certain height from the ellipsoidal reference surface is presented in Equation 22: (22) Fs=Fh+Fh1+(tan a)2 {F_s} = {F_h} + {F_h}\sqrt {1 + {{\left({\tan \,a} \right)}^2}} where F_h, F, tan a denote the area of the parcel with a height from the ellipsoid, the area of the parcel on the surface of the ellipsoid and the average slope of terrain, respectively. The changes caused by the effect of the slope factor on the area of the parcel are presented in Figure 7.

Figure 7

Area values calculated using scaled maps are not exact area values. In order to compute the exact area value of a parcel, different types of area corrections must be added to the area calculated from the scaled map. The impact of these area corrections on the area calculated on the scaled map is shown in Tables 5 and 6 as a percentage.

As seen in Table 9, area corrections calculated based on projection and scale factors have negative values and area corrections calculated based on other factors have positive values. The biggest correction to be added to the parcel area is the area correction caused by the slope factor.