Previous steps

If you would like to return to information from the previous session, please click here

Context

For plotting the spatial locations of monitoring sites, there are a number of potential formats that may be required in order to visualise. These include spatial data in vector format (e.g. ESRI shapefile), spatial points (e.g. coordinates from a hand-held GPS), or raster (i.e. gridded) format.

Most commonly, coastlines will be in a vector format and can be imported using the sp or sf packages. Digital elevation models (DEM) or bathymetry layers are often in raster (or gridded) format and can be imported into R using the raster package.

This wiki page provides a brief instruction to importing of vector data and the procedure is similar for other formats.

Importing coastlines

One source of a fine resolution coasline comes from the Global Self-consistent Hierarchical High-resolution Geography database. One advantage of using this source is that it is able to resolve coastlines of remote archipelagos and unpopulated regions (e.g. Galápagos Islands, regions of the Eastern Tropical Pacific) in a consistent, reproducible fashion.

Accessing this data can be done a number of ways using R. For this example, we will use the downloaded high resolution shapefile.

Once downloaded and trimmed to the region, the GSHHG coastline can be imported into R using the sf (or simple features) using the read_sf() function:

 ## -- import gshhg coastline -- ##
  # point to data locale
    data_locale <- paste0("data_raw/examples/mapping/shp/wio_coastline/")

  # call to wio coastline
    wio_coastline <-
      paste0(data_locale,
             "wio_coastline.",
             coast_res,
             ".shp") %>%
      read_sf()
# Simple feature collection with 1366 features and 6 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 18 ymin: -34.83044 xmax: 63.50208 ymax: 11
# Geodetic CRS:  WGS 84
# # A tibble: 1,366 x 7
   # id    level source parent_id sibling_id     area                                         geometry
   # <chr> <int> <chr>      <int>      <int>    <dbl>                               <MULTIPOLYGON [°]>
 # 1 1         1 WVS           -1          1   2.92e7 (((43.65372 11, 43.67544 10.96236, 43.77717 10.…
 # 2 10        1 WVS           -1         10   5.91e5 (((44.99914 -25.48708, 44.99831 -25.48708, 44.9…
 # 3 183       1 WVS           -1        183   2.52e3 (((55.26581 -21.00122, 55.26711 -21, 55.28128 -…
 # 4 218       1 WVS           -1        218   1.87e3 (((57.54125 -20.00083, 57.54208 -20, 57.54583 -…
 # 5 243       1 WVS           -1        243   1.57e3 (((39.37956 -6.000833, 39.37956 -6.001694, 39.4…
 # 6 313       1 WVS           -1        313   1.03e3 (((43.21956 -11.73503, 43.21956 -11.73331, 43.2…
 # 7 320       1 WVS           -1        320   9.84e2 (((39.67706 -5.002528, 39.67706 -5, 39.67275 -4…
 # 8 560       1 WVS           -1        560   4.47e2 (((39.58708 -7.946694, 39.58708 -7.945028, 39.5…
 # 9 573       1 WVS           -1        573   4.37e2 (((44.20289 -12.1675, 44.20417 -12.16625, 44.20…
# 10 645       1 WVS           -1        645   3.74e2 (((45.03622 -12.72086, 45.03622 -12.72, 45.0462…
# # … with 1,356 more rows

One should note that the geometry type of this object is MULTIPOLYGON and the spatial information (e.g. coordinate reference system (CRS)) and associated data columns (e.g. id) are retained - as well as a geometry column.

Alternatively, we can use the sp package, which produces a SpatialPolygonsDataFrame with similar properties:

       paste0(data_locale,
              "wio_coastline.",
              coast_res,
              ".shp") %>%
        shapefile()
class       : SpatialPolygonsDataFrame
features    : 1366
extent      : 18, 63.50208, -34.83044, 11  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs
variables   : 6
names       :    id, level, source, parent_id, sibling_id,           area
min values  :     1,     1,    WVS,        -1,          1, 0.017926397097
max values  : 99954,     1,    WVS,        -1,     168030,   29220969.727

Once we start working with spatial data, the differences between these two types of objects will become more apparent. But, for now, it is easy enough to switch between the two formats. For example, to visualising the coastline, we will convert the object to a sp object and use plot():

wio_coastline %>% as("Spatial") %>% plot()

Which looks something like this:

Creating zooms

As most coral reef monitoring efforts are conducted in localised areas of the Western Indian Ocean region, we will need to create a zoom or geographic limits to adjust the extent of the coastline represented.

For example, using the fish transect monitoring data from Kenya, we can set the geographic extent of the region by taking the max() and min() values of the dataset:

 ## -- set dimensions for monitoring domain -- ##
  # set country of interest
    country_of_interest <-
      c("Kenya")

  # extract limits from monitoring data
    kenya_limits <-
      kenya_fishes_data %>%
        dplyr::filter(Country %in% country_of_interest) %>%
        dplyr::filter(!Latitude  %>% is.na(),
                      !Longitude %>% is.na()) %>%
        dplyr::select(Station,
                      Latitude,
                      Longitude) %>%
        distinct()

  # set limits
    lim_n <- kenya_limits$Latitude  %>% max(na.rm = TRUE)
    lim_e <- kenya_limits$Longitude %>% min(na.rm = TRUE)
    lim_w <- kenya_limits$Longitude %>% max(na.rm = TRUE)
    lim_s <- kenya_limits$Latitude  %>% min(na.rm = TRUE)

  # create bounds
    kenya_bounds <-
      c(lim_w,
        lim_e,
        lim_s,
        lim_n)

These bounds are then used in ggplot() to constrain the plot. As the coastline object is a sf object, we can use the geom_sf() and coord_sf() to set the limits of the plot. Note that the CRS (i.e. datum = 4326) is specified for WGS84 (i.e. Latitude-Longitude).

  # visualise
    ggplot() +
      geom_sf(fill   = wesanderson::wes_palette("Cavalcanti1")[3],
              colour = "grey75",
              size   = 0.1,
              alpha  = 0.4,
              data   = wio_coastline %>%
                         dplyr::select(id)) +
      geom_sf(aes(colour = Station,
                  size   = abundance),
                 alpha = 0.7,
                 data  = kenya_fishes_sites %>%
                           dplyr::filter(Family %in% taxa_of_interest)) +
      theme_nothing() +
      coord_sf(xlim  = c(kenya_bounds[2], kenya_bounds[1]),
               ylim  = c(kenya_bounds[3], kenya_bounds[4]),
               datum = 4326) +
      facet_grid(. ~ Family) +
      scale_colour_manual(values = c_palette)

Which should look something like this:

Next steps

As the WGS84 coordinate reference system is good for locating positions on the (irregular) curved surface of the globe, it tends to distort coordinate values when represented on a flat surface (i.e. on a map or computer screen). This is largely due to the Longitude values are not equally spaced as one moves further from the equator.

To correct for this, coordinate systems are often projected, which we will find is easily managed here.