DETERMINATION OF AN APPROPRIATE PROJECTION SYSTEM FOR FOREST AREAS IN TURKEY

. Geometrical surfaces such as sphere and ellipsoid are considered as reference surfaces since there is no geometric shape that perfectly represents the earth when translating the earth into a map plane. Hence, on 3D reference surfaces, it is almost impossible to perfectly preserve the angle, direction and area properties and transfer them to a map plane without any deformations. The scaled topographic maps produced in our country under provision of map production regulations are conformal projections that do not preserve area properties but angle and shape properties. Area values calculated by projection coordinates cannot be considered the exact area values therefore, an area reduction is needed. Area values calculated by ignoring this situation in GIS based software do not represent the accurate area values on reference surfaces. The aim of this study is to determine the best area preserving projection for GIS applications in which area values are important. In this study, the real area values of 25 large-extent forest parcels are determined by employing the Danielsen method with geographical coordinates on ellipsoid surface. These parcels are also calculated by using the area-preserving projections available in ArcGIS software and are compared to their real area values.


Introduction
Nowadays, geographic information systems (GIS) -based applications use coordinate data as well as area data. Accuracy of a calculated area depends on the position accuracy of the coordinates used, map scale and selected map projection. The map scale and projection are the most effective ones in accurate computation of an area. Today's technological plotting, computation techniques and GIS software enable the conduction of applications in 1:1 scale, which attaches importance to the accuracy of the used coordinates and selection of a map projection.
Map projection is the reduction of measurements obtained from physical earth surface to a geometrically defined reference surface (such as ellipsoid, sphere etc.) and projection onto a map plane by means of geometrical and mathematical techniques. Projection from a 3D reference surface to a 2D plane (i.e. map) causes deformations in the area, angle and length of original shapes. It is not possible to eliminate all these deformations. Hence, one of the area preserving, angle preserving or length preserving projections is selected to preserve one of these deformations (Pearson, 1990).
forests with terrestrial laser scanning techniques (Guang et al., 2013), to divide endangered species into ecological zones in Canada (Committee on the Status of Endangered Wildlife in Canada [COSEWIC], 2015), to planning European forests (Borkowski, 2016), and to monitor forest inventory of 50 states in conjunction with forest data published in the USDA Forest Service Field Inventory and Analysis Library (Vogt & Smith, 2017).
In our country, boundaries of forest parcels are measured in the Modified Universal Transversal Mercator Projection (MUTM, 3° wide longitude zones) through Global Navigation Satellite Systems (GNSS). Areas of forest parcels are calculated by the Gauss Area Method using MUTM projection coordinate data. The studies carried out by the Regional Directorate of Forestry revealed the fact that growing stock of forests is increasing year by year. With the completion of the forest cadastre, it will be possible to determine the exact forest area with high accuracy. While the forests covered an area of 20.2 million ha area in 1972 in Turkey, this amount increased up to 22.3 million ha in 2015, which makes about 28.6% of the whole country. In parallel to this increase, timber wealth increased from 0.9 billion m^3to 1.6 billion m^3in the same time period (Republic of Turkey General Directorate of Forestry, n.d.). In our country, 1/1000, 1/5000 and 1/10000 scale maps are produced in the Modified UTM projection system, whereas 1/25000, 1/50000, 1/100000 and 1/250000 scale maps are produced in the UTM (6° wide longitude zones) projection. The UTM projection is obtained from the Gauss-Kruger mapping. The Gauss-Kruger projection is a mapping of a reference ellipsoid of the earth onto a plane and some definition of the ellipsoid and various associated constants are useful. In the Gauss-Kruger mapping; x prime axis is the equator and y prime axis is the prime longitude. Turkey is located between 26-45° east longitudes and 36-42° north latitudes. This 19° longitude difference is divided into four 6° zones (zone numbers; 35, 36, 37, 38) in the UTM, whereas into seven 3° zones (central meridians 27°, 30°, 33°, 36°, 39°, 42°, 45°) in the MUTM (See Figure 1).
The UTM system is a conformal projection as it is derived from the Gauss-Kruger mapping. The area deformation in the UTM projection system depends on the size of the area and the distance from the prime longitude. Therefore, geographical coordinates are needed for the accurate determination of areas, rather than UTM coordinates. However, existing geodetic and GIS software are not able to compute areas using geographical coordinates, which forces the users to seek for alternative software and also to make extra calculations. Instead, area preserving projection is usually used to determine areas (Yıldırım & Kaya, 2008;Yıldırım, 2012).
As the forest areas calculated in the conformal projection system increase, they less represent the actual areas. In our country, forests cover a considerable amount of area. They also consist of large-scaled parcels. Hence, it would be more reasonable to determine the exact amount of forest areas with area preserving projections instead of conformal MUTM projection. Thus, selection of area preserving projections in large-extent GIS applications carried out by the General Directorate of Forestry will reduce the deformations in parcel areas.

