Looping over geodist to create a distance matrix

[Jill Adona]

Dear all,

I'm trying to create a matrix from a dataset (since I don't know Mata yet; just playing on my strengths here): I have units with coordinates, and I want to get unit1's distance to itself, to unit2, to unit3, and so on in new variables d_u1u1, d_u1u2, d_u1u3, and so on.

First, looping over observations isn't a recommended practice and I can't get around it either, so I copied the lat and lon coordinates into lat2 and lon2.

levelsof unit_id, local(levels)
foreach 1 of local levels {
    geodist lat lon lat2 lon2, gen(d_`1')
} // THIS WORKED BUT OF COURSE IT JUST CONTAINED DISTANCES TO ITSELF, WHICH IS 0

I'm trying the code below but I get the error that "geodist lat/lon are not variables, nothing to generate"

foreach 1 of local levels {
    foreach x of local lon1 {
        foreach y of local lat1 {
            foreach a of local lon2 {
                foreach b of local lat2 {
                    geodist `y' `x' `b' `a', gen(d_`1')
                }
            }
        }
    }
}

Can anyone help?

[Jill Adona]

Got it! This code worked:

levelsof unit_id, local(levels)
levelsof lon, local(lon1)
levelsof lat, local(lat1)

foreach x of local lon1 {
    foreach y of local lat1 {
        foreach 1 of local levels {
            geodist `y' `x' lat2 lon2 if unit_id!="`1'", gen(d_`1')
        }
    }
}

[Robert Picard]

I'm not completely following your loop logic here but you can use cross to form all pairwise combination of units and then calculate the distance directly without any loops. If you want the results in a wide layout, you can then use reshape to create the variables. Here's an example:

* Example generated by -dataex-. To install: ssc install dataex
clear
input str3 unit_id double(lat lon)
"101"  38.39548681822478 -104.18357761720183
"102"  38.06754283901255 -102.14379541053262
"103" 37.546585177950746 -104.30864410382338
"104"  37.11422746935596 -104.83583446240563
"105" 40.475733076669314 -103.42037523327475
"106"  38.40341958395186 -103.51490262200869
"107" 37.284420367308336 -104.41756057006599
"108"  38.29347179610179 -108.32177649571459
"109" 39.220412647896254 -104.18045895713094
"110" 40.503964043592504 -102.89252782946723
end
save "dataex.dta", replace

rename (unit_id lat lon) =2
cross using "dataex.dta"
geodist lat lon lat2 lon2, gen(d)

drop lat* lon*
format %8.0g d
reshape wide d, i(unit_id) j(unit_id2) string
list, noobs

* spot check distance from "102" to "109"
list unit_id d109 in 2
use "dataex.dta", clear
geodist `=lat[2]' `=lon[2]' `=lat[9]' `=lon[9]'

and the results

. list, noobs

  +-------------------------------------------------------------------------------------------------------------+
  | unit_id      d101      d102      d103      d104      d105      d106      d107      d108      d109      d110 |
  |-------------------------------------------------------------------------------------------------------------|
  |     101         0   182.262   94.8633   153.385   240.119   58.4179   125.029   361.908   91.5769   259.121 |
  |     102   182.262         0   199.211   260.222   289.155   125.696     218.6   541.774   218.661   278.097 |
  |     103   94.8633   199.211         0   66.9736   334.153   117.934   30.6524    362.44   186.138   350.456 |
  |     104   153.385   260.222   66.9736         0   392.889    184.45   41.6615    334.08   240.735   412.393 |
  |     105   240.119   289.155   334.153   392.889         0    230.22   364.684   486.738   153.811   44.8608 |
  |-------------------------------------------------------------------------------------------------------------|
  |     106   58.4179   125.696   117.934    184.45    230.22         0   147.437   420.324   107.551   239.282 |
  |     107   125.029     218.6   30.6524   41.6615   364.684   147.437         0   361.639   215.897   381.096 |
  |     108   361.908   541.774    362.44    334.08   486.738   420.324   361.639         0    374.35   528.007 |
  |     109   91.5769   218.661   186.138   240.735   153.811   107.551   215.897    374.35         0   180.148 |
  |     110   259.121   278.097   350.456   412.393   44.8608   239.282   381.096   528.007   180.148         0 |
  +-------------------------------------------------------------------------------------------------------------+

. 
. * spot check distance from "102" to "109"
. list unit_id d109 in 2

     +-------------------+
     | unit_id      d109 |
     |-------------------|
  2. |     102   218.661 |
     +-------------------+

. use "dataex.dta", clear

. geodist `=lat[2]' `=lon[2]' `=lat[9]' `=lon[9]' 
WGS 1984 ellipsoid(6378137,298.257223563) distance = 218.66141 km

[Robert Picard]

For completeness, here's how to use a loop over observations to do the same thing:

* Example generated by -dataex-. To install: ssc install dataex
clear
input str3 unit_id double(lat lon)
"101"  38.39548681822478 -104.18357761720183
"102"  38.06754283901255 -102.14379541053262
"103" 37.546585177950746 -104.30864410382338
"104"  37.11422746935596 -104.83583446240563
"105" 40.475733076669314 -103.42037523327475
"106"  38.40341958395186 -103.51490262200869
"107" 37.284420367308336 -104.41756057006599
"108"  38.29347179610179 -108.32177649571459
"109" 39.220412647896254 -104.18045895713094
"110" 40.503964043592504 -102.89252782946723
end
save "dataex.dta", replace

forvalues i = 1/`=_N' {
    scalar lat0 = lat[`i']
    scalar lon0 = lon[`i']
    local unit = unit_id[`i']
    geodist lat0 lon0 lat lon, gen(d`unit')

}
format %8.0g d*
list unit_id d*, noobs

and the results:

. list unit_id d*, noobs

  +-------------------------------------------------------------------------------------------------------------+
  | unit_id      d101      d102      d103      d104      d105      d106      d107      d108      d109      d110 |
  |-------------------------------------------------------------------------------------------------------------|
  |     101         0   182.262   94.8633   153.385   240.119   58.4179   125.029   361.908   91.5769   259.121 |
  |     102   182.262         0   199.211   260.222   289.155   125.696     218.6   541.774   218.661   278.097 |
  |     103   94.8633   199.211         0   66.9736   334.153   117.934   30.6524    362.44   186.138   350.456 |
  |     104   153.385   260.222   66.9736         0   392.889    184.45   41.6615    334.08   240.735   412.393 |
  |     105   240.119   289.155   334.153   392.889         0    230.22   364.684   486.738   153.811   44.8608 |
  |-------------------------------------------------------------------------------------------------------------|
  |     106   58.4179   125.696   117.934    184.45    230.22         0   147.437   420.324   107.551   239.282 |
  |     107   125.029     218.6   30.6524   41.6615   364.684   147.437         0   361.639   215.897   381.096 |
  |     108   361.908   541.774    362.44    334.08   486.738   420.324   361.639         0    374.35   528.007 |
  |     109   91.5769   218.661   186.138   240.735   153.811   107.551   215.897    374.35         0   180.148 |
  |     110   259.121   278.097   350.456   412.393   44.8608   239.282   381.096   528.007   180.148         0 |
  +-------------------------------------------------------------------------------------------------------------+

.

[Mike Lacy]

Deleted because my essentially identical response crossed in the ether with Robert Picard's.

[Jill Adona]

Thank you, Robert! I will try that code. For some reason while I'm checking the resulting distances they aren't correct...