17 January 2022

This is the logical riff off of what Jon helped me with to perform some raster math on the bycatch and fishing effort data.

This code generates the CPUE index map in SI Figure S2.

library(raster)
library(ggplot2)
library(dplyr)
library(sp)
#library(leaflet)
library(lubridate)
#library(mapview)
#library(leafsync)
#library(htmlwidgets)
library(RColorBrewer)
library(mapdata)
Loading required package: maps

Attaching package: ‘maps’

The following object is masked from ‘package:purrr’:

    map
fishing <- readRDS("extdata/longline_effort_halfdegree_bins_diana_w_hooks.RDS") %>%
  mutate(season = ifelse(month %in% c(5:9), "breeding", "nonbreeding"))

fish_by_hooks <- fishing %>%
  group_by(season, lonr, latr) %>%
  summarise(sum(hooks)) %>%
  rename(tot_hooks = `sum(hooks)`)
`summarise()` has grouped output by 'season', 'lonr'. You can override using the `.groups` argument.
### read in the effort raster I made
fishing_effort<-raster("extdata/effort.tif")
### read in bycatches
fulmars <- read.csv("bycatch_data/FINAL_nofu_bycatch_downsampled_revised_06122020_newseasons.csv") %>%
  mutate(bycatches=1) %>% #need this for summing points per grid cell later
  filter(!is.na(FisheriesLatDD), !is.na(FisheriesLongDD)) #get rid of NA lat/longs
