[R] PBSmapping and shapefiles

[Denis Chabot]

Hi,

I got no reply to this:
Le 16-Jul-05 à 2:42 PM, Denis Chabot a écrit :

> Hi,
>
> Is there a way, preferably with R, to read shapefiles and transform  
> them in a format that I could then use with package PBSmapping?
>
> I have been able to read such files into R with maptools'  
> read.shape and plot it with plot.Map, but I'd like to bring the  
> data to PBSmapping and plot from there. I also looked at the  
> package shapefile, but it does not seem to do what I want either.
>
> Sincerely,
>
> Denis Chabot
>

but I managed to progress somewhat on my own. Although it does not  
allow one to use shapefiles in PBSmapping, "maptools" at least makes  
it possible to read such files. In some cases I can extract the  
information I want from not-too-complex shapefiles. For instance, to  
extract all the lines corresponding to 60-m isobath in a shapefile, I  
was able to do:

library(maptools)
test <- read.shape("bathy.shp")
test2 <- Map2lines(test)
bathy60 <- subset(test2, test$att.data$Z == 60)

I do not quite understand the structure of bathy60 (list of lists I  
think)
but at this point I resorted to printing bathy60 on the console and  
imported that text into Excel for further cleaning, which is easy  
enough. I'd like to complete the process within R to save time and to  
circumvent Excel's limit of around 64000 lines. But I have a hard  
time figuring out loops in R, coming from a background of  
"observation based" programs such as SAS.

The output of bathy60 looks like this:

[[1]]
            [,1]     [,2]
[1,] -55.99805 51.68817
[2,] -56.00222 51.68911
[3,] -56.01694 51.68911
[4,] -56.03781 51.68606
[5,] -56.04639 51.68759
[6,] -56.04637 51.69445
[7,] -56.03777 51.70207
[8,] -56.02301 51.70892
[9,] -56.01317 51.71578
[10,] -56.00330 51.73481
[11,] -55.99805 51.73840
attr(,"pstart")
attr(,"pstart")$from
[1] 1

attr(,"pstart")$to
[1] 11

attr(,"nParts")
[1] 1
attr(,"shpID")
[1] NA

[[2]]
           [,1]     [,2]
[1,] -57.76294 50.88770
[2,] -57.76292 50.88693
[3,] -57.76033 50.88163
[4,] -57.75668 50.88091
[5,] -57.75551 50.88169
[6,] -57.75562 50.88550
[7,] -57.75932 50.88775
[8,] -57.76294 50.88770
attr(,"pstart")
attr(,"pstart")$from
[1] 1

attr(,"pstart")$to
[1] 8

attr(,"nParts")
[1] 1
attr(,"shpID")
[1] NA

What I need to produce for PBSmapping is a file where each block of  
coordinates shares one ID number, called PID, and a variable POS  
indicating the number of each coordinate within a "shape". All other  
lines must disappear. So the above would become:

PID X Y
1 1 -55.99805 51.68817
1 2 -56.00222 51.68911
1 3 -56.01694 51.68911
1 4 -56.03781 51.68606
1 5 -56.04639 51.68759
1 6 -56.04637 51.69445
1 7 -56.03777 51.70207
1 8 -56.02301 51.70892
1 9 -56.01317 51.71578
1 10 -56.00330 51.73481
1 11 -55.99805 51.73840
2 1 -57.76294 50.88770
2 2 -57.76292 50.88693
2 3 -57.76033 50.88163
2 4 -57.75668 50.88091
2 5 -57.75551 50.88169
2 6 -57.75562 50.88550
2 7 -57.75932 50.88775
2 8 -57.76294 50.88770

I don't know how to do this in R. My algorithm would involve looking  
at the structure of a line, discarding it if not including  
coordinates, and then creating PID and POS for lines with  
coordinates, depending on the content of lines i and i-1. In R?

Thanks in advance,

Denis Chabot