diff --git a/R/data.R b/R/data.R
index 7ecc7f8a..de2fd62e 100644
--- a/R/data.R
+++ b/R/data.R
@@ -956,9 +956,15 @@
 ##' abundance was calculated for seven regions: the Strait of Georgia (SOG),
 ##' West Coast Vancouver Island (WCVI), Queen Charlotte Strait (QCS), Discovery
 ##' Passage (DP), Central Mainland Coast (CMC), Northern Mainland Coast (NMC),
-##' and Haida Gwaii (HG); see Figure 1 of DFO (2022). The estimated abundances
-##' were calculated using Generalised Additive Models, and are included here
-##' (yielding the trends shown in Figure 3 of DFO, 2022). See that figure to see
+##' and Haida Gwaii (HG); see Figure 1 of DFO (2022) for a map showing the
+##' regions. The coastwide estimates are also given (labelled `Coastwide` so the
+##' plot title is capitalised), for which the means are calculated as the sum of
+##' the regions' means and the `low` and `high` values as TODO ???
+##'
+##'
+##' The estimated abundances
+##' were calculated using Generalised Additive Models (GAMs), and are included here
+##' (reproducing the trends shown in Figure 3 of DFO, 2022). See that figure to see
 ##' the amount of data available in each region; these data can be added to
 ##' `pacea` if desired -- please let us know. The final year of data in each
 ##' region is saved in the object `harbour_seals_data_final_year`.
@@ -979,18 +985,22 @@
 ##'   of about 250 days in the GAM, hence some years have more than one value
 ##'   and so we give a date rather than a year}
 ##'   \item{region:}{the region that the estimate corresponds to; one of SOG,
-##'   WCVI, QCS, DP, CMC, NMC, or HG}
+##'   WCVI, QCS, DP, CMC, NMC, HG, or Coastwide}
 ##'   \item{low:}{low end of the estimate of abundance, defined as
-##'   the TODO CHECK}
+##'   the upper value of the 95% confidence interval (calculated on the log
+##'   scale) from the GAM}
 ##'   \item{mean:}{mean estimate of abundance, numbers of seals (in that region)}
-##'   \item{high:}{high end of the estimate of abundance, defined TODO CHECK as
-##'   the TODO CHECK}
+##'   \item{high:}{high end of the estimate of abundance, defined as
+##'    the lower value of the 95% confidence interval (calculated on the log
+##'   scale) from the GAM}
 ##'  }
 ##'
 ##' @examples
 ##' \dontrun{
 ##' harbour_seals
-##' plot(harbour_seals)   # TODO reproduce trends in Figure 3
+##' plot(harbour_seals)
+##' plot(harbour_seals, region = "SOG")
+##' plot(harbour_seals, include_coastwide = FALSE)
 ##' }
 ##' @author Andrew Edwards
 ##' @source Estimates provided by Strahan Tucker, then wrangled and imported using
@@ -1015,18 +1025,17 @@
 "harbour_seals_data_final_year"
 
 #' 200m isobath for separating inshore and offshore Pacific using bathymetry
-#' 
-#' These coordinates roughly follow the continental shelf (200m depth) to provide the inshore and offshore boundary. The line is coarse and therefore ignores canyons. 
-#' 
+#'
+#' These coordinates roughly follow the continental shelf (200m depth) to provide the inshore and offshore boundary. The line is coarse and therefore ignores canyons.
+#'
 #' @format sf multilinestring object
 #'
-#' @examples 
+#' @examples
 #' \dontrun{
 #' isobath_200m
 #' plot(isobath_200m)
 #' plot(bccm_eez_poly, add = TRUE)
 #' }
-#' 
+#'
 #' @source Generated from running `data-raw/roms/isobath_sfline.R`.
 "isobath_200m"
-
diff --git a/R/plot-pacea-harbour-seals.R b/R/plot-pacea-harbour-seals.R
index 90c07327..d244cf1e 100644
--- a/R/plot-pacea-harbour-seals.R
+++ b/R/plot-pacea-harbour-seals.R
@@ -1,18 +1,20 @@
 ##' Plot the Pacific Harbour Seals abundance estimates
 ##'
-##' Default is to plot all seven regions, to mostly replicate
-##' (TODO CHECK) Figure 3 of DFO (2022; see `?harbour_seals`) in a `ggplot`
+##' Default is to plot all seven regions plus the coastwide estimate, to mostly replicate
+##' Figure 3 of DFO (2022; see `?harbour_seals`) in a `ggplot`
 ##' style. If a single region is
-##' indicated selected then a time series plot is shown, in the same style as
+##' indicated then a time series plot is shown, in the same style as
 ##' plots of object class `pacea_biomass` (such as `hake_biomass`).
-##' Note that the values are means and standard errors (not medians and
+##' Note that the values are means and confidence intervals (not medians and
 ##' confidence intervals like other objects); see `?harbour_seals`.
 ##'
 ##' @param obj `pacea_harbour_seals`, as this is of class `pacea_harbour_seals`
 ##'   and the plotting is tailored for this data set.
 ##' @param region which region to plot for a single plot; must be one of SOG,
