Re: st: calculating latitude and longitude of nearby point

[Robert Picard <picard@netbox.com>]

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

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