Re: [R] PBSmapping and shapefiles

From: Denis Chabot <chabotd_at_globetrotter.net>
Date: Fri 22 Jul 2005 - 22:30:04 EST

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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Jul 22 22:37:37 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:33:55 EST