[R-sig-Geo] [R] Plot/manage spatial boundary data

[Roger Bivand]

On Fri, 10 Jun 2005, David Forrest wrote:

> On Fri, 10 Jun 2005, Roger Bivand wrote:
> 
> > On Thu, 9 Jun 2005, David Forrest wrote:
> >
> > > I have some disconnected boundary data from a finite element ocean model
> > > and I'd like to make a plot.
> > >
> > > Maptools looks promising, but since my data is not in a shapefile or a
> > > map, I'm unclear on what the best way to approach the problem.
> >
> > If the line segments are unconnected, you will first need to establish
> > (short) lists saying which segments (in which order and direction) bound
> > each polygon to be filled with colour when plotting. A package you can
> > consider for "rolling your own" spatial rings is sp, which is on CRAN.
> 
> The segments are largely connected -- my domain is bounded by a list of
> open boundary segments, land boundary segments, and several islands.
> library(sp) looks like the proper tool.

If you can identify the open segments, and make Sline objects of them, 
you'll get the bounding box too (bbox() of the output object). That gives 
you the coordinates you need to close them against the study area boundary 
(assumed rectangular in latlong), by rbind() to the input coordinates for 
the missing bits, and copy in coords[1,] at the end to close. Since the 
islands are closed rings anyway, I think this should work. Please let us 
know how you get on.

Roger

> 
> > Once you have the list describing segment membership, order and direction
> > for each ring, building a SpatialRings object is not difficult. The hard
> > bit is going from spaghetti line segments to the list imposing order.
> >
> > Alternatively, the PBSmapping package may have suitable functions for
> > coersing polySet objects into rings. If your coordinates were in the
> > Pacific, I'd say PBSmapping might already have what you need, but your
> > example coordinates are Atlantic.
> >
> > (Could I suggest moving this discussion to R-sig-geo, referenced in the
> > "Spatial" Task View on CRAN (top left corner in navigation bar)?
> 
> Sure.  I just subscribed.
> 
> How do you find the Task Views?  I found the mailing list off of
> http://www.r-project.org/ -- Mailing Lists and the CRAN link brings up
> mirrors.
> 
> Doh! I just answered my own question; these are different:
> 
>    http://www.r-project.org/
>    http://cran.r-project.org/
> 
> I thought they were simple mirrors of the homepage for downloading, but
> the menus are different.
> 
> The page at http://cran.r-project.org/ Task Views / Spatial
> http://cran.r-project.org/src/contrib/Views/Spatial.html looks like just
> what I need.  Thanks again.
> 
> > >
> > > >geom[1:10,]
> > >          lon      lat  depth
> > > 1  -75.42481 35.58192 16.172
> > > 2  -75.40726 35.58567 18.045
> > > 3  -75.41351 35.60312 17.333
> > > 4  -75.38888 35.58959 20.787
> > > 5  -75.39495 35.60706 19.834
> > > 6  -75.36964 35.59370 20.950
> > > 7  -75.37556 35.61159 20.941
> > > 8  -75.35530 35.61660 23.107
> > > 9  -75.34950 35.59800 22.960
> > > 10 -75.33418 35.62194 23.934
> > >
> > > >island1<-c(2,3,4,2)
> > > >water<-c(1,3,5,7,8,10)
> > > > land<-c(1,2,4,6,9,10)
> > > > plot(geom$lon[land],geom$lat[land],pch='.',t='l')
> > >  lines(geom$lon[water],geom$lat[water],pch='.',t='l',col="blue")
> > >  lines(geom$lon[island1],geom$lat[island1],pch='.',t='l',col="green")
> > >
> > > The above is toy-sized: dim(geom) is on the order of 120000,3 and there
> > > are about 30 different islands.  Maptools seems devoted to shapefiles,
> > > and it is unclear how to create 'polylists'.
> > >
> > > Is there a good way to manage and graph data defined on irregular grids?
> > >
> > > Dave
> 
> Dave
> 

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