Applications requiring parcel area information (surface area) in general directorate of forestry
The parcel areas in our country are calculated by Gauss area method using MUTM coordinates. Since the Gauss area method is based on map coordinates, it is applied on planes and hence, is not applied on original ellipsoid surface. Forest area information is used directly in many studies carried out by the General Directorate of Forestry. Direct or indirect studies on forest areas in the central and provincial organizations of the General Directorate of Forestry: -Forest cadastre applications; -Studies on rehabilitation and afforestation areas; -Planning of forest and forest resources; -Production and marketing of wood and non-wood forest products; -Silviculture applications; -Retained carbon quantity in forest areas; -Studies on recreation areas; -Studies on permissions given for forest areas; -Pasture improvement and erosion control studies; -Forest management planning; -Studies on coniferous and deciduous species in coppice and high forests; -Growing stock studies; -Studies on regeneration areas. Calculated forest areas play a significant role in forestry organizations, especially in forestry economy applications. In forestry organizations; plans, practices and outputs are mainly based on forest areas. Considering the extent of forest areas in forest organizations, accurate calculation of forest areas is a necessity to maintain the sustainability of forests and conduct aforementioned studies related to forest management.

Area calculations
Areas, which are calculated from the scale maps generated from the ellipsoid reference surface by the map projection method, have plane geometry. On the other hand, area Figure 1. UTM and MUTM zones for Turkey deformation increases as the area grows. Hence, accurate value of an area should be calculated from the ellipsoid surface by means of the geographical coordinates. Area deformation formulas (F-f) are derived based on each projection type. The magnitude of the deformation in these formulas is usually proportional to the size of the area and the distance to x or y axis (Grossmann, 1976).
In the area deformation formula for MUTM; F is the original ellipsoid area, f is the map projection area, R is the Gauss average radius, and y_m is the average distance from the parcel to the y axis origin (prime meridian). In GIS applications, deformation formulas and geographical coordinates on the ellipsoid surface are not available for area calculation. Therefore, in such applications, an area preserving projection should be preferred, considering the geographical location and size of the study area, instead of conformal projections. In such cases, out of the areapreserving projections, the one that causes the minimum area deformation should be preferred. The most straightforward way to achieve this is to compare the area values on area preserving projections against the area values calculated from the ellipsoid geographic coordinates. In spite of the fact that the area-preserving projections are based on the principle to preserve area information, area deformation is inevitable to some extent (Snyder, 1987;Bugayevskiy & Snyder, 1995).
Many methods have been reported in the literature to calculate the area of any parcel using ellipsoidal geographical coordinates. Some of the most commonly used methods are; Kimerling (1984), Danielsen (1989), Gillissen (1993), Sjöberg (2006), Karney (2013) and Tseng (Tseng et al., 2015). Since the Kimerling method is a spherical solution, the edges are not geodetic line but the great circle in the sphere. For this reason, the accurate ellipsoid is not the solution. The Danielsen, Sjöberg, Karney and Tseng solution methods consider parcel edges as geodetic lines on ellipsoid and calculate the areas from the ellipsoidal area between this line and equator. The Gillissen method uses the Albers Equal Area Conic Projection to calculate the area of each part of the geodetic area divided by secants on the ellipsoid surface.

