importGSHHS: Import Data from a GSHHS Database

In PBSmapping: Mapping Fisheries Data and Spatial Analysis Tools

Import Data from a GSHHS Database

Description

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

Usage

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 ^\circW (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. 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)).

Value

A PolySet with a PolyData attribute.

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^\circ, 360^\circ) or (-180^\circ, 180^\circ).

This version of the R function attempts to stitch together the overlapping edges of gshhs that run from -20^\circ to 360^\circ (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.

Author(s)

Nicholas M. Boers, Staff 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: 2023-10-30

References

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

See Also

In package PBSmapping:
importEvents, importLocs, importPolys

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 documentation built on Nov. 4, 2023, 9:08 a.m.