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

[Armelle ROUYER]

Thank you very much for your help Roger, it really helps and explains me 
how the functions work and generate polygons.
Regards.
Armelle

Roger Bivand a écrit :
> 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.
>>>>
>>>>
>>>>
>>>
>>
>>
> 

-- 
Armelle Rouyer

IFREMER - Station de Lorient
Laboratoire Biologie Halieutique
8 rue Francois Toullec
56100 LORIENT
FRANCE
02 97 87 38 23