idnits 2.17.1 draft-thomson-geopriv-uncertainty-08.txt: Checking boilerplate required by RFC 5378 and the IETF Trust (see https://trustee.ietf.org/license-info): ---------------------------------------------------------------------------- No issues found here. Checking nits according to https://www.ietf.org/id-info/1id-guidelines.txt: ---------------------------------------------------------------------------- No issues found here. Checking nits according to https://www.ietf.org/id-info/checklist : ---------------------------------------------------------------------------- ** The document seems to lack an IANA Considerations section. (See Section 2.2 of https://www.ietf.org/id-info/checklist for how to handle the case when there are no actions for IANA.) Miscellaneous warnings: ---------------------------------------------------------------------------- == The copyright year in the IETF Trust and authors Copyright Line does not match the current year -- The document seems to lack a disclaimer for pre-RFC5378 work, but may have content which was first submitted before 10 November 2008. If you have contacted all the original authors and they are all willing to grant the BCP78 rights to the IETF Trust, then this is fine, and you can ignore this comment. If not, you may need to add the pre-RFC5378 disclaimer. (See the Legal Provisions document at https://trustee.ietf.org/license-info for more information.) -- The document date (August 27, 2013) is 3893 days in the past. Is this intentional? Checking references for intended status: Informational ---------------------------------------------------------------------------- -- Looks like a reference, but probably isn't: '1' on line 1203 -- Looks like a reference, but probably isn't: '2' on line 530 -- Looks like a reference, but probably isn't: '3' on line 530 == Missing Reference: '2d' is mentioned on line 671, but not defined == Missing Reference: '3d' is mentioned on line 671, but not defined -- Looks like a reference, but probably isn't: '0' on line 1203 == Unused Reference: 'RFC3694' is defined on line 1081, but no explicit reference was found in the text -- Obsolete informational reference (is this intentional?): RFC 3825 (Obsoleted by RFC 6225) Summary: 1 error (**), 0 flaws (~~), 4 warnings (==), 7 comments (--). Run idnits with the --verbose option for more detailed information about the items above. -------------------------------------------------------------------------------- 2 GEOPRIV M. Thomson 3 Internet-Draft Microsoft 4 Intended status: Informational J. Winterbottom 5 Expires: February 28, 2014 Andrew Corporation 6 August 27, 2013 8 Representation of Uncertainty and Confidence in PIDF-LO 9 draft-thomson-geopriv-uncertainty-08 11 Abstract 13 The key concepts of uncertainty and confidence as they pertain to 14 location information are defined. Methods for the manipulation of 15 location estimates that include uncertainty information are outlined. 17 Status of This Memo 19 This Internet-Draft is submitted in full conformance with the 20 provisions of BCP 78 and BCP 79. 22 Internet-Drafts are working documents of the Internet Engineering 23 Task Force (IETF). Note that other groups may also distribute 24 working documents as Internet-Drafts. The list of current Internet- 25 Drafts is at http://datatracker.ietf.org/drafts/current/. 27 Internet-Drafts are draft documents valid for a maximum of six months 28 and may be updated, replaced, or obsoleted by other documents at any 29 time. It is inappropriate to use Internet-Drafts as reference 30 material or to cite them other than as "work in progress." 32 This Internet-Draft will expire on February 28, 2014. 34 Copyright Notice 36 Copyright (c) 2013 IETF Trust and the persons identified as the 37 document authors. All rights reserved. 39 This document is subject to BCP 78 and the IETF Trust's Legal 40 Provisions Relating to IETF Documents 41 (http://trustee.ietf.org/license-info) in effect on the date of 42 publication of this document. Please review these documents 43 carefully, as they describe your rights and restrictions with respect 44 to this document. Code Components extracted from this document must 45 include Simplified BSD License text as described in Section 4.e of 46 the Trust Legal Provisions and are provided without warranty as 47 described in the Simplified BSD License. 49 Table of Contents 51 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . 2 52 1.1. Conventions and Terminology . . . . . . . . . . . . . . . 3 53 2. A General Definition of Uncertainty . . . . . . . . . . . . . 3 54 2.1. Uncertainty as a Probability Distribution . . . . . . . . 4 55 2.2. Deprecation of the Terms Precision and Resolution . . . . 6 56 2.3. Accuracy as a Qualitative Concept . . . . . . . . . . . . 7 57 3. Uncertainty in Location . . . . . . . . . . . . . . . . . . . 7 58 3.1. Representation of Uncertainty and Confidence in PIDF-LO . 7 59 3.2. Uncertainty and Confidence for Civic Addresses . . . . . 8 60 3.3. DHCP Location Configuration Information and Uncertainty . 9 61 4. Manipulation of Uncertainty . . . . . . . . . . . . . . . . . 9 62 4.1. Reduction of a Location Estimate to a Point . . . . . . . 9 63 4.1.1. Centroid Calculation . . . . . . . . . . . . . . . . 10 64 4.1.1.1. Arc-Band Centroid . . . . . . . . . . . . . . . . 11 65 4.1.1.2. Polygon Centroid . . . . . . . . . . . . . . . . 11 66 4.2. Conversion to Circle or Sphere . . . . . . . . . . . . . 14 67 4.3. Three-Dimensional to Two-Dimensional Conversion . . . . . 14 68 4.4. Increasing and Decreasing Uncertainty and Confidence . . 15 69 4.4.1. Rectangular Distributions . . . . . . . . . . . . . . 16 70 4.4.2. Normal Distributions . . . . . . . . . . . . . . . . 16 71 4.5. Determining Whether a Location is Within a Given Region . 17 72 4.5.1. Determining the Area of Overlap for Two Circles . . . 18 73 4.5.2. Determining the Area of Overlap for Two Polygons . . 19 74 5. Examples . . . . . . . . . . . . . . . . . . . . . . . . . . 19 75 5.1. Reduction to a Point or Circle . . . . . . . . . . . . . 19 76 5.2. Increasing and Decreasing Confidence . . . . . . . . . . 22 77 5.3. Matching Location Estimates to Regions of Interest . . . 22 78 6. Security Considerations . . . . . . . . . . . . . . . . . . . 22 79 7. Acknowledgements . . . . . . . . . . . . . . . . . . . . . . 23 80 8. Informative References . . . . . . . . . . . . . . . . . . . 23 81 Appendix A. Conversion Between Cartesian and Geodetic 82 Coordinates in WGS84 . . . . . . . . . . . . . . . . 24 83 Appendix B. Calculating the Upward Normal of a Polygon . . . . . 25 84 B.1. Checking that a Polygon Upward Normal Points Up . . . . . 26 85 Authors' Addresses . . . . . . . . . . . . . . . . . . . . . . . 27 87 1. Introduction 89 Location information represents an estimation of the position of a 90 Target. Under ideal circumstances, a location estimate precisely 91 reflects the actual location of the Target. In reality, there are 92 many factors that introduce errors into the measurements that are 93 used to determine location estimates. 95 The process by which measurements are combined to generate a location 96 estimate is outside of the scope of work within the IETF. However, 97 the results of such a process are carried in IETF data formats and 98 protocols. This document outlines how uncertainty, and its 99 associated datum, confidence, are expressed and interpreted. 101 This document provides a common nomenclature for discussing 102 uncertainty and confidence as they relate to location information. 104 This document also provides guidance on how to manage location 105 information that includes uncertainty. Methods for expanding or 106 reducing uncertainty to obtain a required level of confidence are 107 described. Methods for determining the probability that a Target is 108 within a specified region based on their location estimate are 109 described. These methods are simplified by making certain 110 assumptions about the location estimate and are designed to be 111 applicable to location estimates in a relatively small area. 113 1.1. Conventions and Terminology 115 This document assumes a basic understanding of the principles of 116 mathematics, particularly statistics and geometry. 118 Some terminology is borrowed from [RFC3693] and [RFC6280]. 120 Mathematical formulae are presented using the following notation: add 121 "+", subtract "-", multiply "*", divide "/", power "^" and absolute 122 value "|x|". Precedence is indicated using parentheses. 123 Mathematical functions are represented by common abbreviations: 124 square root "sqrt(x)", sine "sin(x)", cosine "cos(x)", inverse cosine 125 "acos(x)", tangent "tan(x)", inverse tangent "atan(x)", error 126 function "erf(x)", and inverse error function "erfinv(x)". 128 2. A General Definition of Uncertainty 130 Uncertainty results from the limitations of measurement. In 131 measuring any observable quantity, errors from a range of sources 132 affect the result. Uncertainty is a quantification of what is known 133 about the observed quantity, either through the limitations of 134 measurement or through inherent variability of the quantity. 136 Uncertainty is most completely described by a probability 137 distribution. A probability distribution assigns a probability to 138 possible values for the quantity. 140 A probability distribution describing a measured quantity can be 141 arbitrarily complex and so it is desirable to find a simplified 142 model. One approach commonly taken is to reduce the probability 143 distribution to a confidence interval. Many alternative models are 144 used in other areas, but study of those is not the focus of this 145 document. 147 In addition to the central estimate of the observed quantity, a 148 confidence interval is succintly described by two values: an error 149 range and a confidence. The error range describes an interval and 150 the confidence describes an estimated upper bound on the probability 151 that a "true" value is found within the extents defined by the error. 153 In the following example, a measurement result for a length is shown 154 as a nominal value with additional information on error range (0.0043 155 meters) and confidence (95%). 157 e.g. x = 1.00742 +/- 0.0043 meters at 95% confidence 159 This result indicates that the measurement indicates that the value 160 of "x" between 1.00312 and 1.01172 meters with 95% probability. No 161 other assertion is made: in particular, this does not assert that x 162 is 1.00742. 164 This document uses the term _uncertainty_ to refer in general to the 165 concept as well as more specifically to refer to the error increment. 167 Uncertainty and confidence for location estimates can be derived in a 168 number of ways. This document does not attempt to enumerate the many 169 methods for determining uncertainty. [ISO.GUM] and [NIST.TN1297] 170 provide a set of general guidelines for determining and manipulating 171 measurement uncertainty. This document applies that general guidance 172 for consumers of location information. 174 2.1. Uncertainty as a Probability Distribution 176 The Probability Density Function (PDF) that is described by 177 uncertainty indicates the probability that the "true" value lies at 178 any one point. The shape of the probability distribution can vary 179 depending on the method that is used to determine the result. The 180 two probability density functions most generally applicable most 181 applicable to location information are considered in this document: 183 o The normal PDF (also referred to as a Gaussian PDF) is used where 184 a large number of small random factors contribute to errors. The 185 value used for the error range in a normal PDF is related to the 186 standard deviation of the distribution. 188 o A rectangular PDF is used where the errors are known to be 189 consistent across a limited range. A rectangular PDF can occur 190 where a single error source, such as a rounding error, is 191 significantly larger than other errors. A rectangular PDF is 192 often described by the half-width of the distribution; that is, 193 half the width of the distribution. 195 Each of these probability density functions can be characterized by 196 its center point, or mean, and its width. For a normal distribution, 197 uncertainty and confidence together are related to the standard 198 deviation (see Section 4.4). For a rectangular distribution, half of 199 the width of the distribution is used. 201 Figure 1 shows a normal and rectangular probability density function 202 with the mean (m) and standard deviation (s) labelled. The half- 203 width (h) of the rectangular distribution is also indicated. 205 ***** *** Normal PDF 206 ** : ** --- Rectangular PDF 207 ** : ** 208 ** : ** 209 .---------*---------------*---------. 210 | ** : ** | 211 | ** : ** | 212 | * <-- s -->: * | 213 | * : : : * | 214 | ** : ** | 215 | * : : : * | 216 | * : * | 217 |** : : : **| 218 ** : ** 219 *** | : : : | *** 220 ***** | :<------ h ------>| ***** 221 .****-------+.......:.........:.........:.......+-------*****. 222 m 224 Figure 1: Normal and Rectangular Probability Density Functions 226 For a given PDF, the value of the PDF describes the probability that 227 the "true" value is found at that point. Confidence for any given 228 interval is the total probability of the "true" value being in that 229 range, defined as the integral of the PDF over the interval. 231 The probability of the "true" value falling between two points is 232 found by finding the area under the curve between the points (that 233 is, the integral of the curve between the points). For any given 234 PDF, the area under the curve for the entire range from negative 235 infinity to positive infinity is 1 or (100%). Therefore, the 236 confidence over any interval of uncertainty is always less than 237 100%. 239 Figure 2 shows how confidence is determined for a normal 240 distribution. The area of the shaded region gives the confidence (c) 241 for the interval between "m-u" and "m+u". 243 ***** 244 **:::::** 245 **:::::::::** 246 **:::::::::::** 247 *:::::::::::::::* 248 **:::::::::::::::** 249 **:::::::::::::::::** 250 *:::::::::::::::::::::* 251 *:::::::::::::::::::::::* 252 **:::::::::::::::::::::::** 253 *:::::::::::: c ::::::::::::* 254 *:::::::::::::::::::::::::::::* 255 **|:::::::::::::::::::::::::::::|** 256 ** |:::::::::::::::::::::::::::::| ** 257 *** |:::::::::::::::::::::::::::::| *** 258 ***** |:::::::::::::::::::::::::::::| ***** 259 .****..........!:::::::::::::::::::::::::::::!..........*****. 260 | | | 261 (m-u) m (m+u) 263 Figure 2: Confidence as the Integral of a PDF 265 In Section 4.4, methods are described for manipulating uncertainty if 266 the shape of the PDF is known. 268 2.2. Deprecation of the Terms Precision and Resolution 270 The terms _Precision_ and _Resolution_ are defined in RFC 3693 271 [RFC3693]. These definitions were intended to provide a common 272 nomenclature for discussing uncertainty; however, these particular 273 terms have many different uses in other fields and their definitions 274 are not sufficient to avoid confusion about their meaning. These 275 terms are unsuitable for use in relation to quantitative concepts 276 when discussing uncertainty and confidence in relation to location 277 information. 279 2.3. Accuracy as a Qualitative Concept 281 Uncertainty is a quantitative concept. The term _accuracy_ is useful 282 in describing, qualitatively, the general concepts of location 283 information. Accuracy is generally useful when describing 284 qualitative aspects of location estimates. Accuracy is not a 285 suitable term for use in a quantitative context. 287 For instance, it could be appropriate to say that a location estimate 288 with uncertainty "X" is more accurate than a location estimate with 289 uncertainty "2X" at the same confidence. It is not appropriate to 290 assign a number to "accuracy", nor is it appropriate to refer to any 291 component of uncertainty or confidence as "accuracy". That is, to 292 say that the "accuracy" for the first location estimate is "X" would 293 be an erroneous use of this term. 295 3. Uncertainty in Location 297 A _location estimate_ is the result of location determination. A 298 location estimate is subject to uncertainty like any other 299 observation. However, unlike a simple measure of a one dimensional 300 property like length, a location estimate is specified in two or 301 three dimensions. 303 Uncertainty in 2- or 3-dimensional locations can be described using 304 confidence intervals. The confidence interval for a location 305 estimate in two or three dimensional space is expressed as a subset 306 of that space. This document uses the term _region of uncertainty_ 307 to refer to the area or volume that describes the confidence 308 interval. 310 Areas or volumes that describe regions of uncertainty can be formed 311 by the combination of two or three one-dimensional ranges, or more 312 complex shapes could be described. 314 3.1. Representation of Uncertainty and Confidence in PIDF-LO 316 A set of shapes suitable for the expression of uncertainty in 317 location estimates in the Presence Information Data Format - Location 318 Object (PIDF-LO) are described in [GeoShape]. These shapes are the 319 recommended form for the representation of uncertainty in PIDF-LO 320 [RFC4119] documents. 322 The PIDF-LO does not include an indication of confidence, but that 323 confidence is 95%, by definition in [RFC5491]. Similarly, the PIDF- 324 LO format does not provide an indication of the shape of the PDF. 326 Absence of uncertainty information in a PIDF-LO document does not 327 indicate that there is no uncertainty in the location estimate. 328 Uncertainty might not have been calculated for the estimate, or it 329 may be withheld for privacy purposes. 331 If the Point shape is used, confidence and uncertainty are unknown; a 332 receiver can either assume a confidence of 0% or infinite 333 uncertainty. The same principle applies on the altitude axis for 334 two-dimension shapes like the Circle. 336 3.2. Uncertainty and Confidence for Civic Addresses 338 Civic addresses [RFC5139] inherently include uncertainty, based on 339 the area of the most precise element that is specified. Uncertainty 340 is effectively defined by the presence or absence of elements -- 341 elements that are not present are deemed to be uncertain. 343 To apply the concept of uncertainty to civic addresses, it is helpful 344 to unify the conceptual models of civic address with geodetic 345 location information. 347 Note: This view is one perspective on the process of geo-coding - 348 the translation of a civic address to a geodetic location. 350 In the unified view, a civic address defines a series of (sometimes 351 non-orthogonal) spatial partitions. The first is the implicit 352 partition that identifies the surface of the earth and the space near 353 the surface. The second is the country. Each label that is included 354 in a civic address provides information about a different set of 355 spatial partitions. Some partions require slight adjustments from a 356 standard interpretation: for instance, a road includes all properties 357 that adjoin the street. Each label might need to be interpreted with 358 other values to provide context. 360 As a value at each level is interpreted, one or more spatial 361 partitions at that level are selected, and all other partitions of 362 that type are excluded. For non-orthogonal partitions, only the 363 portion of the partition that fits within the existing space is 364 selected. This is what distinguishes King Street in Sydney from King 365 Street in Melbourne. Each defined element selects a partition of 366 space. The resulting location is the intersection of all selected 367 spaces. 369 The resulting spatial partition can be considered to represent a 370 region of uncertainty. At no stage does this process select a point; 371 although, as spaces get smaller this distinction might have no 372 practical significance and an approximation if a point could be used. 374 Uncertainty in civic addresses can be increased by removing elements. 375 This doesn't necessarily improve confidence in the same way that 376 arbitrarily increasing uncertainty in a geodetic location doesn't 377 increase confidence. 379 3.3. DHCP Location Configuration Information and Uncertainty 381 Location information is often measured in two or three dimensions; 382 expressions of uncertainty in one dimension only are rare. The 383 "resolution" parameters in [RFC3825] provide an indication of 384 uncertainty in one dimension. 386 [RFC3825] defines a means for representing uncertainty, but a value 387 for confidence is not specified. A default value of 95% confidence 388 can be assumed for the combination of the uncertainty on each axis. 389 That is, the confidence of the resultant rectangular polygon or prism 390 is 95%. 392 4. Manipulation of Uncertainty 394 This section deals with manipulation of location information that 395 contains uncertainty. 397 The following rules generally apply when manipulating location 398 information: 400 o Where calculations are performed on coordinate information, these 401 should be performed in Cartesian space and the results converted 402 back to latitude, longitude and altitude. A method for converting 403 to and from Cartesian coordinates is included in Appendix A. 405 While some approximation methods are useful in simplifying 406 calculations, treating latitude and longitude as Cartesian axes 407 is never advisable. The two axes are not orthogonal. Errors 408 can arise from the curvature of the earth and from the 409 convergence of longitude lines. 411 o Normal rounding rules do not apply when rounding uncertainty. 412 When rounding, the region of uncertainty always increases (that 413 is, errors are rounded up) and confidence is always rounded down 414 (see [NIST.TN1297]). This means that any manipulation of 415 uncertainty is a non-reversible operation; each manipulation can 416 result in the loss of some information. 418 4.1. Reduction of a Location Estimate to a Point 419 Manipulating location estimates that include uncertainty information 420 requires additional complexity in systems. In some cases, systems 421 only operate on definitive values, that is, a single point. 423 This section describes algorithms for reducing location estimates to 424 a simple form without uncertainty information. Having a consistent 425 means for reducing location estimates allows for interaction between 426 applications that are able to use uncertainty information and those 427 that cannot. 429 Note: Reduction of a location estimate to a point constitutes a 430 reduction in information. Removing uncertainty information can 431 degrade results in some applications. Also, there is a natural 432 tendency to misinterpret a point location as representing a 433 location without uncertainty. This could lead to more serious 434 errors. Therefore, these algorithms should only be applied where 435 necessary. 437 Several different approaches can be taken when reducing a location 438 estimate to a point. Different methods each make a set of 439 assumptions about the properties of the PDF and the selected point; 440 no one method is more "correct" than any other. For any given region 441 of uncertainty, selecting an arbitrary point within the area could be 442 considered valid; however, given the aforementioned problems with 443 point locations, a more rigorous approach is appropriate. 445 Given a result with a known distribution, selecting the point within 446 the area that has the highest probability is a more rigorous method. 447 Alternatively, a point could be selected that minimizes the overall 448 error; that is, it minimises the expected value of the difference 449 between the selected point and the "true" value. 451 If a rectangular distribution is assumed, the centroid of the area or 452 volume minimizes the overall error. Minimizing the error for a 453 normal distribution is mathematically complex. Therefore, this 454 document opts to select the centroid of the region of uncertainty 455 when selecting a point. 457 4.1.1. Centroid Calculation 459 For regular shapes, such as Circle, Sphere, Ellipse and Ellipsoid, 460 this approach equates to the center point of the region. For regions 461 of uncertainty that are expressed as regular Polygons and Prisms the 462 center point is also the most appropriate selection. 464 For the Arc-Band shape and non-regular Polygons and Prisms, selecting 465 the centroid of the area or volume minimizes the overall error. This 466 assumes that the PDF is rectangular. 468 Note: The centroid of a concave Polygon or Arc-Band shape is not 469 necessarily within the region of uncertainty. 471 4.1.1.1. Arc-Band Centroid 473 The centroid of the Arc-Band shape is found along a line that bisects 474 the arc. The centroid can be found at the following distance from 475 the starting point of the arc-band (assuming an arc-band with an 476 inner radius of "r", outer radius "R", start angle "a", and opening 477 angle "o"): 479 d = 4 * sin(o/2) * (R*R + R*r + r*r) / (3*o*(R + r)) 481 This point can be found along the line that bisects the arc; that is, 482 the line at an angle of "a + (o/2)". Negative values are possible if 483 the angle of opening is greater than 180 degrees; negative values 484 indicate that the centroid is found along the angle "a + (o/ 485 2) + 180". 487 4.1.1.2. Polygon Centroid 489 Calculating a centroid for the Polygon and Prism shapes is more 490 complex. Polygons that are specified using geodetic coordinates are 491 not necessarily coplanar. For Polygons that are specified without an 492 altitude, choose a value for altitude before attempting this process; 493 an altitude of 0 is acceptable. 495 The method described in this section is simplified by assuming 496 that the surface of the earth is locally flat. This method 497 degrades as polygons become larger; see [GeoShape] for 498 recommendations on polygon size. 500 The polygon is translated to a new coordinate system that has an x-y 501 plane roughly parallel to the polygon. This enables the elimination 502 of z-axis values and calculating a centroid can be done using only x 503 and y coordinates. This requires that the upward normal for the 504 polygon is known. 506 To translate the polygon coordinates, apply the process described in 507 Appendix B to find the normal vector "N = [Nx,Ny,Nz]". This value 508 should be made a unit vector to ensure that the transformation matrix 509 is a special orthogonal matrix. From this vector, select two vectors 510 that are perpendicular to this vector and combine these into a 511 transformation matrix. 513 If "Nx" and "Ny" are non-zero, the matrices in Figure 3 can be used, 514 given "p = sqrt(Nx^2 + Ny^2)". More transformations are provided 515 later in this section for cases where "Nx" or "Ny" are zero. 517 [ -Ny/p Nx/p 0 ] [ -Ny/p -Nx*Nz/p Nx ] 518 T = [ -Nx*Nz/p -Ny*Nz/p p ] T' = [ Nx/p -Ny*Nz/p Ny ] 519 [ Nx Ny Nz ] [ 0 p Nz ] 520 (Transform) (Reverse Transform) 522 Figure 3: Recommended Transformation Matrices 524 To apply a transform to each point in the polygon, form a matrix from 525 the ECEF coordinates and use matrix multiplication to determine the 526 translated coordinates. 528 [ -Ny/p Nx/p 0 ] [ x[1] x[2] x[3] ... x[n] ] 529 [ -Nx*Nz/p -Ny*Nz/p p ] * [ y[1] y[2] y[3] ... y[n] ] 530 [ Nx Ny Nz ] [ z[1] z[2] z[3] ... z[n] ] 532 [ x'[1] x'[2] x'[3] ... x'[n] ] 533 = [ y'[1] y'[2] y'[3] ... y'[n] ] 534 [ z'[1] z'[2] z'[3] ... z'[n] ] 536 Figure 4: Transformation 538 Alternatively, direct multiplication can be used to achieve the same 539 result: 541 x'[i] = -Ny * x[i] / p + Nx * y[i] / p 543 y'[i] = -Nx * Nz * x[i] / p - Ny * Nz * y[i] / p + p * z[i] 545 z'[i] = Nx * x[i] + Ny * y[i] + Nz * z[i] 547 The first and second rows of this matrix ("x'" and "y'") contain the 548 values that are used to calculate the centroid of the polygon. To 549 find the centroid of this polygon, first find the area using: 551 A = sum from i=1..n of (x'[i]*y'[i+1]-x'[i+1]*y'[i]) / 2 553 For these formulae, treat each set of coordinates as circular, that 554 is "x'[0] == x'[n]" and "x'[n+1] == x'[1]". Based on the area, the 555 centroid along each axis can be determined by: 557 Cx' = sum (x'[i]+x'[i+1]) * (x'[i]*y'[i+1]-x'[i+1]*y'[i]) / (6*A) 559 Cy' = sum (y'[i]+y'[i+1]) * (x'[i]*y'[i+1]-x'[i+1]*y'[i]) / (6*A) 561 Note: The formula for the area of a polygon will return a negative 562 value if the polygon is specified in clockwise direction. This 563 can be used to determine the orientation of the polygon. 565 The third row contains a distance from a plane parallel to the 566 polygon. If the polygon is coplanar, then the values for "z'" are 567 identical; however, the constraints recommended in [RFC5491] mean 568 that this is rarely the case. To determine "Cz'", average these 569 values: 571 Cz' = sum z'[i] / n 573 Once the centroid is known in the transformed coordinates, these can 574 be transformed back to the original coordinate system. The reverse 575 transformation is shown in Figure 5. 577 [ -Ny/p -Nx*Nz/p Nx ] [ Cx' ] [ Cx ] 578 [ Nx/p -Ny*Nz/p Ny ] * [ Cy' ] = [ Cy ] 579 [ 0 p Nz ] [ sum of z'[i] / n ] [ Cz ] 581 Figure 5: Reverse Transformation 583 The reverse transformation can be applied directly as follows: 585 Cx = -Ny * Cx' / p - Nx * Nz * Cy' / p + Nx * Cz' 587 Cy = Nx * Cx' / p - Ny * Nz * Cy' / p + Ny * Cz' 589 Cz = p * Cy' + Nz * Cz' 591 The ECEF value "[Cx,Cy,Cz]" can then be converted back to geodetic 592 coordinates. Given a polygon that is defined with no altitude or 593 equal altitudes for each point, the altitude of the result can either 594 be ignored or reset after converting back to a geodetic value. 596 The centroid of the Prism shape is found by finding the centroid of 597 the base polygon and raising the point by half the height of the 598 prism. This can be added to altitude of the final result; 599 alternatively, this can be added to "Cz'", which ensures that 600 negative height is correctly applied to polygons that are defined in 601 a "clockwise" direction. 603 The recommended transforms only apply if "Nx" and "Ny" are non-zero. 604 If the normal vector is "[0,0,1]" (that is, along the z-axis), then 605 no transform is necessary. Similarly, if the normal vector is 606 "[0,1,0]" or "[1,0,0]", avoid the transformation and use the x and z 607 coordinates or y and z coordinates (respectively) in the centroid 608 calculation phase. If either "Nx" or "Ny" are zero, the alternative 609 transform matrices in Figure 6 can be used. The reverse transform is 610 the transpose of this matrix. 612 if Nx == 0: | if Ny == 0: 613 [ 0 -Nz Ny ] [ 0 1 0 ] | [ -Nz 0 Nx ] 614 T = [ 1 0 0 ] T' = [ -Nz 0 Ny ] | T = T' = [ 0 1 0 ] 615 [ 0 Ny Nz ] [ Ny 0 Nz ] | [ Nx 0 Nz ] 617 Figure 6: Alternative Transformation Matrices 619 4.2. Conversion to Circle or Sphere 621 The Circle or Sphere are simple shapes that suit a range of 622 applications. A circle or sphere contains fewer units of data to 623 manipulate, which simplifies operations on location estimates. 625 The simplest method for converting a location estimate to a Circle or 626 Sphere shape is to determine the centroid and then find the longest 627 distance to any point in the region of uncertainty to that point. 628 This distance can be determined based on the shape type: 630 Circle/Sphere: No conversion necessary. 632 Ellipse/Ellipsoid: The greater of either semi-major axis or altitude 633 uncertainty. 635 Polygon/Prism: The distance to the furthest vertex of the polygon 636 (for a Prism, it is only necessary to check points on the base). 638 Arc-Band: The furthest length from the centroid to the points where 639 the inner and outer arc end. This distance can be calculated by 640 finding the larger of the two following formulae: 642 X = sqrt( d*d + R*R - 2*d*R*cos(o/2) ) 644 x = sqrt( d*d + r*r - 2*d*r*cos(o/2) ) 646 Once the Circle or Sphere shape is found, the associated confidence 647 can be increased if the result is known to follow a normal 648 distribution. However, this is a complicated process and provides 649 limited benefit. In many cases it also violates the constraint that 650 confidence in each dimension be the same. Confidence should be 651 unchanged when performing this conversion. 653 Two dimensional shapes are converted to a Circle; three dimensional 654 shapes are converted to a Sphere. 656 4.3. Three-Dimensional to Two-Dimensional Conversion 658 A three-dimensional shape can be easily converted to a two- 659 dimensional shape by removing the altitude component. A sphere 660 becomes a circle; a prism becomes a polygon; an ellipsoid becomes an 661 ellipse. Each conversion is simple, requiring only the removal of 662 those elements relating to altitude. 664 The altitude is unspecified for a two-dimensional shape and therefore 665 has unlimited uncertainty along the vertical axis. The confidence 666 for the two-dimensional shape is thus higher than the three- 667 dimensional shape. Assuming equal confidence on each axis, the 668 confidence of the circle can be increased using the following 669 approximate formula: 671 C[2d] >= C[3d] ^ (2/3) 673 "C[2d]" is the confidence of the two-dimensional shape and "C[3d]" is 674 the confidence of the three-dimensional shape. For example, a Sphere 675 with a confidence of 95% can be simplified to a Circle of equal 676 radius with confidence of 96.6%. 678 4.4. Increasing and Decreasing Uncertainty and Confidence 680 The combination of uncertainty and confidence provide a great deal of 681 information about the nature of the data that is being measured. If 682 both uncertainty, confidence and PDF are known, certain information 683 can be extrapolated. In particular, the uncertainty can be scaled to 684 meet a desired confidence or the confidence for a particular region 685 of uncertainty can be found. 687 In general, confidence decreases as the region of uncertainty 688 decreases in size and confidence increases as the region of 689 uncertainty increases in size. However, this depends on the PDF; 690 expanding the region of uncertainty for a rectangular distribution 691 has no effect on confidence without additional information. If the 692 region of uncertainty is increased during the process of obfuscation 693 (see [I-D.thomson-geopriv-location-obscuring]), then the confidence 694 cannot be increased. 696 A region of uncertainty that is reduced in size always has a lower 697 confidence. 699 A region of uncertainty that has an unknown PDF shape cannot be 700 reduced in size reliably. The region of uncertainty can be expanded, 701 but only if confidence is not increased. 703 This section makes the simplifying assumption that location 704 information is symmetrically and evenly distributed in each 705 dimension. This is not necessarily true in practice. If better 706 information is available, alternative methods might produce better 707 results. 709 4.4.1. Rectangular Distributions 711 Uncertainty that follows a rectangular distribution can only be 712 decreased in size. Since the PDF is constant over the region of 713 uncertainty, the resulting confidence is determined by the following 714 formula: 716 Cr = Co * Ur / Uo 718 Where "Uo" and "Ur" are the sizes of the original and reduced regions 719 of uncertainty (either the area or the volume of the region); "Co" 720 and "Cb" are the confidence values associated with each region. 722 Information is lost by decreasing the region of uncertainty for a 723 rectangular distribution. Once reduced in size, the uncertainty 724 region cannot subsequently be increased in size. 726 4.4.2. Normal Distributions 728 Uncertainty and confidence can be both increased and decreased for a 729 normal distribution. However, the process is more complicated. 731 For a normal distribution, uncertainty and confidence are related to 732 the standard deviation of the function. The following function 733 defines the relationship between standard deviation, uncertainty and 734 confidence along a single axis: 736 S[x] = U[x] / ( sqrt(2) * erfinv(C[x]) ) 738 Where "S[x]" is the standard deviation, "U[x]" is the uncertainty and 739 "C[x]" is the confidence along a single axis. "erfinv" is the 740 inverse error function. 742 Scaling a normal distribution in two dimensions requires several 743 assumptions. Firstly, it is assumed that the distribution along each 744 axis is independent. Secondly, the confidence for each axis is the 745 same. Therefore, the confidence along each axis can be assumed to 746 be: 748 C[x] = Co ^ (1/n) 750 Where "C[x]" is the confidence along a single axis and "Co" is the 751 overall confidence and "n" is the number of dimensions in the 752 uncertainty. 754 Therefore, to find the uncertainty for each axis at a desired 755 confidence, "Cd", apply the following formula: 757 Ud[x] <= U[x] * (erfinv(Cd ^ (1/n)) / erfinv(Co ^ (1/n))) 759 For regular shapes, this formula can be applied as a scaling factor 760 in each dimension to reach a required confidence. 762 4.5. Determining Whether a Location is Within a Given Region 764 A number of applications require that a judgement be made about 765 whether a Target is within a given region of interest. Given a 766 location estimate with uncertainty, this judgement can be difficult. 767 A location estimate represents a probability distribution, and the 768 true location of the Target cannot be definitively known. Therefore, 769 the judgement relies on determining the probability that the Target 770 is within the region. 772 The probability that the Target is within a particular region is 773 found by integrating the PDF over the region. For a normal 774 distribution, there are no analytical methods that can be used to 775 determine the integral of the two or three dimensional PDF over an 776 arbitrary region. The complexity of numerical methods is also too 777 great to be useful in many applications; for example, finding the 778 integral of the PDF in two or three dimensions across the overlap 779 between the uncertainty region and the target region. If the PDF is 780 unknown, no determination can be made. When judging whether a 781 location is within a given region, uncertainties using these PDFs can 782 be assumed to be rectangular. If this assumption is made, the 783 confidence should be scaled to 95%, if possible. 785 Note: The selection of confidence has a significant impact on the 786 final result. Only use a different confidence if an uncertainty 787 value for 95% confidence cannot be found. 789 Given the assumption of a rectangular distribution, the probability 790 that a Target is found within a given region is found by first 791 finding the area (or volume) of overlap between the uncertainty 792 region and the region of interest. This is multiplied by the 793 confidence of the location estimate to determine the probability. 794 Figure 7 shows an example of finding the area of overlap between the 795 region of uncertainty and the region of interest. 797 _.-""""-._ 798 .' `. _ Region of 799 / \ / Uncertainty 800 ..+-"""--.. | 801 .-' | :::::: `-. | 802 ,' | :: Ao ::: `. | 803 / \ :::::::::: \ / 804 / `._ :::::: _.X 806 | `-....-' | 807 | | 808 | | 809 \ / 810 `. .' \_ Region of 811 `._ _.' Interest 812 `--..___..--' 814 Figure 7: Area of Overlap Between Two Circular Regions 816 Once the area of overlap, "Ao", is known, the probability that the 817 Target is within the region of interest, "Pi", is: 819 Pi = Co * Ao / Au 821 Given that the area of the region of uncertainty is "Au" and the 822 confidence is "Co". 824 This probability is often input to a decision process that has a 825 limited set of outcomes; therefore, a threshold value needs to be 826 selected. Depending on the application, different threshold 827 probabilities might be selected. In the absence of specific 828 recommendations, this document suggests that the probability be 829 greater than 50% before a decision is made. If the decision process 830 selects between two or more regions, as is required by [RFC5222], 831 then the region with the highest probability can be selected. 833 4.5.1. Determining the Area of Overlap for Two Circles 835 Determining the area of overlap between two arbitrary shapes is a 836 non-trivial process. Reducing areas to circles (see Section 4.2) 837 enables the application of the following process. 839 Given the radius of the first circle "r", the radius of the second 840 circle "R" and the distance between their center points "d", the 841 following set of formulas provide the area of overlap "Ao". 843 o If the circles don't overlap, that is "d >= r+R", "Ao" is zero. 845 o If one of the two circles is entirely within the other, that is 846 "d <= |r-R|", the area of overlap is the area of the smaller 847 circle. 849 o Otherwise, if the circles partially overlap, that is "d < r+R" and 850 "d > |r-R|", find "Ao" using: 852 a = (r^2 - R^2 + d^2)/(2*d) 853 Ao = r^2*acos(a/r) + R^2*acos((d - a)/R) - d*sqrt(r^2 - a^2) 855 A value for "d" can be determined by converting the center points to 856 Cartesian coordinates and calculating the distance between the two 857 center points: 859 d = sqrt((x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2) 861 4.5.2. Determining the Area of Overlap for Two Polygons 863 A calculation of overlap based on polygons can give better results 864 than the circle-based method. However, efficient calculation of 865 overlapping area is non-trivial. Algorithms such as Vatti's clipping 866 algorithm [Vatti92] can be used. 868 For large polygonal areas, it might be that geodesic interpolation is 869 used. In these cases, altitude is also frequently omitted in 870 describing the polygon. For such shapes, a planar projection can 871 still give a good approximation of the area of overlap if the larger 872 area polygon is projected onto the local tangent plane of the 873 smaller. This is only possible if the only area of interest is that 874 contained within the smaller polygon. Where the entire area of the 875 larger polygon is of interest, geodesic interpolation is necessary. 877 5. Examples 879 This section presents some examples of how to apply the methods 880 described in Section 4. 882 5.1. Reduction to a Point or Circle 884 Alice receives a location estimate from her LIS that contains a 885 ellipsoidal region of uncertainty. This information is provided at 886 19% confidence with a normal PDF. A PIDF-LO extract for this 887 information is shown in Figure 8. 889 890 891 892 -34.407242 150.882518 34 893 894 7.7156 895 896 897 3.31 898 899 900 28.7 902 903 904 43 905 906 907 908 909 911 Figure 8 913 This information can be reduced to a point simply by extracting the 914 center point, that is [-34.407242, 150.882518, 34]. 916 If some limited uncertainty were required, the estimate could be 917 converted into a circle or sphere. To convert to a sphere, the 918 radius is the largest of the semi-major, semi-minor and vertical 919 axes; in this case, 28.7 meters. 921 However, if only a circle is required, the altitude can be dropped as 922 can the altitude uncertainty (the vertical axis of the ellipsoid), 923 resulting in a circle at [-34.407242, 150.882518] of radius 7.7156 924 meters. 926 Bob receives a location estimate with a Polygon shape. This 927 information is shown in Figure 9. 929 930 931 932 933 -33.856625 151.215906 -33.856299 151.215343 934 -33.856326 151.214731 -33.857533 151.214495 935 -33.857720 151.214613 -33.857369 151.215375 936 -33.856625 151.215906 937 938 939 940 942 Figure 9 944 To convert this to a polygon, each point is firstly assigned an 945 altitude of zero and converted to ECEF coordinates (see Appendix A). 946 Then a normal vector for this polygon is found (see Appendix B). The 947 results of each of these stages is shown in Figure 10. Note that the 948 numbers shown are all rounded; no rounding is possible during this 949 process since rounding would contribute significant errors. 951 Polygon in ECEF coordinate space 952 (repeated point omitted and transposed to fit): 953 [ -4.6470e+06 2.5530e+06 -3.5333e+06 ] 954 [ -4.6470e+06 2.5531e+06 -3.5332e+06 ] 955 pecef = [ -4.6470e+06 2.5531e+06 -3.5332e+06 ] 956 [ -4.6469e+06 2.5531e+06 -3.5333e+06 ] 957 [ -4.6469e+06 2.5531e+06 -3.5334e+06 ] 958 [ -4.6469e+06 2.5531e+06 -3.5333e+06 ] 960 Normal Vector: n = [ -0.72782 0.39987 -0.55712 ] 962 Transformation Matrix: 963 [ -0.48152 -0.87643 0.00000 ] 964 t = [ -0.48828 0.26827 0.83043 ] 965 [ -0.72782 0.39987 -0.55712 ] 967 Transformed Coordinates: 968 [ 8.3206e+01 1.9809e+04 6.3715e+06 ] 969 [ 3.1107e+01 1.9845e+04 6.3715e+06 ] 970 pecef' = [ -2.5528e+01 1.9842e+04 6.3715e+06 ] 971 [ -4.7367e+01 1.9708e+04 6.3715e+06 ] 972 [ -3.6447e+01 1.9687e+04 6.3715e+06 ] 973 [ 3.4068e+01 1.9726e+04 6.3715e+06 ] 975 Two dimensional polygon area: A = 12600 m^2 976 Two-dimensional polygon centroid: C' = [ 8.8184e+00 1.9775e+04 ] 978 Average of pecef' z coordinates: 6.3715e+06 980 Reverse Transformation Matrix: 981 [ -0.48152 -0.48828 -0.72782 ] 982 t' = [ -0.87643 0.26827 0.39987 ] 983 [ 0.00000 0.83043 -0.55712 ] 985 Polygon centroid (ECEF): C = [ -4.6470e+06 2.5531e+06 -3.5333e+06 ] 986 Polygon centroid (Geo): Cg = [ -33.856926 151.215102 -4.9537e-04 ] 988 Figure 10 990 The point conversion for the polygon uses the final result, "Cg", 991 ignoring the altitude since the original shape did not include 992 altitude. 994 To convert this to a circle, take the maximum distance in ECEF 995 coordinates from the center point to each of the points. This 996 results in a radius of 99.1 meters. Confidence is unchanged. 998 5.2. Increasing and Decreasing Confidence 1000 Assuming that confidence is known to be 19% for Alice's location 1001 information. This is typical value for a three-dimensional ellipsoid 1002 uncertainty of normal distribution where the standard deviation is 1003 supplied in each dimension. The confidence associated with Alice's 1004 location estimate is quite low for many applications. Since the 1005 estimate is known to follow a normal distribution, the method in 1006 Section 4.4.2 can be used. Each axis can be scaled by: 1008 scale = erfinv(0.95^(1/3)) / erfinv(0.19^(1/3)) = 2.9937 1010 Ensuring that rounding always increases uncertainty, the location 1011 estimate at 95% includes a semi-major axis of 23.1, a semi-minor axis 1012 of 10 and a vertical axis of 86. 1014 Bob's location estimate covers an area of approximately 12600 square 1015 meters. If the estimate follows a rectangular distribution, the 1016 region of uncertainty can be reduced in size. To find the confidence 1017 that he is within the smaller area of the concert hall, given by the 1018 polygon [-33.856473, 151.215257; -33.856322, 151.214973; 1019 -33.856424, 151.21471; -33.857248, 151.214753; 1020 -33.857413, 151.214941; -33.857311, 151.215128]. To use this new 1021 region of uncertainty, find its area using the same translation 1022 method described in Section 4.1.1.2, which is 4566.2 square meters. 1023 The confidence associated with the smaller area is therefore 95% * 1024 4566.2 / 12600 = 34%. 1026 5.3. Matching Location Estimates to Regions of Interest 1028 Suppose than a circular area is defined centered at 1029 [-33.872754, 151.20683] with a radius of 1950 meters. To determine 1030 whether Bob is found within this area, we apply the method in 1031 Section 4.5. Using the converted Circle shape for Bob's location, 1032 the distance between these points is found to be 1915.26 meters. The 1033 area of overlap between Bob's location estimate and the region of 1034 interest is therefore 2209 square meters and the area of Bob's 1035 location estimate is 30853 square meters. This gives the probability 1036 that Bob is less than 1950 meters from the selected point as 67.8%. 1038 Note that if 1920 meters were chosen for the distance from the 1039 selected point, the area of overlap is only 16196 square meters and 1040 the confidence is 49.8%. Therefore, it is marginally more likely 1041 that Bob is outside the region of interest, despite the center point 1042 of his location estimate being within the region. 1044 6. Security Considerations 1045 This document describes methods for managing and manipulating 1046 uncertainty in location. No specific security concerns arise from 1047 most of the information provided. 1049 7. Acknowledgements 1051 Peter Rhodes provided assistance with some of the mathematical 1052 groundwork on this document. Dan Cornford provided a detailed review 1053 and many terminology corrections. 1055 8. Informative References 1057 [Convert] Burtch, R., "A Comparison of Methods Used in Rectangular 1058 to Geodetic Coordinate Transformations", April 2006. 1060 [GeoShape] 1061 Thomson, M. and C. Reed, "GML 3.1.1 PIDF-LO Shape 1062 Application Schema for use by the Internet Engineering 1063 Task Force (IETF)", Candidate OpenGIS Implementation 1064 Specification 06-142r1, Version: 1.0, April 2007. 1066 [I-D.thomson-geopriv-location-obscuring] 1067 Thomson, M., "Obscuring Location", draft-thomson-geopriv- 1068 location-obscuring-03 (work in progress), June 2011. 1070 [ISO.GUM] ISO/IEC, "Guide to the expression of uncertainty in 1071 measurement (GUM) ", Guide 98:1995, 1995. 1073 [NIST.TN1297] 1074 Taylor, B. and C. Kuyatt, "Guidelines for Evaluating and 1075 Expressing the Uncertainty of NIST Measurement Results ", 1076 Technical Note 1297, Sep 1994. 1078 [RFC3693] Cuellar, J., Morris, J., Mulligan, D., Peterson, J., and 1079 J. Polk, "Geopriv Requirements", RFC 3693, February 2004. 1081 [RFC3694] Danley, M., Mulligan, D., Morris, J., and J. Peterson, 1082 "Threat Analysis of the Geopriv Protocol", RFC 3694, 1083 February 2004. 1085 [RFC3825] Polk, J., Schnizlein, J., and M. Linsner, "Dynamic Host 1086 Configuration Protocol Option for Coordinate-based 1087 Location Configuration Information", RFC 3825, July 2004. 1089 [RFC4119] Peterson, J., "A Presence-based GEOPRIV Location Object 1090 Format", RFC 4119, December 2005. 1092 [RFC5139] Thomson, M. and J. Winterbottom, "Revised Civic Location 1093 Format for Presence Information Data Format Location 1094 Object (PIDF-LO)", RFC 5139, February 2008. 1096 [RFC5222] Hardie, T., Newton, A., Schulzrinne, H., and H. 1097 Tschofenig, "LoST: A Location-to-Service Translation 1098 Protocol", RFC 5222, August 2008. 1100 [RFC5491] Winterbottom, J., Thomson, M., and H. Tschofenig, "GEOPRIV 1101 Presence Information Data Format Location Object (PIDF-LO) 1102 Usage Clarification, Considerations, and Recommendations", 1103 RFC 5491, March 2009. 1105 [RFC6280] Barnes, R., Lepinski, M., Cooper, A., Morris, J., 1106 Tschofenig, H., and H. Schulzrinne, "An Architecture for 1107 Location and Location Privacy in Internet Applications", 1108 BCP 160, RFC 6280, July 2011. 1110 [Sunday02] 1111 Sunday, D., "Fast polygon area and Newell normal 1112 computation ", Journal of Graphics Tools JGT, 1113 7(2):9-13,2002, 2002, 1114 . 1116 [Vatti92] Vatti, B., "A generic solution to polygon clipping ", 1117 Communications of the ACM Vol35, Issue7, pp56-63, 1992, 1118 . 1120 [WGS84] US National Imagery and Mapping Agency, "Department of 1121 Defense (DoD) World Geodetic System 1984 (WGS 84), Third 1122 Edition ", NIMA TR8350.2, January 2000. 1124 Appendix A. Conversion Between Cartesian and Geodetic Coordinates in 1125 WGS84 1127 The process of conversion from geodetic (latitude, longitude and 1128 altitude) to earth-centered, earth-fixed (ECEF) Cartesian coordinates 1129 is relatively simple. 1131 In this section, the following constants and derived values are used 1132 from the definition of WGS84 [WGS84]: 1134 {radius of ellipsoid} R = 6378137 meters 1136 {inverse flattening} 1/f = 298.257223563 1138 {first eccentricity squared} e^2 = f * (2 - f) 1139 {second eccentricity squared} e'^2 = e^2 * (1 - e^2) 1141 To convert geodetic coordinates (latitude, longitude, altitude) to 1142 ECEF coordinates (X, Y, Z), use the following relationships: 1144 N = R / sqrt(1 - e^2 * sin(latitude)^2) 1146 X = (N + altitude) * cos(latitude) * cos(longitude) 1148 Y = (N + altitude) * cos(latitude) * sin(longitude) 1150 Z = (N*(1 - e^2) + altitude) * sin(latitude) 1152 The reverse conversion requires more complex computation and most 1153 methods introduce some error in latitude and altitude. A range of 1154 techniques are described in [Convert]. A variant on the method 1155 originally proposed by Bowring, which results in an acceptably small 1156 error, is described by the following: 1158 p = sqrt(X^2 + Y^2) 1160 r = sqrt(X^2 + Y^2 + Z^2) 1162 u = atan((1-f) * Z * (1 + e'^2 * (1-f) * R / r) / p) 1164 latitude = atan((Z + e'^2 * (1-f) * R * sin(u)^3) / 1165 (p - e^2 * R * cos(u)^3)) 1167 longitude = atan(Y / X) 1169 altitude = sqrt((p - R * cos(u))^2 + (Z - (1-f) * R * sin(u))^2) 1171 If the point is near the poles, that is "p < 1", the value for 1172 altitude that this method produces is unstable. A simpler method for 1173 determining the altitude of a point near the poles is: 1175 altitude = |Z| - R * (1 - f) 1177 Appendix B. Calculating the Upward Normal of a Polygon 1179 For a polygon that is guaranteed to be convex and coplanar, the 1180 upward normal can be found by finding the vector cross product of 1181 adjacent edges. 1183 For more general cases the Newell method of approximation described 1184 in [Sunday02] may be applied. In particular, this method can be used 1185 if the points are only approximately coplanar, and for non-convex 1186 polygons. 1188 This process requires a Cartesian coordinate system. Therefore, 1189 convert the geodetic coordinates of the polygon to Cartesian, ECEF 1190 coordinates (Appendix A). If no altitude is specified, assume an 1191 altitude of zero. 1193 This method can be condensed to the following set of equations: 1195 Nx = sum from i=1..n of (y[i] * (z[i+1] - z[i-1])) 1197 Ny = sum from i=1..n of (z[i] * (x[i+1] - x[i-1])) 1199 Nz = sum from i=1..n of (x[i] * (y[i+1] - y[i-1])) 1201 For these formulae, the polygon is made of points 1202 "(x[1], y[1], z[1])" through "(x[n], y[n], x[n])". Each array is 1203 treated as circular, that is, "x[0] == x[n]" and "x[n+1] == x[1]". 1205 To translate this into a unit-vector; divide each component by the 1206 length of the vector: 1208 Nx' = Nx / sqrt(Nx^2 + Ny^2 + Nz^2) 1210 Ny' = Ny / sqrt(Nx^2 + Ny^2 + Nz^2) 1212 Nz' = Nz / sqrt(Nx^2 + Ny^2 + Nz^2) 1214 B.1. Checking that a Polygon Upward Normal Points Up 1216 RFC 5491 [RFC5491] stipulates that polygons be presented in anti- 1217 clockwise direction so that the upward normal is in an upward 1218 direction. Accidental reversal of points can invert this vector. 1219 This error can be hard to detect just by looking at the series of 1220 coordinates that form the polygon. 1222 Calculate the dot product of the upward normal of the polygon 1223 (Appendix B) and any vector that points away from the center of the 1224 Earth from the location of polygon. If this product is positive, 1225 then the polygon upward normal also points away from the center of 1226 the Earth. 1228 The inverse cosine of this value indicates the angle between the 1229 horizontal plane and the approximate plane of the polygon. 1231 A unit vector for the upward direction at any point can be found 1232 based on the latitude (lat) and longitude (lng) of the point, as 1233 follows: 1235 Up = [ cos(lat) * cos(lng) ; cos(lat) * sin(lng) ; sin(lat) ] 1237 For polygons that span less than half the globe, any point in the 1238 polygon - including the centroid - can be selected to generate an 1239 approximate up vector for comparison with the upward normal. 1241 Authors' Addresses 1243 Martin Thomson 1244 Microsoft 1245 3210 Porter Drive 1246 Palo Alto, CA 94304 1247 US 1249 Email: martin.thomson@gmail.com 1251 James Winterbottom 1252 Andrew Corporation 1253 Andrew Building (39) 1254 Wollongong University Campus 1255 Northfields Avenue 1256 Wollongong, NSW 2522 1257 AU 1259 Email: james.winterbottom@andrew.com