Application
Turkey is a large-scale country that extends in the East-West direction and its geographical boundaries do not fit into a single UTM gore (Figure 1). In Turkey, there are 28 Regional directorates of forestry in the General Directorate of Forestry. The Kütahya Regional Directorate of Forestry, which was chosen as study area in this study, is one of these directorates. The Kütahya Regional Directorate of Forestry, which covers an area of approximately 2 million hectares, consists of 6 forestry operation directorates and 40 forest sub-district directorates. The forest area covered by the Kütahya regional directorate of forestry is approximately 708,784.40 hectares, which is about 36% of the total area of the regional directorate of forestry in Turkey ( Figure 2).
In this study, 25 forest parcels from the Kütahya regional directorate of forestry were used, whose areas ranged from 339 ha to 2857 ha (see Figure 2). As a first step, the areas of these forest parcels were calculated by using the MUTM coordinates in the ITRF datum. Then, the area deformation of these parcels were calculated by using the Equation 1 ( m y = 13 km and R = 6370 km). Finally, on the reference ellipsoid (GRS80) surface, the MUTM coordinates of these parcels were converted to geographical coordinates and accurate area values were calculated using the Danielsen method (Ellipsoid Geographical Area) (see Table 1).
The area deformation rates determined by two different methods (see Table 1) are shown in the Figure 3.
As seen in Figure 3, reduction rate from map coordinates (F-f) increased as did the size of the area, whereas the differences from the accurate area were found to vary with respect to the shape of the parcels. Comparisons between the actual area values and areas computed via the (F-f) reduction equations revealed the insufficiency of the reduction equations. The differences between the actual areas of the parcels and areas computed from the reduction equations were found to vary between 1 m 2 and 200 m 2 , which made a total area difference of 1052 m 2 in all parcels. Deformation formulas; it does not match the actual field values when the area grows and is applied  to smooth square or rectangular external shapes (concave polygon). Therefore, in large areas, it is necessary to compute the accurate area values by using the geographical coordinates on the ellipsoid surface, which minimizes the amount of the area deformation determined from maps. Since the existing GIS software do not have the capability to calculate actual area values with reduction equations or geographical coordinates, six different area-preserving projections (see Table 2) implemented in the ArcGIS 10.3 GIS software were used in this study to compute the parcel areas. The projection parameters used for the parcel areas are given in Table 2 (Kennedy & Kopp, 2000).
The areas of the forest parcels given in Table 1 were recalculated by using the area-preserving projection methods given in Table 2. The differences between these areas and Danielsen-derived areas are given in Table 3. Table 3 was used to determine the area-preserving projection that achieved the minimum area deformation value.
The area deformation rates of the area preserving projections are shown in Figure 4.
As shown in Figure 4 and Table 3, the AlbersEAC, LambertEAA and BonneEA were found to achieve the maximum area deformations of 65 m 2 , 71 m 2 and 79 m 2 , respectively. The difference between the actual area of the  forest parcels and AlbersEAC -and LambertEAA -derived parcels were below 30 m 2 at a rate of 96%, whereas the difference with the BonneEA projection was below 30 m 2 at a rate of 92%. The areas of 21 parcels in the Al-bersEAC, 20 parcels in the BonneEA and 20 parcels in the Lambert EAA projection were found to have a maximum of 15 m 2 difference with the actual parcel areas. In addition, 6 parcels in the AlbersEAC, 8 parcels in the Bon-neEA and 6 parcels in the LambertEAA projection were found to have an accuracy of 1 m 2 . Since the AlbersEAC is a conic projection, it is a suitable projection method for Turkey. Conic projections are suitable for countries with a rectangular shape extending in middle latitudes and in the east-west direction. This, of course, enables the production of maps with a single coordinate system and less deformation all around the country. Currently, the MUTM projection is used for forest parcels in Turkey.
In case of using AlbersEAC instead of this projection, the MUTM in Table 1 and the AlbersEAC area deformations in Table 3 are shown in Figure 4. Figure 5 revealed the fact that deformation differences caused by the selection of the MUTM and AlbersEAC projection were found to vary between 1 m 2 and 250 m 2 . The total area difference for all parcels was found about 0.1 ha (1000 m 2 ). All forest areas of Turkey (22 342 935 ha) if we take into consideration; the total area deformation resulting from the projection selection is approximately 77 hectares. This value equals approximately 110 soccer fields.

