How to find "number of neighbors" from a given set of cartesian coordinates

[Dobrokhotov1989]

Hi,

I have a tibble with a set of cartesian coordinates of objects. I need to find a local density around each object. The local density I define as a number of objects reside in the circle around the object of interest. Other words I need to find a number of objects in a red circle minus one (in this example radius = 30) and repeat this for each point on the plot.
Previously I used square instead of circle as approximation but not I need to be more precise. Is there any simple algorithm / function to solve this task?

library(tidyverse)
library(ggforce)

#Sample data
df <- structure(
  list(
    ObjectNumber = 1:83,
    Center_X = c(75.3622047244095, 118.418244406196, 
                              138.2156133829, 151.69918699187, 152.115894039735, 166.151933701657, 
                              170.233890214797, 178.8127090301, 195.203438395415, 197.945378151261, 
                              202.595307917889, 219.089330024814, 228.208888888889, 239.743083003953, 
                              260.227941176471, 4.80110497237569, 4.70967741935484, 4, 10.8529411764706, 
                              12.8484848484848, 14.5706214689266, 10.7, 21.48125, 29.7473684210526, 
                              31.3709677419355, 33.0898876404494, 160.535836177474, 157.873015873016, 
                              74.6802325581395, 100.332764505119, 124.740196078431, 84.4304932735426, 
                              104.28013029316, 128.556451612903, 139.678571428571, 168.125423728814, 
                              168.129370629371, 181.983552631579, 198.579326923077, 223.676975945017, 
                              2.06849315068493, 10.3079584775087, 13.4020100502513, 84.1421800947867, 
                              94.2321428571429, 128.511627906977, 139.585106382979, 167.854237288136, 
                              167.249134948097, 181.57328990228, 198.026378896883, 223.456790123457, 
                              120.709874448592, 155.283625730994, 161.153439153439, 162.259541984733, 
                              184.914285714286, 191.828571428571, 191.511764705882, 189.782805429864, 
                              193.07881773399, 205.176470588235, 204.009411764706, 210.983870967742, 
                              216.93536121673, 219.901098901099, 231.946472019465, 227.904761904762, 
                              234.648910411622, 232.892307692308, 234.834239130435, 235.601286173633, 
                              240.765243902439, 257.485714285714, 259.947692307692, 261.067708333333, 
                              270.232727272727, 273.879518072289, 277.845425867508, 279.330275229358, 
                              285.195599022005, 292.333333333333, 299.894736842105),
    Center_Y = c(3.89763779527559, 
                              22.6006884681583, 61.3122676579926, 11.1517615176152, 85.3973509933775, 
                              43.4861878453039, 70.5298329355609, 7.57859531772575, 77.8080229226361, 
                              27.5546218487395, 11.5923753665689, 23.3002481389578, 289.448888888889, 
                              268.95256916996, 286.632352941176, 203.745856353591, 234.264516129032, 
                              292.7, 185.957219251337, 219.411255411255, 253.189265536723, 
                              269.733333333333, 275.18125, 235.361403508772, 197.322580645161, 
                              223.797752808989, 124.0204778157, 109.320105820106, 54.2093023255814, 
                              27.5546075085324, 10.8480392156863, 243.443946188341, 256.074918566775, 
                              49.7983870967742, 53.1224489795918, 38.3220338983051, 66.9020979020979, 
                              12.3157894736842, 37.7235576923077, 5.7319587628866, 42.7397260273973, 
                              60.3840830449827, 28.7989949748744, 244.530805687204, 258.416666666667, 
                              51.4496124031008, 54.1808510638298, 39.664406779661, 68.3840830449827, 
                              13.6644951140065, 39.0023980815348, 6.41358024691358, 273.777740074652, 
                              154.669590643275, 239.55291005291, 274.834605597964, 287.663492063492, 
                              179.651948051948, 220.174509803922, 262.85520361991, 199.689655172414, 
                              151.244705882353, 241.642352941176, 187.610215053763, 215.041825095057, 
                              258.777472527473, 299.951338199513, 121.095238095238, 170.680387409201, 
                              199.157692307692, 144.823369565217, 262.929260450161, 112.268292682927, 
                              236.848214285714, 286.452307692308, 196.216145833333, 265.141818181818, 
                              172.371485943775, 146.596214511041, 128.651376146789, 213.745721271394, 
                              293.911270983213, 240.784380305603)),
  row.names = c(NA, -83L),
  class = c("tbl_df", "tbl", "data.frame"))

# Sample object 
test_obj <- 7
# Desired radius around object
radius <- 30

