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

[Roger Bivand]

On Mon, 12 Nov 2012, dani jon wrote:

> Hi Roger,
> I apologize for multiple posting.  I thought the list automatically posts
> the questions and  didn't know the questions go through admin before being
> posted.  I'll be patient next time...
>

No administrative moderation takes place on this list - all messages from 
subscribers are circulated immediately, and shown in the archive. All your 
three copies were sent to over 2600 recipients worldwide. Since you 
yourself have chosen to receive digests (once daily), please wait at least 
24 hours for any reply onlist. Please also do use the archives to check 
whether a message has been posted.

This isn't just a response to you, other list members seem to need 
reminding from time to time, certainly about not posting HTML.

Roger

>
>
>
> On Mon, Nov 12, 2012 at 2:40 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>
>> 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<http://www.r-project.org/mail.html#instructions>
>>
>> and posting guide:
>>
>> http://www.r-project.org/**posting-guide.html<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<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<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
>>
>>
>

-- 
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