0

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 Raster with associated points and polygons

Here is the ideal outcome

A count of pixels with points and total raster pixels present broken down by polygon ID

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)

# add data
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)
plot(points, pch=1, add=TRUE)
plot(polys, lwd=3, add=TRUE)

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

1 Answer 1

1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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