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

[Edzer Pebesma]

An example using SpatialPixelsDataFrame, as suggested by Roger, is given
below. It does not work with the current sp 0.9-76 due to a small bug
introduced there, but does work again in 0.9-77 which I just submitted
to CRAN, and might work with older versions.

require(sp)
PD=read.table("pdata.txt",header=T)
x=seq(-5,-1,.2)
y=seq(43,49,.2)
PD$x = x[PD[[1]]]
PD$y = y[PD[[2]]]
coordinates(PD)=~x+y
proj4string(PD)="+proj=longlat +datum=WGS84"
gridded(PD)=TRUE
fullgrid(PD)=TRUE
spplot(PD["Z"], scales = list(draw=TRUE), at = fivenum(PD$Z))

Alternatively you could use image() for plotting. Getting the colors of
the example below is an issue that I didn't look into.

On 02/18/2011 04:00 PM, Armelle ROUYER wrote:
> 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.
>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de