In this script I will explore different packages and tools for producing a small country border (e.g. Portugal) and a local area detail (e.g. Ria de Aveiro or Ria Formosa). It does not intend to be exhaustive, but rather to make a sample script that I often (re)use and might be useful to others. From a short quick map to put in data to higher resolution maps, that may require downloading the data from the internet.
From these bases, we can than add our data or costumize the maps accordingly our taste and desires.

1 Coastline (country border)

We will start by defining the limits of the study area. In this case, we aim to plot just portugal mainland, excluding the islands, so a country code will not be enough on most cases.

1.1 Library maps + ggplot

Excellent for a quick map.

# Establish our two example geographic limits (bbox)
pt.lim = data.frame(ylim=c(36.6, 42), xlim=c(-10, -7.0))
pt.bbox <- st_bbox(c(xmin=pt.lim$xlim[1],
           xmax=pt.lim$xlim[2],
           ymin=pt.lim$ylim[1],
           ymax=pt.lim$ylim[2]))

loc.lim = data.frame(ylim=c(40.6, 40.9), xlim=c(-8.8, -8.4))
loc.bbox <- st_bbox(c(
           xmin=loc.lim$xlim[1],
           xmax=loc.lim$xlim[2],
           ymin=loc.lim$ylim[1],
           ymax=loc.lim$ylim[2]))
require(maps)
ggplot()+
  borders("world", fill="grey90",colour="grey")+
  coord_fixed(ylim=pt.lim$ylim, xlim=pt.lim$xlim)

The resolution is not very high, however, but it is excellent for a quick plot. Below there is a simple version with the country alone.

1.2 Library raster & GADM: global administrative boundaries

Using raster package and GADM library which as the advantage of having altitude (alt), and regional borders (level 1 and 2) for example. It requires to download the data locally.

‘GADM’ is a database of global administrative boundaries.

The same function permits to download data from: ‘worldclim’ is a database of global interpolated climate data. ‘SRTM’ refers to the hole-filled CGIAR-SRTM (90 m resolution). ‘countries’ has polygons for all countries at a higher resolution than the ‘wrld_simpl’ data in the maptools package.

# Using GADM (gives a SpatialPolygonsDataFrame)
require(raster)
pt.gadm <- getData('GADM', country='Portugal', level=0)
# plot(pt)
# plot(getData('GADM', country='Portugal', level=1))

# Make it sf to facilitate and cut
require(sf)
pt.gadm <- sf::st_as_sf(pt.gadm) %>% 
    st_crop(pt.bbox)
 
object.size(pt.gadm)

2097128 bytes

# National map
ggplot()+
  geom_sf(data=pt.gadm)+
  coord_sf(ylim=pt.lim$ylim, xlim=pt.lim$xlim)+
  ggtitle("GADM border")

# Local map, basic
ggplot()+
  geom_sf(data=pt.gadm)+
  coord_sf(ylim=loc.lim$ylim, xlim=loc.lim$xlim)+
  ggtitle("GADM border - local area")

