Abstract
Global geodetic techniques currently can provide the user with worldwide millimeter accuracy. Preservation of this degree of accuracy in derived products is far from straightforward and may leave vast room for trouble in the different steps involved in the collection, storing, processing, analysis, and delivering of geospatial information. This paper is envisioned to serve as a guide for those utilizing map projections, in any possible form of application-cartography, GIS, remote sensing, photogrammetry, etc., to the common (and not so common) causes of error and misconception. This work also explores and questions the validity of some of approximations that are routinely implemented and quantifies the corresponding impact. These include the impact of neglecting meteorological corrections, reduction to ellipsoid and grid scale factors for distances, meridian convergence and arc-to-chord correction for angles, and mixing up with different frames and reference systems, height systems, or deceptively similar map projections. Correct indications are also given for accurately performing geospatial operations such as intersection of lines, determination of minimum point to line distance, and area determination for cadaster, which are often performed with suboptimal accuracy.
1. Introduction
The ever-increasing demands and capabilities of modern technologies and theoretical advances, with space geodetic techniques able to realize the International Terrestrial Reference Frame (ITRF) with accuracies of few millimeters [1], have brought many geodetic challenges that could be overlooked in the past. Now, users have to deal with moving (dynamic) coordinates, diverse local, regional, and global reference systems, as well as epochs to which measurements of different types are referred to, data sources of disparate origin and quality, and a long list of upsetting questions that, if treated improperly, may ruin their work. These issues become even more evident when one uses machine precision, open-source libraries such as GeographicLib [2, 3], such that solutions are no longer affected by numerical truncations. Another point to consider is the fact that in some applications, one can directly work in a 3D earth-centered earth-fixed (ECEF) coordinate system and avoid map projection distortion issues, whereas, for other applications, computations or visualization of data in ECEF can become overly complex and need to be simplified by working on a planar surface.
In this paper, we offer a set of useful facts, common misconceptions and approximations (with quantitative inaccuracies), and potential remedies as well as additional guidance for the users of cartography that may want to incorporate the current part-per-billion geodetic accuracy (i.e., some millimeters per thousand kilometers) into their work. Given that the list of possible topics to be covered is undoubtedly vast, we concentrate on those producing the most frequent misunderstandings (or the largest errors). We deal with the topics, inasmuch as possible, as a set of self-contained sections for the reader to easily identify their possible troublesome issue and the corresponding solution. As a succinct summary of facts and results that can be found in the manuscript, we mention(i)The existence of different terrestrial reference systems and frames (WGS84, ITRF, ETRS89, NATRF2022, etc.), some more accurate than others, some more stable than others, with different realizations for each one (with coordinates easily differentiating at the centimeter level or more) from which the user is normally obliged to select the one required by the work in question.(ii)The necessary correction of field distances by meteorological factors, the reduction to the ellipsoid, projection to a grid, which can easily amount to centimeter, decimeter, and decimeter differences, respectively.(iii)The required correction of field angles by meridian convergence and arc-to-chord angles, which can easily amount to differences of degrees and arc-seconds, respectively.(iv)The existence of different height types, with possible discrepancies up to the order of meters.(v)The existence of different map projections, some of which produce very different results despite their deceptively similar names.(vi)Guidance to proper handling of problems involving different grid zones, transforming between old and new maps, computing geospatial operations as intersections and minimum distances (which are often solved in approximate ways, easily yielding errors of the order of meters), and area determination in cadaster.(vii)New directions in cartography, including the optimization of map projections and the definition of new state plane coordinate systems.
2. Common Issues and Main Corrections
2.1. Reference Surface, System, and Frame
Surprising as it might be given their prolific use, it is not always clear what geodetic or geographic coordinates, namely, latitude and longitude, refer to. While it is obvious that they serve to locate a point on the surface of an ellipsoid, what may not be obscured to the user is the definition and properties of the ellipsoid (it might even be that a spherical Earth model serves to simplify the ellipsoid). Beyond the geometry of the ellipsoid, which can be specified by two parameters, some dynamical parameters such as its angular velocity and associated mass are also conventionally adopted. The reference surface is hereby defined. This is the case, for example, of the Hayford, GRS80, and WGS84 ellipsoids (the latter two having a minuscule discrepancy of 0.1 mm between minor semiaxes, which is negligible for most practical purposes).
However, knowing the particular ellipsoid is not enough unless additional information, such as the location of its center of mass and the orientation of its axes in space, is unambiguously defined. The set of models and parameters (including a particular ellipsoid as reference surface) to do so constitute what is commonly known as a reference system. At present, the International Terrestrial Reference System (ITRS) is the most accurate reference system available and therefore the only reliable support for performing scientific geospatial operations. Being its definition purely theoretical, including the idea that the system corotates with Earth introducing no horizontal net rotations by not being attached to any particular tectonic plate, the question of how it allows users to work in such reference system arises. A realization of the reference system permits one to do so, and in its simplest form, it is accomplished by means of a set of benchmarks with their corresponding coordinates and velocities in the system; this constitutes what is called a datum or reference frame. The International Terrestrial Reference Frame 2014 (ITRF2014) is the newest realization of the ITRS, which for the first time in the history of ITRS realizations includes not only coordinates and velocities for stations but also nonlinear motions [1]. At the time of this writing, it is expected that a new realization, ITRF2020, will be officially released in a few months [4]. Many localized reference frames (e.g., NAD83 (2011) or the forthcoming NATRF2022 for North America) are tied to ITRF at a specific epoch in time (e.g., 2008.00). Table 1 shows some of the most commonly used reference frames along with their corresponding realizations.
Users often find allusions everywhere to the somewhat mysterious WGS84 reference system. This is the U.S. Department of Defense (managed by the NGA) global reference system used for the development and maintenance of the Global Positioning System (GPS). Apart from the few stations comprising the GPS Ground Segment, there are no ground benchmarks to define the system that users can access, leaving the satellites and their broadcast ephemerides as the only possibility to utilize this system. Although the origin of their realizations is a bit different than those of the ITRS, at present the user may consider that the latest realizations of both systems are compatible for practical purposes. Also, compatible with an ITRS realization, say ITRF20xx, is an IGSxx, the corresponding reference frame used by the International GNSS Service (IGS) to derive their products, for example, precise satellite ephemerides, currently given in IGS14 which is fully compatible with ITRF2014.
Now, coordinate variability with time may be an unavoidable fact for scientific applications, but it is clearly not convenient for cartography, let alone for cadastral or land registration, where it is desirable to avoid as much as possible the global variability due to plate tectonics and focus on a particular area of interest. This is the reason for the definition and common use of regional and national reference systems and frames. They include the European Terrestrial Reference System 1989 (ETRS89), defined to be coincident with ITRS at epoch 1989.0 and fixed to the stable part of the Eurasian Plate, and the corresponding reference frames, for which ETRF2000 is conventionally adopted; the Geocentric Datum of Australia of 1994 (GDA94) coincident with ITRS realization ITRF92 at epoch 1994.0; China Geodetic Coordinate System 2000 (CGCS2000), which is coincident with ITRF97 at epoch of 2000.0, or the North American Datum 1983 (NAD83) whose latest realization is the NAD83 (2011) epoch 2010 for the North American Plate. In 2022, NAD83 will be replaced by the North American Terrestrial Reference Frame of 2022 (NATRF2022) which will be defined based on the latest IGS reference frame at an epoch to be determined [5, 6]. Additional reference frames will be developed for the Pacific, Mariana, and Caribbean plates. These frames will be related to IGS through an Euler pole rotation specific to each plate. This effort will also result in a redefinition of the State Plane Coordinate System utilized for many engineering and cadastral applications.
A fact which no user should disregard in principle is that differences between coordinates in the international and a regional frame (e.g., ETRF) are very significant. These differences may reach up to the order of meters; hence, it is the task of the cartographer to properly transform between the corresponding frames (see, e.g., [7] for the case of ETRF). Oftentimes, handy online coordinate conversion calculators may not properly perform these transformations or suffer from numerical transformation. When converting coordinates between reference frames or even within coordinate systems within the same reference frame, one should perform the conversion in reverse to ensure truncation has not occurred. Reference frame transformations may also be specific to a region (e.g., conterminous US), so the user should carefully check that the transformation being utilized is appropriate for their project sites.
Another potential source of confusion is where a user may inadvertently apply a transformation twice. For example, a user may process coordinates for a base station using a postprocessing service such as OPUS (The Online User Positioning Service) by the NGS. This system will provide coordinates in NAD83 (2011) epoch 2010.00 (as well as preview coordinates in NATRF2022) after a conversion from ITRF. When applying this coordinate to the base station coordinates to perform GNSS baseline computations to a rover, one rarely should reapply a reference frame transformation despite the software labelling of the coordinates as “WGS84” coordinates. Otherwise, the reference frame transformation will be doubled from what it should be because it was applied in the postprocessing service and then, again in software computing, the GNSS baselines to the rover positions.
In summary, when providing coordinates, one should be specific on the corresponding reference frame used as well as the epoch the coordinates are referenced to. Most likely, this will be an ITRF, e.g., ITRF2014, for scientific applications or, alternatively, a regional reference frame, such as ETRF2000 for Europe, NAD83 (2011) epoch 2010.00 for North America, and CGCS2000 for China. Furthermore, while a few decades ago it was rare to find the accuracies corresponding to the coordinates of the different benchmarks in official lists, this information is recently being provided by most of the national geodetic institutions and must be taken into account when it comes to calculate the resulting accuracies of the coordinates of derived points; that is, the actual coordinate accuracies of the points used to bring the work into the particular frame must be considered rather than assuming that these coordinates are exact.
2.2. Reduction from Terrain to Ellipsoid
Geodetic measurements (angles, distances, levelling differences, GNSS baselines, etc.) are acquired on the terrain, or more precisely at a short distance off it (i.e., the instrument height). Except for the cases where terrain measurements are themselves the final magnitudes of interest, they have to be transformed or reduced to the corresponding values over the Earth reference surface so that they can be properly used for geospatial operations. Since the 18th century, the ellipsoid of revolution, with slightly different values for different ellipsoid definitions, is the reference surface (or Earth model) for common reference of these measurements. The required operations to refer terrain measurements to the ellipsoid are known as reductions to the ellipsoid. While accurate consideration of all effects playing a part must be left to the expert in geodesy, it is worth outlining here the most significant corrections that have to be considered.
To start with, one should recall that measurements obtained with an electronic distance meter (EDM) are based on the propagation of a light beam (normally of infrared or red wavelength) between the EDM and the reflector and then back again to the EDM. The meteorological conditions existing in the moment of observation are normally different from the ones stored in the EDM as reference for computing the distance based on the manufacturer calibration. More precisely, the actual index of refraction is normally different from the reference index of refraction used by the EDM. A correction to the measured distance in terms of the existing temperature, especially, but also atmospheric pressure and humidity, must be applied before reducing the distance to the ellipsoid. Neglecting this correction may lead to significant errors in the distance, depending on the discrepancies of the on-site temperature, humidity, and atmospheric pressure from the EDM reference values. Approximate relationships between actual and reference differences in temperature, atmospheric pressure, and humidity values and the resulting errors in the measured distance are shown in Table 2 adapted from [8]. They show that neglecting the meteorological correction may easily have an impact on the centimeter level in a distance of several hundred meters. This is also true, e.g., for a terrestrial laser scanner (TLS), since its distance measurement follows the same principle. Note that not all scanners have this capability to apply corrections; however, some long range scanners now include onboard sensors to directly obtain these meteorological parameters.
The instrument user manual normally gives handy formulation to perform this meteorological correction, as well as the possibility to introduce in the field the real meteorological values so that the measurements are automatically corrected. Standard formulation to account for this correction (up to the 1 mm/km level) is also provided by [9]. For even more precise corrections, the user should follow [10, 11].
Once the distance between, say, points and , , has been corrected with meteorological values, it then has to be reduced to the ellipsoid. As it is represented in Figure 1, the main correction, equation (1), brings this distance to the chord going through the ellipsoid . A second correction, equation (2), brings this chord to the normal section onto the ellipsoid , which, except for distances of hundreds of kilometers, is coincident with the distance along the geodesic line (i.e., the shortest line onto the ellipsoid surface) with discrepancies below 1 mm.where is the radius of the normal section of the ellipsoid for the azimuth of the line or simply a mean Earth radius (e.g., 6,371,000 m) and and are ellipsoidal heights (see the subsequent section about heights) for points and , respectively. Alternatively, one can use an approximate method to directly reduce from ground to ellipsoid in terms of a mean elevation and mean geoid height with mean Earth radius [12]:
For the cases where reduction to the ellipsoid is overlooked, we can find significant errors that are very much proportional to the measured distance. They are strongly dependent on height; they are almost completely due to equation (1) and may amount to significant figures. As an example, see the computations in Table 3 and Figure 2 where, among the different results, it is shown that a distance error of 7.8 cm occurs for a distance of 1000 m measured at a height of 500 m above the reference ellipsoid.
In consequence, except for measurements at very minimal heights above the reference surface, at least reduction with equation (1), or equation (3), should be inexcusably done. Amstrong et al. [13] provided a detailed discussion of these factors and the development of low distortion projects to minimize differences between the terrain and the grid of the coordinate system.
Regarding angles, we can also find differences between the angles measured in the field and the ones referred to the ellipsoid. The plumb line or vertical line of force with which the instrument is setup on the terrain is different from the normal line to the ellipsoid surface. They form an angle called vertical deflection, which is usually of a few arc-seconds, but can amount to several tens of arc-seconds esp. in mountainous regions (due to irregular subsurface densities). Some models, such as the EGM2008 [14], permit the computation of this angle given the coordinates of a point on the ellipsoid. Correction of vertical deflection is important for vertical angles but also has an impact on horizontal angles as a second-order correction (possibly of one order of magnitude less). Nonetheless, this correction has to be considered at least whenever the instrument accuracy is of the order of the correction.
2.3. Projection from Ellipsoid to Grid
Map projections or grids are the source of inevitable distortions due to the different character in terms of curvature of a plane and a sphere or ellipsoid (just imagine trying to wrap a ball with some gift papers). Map projections are, however, of undeniable help for visualization as well as a tool for computation (mathematical formulas are much simpler on a plane compared with a curved surface).
Apart from the scale of representation, one has to account for the length distortion. This means that, in spite of the scale information given in the map legend, there is actually no map with a constant scale for the entire area depicted! For a pair of infinitesimally close points and , we can define the linear distortion coefficient as the ratio of the projected distance to the original distance on the ellipsoid surface ds and obtain after some derivations using differential geometry [15] such thatwhere and are the geographic coordinate differences between the infinitesimally close points so that , , and also represent partial derivatives (evaluated in point ) of the functions defining the map projection and with respect to the geographic coordinates latitude and longitude ; and and are the principal radii of curvature of the ellipsoid. For the sake of a better illustration, the reader is invited to use the tool available in [16] for the evaluation of linear distortion coefficients in the most common map projections.
For angle-preserving (so-called conformal) projections, is a function of position on the ellipsoid only (i.e., and finally cancel out after including the expressions for the partial derivatives). An average linear distortion coefficient for the line has to be computed and used for projecting the distance from the ellipsoid onto the grid :
Neglecting the use of the linear distortion coefficient equation (5) may lead to errors of hundreds of mm/km (e.g., 3 cm in just 100 m!).
It has to be noted that sometimes a handy scale coefficient that includes both correction for length grid distortion equation (5) and reduction to the ellipsoid for a mean height in the area equation (3) can be provided.
Referring to angles, they also change when transformed or projected from the ellipsoid surface to a particular map projection. Even if the projection is said to be angle-preserving (conformal), this does not mean that the geodetic north is aligned with the grid north (customary coincident with the -axis). They form an angle called meridian convergence so that geodetic azimuths and grid azimuths differ considerably. An additional minor problem, often overlooked, may be that a geodesic on the ellipsoid does not appear as a straight segment on a grid. We pay special attention to these angular issues in the subsequent sections entitled meridian convergence and the arc-chord problem in grids.
Finally, an expression analogous to equation (4) also exists for the relationship between grid areas and original areas on the ellipsoid. It has to be noted that for conformal projections, the resulting area distortion coefficient is simply the square of the linear distortion coefficient. Obviously, this correction is equally important and has to be taken into account when area matters (e.g., for cadastral purposes).
2.4. Heights
Classically, heights have been treated separately from horizontal coordinates (be these latitude and longitude geographic coordinates or grid coordinates). One should pay particular attention to the height in use, since there are different definitions of height (even for height above sea level) and countries do not always use a unique type of height for its territory, for instance, orthometric heights are used in the U.S. within the North American Vertical Datum of 1988 (NAVD 88) and dynamic heights within the International Great Lakes Datum of 1985 (IGDL 85). The most common types of height are(1)Ellipsoidal , with respect to the particular ellipsoid used as an approximation for the Earth. In consequence, ellipsoidal heights do not have a direct physical meaning but a geometrical one. They are the natural height used for satellite systems (e.g., GNSS).(2)Orthometric , with respect to the geoid (conventional equipotential surface of the Earth gravitational force that would coincide with the mean ocean surface extended below the continents). Its accurate determination entails the usage of gravity values on surface as well as a particular hypothesis about the vertical gradient of gravity in the subsurface until the geoid (the most common hypothesis leads to the technically called Helmert orthometric heights). This type of height is often referred as height above mean sea level and is commanded as the official height type in many countries (e.g., U.S. or some countries in Western Europe). Its relationship with ellipsoidal height is easy in terms of the so-called geoid undulation N, which can be obtained from a geoid model, e.g., EGM2008 [14],(3)Normal , when instead of hypothesizing about the gradient of real gravity, a formula in terms of the so-called normal gravity (i.e., the one generated by the ellipsoid of reference equipped with a mass and rotation rate) is used instead for the computation. It is therefore not referred to the geoid but to the quasigeoid. Similar to equation (6), a simple expression relates this height with the ellipsoid height in terms of the quantity named height anomaly :(4)Dynamic , which is the gravity potential difference with respect to the geoid divided by the value of the normal gravity.
Other slightly different flavors exist for heights (such as normal orthometric) as well as several issues in their practical implementation (e.g., zero surfaces conventionally used may have relative offsets, see, e.g., [17] for observed offsets between national height datums among European countries). Consulting the expert may often be the best advice for the complex cases. In coastal mapping applications, other reference datums for height exist such as mean low water (MLW) or mean higher high water (MHHL). NOAA NGS released a tool, VDatum, to convert between several ellipsoidal, orthometric, and tidal datums as well as several geoid models [18]. For details on these datums, please refer to https://tidesandcurrents.noaa.gov/datum_options.html.
3. Additional Issues and Misconceptions
3.1. The Arc-to-Chord Problem in Grids
One issue to always bear in mind is the fact that straight lines on a grid do not correspond to “straight lines” on the ellipsoid (namely, geodesics) nor on the terrain. The shortest path between two points on the surface of the ellipsoid, a geodesic line segment, appears, in general, as an arc in a map projection. The difference between this arc and the straight segment connecting its ends on a grid is more evident in angle than distance values (for all practical purposes, their length difference is negligible). Figure 3 shows the image of a geodesic segment (curved line), the chord on the grid, the cartographic north (or grid north, normally coincident by definition with the direction of the -axis), and the geodetic north (tangent in to the direction of the meridian).
While the meridian convergence, angle , may be of the order of degrees and is obviously nonnegligible (see next section dedicated to this question), the arc-to-chord correction, angle , which depends on the length, orientation and location of the line, and particular map projection, is considerably smaller and may be safely neglected for many applications, although it is always recommended to exactly compute its value following the corresponding formulation in the grid specifications. The different azimuths appearing in Figure 3 are the cartographic azimuth of the chord, , which is easily measurable in a grid projection given the images of endpoints and , the cartographic azimuth of the geodesic, (not easily accessible on a map projection), which fulfilsand the azimuth of the geodesic on the ellipsoid , which can be obtained by
Notice that the arc-to-chord correction depends on the location of both and , whereas the meridian convergence (see next section) is only point-dependent (on ).
A note of caution is worth here regarding Figure 3 and the corresponding relations in equations (8) and (9): they are only correct for conformal grids. We are depicting together in the same figure, and using their true magnitudes, angles on the ellipsoid and angles on the grid and , which is only possible if the grid is angle-preserving, i.e., conformal.
Just for illustrative purposes, Table 4 gives the arc-to-chord correction in Universal Transverse Mercator projection, Zone 30 N, for different length, orientation, and location of projected lines.
As we can see in Table 4, for lengths up to 10 km, the arc-to-chord correction amounts to a few arc-seconds at most and ordinarily could be neglected, whereas for distances of 100 km, it approaches 60″. Obviously, the situation worsens if this particular projection, which is optimized for middle latitudes, is used in low or high latitudes.
As a consequence of the arc-chord difference, the conservation of large areas in equiareal projections and preservation of angles in conformal projections are not always evident since the conservation occurs for the images of geodesics on the grid, which, as said, are not the straight segments between pin-pointed map locations but curved lines. This problem is in general (i.e., for every region in the grid) hard to avoid, whereas it is solved in practice if we are interested only in some particular areas by customizing a particular projection so that the enclosing geodesics appear as (very nearly) straight lines. The ellipsoidal gnomonic projection [19] is very well-suited to this particular problem. There remain, however, intractable problems by using a grid that should be solved on the surface of the ellipsoid [20].
3.2. Meridian Convergence
The direction of the meridian passing to one point is in general not parallel to the -axis in the grid. These directions form an angle called meridian convergence , which is needed to relate the azimuths in the grid with the azimuths on the ellipsoid surface, see the previous Figure 3 and equations (8) and (9).
Unlike the arc-to-chord correction , which can be usually neglected, the meridian convergence has to be always taken into account, since it can easily reach a few arc degrees. Its value depends on the location of the base point and the particular projection and is easily computed for conformal projections by (e.g., [21] p. 718)where we need the partial derivatives of the defining functions of the grid and .
Recall that Figure 3 and equations (8) and (9) are true for conformal projections so that subtraction of two cartographic azimuths, say and , which gives one angle in the grid, yields the same result as the subtraction of the corresponding geodetic azimuths, in this case and . In other words, meridian convergence cancels out when computing angles and angles (not azimuths) are preserved in conformal projections. We obtain this result by subtracting two expressions of the type of equation (9), one for and one for :
The preservation of angles on the ellipsoid to angles between chords in the grid is also true inasmuch as the arc-chord corrections can be neglected (see previous section). Otherwise, we will have to take into account the particular arc-chord values:
Just for illustrative purposes, Table 5 provides the meridian converge angle in Universal Transverse Mercator projection, Zone 30 N, for different locations of base points.
3.3. Geocentric Cartesian Coordinates (ECEF Coordinates)
The advent of spatial era has brought into play an additional natural system of 3D coordinates for positioning spatial aircrafts or expressing the measurements from and to them: earth-centered earth-fixed (ECEF) Cartesian coordinates . These and should not be confused with grid coordinates (approximately in the east-west and north-south directions) or horizontal coordinates in a local reference system inasmuch as ECEF does not represent the vertical component. Rather represents the distance along the axis from the center of mass of the Earth through the North Pole while the and axes are orthogonal and intersect the Earth at the equator. All of them (, , and ) are partly related to the horizontal plane as well as to the vertical direction.
There is a clear and straightforward geometrical relationship between these geocentric coordinates (, , and ) and geographic coordinates (latitude , longitude , and ellipsoidal height ):where is the radius of curvature of the ellipsoid in the east-west direction and a and b are the major and minor semiaxes of the ellipsoid. For the inverse computation, geographic coordinates in terms of Cartesian coordinates, there are a vast number of approximate methods, iterative methods, and also closed-form solutions (e.g., [22]).
3.4. Use Correct Formulation
Trivial as it may seem, the issue of applying the proper formulas causes a lot of trouble. One particularly common case is the confusion between Mercator and so-called Web Mercator projections. The first one is the famous conformal projection devised in the 16th century by the Flemish cartographer bearing its name, whereas the second one is a simplification commonly used for web mapping applications, such as Google Maps, that require a quick representation on-screen depending on the zooming value. The simplification, which apart from some scale coefficients, basically entails using the formula of the spherical isometric latitude instead of the correct, but more complex isometric latitude formula for an ellipsoid, makes the projection slightly nonconformal, which cannot be fully appreciated in the view displayed on the screen but can definitely be observed in computations. This projection is often termed as Pseudo-Mercator projection, which unfortunately makes things even less clear for the user.
Between the community of cartographers, it is quite common the usage of EPSG codes [23] to define the reference system, frame, and grid used for a particular application. Some of them may be deceivingly similar for a nonexperienced user but a careful look should reveal that, e.g., EPSG3857 coordinate reference system is “Pseudo-Mercator” using “spherical development of ellipsoidal coordinates,” i.e., Web Mercator with radius equal to the major semiaxis of WGS84 ellipsoid, whereas EPSG3395 is the standard Mercator grid for a spherical reference surface centered in latitude 0° and longitude 0°, which is often referred as World Mercator [24, 25]. The differences between them may amount to “errors of 0.7 percent in scale and differences in northing of up to 43 km in the map” [25].
Web Mercator (EPSG3857) and World Mercator (EPSG3395) projections are also different from Universal Transverse Mercator (UTM, EPSG258xx, where xx represents the zone number, e.g., EPSG25830), a particularization of the Transverse Mercator (none of the above two!) for zones of a width of 6 degrees in longitude, with a scale factor of 0.9996 in the central meridian (which diminishes the average scale distortion in the zone of interest) a false easting of 500,000 m and in the case of the southern hemisphere also a false northing of 1,000,000 m. The confusion between any of these different projections may make a complete mess of the task at hand!
3.5. Different Grid Zones
All successful alternatives to solve a cartographic problem (e.g., spatial intersection and buffering) need the unambiguous definition and use of a unique mathematical surface be it the ellipsoid surface, along with the tools of geometrical geodesy, or a particular grid along with the tools applicable in map projections. While trying to solve a particular problem in the ellipsoid of reference may be challenging due to the complex formulation required for this curved surface, the analysis of the problem in a grid involves only the use of planar geometry and consideration of the corresponding grid distortions. When the problem at hand comprises points located in different grid zones in a projection with zone division, such as UTM, one should choose one common zone (advisably the one that has more data or contains a larger portion of the project) to project all points and solve the problem therein. The only costs to pay are (1) the increase in distortion values, which, if are accurately known and accounted for, suppose no problem for computing the final result, and (2) the need to use an algorithm with extended precision, e.g., [26] for UTM. For problems involving very large areas, however, a projection with zone division should be abandoned altogether in favour of a more convenient projection (see next section).
3.6. The Best Projection
In many occasions, the particular map projection to be used is prescribed and unchangeable (e.g., official map projections in Europe [27]). In some other cases, it can be user-selected either for mathematical convenience (e.g., to obtain geodesics represented as straight lines [19]) or for the sake of a better interpretability (e.g., the snake grid [28]). The expert has to evaluate the convenience of applying one of the thousands of existing map projections or deriving yet a new one.
A general recapitulation of some of the most used projections can be the following (largely adapted from [29]):(i)For the representation of the entire world (atlases), Robinson and Winkel-Tripel projections, as have been customary used by the National Geographic Society, or the recently presented Equal Earth projection [30].(ii)For topographic maps, conformal projections such as Mercator (in low latitudes), Transverse Mercator or Lambert Conformal Conic (in mid latitudes), and stereographic projection (high latitudes).(iii)For equal area representations, Lambert Azimuthal Equal Area and Albers Equal Area Conic projections.
We will not extend further on this question and simply refer the reader to the excellent manuals written on the topic, e.g., [31, 32].
3.7. Accurate Geospatial Operations in Spatial Database Management Systems
Geospatial operations such as the intersections of geometries and computation of areas of influence (buffering) entails solving some basic geometric calculations such as (1) calculation of a point from a starting point plus azimuth and distance (also known as the direct problem of geodesy), (2) determination of distance and azimuth between two points (also known as the inverse problem of geodesy), (3) identifying the location of line intersection, and (4) calculation of the minimum distance from a point to a line. These problems, which need to be solved on the surface of the ellipsoid, have often been resolved by making use of auxiliary projections, which entails some of the aforementioned problems (e.g., difference between a straight line in the grid and the geodesic line), therefore producing solutions of low accuracy. In recent years, a trend has emerged that all computations should be performed without any detectable accuracy degradation, that is within accuracies close to machine precision [2, 19, 26]. In consequence, the most powerful spatial analysis software solutions (Oracle, Google Earth Engine, PostgreSQL, and PostGIS) have started to implement some of these geometric calculations directly on the ellipsoid, namely, the direct and inverse problems of geodesy, while the intersection of lines and minimum distance from a point to a line still lack an adequate treatment such that significant errors (up to hundreds of meters or more in the worst cases) can be committed [33]. These issues show how necessary it is to fully implement close to machine precision geometric computations on the surface on the ellipsoid for a reliable spatial database management system.
3.8. Old Maps and New Maps
This problem of combining information from maps of different vintages involves, in general, several of the issues mentioned above, for example, different reference systems, different map projections, and different grid zones, as well as the inherent inaccuracies to each of the map elaborations and levels of detail.
Transformations are in general performed in geodetic coordinates (i.e., latitude and longitude for the reference ellipsoid), whereas rule-of-thumb solutions (e.g., constant and shifts between UTM coordinates in ED50 and ETRS89 systems) are normally not advisable. The disciplines in which this question arises range from cadastral and cartographic updating to remote sensing or photogrammetry, see, e.g., [34–36]. For the case of transformation between old and new frames and systems, distortion grids are normally used, see, e.g., [37]. Mapping agencies have developed official tools for this purpose, e.g., NOAA NGS has developed a converter between the reference system NAD83 and the old NAD27 [38]. The first implementation was called NADCON, which has now been integrated into the NGS Coordinate Conversion and Transformation Tool (NCAT). NCAT also enables conversions between several realizations of NAD83. Additional tools have been developed to model tectonic movement such that coordinates can be estimated for other epochs than those used for the acquisition. Examples would be the Horizontal Time Dependent Positioning Service offered by NGS [39] or the Land Information New Zealand coordinate converter [40]. The new reference frames (e.g., NATRF2022) for the U.S. will incorporate these models in the supplemental online tools. Despite these many tools and resources available, the best advice to be given, due to the number of different issues involved, is here more than ever to consult the expert.
3.9. Accurate Area Determination in Cadaster
In contrast to a theoretically well-defined polygon on the ellipsoid surface, the legal definition of a parcel is made at the actual height of the land, which is the area of interest for everyday use: land cultivation, cadastral income taxes, agricultural subsidies, etc. While it is common to assume that this area can be regarded the same as the one projected on the horizontal plane at the mean elevation of the parcel [41, 42], the differences between the horizontal area at the average elevation of the parcel and the corresponding area on the surface of the ellipsoid can be substantial, easily reaching hundreds of ppm. The interested reader is referred to [42] for a detailed analysis of the different possible sources of errors affecting the resulting area (due to ignoring elevation, disregarding the map distortion factor, the simplification of geodesics as straight lines, rounding effects, etc.) as well as for the description of an accurate method for parcel area computation.
Finally, it is also worth mentioning that, as [43] demonstrated, even though the parcel area is usually rounded to a square meter, the correctness up to this order cannot be guaranteed for parcels with large or complex boundaries if point coordinates are determined up to the centimeter level, a fact that should be acknowledged when defining the technical requirements for cadastral works.
3.10. New Directions in Cartography
Far from being a science with unshakable principles, cartography is constantly evolving, whether to preserve the unprecedented level of accuracy achievable with today’s technology, to meet the fast-response demands of online viewing on the Internet, or to facilitate the calculations of a global, ever growing, community of users (of which professional cartographers and geodesists make up only a tiny fraction).
It is virtually impossible to account for all the recent advances, but some areas worth noting in which considerable effort is being made include(i)Web mapping applications: an adaptation of the classic Mercator projection with the focus on fast online visualization, which is known as Web Mercator, has been proposed in the last years. Although it has brought additional challenges and discussion [44], it is currently used by the majority of online map providers (Google Maps, Bing Maps, OpenStreetMap, etc.), which face the inevitable challenge of efficiently handling the immense volumes of available geospatial data [45].(ii)Augmented reality (AR), virtual reality (VR), and 3D maps: several challenges are yet to be resolved, among them the lack of standards for user interaction with the map [46] and an even greater improvement of the user experience for scientific and societal applications [47]. Also note that most of these works are being done in basic local 3D Cartesian coordinate systems and not really considering distortion, projections, and the corresponding coordinate systems when bringing in input data. While this could be a successful approach for small sites, certainly it is not for large areas.(iii)Visualization of uncertainty in geographic data: a correct interpretation of results is often strongly dependent on the uncertainty of the values depicted. Among the most successful proposals to represent uncertainty on bivariate maps using visual variables, one may find boundary fuzziness and color lightness [48]. The visualization of uncertainties in 3D models can be performed by blending multiple colors with shaders as in [49], which also allows analyzing different error components (e.g., angular, range, horizontal or vertical components).(iv)Optimization of map projections in terms of ground-to-grid distortions, which can be more meaningful to the user than the usual definition in terms of ellipsoid-to-grid distortions, including the definition of low distortion projections (LDPs). The reader is referred to more specific sources, such as [13, 50–56].(v)The forthcoming release of the 2022 United States State Plane Coordinate System (SPCS2022) deserves special attention, see, e.g., [54, 57], as part of the transition from the North American Datum of 1983 (NAD83) to the 2022 Terrestrial Reference Frames.
4. Conclusion
Map projections are a useful tool not only for visualization purposes but also to facilitate geospatial operations that become much simpler and computationally efficient on a plane than on the Earth reference surface. Unfortunately, the simplification comes at the cost of the unavoidable appearance of distortions. Beyond the well-known, and visually evident, deformations in small-scale maps, this paper presents a more thorough analysis of distortion in maps including numerical quantification of common approximations and misconceptions from the fact that any map cannot have a constant scale in its entire domain due to unsolvable issues with a map projection.
Data Availability
The data used in this study are available upon request to the corresponding author.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The first author wants to thank Fernando Sánchez Miramón at Institut Cartogràfic Valencià for fruitful discussions at the crossroads between geodesy and cartography.