Geodata: Calculating distance from a point to a border

[Jason DeBacker]

Dear Stata Users,
I am interested in using a geographic regression discontinuity design where there is a discontinuity at a U.S. state border. My problem is that I do not know how to calculate the shortest distance from a point to a polygon. In ArcGIS, I see that there is a tool called “Near” that can determine the distance from a point to the nearest point on a polygon – this is exactly what I would like to do, however I cannot use ArcGIS and need to do the task in Stata (or SAS).

After much searching, I have not found a similar function in Stata to the ArcGIS “Near” function. I have seen the nearest.ado function, which is an efficient way to calculate the nearest point to another given a set of points. I need to calculate the distance to the border of just one state, so it might be very feasible to create a set of points that are on the border of this state (say a point at each half-mile along the border) and use nearest.ado to find the nearest border point to each non-border point. Is there a way to take a shapefile (e.g., the Census TIGER/Line state shapefile) and create dataset with the set of points along the border of a state?

Any other suggestions on how to calculate the distance from a point to a border?

Many thanks for your help!

Sincerely,
Jason

[Robert Picard]

I've done something like this in the past. I had a set of grid points spaced 1 mile apart in Texas and I needed to find, for each grid point, the closest point on the Mexican border. Here's the code that processes the Mexico border shapefile:

* There are more than one polygon. Identify each polygon. The first one
* represent mainland Mexico. 

    use "mexico_coor.dta", clear
    
    rename _X lon
    rename _Y lat

    // a new polygon starts with missing values for lat/lon
    gen poly = sum(lon == .)
    gen id = _n
    sort poly id
    by poly: assert lon == . if _n == 1
    by poly: assert lon[2] == lon[_N]
    by poly: assert lat[2] == lat[_N]
    
    sum
    
    
* Drop a point that is an obvious error in the data

    list in 319378/319381
    drop in 319380

    
* We will need to measure the distance between points in Texas to the
* Mexican border. Make sure that the line segments are short. Target
* 0.1 miles between each point. Create new points for polygon segments
* that are too long

    // distance between consecutive points
    gen double lat2 = lat[_n+1]
    gen double lon2 = lon[_n+1]
    geodist lat lon lat2 lon2, gen(d) mi
    
    // number of points needed to keep the segments at < 0.1 mi
    gen n = 1 + int(d/.1)
    
    expand n if !mi(d)
    sort id
    
    // interpolate linearly each intermediary points
    by id: gen double lat3 = lat + (_n-1) * (lat2 - lat) / n
    by id: gen double lon3 = lon + (_n-1) * (lon2 - lon) / n
    
    // use these new coordinates
    drop lat lon
    rename lat3 lat
    rename lon3 lon
    replace id = _n
    
    // double check; ignore missing distances (separate polygons)
    keep poly id lat lon
    gen double lat2 = lat[_n+1]
    gen double lon2 = lon[_n+1]
    geodist lat lon lat2 lon2, gen(d) mi
     assert round(d,.1) <= .1 if !mi(d)


* Clean-up and save

    keep poly id lat lon

and here's the code to find the nearest border point

    use "mexico/mexican_border_points.dta"
    keep if poly == 1    // keep only mainland Mexico
    tempfile mexico
    qui save "`mexico'"
    
    use "tx_mile_grid_inside.dta"
    
    geonear gridid lat lon using "`mexico'", ///
        n(id lat lon) mi long near(1)

Note that both geodist and geonear are available from SSC.

[Jason DeBacker]

Robert,
Thank you! This method looks like it'll work well for my purposes - I'll give it a go.

Sincerely,
Jason

[Jesse Wursten]

Robert Picard