Given a line on the earth's surface, how do I plot a line perpendicular to it?

Given a line on the earth's surface, how do I plot a line perpendicular to it?

Apologies if this is a very simple question. I thought this would be a straightforward task, but it is proving counter-intuitive.

I start with the blue line in the figure below (see link -- I'm unable to upload figures yet). I found a perpendicular line by calculating the gradient of the blue line (m), then plotting another line (green) with gradient -1/m. When I plot the lines in Matlab (using 'plot' and 'axis equal'), they look perpendicular, as expected.

http://imgur.com/7qMkx

However, when I export these lines to Google Earth (using the KML Toolbox), they no longer look perpendicular (see link below; the shorter line is the blue line from the previous figure).

http://imgur.com/ncJQ7

I understand that strange things happen on curved surfaces, but I thought the lines should at least look perpendicular locally. I suspect this has something to do with the projection in Google Earth -- in particular, the fact that the grid cells appear to have roughly similar side lengths, yet the longitudinal edge has length = 1 degree, whereas the latitudinal edge has length = 0.5 degree.

So, in summary:

• is my method of finding a perpendicular line valid on a curved surface? (i.e., plotting a line with gradient -1/m)
• in the Google Earth image, do the perpendicular lines look as expected, or is something strange going on?

UPDATE:

To provide more context: I am looking at radar data taken from a plane. The multi-colored area is the 'swath', where observations have been recorded. The blue line I started with in the explanation above is parallel to the swath: this is the aircraft's flight line (the plane was moving in the roughly south-west direction). The radar looks in the direction orthogonal to the flight line, on the left. I am trying to draw a line perpendicular to the flight line; this should be the direction the radar is looking, and should cut the swath neatly. As you can see, this is not the case. • Actually, Google portrays the angle correctly, whereas your Matlab plot does not... – whuber Jul 16 '12 at 20:23
• Thanks for the reply, whuber. Does that mean the line I have calculated is not, in fact, perpendicular to the original? Or does it mean that the way I view the lines in Matlab is somehow wrong? I used the 'axis equal' command in Matlab to convince myself that the lines are perpendicular, but I expect this would not properly represent the curvature of the earth. – Paul Keating Jul 16 '12 at 20:37
• Your first interpretation is correct, Paul. Using (lon, lat) as if they were Cartesian coordinates distorts the directions so grievously that most angles are visibly distorted. Therefore, using the Cartesian formula -1/m won't even get you close to a right angle except in a couple of lucky directions (or at the Equator). To provide an alternative method, it would help to know what kind of "line" you want (geodesic or rhumb line, for instance), your purpose in plotting it, and how long it would need to be. (Short lengths allow simple approximations to be used.) – whuber Jul 16 '12 at 20:59
• BTW, you can learn a lot about these issues by reading some of the threads on our site about conformal projections. – whuber Jul 16 '12 at 21:01
• Thanks very much, whuber. Very interesting. I have provided more context in the problem description, above. The length I need is relatively short (the width of the swath, shown above, which is a few tenths of a degree) so any suggestions on valid approximations would be much appreciated. Thanks also for the link -- I'll check it out. – Paul Keating Jul 16 '12 at 21:18

An elegant principle provides a simple answer:

All points on a smooth curved surface are flat at a sufficiently large scale.

This means that after affine change of coordinates (usually involving just a rescaling of one of them), we can use formulas of Euclidean geometry, such as the Pythagorean Theorem for computing distances and the negative-reciprocal-slope formula for finding perpendiculars.

With latitude and longitude coordinates on the sphere (away from the poles, where longitude becomes singular), all we need to do is rescale the east-west direction to reflect the shrinking length of a degree of longitude as one approaches the poles. With a spherical model of the earth, that shrinkage is given by the cosine of the latitude. This is merely a change in the plot's aspect ratio, nothing more.

This works for regions that extend no more than a few degrees of latitude north-south and do not approach either pole.

Therefore, all you have to do is:

1. Multiply all longitudes by the cosine of a typical latitude.

2. Compute the perpendicular line.

3. Undo the coordinate adjustment.

For example, suppose the plane's track took it from (lon, lat) = (-78, 40) to (-79, 41). We may take a typical latitude to lie between 40 and 41, such as 40.5.

Step 1 The adjusted coordinates are (-78 * cos(40.5), 40) = (-59.31167, 40) and (-79 * cos(40.5), 41) = (-60.07207, 41).

Step 2 The question proposes doing this using a negative-reciprocal-slope method. That would be correct, but it will fail in some cases (where the slope is infinite). It is more general and more powerful to use vector arithmetic. Here's how the calculation goes.

The direction vector for the flight's path is the displacement from its beginning to its end,

v =  (-60.07207, 41) - (-59.31167, 40)
=  (-0.7604, 1.0).

Turning any vector (x,y) at right angles clockwise produces (y,-x), whence a perpendicular direction to the right is

w = (1.0, 0.7604).

According to the Pythagorean Theorem, the length of this vector is the square root of the sum of the squares of its coefficients,

|w| = sqrt(1^2 + 0.7604^2) = 1.256268

Let's move, say, 0.2 degrees along this vector from the starting point of the plane's flight. The start is at (-59.31167, 40) and the displacement is 0.2 / |w| times w, ending up at

(-59.31167, 40) + 0.2 / 1.256268 * (1.0, 0.7604) = (-59.15247  40.12106).

Step 3 To undo the adjustment, divide the first coordinates of any resulting points by the same cosine used in Step 1:

(-59.15247/cos(40.5), 40.12106) = (-77.79064, 40.12106)

If you plot these points using a 1:1 aspect ratio, the angle will appear to be obtuse rather than a right angle. But if you change the aspect ratio to 1:cos(40.5) (about 4:3), the angle will correctly appear to be 90 degrees. When you plot the points using any conformal projection--including Google's Mercator--the angle will also be correct.

• This is fantastic, whuber. Thanks very much. It's a very clear explanation of where I am going wrong, and how I can fix the problem. I've ended up using Matlab's 'azimuth' and 'reckon' functions to get the perpendicular line, but it's good to see a more detailed, first-principles explanation of how to do this. Many thanks. – Paul Keating Jul 17 '12 at 13:41