Notice: On March 31, it was **announced** that Statalist is moving from an email list to a **forum**. The old list will shut down on April 23, and its replacement, **statalist.org** is already up and running.

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

From |
Robert Picard <picard@netbox.com> |

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

**Follow-Ups**:**Re: st: calculating latitude and longitude of nearby point***From:*"Dimitriy V. Masterov" <dvmaster@gmail.com>

**References**:**st: calculating latitude and longitude of nearby point***From:*"Dimitriy V. Masterov" <dvmaster@gmail.com>

**Re: st: calculating latitude and longitude of nearby point***From:*Scott Merryman <scott.merryman@gmail.com>

**RE: st: calculating latitude and longitude of nearby point***From:*Nick Cox <n.j.cox@durham.ac.uk>

**Re: st: calculating latitude and longitude of nearby point***From:*"Dimitriy V. Masterov" <dvmaster@gmail.com>

- Prev by Date:
**RE: st: Select the last 4 numbers out of 8...** - Next by Date:
**st: Merging based on dates** - Previous by thread:
**Re: st: calculating latitude and longitude of nearby point** - Next by thread:
**Re: st: calculating latitude and longitude of nearby point** - Index(es):