# calculate distance of points from centroid

I have a spatialpoints object, and a spatialpolygonsdataframe, composed of rectangular cells, and which the spatialpoints fall into. What I would like is to determine which polygon in the spatialpolygonsdataframe each point from the spatialpoints falls into, the distance each point is from the centroid of the respective polygon, and the angle that the point is in relationship to the centroid (with the centroid acting as 0,0).

I know, for example, I could use geosphere::bearing(), raster::pointDistance(), and one of the many centroid functions, but I am not sure how to combine all this into working code.

Here's some sample code. SP is the spatialpoints object, and rintersect3 is the spatialpolygonsdataframe.

owin2Polygons <- function(x, id="1") {
stopifnot(spatstat::is.owin(x))
x <- spatstat::as.polygonal(x)
closering <- function(df) { df[c(seq(nrow(df)), 1), ] }
pieces <- lapply(x\$bdry,
function(p) {
sp::Polygon(coords=closering(cbind(p\$x,p\$y)),
hole=spatstat.utils::is.hole.xypolygon(p))  })
z <- sp::Polygons(pieces, id)
return(z)
}

owin2SP <- function(x) {
stopifnot(spatstat::is.owin(x))
y <- owin2Polygons(x)
z <- sp::SpatialPolygons(list(y))
return(z)
}

test.verts <- list()
test.verts\$x <- c(518, 594, 421, 286, 104, 129, 315)
test.verts\$y <- c(273, 608, 1202, 1201,  620,  285,  135)
test.owin <- spatstat::owin(poly = test.verts)
ppp <- spatstat::ppp(x = c(344, 276, 535, 506, 459), y = c(194, 865, 551, 272, 438), window = test.owin)

SP <- maptools::as.SpatialPoints.ppp(ppp)

w3 <-  spatstat::Window(ppp)

#below converts window, and overlays a grid with 10x10 cells
spatpoly3 <- owin2SP(w3) #internal function

bbow3 <- spatstat::as.owin(spatstat::boundingbox(w3))
spatpoly3 <- sf::st_as_sf(spatpoly3)

#---Alphanumeric Grid:-------
# make grid
ng1 <- sf::st_make_grid(spatpoly3, n = c(5, 5))
ng1 <- sf::st_sf('ID' = seq(length(ng1)), ng1)

# rasterise grid (ID becomes cell value)
nr1 <- fasterize::raster(ng1, ncol = 5, nrow = 5)
nr1[] <- ng1\$ID

# get X and Y from row and col numbers, convert Y to letters

ng1\$X <- raster::colFromCell(nr1, ng1\$ID) # or seq(ncell(nr1))
ng1\$Y <- rev(sapply(raster::rowFromCell(nr1, ng1\$ID), excel_col))

# combine for label
ng1\$LAB <- paste0(ng1\$Y, ng1\$X)

# trim back to polygon
clip2 <- sf::st_join(lwgeom::st_make_valid(ng1), sf::st_as_sf(spatpoly3), join = st_intersects)

geo <- as(clip2, "Spatial")
spatpoly3 <- owin2SP(w3) #internal function
spatpoly3 <- rgeos::gBuffer(spatpoly3, byid=TRUE, width=0)
geo <- rgeos::gBuffer(geo, byid=TRUE, width=0)

CRS.new <- sp::CRS("+init=epsg:27700")

#proj4string(geo) <- CRS.new
sp::proj4string(spatpoly3) <- CRS.new

geo <- sp::spTransform(geo, sp::CRS("+init=epsg:27700"))

#Finished product; leaf with alphanumerically labelled grid
rintersect3 <- raster::intersect(geo, spatpoly3)

sp::proj4string(SP) <- CRS.new

points_over3 <- GISTools::poly.counts(SP, rintersect3)

the following worked out for me:

library(rgdal)
library(geosphere)
library(rgeos)
# crs projection for standardizing:
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

SP <- spTransform(SP, crs.geo)
rintersect3 <- spTransform(rintersect3, crs.geo)
rintersect3\$info <- 1:nrow(rintersect3) # for matching row number of each polygon
# adding data to the SPoint object:
SP <- SpatialPointsDataFrame(coords = SP@coords, data = data.frame(id = 1:nrow(SP)), proj4string = crs.geo)
# the following to be calculated:
SP\$polyg <- as.integer(NA)
SP\$distance <- as.numeric(NA)
SP\$angle <- as.numeric(NA)

for (i in 1:nrow(SP)) {

point <- SP[i,]
# overlaying point and polygon:
aa<-over(point,rintersect3[,'info']) # row number of the polygon over which the point lays.
SP@data[i,'polyg'] <- as.integer(aa)
polyg <- rintersect3[as.integer(aa),] # polygon where the point is located
# calculating distances (in meters) between the point and the polygon's centroid:
dist <- distm (point, gCentroid(polyg), fun = distHaversine)
SP@data[i,'distance'] <- as.numeric(dist)
# calculating angle -- use trigonometric properties of angle sucha that tan(x)^1 = opposite/adjacent
sides <- abs(point@coords-gCentroid(polyg)@coords)
angle<-atan(sides[1]/sides[2])*180/pi
SP@data[i,'angle'] <- angle

}
>>SP@data
id polyg distance    angle
1   1  9054 27230.52 25.95953
2   2  9054 27102.47 28.30076
3   3  9130 25359.85 77.90524
4   4  9338 31300.92 40.28202
5   5  9340 19369.56 70.48203
6   6  9339 28448.07 39.80557
7   7  9340 24921.74 84.69489
8   8  9132 25627.00 45.00000
9   9  9132 25627.00 45.00000
10 10  9131 23645.29 30.93837

references for calculating the angle: here and here. Important: the angle might change depending on your preferred projection.