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