[GIS] Adding a grid and selecting random cells from spatial polygon in r

ggmaprrasterrgdalsp

I have a fairly similar question to this: Similar Question
I am working with a map of Denmark or more a subset of a Denmark Map.

This file comes from the Global administrative areas site GADM site

require (raster)
require (rgeos)
require (rgdal)

DK1 <- getData('GADM', country='DNK', level=0)

#create an extent for clipping in raster package 
ext=  extent (8.41, 9.04, 56.92, 57.16)   
clipe <- as(ext, "SpatialPolygons") 
proj4string(clipe) <- CRS(proj4string(DK1)) 
cropDK <- SpatialPolygonsDataFrame(clipe, data.frame(x = 1), match.ID = FALSE) 

#clip the shapefile using rgeos package
DK_reg1<- gIntersection(DK1, cropDK) 
plot (DK_reg1)

enter image description here
I would like to create a grid of 1000 equally sized cells, plot it and then select and highlight 100 random cells. I seem to be having some trouble because of the projection and undefined units. Does anyone have any tips on how to fix these issues? 

#Grid the clipped shapefile to 1000 equal size units. 
proj4string (DK_reg1)

bb <- bbox(DK_reg1)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration); 1 m = 3.28084 ft
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd

sp_grd <- SpatialGridDataFrame(grd,data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(DK_reg1)))

Finnally is it possible to add this gridded shape file to a ggmap raster showing the main roads? 

library (ggmap)
library(mapproj)
map <- get_map(location=c(lon=8.65, lat=56.955665), zoom =10)
ggmap(map, extent = "normal")
  plot (DK_reg1, add=T)

enter image description here

Best Answer

Here are some ideas

library(raster)
library(rgeos) 

DK1 <- getData('GADM', country='DNK', level=0)

# extent coordinates 
e <- c(8.41, 9.04, 56.92, 57.16)   
DK_reg1 <- crop(DK1, extent(e))

# resolution to get about 1000 cells
res <- sqrt((e[2]-e[1]) * (e[4]-e[3])) / sqrt(1000)

r <- raster(extent(e), res=res, crs=crs(DK1))

# rasterize if you only want to extract values within the area of the polygons
r <- rasterize(DK_reg1, r)

s <- sampleRandom(r, 100, cells=TRUE, xy=TRUE)

# display 
r[s[,1]] <- 2
plot(r)


# add to Google map
library (dismo)
g <- gmap(r)
plot(g, interpolate=TRUE)
points(Mercator(s[, c('x', 'y')]), col='red', pch=20)
Related Question