-##'   WCVI, QCS, DP, CMC, NMC, or HG. If `NULL` (the default) then do a panel
+##'   WCVI, QCS, DP, CMC, NMC, HG, or coastwide. If `NULL` (the default) then do a panel
 ##'   plot showing all seven regions.
+##' @param include_coastwide logical, whether to include the coastwide population in the
+##'   eight panel plot (keep `region = NULL`).
 ##' @param value the column to plot as a solid line, in this case always the
 ##'   default `mean`
 ##' @param y_max maximum y value for certain types of plot (use this if you get
@@ -35,14 +37,14 @@
 ##' @author Andrew Edwards
 ##' @examples
 ##' \dontrun{
-##' plot(hake_biomass)  TODO
-##' plot(hake_biomass,
-##'      xlim = c(lubridate::dmy(01011950),
-##'               lubridate::dmy(01012040))) # to expand x-axis
+##' plot(harbour_seals)
+##' plot(harbour_seals, region = "SOG")
+##' plot(harbour_seals, include_coastwide = FALSE)
 ##' }
 plot.pacea_harbour_seals <- function(obj,
                                      region = NULL,
                                      value = "mean",
+                                     include_coastwide = TRUE,
                                      xlab = "Year",
                                      ylab = attr(obj, "axis_name"),
                                      y_tick_by = 1000,
@@ -64,8 +66,8 @@ plot.pacea_harbour_seals <- function(obj,
   stopifnot("value must be a column of the object in the first argument" =
             value %in% names(obj))
 
-  stopifnot("region must be one of SOG, WCVI, QCS, DP, CMC, NMC, HG" =
-            region %in% c("SOG", "WCVI", "QCS", "DP", "CMC", "NMC", "HG"))
+  stopifnot("region must be one of SOG, WCVI, QCS, DP, CMC, NMC, HG, Coastwide" =
+            region %in% c("SOG", "WCVI", "QCS", "DP", "CMC", "NMC", "HG", "Coastwide"))
 
   if(!is.null(region)){
     region_choice <- region     # Else region == region in next line does not work
@@ -94,18 +96,22 @@ plot.pacea_harbour_seals <- function(obj,
                   x_tick_extra_years = x_tick_extra_years,
                   start_decade_ticks = start_decade_ticks)
   } else {
-    ggplot(data = obj,
+    if(!include_coastwide){
+      obj <- dplyr::filter(obj, region != "Coastwide")
+    }
+
+    ggplot2::ggplot(data = obj,
            aes(x = date,
                y = mean)) +
-      geom_ribbon(aes(ymin = low,
+      ggplot2::geom_ribbon(aes(ymin = low,
                       ymax = high),
                   fill = uncertainty_shade_col) +
-      geom_line(col = uncertainty_line_col) +
+      ggplot2::geom_line(col = uncertainty_line_col) +
       # Add back in if data get included
       # geom_point(data = harbour_seals_data,
       #        aes(x = make_date(year),
       #           y = mean)) +
-      facet_wrap(~region,
+      ggplot2::facet_wrap(~region,
                  scales = "free_y",
                  ncol = 2)
   }
diff --git a/data-raw/mammals/harbour-seals/harbour-seals.R b/data-raw/mammals/harbour-seals/harbour-seals.R
index e0fc00d4..99030e7d 100644
--- a/data-raw/mammals/harbour-seals/harbour-seals.R
+++ b/data-raw/mammals/harbour-seals/harbour-seals.R
@@ -78,9 +78,9 @@ p <- ggplot(data = harbour_seals_gam_new,
              ncol = 2)
 p
 
-#Looks like matches SAR figure. Could add data, but TODO refer people to
+#Looks like matches SAR figure. Could add data, but refer people to
 # original SAR in the help file.
-#  TODO in help clarify what these are, as caption talks about se only
+
 
 class(harbour_seals_gam_new) <- c("pacea_harbour_seals",
                               class(harbour_seals_gam_new))
@@ -98,6 +98,21 @@ harbour_seals <- harbour_seals_gam_new %>%
          mean,
          high)                  # not keeping year as duplicated for some years
 
+# Calculate coastwide numbers and uncertainty
+
+coastwide <- dplyr::summarise(group_by(harbour_seals, date),
+                              low = sum(low),          # TODO this and high are
+                                        # wrong, emailed Kurtis
+                              mean = sum(mean),
+                              high = sum(high)) %>%
+  dplyr::mutate(region = "Coastwide") %>%
+  dplyr::relocate(date, region)
+
+stopifnot(names(harbour_seals) == names(coastwide))
+
+harbour_seals <- rbind(harbour_seals,
+                       coastwide)
+
 # Also need the final date of the data, to then colour code the
 # figures.
 
diff --git a/data/harbour_seals.rda b/data/harbour_seals.rda
index 8e400559..90d58e44 100644
Binary files a/data/harbour_seals.rda and b/data/harbour_seals.rda differ
diff --git a/man/harbour_seals.Rd b/man/harbour_seals.Rd
index 41fe6d52..383c3168 100644
--- a/man/harbour_seals.Rd
+++ b/man/harbour_seals.Rd
@@ -11,12 +11,14 @@ A tibble also of class `pacea_harbour_seals` with columns:
   of about 250 days in the GAM, hence some years have more than one value
   and so we give a date rather than a year}
   \item{region:}{the region that the estimate corresponds to; one of SOG,
-  WCVI, QCS, DP, CMC, NMC, or HG}
+  WCVI, QCS, DP, CMC, NMC, HG, or Coastwide}
   \item{low:}{low end of the estimate of abundance, defined as
-  the TODO CHECK}
+  the upper value of the 95% confidence interval (calculated on the log
+  scale) from the GAM}
   \item{mean:}{mean estimate of abundance, numbers of seals (in that region)}
-  \item{high:}{high end of the estimate of abundance, defined TODO CHECK as
-  the TODO CHECK}
+  \item{high:}{high end of the estimate of abundance, defined as
+   the lower value of the 95% confidence interval (calculated on the log
+  scale) from the GAM}
  }
 }
 \source{
@@ -41,9 +43,15 @@ pupping, and proportion of seals hauled out during surveys. Estimated
 abundance was calculated for seven regions: the Strait of Georgia (SOG),
 West Coast Vancouver Island (WCVI), Queen Charlotte Strait (QCS), Discovery
 Passage (DP), Central Mainland Coast (CMC), Northern Mainland Coast (NMC),
-and Haida Gwaii (HG); see Figure 1 of DFO (2022). The estimated abundances
-were calculated using Generalised Additive Models, and are included here
-(yielding the trends shown in Figure 3 of DFO, 2022). See that figure to see
+and Haida Gwaii (HG); see Figure 1 of DFO (2022) for a map showing the
+regions. The coastwide estimates are also given (labelled `Coastwide` so the
+plot title is capitalised), for which the means are calculated as the sum of
+the regions' means and the `low` and `high` values as TODO ???
+
+
+The estimated abundances
+were calculated using Generalised Additive Models (GAMs), and are included here
+(reproducing the trends shown in Figure 3 of DFO, 2022). See that figure to see
 the amount of data available in each region; these data can be added to
 `pacea` if desired -- please let us know. The final year of data in each
 region is saved in the object `harbour_seals_data_final_year`.
@@ -61,7 +69,9 @@ https://www.dfo-mpo.gc.ca/csas-sccs/Publications/SAR-AS/2022/2022_034-eng.html
 \examples{
 \dontrun{
 harbour_seals
-plot(harbour_seals)   # TODO reproduce trends in Figure 3
+plot(harbour_seals)
+plot(harbour_seals, region = "SOG")
+plot(harbour_seals, include_coastwide = FALSE)
 }
 }
 \author{
diff --git a/man/plot.pacea_harbour_seals.Rd b/man/plot.pacea_harbour_seals.Rd
index 3c077cdb..ae096b23 100644
--- a/man/plot.pacea_harbour_seals.Rd
+++ b/man/plot.pacea_harbour_seals.Rd
@@ -8,6 +8,7 @@
   obj,
   region = NULL,
   value = "mean",
+  include_coastwide = TRUE,
   xlab = "Year",
   ylab = attr(obj, "axis_name"),
   y_tick_by = 1000,
@@ -30,12 +31,15 @@
 and the plotting is tailored for this data set.}
 
 \item{region}{which region to plot for a single plot; must be one of SOG,
-WCVI, QCS, DP, CMC, NMC, or HG. If `NULL` (the default) then do a panel
+WCVI, QCS, DP, CMC, NMC, HG, or coastwide. If `NULL` (the default) then do a panel
 plot showing all seven regions.}
 
 \item{value}{the column to plot as a solid line, in this case always the
 default `mean`}
 
+\item{include_coastwide}{logical, whether to include the coastwide population in the
+eight panel plot (keep `region = NULL`).}
+
 \item{xlab}{x-axis label}
 
 \item{ylab}{y-axis label, the default is an attribute of the `pacea_index`
@@ -79,20 +83,19 @@ plot of the time series as median with bars showing uncertainty (if
   `low` and `high` are columns of `obj) to the current device; returns nothing.
 }
 \description{
-Default is to plot all seven regions, to mostly replicate
-(TODO CHECK) Figure 3 of DFO (2022; see `?harbour_seals`) in a `ggplot`
+Default is to plot all seven regions plus the coastwide estimate, to mostly replicate
+Figure 3 of DFO (2022; see `?harbour_seals`) in a `ggplot`
 style. If a single region is
-indicated selected then a time series plot is shown, in the same style as
+indicated then a time series plot is shown, in the same style as
 plots of object class `pacea_biomass` (such as `hake_biomass`).
-Note that the values are means and standard errors (not medians and
+Note that the values are means and confidence intervals (not medians and
 confidence intervals like other objects); see `?harbour_seals`.
 }
 \examples{
 \dontrun{
-plot(hake_biomass)  TODO
-plot(hake_biomass,
-     xlim = c(lubridate::dmy(01011950),
-              lubridate::dmy(01012040))) # to expand x-axis
+plot(harbour_seals)
+plot(harbour_seals, region = "SOG")
+plot(harbour_seals, include_coastwide = FALSE)
 }
 }
 \author{