basic projection question

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 question

Trying 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) is

d = 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

Gerald I. Evenden (gie@charon.er.usgs.gov) writes on 7 Jul 93:

You have totally lost me.

and you were the only reply :frowning:

Let me re-phrase the question. Forget about spheres
and ellipsoids and earth for a second. I am coding up
an algorithm that was written for points in a plane
(call it x-y). A line extends through two points,
(x1,y1) and (x2,y2), both of which are in the first
quadrant of a Cartesian coordinate system. This line
crosses the y-axis at some point. The angle 't'
between the line and the axis is what I am interested
in.

     | * (x2, y2)
     | /
     | /
     | * (x1, y2)
     | /
     | /
     -----------------------
     | /
     | /
     |t /
     | /
     |/
     |
    /|
   / |

Calculation of the angle is straightforward in the
x-y plane. Now enter the complexities of geo
referenced data, an funny-shaped earth, etc. How would
this angle be calculated using GRASS library
functions? (Assuming that the data is in one of
the supported coordinate systesms.)

I suspect that most users of this program that I'm
writing will never have anything by x-y
(unreferenced) data, but I was hoping to do what
is right-n-proper by including support for lat-long
data.

--Darrell

In info.grass.programmer you write:

Gerald I. Evenden (gie@charon.er.usgs.gov) writes on 7 Jul 93:

You have totally lost me.

and you were the only reply :frowning:

Let me re-phrase the question. Forget about spheres
and ellipsoids and earth for a second. I am coding up
an algorithm that was written for points in a plane
(call it x-y). A line extends through two points,
(x1,y1) and (x2,y2), both of which are in the first
quadrant of a Cartesian coordinate system. This line
crosses the y-axis at some point. The angle 't'
between the line and the axis is what I am interested
in.

    | * (x2, y2)
    | /
    | /
    | * (x1, y2)
    | /
    | /
    -----------------------
    | /
    | /
    |t /
    | /
    |/
    |
   /|
  / |

Calculation of the angle is straightforward in the
x-y plane. Now enter the complexities of geo
referenced data, an funny-shaped earth, etc. How would
this angle be calculated using GRASS library
functions? (Assuming that the data is in one of
the supported coordinate systesms.)

I suspect that most users of this program that I'm
writing will never have anything by x-y
(unreferenced) data, but I was hoping to do what
is right-n-proper by including support for lat-long
data.

--Darrell

Darrell, what Gerry says is important to understand. The definition of
a line between two points on a plane is unique, but not on on a curved
surface like a sphere or ellipsoid. The "natural" definition is the
line that is the shortest distance between the two points (what is
called a geodesic). If you define lines this way you get a different
angle than if you define the lines some other way- for example
rhumblines (these are lines that appear to be straight on a mercator
projection). Rhumblines are good for navigating based on angles to the
north pole, but they are oddly shaped when drawn on the surface of the
earth and the angle between rhumblines that leave one common points and
join two other points will not be the same as the geodesics which
connect these three same points.

Just code this for planimetric projections and document that it
isn't right for lat/lon. When someone provides a routines for
these calculations perhaps the routine can be included in your code.
--
Michael Shapiro shapiro@zorro.cecer.army.mil
U.S. Army CERL (217) 373-7277
P.O. Box 9005
Champaign, Ill. 61826-9005