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 (-zx, -zy, 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 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.
Best Answer
QGIS >= 2.18
New function (after this question was posted)
azimuth()has become available, and made this calculation easier.If we have two point layers, one of which is your samples layer and the other is the specific third point. In the example below they correspond to:
TargetPoints_4326: sample points. Only unique id ("id") field is required. I have only three ("id"=1,2,3) pink-colored points in this example.ObsPoints_4326: your observation points to measure the angle. I have three overlapping points (green color). This layer has to have "EPSG_code" field to define the CRS you want to work, along with "Start" and "End" fields as pointers to the"id"field of the sample points layer.Then open the attribute table of your observation points layer (e.g.
ObsPoints_4326) and create a new field ("angle") with following expression:Output ("angle") field is measured anti-clockwise (e.g. first row on the attached picture shows angle= 27.185). If you want clockwise- measurement, please edit the above expression to switch "Start" and "End".
Note: This example used Lambert Conformal Conic (EPSG:3034) which will be good if your points are spread across large area in Europe. If your study area is local, UTM would do.