2

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)
1

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.

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.