[R-sig-Geo] PBSmapping to SpatialPolygonsDF
Virgilio Gómez-Rubio
v.gomezrubio at imperial.ac.uk
Mon Mar 26 22:53:25 CEST 2007
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?
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
> >
> >
>
More information about the R-sig-Geo
mailing list