Creating a distance matrix from geocoordinates

[Jan Grbs]

Hi,

I have a data set of about 2500 hospitals and their geocoordinates (longitude, latitude). In order to measure the number of hospitals within a 10 km radius of each hospital, I want to create a distance matrix. I found some code on the internet using spatwmat to create a distance matrix.

However, I do constantly receive error messages, e.g. "matrix operation not found". Can anyone help me on this or suggest a different way to create such a distance matrix?

Many thanks in advance.
Jan

[Mike Lacy]

-geonear- has worked very well for me. I don't know if it will return the matrix you want, but it will identify all observations within a specified distance of each observation. It's amazingly fast, I might add. See -ssc describe geonear-.

[Jan Grbs]

Thanks Mike! The geonear function sounds good. However, I am not very familiar with Stata and thus having some problems using the function. So my data set looks like this (2444 rows in total):

Primary Key Longitude Latitude
2600000660 46829723 9855609
26010002399 54790663 9426379
... .... ...

What I would like to achieve is calculating the number of hospitals within 10km, 30km and 50km distance for each Primary Key and saving the data in newly generated variables c1, c2 and c3. I would be very thankful if you could help me with that.

Many thanks,
Jan

[Mike Lacy]

For your case, the "neighbor file" and the "base file" would be the same Stata data file, i.e. "hospitals.dta"
I presume you don't want to count each subject hospital in your totals.

So, with hospitals.dta resident in Stata, and using the variable names you show, I think you would want:

geonear PrimaryKey Latitude Longitude using "hospitals.dta", ///
   wide neighbors(PrimaryKey Latitude Longitude) ignoreself radius(10)

I have not used -geonear- for exactly this purpose, but my understanding is that this would return in variables named nid1,...,nidMax the PrimaryKey of each of the hospitals within 10km of the subject hospital. It would also return in the variables km_to_nid1, ... km_to_nidMax to identify the distances. (Max is my notation for whatever the maximum number of hospitals is that is found for among all your subject hospitals.) Some of your subject hospitals would have only (say) 1, while others might have 3 or 4. The higher numbered neighbor variables would be missing for hospitals with fewer neighbors, so you can count the number of neighbors, your desired value, with:

egen NHospitalsWithin10km = rownomiss(km_to*)

You can then drop all the nid* and km* variables, and repeat the preceding command at 30km, etc.

I have not tried this out, but my code should be close.

[Robert Picard]

If your data is really as you show (and note our strong preference for data examples generated using dataex, see item 12.2 of the FAQ), then you can't use geonear because these are not latitudes and longitudes. These are most likely (x,y) cartesian coordinates generated using some map projection. Because of the inevitable distortions needed to fit a spherical surface on a 2d map, cartesian distances will be inaccurate. If you still want to calculate the distances, you have the problem of finding the unit of measurement. Anyway, the number of point is small enough that you can use a brute force approach to calculate the distance between all points.

* Example generated by -dataex-. To install: ssc install dataex
clear
input double(PrimaryKey Longitude Latitude)
 2600000660 46829723 9855609
26010002399 54790663 9426379
12345677888 54799999 9426666
end
format %15.0g PrimaryKey
save "temp_statalist.dta", replace

* form all pairwise combination of points
rename * *0
cross using "temp_statalist.dta"

* cartesian distance
gen distance = sqrt((Longitude-Longitude0)^2 + (Latitude-Latitude0)^2)

* so what's the unit of distance?
sort PrimaryKey0 d PrimaryKey
list, sepby(PrimaryKey0)

Alternatively, you can use GIS software to convert these points back to lat/lon or you can find another source that is already in lat/lon.

[Jan Grbs]

Thanks a lot! I managed to solve the problem by using the geonear function. And yes, I had to put my data in the right format first as geonear requires the values for longitude and latitude to be between -90 and 90 I think.

[Kees Maat]

Dear Robert (or someone else),

I was looking for a way to create just a distance matrix from coordinates, so cross was very helpful for that. But I have an additional question:

I used simply 1, 2 and 3 for the PrimaryKey0 and PrimaryKey and sorted the data as a distance matrix,
so sort PrimaryKey0 PrimaryKey. I then get:

1       1        0  
1       2        7972503  
1       3        7981810  
2       1        7972503  
2       2        0  
2       3        9340.41  
3       1        7981810  
3       2        9340.41  
3       3        0

But what I need (in txt or csv or Excel) is a real matrix with distances between all locations:

         1        2         3
1        0        7972503   7981810  
2        7972503  0         9340.41  
3        7981810  9340.41   0

Do you know how to do that?

Thanks, Kees

[Kees Maat]

The answer on my previous question is reshape wide.

[Diti Roy]

I am trying to perform geo near. I want to look at the border districts' effect (spillover effect) i.e I want to find the non-treated neighboring districts for my treated districts. However, the treatment of the district varies over time. In other words, districts come in for treatment and go out over time. Now, I have got two datasets with the centroid values - one includes only treated districts and the other non-treated districts. After which I perform the geo near using the following codes with 2000 km?

