[R] PBSmapping and shapefiles

[Roger Bivand]

On Fri, 22 Jul 2005, Denis Chabot wrote:

> 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?
> 

The way to do it would be to manipulate the Map$Shapes object directly, 
making sure that the Map$att.data records stay associated with them. If 
you could send me (off-list) a small sample shapefile, I'll see if there 
are obvious solutions - they'll probably be through the "sp" package.

Roger Bivand

> Thanks in advance,
> 
> Denis Chabot
> 
> ______________________________________________
> R-help at 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
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no