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

  • P.S. if anyone knows how to grid the ggmap image. That would also be helpful. – I Del Toro Sep 13 '14 at 12:39
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)
  • 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 8: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 8: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 9:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.