[R-sig-Geo] Creating another grid from n*n grid

[Roger Bivand]

On Sat, 12 Jan 2008, Takatsugu Kobayashi wrote:

> Hi,
>
> (As an exercise)I am trying to create a SpatialPolygons object from n*n
> grid. By reading the material (from vignette{'sp'), I managed to convert
> grid to polygons. Now I would like plot this polygons, and here is what
> I have done so far:
>
> Step1: Create a 9*9 grid
> n<-9 # number of cells in each axis
> x <- 1:n+1
> y <- 1:n+1
> r.grid <- makeGrid(x,y)

Trying to reproduce your case, there is a problem:

> library(sp)
> n<-9
> x <- (1:n)+1
> y <- (1:n)+1
> r.grid <- makeGrid(x,y)
Error: could not find function "makeGrid"

(It seems to be in PBSmapping).
Doing things directly might be easier:

library(sp)
grd <- GridTopology(c(2,2), c(1,1), c(9,9))
SP <- as(grd, "SpatialPolygons")
plot(SP, axes=TRUE)

See ?GridTopology and ?as.SpatialPolygons.GridTopology, the latter is the 
first line in help(package=sp) - though other names in that listing are 
less intuitive.

Looking at:

getMethod("coerce", c("GridTopology", "SpatialPolygons"))

sp:::as.SpatialPolygons.GridTopology

shows what is going on, and:

sp:::as.SpatialPolygons.PolygonsList

and

SpatialPolygons

the actual code building the SpatialPolygons object. In your exercise, you 
could say debug(SpatialPolygons), which I think would show that the extra 
list() is unnecessary. See if you can get your exercise to give the same 
results as coercing from a GridTopology object to a SpatialPolygons 
object, this is a very typical and effective way to find out how things 
work.

There was correspondence on the list last year with sample code for 
coercing a PolySet object to SpatialPolygons - would it be useful to put 
such a function into the maptools package?

Roger

>
> Step2: Convert the grid to SpatialPolygon
> Srs <- list()
> for (i in 1:n^2)
> {
> x.coord <- c(r.grid$X[(1+4*(i-1)):(1+4*(i-1)+3)],r.grid$X[(1+4*(i-1))])
> y.coord <- c(r.grid$Y[(1+4*(i-1)):(1+4*(i-1)+3)],r.grid$Y[(1+4*(i-1))])
> subPoly <- Polygon(cbind(x.coord, y.coord),hole=FALSE)
> Srs[[i]] <- Polygons(list(subPoly), paste("SP",i,sep=''))
> }
>
> SpP <- SpatialPolygons(list(Srs),1:n^2)
> Error in is.vector(X) : cannot get a slot ("Polygons") from an object of
> type "list"
>
> I am trying to what the following message tells me. class(Srs) returns
> "list" as class and Src has an object of class "Polygons" as follows:
>
> [[81]]
> An object of class ?Polygons?
> Slot "Polygons":
> [[1]]
> An object of class ?Polygon?
> Slot "labpt":
> [1] 9.5 9.5
>
> Slot "area":
> [1] 1
>
> Slot "hole":
> [1] FALSE
>
> Slot "ringDir":
> [1] 1
>
> Slot "coords":
> x.coord y.coord
> [1,] 9 9
> [2,] 9 10
> [3,] 10 10
> [4,] 10 9
> [5,] 9 9
>
> Slot "plotOrder":
> [1] 1
>
> Slot "labpt":
> [1] 9.5 9.5
>
> Slot "ID":
> [1] "SP81"
>
> Slot "area":
> [1] 1
>
> I am sure I am doing fundamental mistakes here. As I am thinking to
> increase n later, I wouldn't want to type Srs[[1]..... by hand.
>
> Also, I am going to convert this polygons into shapefile to get an
> polygon-adjacency matrix through GeoDA. I have reading a few threads
> about a polygon adjacency matrix in R, but it seems GeoDA is a way to
> go? For my case, I would assign 1 if 2 polygons share multiple
> coordinates and store neighbors in list if this is a good way to do in R.
>
> Sorry for this very fundamental question, but I appreciate help.
>
> Thank you.
>
> Taka
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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