I have a hard time trying to describe the shapefile exactly, but I am trying to generate a grid over the Gulf of St. Lawrence, such that only the gridsquares which lie over water are included. I would like this to extend to the narrow southern strait at (-60, 47)

I'm really new to GIS in R, and below is my attempt at brute forcing it by generating rows of gridsquares and specifying the bounding box for each line individually. I used st_make_grid from the sf library to make the 5 rows. Each gridsquare is 0.25x0.25 degrees.

enter image description here

I feel like there must be a more elegant solution to this... if anyone knows one.

2 Answers 2


Maybe this isn't an elegant solution, but works:

library(maptools) # for example data


r <- raster() # sample raster
res(r) <- c(5,5) # example resoution (I used a coarse one for speed up computation)

r[] <- 1:ncell(r) # this could be any value 

usa <- wrld_simpl[wrld_simpl$ISO3 == 'USA',] # example country

re <- raster::extract(r,usa, weights = T, normalizeWeights=F, cellnumber = T) # here you extract cell number intersected by usa, I put weights to filter

sup <- re[[1]][,'cell'][which(re[[1]][,'weight'] == 1)] # identify only cells with 100% land coverage, delteing [which(re[[1]][,'weight'] == 1)] you'll have sea-only pixels

r[sup] <- NA # remvoe values for this cells

poly <- as(r,'SpatialPolygonsDataFrame') # to grid

plot(usa);plot(poly,border='red',add=T) # plot

enter image description here


It's hard to define clearly what inside and outside means, because you probably don't want lakes or rivers - but at any rate creating the polygons from a raster mask give some control.

Here I use sf version of Canada/USA states and produce all the polygons, then filter.

x <- rbind(rnaturalearth::ne_states( "United States of America", returnclass = "sf"), 
       rnaturalearth::ne_states( "Canada", returnclass = "sf"))
## grid template on a specific extent
r <- raster(extent(-75, -50, 42, 53), res = 0.25, crs = "+init=epsg:4326")

## fasterize
g <- fasterize(x, r)

## fast conversion pixel to polygons
p <- spex::polygonize(!is.na(g))
## layer is whether we were in or out of land
plot(subset(p, !layer)$geometry)
plot(x$geometry, add = TRUE)

enter image description here


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.