importGSHHS() R function from [PBSmapping]

Import Data from a GSHHS Database

Import data from a GSHHS database and convert data into a PolySet

with a PolyData attribute. The database was originally called Global Self-consistent, Hierarchical, High-resolution Shoreline

(GSHHS, Wessel and Smith 1996), but Shoreline was subsequently expanded to include more Geography (GSHHG).


importGSHHS(gshhsDB, xlim, ylim, maxLevel=4, n=0, useWest=FALSE)

Arguments

gshhsDB: character -- path name to binary GSHHS database. If unspecified, looks for gshhs_f.b in the root of the PBSmapping library directory.
xlim: numeric -- range of X-coordinates (for clipping). The range should be between 0 and 360, starting at the Greenwich meridian and wrapping eastward around the globe.
ylim: numeric -- range of Y-coordinates (for clipping).
maxLevel: numeric -- maximum level of polygons to import: 1 (land), 2 (lakes on land), 3 (islands in lakes), or 4 (ponds on islands); ignored when importing lines.
n: numeric -- minimum number of vertices that must exist in a line/polygon in order for it to be imported.
useWest: logical -- if TRUE, convert the X-coordinates (longitude) to $deg$ W (western hemisphere -180 to 0, i.e., west of the Greenwich meridian).

Details

This routine requires a binary GSHHG (Global Self-consistent, Hierarchical, High-resolution Geography) database file. The GSHHG database has been released in the public domain.

(www.soest.hawaii.edu/pwessel/gshhg)

At the time of writing, the most recent binary database was the archive file called gshhg-bin-2.3.7.zip.

The archive contains multiple binary files that contain geographical coordinates for shorelines (gshhs), rivers (wdb_rivers), and borders (wdb_borders). The latter two come from World DataBank II (WDBII). The five resolutions available are: full (f), high (h), intermediate (i), low (l), and coarse (c).

This routine returns a PolySet object with an associated PolyData attribute. The attribute contains four fields: (a) PID, (b) SID, (c) Level, and (d) Source. Each record corresponds to a line/polygon in the PolySet. The Level indicates the line's/polygon's level (1=land, 2=lake, 3=island, 4=pond). The Source identifies the data source (1=WVS, 0=CIA (WDBII)).

Returns

A PolySet with a PolyData attribute.

References

Wessel, P., and Smith, W.H.F. (1996) A Global Self-consistent, Hierarchical, High-resolution Shoreline Database. J. Geophys. Res. 101 8741-8743.

(www.soest.hawaii.edu/pwessel/gshhg/Wessel+Smith_1996_JGR.pdf)

Author(s)

Nicholas M. Boers, Senior Software Engineer

Jobber, Edmonton AB

Maintainer: Rowan Haigh, Program Head -- Offshore Rockfish

Pacific Biological Station (PBS), Fisheries & Oceans Canada (DFO), Nanaimo BC

locus opus: Offsite, Vancouver BC

Last modified Rd: 2024-09-25

Note

The function calls a C routine, also called importGSHHS, which returns a set of map coordinates that is not always predictably laid out. This issue stems from how the world is divided at the Greenwich meridian and at the International Date Line. The unpredictability occurs when user-specified X-limits span either of the longitudinal meridians -- (0 $deg$ , 360 $deg$ ) or (-180 $deg$ , 180 $deg$ ).

This version of the R function attempts to stitch together the overlapping edges of gshhs that run from -20 $deg$ to 360 $deg$ (see example map 5 below). At present, no attempt has been made to deal with the overlap at the International Date Line where Russia overlaps the Aleutian Islands of Alaska. To some extent, the C-code can deal with this, but not in all cases.

Therefore, the user will likely experience some limitations when using importGSHHS. The solution is to import the whole dataset with this function using xlim=c(0,360), and then apply the function refocusWorld

with user-desired X-limits. The Y-limits are generally not problematic unless the user wants to focus on either pole.

Examples


## Not run:

