Determine If One Coordinate Is in Radius of Another

Determine if one coordinate is in radius of another

This post shows how to do this in SQL Server.

And here is how to do it in MySQL:

SELECT ((ACOS(SIN($lat * PI() / 180) * SIN(lat * PI() / 180) + 
         COS($lat * PI() / 180) * COS(lat * PI() / 180) * COS(($lon - lon) * 
         PI() / 180)) * 180 / PI()) * 60 * 1.1515) AS distance 
FROM members 
HAVING distance<='10' ORDER BY distance ASC

Formula to determine whether a radius is within another radius

Let circle radii are R and r.

Calculate distance d between circle centers using haversine formula.

Compare d with radii:

d > R + r:  circles don't intersect
Abs(R-r) <= d <= R + r: circles do intersect
Abs(R-r) > d : one circle lies inside another

Check if point is inside a circle

Function to calculate the distance between two coordinates (converted to C# from this answer):

double GetDistance(double lat1, double lon1, double lat2, double lon2) 
{
    var R = 6371; // Radius of the earth in km
    var dLat = ToRadians(lat2-lat1);
    var dLon = ToRadians(lon2-lon1); 
    var a = 
        Math.Sin(dLat/2) * Math.Sin(dLat/2) +
        Math.Cos(ToRadians(lat1)) * Math.Cos(ToRadians(lat2)) * 
        Math.Sin(dLon/2) * Math.Sin(dLon/2);

    var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1-a)); 
    var d = R * c; // Distance in km
    return d;
}

double ToRadians(double deg) 
{
    return deg * (Math.PI/180);
}

If the distance between the two points is less than the radius, then it is within the circle.

Checking whether coordinates fall within a given radius

Updated to more scalable solution:

The previous answer (still included below) is not suited for large data sets. The reason is that we need to compute the distance for each pair of shops and bus. Therefore both the memory and computation scale as O(N*M) for N shops and M buses. A more scalable solution uses a data structure such as a KD-Tree to perform nearest neighbor search for each shop. The advantage here is that the computational complexity becomes O(M*logM) for building the KD-Tree for the bus stops and O(N*logM) for searching the nearest neighbor for each shop.

To do this, we can use nn2 from the RANN package. The complication here is that nn2 deals only with Euclidean distances and does not know anything about lat/long. Therefore, we need to convert the lat/long coordinates to some map projection (i.e. UTM) in order to use it correctly (i.e., in order to compute the Euclidean distance between shops and bus stops correctly).

Note: The following borrows heavily from Josh O'Brien's solutions for determining the UTM zone from a longitude and for converting lat/long to UTM, so he should take a bow.

## First define a function from Josh OBrien's answer to convert
## a longitude to its UTM zone
long2UTM <- function(long) {
  (floor((long + 180)/6) %% 60) + 1
}

## Assuming that all points are within a zone (within 6 degrees in longitude),
## we use the first shop's longitude to get the zone.
z <- long2UTM(shops[1,"long"])

library(sp)
library(rgdal)

## convert the bus lat/long coordinates to UTM for the computed zone
## using the other Josh O'Brien linked answer
bus2 <- bus
coordinates(bus2) <- c("long", "lat")
proj4string(bus2) <- CRS("+proj=longlat +datum=WGS84")

bus.xy <- spTransform(bus2, CRS(paste0("+proj=utm +zone=",z," ellps=WGS84")))

## convert the shops lat/long coordinates to UTM for the computed zone
shops2 <- shops
coordinates(shops2) <- c("long", "lat")
proj4string(shops2) <- CRS("+proj=longlat +datum=WGS84")

shops.xy <- spTransform(shops2, CRS(paste0("+proj=utm +zone=",z," ellps=WGS84")))

library(RANN)

## find the nearest neighbor in bus.xy@coords for each shops.xy@coords
res <- nn2(bus.xy@coords, shops.xy@coords, 1)
## res$nn.dist is a vector of the distance to the nearest bus.xy@coords for each shops.xy@coords
## res$nn.idx is a vector of indices to bus.xy of the nearest bus.xy@coords for each shops.xy@coords
shops$Bus_Stop <- res$nn.dists <= 500
shops$Bus_ID <- ifelse(res$nn.dists <= 500, bus[res$nn.idx,"Bus_Stop_ID"], NA)

Although more complicated, this approach is much better suited for realistic problems where you may have large numbers of shops and bus stops. Using the same supplied data:

print(shops)
##  Shop_ID       lat     long Bus_Stop Bus_ID
##1       1 -34.03935 18.61796     TRUE      A
##2       2 -33.92782 18.41052    FALSE   <NA>

---

You can do this using the package geosphere. Here, I'm assuming that your first data frame is named bus, and your second data frame is named shops:

library(geosphere)
g <- expand.grid(1:nrow(shops), 1:nrow(bus))
d <- matrix(distGeo(shops[g[,1],c("long","lat")], bus[g[,2],c("long","lat")]),
            nrow=nrow(shops))
