Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down at the end of May, and its replacement, statalist.org is already up and running.


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

Re: st: calculating latitude and longitude of nearby point


From   "Dimitriy V. Masterov" <dvmaster@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: calculating latitude and longitude of nearby point
Date   Thu, 9 Jun 2011 10:06:09 -0400

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/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index