useWest=FALSE
useVers=c("2.2.0","2.2.3","2.3.0","2.3.4") # GSHHG versions
mapswitch = 5
for (i in c("land","rivers","borders"))
  if (exists(i)) eval(parse(text=paste0("rm(",i,")")))
switch( mapswitch,
 # 1. Canada------------------------------------------------
      {vN=4; useWest=T; xlim=c(-150,-50)+360;ylim=c(40,75)},
 # 2. NW Canada & America-----------------------------------
      {vN=4; useWest=T;xlim=c(-136,-100)+360;ylim=c(40,75)},
 # 3. Black Sea (user Ivailo)-------------------------------
      {vN=4; xlim=c(27.5, 34.3); ylim=c(40.9, 46.7)},
 # 4. W Europe, NW Africa (user Uli)------------------------
      {vN=4; xlim=c(-20,10); ylim=c(20,50)},
 # 5. W Europe + Iceland------------------------------------
      {vN=4; xlim=c(-25, 20); ylim=c(40, 68)},
 # 6. New Zealand-------------------------------------------
      {vN=4; xlim=c(163, 182); ylim=c(-48,-34)},
 # 7. Australia---------------------------------------------
      {vN=4; xlim=c(112,155); ylim=c(-44,-10)},
 # 8. Japan-------------------------------------------------
      {vN=4; xlim=c(127,148); ylim=c(30,47)},
 # 9. Central America---------------------------------------
      {vN=4; useWest=T; xlim=c(-95,-60)+360;ylim=c(-10,25)},
 #10. North Pacific-----------------------------------------
      {vN=4; useWest=T; xlim=c(150,220); ylim=c(45,80)},
 #11. Pacific Ocean-----------------------------------------
      {vN=4; xlim=c(112,240); ylim=c(-48,80)},
 #12. North Atlantic (world coordinates)--------------------
      {vN=4; xlim=c(285,360); ylim=c(40,68)},
 #13. North Atlantic (western hemisphere coordinates)-------
      {vN=4; xlim=c(-75,0); ylim=c(40,68)},
 #14. Atlantic Ocean----------------------------------------
      {vN=4; xlim=c(285,380); ylim=c(-50,68)},
 #15. Northern hemisphere-----------------------------------
      {vN=4; xlim=c(-180,180); ylim=c(0,85)},
 #16. Asia--------------------------------------------------
      {vN=4; xlim=c(0,180); ylim=c(0,80)},
 #17. North America-----------------------------------------
      {vN=4; xlim=c(-180,0); ylim=c(0,80)},
 #18. International date line-------------------------------
      {vN=4; xlim=c(45,315); ylim=c(0,80)},
 #19. Indian Ocean------------------------------------------
      {vN=4; xlim=c(20,130); ylim=c(-40,40)},
 #20. Moose County ("400 miles north of everywhere")--------
      {vN=4; xlim=c(272.5,280.5); ylim=c(43,47.5)}
)
db=paste0("gshhg-bin-",useVers[vN])        # database version folder
gshhg   = paste0("C:/Ruser/GSHHG/",db,"/") # directory with binary files
land    = importGSHHS(paste0(gshhg,"gshhs_i.b"),
          xlim=xlim,ylim=ylim,maxLevel=4,useWest=useWest)
rivers  = importGSHHS(paste0(gshhg,"wdb_rivers_i.b"),
          xlim=xlim,ylim=ylim,useWest=useWest)
borders = importGSHHS(paste0(gshhg,"wdb_borders_i.b"),
          xlim=xlim,ylim=ylim,useWest=useWest,maxLevel=1)
if(exists("land")){
  plotMap(land,xlim=xlim-ifelse(useWest,360,0),ylim=ylim,
    col="lemonchiffon",bg="aliceblue")
  if(!is.null(rivers)) addLines(rivers,col="blue")
  if(!is.null(borders)) addLines(borders,col="red",lwd=2)
}
## End(Not run)

PBSmapping package Read PDF manual

Maintainer: Rowan Haigh
License: GPL (>= 2)
Last published: 2024-09-25

Useful links

importGSHHS function