ggplot() +
  geom_point(data = df,
             aes(x = Center_X,
                 y = Center_Y)) +
  geom_circle(data = df[test_obj,],
              aes(x0 = Center_X,
                  y0 = Center_Y,
                  r = radius),
              color = "red") +
  annotate(geom = "rect",
           xmin = df[[test_obj, "Center_X"]] - radius,
           xmax = df[[test_obj, "Center_X"]] + radius,
           ymin = df[[test_obj, "Center_Y"]] - radius,
           ymax = df[[test_obj, "Center_Y"]] + radius,
           alpha = 0,
           color = "blue") +
  theme_classic()

# To count number of neighbors in a square around object, I came up with the 
# following approach

count_neighbors <- function(x, y, Center_X, Center_Y, radius) {
  sum(Center_X >= x-radius & Center_X <= x+radius &
        Center_Y >= y-radius & Center_Y <= y+radius) - 1
}

df %>% mutate(density = map2_dbl(Center_X,
                                 Center_Y,
                                 ~count_neighbors(.x,
                                                  .y,
                                                  Center_X,
                                                  Center_Y,
                                                  radius)))
#> # A tibble: 83 x 4
#>    ObjectNumber Center_X Center_Y density
#>           <int>    <dbl>    <dbl>   <dbl>
#>  1            1     75.4     3.90       1
#>  2            2    118.     22.6        4
#>  3            3    138.     61.3       10
#>  4            4    152.     11.2        5
#>  5            5    152.     85.4        5
#>  6            6    166.     43.5        9
#>  7            7    170.     70.5        5
#>  8            8    179.      7.58       5
#>  9            9    195.     77.8        3
#> 10           10    198.     27.6       10
#> # ... with 73 more rows

[rata]

You could translate this into a spatial problem and use the amazing package sf . However, the results from your function and sf's output seems to differ, I'm not sure why.

library(sf)
df_sf <- st_as_sf(df, coords = c("Center_X", "Center_Y"), remove = FALSE)


df_sf <- df_sf %>%
  mutate(
    density = map_int(st_is_within_distance(df_sf, dist = 30), ~length(.x))
  )

df_sf

Simple feature collection with 83 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 2.068493 ymin: 3.897638 xmax: 299.8947 ymax: 299.9513
CRS:           NA
# A tibble: 83 × 5
   ObjectNumber Center_X Center_Y            geometry density
 *        <int>    <dbl>    <dbl>             <POINT>   <int>
 1            1     75.4     3.90  (75.3622 3.897638)       1
 2            2    118.     22.6  (118.4182 22.60069)       4
 3            3    138.     61.3  (138.2156 61.31227)       7
 4            4    152.     11.2  (151.6992 11.15176)       4
 5            5    152.     85.4  (152.1159 85.39735)       6
 6            6    166.     43.5  (166.1519 43.48619)       8
 7            7    170.     70.5  (170.2339 70.52983)       6
 8            8    179.      7.58 (178.8127 7.578595)       6
 9            9    195.     77.8  (195.2034 77.80802)       4
10           10    198.     27.6  (197.9454 27.55462)       8
# … with 73 more rows

ggplot(df_sf) +
  geom_sf(aes(color = density)) +
  scale_color_gradient2(low = "blue",high = "red",midpoint = 5)

[Dobrokhotov1989]

@rata ,
Thank you! Indeed, sf library works well for this task.
The only issue I faced at the beginning is that st_is_within_distance() does not care about multiple images in the same tibble. To solve this, I've just nested data by image number.
The fact that results from sf and my function are different is expected: as square has large area, the number of neighbors fit in this area might be slightly bigger.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(ggforce)

