Approximate distance between two points (longitude / latitude) without Haversine

Question

How can I judge from the long/lat coordinates of a point A alone whether it is possibly/certainly not within a short distance of x meters from another coordinate B, using R (without applying geosphere::distHaversine or any other advanced calculation)?

"Certainly" in this case means that the long/lat of A tells me that it cannot be within the critical distance to B, while "possibly" means the complementary case.

Example:

x <- 100
point_A <- c("long"=8.678180,"lat"=50.114390)
point_B <- data.frame("id"=c(1,2,3,4),"long"=c(8.678459,8.618162,8.610878,7.208705),
"lat"=c(50.114124,50.226831,50.230530,50.799290))
point_B_certainly_not_in_distance <- subset(point_B, ...) # wanted
point_B_possibly_in_distance <- subset(point_B, ...) # wanted

In this example, only coordinate with id 1 should be identified as possibly within distance, whereas id 2, 3 and 4 should be identified as certainly not in distance.

To clarify the purpose: Imagine that point_B has >10^6 rows. In that case, a simplified "pre-subsetting" before applying the Haversine formula makes great sense.

Answer 1

1 degree of latitude is the same distance (1/360 circumference) everywhere on the globe, so you can cut on latitude difference. Two points that have a latitude difference of 1 degree can't be nearer than 1/360 of the circumference.

Having done that, you're left with points in a narrow latitude band, you could work out the radius of that band, approximate it as a circle, map the points to the circle (a sine and a cosine and a multiply or two), and then use the circle distance as a threshold, although points nearer the poles could have under-threshold great circle distances but over-threshold approximate latitude band distances. So the threshold may need some correction which would depend on latitude. Now its getting complex, and the benefit isn't clear.

Another method with a couple of sines and cosines would be to compute the 3d position of a point, then use Pythagoras to compute the straight-line Euclidean distance between points through the globe rather than the great circle distance. This distance will always be less than the surface distance, so if Euclidean distance > threshold then the great circle distance will be greater than the threshold. Drop the point.

In any event you don't need to compute the distance, which involves computing an arcsin (in the formula on Wikipedia), you only need compute the angle, which is a couple or so sines and cosines. Work out the angle threshold for your distance threshold and then you only do arcsin on the points in your threshold to get the distance.

Answer 2

I had a similar issue recently (150 point_As and 700,000 point_Bs) solved using the point in polygon concept and the sp and rgeos packages. Basically, you can turn your point_A into a ring of radius x and then 'subset' the point_B points which fall within that ring.

When you have your data in the right format, all you need is the %over% command from rgeos.

First, convert your points to spatialPoints using the sp package and set an appropriate coordinate reference (CRS). I just picked the WGS 84 projection as a 'standard' and appears fine for this:

 library(dplyr)
    library(sp)

x <- 13500  ##setting this as the radius of interest

point_A <- data.frame("long"=8.678180,"lat"=50.114390) %>% SpatialPoints(proj4string = CRS("+init=epsg:4326"))
point_B <- data.frame("id"=c(1,2,3,4),"long"=c(8.678459,8.618162,8.610878,7.208705),
                      "lat"=c(50.114124,50.226831,50.230530,50.799290))
point_B <- SpatialPointsDataFrame(coords=point_B[,c("long","lat")],data=point_B,proj4string = CRS(proj4string(point_A)))

Next, you will need to project the coordinates to a planar reference system.

A projection is a formula used to convert long/lat coordinates into a flat coordinate system that you can use on paper or a computer screen. It's usually done from a geographic coordinate system (long/lat).

I just picked ESPG 2192 here as it seemed close enough (Western Europe). You can use sites such as http://spatialreference.org to find the best projection to use for your needs.

After projecting, you can run rgeos::gBuffer to create the ring around your point. After that, everything should be converted back to the original long/lat CRS.

library(rgeos) 
point_A <- spTransform(point_A,CRS("+init=epsg:2192")) 
Point_A_buffer <- gBuffer(point_A,width = x,quadsegs = 100)
point_A <- spTransform(point_A,CRS("+init=epsg:4326")) 
Point_A_buffer <- spTransform(Point_A_buffer,CRS("+init=epsg:4326")) 

Then the %over% operator from sp will indicate which points are encompassed by the ring. From there you can subset as intended.

point_in_zone <- point_B %over% Point_A_buffer %>% data.frame("id"=point_B@data[,"id"],"test"=.)

> point_in_zone
  id in_radius
1  1         1
2  2         1
3  3        NA
4  4        NA

point_B_possibly_in_distance <- point_B[!is.na(point_in_zone$test),"id"]
point_B_certainly_not_in_distance <- point_B[is.na(point_in_zone$test),"id"]

You can verify the calculation has worked as intended through a plot or through checking the great circle distances on a subset.

library(leaflet)

leaflet(Point_A_buffer) %>% addTiles %>% addMarkers(data=point_A) %>% 
  addCircleMarkers(data=point_B,radius = 2,fillOpacity=1) %>% 
  addPolygons(fillColor = "red",weight=1,
              color="red") 

Answer 3

My original answer here was not correct; sorry. Here comes my solution, well-tested for my use case. Expecting a vector of decimal latitude values and a distance in meters, calculateProperties returns the decimal long/lat equivalent of one such distance at the location with the narrowest latitude. This means you can use these two values to filter existing long/lat points and include at least all points which are definitely within the range (i.e. making a square subset).

decimalDegreesToDMS <- function(degree){
  d <- as.integer(degree)
  m <- as.integer((degree - d) * 60)
  s <- (degree - d - m/60) * 3600
  return(
    list("d"=d,"m"=m,"s"=s)
  )
}

DMStoAngle <- function(d,m,s){
  return( d + m/60 + s/3600 )
}

# http://formulas.tutorvista.com/math/degrees-to-radians-formula.html
angleToRadian <- function(angle){
  return( angle * pi / 180 )
}

calculateProperties <- function(latitude_vector,distance){
  if(abs(max(latitude_vector)) >= abs(min(latitude_vector))){
    narrowest_lat <- max(latitude_vector)
  } else {
    narrowest_lat <- min(latitude_vector)
  }

  # https://gis.stackexchange.com/questions/142326/calculating-longitude-length-in-miles
  dms <- decimalDegreesToDMS( narrowest_lat )
  angle <- DMStoAngle( dms$d, dms$m, dms$s )
  radians <- angleToRadian(angle)
  x_per_lat_degree <- distance / (111.32 * 10^3)
  x_per_long_degree <- distance / (40075 * 10^3 * cos( radians ) / 360)

  return( list("narrowest_lat"=narrowest_lat,
               "x_per_lat_degree"=x_per_lat_degree,
               "x_per_long_degree"=x_per_long_degree) )
}

# subset existing points by meter distance
points <-data.frame("longitude"=c(8.678180,8.678459,8.618162,8.610878,7.208705),
         "latitude"=c(50.114390,50.114124,50.226831,50.230530,50.799290))
cp <- calculateProperties(points[,"latitude"],1000)
x_per_long_degree <- cp[["x_per_long_degree"]]
x_per_lat_degree <- cp[["x_per_lat_degree"]]
square_subset <- subset(points, longitude < points[1,"longitude"]+x_per_long_degree
                        & longitude > points[1,"longitude"]-x_per_long_degree
                        & latitude < points[1,"latitude"]+x_per_lat_degree
                        & latitude > points[1,"latitude"]-x_per_lat_degree)

(Please note that this does not cover the original question completely, but it is easy to solve the equation for distance.)