Spatial Join in R

Question

I'm fairly new to using R for spatial work so bear with me. I have found questions similar to this but nothing that answers this question specifically enough for me to implement it.

I have point data in a spatialPointDataFrame and polygon data in a SpatialPolygons class. I am trying to do a spatial join where the polygons are given the sum of one of the attributes of all the points that fall within the polygon's boundary. This seemed simple enough, but I can't find documentation on how to do it.

I have tried using sp:over for this and it doesn't seem to work.

I have referenced and it works with one polygon set but not the other. Basically the difference is one polygon set is counties and one is ZCTA's. Other than that they should be the same. But, the code works fine for counties but gives

Error in [[<-.data.frame(*tmp*`, name, value = c(122221L, 1465574L,  : 
  replacement has 263 rows, data has 408 

anytime I try to use ZCTA's or a census tracts polygon set. I thought that it may be broken geometry problem so I ran all polygons sets through Repair geometry in arcmap to hopefully eliminate that potential problem.

srdfcounties = SpatialPolygonsDataFrame( centTrue, 
                   data = data.frame(row.names=paste("0", (1:length(centTrue))),
                                                           PIDS=1:length(centTrue)))

srdfZCTA =  SpatialPolygonsDataFrame(centFalse,
                   data = data.frame(row.names=paste("0", (1:length(centFalse))),
                                                           PIDS=1:length(centFalse)))

tail(srdfZCTA@data)
head(srdfcounties@data)
tail(busiPTS@data)
plot(srdfcounties)

points(busiPTS, pch=20)

ptsZCTA.poly <- point.in.poly(busiPTS, srdfZCTA)
ptsCounty.poly <- point.in.poly(busiPTS, srdfcounties)

srdfcounties$sales_vol <- tapply(ptsCounty.poly@data$SALES_VOL,      
                                 ptsCounty.poly@data$PIDS, FUN=sum)
spplot(srdfcounties, zcol='sales_vol')

srdfZCTA$sales_vol <- tapply(pts.poly@data$SALES_VOL, pts.poly@data$PIDS, FUN=sum)
spplot(srdfZCTA, zcol='sales_vol')

EDIT2:

I figured it out. Just to provide closure and insight for anyone else who may have this same problem I will post the working code soon.

Answer 1

This does not seem to be an issue with the function. From the R code, included in your zip file, it looks like there was some confusion in your naming convention that caused a syntax error by calling the wrong objects in tapply. You were also attempting the add the tapply results to the points and not to the source polygon object, which the values were aggregated to

That said, methodologically it would be much safer to create a data.frame from the tapply results and then use merge to join back to the polygon object rather than trying to pipe the results directly to an object. If you think about it, if not all polygons have points that occur in them, then the resulting aggregation will not have the same dimension as the source polygons. Therefore, you cannot just pipe the shorter vector to the polygons and have them match.

Based on your naming conventions and code here is what I got to work.

Point in polygon using "busiPTS" (points) and "srdfZCTA" (tessellated polygons). This results in a SpatialPointsDataFrame object with the polygon IDS stored as the PIDS column.

ptsZCTA.poly <- point.in.poly(busiPTS, srdfZCTA)
  plot(srdfZCTA)
    plot(ptsZCTA.poly, pch=20, cex = 0.5, add=TRUE)

Here we run tapply to aggregate the sum of SALES_VOL to each polygon in "srdfZCTA". This results are an array class object representing the sum for each polygon with intersecting points. We can coerce the array into a data.frame and use names() to include the PIDS values (loosing the names (PIDS) is why we do not directly coerce tapply to a data.frame while calling the function).

zcta.sales.vol <- tapply(ptsZCTA.poly@data$SALES_VOL, 
                     ptsZCTA.poly@data$PIDS, FUN=sum) 
zcta.sales.vol <- data.frame(PIDS=names(zcta.sales.vol), 
                    sales.vol=as.vector(zcta.sales.vol))

Now that we have a data.frame representing sums of SALES_VOL for each polygon that has intersecting points we can use merge to join the results back to the tessellated polygons "srdfZCTA". Polygons that do not have intersecting points will be assigned a NA.

srdfZCTA <- merge(srdfZCTA, zcta.sales.vol, by="PIDS")
  str(srdfZCTA@data)                    
  spplot(srdfZCTA, zcol='sales.vol')   

Because they should be the same as the number of values zcta.sales.vol, we can verify results by looking at the number of non-NA values contained in "srdfZCTA". You will see that of the 408 polygons in "srdfZCTA", 254 of them contain values in the sales.vol column. That matches the number of aggregated values in the "zcta.sales.vol" object created by tapply.

dim(srdfZCTA)[1]
dim(zcta.sales.vol)[1]
dim( srdfZCTA@data[!is.na(srdfZCTA@data$sales.vol),])[1]

I did notice that in your code you had a line of code busiPTS[!is.na(busiPTS@data$SALES_VOL),] that, if there were NA values, would have caused you considerable trouble. When removing NA values from a spatial object you cannot just operate on the @data slot. The removed values will not propagate through the entire sp object leaving you with a mismatch between the number of observations in @data and the rest of the object. I provide a function "sp.na.omit" in spatialEco the deals with this issue.

Answer 2

plot(A)
points(B)
# Overlay points and extract just the code column: 
a.data <- over(A, B[,"code"])
# Add that data back to A:
A$bcode <- a.data$code

From here