Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am looking for the most accurate way of calculating the view angle (Azimuth, Elevation) from a reference point (my location) A to another point B. I have the Lat, Lon, and Altitude for each point.

Currently I am converting to local ENU coordinates using the Wikipedia formulas for ECEF to ENU and using that to calculate the angles.

Az = atan2(E,N)
El = atan(U / sqrt(E^2 + N^2)

Is there a more accurate way of calculating these Angles? Maybe a different coordinate system I should be projecting to? Is this the correct type of Angle, should a rhumb line azimuth or something else be used instead?

I have also looked at using Vincenty's formula for calculating the azimuth but it does not take into account the altitude of each point. A change in Altitude using ENU seems to slightly change the Azimuth which seems correct(?). When comparing Vincenty Azimuth with the ENU Azimuth on the surface (Altitude 0), the answers are extremely close but start to differ if the points are not on the surface.

share|improve this question

1 Answer 1

up vote 4 down vote accepted

For accurate calculations, convert (lat, lon, elevation) directly to earth-centered (x,y,z). (If you don't do this, you need to retain additional information about the local normal ["up"] directions in order to compute angles accurately at nonzero elevations.)

Elevation

Given two points (x,y,z) and (x',y',z') in an earth-centered coordinate system, the vector from the first to the second is (dx,dy,dz) = (x'-x, y'-y, z'-z), whence the cosine of the angle made to the normal at (x,y,z) is the inner product of the unit length versions of those vectors:

Cos(elevation) = (x*dx + y*dy + z*dz) / Sqrt((x^2+y^2+z^2)*(dx^2+dy^2+dz^2))

Obtain its principal inverse cosine. Subtract this from 90 degrees if you want the angle of view relative to a nominal horizon. This is the "elevation."

Azimuth

A similar calculation obtains the local direction of view ("azimuth"). We need a level vector (u,v,w) that points due north. One such vector at the location (x,y,z) is (-z*x, -z*y, x^2+y^2). (The inner product of these two vectors is zero, proving it is level. Its projection onto the Equatorial plane is proportional to (-x,-y) which points directly inward, making it the projection of a north-pointing vector. These two calculations confirm that this is indeed the desired vector). Therefore

Cos(azimuth) = (-z*x*dx - z*y*dy + (x^2+y^2)*dz) / Sqrt((x^2+y^2)(x^2+y^2+z^2)(dx^2+dy^2+dz^2))

We also need the sine of the azimuth, which is similarly obtained once we know a vector pointing due East (locally). Such a vector is (-y, x, 0), because it clearly is perpendicular to (x,y,z) (the up direction) and the northern direction. Therefore

Sin(azimuth) = (-y*dx + x*dy) / Sqrt((x^2+y^2)(dx^2+dy^2+dz^2))

These values enable us to recover the azimuth as the inverse tangent of the cosine and sine.


Example

A pilot in an airplane flying flying west at 4000 meters, located at (lat, lon) = (39, -75), sees a jet far ahead located at (39, -76) and flying at 12000 meters. What is are the angles of view (relative to the level direction at the pilot's location)?

The XYZ coordinates of the airplanes are (x,y,z) = (1285410, -4797210, 3994830) and (x',y',z') = (1202990, -4824940, 3999870), respectively (in the ITRF00 datum, which uses the GRS80 ellipsoid). The pilot's view vector therefore is (dx,dy,dz) = (-82404.5, -27735.3, 5034.56). Applying the formula gives the cosine of the view angle as 0.0850706. Its inverse cosine is 85.1199 degrees, whence the elevation is 4.88009 degrees: the pilot is looking up by that much.

A north-pointing level vector is (-5.13499, 19.1641, 24.6655) (times 10^12) and an east-pointing level vector is (4.79721, 1.28541, 0) (times 10^6). Therefore, applying the last two formulas, the cosine of the azimuth equals 0.00575112 and its sine equals -0.996358. The ArcTangent function tells us the angle for the direction (0.00575112, -0.996358) is 270.331 degrees: almost due west. (It's not exactly west because the two planes lie on the same circle of latitude, which is curving toward the North pole: see Why is the 'straight line' path across continent so curved? for an extended explanation.)

By the way, this example confirms we got the orientation correct for the azimuth calculation: although it was clear the east-pointing vector was orthogonal to the other two vectors, until now it was not plain that it truly points east and not west.

share|improve this answer
    
Forgive me if I am missing something but I am looking for the Azimuth Angle and Elevation Angle from A to B where the Azimuth Angle is 0 when pointing directly North and positive when pointing East, and the Elevation Angle is negative when pointing down and positive when pointing up. It seems like this provides Elevation angle but how do I calculate the Azimuth? –  user66332 Apr 23 '13 at 17:24
    
Also, in your example, xyz and x'y'z' are the same. –  user66332 Apr 23 '13 at 17:33
    
Re the first comment: I have added the azimuth calculation. Re the second comment: thank you for catching that! (It was a manual copying error.) I have fixed it. –  whuber Apr 23 '13 at 17:53
    
Are you sure the North vector is calculated correctly? When I compare your East vector to the East vector calculated from the Wiki page, the vectors are the same. They are not the same for the North vector. Wiki north unit vector = -0.16288010267506386, 0.6078768187253738, 0.7771459614569709 Your north unit vector = 0.16221935131335216, 0.6054108610722955, 0.779206372763453 –  user66332 Apr 24 '13 at 19:22
1  
ok i got the answer for whoever wants to do this (convert lat/lon to cartesian) in the future. uses trac.osgeo.org/proj/wiki/man_cs2cs : cs2cs +proj=latlong +datum=WGS84 +units=m +to +proj=geocent +datum=WGS84… not knowing the GIS vocabulary is painful –  mga Aug 5 at 16:30

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.