coordinates(fish_by_hooks) <- ~lonr+latr
proj4string(fish_by_hooks) <- CRS("+proj=longlat +datum=WGS84")
fish_by_hooks_projected <- spTransform(fish_by_hooks, CRS("+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))
### alternative raster
rast <- raster(resolution=1, ext=extent(fish_by_hooks), crs=proj4string(fish_by_hooks)) #create
fishing_effort2 <- rasterize(fish_by_hooks, rast, 'tot_hooks', fun=sum)

# note: this is not in the aea projection.
### make bycatch spatial points, rasterize to match effort raster and sum catches per grid cell
coordinates(fulmars) <- ~FisheriesLongDD+FisheriesLatDD
proj4string(fulmars) <- CRS("+proj=longlat +datum=WGS84")
fulmars_projected <- spTransform(fulmars, CRS("+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

fulmar_catches <- rasterize(fulmars, fishing_effort, 'bycatches', fun=sum, background=0) #use background=0 to make sure you get zeros (instead of NAs) later on when you do raster math since you want the answer "0" not "NA" when there was no bycatch but there was fishing effort
fulmar_catches2 <- fulmar_catches #this one is just for plotting not for math
fulmar_catches2[fulmar_catches2 == 0] <- NA #convert zeroes to NAs just for plotting not for math
# Multiply out the number of bycatches to expand from 2% (for which we have samples) to 100%
fulmar_CPUE <- 1000000*((100*fulmar_catches)/2)/fishing_effort #calculate bycatch per unit effort per million hooks

# convert raster to data frame
test_spdf <- rasterToPoints(fulmar_CPUE)
test_df <- as.data.frame(test_spdf)
colnames(test_df) <- c("x", "y", "value")

# fix the longitude values
cpue.df <- test_df %>%
  mutate(lon = ifelse(x > 0, (360-x)*-1, x)) %>%
  mutate(lat = y)
# set limits
xmax = max(cpue.df$lon)
xmin = min(cpue.df$lon)
ymax = max(cpue.df$lat)
ymin = min(cpue.df$lat)

ylim = c(ymin, ymax)
xlim = c(xmin, xmax)
# get regional polygons
reg = map_data("world2Hires")
regions = subset(reg, region %in% c('USA', 'Canada')) %>%
  filter(lat > 47 & lat < 65, long < 240)

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

# set map limits
#lons = c(-188, -140) #this needs to extend beyond the date line
#lats = c(50, 61.75)

region.zoom <- regions %>%
  filter(long < -145, long > -181, lat > 50, lat < 62)
ggplot() +  
  geom_polygon(data = region.zoom, aes(x = long, y = lat, group = group), 
               fill = "grey38", color = NA) +
  geom_tile(data=cpue.df, aes(x=lon, y=lat, fill=log(value)), alpha=0.7) + 
  #geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
               #fill=NA, color="grey50", size=0.25) +
  xlim(xlim) +
  ylim(ylim) +
  coord_map(projection = "albers", par=55,62) +
  theme_minimal() +
  scale_fill_gradientn(colors = c("lemonchiffon", "gold","orange", "firebrick")) +
  labs(
    y = "Latitude",
    x = "Longitude",
    fill = "CPUE index"
  ) +
  theme(
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10))
  )

ggsave("pdf_outputs/nofu_bycatch_cpue_scaled.pdf")
Saving 5.98 x 3.7 in image

LS0tCnRpdGxlOiAiZmlzaGluZyBlZmZvcnQgdnMuIGJ5Y2F0Y2ggLSBDUFVFIGluZGV4IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgoxNyBKYW51YXJ5IDIwMjIKClRoaXMgaXMgdGhlIGxvZ2ljYWwgcmlmZiBvZmYgb2Ygd2hhdCBKb24gaGVscGVkIG1lIHdpdGggdG8gcGVyZm9ybSBzb21lIHJhc3RlciBtYXRoIG9uIHRoZSBieWNhdGNoIGFuZCBmaXNoaW5nIGVmZm9ydCBkYXRhLgoKVGhpcyBjb2RlIGdlbmVyYXRlcyB0aGUgQ1BVRSBpbmRleCBtYXAgaW4gU0kgRmlndXJlIFMyLgoKYGBge3IgbG9hZC1wYWNrYWdlc30KbGlicmFyeShyYXN0ZXIpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShkcGx5cikKbGlicmFyeShzcCkKbGlicmFyeShsdWJyaWRhdGUpCmxpYnJhcnkoUkNvbG9yQnJld2VyKQpsaWJyYXJ5KG1hcGRhdGEpCmBgYAoKCmBgYHtyIGxvYWQtZmlzaGluZy1lZmZvcnQtZGF0YX0KZmlzaGluZyA8LSByZWFkUkRTKCJleHRkYXRhL2xvbmdsaW5lX2VmZm9ydF9oYWxmZGVncmVlX2JpbnNfZGlhbmFfd19ob29rcy5SRFMiKSAlPiUKICBtdXRhdGUoc2Vhc29uID0gaWZlbHNlKG1vbnRoICVpbiUgYyg1OjkpLCAiYnJlZWRpbmciLCAibm9uYnJlZWRpbmciKSkKCmZpc2hfYnlfaG9va3MgPC0gZmlzaGluZyAlPiUKICBncm91cF9ieShzZWFzb24sIGxvbnIsIGxhdHIpICU+JQogIHN1bW1hcmlzZShzdW0oaG9va3MpKSAlPiUKICByZW5hbWUodG90X2hvb2tzID0gYHN1bShob29rcylgKQoKYGBgCgpgYGB7ciBlZmZvcnQtcmFzdGVyfQojIyMgcmVhZCBpbiB0aGUgZWZmb3J0IHJhc3RlciBJIG1hZGUKZmlzaGluZ19lZmZvcnQ8LXJhc3RlcigiZXh0ZGF0YS9lZmZvcnQudGlmIikKYGBgCgoKYGBge3IgbG9hZC1ieWNhdGNoLWRhdGF9CiMjIyByZWFkIGluIGJ5Y2F0Y2hlcwpmdWxtYXJzIDwtIHJlYWQuY3N2KCJieWNhdGNoX2RhdGEvRklOQUxfbm9mdV9ieWNhdGNoX2Rvd25zYW1wbGVkX3JldmlzZWRfMDYxMjIwMjBfbmV3c2Vhc29ucy5jc3YiKSAlPiUKICBtdXRhdGUoYnljYXRjaGVzPTEpICU+JSAjbmVlZCB0aGlzIGZvciBzdW1taW5nIHBvaW50cyBwZXIgZ3JpZCBjZWxsIGxhdGVyCiAgZmlsdGVyKCFpcy5uYShGaXNoZXJpZXNMYXRERCksICFpcy5uYShGaXNoZXJpZXNMb25nREQpKSAjZ2V0IHJpZCBvZiBOQSBsYXQvbG9uZ3MKCmBgYAoKCmBgYHtyIHByb2plY3QtZmlzaGluZy1lZmZvcnR9CmNvb3JkaW5hdGVzKGZpc2hfYnlfaG9va3MpIDwtIH5sb25yK2xhdHIKcHJvajRzdHJpbmcoZmlzaF9ieV9ob29rcykgPC0gQ1JTKCIrcHJvaj1sb25nbGF0ICtkYXR1bT1XR1M4NCIpCmZpc2hfYnlfaG9va3NfcHJvamVjdGVkIDwtIHNwVHJhbnNmb3JtKGZpc2hfYnlfaG9va3MsIENSUygiK3Byb2o9YWVhICtsYXRfMT01NSArbGF0XzI9NjUgK2xhdF8wPTUwICtsb25fMD0tMTU0ICt4XzA9MCAreV8wPTAgK2VsbHBzPUdSUzgwICtkYXR1bT1OQUQ4MyArdW5pdHM9bSArbm9fZGVmcyIpKQoKYGBgCgpgYGB7cn0KIyMjIGFsdGVybmF0aXZlIHJhc3RlcgpyYXN0IDwtIHJhc3RlcihyZXNvbHV0aW9uPTEsIGV4dD1leHRlbnQoZmlzaF9ieV9ob29rcyksIGNycz1wcm9qNHN0cmluZyhmaXNoX2J5X2hvb2tzKSkgI2NyZWF0ZQpmaXNoaW5nX2VmZm9ydDIgPC0gcmFzdGVyaXplKGZpc2hfYnlfaG9va3MsIHJhc3QsICd0b3RfaG9va3MnLCBmdW49c3VtKQoKIyBub3RlOiB0aGlzIGlzIG5vdCBpbiB0aGUgYWVhIHByb2plY3Rpb24uCmBgYAoKCmBgYHtyfQojIyMgbWFrZSBieWNhdGNoIHNwYXRpYWwgcG9pbnRzLCByYXN0ZXJpemUgdG8gbWF0Y2ggZWZmb3J0IHJhc3RlciBhbmQgc3VtIGNhdGNoZXMgcGVyIGdyaWQgY2VsbApjb29yZGluYXRlcyhmdWxtYXJzKSA8LSB+RmlzaGVyaWVzTG9uZ0REK0Zpc2hlcmllc0xhdERECnByb2o0c3RyaW5nKGZ1bG1hcnMpIDwtIENSUygiK3Byb2o9bG9uZ2xhdCArZGF0dW09V0dTODQiKQpmdWxtYXJzX3Byb2plY3RlZCA8LSBzcFRyYW5zZm9ybShmdWxtYXJzLCBDUlMoIitwcm9qPWFlYSArbGF0XzE9NTUgK2xhdF8yPTY1ICtsYXRfMD01MCArbG9uXzA9LTE1NCAreF8wPTAgK3lfMD0wICtlbGxwcz1HUlM4MCArZGF0dW09TkFEODMgK3VuaXRzPW0gK25vX2RlZnMiKSkKCmZ1bG1hcl9jYXRjaGVzIDwtIHJhc3Rlcml6ZShmdWxtYXJzLCBmaXNoaW5nX2VmZm9ydCwgJ2J5Y2F0Y2hlcycsIGZ1bj1zdW0sIGJhY2tncm91bmQ9MCkgI3VzZSBiYWNrZ3JvdW5kPTAgdG8gbWFrZSBzdXJlIHlvdSBnZXQgemVyb3MgKGluc3RlYWQgb2YgTkFzKSBsYXRlciBvbiB3aGVuIHlvdSBkbyByYXN0ZXIgbWF0aCBzaW5jZSB5b3Ugd2FudCB0aGUgYW5zd2VyICIwIiBub3QgIk5BIiB3aGVuIHRoZXJlIHdhcyBubyBieWNhdGNoIGJ1dCB0aGVyZSB3YXMgZmlzaGluZyBlZmZvcnQKZnVsbWFyX2NhdGNoZXMyIDwtIGZ1bG1hcl9jYXRjaGVzICN0aGlzIG9uZSBpcyBqdXN0IGZvciBwbG90dGluZyBub3QgZm9yIG1hdGgKZnVsbWFyX2NhdGNoZXMyW2Z1bG1hcl9jYXRjaGVzMiA9PSAwXSA8LSBOQSAjY29udmVydCB6ZXJvZXMgdG8gTkFzIGp1c3QgZm9yIHBsb3R0aW5nIG5vdCBmb3IgbWF0aApgYGAKCgpgYGB7ciBmdWxtYXItY3B1ZX0KIyBNdWx0aXBseSBvdXQgdGhlIG51bWJlciBvZiBieWNhdGNoZXMgdG8gZXhwYW5kIGZyb20gMiUgKGZvciB3aGljaCB3ZSBoYXZlIHNhbXBsZXMpIHRvIDEwMCUKZnVsbWFyX0NQVUUgPC0gMTAwMDAwMCooKDEwMCpmdWxtYXJfY2F0Y2hlcykvMikvZmlzaGluZ19lZmZvcnQgI2NhbGN1bGF0ZSBieWNhdGNoIHBlciB1bml0IGVmZm9ydCBwZXIgbWlsbGlvbiBob29rcwoKIyBjb252ZXJ0IHJhc3RlciB0byBkYXRhIGZyYW1lCnRlc3Rfc3BkZiA8LSByYXN0ZXJUb1BvaW50cyhmdWxtYXJfQ1BVRSkKdGVzdF9kZiA8LSBhcy5kYXRhLmZyYW1lKHRlc3Rfc3BkZikKY29sbmFtZXModGVzdF9kZikgPC0gYygieCIsICJ5IiwgInZhbHVlIikKCiMgZml4IHRoZSBsb25naXR1ZGUgdmFsdWVzCmNwdWUuZGYgPC0gdGVzdF9kZiAlPiUKICBtdXRhdGUobG9uID0gaWZlbHNlKHggPiAwLCAoMzYwLXgpKi0xLCB4KSkgJT4lCiAgbXV0YXRlKGxhdCA9IHkpCgpgYGAKCmBgYHtyfQojIHNldCBsaW1pdHMKeG1heCA9IG1heChjcHVlLmRmJGxvbikKeG1pbiA9IG1pbihjcHVlLmRmJGxvbikKeW1heCA9IG1heChjcHVlLmRmJGxhdCkKeW1pbiA9IG1pbihjcHVlLmRmJGxhdCkKCnlsaW0gPSBjKHltaW4sIHltYXgpCnhsaW0gPSBjKHhtaW4sIHhtYXgpCmBgYAoKCgpgYGB7ciBnZXQtY29hc3RsaW5lLXBvbHlnb259CiMgZ2V0IHJlZ2lvbmFsIHBvbHlnb25zCnJlZyA9IG1hcF9kYXRhKCJ3b3JsZDJIaXJlcyIpCnJlZ2lvbnMgPSBzdWJzZXQocmVnLCByZWdpb24gJWluJSBjKCdVU0EnLCAnQ2FuYWRhJykpICU+JQogIGZpbHRlcihsYXQgPiA0NyAmIGxhdCA8IDY1LCBsb25nIDwgMjQwKQoKIyBjb252ZXJ0IGxhdCBsb25ncwpyZWdpb25zJGxvbmcgPSAoMzYwIC0gcmVnaW9ucyRsb25nKSotMQoKIyBzZXQgbWFwIGxpbWl0cwojbG9ucyA9IGMoLTE4OCwgLTE0MCkgI3RoaXMgbmVlZHMgdG8gZXh0ZW5kIGJleW9uZCB0aGUgZGF0ZSBsaW5lCiNsYXRzID0gYyg1MCwgNjEuNzUpCgpyZWdpb24uem9vbSA8LSByZWdpb25zICU+JQogIGZpbHRlcihsb25nIDwgLTE0NSwgbG9uZyA+IC0xODEsIGxhdCA+IDUwLCBsYXQgPCA2MikKYGBgCgoKYGBge3IgcGxvdC1jcHVlLWluZGV4fQojIG1ha2UgcGxvdCBvZiBieWNhdGNoIHBlciBmaXNoaW5nIGVmZm9ydCAKZ2dwbG90KCkgKyAgCiAgZ2VvbV9wb2x5Z29uKGRhdGEgPSByZWdpb24uem9vbSwgYWVzKHggPSBsb25nLCB5ID0gbGF0LCBncm91cCA9IGdyb3VwKSwgCiAgICAgICAgICAgICAgIGZpbGwgPSAiZ3JleTM4IiwgY29sb3IgPSBOQSkgKwogIGdlb21fdGlsZShkYXRhPWNwdWUuZGYsIGFlcyh4PWxvbiwgeT1sYXQsIGZpbGw9bG9nKHZhbHVlKSksIGFscGhhPTAuNykgKyAKICAjZ2VvbV9wb2x5Z29uKGRhdGE9T1IsIGFlcyh4PWxvbmcsIHk9bGF0LCBncm91cD1ncm91cCksIAogICAgICAgICAgICAgICAjZmlsbD1OQSwgY29sb3I9ImdyZXk1MCIsIHNpemU9MC4yNSkgKwogIHhsaW0oeGxpbSkgKwogIHlsaW0oeWxpbSkgKwogIGNvb3JkX21hcChwcm9qZWN0aW9uID0gImFsYmVycyIsIHBhcj01NSw2MikgKwogIHRoZW1lX21pbmltYWwoKSArCiAgc2NhbGVfZmlsbF9ncmFkaWVudG4oY29sb3JzID0gYygibGVtb25jaGlmZm9uIiwgImdvbGQiLCJvcmFuZ2UiLCAiZmlyZWJyaWNrIikpICsKICBsYWJzKAogICAgeSA9ICJMYXRpdHVkZSIsCiAgICB4ID0gIkxvbmdpdHVkZSIsCiAgICBmaWxsID0gIkNQVUUgaW5kZXgiCiAgKSArCiAgdGhlbWUoCiAgICBheGlzLnRpdGxlLnggPSBlbGVtZW50X3RleHQobWFyZ2luID0gbWFyZ2luKHQgPSAxMCkpLAogICAgYXhpcy50aXRsZS55ID0gZWxlbWVudF90ZXh0KG1hcmdpbiA9IG1hcmdpbihyID0gMTApKQogICkKCiMgTWFrZSBTSSBGaWd1cmUgMiB3aXRoIENQVUUgaW5kZXgKIyBtaXNtYXRjaCBiZXR3ZWVuIGZpc2hpbmcgZWZmb3J0IGFuZCBieWNhdGNoICUgbWFrZXMgdGhlIG51bWVyaWMgc2NhbGUgbWVhbmluZ2xlc3MsIGJ1dCB0aGUgcmVsYXRpdmUgYW1vdW50cyBjb3JyZWN0Cmdnc2F2ZSgicGRmX291dHB1dHMvbm9mdV9ieWNhdGNoX2NwdWVfc2NhbGVkLnBkZiIpCmBgYAoKCg==