Outline method

Compute the full distance matrix and then expand it as a dataframe with three columns, {i, j, distance}.

Then merge with the input data on i and j so you have rows with columns:

{i, j, distance, year_i, year_j, territory_i, territory_j}

then applying over groups defined by i (eg use split to break the dataframe into one dataframe for each i) filter for year_i == year_j and territory_i != territory_j. Sort by distance and take the first 1 or first 3.

Alternatively, use the FNN package with the knnx.index function, but you'll probably have to loop over each point, select the other points in the same year and different territory, and use knn.index to find the 1 or 3 nearest neighbours.

I don't think we have a sample dataset sufficient to test at the moment, but these approaches should get you there. As a beginner in R you may need to step back to master some of the basics like dataframe subsetting, using lapply to loop over lists and so on.

Full Solution with Sample Data and Code

First lets make a sample data set with n rows, with random locations, years, and territories:

set.seed(123) library(sf) n=100 df = st_as_sf(data.frame(x=runif(n), y=runif(n), year = sample(2000:2004, n, TRUE), territory = sample(LETTERS[1:12], n, TRUE) ), coords=1:2) 

Let's get the full distance matrix. This is the inefficientest part of this, since it has to compute all NxN distances. Proper nearest neighbour codes can improve by making spatial indexes, but if you're not doing this with a Big Data set or doing it a zillion times with the same points then this should be okay:

dm = st_distance(df) 

At this point:

> class(dm) [1] "matrix" "array" > dim(dm) [1] 100 100 

Next turn that NxN matrix into a data frame with index columns for the rows of df and the corresponding distance:

ijd = data.frame(expand.grid(i=1:n, j=1:n)) ijd$distance = c(dm) 

At this point:

> class(ijd) [1] "data.frame" > dim(ijd) [1] 10000 3 > head(ijd) i j distance 1 1 1 0.0000000 2 2 1 0.5675434 3 3 1 0.1647495 4 4 1 0.6929705 

Now put the year and territory into that data frame by getting the ith and jth rows from df:

ijd$year.i = df$year[ijd$i] ijd$year.j = df$year[ijd$j] ijd$territory.i = df$territory[ijd$i] ijd$territory.j = df$territory[ijd$j] 

At this point:

> head(ijd) i j distance year.i year.j territory.i territory.j 1 1 1 0.0000000 2003 2003 I I 2 2 1 0.5675434 2000 2003 A I 3 3 1 0.1647495 2002 2003 E I 4 4 1 0.6929705 2003 2003 L I 5 5 1 0.6633056 2002 2003 E I 

Filter out differing territories and keep matching years. The "differing territories" criterion ensures a point can't be a nearest neighbour of itself.

ijd = ijd[ijd$year.i == ijd$year.j,] ijd = ijd[ijd$territory.i != ijd$territory.j,]

At this point:

> head(ijd) i j distance year.i year.j territory.i territory.j 4 4 1 0.6929705 2003 2003 L I 7 7 1 0.3958940 2003 2003 E I 8 8 1 0.6049048 2003 2003 G I 9 9 1 0.3247383 2003 2003 D I 

(Note, I've not tested this if there are no points matching the above criteria for a given point)

Split into data frames for each i point.

dg = split(ijd, ijd$i) 

At this point:

> class(dg) [1] "list" > length(dg) [1] 100 > head(dg[[23]]) i j distance year.i year.j territory.i territory.j 1023 23 11 0.7027936 2001 2001 I A 1723 23 18 0.8808717 2001 2001 I C 2923 23 30 0.6241636 2001 2001 I C 3123 23 32 0.6396957 2001 2001 I K 

Now we want to sort those by distance and get the top 3, unless there's less than 3 in which case return the lot. Let's write a function that does that given a data frame so we can test this before trying:

n3 = function(d){ d = d[order(d$distance),] d[1:min(c(nrow(d),3)),] } 

And let's apply it:

dn = lapply(dg,n3) 

At this point:

> length(dn) [1] 100 > class(dn) [1] "list" > dn[[23]] i j distance year.i year.j territory.i territory.j 6423 23 65 0.1805550 2001 2001 I C 6723 23 68 0.1924738 2001 2001 I B 9823 23 99 0.2330713 2001 2001 I B 

This returns a list of data frames like this:

> dn[[1]] i j distance year.i year.j territory.i territory.j 5601 1 57 0.1624860 2003 2003 I B 3701 1 38 0.1994243 2003 2003 I G 7301 1 74 0.3222224 2003 2003 I J > dn[[2]] i j distance year.i year.j territory.i territory.j 8602 2 87 0.2080085 2000 2000 A J 6002 2 61 0.2095182 2000 2000 A C 9302 2 94 0.3158656 2000 2000 A C 

The first one showing that point 1's nearest neighbours matching the criteria are 57, 38 and 74, point 2's nearest neighbours are 87, 61 and 94. The list has 100 elements (unless a point has no neighbours under the matching criteria, because I think they will just drop out, untested...). Note the years match and the territories dont.

If you want to drop all the other columns and make a simple data frame then:

nnij = do.call(rbind, dn)[,c("i","j")] 

giving:

> head(nnij, 10) i j 1.5601 1 57 1.3701 1 38 1.7301 1 74 2.8602 2 87 2.6002 2 61 2.9302 2 94 3.1803 3 19 3.6203 3 63 3.1103 3 12 4.8804 4 89 >