Minimum distance from point to linestring on earth's surface

General Usage
,
1

Say I have a shapefile representing the highway network in South America and want to find the minimum geographic distance from that network to another point on the continent. Is there no simple way to do this with the available geospatial packages ?

I understand that I can use nearestPoints() in LibGEOS.jl to find the point in the network closest to the point of interest under a naive Eucludean interpretation of latitude and longitude, and then use another package to compute the distance between the two points in a coordinate reference system such as WGS 1984. But technically that’s not quite right since the point returned by nearestPoints() might not actually be the nearest point in WGS 1984.

In R, the sf package provides st_distance(), which at least computes great circle distances between various kinds of geometry objects.

1 Like
2 
"""
    haversine_distance_km(lon1, lat1, lon2, lat2)

Calculate geodesic distance between two points using the Haversine formula.
"""
function haversine_distance_km(lon1::Float64, lat1::Float64, lon2::Float64, lat2::Float64)
    dlat = deg2rad(lat2 - lat1)
    dlon = deg2rad(lon2 - lon1)
    
    a = sin(dlat/2)^2 + cos(deg2rad(lat1)) * cos(deg2rad(lat2)) * sin(dlon/2)^2
    c = 2 * atan(sqrt(a), sqrt(1-a))
    
    return EARTH_RADIUS_KM * c
end

This ought to work given the easternmost, e.g., node on the network.

1 Like
3

The GMT.jl module mapproject can do that. Unfortunately, being a such multi-functionality program I have not yet translated the terse GMT syntax to the verbose version of GMT.jl. The work is done by option -G and in GMT.jl it would be something like

using GMT

D = gmtread("the_shape_file");
mapproject(D, G=(lon,lat))      # Here (lon,lat) is your "point on the continent"
1 Like
4

@technocrat, this is the distance between two points, right? Not the minimum distance between a point and a network/MultiLineString.

5

It is two points. That’s why I suggested that one node on the network would need to be chosen. Am I missing something in what you’re trying to do?

6

I’m not sure what a MultiLineString is. But in order to find the distance between, e.g., a line segment and a point, pointwise distances are not sufficient. You also need to determine the actual closest point along the segment.

The missing part is that the closest point might not be an endpoint.

7

As I said, GMT does it. Finds the distance and the closest point.

8

Yes. There’s a non-trivial algorithmic step in picking the nearest point in the network and doing that properly in spherical/ellipsoidal coordinates.

9

Thank you @joa-quim. I tried this. My actual shapefile has 542 entries, each represented with a vector of coordinate pairs. E.g., the first entry has 94 coordinate pairs. Do I understand right that mapproject() is augmenting this data structure with a 3rd column, which holds the distance to each coordinate pair?

If so, then this is not quite the correct algorithm, right? The minimum of the distances to the vertices is not necessarily the minimum distance to the network.

10

Sorry, it’s not option G it’s option L, which in GMT.jl is aliased to dist2line. (read the linked description of that option)

I know what you are after (I needed that operation several times before) and I’m sure GMT does it with ellipsoidal accuracy if needed. See this small example that is computing the closest point along the great circle line between [0 0] and [1 1] from point at [0.5 0.5]

julia> mapproject([0.5 0.5], dist2line=[0.0 0; 1 1])

1×5 GMTdataset{Float64, 2}
 Row │ col.1  col.2    col.3     col.4     col.5
─────┼───────────────────────────────────────────
   1 │   0.5    0.5  4.46581  0.499972  0.500029

The results are in meters and col3=dist, col4=lon_along_line; col5=lat_along_line

In your real case the command should be

mapproject([lon lat], dist2line=D)

and since you have 542 segments, you’ll have 542 rows. You can test it with a single segment (the first) with

mapproject([lon lat], dist2line=D[1])
1 Like