# Using R to combine information in raster and point vector according to separate polygon vector associated with each

This builds off a previous question I had about counting the pixels in a raster that contains at least 1 point (Using R to count number of pixels in raster file for which there are at least one point object present) which was perfectly answered by @dbaston.

I guess my real problem is a bit more complex and I am curious how to arrive at a similar solution with an added layer of complexity, where instead of just counting up the pixels with points for the entire raster, the counts of pixels with points is made separately for a number of uniquely ID'ed polygons.

Here is a visual representation of the problem

Here is the ideal outcome

Reproducible Example

``````library(raster)
set.seed(123)
#===========================
# POINTS
#===========================
# create points
points <- structure(list(longitude = runif(10), latitude =runif(10)), .Names = c("longitude", "latitude"), class = "data.frame", row.names = c(NA, -10L))

# set coordinates for spatial data frame
coordinates(points) <- cbind(points\$longitude , points\$latitude)

# assign a CRS to spatial object
proj4string(points) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#===========================
# RASTER WITH DATA
#===========================
# create raster
rast <- raster(xmn=0, xmx=1, ymn=0, ymx=1, res=0.25)

rast\$data <- rnorm(length(rast), mean = 10, sd=1)

#===========================
# POLYGONS
#===========================
square <- rbind(c(1, 0.5, 0.5, 0.5, 0.5, 1,
1, 1, 1, 0.5, 1, 0.5),
c(0.5, 0, 0, 0, 0, 0.5,
0.5, 0.5, 0.5, 0, 0.5, 0))

ID <- c("Poly1", "Poly2")

polys <- SpatialPolygons(list(
Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID[1]),
Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID[2])
))

#===========================
# VISUALIZE
#===========================
plot(rast)

``````

That is where I am stuck. I am trying to learn how to work with spatial data in R.

The approach that worked for me was very straightforward once I figured it out.

Approach

(1) Use the rasterize function in package raster to convert the points to pixels that match the resolution of the original raster. If only presence/absence is of interest, then the 'field' element in the rasterize function can be set to 1; this returns either a 1 for pixels that have one or more points present or an NA for pixels with an absence of points.

``````points_to_raster <- rasterize(x=points, y=rast, field=1)
``````

(2) Use the extract function in package raster polys vector to extract presence/absence data from the rasterized point layer and store as a dataframe containing a record for each pixel with the associated poly ID and a 1 if a point was present.

``````df.points_to_raster <- raster::extract(points_to_raster, polys, df=TRUE, weights=TRUE)
``````

(3) Use dplyr to summarize the data to match the desired output.

``````# Load library
library(dplyr)

# Convert NA's to 0's for counting the total number of pixels
df.points_to_raster[is.na(df.points_to_raster)] <- 0

# Group by poly ID, add up the number of 1's to get the count of pixels with points in each poly, and count the total pixels present in each poly
df.points_to_raster %>%
group_by(ID)%>%
dplyr::summarize(`Count of pixels with points` = sum(layer), `Total pixels` = length(layer))
``````