Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

# Re: st: calculating latitude and longitude of nearby point

 From Robert Picard To statalist@hsphsun2.harvard.edu Subject Re: st: calculating latitude and longitude of nearby point Date Thu, 9 Jun 2011 13:46:12 -0400

```I had to do this recently so I implemented a quick and dirty port of
the equations from:

http://www.movable-type.co.uk/scripts/latlong.html

Robert

* --------------------- begin example ---------------------
* This is the formula for calculating the coordinates of a
* destination point given distance and bearing from start
* point. Note Excel's ATAN2(x,y) is specified atan2(y,x) in
* Stata.
*
*lat2: =ASIN(SIN(lat1)*COS(d/R) + COS(lat1)*SIN(d/R)*COS(brng))
*lon2: =lon1 + ATAN2(COS(d/R)-SIN(lat1)*SIN(lat2), SIN(brng)*SIN(d/R)*COS(lat1))
* (source: http://www.movable-type.co.uk/scripts/latlong.html)

capture program drop destpoint
program define destpoint

args dlat1 dlon1 d b lat2 lon2

local R 6371

tempname lat1 lon1 brng
gen double `lat1' = `dlat1' * _pi / 180
gen double `lon1' = `dlon1' * _pi / 180
gen double `brng' = `b' * _pi / 180

gen double `lat2' = asin(sin(`lat1') * cos(`d'/`R') + ///
cos(`lat1') * sin(`d'/`R') * cos(`brng'))

gen double `lon2' = `lon1' + atan2(sin(`brng') * sin(`d'/`R') *
cos(`lat1'), ///
cos(`d'/`R') - sin(`lat1') * sin(`lat2'))

qui replace `lat2' = `lat2' * 180 / _pi
qui replace `lon2' = `lon2' * 180 / _pi

end

* Create 360 neighbor points

clear
set obs 360
gen nid = _n

// set reference point at 0,0. Use a radius of 1km
gen double clat = 0
gen double clon = 0
gen double bearing = _n - 1

destpoint clat clon 10 bearing nlat nlon

// -geodist- is available from SSC
geodist clat clon nlat nlon, gen(d) sphere
assert abs(10 - d) < 1e-10
drop d

* Make a nice plot

scatter nlat nlon, ylabel(none) xlabel(none) aspect(1)

* --------------------- begin example ---------------------

On Thu, Jun 9, 2011 at 10:06 AM, Dimitriy V. Masterov
<dvmaster@gmail.com> wrote:
> 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
> 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/
>

*
*   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/
```