Re: st: calculating latitude and longitude of nearby point

["Dimitriy V. Masterov" <dvmaster@gmail.com>]

I think I may have figured out how to do this using plane geometry,
which does not seem like a terrible approximation. My trig is a bit
rusty, so I could not get the coordinates to come out using the
formulas Scott referred to.

Here's the code. I tried to verify the results using geodist. I am off
by 1.6 meters. The two approaches do define the angles differently.

******************************************
clear;
set obs 1;

/* Find the lat/lon of a point one mile north of Stata HQ */
gen lat  = 30.55233;
gen lon = -96.24564;
gen d = 1.609344; // 1 mile in km

/* This method seems to work */
/* The constants 110540 and 111320 reflect the earth's oblateness
(polar and equatorial circumferences are different). */
gen theta = _pi/2; // This formula measures angles counterclockwise
from due east, so N = 90 deg or pi/2 in radians
gen lat_new  = lat  + (d*sin(theta))/110540;
gen lon_new  = lon + (d*cos(theta))/(111320*cos(_pi*lat/180));

geodist lat lon lat_new lon_new, gen(distance);

list;
drop *_new distance theta;


/* This does not work as implemented here, formulas from
http://www.movable-type.co.uk/scripts/latlong.html */
gen theta = 0; // North is 0 deg/radians, and pi/2 or 90 deg is E in
radians (different from above)
gen lat_new = asin(sin(lat*_pi/180)*cos(d/6371) +
cos(lat*_pi/180)*sin(d/6371)*cos(theta));
gen lon_new = lon +
atan2(sin(theta)*sin(d/6371)*cos(_pi*lat/180),cos(d/6371)-sin(_pi*lat/180)*sin(_pi*lat_new/180));

geodist lat lon lat_new lon_new, gen(distance);

list;

******************************************
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/