I am trying to create a grid over a polygon from which I can choose a random subset of cells for sampling. The problem occurs when I crop the grid, which I am converting to polygons, by the polygon shapefile. For some reason this results in a polygon the shape of the cropping polygon but with lines across it. After opening the resulting shapefile in QGIS I can see that the lines are "edges" of the polygon, because when I get rid of the outline of the polygon, the lines go away.
Can you explain what error, if any, I am making, or do you have an alternate means of generating a grid over a polygon from which a subset of grid cells can be selected (and ideally turned in a shapefile)?
Here is some sample code:
(California boundary found at https://data.ca.gov/dataset/ca-geographic-boundaries)
require(rgdal)
require(sp)
require(raster)
require(rgeos)
# Load shapefile
ca <- readOGR(dsn = "...", layer = "CA_State_TIGER2016")
ca.utm <- spTransform(ca, CRSobj = "+proj=utm +zone=11 ellps=WGS84")
# Calculate number of grid cells necessary
x.cells <- ceiling((ca.utm@bbox[1,2]-ca.utm@bbox[1,1])/10000)
y.cells <- ceiling((ca.utm@bbox[2,2]-ca.utm@bbox[2,1])/10000)
# Generate grid
ca.grid <- raster(extent(ca.utm), nrow=y.cells,
ncol=x.cells, crs=ca.utm@proj4string)
ca.grid.poly <- as(as(ca.grid, "SpatialPixels"),"SpatialPolygons")
ca.grid.utm <- spTransform(ca.grid.poly, CRSobj = "+proj=utm +zone=11 ellps=WGS84")
# Crop grid to polygon
ca.grid.inter <- gIntersection(ca.grid.utm, ca.utm)
# or
ca.grid.crop <- crop(ca.grid.utm, ca.utm)
# both gIntersection and crop give same result as far as I can tell
Here is a screen shot:
plot(ca.grid.inter)
