Date: Wed, 7 Jul 93 02:37:21 -0500

Message-Id: <9307070737.AA12583@bushland.ecn.purdue.edu>

From: Darrell McCauley <mccauley@ecn.purdue.edu>

To: grassp-list@max.cecer.army.mil

Subject: basic projection questionTrying to make my programs more lat-long friendly...

G_begin_distance_calculations();

I need to calculate the angle of a line made between two

points (x1,y1) and (x2,y2).

First, you need to define what you mean by a "line" between two point

when dealing with geographic coordinate: a rhumb line (constant

azimuth) or geodesic (shortest distance).

The distance between the points (the hypotenuse in

geometric/planimetric terms) isd = G_distance (x1, y1, x2, y2);

In a planimetric projection, the angle (in degrees) would be:

What do you mean by "planimetric projection?" Projections have nothing

to do with rhumb or geodesic distances---if that's what we are dealing

with.

angle = 180.0 / 3.14159 * acos ( (x1-x2) /d )

but with lat-long I don't believe that this is the case

(because the distance between two longitudes is different,

depending upon what latitude you are at). Would the

correct way to compute this angle be:angle = 180.0 / 3.14159

* acos ((G_distance(x1,y1,x2,y1)+G_distance(x1,y2,x2,y2)) /2.0

/ d);(that is, averaging the two distances at lat1 and lat2 to come

up with the numerator for the inverse cosine quotient).

This approach seems sort of naive but I cannot think of

anything better at the moment.

You have totally lost me.

Is this correct approach? Is there a library function that I am

missing which would do what I want?--Darrell

I am not sure that I understand your problem. But you seem to be needing

the formulary for either a rhumbline or geodesic and I suspect the latter.

The program geod distributed with the PROJ.4 system will give you

forward and inverse geodesic computations with forward and back azimuths

at each point and great circle distances for both spherical and elliptical

Earths. The elliptical formulary is too complex to post but the spherical

formulae are basically:

distance = 2*R*asin(sqrt(sin^2((lat2-lat1)/2) + cos(lat1)*cos(lat2)*

sin^2((lon2-lon1)/2)))

Azi1-2 = atan2(cos(lat2)*sin(lon2-lon1), cos(lat1)*sin(lat2)-

sin(lat1)*cos(lat2)*cos(lon2-lon1))

where lat1,lon1 and lat2,lon2 are geographic coordinates of points,

R is Earth's radius. To get back azimuth from point 2 to point 1,

use second equation and with reversed coodinantes. For more details

and verification of the above see Snyder's "Map Projections---A Working

Manual", USGS Prof.Paper 1395, p. 30. For rhumb line computations, see

p. 46 of same publication.

Gerald (Jerry) I. Evenden Internet: gie@charon.er.usgs.gov

voice: (508)563-6766 Postal: P.O. Box 1027

fax: (508)457-2310 N.Falmouth, MA 02556-1027