test_data <- structure(list(ObjectNumber = c(1, 12, 15, 16, 19, 23, 24, 26, 
                                             37, 41, 42, 53, 55, 58, 60, 67, 73, 76, 79, 86, 87, 94, 98, 101, 
                                             109, 112, 1, 2, 3, 4, 21, 24, 28, 29, 35, 36, 37, 45, 48, 49, 
                                             50, 52, 58, 63, 66, 67, 68, 70, 77, 78, 85, 86, 87, 93, 95, 98, 
                                             100, 103, 106, 107, 113, 116, 117, 118, 124, 129, 134, 135, 141, 
                                             149, 150, 152, 153, 161, 167, 171, 172, 178, 182, 187, 192, 194, 
                                             195, 205, 208, 212, 213, 215, 225, 226, 227, 228),
                            ImageNumber = c(1, 
                                            1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                            1, 1, 1, 1, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 
                                            231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 
                                            231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 
                                            231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 
                                            231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 231, 
                                            231, 231, 231, 231),
                            Location_Center_X.nuc = c(2.7, 10.78, 16.5, 
                                                      18.55, 18.97, 29.33, 31.96, 37.65, 57.55, 68.33, 68.42, 91.45, 
                                                      97.53, 96.94, 100.7, 110.67, 125.04, 131.77, 144.7, 152.67, 156.71, 
                                                      171.42, 182.47, 184.09, 188.84, 195.69, 3.34, 2.29, 5.03, 7.46, 
                                                      11.12, 11.59, 19.03, 17.41, 19.12, 25.01, 23.44, 32.16, 33.52, 
                                                      37.06, 36.52, 43.12, 45.1, 49.66, 52.97, 57.93, 55.05, 58.03, 
                                                      62.17, 62.33, 69.31, 66.63, 70.4, 77.31, 76.77, 80.33, 86.56, 
                                                      87.18, 90.12, 91.2, 95.23, 97.14, 95.75, 96.49, 104.34, 102.33, 
                                                      109.73, 113.62, 114.56, 121.83, 123.6, 124.52, 130.46, 134.06, 
                                                      137.46, 145.4, 141.91, 151.23, 157.52, 161.62, 166.79, 168.69, 
                                                      168.55, 173.45, 176.47, 180.37, 182.56, 184.23, 194.7, 195.26, 
                                                      195.93, 196.81), 
                            Location_Center_Y.nuc = c(124.2, 9.07, 39.69, 
                                                      88.32, 107.83, 130.32, 16.64, 163.7, 108.46, 79.67, 136.2, 147.17, 
                                                      57.75, 82.03, 14.55, 33.63, 132.93, 41.27, 70.28, 165.24, 94.59, 
                                                      123.95, 169.95, 59.19, 78.74, 122.71, 25.82, 80.05, 128.69, 153.43, 
                                                      2.65, 60.69, 77.48, 114.76, 29.65, 95.9, 173.43, 16.66, 44.9, 
                                                      154.85, 184.81, 116.33, 59.76, 99.47, 19.71, 134.08, 156.6, 189.99, 
                                                      72.41, 110.83, 33.88, 79.92, 164.67, 184.14, 106.82, 143.57, 
                                                      130.66, 79.56, 35.85, 98.25, 55.28, 157.38, 185.17, 119.71, 139.45, 
                                                      12.83, 95.28, 74.37, 122.15, 151.17, 182.18, 102.21, 165.77, 
                                                      127.62, 199.93, 14.21, 103.97, 72.85, 38.29, 112.84, 143.04, 
                                                      81.42, 172.71, 19.45, 199.71, 68.15, 45.52, 8.35, 84.83, 134.14, 
                                                      29.76, 172.72)),
                       row.names = c(NA, -92L),
                       class = c("tbl_df", "tbl", "data.frame"))

radius = 50

df_sf <- test_data %>%
  group_by(ImageNumber) %>%
  nest() %>%
  mutate(data = map(data, ~st_as_sf(.x, coords = c("Location_Center_X.nuc", "Location_Center_Y.nuc"), remove = FALSE)),
         data = map(data, .f = function(x){
           x <- x %>%
             mutate(neighbours = st_is_within_distance(x, dist = radius),
                    density_circle = map_int(neighbours, ~length(.x)-1L)) %>%
             select(-neighbours)
           return(x)
         })) %>%
  unnest(cols = data)

test_obj = 85

ggplot() +
  geom_point(data = df_sf %>% filter(ImageNumber == 231),
             aes(x = Location_Center_X.nuc,
                 y = Location_Center_Y.nuc)) +
  geom_circle(data = df_sf %>% filter(ImageNumber == 231 & ObjectNumber == test_obj),
              aes(x0 = Location_Center_X.nuc,
                  y0 = Location_Center_Y.nuc,
                  r = radius),
              color = "red") +
  theme_classic()

df_sf %>% filter(ImageNumber == 231 & ObjectNumber == test_obj)
#> # A tibble: 1 x 6
#> # Groups:   ImageNumber [1]
#>   ImageNumber ObjectNumber Location_Center_X.~ Location_Center_Y~ density_circle
#>         <dbl>        <dbl>               <dbl>              <dbl>          <int>
#> 1         231           85                69.3               33.9             10
#> # ... with 1 more variable: geometry <POINT>
count_neighbors <- function(x, y, xVal, yVal, radius) {
  sum(xVal >= x-radius & xVal <= x+radius &
        yVal >= y-radius & yVal <= y+radius) - 1
}

df_sf <- df_sf %>%
  mutate(
    density_square = map2_dbl(Location_Center_X.nuc,
                           Location_Center_Y.nuc,
                           ~count_neighbors(.x,
                                            .y,
                                            Location_Center_X.nuc,
                                            Location_Center_Y.nuc,
                                            radius))
  )

df_sf %>%
  filter(ImageNumber == 231) %>%
  ggplot(aes(x = density_circle,
             y = density_square)) +
  geom_point() +
  geom_abline(slope = 1) +
  theme_classic()

Created on 2021-09-13 by the reprex package (v2.0.1)