surveyjoin: a standardized database of scientific trawl surveys in the Northeast Pacific Ocean

Data sources

Scientists and government agencies contributing to the management of sustainable fisheries collect a wide variety of data types from a wide variety of data sources, many of which could be combined across boundaries. Because trawl surveys are one of the more common types of fisheries data collected around the world, we initially limited our focus to creating a database of trawl survey data. Bottom trawls use nets towed on or near the ocean bottom at predefined locations and depths, and are designed to sample benthic fish and invertebrates. In many regions around the world, such surveys are used by governments and industry to monitor the sustainability of fisheries. Total catches for a particular species of interest can be converted into catch-per-unit-effort (CPUE, in kg per ha), and after accounting for sources of variability (e.g., habitat, environmental, random spatial variation), data can be used to generate relative estimates of total biomass in the survey domain. Time series of relative biomass density estimates can then be used as inputs into integrated population models (i.e., fisheries stock assessments).

While bottom trawl surveys have been designed to monitor temporal changes in groundfish biomass, trawl survey gear is not very selective. As a result, species encountered in these surveys include both commercial groundfish species of interest, as well as benthic invertebrates and non-commercial species. Additionally, the bottom trawl surveys are also used as platforms to collect other data critical for stock assessment and ecosystem based fisheries management, including individual lengths, weights, diet, fecundity, and age information. In situ oceanographic data (temperature, oxygen, salinity) may also be collected to better inform relationships between species and the environments they inhabit.

In total, surveyjoin brings together catch-per-unit-effort, haul, and in-situ bottom temperature data from 14 distinct bottom trawl surveys from across the west coast of North America (Table 1; Fig. 1). This includes five surveys in Alaska ecosystems conducted by NOAA Fisheries Alaska Fisheries Science Center since the mid-1980s (the eastern Bering Sea shelf, northern Bering Sea shelf, the eastern Bering Sea slope, the Aleutian Islands, and the Gulf of Alaska), four surveys in waters of British Columbia conducted by Fisheries and Oceans Canada since the early-2000s (Hecate Strait and Queen Charlotte Sound generally sampled in odd years, and west coasts of Haida Gwaii and Vancouver Island sampled in even years), and five surveys on the U.S. West Coast conducted by the NOAA Fisheries Northwest Fisheries Science Center. A full list of these surveys may be accessed in surveyjoin with the function get_survey_names().

