Are there differences between the spatial distributions of bycatch from the different colonies during the breeding and nonbreeding seasons? I’ll use BA values to compare among the pairwise comparisons.
library(tidyverse)
[30m── [1mAttaching packages[22m ─────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──[39m
[30m[32m✓[30m [34mggplot2[30m 3.3.2 [32m✓[30m [34mpurrr [30m 0.3.3
[32m✓[30m [34mtibble [30m 3.0.1 [32m✓[30m [34mdplyr [30m 1.0.4
[32m✓[30m [34mtidyr [30m 1.1.2 [32m✓[30m [34mstringr[30m 1.4.0
[32m✓[30m [34mreadr [30m 1.3.1 [32m✓[30m [34mforcats[30m 0.4.0[39m
package ‘ggplot2’ was built under R version 3.6.2package ‘tibble’ was built under R version 3.6.2package ‘tidyr’ was built under R version 3.6.2[30m── [1mConflicts[22m ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31mx[30m [34mtidyr[30m::[32mextract()[30m masks [34mraster[30m::extract()
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mpurrr[30m::[32mis_null()[30m masks [34mtestthat[30m::is_null()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()
[31mx[30m [34mdplyr[30m::[32mmatches()[30m masks [34mtidyr[30m::matches(), [34mtestthat[30m::matches()
[31mx[30m [34mdplyr[30m::[32mselect()[30m masks [34mraster[30m::select()[39m
library(dplyr)
library(adehabitatHR)
Loading required package: deldir
deldir 0.1-25
Loading required package: ade4
Loading required package: adehabitatMA
Registered S3 methods overwritten by 'adehabitatMA':
method from
print.SpatialPixelsDataFrame sp
print.SpatialPixels sp
Attaching package: ‘adehabitatMA’
The following object is masked from ‘package:raster’:
buffer
Loading required package: adehabitatLT
Loading required package: CircStats
Loading required package: MASS
Attaching package: ‘MASS’
The following object is masked from ‘package:dplyr’:
select
The following objects are masked from ‘package:raster’:
area, select
Loading required package: boot
Attaching package: ‘adehabitatLT’
The following object is masked from ‘package:dplyr’:
id
library(ggplot2)
library(magrittr)
Attaching package: ‘magrittr’
The following object is masked from ‘package:purrr’:
set_names
The following object is masked from ‘package:tidyr’:
extract
The following objects are masked from ‘package:testthat’:
equals, is_less_than, not
The following object is masked from ‘package:raster’:
extract
library(sp)
library(raster)
library(marmap)
Attaching package: ‘marmap’
The following object is masked from ‘package:oce’:
plotProfile
The following object is masked from ‘package:raster’:
as.raster
The following object is masked from ‘package:grDevices’:
as.raster
# assignments
nofu_w_meta <- read_csv("bycatch_data/FINAL_nofu_bycatch_downsampled_revised_06122020_newseasons.csv")
Parsed with column specification:
cols(
SAMPLE_ID = [31mcol_character()[39m,
NMFS_DNA_ID = [31mcol_character()[39m,
collection = [31mcol_character()[39m,
PofZ = [32mcol_double()[39m,
z_score = [32mcol_double()[39m,
TemporalAvail = [32mcol_double()[39m,
SpatialAvail = [32mcol_double()[39m,
FisheriesCollMonth = [32mcol_double()[39m,
FisheriesCollSeason = [31mcol_character()[39m,
FisheriesLatDD = [32mcol_double()[39m,
FisheriesLongDD = [32mcol_double()[39m,
MissingSBND = [32mcol_double()[39m
)
# summarise those data a bit...
nofu_w_meta %>%
filter(!is.na(FisheriesLatDD)) %>% # remove entries that don't have Lat/Lon info
filter(!is.na(FisheriesCollSeason)) %>% # remove entries that don't have temporal info
#group_by(collection, FisheriesCollSeason) %>% # how many birds per colony/season
#group_by(FisheriesCollSeason) %>%
tally()
# are all data >90% likelihood?
nofu_w_meta %>%
filter(PofZ > 0.90) %>%
dplyr::select(SAMPLE_ID) %>%
unique()
# colony locations
colony_locs <- read_csv("bycatch_data/colony_lat_lon.csv")
Parsed with column specification:
cols(
Colony = [31mcol_character()[39m,
lon = [32mcol_double()[39m,
lat = [32mcol_double()[39m
)
# how many during which seasons?
nofu_w_meta %>%
group_by(FisheriesCollSeason) %>%
tally()
# arrange assignment and fishery data for spatial analysis
new_spatial_df <- nofu_w_meta %>%
#select(-TAG_NUMBER, -NMFS_DNA_ID, -OriginColony, -ColonyLatDD, -ColonyLongDD, -PofZ, -z_score) %>%
dplyr::select(NMFS_DNA_ID, collection, FisheriesCollSeason, FisheriesLatDD, FisheriesLongDD) %>%
rename(OriginColony = collection) %>%
filter(!is.na(FisheriesLatDD), FisheriesLatDD>50) %>% #get rid of blanks, keep only alaska data
filter(!is.na(FisheriesCollSeason)) %>% #get rid of blanks
mutate(ColonyPhenology = paste(OriginColony, FisheriesCollSeason, sep="_")) #make colony-phenology grouping fieldg
NOFU <- new_spatial_df %>%
dplyr::select(OriginColony, FisheriesLongDD, FisheriesLatDD, ColonyPhenology) #keep only these fields for kerneling
972 birds at this stage (29 birds were caught below 50 degrees N latitude)
plot(NOFU$FisheriesLongDD, NOFU$FisheriesLatDD, asp = 1)
Because of the date line, these birds at +176 are just on the other side of 180 degrees…
NOFU %>%
filter(FisheriesLongDD > 150) #%>%
# group_by(ColonyPhenology) %>%
# tally() %>%
# mutate(prop = n/sum(n))
I don’t need to mess with the longitude values in the new projection.
nofu_dat_fixed <- NOFU %>%
#mutate(long = ifelse(FisheriesLongDD > 150, -FisheriesLongDD, FisheriesLongDD)) %>% # 36 birds like this
mutate(lat = FisheriesLatDD) %>%
mutate(long = FisheriesLongDD)
# and save a version of that df that is still in df format
nofu_data_frame <- NOFU %>%
#mutate(long = ifelse(FisheriesLongDD > 150, -FisheriesLongDD, FisheriesLongDD)) %>%
mutate(lat = FisheriesLatDD) %>%
mutate(long = FisheriesLongDD)
# plot again
plot(nofu_dat_fixed$long, nofu_dat_fixed$lat, asp = 1)
NA
NA
### define projections
input_proj <- "+proj=longlat +datum=WGS84" ###the input projection of track points
desired_proj <- "+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" ###the desired projection for track points
### convert data to a spatial object and reproject to desired coordinate system
coordinates(nofu_dat_fixed) = ~long+lat ### define the x and y spatial fields
proj4string(nofu_dat_fixed) <- CRS(input_proj) ### define the projection
nofu_dat2 <- spTransform(nofu_dat_fixed, CRS(desired_proj)) ###transform to desired projection
# another quick plot
plot(nofu_dat2, pch = 19, cex = 0.5)
# I need this in df form again - not spatial df
NOFU_df2 <- new_spatial_df %>%
#mutate(long = ifelse(FisheriesLongDD > 150, -FisheriesLongDD, FisheriesLongDD)) %>%
mutate(lat = FisheriesLatDD) %>%
mutate(long = FisheriesLongDD)
Add bathymetric map
# antimeridian region
aleu <- getNOAA.bathy(175, -145, 50, 65, resolution = 4,
antimeridian = TRUE)
Querying NOAA database ...
This may take seconds to minutes, depending on grid size
Building bathy matrix ...
# Make it a raster
bathy <- as.raster(aleu)
# Create a xyz table for ggplot
bath<-fortify(aleu)
str(bathy)
Formal class 'RasterLayer' [package "raster"] with 12 slots
..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots
.. .. ..@ name : chr ""
.. .. ..@ datanotation: chr "FLT4S"
.. .. ..@ byteorder : chr "little"
.. .. ..@ nodatavalue : num -Inf
.. .. ..@ NAchanged : logi FALSE
.. .. ..@ nbands : int 1
.. .. ..@ bandorder : chr "BIL"
.. .. ..@ offset : int 0
.. .. ..@ toptobottom : logi TRUE
.. .. ..@ blockrows : int 0
.. .. ..@ blockcols : int 0
.. .. ..@ driver : chr ""
.. .. ..@ open : logi FALSE
..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots
.. .. ..@ values : int [1:134775] 44 58 41 28 22 16 29 31 31 30 ...
.. .. ..@ offset : num 0
.. .. ..@ gain : num 1
.. .. ..@ inmemory : logi TRUE
.. .. ..@ fromdisk : logi FALSE
.. .. ..@ isfactor : logi FALSE
.. .. ..@ attributes: list()
.. .. ..@ haveminmax: logi TRUE
.. .. ..@ min : int -7380
.. .. ..@ max : int 4537
.. .. ..@ band : int 1
.. .. ..@ unit : chr ""
.. .. ..@ names : chr ""
..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots
.. .. ..@ type : chr(0)
.. .. ..@ values : logi(0)
.. .. ..@ color : logi(0)
.. .. ..@ names : logi(0)
.. .. ..@ colortable: logi(0)
..@ title : chr(0)
..@ extent :Formal class 'Extent' [package "raster"] with 4 slots
.. .. ..@ xmin: num 175
.. .. ..@ xmax: num 215
.. .. ..@ ymin: num 50
.. .. ..@ ymax: num 65
..@ rotated : logi FALSE
..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
.. .. ..@ geotrans: num(0)
.. .. ..@ transfun:function ()
..@ ncols : int 599
..@ nrows : int 225
..@ crs :Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
..@ history : list()
..@ z : list()
The projection for the bathy raster appears to be: “+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0”
I need the blank grid and the nofu df in the same projection
** Need new blank grid with the appropriate projection:
# new desired projection
new_proj <- "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154"
### convert data to a spatial object and reproject to desired coordinate system
#coordinates(nofu_dat_fixed) = ~long+lat ### define the x and y spatial fields
#proj4string(nofu_dat_fixed) <- CRS(input_proj) ### define the projection
nofu_dat3 <- spTransform(nofu_dat_fixed, CRS(new_proj)) ###transform to desired projection
### build a grid to use for the KUDs
ext<-1.3*extent(nofu_dat3) #define the spatial extent as 1.3X the extent of your data
grid_blank2 <- raster(ext=ext, crs=new_proj, resolution=3000) #make a grid using the extent, define the projection, and set the cell size (resolution)
vals <- 1:ncell(grid_blank2) #give the grid dummy values
grid_blank2 <- setValues(grid_blank2, vals) #give the grid dummy values
grid_blank2 <- as(grid_blank2, "SpatialPixelsDataFrame") #convert grid to Spatial Pixels Data Frame (what kernelUD needs)
this mapping from: https://hansenjohnson.org/post/bathymetric-maps-in-r/
library(oce)
library(ocedata)
package ‘ocedata’ was built under R version 3.6.2
data("coastlineWorldFine")
# convert bathymetry
bathyLon = as.numeric(rownames(aleu))
bathyLat = as.numeric(colnames(aleu))
bathyZ = as.numeric(aleu)
dim(bathyZ) = dim(aleu)
# define plotting region
mlon = mean(NOFU_df2$long)
mlat = mean(NOFU_df2$lat)
span = 2000
col_comparisons_kernel <- function(colony_phenology1, colony_phenology2){
map_dat1 <- NOFU_df2 %>%
filter(ColonyPhenology %in% c(colony_phenology1, colony_phenology2))
# Make a SPDF
spdf <- SpatialPointsDataFrame(coords=cbind(map_dat1$long,
map_dat1$lat),
data=data.frame(id=map_dat1$ColonyPhenology),
proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
# convert my UDs into the right projection to match up with the bathymap
# Project data into laea
spdf.t <- spTransform(spdf,
CRS("+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154"))
# calculate kernelUD
ud <- kernelUD(spdf.t, h = 30000, grid=grid_blank2)
tmp_ud <- getverticeshr(ud, percent=50, standardize=T)
# transform to the same projection as the bathymetry
tmp_ud2 <- spTransform(tmp_ud, CRS = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
# fortify for ggplotting
uds_df <- fortify(tmp_ud2) %>%
mutate(lon = ifelse(long < 0, -long, long))
# colors and formatting
col1 <-uds_df %>%
filter(id == colony_phenology1)
col2 <- uds_df %>%
filter(id == colony_phenology2)
return(list
(col1,
col2))
}
Apply that function to Chagulak and StMatthew, breeding season:
comp1
[[1]]
[[2]]
NA
# and just the locations for st matt and chagulak
matt_chag_locs <- colony_locs %>%
filter(Colony %in% c("Chagulak", "StMatthew"))
matt_chag_locs$Colony <- as.factor(matt_chag_locs$Colony)
colors <- c("purple", "dodgerblue")[matt_chag_locs$Colony]
chag <- comp1[[1]]
stmatt <- comp1[[2]]
Plot it up
# set up the output
pdf("pdf_outputs/chagulak_breeding.pdf", width = 5.5, height = 4.5)
plot(coastlineWorldFine, clon = -168, clat = 56.8, span = 825,
projection="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154", col = 'lightgrey')
# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
levels = c(-500, -1000, -2500),
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
col = 'darkgray')
# add depth legend
legend("bottomleft", seg.len = 3, cex = 0.7,
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
legend = c("-500", "-1000", "-2500"),
box.col = "white", col = 'darkgray', title = "Depth (m)", bg = "white")
# add map data
#mapPolygon(longitude = stmatt$long, latitude = stmatt$lat, lty = 0, col = alpha("dodgerblue", 0.2))
mapPolygon(longitude = chag$long, latitude = chag$lat, lty = 0, col = alpha("purple", 0.2))
# add colony legend
# legend("bottom", seg.len = 3, cex = 0.7,
# legend = c("lightblue", "lightgreen"),
# box.col = "white", col = 'darkgray', title = "Breeding colony", bg = "white")
# add colony locations
mapPoints(longitude = matt_chag_locs$lon, latitude = matt_chag_locs$lat, pch = 23, col = colors, cex = 1.2, lwd = 2)
dev.off()
null device
1
Apply that function to Chagulak and StMatthew, nonbreeding season:
colony_phenology1 <- "Chagulak_nonbreeding"
colony_phenology2 <- "StMatthew_nonbreeding"
comp2 <- col_comparisons_kernel(colony_phenology1, colony_phenology2)
Regions defined for each Polygons
nob.chag <- comp2[[1]]
nob.stmatt <- comp2[[2]]
# set up the output
#plot(coastlineWorldFine, clon = -168, clat = mlat, span = 1100,
#pdf("pdf_outputs/stmatt_nonbreeding.pdf", width = 5.5, height = 4.5)
pdf("pdf_outputs/chagulak_nonbreeding.pdf", width = 5.5, height = 4.5)
plot(coastlineWorldFine, clon = -168, clat = 56.8, span = 825,
projection="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154", col = 'lightgrey')
# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
levels = c(-500, -1000, -2500),
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
col = 'darkgray')
# add depth legend
legend("bottomleft", seg.len = 3, cex = 0.7,
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
legend = c("-500", "-1000", "-2500"),
box.col = "white", col = 'darkgray', title = "Depth (m)", bg = "white")
# add map data
#mapPolygon(longitude = nob.stmatt$long, latitude = nob.stmatt$lat, col = alpha("dodgerblue", 0.2), lty = 0)
mapPolygon(longitude = nob.chag$long, latitude = nob.chag$lat, lty = 0, col = alpha("purple", 0.2))
# add colony locations
mapPoints(longitude = matt_chag_locs$lon, latitude = matt_chag_locs$lat, pch = 23, col = colors, cex = 1.2, lwd = 2)
dev.off()
null device
1
Apply that function to Pribilofs and Semidis, breeding season:
colony_phenology1 <- "Semidi_breeding"
colony_phenology2 <- "Pribilof_breeding"
comp1 <- col_comparisons_kernel(colony_phenology1, colony_phenology2)
Regions defined for each Polygons
# and just the locations for pribs and semidis
semidi_prib_locs <- colony_locs %>%
filter(Colony %in% c("Pribilof", "Semidi"))
semidi_prib_locs$Colony <- as.factor(semidi_prib_locs$Colony)
colors2 <- c("red", "darkorange")[semidi_prib_locs$Colony]
# select-the-output-for-mapping
semidi.b <- comp1[[1]]
prib.b <- comp1[[2]]
Plot it up
Nonbreeding season: Pribs and Semidis
colony_phenology1 <- "Semidi_nonbreeding"
colony_phenology2 <- "Pribilof_nonbreeding"
comp1 <- col_comparisons_kernel(colony_phenology1, colony_phenology2)
Regions defined for each Polygons
# select-the-output-for-mapping
semidi.nob <- comp1[[1]]
prib.nob <- comp1[[2]]
Plot it up
# set up the output
#pdf("pdf_outputs/prib_semidi_nonbreeding_map.pdf", width = 6.5, height = 5.5)
pdf("pdf_outputs/prib_nonbreeding.pdf", width = 5.5, height = 4.5)
plot(coastlineWorldFine, clon = -168, clat = 56.8, span = 825,
#plot(coastlineWorldFine, clon = -168, clat = mlat, span = 1100,
projection="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154", col = 'lightgrey')
# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
levels = c(-500, -1000, -2500),
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
col = 'darkgray')
# add depth legend
legend("bottomleft", seg.len = 3, cex = 0.7,
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
legend = c("-500", "-1000", "-2500"),
box.col = "white", col = 'darkgray', title = "Depth (m)", bg = "white")
# add map data
#mapPolygon(longitude = semidi.nob$long, latitude = semidi.nob$lat, lty = 0, col = alpha("darkorange", 0.2))
mapPolygon(longitude = prib.nob$long, latitude = prib.nob$lat, lty = 0, col = alpha("red", 0.2))
# add colony legend
# legend("bottom", seg.len = 3, cex = 0.7,
# legend = c("lightblue", "lightgreen"),
# box.col = "white", col = 'darkgray', title = "Breeding colony", bg = "white")
# add colony locations
mapPoints(longitude = semidi_prib_locs$lon, latitude = semidi_prib_locs$lat, pch = 23, col = colors2, cex = 1.2, lwd = 2)
dev.off()
null device
1
Next set of comparisons: St. Matthew and the Pribilof Islands
colony_phenology1 <- "Pribilof_breeding"
colony_phenology2 <- "StMatthew_breeding"
comp3 <- col_comparisons_kernel(colony_phenology1, colony_phenology2)
Regions defined for each Polygons
# colors and formatting
prib.b <- comp3[[1]]
matt.b <- comp3[[2]]
# select those colonies
# keep the color consistent for St Matt's
# and just the locations for st matt and chagulak
p.m_locs <- colony_locs %>%
filter(Colony %in% c("Pribilof", "StMatthew"))
p.m_locs$Colony <- as.factor(p.m_locs$Colony)
colors <- c("red", "dodgerblue")[p.m_locs$Colony]
Plot it up
# set up the output
#plot(coastlineWorldFine, clon = -168, clat = mlat, span = 1100,
pdf("pdf_outputs/prib_breeding.pdf", width = 5.5, height = 4.5)
plot(coastlineWorldFine, clon = -168, clat = 56.8, span = 825,
projection="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154", col = 'lightgrey')
# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
levels = c(-500, -1000, -2500),
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
col = 'darkgray')
# add depth legend
legend("bottomleft", seg.len = 3, cex = 0.7,
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
legend = c("-500", "-1000", "-2500"),
box.col = "white", col = 'darkgray', title = "Depth (m)", bg = "white")
# add map data
mapPolygon(longitude = prib.b$long, latitude = prib.b$lat, col = alpha("red", 0.2), lty = 0)
#mapPolygon(longitude = matt.b$long, latitude = matt.b$lat, lty = 0, col = alpha("dodgerblue", 0.2))
# add colony locations
mapPoints(longitude = p.m_locs$lon, latitude = p.m_locs$lat, pch = 23, col = colors, cex = 1.2, lwd = 2)
dev.off()
null device
1
nonbreeding:
colony_phenology1 <- "Pribilof_nonbreeding"
colony_phenology2 <- "StMatthew_nonbreeding"
comp4 <- col_comparisons_kernel(colony_phenology1, colony_phenology2)
Regions defined for each Polygons
# colors and formatting
prib.nob <- comp4[[1]]
matt.nob <- comp4[[2]]
plot(coastlineWorldFine, clon = -168, clat = mlat, span = 1100,
projection="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154", col = 'lightgrey')
# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
levels = c(-500, -1000, -2500),
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
col = 'darkgray')
# add depth legend
legend("bottomleft", seg.len = 3, cex = 0.7,
lwd = c(1.5, 1, 1),
lty = c(3, 1, 3),
legend = c("-500", "-1000", "-2500"),
box.col = "white", col = 'darkgray', title = "Depth (m)", bg = "white")
# add map data
mapPolygon(longitude = prib.nob$long, latitude = prib.nob$lat, col = alpha("red", 0.2), lty = 0)
mapPolygon(longitude = matt.nob$long, latitude = matt.nob$lat, lty = 0, col = alpha("dodgerblue", 0.2))
# add colony locations
mapPoints(longitude = p.m_locs$lon, latitude = p.m_locs$lat, pch = 23, col = colors, cex = 1.2, lwd = 2)
Colony size map Quick and dirty way: list of colony sizes
colony_pops <- colony_locs %>%
mutate(pop = ifelse(Colony == "Chagulak", 500000, 450000)) %>%
mutate(pop = ifelse(Colony == "Pribilof", 79700, pop)) %>%
mutate(pop = ifelse(Colony == "Semidi", 440000, pop))
# add the % bycatch per colony
bycatch_prop <- c(.14, .36, .23, .27) # the order is Chag/Semidi/Prib/StM
by_percent <- as_tibble(bycatch_prop) %>%
rename(by_prop = value)
colony_perc <- colony_pops %>%
mutate(pop_perc = pop/(sum(pop))) %>% # and then make that a percentage
bind_cols(by_percent)
#colony_pops$pop <- factor(colony_pops$pop)
colony_locs$Colony <- as.factor(colony_locs$Colony)
#colors <- c("purple", "red", "goldenrod", "dodgerblue")[colony_locs$Colony]
# # easiest to grab this from Jessie's table
# colony <- c("Chagulak", "Pribilof", "StMatthew", "Semidi")
#
#
# by_prop <- as_data_frame(bycatch_prop) %>%
# rename(by_prop = value)
#
# as_data_frame(colony) %>%
# rename(colony = value) %>%
# bind_cols(by_prop)
library(ggpattern)
Attaching package: ‘ggpattern’
The following objects are masked from ‘package:ggplot2’:
flip_data, flipped_names, gg_dep, has_flipped_aes, remove_missing, should_stop,
waiver
Straight colony map
colors
[1] "purple" "darkorange" "red" "dodgerblue"
Warning messages:
1: package ‘raster’ was built under R version 3.6.2
2: package ‘sp’ was built under R version 3.6.2
#pdf("pdf_outputs/August28_test_disproportionate_map.pdf")
plot(coastlineWorldFine, clon = -167, clat = mlat, span = 1100,
projection="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154", col = 'lightgrey')
# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
levels = c(-500, -1000, -2500),
lwd = c(1, 1, 0.5),
lty = c(3, 1, 3, 1, 3),
col = 'darkgray')
# add depth legend
legend("bottomleft", seg.len = 3, cex = 0.7,
lwd = c(1, 1, 0.5),
lty = c(3, 1, 3, 1, 3),
legend = c("-500", "-1000", "-2500"),
box.col = "white", col = 'darkgray', title = "Depth (m)", bg = "white")
# add map data
# add colony locations
# mapPoints(longitude = colony_locs$lon, latitude = colony_locs$lat, pch = 16, col = alpha(colors, 0.5), cex = 0.00001* (colony_pops$pop), lwd = 2) # relative population sizes
mapPoints(longitude = colony_locs$lon, latitude = colony_locs$lat, pch = 1, col = alpha("dodgerblue",0.7), cex = 11*(colony_perc$pop_perc), lwd = 2) # relative population sizes
mapPoints(longitude = colony_locs$lon, latitude = colony_locs$lat, pch = 1, col = "red", cex = 11*(colony_perc$by_prop), lwd = 2) # relative bycatch amounts
#dev.off()
I need to recalculate the BA values for all of my colony comparisons. Fortunately, I don’t think I need to do so using the p-value randomized framework because none of those p-values were remotely close to statistically significant.
source("R/nofu-functions.R")
build blank grid for use in the kernel overlap function
### define projections
input_proj <- "+proj=longlat +datum=WGS84" ###the input projection of track points
desired_proj <- "+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" ###the desired projection for track points
### convert data to a spatial object and reproject to desired coordinate system
coordinates(NOFU) = ~FisheriesLongDD+FisheriesLatDD ### define the x and y spatial fields
proj4string(NOFU) <- CRS(input_proj) ### define the projection
NOFU2 <- spTransform(NOFU, CRS(desired_proj)) ###transform to desired projection
### build a grid to use for the KUDs
ext<-1.3*extent(NOFU2) #define the spatial extent as 1.3X the extent of your data
grid_blank <- raster(ext=ext, crs=desired_proj, resolution=3000) #make a grid using the extent, define the projection, and set the cell size (resolution)
vals <- 1:ncell(grid_blank) #give the grid dummy values
grid_blank <- setValues(grid_blank, vals) #give the grid dummy values
grid_blank <- as(grid_blank, "SpatialPixelsDataFrame") #convert grid to Spatial Pixels Data Frame (what kernelUD needs)
Basically, the function I want is to enter the colony names and season of interest, and have the rest wrapped up in a function that leans on Abram’s fct.
colony_season_BA <- function(colony1, colony2, season, nofu_spatial_df){
# subset the data by colony comparisons
comparison_df <- nofu_spatial_df %>%
filter(OriginColony == colony1 | OriginColony == colony2) %>%
filter(FisheriesCollSeason == season)
comparison_BA <- kernalOverlapBA_p(tracks = comparison_df, tripid="NMFS_DNA_ID", groupid = "ColonyPhenology",
lon="FisheriesLongDD",lat="FisheriesLatDD",
colonyLon=-169.6449, colonyLat=56.66225,
its=2, h=30000, ud.grid = grid_blank, Plot=F)
return(comparison_BA)
}
test that out for the St Matthew/Prib breeding and nonbreeding
colony_season_BA("StMatthew", "Pribilof", "breeding", new_spatial_df)
[1] "iteration: 1 of 2"
[1] "iteration: 2 of 2"
colony_season_BA("StMatthew", "Pribilof", "nonbreeding", new_spatial_df)
colony_season_BA("StMatthew", "Chagulak", "nonbreeding", new_spatial_df)
colony_season_BA("StMatthew", "Chagulak", "breeding", new_spatial_df)
colony_season_BA("Pribilof", "Chagulak", "breeding", new_spatial_df)
colony_season_BA("Pribilof", "Chagulak", "nonbreeding", new_spatial_df)
colony_season_BA("Semidi", "Chagulak", "breeding", new_spatial_df)
colony_season_BA("Semidi", "Chagulak", "nonbreeding", new_spatial_df)
colony_season_BA("Semidi", "StMatthew", "breeding", new_spatial_df)
colony_season_BA("Semidi", "StMatthew", "nonbreeding", new_spatial_df)
colony_season_BA("Semidi", "Pribilof", "breeding", new_spatial_df)
colony_season_BA("Semidi", "Pribilof", "nonbreeding", new_spatial_df)