Geographic math question

P Kishor punk.kish at GMAIL.COM
Thu Jun 28 18:17:51 EDT 2007


On 6/28/07, M5 <m5 at renefournier.com> wrote:
> Not being very strong at math, I have a little problem that I'm not sure how
> to solve. Maybe someone can help me.
>
> Basically, given a point (latitude, longitude) and a radius (100 meters)
> (think circle), I need to compute an equivalent square: That is, two points
> that would correspond to two corners of the square.
>
> From: 51, -114 100 meters
> To: 51.005, -114.005 NE corner
>  49.995, -113.995 SW corner
>
> Now, the above is not really accurate, of course, since the earth is
> spherical (well, at least most people think so), and I would like this
> computation to run in MySQL query, e.g.:
>
> UPDATE places SET ne_latitude = (*), ne_longitude = (*), sw_latitude = (*),
> sw_longitude = (*)
>
> In the above table, there are already three columns with the centre latitude
> and longitude and radius. Any ideas? Thanks.
>

use the great circle method for calculating the distance between two
lat/lon pairs. See below for an Oracle-specific proc I wrote a while
back. Modify it for your database, and adapt it for your problem.

============================
/*
The great_circle_distance function calculates distance between two points
given two sets of lon/lat and the radius of the earth. I call it lon/lat
instead of the more colloquial lat/lon to reflect that lon is actually x coord
and lat is y coord. Easy to forget this.

This function will return the distance in whatever units the earth radius is
provided. FYI, the approx. spherical radius of earth in miles is 3956, and in
kms is 6367. In any case, keep in mind that the radius of the
earth varies from equator to pole because the earth is not a sphere but an
oblate spheroid. It is possible to approximate adjustment to the radius of the
earth based on the latitude. For that, two radii would be required. The
calculation would be like so

earth_radius_equator = radius at equator
earth_radius_pole    = radius at pole
lat                  = latitude at which the radius is desired
earth_radius =  earth_radius_equator -
                ((earth_radius_equator - earth_radius_pole) * SIN(lat))

Suit yerself.

Puneet Kishor, punkish at eidesis.org
Aug, 2006
*/

CREATE OR REPLACE FUNCTION great_circle_distance(
  lon1         NUMBER,
  lat1         NUMBER,
  lon2         NUMBER,
  lat2         NUMBER,
  earth_radius NUMBER,
)

RETURN NUMBER IS

-- change this to suit yourself
pi      CONSTANT NUMBER := 3.14159;

/*
Most computers require the arguments of trigonometric functions to be
expressed in radians. To convert dms lon/lat to radians, multiply by dec2rad
constant declared above
*/
dec2rad CONSTANT NUMBER := pi / 180;

a     NUMBER;
c     NUMBER;
d_lat NUMBER;
d_lon NUMBER;

BEGIN

  d_lon := abs( (lon2 * dec2rad) - (lon1 * dec2rad) );
  d_lat := abs( (lat2 * dec2rad) - (lat1 * dec2rad) );

  a :=  POWER(SIN(d_lat / 2), 2) +
        (
          COS( (lat1 * dec2rad) ) *
          COS( (lat2 * dec2rad) ) *
          POWER(SIN(d_lon / 2), 2)
        );

  -- distance in radians
  c := 2 * ATAN2(SQRT(a), SQRT(1-a));

  RETURN ROUND(earth_radius * c, 1);

END great_circle_distance;
=================================

-- 
Puneet Kishor http://punkish.eidesis.org/
Nelson Inst. for Env. Studies, UW-Madison http://www.nelson.wisc.edu/
Open Source Geospatial Foundation http://www.osgeo.org/education/
S&T Policy Fellow, National Academy of Sciences http://www.nas.edu/
---------------------------------------------------------------------
collaborate, communicate, compete
=====================================================================



More information about the mapserver-users mailing list