shops$Bus_Stop <- apply(d, 1, function(x) any(x <= 500))
shops$Bus_ID <- bus[apply(d, 1, function(x) {
                                  c <-which(x <= 500)
                                  if(length(c)==0) NA else c[1]
                                }), "Bus_Stop_ID"]
print(shops)
##  Shop_ID       lat     long Bus_Stop Bus_ID
##1       1 -34.03935 18.61796     TRUE      A
##2       2 -33.92782 18.41052    FALSE   <NA>

Notes:

1. We first use expand.grid to enumerate all pair combinations of shops and bus stops. These are ordered by shops first.
2. We then compute the distance matrix d using geosphere::distGeo. Note here that the input expects (lon, lat) coordinates. distGeo returns distances in meters. The resulting d matrix is now(shops) by now(bus) so that each row gives the distance from a shop to each bus stop.
3. We then see if there is a bus stop within 500 meters of each shop by applying the function any(x <= 500) for each row x in d using apply with MARGIN=1.
4. Similarly, we can extract the column of d (corresponding to the row in bus) for the first shop that is within 500 meters using which instead of any in our applied function. Then use this result to select the Bus_Stop_ID from bus.

By the way, we don't have to apply the condition x <= 500 twice. The following will also work:

shops$Bus_ID <- bus[apply(d, 1, function(x) {
                                  c <-which(x <= 500)
                                  if(length(c)==0) NA else c[1]
                                }), "Bus_Stop_ID"]
shops$Bus_Stop <- !is.na(shops$Bus_ID)

and is more efficient.

Data:

bus <- structure(list(Bus_Stop_ID = structure(1:2, .Label = c("A", "B"
), class = "factor"), lat = c(-34.04199, -33.92312), long = c(18.61747, 
18.44649)), .Names = c("Bus_Stop_ID", "lat", "long"), class = "data.frame",  row.names = c(NA, 
-2L))

shops <- structure(list(Shop_ID = 1:2, lat = c(-34.03935, -33.92782), 
long = c(18.617964, 18.41052), Bus_ID = structure(c(1L, NA
), .Label = c("A", "B"), class = "factor"), Bus_Stop = c(TRUE, 
FALSE)), .Names = c("Shop_ID", "lat", "long", "Bus_ID", "Bus_Stop"
), row.names = c(NA, -2L), class = "data.frame")

Determine if point intersects 35km radius around another point? Possible in Linq?

Sure. Assume you have a function that computes the Haversine distance between two Positions (consisting of a latitude and longitude coordinate). If you don't you can find one here. Then simply use the function as the selector in a Where clause. If using LINQ to SQL, you'll need to materialize them to your Position objects so that the you can use the Haversine function on them as LINQ to objects; there isn't a translation to SQL, though you could probably create a table-valued function that does the same thing if you really don't want to return all the points first.

var origin = new Position( 47.6, 122.3 );
var close = positions.Where( p => Haversine.Distance( origin, p, DistanceType.Km ) <= 35 );

Finding all points in certain radius of another point

You're still going to have to iterate through every point, but there are two optimizations you can perform:

1) You can eliminate obvious points by checking if x1 < radius and if y1 < radius (like Brent already mentioned in another answer).

2) Instead of calculating the distance, you can calculate the square of the distance and compare it to the square of the allowed radius. This saves you from performing expensive square root calculations.

This is probably the best performance you're gonna get.

How to check if a certain coordinates fall to another coordinates radius using PHP only

Thanks for the help. Below is an example function that takes two sets of longitude and latitude co-ordinates and returns the distance between the two.

function getDistance($latitude1, $longitude1, $latitude2, $longitude2) {  
  $earth_radius = 6371;

  $dLat = deg2rad($latitude2 - $latitude1);  
  $dLon = deg2rad($longitude2 - $longitude1);  

  $a = sin($dLat/2) * sin($dLat/2) + cos(deg2rad($latitude1)) * cos(deg2rad($latitude2)) * sin($dLon/2) * sin($dLon/2);  
  $c = 2 * asin(sqrt($a));  
  $d = $earth_radius * $c;  

  return $d;  
}

$distance = getDistance(56.130366, -106.34677099999, 57.223366, -106.34675644699);
if ($distance < 100) {
  echo "Within 100 kilometer radius";
} else {
  echo "Outside 100 kilometer radius";
}

How to check if coordinate inside certain area Python

from recommendation of @user1753919 in his/her comment, I got the answer here: Haversine Formula in Python (Bearing and Distance between two GPS points)

final code:

from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

center_point = [{'lat': -7.7940023, 'lng': 110.3656535}]
test_point = [{'lat': -7.79457, 'lng': 110.36563}]

lat1 = center_point[0]['lat']
lon1 = center_point[0]['lng']
lat2 = test_point[0]['lat']
lon2 = test_point[0]['lng']

radius = 1.00 # in kilometer

a = haversine(lon1, lat1, lon2, lat2)

print('Distance (km) : ', a)
if a <= radius:
    print('Inside the area')
else:
    print('Outside the area')

Thanks

---

Related Topics