4

I am working with a dataset that is similar to this one:

library(raster)
library(sf)

# Define CRSs
crs.latlon <- "+init=epsg:4326"
crs.sirgas <- "+init=epsg:5880"

# define edges in somewhere within brazil
(spbb <- matrix(c(-53,-44,-25,-20), nrow=2, byrow=F))

# create spatial points in lat/lon
spbb <- SpatialPoints(spbb, CRS(crs.latlon))

# convert to sirgas 2000
(spbb <- spTransform(spbb, CRS(crs.sirgas)))

# make a raster extent object rounded to km:
ext <- extent(5200000, 6046000, 7234000, 7756000)

# create raster
r <- raster(ext, ncol=252, nrow=156, crs=crs.sirgas) # ~3 km resolution
r[] <- runif(ncell(r)) * 10

# create station points
pts <- sampleRandom(r, 500, na.rm=TRUE, sp=TRUE, cells=T) #%>%
#plot(st_geometry(sf::st_as_sfc(pts)), col=sf.colors(4), axes = TRUE, cex=2, pch=19)

#check it
plot(r)
plot(pts, add=T, col="red", pch=19)

For each point in pts, I would like to find the n (i.e. 5) closest surrounding ones. The resulting data frame would be something roughly similar to:

id      nearest     distance
pt1       pt3         14607
pt1       pt204       8540
pt1       pt301       6306
pt1       pt455       4956
pt1       pt337       2145

And so on for each point in the dataset.

What is the fastest way to achieve that?

  • 1
    The knn function in the spatialEco package will do exactly this, with the added ability to expand the problem into multivariate space (ie., add covariates in addition to distance). – Jeffrey Evans Feb 28 at 18:53
3

The geos functions provided in sf, combined with some matrix operations, will do the trick. First, convert and transform your points.

library(sf)

spbb                             # define your points here

spbb <- st_as_sf(spbb)           # convert to simple features
spbb <- st_transform(spbb, 5880) # transform to your desired proj (unit = m)

Then, use st_distance to create a matrix of distances between each point in the data frame. Replace the diagonal with NAs (because distance of a point to itself is 0).

dist_matrix   <- st_distance(spbb, spbb)           # creates a matrix of distances
diag(dist_matrix) <- NA                            # replaces 0s with NA

Lastly, use some matrix operations to get the smallest value of each row.

spbb$distance <- rowMins(dist_matrix)              # get the dist of nearest element
spbb$nearest  <- rowMins(dist_matrix, value = T)   # get the index of the nearest element
| improve this answer | |
  • Cool, I look forward to testing your suggestion! Which package is the function rowMins from? – thiagoveloso Feb 27 at 18:43
1

As I am not completely understanding what the exact output is you want (you want new "features" with every time the point and the five nearest (=500 new features), or do you want just to eliminate the points that are never the nearest, or just a dataframe with indicating the index of the 5 nearest point for every point, or....?), I will tell you how it is possible to calculate the distance between every point:

You can convert your pts to an sf object using st_as_sf(pts) and then calculate the distance between the points using:

dist = as.data.frame(st_distance(pts,pts))

Then you can extract for every point the 5 nearest ones. If you specify what the exact output is you want, I can perhaps give a more precise answer.

| improve this answer | |
0

I ended up writing a code that works as expected. Starting from the example provided in the original post, the following should work:

# convert points to sf objects
pts <- st_as_sf(pts)

# creates a matrix of distances between the points
dist_matrix <- st_distance(pts)

# replaces 0s with NA
diag(dist_matrix) <- NA

# convert matrix to data frame and set column and row names
dist_matrix <- data.frame(dist_matrix)
names(dist_matrix) <- pts$cell
rownames(dist_matrix) <- pts$cell

# find the 5 nearest stations and create new data frame
library(tidyr)
library(dplyr)

near <- dist_matrix %>% 
  mutate(ID=rownames(.)) %>% 
  gather('closest','dist',-ID) %>% 
  filter(!is.na(dist)) %>% 
  group_by(ID) %>% 
  arrange(dist) %>% 
  slice(1:5) %>% 
  mutate(dist_rank=1:5)

Hope this is useful to someone else.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.