[R-sig-Geo] PBSmapping to SpatialPolygonsDF

Roger Bivand Roger.Bivand at nhh.no
Tue Mar 27 08:13:35 CEST 2007


On Mon, 26 Mar 2007, Virgilio Gómez-Rubio wrote:

> Hi,
> 
> I am also interested in this convertion function. In fact, yesterday
> night I wrote some code to convert from PolySet to SP, but Roger's code
> is much better than mine... 
> 
> If nobody is already doing it, can I take Roger's code, add a few checks
> and make a function to be included in, say, maptools? 
> 

Yes, please do - I'm sure more users are interested too.

Roger

> Best regards,
> 
> Virgilio
> 
> El lun, 26-03-2007 a las 16:21 -0400, Andrew Niccolai escribió:
> > Roger!!! Thank you very much.  I needed to adjust your code only to account
> > for the fact that my tesselated object did not have a SID column and also I
> > needed to add a step to "repeat" the first set of XY cords for each polygon
> > to close the loop.  Other than that, it worked like a charm out of the gate.
> > 
> > Thank you.
> > 
> > Here are my changes to Roger's code in case others have data sets w/o closed
> > loops.
> > 
> > ##--- convert from Polyset in PBSmapping to SPDF
> > ## example(PolySet2SpatialPolygons ---written by R.Bivand)
> > res0 <- split(cb.tess, cb.tess$PID)
> > for(i in 1:length(res0)){
> > res0[[i]] <- rbind(res0[[i]],res0[[i]][1,])
> > }
> > outPolygons <- vector(mode="list", length=length(res0)) 
> > for (i in seq(along=outPolygons)) {
> > outPolygons[[i]] <- Polygons(lapply(res0[[i]], function(x)
> > Polygon(cbind(res0[[i]]$X, res0[[i]]$Y))), ID=as.character(i))
> > } 
> > outSP <- SpatialPolygons(outPolygons)
> > plot(outSP)
> > 
> > Andrew Niccolai
> > Doctoral Candidate
> > Yale School of Forestry
> > 203-432-5144
> >  
> > 
> > -----Original Message-----
> > From: Roger Bivand [mailto:Roger.Bivand at nhh.no] 
> > Sent: Monday, March 26, 2007 3:34 PM
> > To: Andrew Niccolai
> > Cc: r-sig-geo at stat.math.ethz.ch
> > Subject: Re: PBSmapping to SpatialPolygonsDF
> > 
> > On Mon, 26 Mar 2007, Andrew Niccolai wrote:
> > 
> > > Is there a pre-existing function or a clever way to pull a PBSmapping
> > object
> > > of class "Polyset", "data.frame" into a SpatialPolygonsDataFrame?
> > > 
> > > I found the SpatialPolygons2PolySet function and was hoping that there is
> > a
> > > way to go the other direction.
> > > 
> > > I have created a voronoi tessellation for a set of 627 points and then
> > > clipped that by the convex hull of those points.  The object looks like
> > > this:
> > > 
> > > > class(cb.tess)
> > > [1] "PolySet"    "data.frame"
> > > > names(cb.tess)
> > > [1] "PID" "POS" "X"   "Y"  
> > > > head(cb.tess)
> > >   PID POS        X        Y
> > > 1   1   1 945062.2 860290.8
> > > 2   1   2 945073.5 860273.0
> > > 3   1   3 945062.7 860272.6
> > > 4   1   4 945062.5 860272.7
> > > 5   1   5 945055.9 860277.4
> > > 6   1   6 945057.8 860291.4
> > > > tail(cb.tess)
> > >   PID POS        X        Y
> > > 1 626   9 944972.3 859813.1
> > > 2 627   1 945009.7 859794.3
> > > 3 627   2 944993.8 859769.0
> > > 4 627   3 944985.3 859763.3
> > > 5 627   4 944990.9 859804.2
> > > 6 627   5 945008.2 859804.6
> > > 
> > > I would like to slip back into the world of SPDF if possible.
> > > 
> > > Thanks in advance for help, suggestions or code.
> > 
> > Not yet a function, but if you could try this and polish it a bit, I think 
> > it's workable:
> > 
> > example(SpatialPolygons2PolySet)
> > # gives the input PolySet nor_coast_poly_PS
> > res0 <- split(nor_coast_poly_PS, nor_coast_poly_PS$PID)
> > res1 <- lapply(res0, function(x) split(x, x$SID))
> > outPolygons <- vector(mode="list", length=length(res0))
> > for (i in seq(along=outPolygons)) {
> >    outPolygons[[i]] <- Polygons(lapply(res1[[i]], function(x) 
> >    Polygon(cbind(x$X, x$Y))), ID=as.character(i))
> > }
> > outSP <- SpatialPolygons(outPolygons)
> > plot(outSP)
> > 
> > # get the projection/zone from attr(nor_coast_poly_PS, "projection") and
> > # attr(nor_coast_poly_PS, "zone") if UTM, use in constructing
> > # the SpatialPolygons object. Now the ID= in Polygons() is just the loop 
> > # index - should be the PID values if they are the key value for adding
> > # in a data frame.
> > 
> > Hope this helps,
> > 
> > Roger
> > 
> > > 
> > > cheers
> > > 
> > > Andrew Niccolai
> > > Doctoral Candidate
> > > Yale School of Forestry
> > > 203-432-5144
> > > 
> > > 
> > 
> 

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




More information about the R-sig-Geo mailing list