geonear district x_stub y_stub using "Non- treated districts with centroids.dta", n(district1 x_stub1 y_stub1) ign long within(2000) near(2)

However, I am getting the following syntax.

nbor latitude var x_stub1 must be between -90 and 90
r(198);

Can anybody please help with these concerns?

[Andrew Musau]

Check for missing values. Otherwise, post a reproducible example.

[Diti Roy]

there is no missing value

[Andrew Musau]

As I said in #10, we need a data example that replicates the issue. It does not have to be your full or actual dataset.

[Diti Roy]

I have several years of panel data with nth districts. My treatment varies over time (i.e some of the districts come in and go out of the treatment).To solve this, I have merged the latitudes and the longitudes (generated from shapefile by using the gencentroids(stub) command) with the dataset containing the treatment and year variable in the panel data by using districts_id- from which I created two files to perform geonear - one with the banned districts and the other with the non-ban districts.


I have not pasted the entire tabulated dataset

ssc install dataex

. tab x_stub

x-coordinat |
e of area |
centroid | Freq. Percent Cum.
------------+-----------------------------------
69.92829 | 20 2.81 2.81
69.93822 | 20 2.81 5.63
70.36705 | 20 2.81 8.44
70.76935 | 20 2.81 11.25
70.90977 | 5 0.31 11.56
71.15736 | 20 2.81 14.38
71.47205 | 5 0.31 14.69
71.53689 | 20 2.81 17.50
71.79604 | 20 2.81 20.31
71.92759 | 20 2.81 23.13
72.2083 | 5 0.31 23.44
72.22892 | 20 2.81 26.25
72.25824 | 20 2.81 29.06
72.6628 | 18 2.38 31.44
72.78568 | 5 0.31 31.75



tab y_stub

y-coordinat |
e of area |
centroid | Freq. Percent Cum.
------------+-----------------------------------
8.296923 | 17 1.06 1.06
8.603258 | 4 0.25 1.31
8.835245 | 17 1.06 2.38
9.106371 | 4 0.25 2.63
9.438544 | 4 0.25 2.88
9.619133 | 17 1.06 3.94
9.794156 | 4 0.25 4.19
10.01389 | 4 0.25 4.44
10.13618 | 17 1.06 5.50
10.47388 | 4 0.25 5.75
10.64858 | 17 1.06 6.81
10.79416 | 4 0.25 7.06
10.868 | 17 1.06 8.13



.use "treated districts with centroids.dta", clear
.geonear district x_stub y_stub using "Non- treated districts with centroids.dta", n(district1 x_stub1 y_stub1) ign long within(2000) near(2)

nbor latitude var x_stub1 must be between -90 and 90
r(198);

[Andrew Musau]

A data example means the raw data, not the results of a tabulation. You can, e.g., copy and paste the results of

use "treated districts with centroids.dta", clear
contract district x_stub y_stub
sort district
dataex district x_stub y_stub in 1/100
use "Non- treated districts with centroids.dta", clear
contract district1 x_stub1 y_stub1
sort district1
dataex district1 x_stub1 y_stub1 in 1/100

[Diti Roy]

Thanks. Here it is

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(district_id x_stub y_stub)
  1  77.57578 14.484938
  2  79.01996  13.46324
  3  78.78391 14.461876
  4   82.0033 17.190853
  5  80.07764 16.291048
  6  80.79324 16.568815
  7  78.01077 15.522762
  8  79.74493  14.42524
  9  79.51898 15.611388
 10  83.66933 18.534872
 11  82.71606 17.879093
 12  81.40018 16.892174
 22  86.95314  25.04594
 23  84.60091 26.878845
 24  86.06253 26.107756
 25  84.92094 24.800936
 26  86.29471  25.16689
 27  85.40443  26.18871
 28  85.32538  25.35474
 29  87.53761  25.91772
 30     86.77  25.99135
 31  84.55655 26.146526
 32  84.03043  25.17149
 39  72.22892  22.71804
 40 71.157364 21.376474
 41  71.92759 24.218157
 42   73.2238  21.75263
 43  71.79604 21.631826
 44  73.71906 20.782715
 45 72.662796 23.203684
 46 69.938225  22.25322
 47  70.36705 21.291815
 48 72.901024 22.702017
 49  69.92829  23.64138
 50  72.25824   23.6403
 51  73.84151  22.89135
 52  70.76935 22.250036
 53  73.14674 23.633554
 54  73.34554  21.23379
 55  71.53689 22.785044
 56  73.53912 22.231274
 57  73.12912  20.58801
 58  77.11015  30.38144
 59  77.03268  28.17054
 60   75.5223  29.28937
 61  76.31522 29.425783
 62  76.75643  29.74056
 63  76.15585  28.32321
 64  76.68118  28.83104
 85  86.40241  23.74969
 86  85.58877  24.11807
 87  84.10754 24.010754
 88  84.89087 23.071264
 89    87.245 24.493023
 90  85.85558  22.54394
 91   77.4726 12.913547
 92  74.82811 16.113764
 93  76.44615 15.046112
 94  77.22865 17.942135
 95  75.81337 16.560215
 96  75.68145 13.428985
 97  76.41669 14.223225
 98  75.79681 12.323008
 99  75.41306 15.144032
