I had a similar issue recently (150 point_As and 700,000 point_Bs) solved using the
point in polygon concept and the
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
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:
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),
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.
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"))
%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"=.)
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.
leaflet(Point_A_buffer) %>% addTiles %>% addMarkers(data=point_A) %>%
addCircleMarkers(data=point_B,radius = 2,fillOpacity=1) %>%
addPolygons(fillColor = "red",weight=1,