All data contained in this database and retrievable through the surveyjoin R package are the final, validated survey data that are publicly accessible soon after surveys are completed. Each survey listed in Table 1 collects station-level catch and CPUE, haul, and in situ environmental data. Raw data from each organization are also available from different sources; data collected by NOAA Fisheries’ Alaska Fisheries Science Center are available through the Fisheries One Stop Shop (FOSS) platform (https://www.fisheries.noaa.gov/foss/f?p=215:200), data collected by NOAA Fisheries’ Northwest Fisheries Science Center are available online (https://www.webapps.nwfsc.noaa.gov/data/map), and data collected by Fisheries and Oceans Canada is available online (https://open.canada.ca/data/en/dataset/a278d1af-d567-4964-a109-ae1e84cbd24a/resource/9992b5d6-fa0d-4dd5-a207-dd1c20ee3f6d).

While the gear and sampling protocols are largely consistent across these surveys, some differences in gear design and sampling protocols exist, which may influence catchability and selectivity (Table 1). As an example, the eastern and northern Bering Sea Shelf surveys use a smaller footrope and slightly different net design compared to other surveys in Alaska (Stauffer GD (compiler), 2004). However, ongoing research to calibrate catch data between the Bering shelf and slope surveys aims to standardize catches for some common species; this region faces increased challenges because it is experiencing borealization (Litzow et al., 2024) and surveys are having to adapt to changing environmental conditions and species distributions (Vilas et al., 2024). In the absence of such calibrations or survey redesign, catchability may also be estimated statistically based on spatiotemporal proximity (e.g., by including survey as a factor predictor).

Throughout the history of these surveys, there have been minimal changes in the protocols regarding gear construction, configuration, towing speeds, and survey design within an individual survey (Vilas et al., 2024). Without changing these concepts, technological advances (such as the development of net mensuration equipment, on-bottom sensors, and satellite-based global positioning systems) have enhanced the accuracy of effort measurement. Moreover, changes in species identification and categorization protocols, although not universally applied across all surveys, have also influenced the data collection process (Stevenson & Hoff, 2009).

Package development

While each of the individual surveys in our dataset collects hundreds (or thousands) of species, we only included a subset in our initial version of the surveyjoin package. Each of the surveys (Table 1) can be assigned to one of four major geographical areas: Eastern Bering Sea, Gulf of Alaska, British Columbia, and west coast of California/Oregon/Washington states. Within each of the four regions, we first filtered and identified species that occurred in more than 5% of tows. In the spirit of cross-border data standardization, we then filtered that list to identify those occurring in two or more of our four geographic areas, and finally queried all data from those species across all regions. These steps resulted in nearly 3.3 million observations of 55 species. To help future efforts in using these data, we include the common name, scientific name, and Integrated Taxonomic Information System (ITIS; https://itis.gov/) identification for each species (Table 2); these may be accessed in surveyjoin with the function get_species().

Trawl datasets were standardized across surveys by converting the raw data from each survey into the same units (effort in hectares, catch count in numbers, and weight in kilograms). In addition to catch and effort, we identified common variables across surveys often used in analyses of these data including location (start and end longitude and latitude in decimal degrees), date, and the in situ bottom temperature associated with each haul. Finally, because many analyses using the surveyjoin database will be focused on the estimation of population trends, we also include the design grids associated with each of the surveys as data objects. Each of these objects are included in the package as a dataframe (each grid consisting of >10,000 rows), with cell centroids and cell areas. Examples of using these grids with model predictions are included in the case studies below. In developing the surveyjoin R package, we are following best practices in code development on GitHub (e.g., versioning, releases, code review), and to facilitate interpretation across regions and surveys, we link to survey-specific metadata. As our data originates from multiple agencies and regions, we are also extending these practices to versioning individual datasets. Tracking data source versions should facilitate transparency and reproducibility for end users, and allow changes to be compared over time.

Case studies

To demonstrate the utility of the surveyjoin package, we developed three illustrative case studies focused on different statistical models and questions. Each of these models can be viewed as a different type of geostatistical model, where raw observations are associated with a unique haul identifier and latitude and longitude coordinates. Data and code to replicate our analysis and figures is on our GitHub repository (https://github.com/DFO-NOAA-Pacific/surveyjoin-paper).

Case study 1

For a first case study, a geospatial index standardization (Thorson et al., 2015) was applied across international borders to develop a coastwide index of Pacific hake (Merluccius productus, also known as Pacific whiting). Pacific hake are a fast growing, semi-pelagic species found in the northeast Pacific Ocean off the coasts of Canada and the U.S. (ranging from Alaska to California) and, by volume, represent one of the largest fisheries on the west coast of North America. Pacific hake are managed under the bilateral Pacific Whiting Agreement and an estimate of their population status is updated annually by a team of scientists from NOAA Fisheries and Fisheries and Oceans Canada. The assessment model implemented for Pacific hake can be described as a Bayesian age-structured population model based on the Stock Synthesis modeling framework (Methot & Wetzel, 2013), with data inputs including absolute indices of abundance from an acoustic survey (Grandin et al., 2024). Similar indices of hake abundance have not yet been developed from trawl surveys, as Pacific hake are semi-pelagic and may not be well sampled by bottom trawl gear (most fish sampled by trawl gear are mature, through the survey does occasionally catch immature individuals). However, information from bottom trawl gear may help provide a clearer picture of population size for the semi-pelagic species, as each survey samples different vertical portions of the water column (Monnahan et al., 2021). Our objectives in this analysis are to (1) construct an index using trawl data from the surveyjoin package and compare it to the acoustic survey, and (2) demonstrate how indices may be used to evaluate distributional shifts.

The surveyjoin package was used to query all haul and catch data from British Columbia and the West Coast of the U.S. from 2003 to 2023. Data from surveys in the Gulf of Alaska were not included because observations of Pacific hake are sparse. Since 2003, there have been an average of 673 hauls per year in the Gulf of Alaska, with an average of 23.5 occurrences of hake per year (ranging from 1 to 69). Bottom trawl surveys along the West Coast of the U.S. have been completed annually, except in 2020 due to the COVID-19 pandemic (Keller, Wallace & Methot, 2017). Surveys in British Columbia are stratified into four regions, with two regions usually sampled in odd years (Hecate Strait and Queen Charlotte Sound) and two in even years (West Coast Vancouver Island and West Coast Haida Gwaii) as discussed above (also see Sinclair et al., 2003; Anderson, Keppel & Edwards, 2019). For this analysis, we used a spatiotemporal Generalized Linear Mixed Model (GLMM) using the sdmTMB R package (Anderson et al., 2025; R Core Team, 2024). The sdmTMB package combines the Stochastic Partial Differential Equation (SPDE; Lindgren, Rue & Lindström, 2011; Lindgren & Rue, 2015) approach with fast maximum likelihood estimation using Template Model Builder (Kristensen et al., 2016). Given the zero-inflated and highly skewed catch data relative to most distribution families, CPUE was modeled with a delta-Gamma model (Pennington, 1983). Presence-absence was modeled using a Bernoulli distribution and positive catch rate as a Gamma distribution.

For both the presence-absence and positive model components, the general form of the spatiotemporal GLMM be expressed as

$$u_t=f^{-1}\left(Xb+\mathit{\omega }+{\mathrm{ϵ}}_t\right)$$where ${\mathit{u}}_t$ represents the predicted occurrence or density in link space at all locations $s$ in time $t$, $f^{-1}()$ is the inverse link function (e.g., logit or log), $\mathbf{X}$ represents a matrix of fixed-effects coefficients (year effects) with estimated coefficients $\mathbf{b}$. We separate the spatial variation $\mathbf{\text{ω}}\sim MVN\left(0,{\mathrm{\Sigma }}_{\omega }\right)$ from the year-to-year spatiotemporal variation ${\mathrm{ϵ}}_t$, where the spatial component represents a spatial intercept (treated as a Gaussian Markov random field) and the spatiotemporal component represents spatially correlated temporal deviations from $\mathbf{\text{ω}}$. Because of the even-odd year sampling and checkerboard pattern in the Canadian survey data, the spatiotemporal fields were modeled as an AR(1) process. Year effects were included as predictors (factors) in both model components, and did not include other predictors. However, exploratory versions of the model that included survey catchability effects resulted in models not converging. Future work could consider gear-level catchability effects (e.g., Davidson et al., 2025) or priors on catchability effects (e.g., Monnahan et al., 2021). Convergence was assessed by assuring that the Hessian matrix was positive definite and all absolute log likelihood gradients were <0.001. Residuals were checked via the DHARMa package (Hartig, 2024).

After fitting a spatiotemporal model, the second step in geostatistical index standardization is making predictions to a gridded surface and computing a weighted sum of predicted biomass across grid cells by year (Thorson et al., 2015). Because our model fitting process was restricted to surveys from British Columbia and the West Coast of the U.S., we used the combined survey grids from these regions for the prediction grid. The “epsilon” bias correction in sdmTMB was implemented, which accounts for non-linear transformations when dealing with random effects (Thorson & Kristensen, 2016). Estimated hake biomass time series indicates that the population declined over the 2003–2015 period, and has since increased slightly (Fig. 2). In contrast, the acoustic survey over the same period is variable and without any clear trends. The same prediction grid used to generate coastwide biomass indices may be further subdivided to quantify regional trends. As an example, we generated three regional indices, representing biomass in Canadian waters and biomass north and south of 42° latitude (this breakpoint corresponds to the northern border of California, but also represents a biogeographic break in the California Current; Sivasundar & Palumbi, 2010). Comparing recent trends in these indices suggests that over the last 5 years, biomass in Canada has remained relatively unchanged while biomass in the U.S. has increased slightly (Fig. 2).

Case study 2

In our second case study, we used data from the surveyjoin package to characterize variation in groundfish community structure among species with overlapping ranges. The diversity and composition of groundfish communities are known to be structured by oceanography (e.g., currents), bathymetry (e.g., shelf width), environmental conditions (e.g., temperature), habitat types (e.g., substrate), and fishing pressure (Howard et al., 2021). Understanding how groundfish communities vary spatially can add to our understanding of these dynamics, as well as the resulting influences on food web dynamics and connectivity of stocks among management regions (e.g., across national borders). The ability to describe groundfish communities across this broad spatial range also facilitates the creation of indicators to rapidly detect community shifts in the future, providing an early warning tool for fishery managers and members of the fishing industry to adapt to environmental change (Thorson, Pinsky & Ward, 2016; Fredston et al., 2021; Thompson et al., 2023).

The surveyjoin package was first used to query survey data from 2003 to 2023 for a subset of groundfish with similar ranges across the northeast Pacific that are well-sampled by the bottom trawl surveys. Species included arrowtooth flounder (Atheresthes stomias), Pacific halibut (Hippoglossus stenolepis), flathead sole (Hippoglossoides elassodon), Pacific cod (Gadus macrocephalus), rex sole (Glyptocephalus zachirus), Pacific ocean perch (Sebastes alutus), sablefish (Anoplopoma fimbria), Dover sole (Microstomus pacificus), English sole (Parophrys vetulus), and lingcod (Ophiodon elongatus).

To understand the variation in groundfish community structure, we constructed a joint (multispecies) species distribution model that included an intercept and an independent spatial random field for each species (with shared range or decorrelation parameter, controlling the distance at which two points are functionally independent). We used an isotropic Matérn correlation function defined by spatial coordinates in latitude-longitude, using a 0.5 degree cutoff to define the SPDE mesh. A Tweedie distribution (log link) was used to represent observation error for each species. Model fitting was done using the tinyVAST package (Thorson et al., 2025). After evaluating convergence and diagnostics, as in Case Study 1, we used the fitted model to predict population density for each species across the survey domain.

Using these predictions of population density, we performed a cluster analysis to summarize variation in species’ biomass across the domain and illustrate how community structure varies over space. Ward’s dissimilarity for log-density between each pair of locations (fastcluster R package; Müllner, 2013) was calculated, and then hierarchical clustering was used to identify clusters. Using a value of K = 3 clusters, we then calculated the average log-density for each species.

Cluster analyses identified breakpoints between the three groundfish community components across the domain: (1) near Cape Mendocino, California, in the south and (2) the intersection of the Alaska peninsula and the Aleutian Islands in the north (Fig. 3). The southern community cluster was characterized by relatively higher biomass estimates of Dover sole, sablefish, rex sole, English sole, and lingcod. The central community ranging from Washington state north through British Columbia and the Gulf of Alaska is characterized by relatively higher biomass estimates of arrowtooth flounder, Pacific halibut, sablefish, and Pacific cod. The third community, located in the Aleutian Islands, had relatively equal representation across species, with slightly lower abundances of the southern species (e.g., English sole, lingcod). Flatfish species clustered together by similar life-history traits, including small-mouthed flatfishes (English sole, rex sole, and Dover sole) and larger piscivorous flatfish (arrowtooth flounder and Pacific halibut). Some of the patterns among species can be explained by ontogenetic habitat shifts, such as sablefish moving to the edge of the continental shelf down the continental slope as they mature, thus becoming less available to the regional surveys that are restricted to the shelf (e.g., in the Aleutian Islands and western Gulf of Alaska).

Case study 3

For our third case study, we demonstrated how trawl survey datasets can be input into a coastwide model to quantify changes in time and space at fine spatial scales. This case study focused on sablefish, which is one of the most commercially important species in the northeast Pacific and for which there has been recent interest in integrating science and management coastwide (Kapur et al., 2024). Models with spatial covariates or trends are becoming increasingly used in fisheries and ecology, and are also referred to as spatially varying coefficient (SVC) models (Hastie & Tibshirani, 1993; Barnett, Ward & Anderson, 2021; Thorson et al., 2023). Models with SVCs can be thought of as an extension of simple linear regression, where the estimated field $\mathit{\omega }$ represents an estimated spatial intercept and the field $\mathbf{\text{ς}}$ represents a spatially varying slope; for this case study, we are interested in long term temporal trends. We can extend the notation in Case Study 1 to represent this model as

$${\mathit{u}}_t=f^{-1}\left(\mathbf{X}\mathbf{b}+\mathit{\omega }+{\mathbf{X}}^{SVC}\mathbf{\text{ς}}\right)$$

We use ${\mathbf{X}}^{SVC}$ as a design matrix for variables associated with the spatially varying coefficient (here ${\mathbf{X}}^{SVC}$ just includes a continuous vector corresponding to years). When the covariate of interest is time, values of the intercept field $\mathit{\omega }$ that are below average correspond to fine scale locations with consistently below average biomass density, and similarly, values of the slope field $\mathbf{\text{ς}}$ that are below average can be used to identify fine scale locations where average annual trends in biomass density are declining. Such approaches may be particularly useful in situations where population trends are patchy or exhibit spatial gradients within a large domain (Davidson et al., 2025).

To construct the coastwide model of sablefish, we combined trawl survey data (2003–2023) from three regions in the surveyjoin package: the West Coast of the U.S., British Columbia, and the Gulf of Alaska. We used a cutoff distance of 25 km to construct the SPDE mesh (n = 944 vertices). Total CPUE was modeled using a Poisson-link delta-lognormal family (Thorson, 2018). We followed the same approaches in previous case studies to evaluate model convergence and diagnostics.

We then used the fitted model to predict on to the grids associated with each survey. Since a Poisson-link function was used, the main and spatial effects from each linear predictor were on the same scale (log) and could be summed into a single intercept field. The main and spatially varying coefficient effects of year were similarly summed into a single slope field. Estimates from the coastwide model of sablefish demonstrate that the West Coast of the U.S. and offshore regions in the Gulf of Alaska appear to consistently have above average biomass density for sablefish, while waters in British Columbia appear slightly below average (Fig. 4). With the exception of southern California, estimates of biomass were trending positively for nearly all of the coastline (strongest positive trends occurred in British Columbia and the Gulf of Alaska; Fig. 4).