500论坛币求stata代码，关于坐标距离计算

（1）有3个基准经纬度坐标，分别是（110.23 36.23；111.33 33.45；116.80 40.35）
（2）有9个id的坐标，如示例数据。
（3）在第一个基准坐标（110.23 36.23）半径5km之内，或者，在第二个基准坐标（111.33 33.45）半径3km之内，或者在第三个基准坐标半径1km之内，wanted=1，否则wanted=0.
（4）实际数据量，id有几万条，基准坐标有几千个。请问代码如何编写呢？

1. * Example generated by -dataex-. For more info, type help dataex
   
2. clear
   
3. input float(id lon lat) str3 wanted是否在范围内
   
4. 1 116.78294 39.60197 "？"
   
5. 2 115.97266 40.44241 "？"
   
6. 3 116.11597 39.89557 "？"
   
7. 4 116.49863 39.80755 "？"
   
8. 5  116.8082  40.3568 "？"
   
9. 6 116.13786 39.71119 "？"
   
10. 7 116.64865 40.12924 "？"
    
11. 8 115.91908 40.36859 "？"
    
12. 9 116.66607 39.64987 "？"
    
13. end

复制代码

diannaoasd 发表于 2025-7-26 08:27
stata提示：refcoord not foundinvalid syntax

The code is fine, but there seems a problem with the "复制代码" button.
Don't hit that button. To copy the code, highlight it, right-click, and select 'Copy.'

Alternatively, use the attached do file: dist_20250724.do (1.03 KB)

2025-7-26 23:24:30 上传

stata_code

Once the code is executed, the output will appear as follows:

Try the following code to compute distance:

1. clear all
   
2. input float(id lon lat)
   
3. 1 116.78294 39.60197
   
4. 2 115.97266 40.44241
   
5. 3 116.11597 39.89557
   
6. 4 116.49863 39.80755
   
7. 5 116.80820 40.35680
   
8. 6 116.13786 39.71119
   
9. 7 116.64865 40.12924
   
10. 8 115.91908 40.36859
    
11. 9 116.66607 39.64987
    
12. end
    
13. 
    
14. scalar lon1 =110.23 //point1
    
15. scalar lat1 = 36.23
    
16. 
    
17. scalar lon2 =111.33 //point2
    
18. scalar lat2 = 33.45
    
19. 
    
20. scalar lon3 =116.80 //point3
    
21. scalar lat3 = 40.35
    
22. 
    
23. scalar earthradius =6371
    
24. scalar radians =c(pi) /180
    
25. 
    
26. forvalues i =1/3 {
    
27.     generate double dist`i'= ///
    
28.         2 *earthradius * ///
    
29.         asin(sqrt((1 - ///
    
30.         cos((lat`i' -lat) *radians) + ///
    
31.         cos( lat`i'       *radians) * ///
    
32.         cos(         lat  *radians) * ///
    
33.     (1 -cos((lon`i' -lon) *radians)))/2))
    
34. }
    
35. 
    
36. label variable dist1 "distance from point1 in km"
    
37. label variable dist2 "distance from point2 in km"

复制代码

See also the Haversine formula: https://en.wikipedia.org/wiki/Haversine_formula

Here is an updated version:

1. clear all
   
2. input float(id lon lat wanted)
   
3. 1  116.78294 39.60197 0
   
4. 2  115.97266 40.44241 0
   
5. 3  116.11597 39.89557 0
   
6. 4  116.49863 39.80755 0
   
7. 5  116.80820 40.35680 0
   
8. 6  116.13786 39.71119 0
   
9. 7  116.64865 40.12924 0
   
10. 8  115.91908 40.36859 0
    
11. 9  116.66607 39.64987 0
    
12. end
    
13. 
    
14. matrix refcoord =110.23, 36.23, 5 \ ///
    
15.                  111.33, 33.45, 3 \ ///
    
16.                  116.80, 40.35, 1
    
17. matrix colname refcoord ="lon" "lat" "radius"
    
18. //this refcoord matrix may be imported from external files
    
19. 
    
20. scalar earthradius =6371
    
21. scalar radians =c(pi) /180
    
22. 
    
23. forvalues i =1/`=rowsof(refcoord)' {
    
24.     loc loni =el(refcoord,`i',1)
    
25.     loc lati =el(refcoord,`i',2)
    
26.     loc radi =el(refcoord,`i',3)
    
27. 
    
28.     capture drop dist
    
29.     generate double dist= ///
    
30.         2 *earthradius * ///
    
31.         asin(sqrt((1 - ///
    
32.         cos((`lati' -lat) *radians) + ///
    
33.         cos( `lati'       *radians) * ///
    
34.         cos(         lat  *radians) * ///
    
35.     (1 -cos((`loni' -lon) *radians)))/2))
    
36. 
    
37.     replace wanted =`i' if wanted==0 & dist<=`radi'
    
38. }

复制代码

kirin_guess 发表于 2025-7-24 14:51
Here is an updated version:

stata提示：refcoord not foundinvalid syntax