[R-sig-Geo] convert a PBSmapping grid to shapefile

[Roger Bivand]

On Fri, 18 Feb 2011, Armelle ROUYER wrote:

> Thanks for your answer Roger but it still doesn't work...
> The problem is with mo, has you said it needs more study and I don't how to 
> solve the problem.
>
> Here is the problem :
>
> library(PBSmapping)
>
> PS0=makeGrid(x=seq(-5,-1,.2),y=seq(43,49,.2),projection="LL")
> PD=read.table("pdata.txt",header=T) # text file attach
>
> brks=c(fivenum(PD$Z))
> cols=c("white","yellow","orange","red","darkred")
> PD0=makeProps(PD,breaks=brks,"col",cols)
>
> plotPolys(PS0, polyProps=as.PolyData(PD0), projection="LL",bg="lightblue")
> sp1 <- PolySet2SpatialPolygons(PS0)
> mo <- match(as.character(PD0$PID), row.names(sp1))
> sp2 <- SpatialPolygonsDataFrame(sp1, data=data.frame(col=PD0$col[mo],
>  stringsAsFactors=FALSE))
> plot(sp2, col=sp2$col, axes=TRUE)
>
>
> Many thanks in advance if you find a solution.

The underlying problem is that makeGrid() understands PolySet objects in 
its own way, with the interaction of PID and SID being a top level 
identifier. In ?PolySet, it looks as though the PID is the top level 
identifier, and SIDs within a given PID are its component parts, not 
independent top-level objects. This is how PolySet2SpatialPolygons() sees 
things too, so the multiple SIDs in a single PID are treated as member 
component polygons of a multipolygon (Polygons) object. So we need to 
disambiguate:

library(PBSmapping)
PS0=makeGrid(x=seq(-5,-1,.2),y=seq(43,49,.2),projection="LL")
PS1 <- as.PolySet(data.frame(PID=as.integer((PS0$PID*100)+PS0$SID),
   SID=1L, POS=PS0$POS, X=PS0$X, Y=PS0$Y), projection="LL")
# Be careful to ensure uniqueness in the new PID!!
PD=read.table("pdata.txt",header=T)
# text file attach
PD1 <- data.frame(PID=as.integer((PD$PID*100)+PD$SID), Z=PD$Z)
# Be careful to ensure uniqueness in the new PID!!
brks=c(fivenum(PD$Z))
cols=c("white","yellow","orange","red","darkred")
PD0=makeProps(PD1,breaks=brks,"col",cols)
plotPolys(PS1, polyProps=as.PolyData(PD0), projection="LL",bg="lightblue")
library(maptools)
sp1 <- PolySet2SpatialPolygons(PS1)
mo <- match(as.character(PD0$PID), row.names(sp1))
Z <- rep(as.numeric(NA), length(row.names(sp1)))
Z[mo] <- PD0$Z
col <- rep(as.character(NA), length(row.names(sp1)))
col[mo] <- PD0$col
df <- data.frame(Z=Z, col=col, stringsAsFactors=FALSE,
   row.names=row.names(sp1))
sp2 <- SpatialPolygonsDataFrame(sp1, data=df)
sp3 <- sp2[!is.na(sp2$Z),]
plot(sp3, col=sp3$col, bg="lightblue", axes=TRUE)

> And if it is not possible I will do the same using only the sp package...

You are right that for this data set, representation as a 
SpatialPixelDataFrame object might be best.

Roger

>
> Regards.
> Armelle Rouyer
>
>
>
> Roger Bivand a écrit :
>> On Thu, 17 Feb 2011, Armelle ROUYER wrote:
>> 
>>> Dear all,
>>> 
>>> I used the PBSmapping package to make a grid with the abundance of fishes. 
>>> Each cell has an allocated value. I want to convert this grid into a 
>>> shapefile in order to export it on ArcGIS.
>>> 
>>> So I tried the function "PolySet2SpatialPolygons". The grid is correctly 
>>> transformed into polygons but the allocated value is ignored. Then I tried 
>>> the function "SpatialPolygonsDataFrame" but it didn't work properly (only 
>>> 13 values are associated with the 50 cells).
>>> 
>>> I am not sure that I am using the right functions. So if you have any idea 
>>> or suggestion thank you for your help.
>> 
>> library(PBSmapping)
>> library(maptools)
>> gt <- GridTopology(c(0.5, 0.5), c(1, 1), c(5, 10))
>> sp <- as(gt, "SpatialPolygons")
>> proj4string(sp) <- CRS("+proj=longlat")
>> PS0 <- SpatialPolygons2PolySet(sp)
>> PD0 <- as.PolyData(data.frame(PID=1:50, col=rep(rainbow(10), each=5),
>>   stringsAsFactors=FALSE))
>> plotPolys(PS0, polyProps=PD0, projection="LL")
>> # the above to set up a PolySet opject, PolyData object, and plot them
>> sp1 <- PolySet2SpatialPolygons(PS0)
>> mo <- match(as.character(PD0$PID), row.names(sp1))
>> # check mo - this is an easy match, in practice it may need more study
>> sp2 <- SpatialPolygonsDataFrame(sp1, data=data.frame(col=PD0$col[mo],
>>   stringsAsFactors=FALSE))
>> plot(sp2, col=sp2$col, axes=TRUE)
>> 
>> seems to work. You match on the PolyData PID column and the row.names() of 
>> the SpatialPolygons object, which are the unique PolySet PID values. In 
>> this case, and probably in general, the PIDs in the PolyData and PolySet 
>> match by design. If there isn't an exact match, you'll need to construct an 
>> object of the correct length manually.
>> 
>> Roger
>> 
>>> 
>>> Regards.
>>> 
>>> 
>>> 
>> 
>
>

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