Extracting GSHHS global coastlines
On Thu, 9 Aug 2012, r-sig-geo.20.trevva at spamgourmet.com wrote:
Dear R-sig-geo,
I seem to be having problems with the GSHHS interface that I was
wondering if you could help with? I would like to extract the global,
low-resolution coastline from the database. Here's my script:
library(maptools)
xlims <- c(0,360)
ylims <- c(-90,90)
gshhs.file <- system.file("share/gshhs_c.b", package="maptools")
world <- Rgshhs(gshhs.file, xlim=xlims, ylim=ylims,level=1,checkPolygons=TRUE)
Which results in the following output:
Data are polygon data
Loading required package: rgeos
Loading required package: stringr
Loading required package: plyr
rgeos: (SVN revision 343M)
GEOS runtime version: 3.3.2-CAPI-1.7.2
Polygon checking: TRUE
Polygon 5 is Antarctica
area 13956598
dropping south edge to -90
Rgshhs: clipping 4 of 785 polygons ...
Error in slot(gI, "polygons") :
no slot of name "polygons" for this object of class "SpatialCollections"
Using getRgshhsMap() instead results in a different error:
world <- getRgshhsMap(gshhs.file, xlim=xlims, ylim=ylims,level=1)
Data are polygon data Polygon 5 is Antarctica area 13956598 dropping south edge to -90 Rgshhs: clipping 4 of 785 polygons ... Error in slot(gI, "polygons") : no slot of name "polygons" for this object of class "SpatialCollections" Data are polygon data Polygon is Antarctica area 50577991 dropping south edge to -90 Error in polys[[which(chosen_0 == (Antarctica - 1))]] : attempt to select less than one element
Am I doing something wrong? Or is this a bug?
Either of:
world <- Rgshhs(gshhs.file, level=1,checkPolygons=TRUE)
or
xlims <- c(-180,360)
world <- Rgshhs(gshhs.file, xlim=xlims, ylim=ylims, level=1,
checkPolygons=TRUE)
work. Clipping on Greenwich seems to cause problems with a sliver,
generating a single point (a degenerate polygon) and two polygons for Kent
and East Anglia. I've committed a fix to R-Forge SVN, but the output is
unattractive, with polygon splitting on 180 degrees causing trouble. A few
islands seem to flip from +180 to -180, creating havoc. In addition,
Antarctica is -180, 180, so looks bad. Consider using the wrld_simpl
dataset in maptools, or similar in the rworldmaps package instead if you
need country boundaries, or merge them out with:
all_wrld <- unionSpatialPolygons(wrld_simpl, rep("land",
nrow(wrld_simpl)), threshold=5)
which uses threshold= to filter out boundary slivers. Maybe even easier is
the maps package and map2SpatialPolygons() in maptools if all you need is
a background world land map.
Roger
Am running maptools 0.8-14 with sp 0.9-99 and rgeos 0.2-6 on R 2.14.1 Best wishes, Mark
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Department of Economics, NHH Norwegian School of Economics, 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