diff --git a/.github/workflows/bookdown.yaml b/.github/workflows/bookdown.yaml index 5b17d648..d3cb5282 100644 --- a/.github/workflows/bookdown.yaml +++ b/.github/workflows/bookdown.yaml @@ -4,8 +4,7 @@ on: push: branches: - master - - brandon_url_changes - - andytest1 + - dev pull_request: branches: diff --git a/.gitignore b/.gitignore index 3633393c..0386e3f1 100644 --- a/.gitignore +++ b/.gitignore @@ -35,5 +35,4 @@ rsconnect/ # bookdown files _book/ _bookdown_files/ -packages.bib _book/*.json diff --git a/R/stored_scripts/species_density_analysis.R b/R/stored_scripts/species_density_analysis.R index 40c47024..9b10fc5b 100644 --- a/R/stored_scripts/species_density_analysis.R +++ b/R/stored_scripts/species_density_analysis.R @@ -1,356 +1,356 @@ -# species density analysis - -# -# ```{r, echo = T, message= F, warning=F, include=T, eval = T} -library(ecodata) -library(maps) -library(mapdata) -library(ks) -library(marmap) -library(raster) -library(geosphere) -library(dplyr) - -# plot ks plots - -#----STUFF TO SET----------------- -data.dir <- here::here("data") -gis.dir <- here::here("gis") - -# Gives most recent period. 2015-2019 input is 2016-2018 data -rminyr <- 2015 -rmaxyr <- 2019 - -# tlevel is density of color for KD contours areas -tlevel=75 # move later -color_b="blue" -color_r="orange3" -color_r="tomato3" - -# Code to read in strata and compute areas. Or read from cache. - -# readin in strata.shp and compute areas of strata -# TrawlStrata <- raster::shapefile(file.path(gis.dir,"BTS_Strata.shp")) -# AREA <- geosphere::areaPolygon(TrawlStrata, r=6371000)/10^6 - -# for array of strata and area and make into dataframe -# stratareas <- cbind(TrawlStrata@data$STRATA, AREA) -# colnames(stratareas) <- c("STRATA","AREA") -# stratareas <- data.frame(stratareas) -load(file.path(data.dir, "StratAreas.Rdata")) - -# Query bathymetry or load from cache -# getNOAA.bathy(lon1 = -77, lon2 = -65, lat1 = 35, lat2 = 45, -# resolution = 10) -> nesbath -load(file.path(data.dir, "nesbath.Rdata")) - -#Load raw survey data -load(file.path(data.dir, "Survdat.RData")) - -# MUST run addTrans function -# color transparency -addTrans <- function(color,trans){ - # This function adds transparancy to a color. - # Define transparancy with an integer between 0 and 255 - # 0 being fully transparant and 255 being fully visable - # Works with either color and trans a vector of equal length, - # or one of the two of length 1. - - if (length(color)!=length(trans)&!any(c(length(color), - length(trans))==1)) - stop("Vector lengths not correct") - if (length(color)==1 & length(trans)>1) color <- rep(color,length(trans)) - if (length(trans)==1 & length(color)>1) trans <- rep(trans,length(color)) - - num2hex <- function(x) - { - hex <- unlist(strsplit("0123456789ABCDEF",split="")) - return(paste(hex[(x-x%%16)/16+1],hex[x%%16+1],sep="")) - } - rgb <- rbind(col2rgb(color),trans) - res <- paste("#",apply(apply(rgb,2,num2hex),2,paste,collapse=""),sep="") - return(res) -} - -plot_kd <- function(species, season, exclude_years){ - - # stata to use - # offshore strata to use - CoreOffshoreStrata <- c(seq(1010,1300,10),1340, seq(1360,1400,10),seq(1610,1760,10)) - - # inshore strata to use, still sampled by Bigelow - CoreInshore73to12 <- c(3020, 3050, 3080 ,3110 ,3140 ,3170, 3200, 3230, - 3260, 3290, 3320, 3350 ,3380, 3410 ,3440) - # combine - strata_used <- c(CoreOffshoreStrata,CoreInshore73to12) - - survdat <- survdat %>% - dplyr::select(c(CRUISE6,STATION,STRATUM,SVSPP,YEAR, - SEASON,LAT,LON,ABUNDANCE,BIOMASS)) %>% - filter(SEASON == season, - STRATUM %in% strata_used) %>% # delete record form non-core - #strata and get unique records, - # should be one per species - distinct() %>% - # add field with rounded BIOMASS scaler used to adjust distributions - mutate(LOGBIO = round(log10(BIOMASS * 10+10))) - - # trim the data....to prepare to find stations only - survdat_stations <- survdat %>% - dplyr::select(CRUISE6, STATION, STRATUM, YEAR) %>% - distinct() - - # make table of strata by year - numtowsstratyr <- table(survdat_stations$STRATUM,survdat_stations$YEAR) - - # find records to keep based on core strata - rectokeep <- stratareas$STRATA %in% strata_used - - # add rec to keep to survdat - stratareas <- cbind(stratareas,rectokeep) - - # delete record form non-core strata - stratareas_usedonly <- stratareas[!stratareas$rectokeep=="FALSE",] - - areapertow=numtowsstratyr - - #compute area covered per tow per strata per year - for(i in 1:50){ - areapertow[,i]=stratareas_usedonly$AREA/numtowsstratyr[,i] - } - - # change inf to NA and round and out in DF - areapertow[][is.infinite(areapertow[])]=NA - areapertow=round(areapertow) - areapertow=data.frame(areapertow) - colnames(areapertow) <- c("STRATUM","YEAR","AREAWT") - areapertow$STRATUM <- as.numeric(as.character(areapertow$STRATUM)) - areapertow$YEAR <- as.numeric(as.character(areapertow$YEAR)) - - survdat <- survdat %>% - inner_join(.,areapertow, by= c("STRATUM","YEAR")) %>% - dplyr::rename(AREAPERTOW = AREAWT) - - # add col to survdat for PLOTWT - survdat$PLOTWT <- NA - survdat$PLOTWT <- ceiling(survdat$AREAPERTOW/1000*survdat$LOGBIO/9) - - if (!is.null(exclude_years)){ - sdat <- survdat %>% filter(!YEAR %in% exclude_years) - } else { - sdat <- survdat - } - - # read species list - sps <- ecodata::species_groupings %>% filter(!is.na(SVSPP)) %>% - dplyr::select(COMNAME, SVSPP) - sps <- sps[!duplicated(sps),] - numsps <- nrow(sps) - - # graph par - par(mar = c(0,0,0,0)) - par(oma = c(0,0,0,0)) - - # index 1:numsps, or by species record number for one species, i.e.25:25 - tspe <- sps %>% filter(COMNAME == species) - - # start map - map("worldHires", xlim=c(-77,-65),ylim=c(35,45), fill=T,border=0,col="gray") - map.axes() - - plot(nesbath,deep=-200, shallow=-200, step=1,add=T,lwd=1,col="gray50",lty=2) - - - # for base period, 1970 to 1979, find call lons for species and by biomass weighting - minyr=1969;maxyr=1980 - clons1 = - sdat$LON[(sdat$YEAR>minyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEAR +# ```{r, echo = T, message= F, warning=F, include=T, eval = T} +library(ecodata) +library(maps) +library(mapdata) +library(ks) +library(marmap) +library(raster) +library(geosphere) +library(dplyr) + +# plot ks plots + +#----STUFF TO SET----------------- +data.dir <- here::here("data") +gis.dir <- here::here("gis") + +# Gives most recent period. 2015-2019 input is 2016-2018 data +rminyr <- 2015 +rmaxyr <- 2019 + +# tlevel is density of color for KD contours areas +tlevel=75 # move later +color_b="blue" +color_r="orange3" +color_r="tomato3" + +# Code to read in strata and compute areas. Or read from cache. + +# readin in strata.shp and compute areas of strata +# TrawlStrata <- raster::shapefile(file.path(gis.dir,"BTS_Strata.shp")) +# AREA <- geosphere::areaPolygon(TrawlStrata, r=6371000)/10^6 + +# for array of strata and area and make into dataframe +# stratareas <- cbind(TrawlStrata@data$STRATA, AREA) +# colnames(stratareas) <- c("STRATA","AREA") +# stratareas <- data.frame(stratareas) +load(file.path(data.dir, "StratAreas.Rdata")) + +# Query bathymetry or load from cache +# getNOAA.bathy(lon1 = -77, lon2 = -65, lat1 = 35, lat2 = 45, +# resolution = 10) -> nesbath +load(file.path(data.dir, "nesbath.Rdata")) + +#Load raw survey data +load(file.path(data.dir, "Survdat.RData")) + +# MUST run addTrans function +# color transparency +addTrans <- function(color,trans){ + # This function adds transparancy to a color. + # Define transparancy with an integer between 0 and 255 + # 0 being fully transparant and 255 being fully visable + # Works with either color and trans a vector of equal length, + # or one of the two of length 1. + + if (length(color)!=length(trans)&!any(c(length(color), + length(trans))==1)) + stop("Vector lengths not correct") + if (length(color)==1 & length(trans)>1) color <- rep(color,length(trans)) + if (length(trans)==1 & length(color)>1) trans <- rep(trans,length(color)) + + num2hex <- function(x) + { + hex <- unlist(strsplit("0123456789ABCDEF",split="")) + return(paste(hex[(x-x%%16)/16+1],hex[x%%16+1],sep="")) + } + rgb <- rbind(col2rgb(color),trans) + res <- paste("#",apply(apply(rgb,2,num2hex),2,paste,collapse=""),sep="") + return(res) +} + +plot_kd <- function(species, season, exclude_years){ + + # stata to use + # offshore strata to use + CoreOffshoreStrata <- c(seq(1010,1300,10),1340, seq(1360,1400,10),seq(1610,1760,10)) + + # inshore strata to use, still sampled by Bigelow + CoreInshore73to12 <- c(3020, 3050, 3080 ,3110 ,3140 ,3170, 3200, 3230, + 3260, 3290, 3320, 3350 ,3380, 3410 ,3440) + # combine + strata_used <- c(CoreOffshoreStrata,CoreInshore73to12) + + survdat <- survdat %>% + dplyr::select(c(CRUISE6,STATION,STRATUM,SVSPP,YEAR, + SEASON,LAT,LON,ABUNDANCE,BIOMASS)) %>% + filter(SEASON == season, + STRATUM %in% strata_used) %>% # delete record form non-core + #strata and get unique records, + # should be one per species + distinct() %>% + # add field with rounded BIOMASS scaler used to adjust distributions + mutate(LOGBIO = round(log10(BIOMASS * 10+10))) + + # trim the data....to prepare to find stations only + survdat_stations <- survdat %>% + dplyr::select(CRUISE6, STATION, STRATUM, YEAR) %>% + distinct() + + # make table of strata by year + numtowsstratyr <- table(survdat_stations$STRATUM,survdat_stations$YEAR) + + # find records to keep based on core strata + rectokeep <- stratareas$STRATA %in% strata_used + + # add rec to keep to survdat + stratareas <- cbind(stratareas,rectokeep) + + # delete record form non-core strata + stratareas_usedonly <- stratareas[!stratareas$rectokeep=="FALSE",] + + areapertow=numtowsstratyr + + #compute area covered per tow per strata per year + for(i in 1:50){ + areapertow[,i]=stratareas_usedonly$AREA/numtowsstratyr[,i] + } + + # change inf to NA and round and out in DF + areapertow[][is.infinite(areapertow[])]=NA + areapertow=round(areapertow) + areapertow=data.frame(areapertow) + colnames(areapertow) <- c("STRATUM","YEAR","AREAWT") + areapertow$STRATUM <- as.numeric(as.character(areapertow$STRATUM)) + areapertow$YEAR <- as.numeric(as.character(areapertow$YEAR)) + + survdat <- survdat %>% + inner_join(.,areapertow, by= c("STRATUM","YEAR")) %>% + dplyr::rename(AREAPERTOW = AREAWT) + + # add col to survdat for PLOTWT + survdat$PLOTWT <- NA + survdat$PLOTWT <- ceiling(survdat$AREAPERTOW/1000*survdat$LOGBIO/9) + + if (!is.null(exclude_years)){ + sdat <- survdat %>% filter(!YEAR %in% exclude_years) + } else { + sdat <- survdat + } + + # read species list + sps <- ecodata::species_groupings %>% filter(!is.na(SVSPP)) %>% + dplyr::select(COMNAME, SVSPP) + sps <- sps[!duplicated(sps),] + numsps <- nrow(sps) + + # graph par + par(mar = c(0,0,0,0)) + par(oma = c(0,0,0,0)) + + # index 1:numsps, or by species record number for one species, i.e.25:25 + tspe <- sps %>% filter(COMNAME == species) + + # start map + map("worldHires", xlim=c(-77,-65),ylim=c(35,45), fill=T,border=0,col="gray") + map.axes() + + plot(nesbath,deep=-200, shallow=-200, step=1,add=T,lwd=1,col="gray50",lty=2) + + + # for base period, 1970 to 1979, find call lons for species and by biomass weighting + minyr=1969;maxyr=1980 + clons1 = + sdat$LON[(sdat$YEAR>minyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEARrminyr & sdat$YEAR