100   76.8756 17.054512
101  76.10504 12.978118
102  78.01498 13.370895
103  76.79259 12.604687
104  76.75268 12.079443
105  74.60293 14.774975
106  76.63239  15.87672
107 75.326836 14.036293
108  75.08897 13.120307
109  76.94049 13.513035
110  76.47502  9.438544
111  75.46492 12.129185
112  76.56656  10.01389
113  76.90978  9.106371
114  76.92366  9.794156
115 75.952866 11.549033
116  76.15894 11.137506
117  76.56069 10.794156
118  76.31778 10.473877
119  77.00609  8.603258
183  86.78558 21.288477
184  84.08356  20.26519
185  83.34911  20.67327
186  86.14165  20.55069
187  85.16365  20.92387
188   84.4937 19.466806
189  82.93784 20.060373
190  85.70293 21.532385
191  82.68578 18.988888
192  86.40359  21.90428
193  85.40851 20.060816
194  84.01537   21.4042
195  84.52103   22.0683
196  74.89443 31.517076
197  75.11371  30.18027
198  74.59668 30.522146
end

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(district_id1 x_stub1 y_stub1)
  1  77.57578 14.484938
  2  79.01996  13.46324
  3  78.78391 14.461876
  4   82.0033 17.190853
  5  80.07764 16.291048
  6  80.79324 16.568815
  7  78.01077 15.522762
  8  79.74493  14.42524
  9  79.51898 15.611388
 10  83.66933 18.534872
 11  82.71606 17.879093
 12  81.40018 16.892174
 13  92.67226 24.683865
 14  92.57874  26.71558
 15  90.23173 26.291393
 17  94.96799  27.42658
 18  93.13989 26.046877
 19  93.02378   25.3485
 20  92.68436  26.22977
 21  94.25702  26.72941
 22  86.95314  25.04594
 23  84.60091 26.878845
 24  86.06253 26.107756
 25  84.92094 24.800936
 26  86.29471  25.16689
 27  85.40443  26.18871
 28  85.32538  25.35474
 29  87.53761  25.91772
 30     86.77  25.99135
 31  84.55655 26.146526
 32  84.03043  25.17149
 33  81.26424 19.232664
 34  82.28456 22.336386
 35  81.10991   21.2673
 36  83.55159 22.384794
 37  82.21187  21.00315
 38  82.93346 23.386633
 44  73.71906 20.782715
 45 72.662796 23.203684
 58  77.11015  30.38144
 59  77.03268  28.17054
 60   75.5223  29.28937
 61  76.31522 29.425783
 62  76.75643  29.74056
 63  76.15585  28.32321
 64  76.68118  28.83104
 65  76.67768   31.3584
 66  76.46109  32.70605
 67  76.33959   31.9426
 68  78.40654 31.572264
 69  77.39387 31.888023
 70  77.75434 32.413715
 71  77.65334  31.17899
 72  77.00079  31.61929
 73  76.91825  31.03793
 74   77.4459  30.65933
 75  75.12834  33.75286
 76   74.4601  34.38276
 77  75.84465 33.417263
 78  74.82651  32.69138
 79  75.57309 32.555298
 80  78.02759 34.283478
 81  74.30234  33.64733
 82  74.38713  33.20358
 83  74.87225  34.08359
 84  75.07751  33.02404
 85  86.40241  23.74969
 86  85.58877  24.11807
 87  84.10754 24.010754
 88  84.89087 23.071264
 89    87.245 24.493023
 90  85.85558  22.54394
 91   77.4726 12.913547
 92  74.82811 16.113764
 93  76.44615 15.046112
 94  77.22865 17.942135
 95  75.81337 16.560215
 96  75.68145 13.428985
 97  76.41669 14.223225
 98  75.79681 12.323008
 99  75.41306 15.144032
100   76.8756 17.054512
101  76.10504 12.978118
102  78.01498 13.370895
103  76.79259 12.604687
104  76.75268 12.079443
105  74.60293 14.774975
106  76.63239  15.87672
107 75.326836 14.036293
108  75.08897 13.120307
109  76.94049 13.513035
110  76.47502  9.438544
111  75.46492 12.129185
112  76.56656  10.01389
113  76.90978  9.106371
114  76.92366  9.794156
115 75.952866 11.549033
116  76.15894 11.137506
117  76.56069 10.794156
118  76.31778 10.473877
end