# Local map with a better look:
ggplot()+
  geom_sf(data=pt.gadm, fill="white", col="grey30")+
  coord_sf(ylim=loc.lim$ylim, xlim=loc.lim$xlim)+
  theme_minimal()+
  theme(panel.background = element_rect(fill = 'lightblue'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())+
  xlab("Longitude")+ylab("Latitude")+
  ggtitle("Ria de Aveiro, Portugal (GADM border)")

# Altitude data (gives a raster list):
alt.gadm <- getData("alt", country = "Portugal")[[1]]
plot(alt.gadm, main="Elevation") 

slope.raster <- terrain(alt.gadm, opt='slope')
aspect.raster <- terrain(alt.gadm, opt='aspect')
hill.raster <- hillShade(slope.raster, aspect.raster, 40, 270)
hill.m <- rasterToPoints(hill.raster)
hill.df <-  data.frame(hill.m)
colnames(hill.df) <- c("lon", "lat", "hill")

ggplot() +
  geom_raster(data = hill.df, aes(lon, lat, fill = hill), alpha = .45) +
  scale_fill_gradientn(colours = grey.colors(100)) +
  geom_sf(data = pt.gadm, fill = NA) +
  coord_sf(ylim=pt.lim$ylim, xlim=pt.lim$xlim) +
  theme_bw()+
  ggtitle("Relief map")

Note that the geographic aspect ratio has been preserved in all cases and the presence/absence of the border countries (i.e. Spain).

1.3 Library rnaturalearth

natural earth “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.”

example/source/adapted

# package rnaturalearth 
require(rnaturalearth)
pt.rnat <- ne_countries(country="portugal",scale = "large", returnclass = "sf") %>% 
  st_crop(pt.bbox)

object.size(pt.rnat)

41040 bytes

pt.stat <- ne_states(country = 'portugal', returnclass = "sf") %>% 
  st_crop(pt.bbox)
object.size(pt.stat)

140360 bytes

# Country map  
ggplot() +
    geom_sf(data = pt.rnat,fill="lightblue",colour="grey")+
  geom_sf(data = pt.stat,colour="blue")+
  ggtitle("rnaturalearth border (blue is ne_countries and red is ne_coastline")

# Local map  
ggplot() +
  geom_sf(data = pt.rnat,fill="grey90",colour="grey")+
  coord_sf(ylim=loc.lim$ylim, xlim=loc.lim$xlim)+
  ggtitle("rnaturalearth border - local area")

Note that rnaturalearth has several datasets that can be downloaded and added to the map.

1.4 Library mapdata, worldHires: CIA World Data Bank II

Mapdata is a supplement to maps package, providing the larger and/or higher-resolution databases. Comes with worldHires (“CIA World Data Bank II”) with world coastlines and national boundaries.

worldLores contains exactly the same polygon information, but the boundaries have been reduced to a very low resolution.

require(mapdata)
# get regional polygons
pt.hires = map_data('worldHires','Portugal')

## Alternative
# reg = map_data("world2Hires")
# reg = subset(reg, region %in% c('Portugal'))
# # convert lat longs
# reg$long = (360 - reg$long)*-1


# Get a data frame with longitude, latitude, and size of bubbles (a bubble = a city)
# https://www.r-graph-gallery.com/330-bubble-map-with-ggplot2.html
cities <- world.cities %>% filter(country.etc=="Portugal") %>% 
filter(name %in% c("Coimbra", "Aveiro", "Porto", "Setubal", "Faro", "Tavira","Lisboa","Matosinhos")) 

# convert lat longs
# reg$long = (360 - reg$long)*-1

# make plot
ggplot()+
  # add coastline
  geom_polygon(data = pt.hires, aes(x = long, y = lat, group = group), fill = "antiquewhite", color = "grey")+
  geom_point(data=cities, aes(long, lat))+
  geom_text(data=cities, aes(long, lat, label=name), size=2, angle=45, hjust=0.1, vjust=0.1)+
  # configure projection and plot domain
  coord_fixed(xlim = pt.lim$xlim, ylim = pt.lim$ylim)+
  theme(panel.background = element_rect(fill = "aliceblue"))+
  ggtitle("WorldHires")

ggplot()+
  # add coastline
  geom_polygon(data = pt.hires, aes(x = long, y = lat, group = group), fill = "antiquewhite", color = "grey")+
  coord_fixed(xlim = loc.lim$xlim, ylim = loc.lim$ylim)+
  theme(panel.background = element_rect(fill = "aliceblue"))+
  ggtitle("WorldHires")

1.5 Library sf, with NOAA Global Self-consistant Hierarchical High-resolution (GSHHG) 

  • Needs to be downloaded from the internet to a local computer.

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

  • It is the coastline used in bathymetry GEBCO products, thus expected to have a optimal resolution for global scale data.

See also

# Shoreline data for us detail
pt.gshhg.l1 <- sf::read_sf(          "/Users/martarufino/Downloads/gshhg-shp-2/GSHHS_shp/f/GSHHS_f_L1.shp") %>% st_crop(pt.bbox)
# plot(pt.gshhg.l1["id"])

object.size(pt.gshhg.l1)

267072 bytes

# National map
ggplot()+
  geom_sf(data=pt.gshhg.l1)+
  coord_sf(ylim=pt.lim$ylim, xlim=pt.lim$xlim)+
  ggtitle("GSHHS_f_L1 border")

# Local map, basic
ggplot()+
  geom_sf(data=pt.gshhg.l1)+
  coord_sf(ylim=loc.lim$ylim, xlim=loc.lim$xlim)+
  ggtitle("GSHHS_f_L1 border - local area")

1.6 Interactive map

For the detail of estuarine larger areas we should add the estuaries from the esturine databse. I have produced a estuary visualization rpub which can be found in my profile.

# interactive map
pt.gshhg.l1 %>% 
mapview::mapview(map.types = "Esri.OceanBasemap",
                 alpha.regions = 0.3,
                 col.regions=NA,
                 color="red",
                 legend = FALSE)+
pt.gadm %>% 
mapview::mapview(alpha.regions = 0.3,
                 col.regions=NA,
                 color="blue",
                 legend = FALSE)

In this example we can see how the gshhg is totally overlaied by gadm.

2 Bathymetric maps

2.1 Library marmap, NOAA-ETOPO1 data

Bathymetry lines (depth) are mainly from two sources: NOAA-hosted ETOPO1 and GEBCO(30 second and 1 minute resolutions) (download the data from the site and use ‘readGEBCO.bathy’, using the ncdf4). R can access and plot both easily.

Let’s start by a simplest case, using NOAA.

# https://www.molecularecologist.com/2015/07/marmap/
require(marmap)

# get bathymetry data
pt.baty = getNOAA.bathy(
                  lon1 = pt.lim$xlim[1], 
                  lon2 = pt.lim$xlim[2], 
                  lat1 = pt.lim$ylim[1], 
                  lat2 = pt.lim$ylim[2], resolution = 1)
# plot(pt.baty)
# autoplot(pt.baty, geom=c("r", "c"), colour="grey", size=0.1) + scale_fill_etopo()

# plot it:
ggplot(pt.baty, aes(x=x, y=y)) +
  coord_quickmap() +
  geom_raster(aes(fill=z)) +
  scale_fill_etopo() +
  geom_contour(aes(z=z),
    breaks=c(0, -10, -20, -50, -100, -200, -1000), colour="grey", size=0.2) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  xlab("Longitude")+ylab("Latitude")+
  ggtitle("NOAA-ETOPO marmap")

# detailed map
ggplot(pt.baty, aes(x=x, y=y)) +
  coord_quickmap() +
  geom_raster(aes(fill=z)) +
  scale_fill_etopo() +
  geom_contour(aes(z=z),
    breaks=c(0, -10, -20, -50, -100, -200, -1000),
    colour="grey", size=0.2) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  xlab("Longitude")+ylab("Latitude")+
  ggtitle("NOAA-ETOPO marmap")+
  coord_fixed(xlim = loc.lim$xlim, ylim = loc.lim$ylim)

And now, just contours.

# Now, a basic plot with simple bahymetric contours

# convert bathymetry to data frame
pt.baty2 = fortify.bathy(pt.baty) 

# get regional polygons
reg = map_data("world2Hires")
reg = subset(reg, region %in% c('Portugal'))
# convert lat longs
reg$long = (360 - reg$long)*-1

# make plot
ggplot()+
  # add 5m contour
  geom_contour(data = pt.baty2, 
               aes(x=x, y=y, z=z),
               breaks=c(-5),size=c(0.1), 
               colour="lightsteelblue4")+
  # add 10m contour
  geom_contour(data = pt.baty2, 
               aes(x=x, y=y, z=z),
               breaks=c(-10), size=c(0.2),
               colour="lightsteelblue3")+
  # add 20m contour
  geom_contour(data = pt.baty2, 
               aes(x=x, y=y, z=z),
               breaks=c(-20), size=c(0.3),
               colour="lightsteelblue2")+
  coord_fixed(xlim = pt.lim$xlim, ylim = pt.lim$ylim)+
  # add coastline
  geom_polygon(data = reg, aes(x = long, y = lat, group = group), fill = "antiquewhite", color = "grey")+
  theme(panel.background = element_rect(fill = "aliceblue"))+
  ggtitle("NOAA-etopo marmap")+
  ylab("Latitude")+xlab("Longitude")

2.2 library marmap, GEBCO

For Gebco data, we need to first download the data from GEBCO site and than read it (easily) in r. the resolution is 15 arc-second intervals.

# global terrain model for ocean and land at 15 arc-second intervals
pt.gebco <- readGEBCO.bathy("/Users/martarufino/Google Drive/1_marta/trabalho/IPMA/campanhas biva/dados originais/shapes/GEBCO_2021/gebco_2020_n42.0_s36.6_w-10.0_e-7.0.nc")

# also using directly the raster package:
# plot(raster("/Users/martarufino/Google Drive/1_marta/trabalho/IPMA/campanhas biva/dados originais/shapes/GEBCO_2021/gebco_2020_n42.0_s36.6_w-10.0_e-7.0.tif"))

# plot(pt.gebco)
# autoplot(pt.gebco, geom=c("r", "c"), colour="grey", size=0.1) + scale_fill_etopo()

# plot it:
ggplot(pt.gebco, aes(x=x, y=y)) +
  coord_quickmap() +
  geom_raster(aes(fill=z)) +
  scale_fill_etopo() +
  geom_contour(aes(z=z),
    breaks=c(0, -10, -20, -50, -100, -200, -1000), colour="grey", size=0.2) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  xlab("Longitude")+ylab("Latitude")+
  ggtitle("GEBCO marmap")

# detailed map
ggplot(pt.gebco, aes(x=x, y=y)) +
  coord_quickmap() +
  geom_raster(aes(fill=z)) +
  scale_fill_etopo() +
  geom_contour(aes(z=z),
    breaks=c(0, -10, -20, -50, -100, -200, -1000), colour="red") +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  xlab("Longitude")+ylab("Latitude")+
  ggtitle("GEBCO marmap")+
  coord_fixed(xlim = loc.lim$xlim, ylim = loc.lim$ylim)

## now, if we want to take 'just' the coast line:
# pt.gebco <- fortify(pt.gebco)
# pt.gebco <- pt.gebco[pt.gebco$z==0,]