ggplot world maps

A basic coastline can be obtained quickly using just ggplot. “This is a quick and dirty way to get map data (from the maps package) on to your plot. This is a good place to start if you need some crude reference lines, but you’ll typically want something more sophisticated for communication graphics.” Note that maps package is based on rnaturalearth, as seen bellow.

ggplot()+
  borders("world", fill="grey90",colour="grey")

# a part of it (note that the country polygons are not cut this way):
ggplot()+
  borders("world", 
          xlim = c(-77, -70), 
          ylim = c(36, 41),
          fill="grey90",colour="grey")

# for that we should add  line of coord_fixed:
ggplot()+
  borders("world", 
          fill="grey90",colour="grey")+
  coord_fixed(xlim = c(-77, -70), 
          ylim = c(36, 41))

Natural Earth

The site Natural Earth and the respective library ‘rnaturalworld’ is probably the most commonly used geographic features in r at the moment. It is extremly pratical to use and well integrated in r, namely with geographical data packages such as ‘sf’ and ‘sp’.

“Natural Earth is a public domain map dataset available at 1:10m, 1:50m, and 1:110 million scales. Featuring tightly integrated vector and raster data, with Natural Earth you can make a variety of visually pleasing, well-crafted maps with cartography or GIS software.”

# package rnaturalearth 
require(rnaturalearth)
world <- ne_countries(scale = "large", returnclass = "sf")
world_map <- world %>% 
  dplyr::group_by(continent) %>% 
  dplyr::select(continent) %>% 
  st_buffer(dist=.001) 
  
ggplot(data = world) +
    geom_sf() +
    geom_sf(data=world_map, col=2, alpha=.8)+
    coord_sf(xlim = c(-77, -70), ylim = c(36, 41), expand = FALSE)+
    borders("world",alpha=.3, fill="lightblue",colour="grey")

# # package maps example
# require(maps)
# 
# ia <- map_data("world","portugal")
# mid_range <- function(x) mean(range(x))
# seats <- do.call(rbind, lapply(split(ia, ia$subregion), function(d) {
#   data.frame(lat = mid_range(d$lat), long = mid_range(d$long), subregion = unique(d$subregion))
# }))
# 
# ggplot(ia, aes(long, lat)) +
#   geom_polygon(aes(group = group), fill = NA, colour = "grey60") +
#   geom_text(aes(label = subregion), data = seats, size = 2, angle = 45)

Global Self-consistant Hierarchical High-resolution (Shorelines, gshhg)

The World Vector Shorelines gshhg produced by NOAA, is a high-resolution geography data set, amalgamated from two databases: World Vector Shorelines (WVS) and CIA World Data Bank II (WDBII). The former is the basis for shorelines while the latter is the basis for lakes, although there are instances where differences in coastline representations necessitated adding WDBII islands to GSHHG. The WDBII source also provides political borders and rivers. GSHHG data have undergone extensive processing and should be free of internal inconsistencies such as erratic points and crossing segments. The shorelines are constructed entirely from hierarchically arranged closed polygons.

GSHHG combines the older GSHHS shoreline database with WDBII rivers and borders, available in either ESRI shapefile format or in a native binary format. Geography data are in five resolutions: crude(c), low(l), intermediate(i), high(h), and full(f). Shorelines are organized into four levels: boundary between land and ocean (L1), boundary between lake and land (L2), boundary between island-in-lake and lake (L3), and boundary between pond-in-island and island (L4). Datasets are in WGS84 geographic (simple latitudes and longitudes; decimal degrees).

See also

# lets get a bbox to make a smaller sized objects and facilitate the examples:
geo.box <- c(xmin=-74.5, xmax=-68, ymin=38, ymax=42.5)

# Shoreline data for us detail
gshhg.l1 <- sf::read_sf("D:/shapes/gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp") %>% st_crop(geo.box)
plot(gshhg.l1["id"])

# now if we want to see how the estuaries fit in (see the other rpub):
estuary <- sf::read_sf("C:/Users/lioca/Google Drive/1_marta/trabalho/FCUL/LifeLine shared/kml/estuarios/estuary_rgdal.kml") %>% st_crop(geo.box)
plot(estuary["Name"])

# interactive map
world_map %>% st_crop(geo.box) %>% 
  mapview::mapview(map.types = "Esri.OceanBasemap", 
                   legend = FALSE,
                   col.regions=NA,
                   color="blue")+
estuary %>% 
mapview::mapview(zcol = "Name",
                 col.regions="green",
                 alpha.regions = 0.3,
                 legend = FALSE)+
gshhg.l1 %>% 
mapview::mapview(zcol = "id",
                 alpha.regions = 0.3,
                 col.regions=NA,color="red",
                 legend = FALSE)

In this example we can see how the gshhg has a better resolution, but naturally it is much bigger, thus we should choose the optimal one for each case.

OCE packages coastline

# example from http://emit.phys.ocean.dal.ca/~kelley/oar-internal/oar.pdf

require(oce)
# topo2 dataset in the ocedata package
data(coastlineWorldMedium, package="ocedata")
cwlon <- coastlineWorldMedium[["longitude"]]
cwlat <- coastlineWorldMedium[["latitude"]]
# polygon(cwlon, cwlat, col='gray')

# make a spatial 'sf' object:
kk <- data.frame(long = cwlon, lat = cwlat) %>% 
          na.omit() %>% 
          sf::st_as_sf(coords=c("long","lat")) %>% 
          st_set_crs(4326) %>% 
  st_crop(geo.box)

# interactive map of best coastline with estuaries:
gshhg.l1 %>% 
  mapview::mapview(map.types = "Esri.OceanBasemap",zcol = "id",alpha.regions = 0.3,legend = FALSE)+
kk %>% mapview::mapview(legend = FALSE, col.regions=3)
# Cleraly worst - besides not having the polygons correctly extracted for the moment.

# To plot the route, as in Figure 2.27, we may use map projection capabilities of the oce package (see Chapter 3).
# library(oce)
# data(coastlineWorldFine, package="ocedata")
# mapPlot(coastlineWorldFine, projection="+proj=merc", col='lightgray', longitudelim=range(lon), latitudelim=range(lat))
# # Adding the coastline and international border provides context, but necessitates the redrawing of the path and grid.
# mapLines(lon, lat, lwd=5)

Yet to explore: