Direct georeferencing with correction of map projection distortions for active imaging Zdeněk Švec Department of Geomatics, Czech Technical University in Prague sveczde1@gmail.com Abstract In aerial photogrammetry, the Cartesian coordinate system for the description of object space is commonly used. In contrast, many projects have to be processed in specific map projection and vertical datum. In that space, some geometric de- formations exist. There are some compensative methods for active and passive sensors. In the case of active sensors, decomposition and the correction of ob- servation vector for each ground point can be used. We obtain height, horizontal distance and horizontal angle in this process. All of these values should be cor- rected for precise georeferencing. The contribution deals with the derivation of the corrections and gets some theoretical values from the area of the Czech Republic. Esspetially in the case of high flying heights the corrections may gain values even in order of meters. Keywords: Sensor Orientation, National coordinates, Mapping, LiDAR, GNSS/IMU Introduction The task of georeferencing in the field of aerial topographic survey is the determination of geometric relations between captured data and the real world [7]. It includes two consecutive procedures:the determination of the exterior orientation parameters (EOP) and the restitution scene from EOP and observed data [9]. EOP comprise of the spatial position of the sensor projection center and sensor attitude at the time of the observation. EOP for passive sensors can be gathered indirectly (by the measurement of image coordinates of ground control points) or directly from records of the on-the-board navigation system. On the contrary, active imaging is dependent (due to the sequential measurement principle andthe motion of the carrier vehicle) on direct methods. Referring to the direct georeferencing EOP is typically measured by the GNSS and inertial measurement unit (IMU). GNSS provides absolute position with sufficient frequency (at least 1Hz) and IMU sensor attitude and acceleration. Final trajectory is deduced by combining GNSS and IMU observations [1]. The Aerial survey product is mostly provided in “national coordinates”. It means the co- ordinate system realised by a combination of the national geodetic datum with a national map projection with the associated vertical system. The model for direct georeferencing is designed for cartesian space but national coordinates do not fulfill the condition and therefore cause various geometrical distortions. There are two ways to obtain accuratedata in national coordinates (See Fig. 1): restitution of the scene in cartesian (usually Earth-fixed) frame and Geoinformatics FCE CTU 15(1), 2016, doi:10.14311/gi.15.1.3 35 http://orcid.org/0000-0002-6414-6696 http://dx.doi.org/10.14311/gi.15.1.3 http://creativecommons.org/licenses/by/4.0/ Z. Švec: Direct georeferencing with correction of map projection transforming the scene to national coordinates or restitution of the scene in national coordi- nates with the corrected observation vector. Further text refers to active imaging, especially to georeferencing of LiDAR data. Figure 1: Accurate georeferencing in national coordinates for LiDAR. Direct Georeferencing for active imaging Reference frames and EOP transformation The EOP should be transformed into desired reference frame. Overview of reference frames purveys Tab 1. Detailed description of transformations can be found in [9] and [6]. Brief summary contain Fig. 2. Table 1: Overview of the reference frames frame description e Earth-centered earth fixed frame realized by International Terrestrial Frame n Eccentric earth fixed frame of national ellipsoid l Tangent frame of national ellipsoid p Projection frame established by national map projection b Body frame realised by acceleromerers of IMU Model for georeferencing of LiDAR data According to [10], the coordinates of ground point XG can be (in cartesian frame) expressed as XG = XEO + REORscanXrange = XEO + Xdg, (1) where XEO and REO are sensor position vector and the rotation matrices formed by angular EOP, respectively, Rscan and Xrange are scan angle matrix and range vector, respectively. The Geoinformatics FCE CTU 15(1), 2016 36 Z. Švec: Direct georeferencing with correction of map projection Figure 2: Schema of EOP transformation second term forms the observation vector for direct georeferencing Xdg. If the georeferencing is carried out in national coordinates (p-frame), XEO can be given by exact formulas (usually provided by state mapping authority), contrariwise Xdg is skewed due to datum scale change and p-frame geometry. We may assume that the correct position of ground point in p-frame X ptrans G is given by georeferencing in e-frame and subsequently rigorous transformation to p-frame. Then we gain the correct observation vector in p-frame Xṕdg as X ṕ dg = X ptrans G −X p EO. (2) Our task is to modify vector Xpdg to X ṕ dg. It will be to apply processes published in [10] as a practical approach. It involves some simplifications and refers to conformal map projection. Correction of p-frame distortions According to [10], the method consists of four subsequent steps (Fig. 3). Datum scale distortion. The Cartesian e-frame and n-frame have a different scale (if they have not used the same datum). Hence the length of the observation vector is different as well and it should be multiplied by scale factor mdatum. Decomposition of Xpdg to height component Zdg (it is always negative), horizontal distance D and horizontal angle ϕ (Fig. 4). It will be called “spatial observations”. Application of map projection corrections to spatial observations. The Earth curva- ture correction hec is added to Zdg, D is transformed to the geodesic distance S and to the projected length D´, skew-to-normal correction ζ and arc-to-chord correction δ Geoinformatics FCE CTU 15(1), 2016 37 Z. Švec: Direct georeferencing with correction of map projection is added to horizontal angle (normal-section-to-geodesic correction ξ is not taken into account, because it is numerically insignificant) Composition of map-projected observations to Xṕdg Figure 3: Sequence of correction of observation vector in p-frame IMU works based on Newtonian laws and its Z-axis is aligned with plumb line. If the vertical datum is not gravity-related, the gravimetric correction is required to be added in observations [5]. Decomposition of Xpdg is executed by simple equations Zdg = Zdg,D = √ X2dg + Y 2 dg,ϕ = tg Xdg Ydg · (3) The Earth curvature correction Height and length component are deduce from geometry in Fig. 5 and Fig. 6. hec can be expressed as hec = cos α |GP| = cos α(|OF| 1 cos α −|OG|), (4) G and F can be approximately regarded as locating at a same arc with the radius R + HG, then hec = cos α[(R + HG) 1 cos α − (R + HG)], (5) due to the small value of α we can simplify cos α≈1 − 12α 2 hec≈ α2 2 (R + hG), (6) α = atan D R + HS −Zdg ≈ D R + HS −Zdg , (7) where HS is the sensor projection centre. By combining (6), (7) we can obtain Geoinformatics FCE CTU 15(1), 2016 38 Z. Švec: Direct georeferencing with correction of map projection X axis Y axis Z axis G S D Xdg Ydg Zdg Φ Figure 4: Decomposition of observation vector according to [10] hec = D2(R + HG) 2(R + HS + Zdg )2 , (8) HG = HS + Zdg, (9) hec = D2 2(R + HS + Zdg ) · (10) Correction of length component The first step is the conversion of the horizontal distance D to the geodetic distance S. S is always much shorter than the radius of reference sphere, therefore we can approximate the geodesic line by a circular arc. By using (7) we can form a relationship S = Rα = Ratan D R + hS + Zdg , (11) from the definition of the map projection scale factor we obtain the equation for calculation projected length D́ = mS = mRatan D R + hS + Zdg , (12) where m stands for the map projection length factor. For this purpose it is sufficient to compute m in one point (preferably in the position of sensor projection centre). Geoinformatics FCE CTU 15(1), 2016 39 Z. Švec: Direct georeferencing with correction of map projection D HG hec D′ HS Zdg HG S F P G G′ map projection plane reference sphere α Xdg X′dg R G0 Figure 5: Map projection distortion. G´ is the coresponding projected point of the G accord- ing to [10]. Correction of horizontal angle Since the normal-section-to-geodesic correction distinguishes negligible values, it should be applied to skew-to-normal and arc-to-chord corrections. Skew-to-normal corrections express the angle between the directions of the spatial straight line and its corresponding normal section. Regardless to used coordinate system it is given by [4] ζ = hG 2ρm e2 sin(2α)cos2ϕG, (13) where ϕG and hG are the geodetic latitude and ellipsoidal height of the ground point re- spectively, ρm represents the curvature radius in the meridian plane and α is azimuth of observation vector. For practical calculations hG should be replace by hS + Zdg. Arc-to-chord correction is angle between tangent of projected geodesic and its corresponding chord line. Computation is individual for each map projection. Experiment The asset of the above-mentioned procedure was proven in [10]. Since the real LiDAR data is distorted by variety of random and systematic errors, the simulated data was used for the experiment. It was an applied method based on (2). It compared coordinates of ground Geoinformatics FCE CTU 15(1), 2016 40 Z. Švec: Direct georeferencing with correction of map projection hec α α O FP G Figure 6: Detail for deriving the Earth curvature correction. O stands for centre of the reference sphere (rotated 90 degrees clockwise). points restituted in e-frame and transformed to p-frame (assume as the “correct” data) and the coordinates of ground points restituted in p-frame with application of correction. Most differences gain sub-millimeter values (Tab 2.) Table 2: Experiment published in [10]. Error statistic of simulated data (mdatum= 1,00005; Krassovsky ellipsoid, UTM projection) Relative flight correction mean error [mm] σ [mm] max error [mm] height plane height plane height plane height 500m no 152.9 16.1 56.9 5.6 258.6 25 yes 0.2 0 0.1 0 0.3 0 2000m no 611.7 -41.9 227.7 89 1043.2 -254.9 yes 0.6 0.3 0.3 0.2 1.1 -0.4 8000m no 2446.8 -1871.1 914.9 1224.2 4313.5 -5278.4 yes 2.7 0 1.2 3.6 5.2 -7.2 Situation in the Czech Republic The Earth curvature correction is area-independent thus we will discuss the other components in the most frequent national coordinates in the area of the Czech Republic: S-JTSK (Datum of Uniform Trigonometric Cadastral Network) / Bpv (Baltic Vertical Datum - After Adjust- ment) and UTM 33(34)N / ellipsoidal height (GRS80). Regardless to a projection used, we can compute the skew-to-normal correction as [3] ζ́́ = 0, 108 hG 1000 sin (2α) cos2ϕG· (14) As the correction attains a value of 0,07´´ in extreme cases (Sněžka Mountain, azimuth 45°), we can disrespect it. Geoinformatics FCE CTU 15(1), 2016 41 Z. Švec: Direct georeferencing with correction of map projection S-JTSK/Bpv Change of the datum scale between WGS 84 to S-JTSK is -8.750 ppm, i.e if the observation vector is 1 km length it cause its shortening by 8.75 mm. This factor should be taken into account esspetialy in the case of high survey flights. S-JTSK use double conformal conical Krovak projection (EPSG code 5514). Local scale is given by [2] m = αR cosU N cos ϕ , (15) where the first part marks constants α = 1.000597, N radius of the cross section and U spherical latitude. The projection has 2 undistorted parallels. The local scale causes distortion in the range from -10cm/km to 14cm/km.The difference between D and D´ in some selected locations shows Tab.3 Arc-to-chord correction is given (after removing high-order terms) as [8] δ = (D́G − D́S) [2KS RG RS + KG RS RG ], (16) Ki = sin S0 − sin Si 6 sin S0 ≈5.314510−9∆Ri, (17) ∆Ri = Ri −R0,R0 = 1298039m,Ri = √ X2i + Y 2 i ,D́i = atan Yi Xi , (18) where R means distance of the point from krovak´s projection origin, D́i angle beetween X-axis and the point form projection origin and Si cartographic parallel. The undistorted parallel has been chosen as S0 = 78◦30́ . For calculation parameters of the ground point, we can use the rough position obtained by the uncorrected vector Xpdg.The magnitude of the correction depends on the length of the horizontal distance, orientation of Xpdg and distance from undistorted cartographical parallels. It attains a size of few (in extreme case up to 10) arcseconds. UTM/elipsoidal height (GRS80) For practical applications, elipsoid GRS80 and elipsoid WGS 84 are identical, thus the mdatum is the same. The local scale of each UTM strip is given by this simplified formula: m = (1 + λ2 2 cos2ϕ)m0, (19) where λ is longitude (measured from central meridian of the strip), ϕ is latitude and m0 thescale factor of the central meridian. For the reduction of the scale distortion at the boundaries of the strips m0 = 0, 9996 has been chosen. It follows the distortion -40cm/km at the central meridian and +17cm/km on the boundary of the strip. The length distortion of both projections shows Fig. 7. After some simplification and removing high-order terms we can compute arc-to-cord correc- tion as δ = −Ydg (3XS + Xdg ) 6m20R2 · (20) Geoinformatics FCE CTU 15(1), 2016 42 Z. Švec: Direct georeferencing with correction of map projection Figure 7: Scale distortion of UTM 33N (up), S-JTSK (down). Unit of the values is cm/km. Geoinformatics FCE CTU 15(1), 2016 43 Z. Švec: Direct georeferencing with correction of map projection Table 3: Difference between D and D´ in some selected locations. Location H D S-JTSK UTM 33N [m] [m] m D’[m] ∆D [m] m D’[m] ∆D [m] Jizera 1166 1000 1.00006 999.877 -0.123 0.99962 999.437 -0.563 Zruč nad Sázavou 400 1000 0.9999 999.837 -0.163 0.9996 999.537 -0.463 Znojmo 300 1000 0.99994 999.893 -0.107 0.9997 999.653 -0.347 Conclusion Mathematically, the only rigid way of direct georeferencing is to restitute the scene in an e-frame. However, after applying the above mentioned corrections, the residuals of the p- frame deformations are negligible in comparison with the noise of GNSS/IMU measurements. Although some formulas are simplified, remaining residuals are below 1 cm even by the relative flying height 8000m. The least impacted have angular corrections, and it can be ignored in the case of low relative flying height (up to 1000 m). One of the main arguments for the choice of frames for georeferencing can be computational costs. The most demanding step comprises transformation of EOP to national coordinates (approximately four times higher computation time than the transformation of ground points). However, the GNSS/IMU record data with a frequency up to several hundred Hz and the state- of-the-art LiDAR systems emit pulses with frequency even 500 kHz [10]. Therefore the number of ground points is much higher than the number of EOP observation and transformation of ground points cause consequently higher computation effort. Georeferencing in national coordinates is strictly recommended to verify whether software take into account the p-frame corrections. References [1] Antonio Angrisano. “GNSS/INS Integration Methods”. PhD thesis. Universita degli studi di Napoli, 2010. url: http : / / www . ucalgary . ca / engo _ webdocs / other / Angrisano_PhD%20Thesis_ENG2.pdf. [2] Petr Buchar. Základy matematické kartografie. (in Czech). ČVUT v Praze, 2002. [3] Miloš Cimbálník and Leoš Mervart. Vyšší Geodézie 1. (in Czech). ČVUT v Praze, 1996. [4] Rodney E. Deakin. Traverse Computation on the Ellipsoid and on the UniversalTrans- verse Mercator Projection. RMIT University, Mar. 2010. url: http://www.mygeodesy. id.au/documents/Trav_Comp_V2.1.pdf. [5] Dorota Grejner-Brzezinska, Charles Toth, and Yudan Yi. “On improving navigation accuracy of GPS/INS systems”. In: Photogrammetric engineering & remote sensing 71.4 (2005), pp. 377–389. url: http://eserv.asprs.org/PERS/2005journal/apr/ 2005_04_377-389.pdf. [6] B. Hofmann-Wellenhof, H. Lichtenegger, and J. Collins. Global Positioning System The- ory and Practice. 5th ed. Springer, 2001. doi: 10.1007/978- 3- 7091- 6199- 9. url: http://eserv.asprs.org/PERS/2005journal/apr/2005_04_377-389.pdf. Geoinformatics FCE CTU 15(1), 2016 44 http://www.ucalgary.ca/engo_webdocs/other/Angrisano_PhD%20Thesis_ENG2.pdf http://www.ucalgary.ca/engo_webdocs/other/Angrisano_PhD%20Thesis_ENG2.pdf http://www.mygeodesy.id.au/documents/Trav_Comp_V2.1.pdf http://www.mygeodesy.id.au/documents/Trav_Comp_V2.1.pdf http://eserv.asprs.org/PERS/2005journal/apr/2005_04_377-389.pdf http://eserv.asprs.org/PERS/2005journal/apr/2005_04_377-389.pdf http://dx.doi.org/10.1007/978-3-7091-6199-9 http://eserv.asprs.org/PERS/2005journal/apr/2005_04_377-389.pdf Z. Švec: Direct georeferencing with correction of map projection [7] Klaus Legat. “Approximate direct georeferencing in national coordinates”. In: ISPRS Journal of Photogrammetry and Remote Sensing 60 (2006), pp. 239–255. doi: 10.1016/ j.isprsjprs.2006.02.004. [8] Jan Skaloud. “Direct georeferencing in aerial photogrammetric mapping”. In: Pho- togrammetric Engineering and Remote Sensing 68.3 (2002), pp. 207–210. url: http: //worldcat.org/issn/00991112. [9] Jan Skaloud and Klaus Legat. “Theory and reality of direct georeferencing in national coordinates”. In: International Journal of Photogrammetry and Remote Sensing 63 (Mar. 2008), pp. 272–282. doi: 10.1016/j.isprsjprs.2007.09.002. [10] Yongjun Zhang and Xiang Shen. “Direct georeferencing of airborne LiDAR data in national coordinates”. In: ISPRS Journal of Photogrammetry and Remote Sensing 84 (Oct. 2013), pp. 43–51. doi: 10.1016/j.isprsjprs.2013.07.003. Geoinformatics FCE CTU 15(1), 2016 45 http://dx.doi.org/10.1016/j.isprsjprs.2006.02.004 http://dx.doi.org/10.1016/j.isprsjprs.2006.02.004 http://worldcat.org/issn/00991112 http://worldcat.org/issn/00991112 http://dx.doi.org/10.1016/j.isprsjprs.2007.09.002 http://dx.doi.org/10.1016/j.isprsjprs.2013.07.003