[R-sig-Geo] Extracting GSHHS global coastlines

[Roger Bivand]

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