David Torres<[email protected]> :
The key question to answer yourself before you adopt any of the many
strategies elucidated on the list is what you are going to do with a
separate list of identifiers for each observation on your original
dataset. Having this answer in hand, think about how the data should
be organized to do that efficiently. Then get the lists together.
Here is one way of doing the last in your original example (note I
changed only one thing, an equality test to an inequality), but I
cannot see the utility of this exercise, so it is unlikely to be what
you actually want:
g long rad10=.
g double dist1=.
g id1=""
merge using univpoints00_0
forv i=1/`np' {
qui vincenty `=latitude[`i']' `=longitude[`i']' sch_lat sch_lon, hav(dis)
g within10=dis<10
su within10, meanonly
qui replace rad10=r(sum) in `i'
su dis, meanonly
qui replace dist1=r(min) in `i'
qui levelsof unitid if dis<10, loc(vs)
foreach v of loc vs {
qui replace id1="`v' "+id1 in `i'
}
drop dis within10
}
On Sun, Aug 30, 2009 at 10:15 AM, David Torres<[email protected]> wrote:
> Austin, I've looked over the code you offered and it seems that what I'm
> getting is just the xy coordinates rather than what I need. There's got to
> be a way, somehow, to have stata return a list of all the unitids that are
> within the specified radius. When I run the code I have, all I get is a
> return of the unitid on the nearest school to my 65000 tracts. However,
> within the rad10 variable, several tracts have multiple schools within my
> radius of 10 miles. I need to get all of these unitids regardless of
> whether I end up with wide or long data; I can always reshape later.
>
> Is my question clearer?
>
> Well, any help you or anyone else can give is greatly appreciated. In the
> meantime, I'll keep working at it.
>
> Thanks
> David D. Torres
>
>
> Quoting Austin Nichols <[email protected]>:
>
>> David Torres<[email protected]> :
>> This is in reference to e.g.
>> http://stata.com/statalist/archive/2009-07/msg00083.html
>> http://stata.com/statalist/archive/2007-01/msg00098.html
>> I believe.
>>
>> Note that standard Stata commands like -summarize- and -replace-
>> inside the loop do all the work of saving characteristics of the
>> nearest, or k nearest, observations (and whatever else, such as mean
>> of Z over those with distance less than 10 units). No program like
>> -vincenty- (on SSC) does that; a -replace- command does it.
>> -vincenty- does not return "the min, max, mean, or whatever" and it
>> appears that -globdist- (which I found with -findit-) does only part
>> of what -vincenty- does (and the -globdist- help file contains the
>> caveat that "Testing of this program has been limited"). In other
>> words, add a -replace- to put the result you want in a variable.
>>
>> Maybe you want something like:
>>
>> g long rad10=.
>> g double dist1=.
>> g id1=""
>> g xy1=""
>> merge using univpoints00_0
>> forv i=1/`np' {
>> qui vincenty `=latitude[`i']' `=longitude[`i']' sch_lat sch_lon, hav(dis)
>> g within10=dis<10
>> su within10, meanonly
>> qui replace rad10=r(sum) in `i'
>> su dis, meanonly
>> qui replace dist1=r(min) in `i'
>> qui levelsof unitid if dis==dist1[`i'], loc(vs)
>> foreach v of loc vs {
>> qui replace id1="`v' "+id1 in `i'
>> su sch_lat if unitid==`v', meanonly
>> loc y1=r(mean)
>> su sch_lon if unitid==`v', meanonly
>> loc x1=r(mean)
>> qui replace xy1="(`x1',`y1') "+id1 in `i'
>> }
>> drop dis within10
>> }
>>
>> On Sat, Aug 29, 2009 at 4:26 PM, David Torres<[email protected]> wrote:
>>>
>>> I am trying to calculate multiple distances from one set of xy
>>> coordinates
>>> to another
>>> using vincenty. In addition to getting those distances, however, I need
>>> to
>>> somehow list
>>> the ids for the xy coordinates used to calculate the distance.
>>>
>>> Here is the piece of code I'm working with:
>>>
>>> code excerpt:
>>>
>>> g long rad10=.
>>> g double dist1=.
>>> g id1=""
>>> merge using univpoints00_0
>>>
>>> forv i=1/`np' {
>>> qui vincenty `=latitude[`i']' `=longitude[`i']' sch_lat sch_lon, hav(dis)
>>> g within10=dis<10
>>> su within10, meanonly
>>> qui replace rad10=r(sum) in `i'
>>> su dis, meanonly
>>> qui replace dist1=r(min) in `i'
>>> qui levelsof unitid if dis==dist1[`i'], loc(vs)
>>> foreach v of loc vs {
>>> qui replace id1="`v' "+id1 in `i'
>>> }
>>> drop dis within10
>>> }
>>>
>>> end excerpt
>>>
>>> As you can see, I'm getting a count of all sch_lat and sch_lon
>>> coordinates
>>> within a 10
>>> mile radius of the other set of latitude and longitude points. Then I
>>> return only the
>>> minimum distance. Will vincenty allow for returning more than just the
>>> min,
>>> max, mean,
>>> or whatever? Or is globdist a better choice for returning multiple
>>> distances and the ids
>>> attached to the point the distance to which is being calculated?
>>>
>>> Thanks
>>>
>>> David Torres
>>> UMich
*
* 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/
