1

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

Community
  • 1
I Del Toro
  • 757
  • 1
  • 10
  • 20

1 Answers1

1

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)
Robert Hijmans
  • 10,683
  • 25
  • 35
  • Thanks! Do you know how to get the size of each cell? res= 0.01229634. Is that a decimal degree? How can I convert that to square meters? – I Del Toro Sep 14 '14 at 08:03
  • 1
    If you want everything in meters then use spTransform to reproject to a metric coordinate system (commonly a UTM zone or a Danish national grid system) and then work exclusively with that until you have to flip it back to lat-long. – Spacedman Sep 14 '14 at 08:20
  • @Spacedman Thanks for this suggestion but I'm not sure if spTransfrom would work on the final raster to convert the grid cells to UTMs: UTM_Map <- spTransform(r,CRS(" +proj=utm +zone=32 +ellps=WGS84+datum=WGS84 +units=m +no_defs +towgs84=0,0,0")) – I Del Toro Sep 15 '14 at 09:54