Re: [R] PBSmapping and shapefiles

From: Roger Bivand <Roger.Bivand_at_nhh.no>
Date: Sun 24 Jul 2005 - 02:00:01 EST


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@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@nhh.no

______________________________________________
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 Sun Jul 24 02:07:01 2005

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