Conclusions
In forestry operation directorates, forest areas play an important role in many studies carried out. In the Forest management, plans, practices and outputs are largely shaped by forest areas. When considering the size of forest areas in the forest general directorate, the accurate value of forest areas in terms of forestry management should be calculated with high precision. In addition, forest areas are very important for sustainable forest management and for the healthy implementation of plans and activities. All mapping applications including forest parcels in our country are produced in MUTM projection.
Projections can not represent the original surface (ellipsoid and sphere) perfectly. The angles (conform), distances and areas of shapes on an original surface are distorted when they are projected onto a map. Projection is done by preserving at least one of these three deformations. The MUTM projection is a conformal projection that is more deformed in terms of length and area. The area deformation in the MUTM projection varies according to the size and shape of the area, and to the distance from the prime longitude. Area deformation can be minimized by using geographical coordinates on ellipsoid surface instead of MUTM coordinates, by using projection area reduction equations (F-f) or by using one of the area-preserving projections. In this study, 25 forest parcels with areas ranging between 339 ha and 2857 ha were used to decide which method determines the parcel areas with minimum deformation. The accurate area values of these parcels were calculated by Danielsen method using the geographical coordinates on ellipsoid surface.
For forest application parcels, it can be seen that there are differences up to 200 m 2 when comparing the reduction values calculated from the area reduction equations and the accurate area values in the MUTM projection. The sum of these differences in all parcels was found 1052 m 2 .
Area reduction equations are insufficient to determine the areas of large and concave-shaped parcels in the MUTM projection. Hence, exact area values should be computed with geographical coordinates. However, existing GIS and CAD-based geodetic software are not capable of computing areas with geographical coordinates, hence; additional software and workload are needed for such purposes. Workload can be reduced by using area-preserving projections.
In this study, the area values calculated with six different area-preserving projections available on ArcGIS 10.3 GIS software were compared against the areas computed using the Danielsen method with geographical coordinates to determine which projection provides the minimum area deformation for Turkey. The maximum deformations of 65 m 2 , 71 m 2 and 79 m 2 were found with the AlbersE-AC, LambertEAA and BonneEA projections, respectively. The AlbersEAC and LambertEAA projections were found to be the ones that achieved the minimum deformations (360 m 2 ) in terms of the total area of all parcels. Since the AlbersEAC is a conic projection, it is suitable for Turkey. It is a well-known fact that the conic projections are suitable for countries in middle latitudes with a rectangular shape and lying in the east west direction. This also enables the production of maps with a single coordinate system, minimizing area deformations.
In our country, it was observed that the area deformation differences resulting from the selection of the MUTM and AlbersEAC projections varied between 1 m 2 and 250 m 2 . For all application forest parcels (28840 ha), this difference was approximately 1000 m 2 (0.1 ha). Considering all forest parcel areas in Turkey, the total deformation area was found to be about 77 hectares.
Since the AlbertEAC projection is not used in Turkey, area deformation differences are expected to be very high, considering the total forest area in Turkey. These differences can be considered too much for forestry departments with very high economical outcomes.
Forest areas are calculated on the elipsoid surface using the geographical Coordinates. Parcel corners have different elevation values in areas calculated with geographical coordinates on ellipsoid surface. In this context, accurate area value needs to be analysed and investigated by considering the topography of the area. Plenty of attention should be paid to interpret the actual area in sloping terrains. This was ignored in this study.
The optimum projection selection is important in GIS applications where accurate area information is needed. In addition, the users should pay attention to choose the optimum projection parameters (central meridian, initial longitude or central meridian, one or two standard parallels) in the areas close to study area to minimize the deformations. In this work, the errors due to the position and scale accuracy of the coordinates were ignored. In large forest area GIS applications, accuracy will be further reduced when the effects of errors resulting from the production of coordinates are taken into consideration. The Albers EAC projection proposed for Turkey may not be suitable for other countries. For other countries, one of the area preserving projection types according to the geographical position and size of the country should be selected first. Thus, area values calculated using the geographical coordinates of the parcels should be determined by comparing them with the results obtained using the danielsen method.