Search code examples
rr-rastercentroid

Delete cells of grid whose centroid is not over a raster

I have two problems.

Problem 1: I would like to delete all cells of the grid whose centroid is not located over the raster. I'm not even sure if I'm working with the right "types of objects" (RasterLayer, SpatialPixels, etc.).

See example with dummy data below:

# Load package 
library(raster)

# Create raster and define coordinate reference system
ras <- raster(nrows = 100, ncol = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)
proj4string(ras) <- CRS("+init=epsg:32198")

# Generate random values
val <- sample(x = 1:100, size = ncell(ras), replace = T)
values(ras) <- val

# Create effort grid
xym <- matrix(c(-30,130,130,-30,-30,-30,130,130), nrow = 4, ncol = 2) 
p <- Polygon(xym)
ps <- Polygons(list(p), 1)
sps <- SpatialPolygons(list(ps))
proj4string(sps) <- CRS("+init=epsg:32198")
data <- data.frame(f = 99.9)
spdf <- SpatialPolygonsDataFrame(sps, data)
ptsreg <- spsample(spdf, 50, type = "regular")
grid <- SpatialPixels(ptsreg)

# Plot raster over grid
plot(grid)
plot(ras, add = T)

Problem 2: Is there another way to create an effort grid? My code works but I'm pretty sure there's a simpler way?

Also... In this example, I plotted the grid first and added the raster. If I do it the other way (raster first and grid second), the resulting plot doesn't show the whole extent of the grid.

enter image description here

How can I plot the grid over the raster but still show the entire grid? Such as: enter image description here

Solution

  • For your first question, here is an option. Try:

    cont <- as(extent(ras), "SpatialPolygons") 
proj4string(cont ) <- CRS("+init=epsg:32198")
grid.small <- grid[!is.na(grid%over%cont)]

    You first make and polygon with the extent of your raster, then you subsample your grid with the %over% function. I think it gives you the result you are looking for.

    For your second question, it isn't the best, but try:

    plot(grid)
plot(ras,add=T)
plot(grid,add=T)

    Simple and works. Normally, 

    plot(ras,ext=extent(grid))
plot(grid,add=T)

    Should work but it's not and I couldn't figure out why...

    As for everything in R, there is always multiple ways to do one thing. If your try what I propose and it's not efficient with your real data, then it could worth trying to find a other solution.