[R-sig-Geo] Overlay, sample, and extract

[Roger Bivand]

Patience is a virtue - this is your third identical posting on this topic. 
Please do respect the time list members give. Please also read the list 
instructions:

http://www.r-project.org/mail.html#instructions

and posting guide:

http://www.r-project.org/posting-guide.html

and note that you should only post plain text, not HTML, content. If you 
look at the archived version of your posting, you see:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016684.html

which treats your posting as an attachment. We strip HTML to reduce carbon 
footprint of the list (fewer bytes to store and transmit), and to reduce 
danger from included content. Please play by the rules ...

A slightly tired list admin.

On Mon, 12 Nov 2012, dani jon wrote:

> Dear List,
> I am interested in overlaying a regular polygon grid of 8 by 8 into a
> spatial surface. Essentially, dividing the spatial surface into 64 polygon
> grids and would like to take a sample of 10 points (pixels) with
> corresponding xy coordinate and values for each of the 64 grids.  Here is
> my attempt.
>
> library(maptools)
> library(raster)
> library(sp)
> library(gstat)
> library(maps)
>
> ## Creating spatial surface
>
> xg = seq(1, 100,, 100)
> yg = seq(1, 100,, 100)
> xy = expand.grid(xg, yg)
> names(xy) <- c("x","y")
> g.dummy <- gstat(formula=z~1, locations=~x+y, dummy=T, beta=1,
> model=vgm(psill=0.025,model="
> Exp",range=15), nmax=20)
> yy <- predict(g.dummy, newdata=xy, nsim=1)
> gridded(yy) = ~x+y
> coor.yy <- coordinates(yy[1])
> ytrue <-yy[1]@data
> yyspdf<- SpatialPointsDataFrame(coords=coor.yy, data=ytrue)
> yyspdf[1:10,]
> coordinates     sim1
> 1       (1, 1) 1.305260
> 2       (2, 1) 1.287697
> 3       (3, 1) 1.219094
> 4       (4, 1) 1.333913
> 5       (5, 1) 1.268508
> 6       (6, 1) 1.263417
> 7       (7, 1) 1.243237
> 8       (8, 1) 1.276870
> 9       (9, 1) 1.231528
> 10     (10, 1) 1.162322
>
> ## creating regular polygon grid (8x8)
>
> grd <- GridTopology(c(6.3,6.3), c(12.5,12.5), c(8,8))
> SpP_grd <-as.SpatialPolygons.GridTopology(grd)
> plot(SpP_grd)
> text(coordinates(SpP_grd), sapply(slot(SpP_grd, "polygons"), function(i)
> slot(i, "ID")), cex=0.5)
> ids <- sapply(slot(SpP_grd, 'polygons'), function(i) slot(i, 'ID'))
> trdata <- data.frame(A=rep(c(1,2,3,4,5,6,7,8), 8),
> B=rep(c(1,2,3,4,5,6,7,8), each=8),
> row.names=sapply(slot(SpP_grd, "polygons"), function(i) slot(i, "ID")))
> SpPDF <- SpatialPolygonsDataFrame(SpP_grd, trdata)
> SpPDF at data = cbind(SpPDF at data,ids)
>
>> SpPDF at data[1:10,]
>    A B ids
> g1  1 1  g1
> g2  2 1  g2
> g3  3 1  g3
> g4  4 1  g4
> g5  5 1  g5
> g6  6 1  g6
> g7  7 1  g7
> g8  8 1  g8
> g9  1 2  g9
> g10 2 2 g10
>
> ## Overlay
> image(yy[1])
> plot(SpPDF, border="magenta", add=TRUE)
> text(coordinates(SpP_grd), sapply(slot(SpP_grd, "polygons"), function(i)
> slot(i, "ID")), cex=0.5)
>
> My question is how to overlay yy[1] and SpPDF and take 10 samples within
> each polygon grid (g1-g64) and extract the values and coordinates of
> yy[1]?  So I will have an output with a sample of 640
>
> x  y   sim1        ids
> 1  2   1.305260   g1
> 1  3   1.02567     g1
> .
> .
> 95 95 1.3487    g64
>
> Edzer  suggested something likes this, but it does not exactly  640 samples
> (10*64) as I wanted.
>
> x = do.call(rbind, lapply(1:10, function(x) spsample(yy, cellsize = 12.5,
> type = "stratified")))
>
> Thanks in advance
> Best regards
> Sam
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no