Habitat Analysis

#library(rremat)
library(ggplot2)
library(dplyr)
library(RColorBrewer)
library(knitr)
library(scales)
library(tidyr)
opts_chunk$set(fig.width = 12, fig.height = 8, echo = FALSE, error = FALSE)

Habitat definitions

Water volume in the estuary is distributed such that the upper layers contain the most habitat. Both standard and logarithmic color scales are shown.

data(volumes)
cellvol = prod(attr(volumes, "resolution")[c("xres", "yres", "zres")])
volumes["volume"] = cellvol*volumes$count
for(h in c("littoral", "limnetic", "epibenthic", "sublimnetic", "profundal"))
  volumes[paste0("volume.", h)] = cellvol*volumes[[paste0("count.", h)]]

ggplot(volumes[volumes$count > 0,], aes(x = dist, y = elev, fill = volume)) + 
  geom_raster() + scale_fill_gradientn(colours = brewer.pal(n=6,name = "Set1")) +
  scale_y_continuous("elevation above NAVD29 (m)", labels = comma) +
  xlab("Distance from river mouth (m)")
#ggsave(paste0("C:/Users/Michael/Desktop/john_plots/", "volumedistribution.png"),
  width = 10, height = 8)
ggplot(volumes[volumes$count > 0,], aes(x = dist, y = elev, 
  fill = log10(volume))) + geom_raster() + 
  scale_fill_gradientn(colours = brewer.pal(n=6,name = "Set1")) 

Habitat quality for juvenile salmonids is defined according to the specifications outlined by the WQ Parameter and Habitat Productivity Technical Group (10/6/2015):

Temperature:

| optimal growth | positive growth (suitable) | no or negative growth (stressful) | unsuitable | |:---------------:|:---------------------------------:|:---------------------------------:|:--------------:| | $14-18^\circ$ C | $< 14^\circ$ C or $18-21^\circ$ C | $21--25^\circ$ C | $> 25^\circ$ C |

Dissolved Oxygen:

| function without impairment | suitable | limited by lack of oxygen | unsuitable | |:---------------------------:|:----------:|:-------------------------:|:----------:| | $> 6$ mg/L | $4-6$ mg/L | $3-4$ mg/L | $< 3$ mg/L |

Salinity:

| ion importing | ion neutral | induced ion exporting | ion exporting | |:-------------:|:-----------:|:---------------------:|:-------------:| | $< 10$ ppt | $10-15$ ppt | $15-28$ ppt | $> 28$ ppt |

Overall habitat is assessed to fall into one of 13 categories:

| category | Description | |:---------------------------------|:----------------------------------------------------------------------------------------| | Optimal | all of: optimal growth (T), ion importing (S), function without impairment (DO) | | Unsuitable | any one of: unsuitable (T), ion-exporting (S), unsuitable (DO) | | Sub-optimal (T) | positive growth (T), (S) and (DO) as in "optimal" category | | Stressful (T) | no or negative growth (T), (S) and (DO) as in "optimal" category | | Sub-optimal (DO) | either suitable or limited by lack of oxygen (DO), (T) and (S) as in "optimal" category | | Sub-optimal (S) | either ion neutral or induced ion exporting (S), (T) and (DO) as in "optimal" category | | Sub-optimal (T, DO) | (S) as in "optimal" category | | Sub-optimal (T, S) | (DO) as in "optimal" category | | Sub-optimal (DO, S) | (T) as in "optimal" category | | Stressful sub-optimal (T, DO) | (S) as in "optimal" category | | Stressful sub-optimal (T, S) | (DO) as in "optimal" category | | Sub-optimal (T, S, DO) | Suitable (T), (S) and (DO) as in "sub-optimal" categories | | Stressful sub-optimal (T, S, DO) | Stressful (T), (S) and (DO) as in "sub-optimal" categories |

Habitat quality can additionally be assessed based on water depth:

| Category | Description | Criteria | |:--------------------------|:--------------------------------------------------------------------------------------------------------------|:--------------------------------------:| | Shallow/shoal (Littoral) | high prey availability; high avian predation risk | Depth < 2m and \max{depth} < 2m | | Surface (Limnetic) | lower prey availability; high avian predation risk | $depth < 2m$ and $\max{depth} \geq 2m$ | | Subsurface (Limnetic) | high prey availability where sediment <5m; low prey availability where depth >5m); lower avian predation risk | $2m \leq depth < 5m$ | | Stagnant (Profundal) zone | unsuitable due to hypoxia/anoxia; low avian predation risk | $depth \geq 5m$ |

ta.windows = list(
  "optimal" = expression(x >= 14 & x <= 18),
  "suitable" = expression(x < 14 | (x >= 18 & x <= 21)),
  "stressful" = expression(x > 21 & x <= 25),
  "unsuitable" = expression(x > 25)
)
ta_hsd = make_hsd(ta.windows, names(ta.windows))

oa.windows = list(
  "no impairment" = expression(x >= 6),
  "suitable" = expression(x >= 4 & x < 6),
  "limited" = expression(x >= 3 & x < 4),
  "unsuitable" = expression(x < 3)
)
oa_hsd = make_hsd(oa.windows, names(oa.windows))

sa.windows = list(
  "ion importing" = expression(x < 10),
  "ion neutral" = expression(x >= 10 & x <= 15),
  "induced ion exporting" = expression(x > 15 & x <= 28),
  "ion exporting" = expression(x > 28)
)
sa_hsd = make_hsd(sa.windows, names(sa.windows))


hab_cat = function(xs){ # xs = c(ta, oa, sa)
  if(any(xs %in% c("unsuitable", "ion exporting")))
    "unsuitable"
  else if(all(xs %in% c("optimal", "no impairment", "ion importing")))
    "optimal"  
  else if(xs[[1]] == "suitable" && all(xs[2:3] %in% c("no impairment", "ion importing")))
    "sub-optimal (T)"
  else if(xs[[1]] == "stressful" && all(xs[2:3] %in% c("no impairment", "ion importing")))
    "stressful (T)"
  else if(xs[[2]] %in% c("suitable", "limited") && all(xs[c(1,3)] %in% c("optimal", "ion importing")))
    "sub-optimal (DO)"
  else if(all(xs[1:2] %in% c("optimal", "no impairment")) && xs[[3]] %in% c("ion neutral", "induced ion exporting"))
    "sub-optimal (S)"
  else if(all(xs[1:2] %in% c("suitable", "limited")) && xs[[3]] %in% c("ion importing"))
    "sub-optimal (T, DO)"
  else if(all(xs[c(1,3)] %in% c("suitable", "ion neutral", "induced ion exporting")) && xs[[2]] %in% c("no impairment"))
    "sub-optimal (T, S)"
  else if(all(xs[2:3] %in% c("suitable", "limited", "ion neutral", "induced ion exporting")) && xs[[1]] %in% c("optimal"))  
    "sub-optimal (DO, S)" 
  else if(xs[[1]] %in% c("stressful") && xs[[2]] %in% c("suitable", "limited") && xs[[3]] %in% c("ion importing"))
    "stressful sub-optimal (T, DO)"
  else if(all(xs[c(1, 3)] %in% c("stressful", "ion neutral", "induced ion exporting")) && xs[[2]] %in% c("no impairment"))
    "stressful sub-optimal (T, S)"
  else if(all(xs[2:3] %in% c("suitable", "limited", "ion neutral", "induced ion exporting")) && xs[[1]] %in% c("suitable"))
    "sub-optimal (T, DO, S)"
  else if(all(xs[2:3] %in% c("suitable", "limited", "ion neutral", "induced ion exporting")) && xs[[1]] %in% c("stressful"))
    "stressful sub-optimal (T, DO, S)"
  else
    stop("no category for ", paste(c("ta", "oa", "sa"), xs, sep = " = ", collapse = ", "))
}

tot_hab = make_hsi(list(ta_hsd, oa_hsd, sa_hsd), hab_cat)

Assessing habitat quality

data(grids)
habgrids = grids
habgrids["ta.qual"] = sapply(habgrids$ta, ta_hsd)
habgrids["sa.qual"] = sapply(habgrids$sa, sa_hsd)
habgrids["oa.qual"] = sapply(habgrids$oa, oa_hsd)
habgrids["habitat"] = factor(apply(grids[c("ta", "oa", "sa")], 1, tot_hab),
  levels = c("optimal", "sub-optimal (T)", "sub-optimal (DO)", 
    "sub-optimal (S)", "sub-optimal (T, DO)", "sub-optimal (T, S)",
    "sub-optimal (DO, S)", "sub-optimal (T, DO, S)", "stressful (T)", 
    "stressful sub-optimal (T, DO)", "stressful sub-optimal (T, S)", 
    "stressful sub-optimal (T, DO, S)", "unsuitable"))
habgrids["ta.qual"] = factor(habgrids$ta.qual, levels = names(ta.windows))
habgrids["sa.qual"] = factor(habgrids$sa.qual, levels = names(sa.windows))
habgrids["oa.qual"] = factor(habgrids$oa.qual, levels = names(oa.windows))
# add closure meta data
data(closures)
habgrids["code"] = "O"
habgrids["days.since.closure"] = NA
for(d in unique(habgrids$date)){
  ind = which(d >= closures$start & d <= closures$end)
  if(length(ind) > 0){
    habgrids[habgrids$date == d, "code"] = as.character(closures[ind, "code"])
    if(as.character(closures[ind, "code"]) == "C")
      habgrids[habgrids$date == d, "days.since.closure"] = 
        habgrids[habgrids$date == d, "date"] - 
        closures$start[max(which(d >= closures$start))]
  }
}
habgrids["code"] = factor(habgrids$code)
# add ctd meta data
data(ctdmeta)
gridid = as.character(interaction(habgrids$date, habgrids$id))
metaid = as.character(interaction(strftime(ctdmeta$start, format = "%Y-%m-%d", 
  tz = "US/Pacific"), ctdmeta$id))
for(i in seq(nrow(ctdmeta))){
  habgrids[gridid == metaid[i], "numcasts"] = ctdmeta[i, "numcasts"]  
}

Habitat can be summarized in a variety of different ways: we can consider the estuary as whole, or divide it into segments based on distance from the river mouth (lower reach, middle reach, upper reach) or vertical strata (e.g. by depth categories). We can also group transects by estuary mouth condition (i.e. open, closed, perched), by season, or by inflow (e.g. high or low inflow).

# strip out incomplete grids
habgrids = filter(habgrids, numcasts == 12)
# total habitat
overall = summarize_by_strata(habgrids, 
  stratcols = c("date", "id", "habitat", "code", "days.since.closure"), 
  summexpr = c(volume = ~sum(volume)))

# stratify by depth (above/below 2m) and by depth qual
bydepthqual = summarize_by_strata(habgrids, stratcols = c("date", "id", "code", 
  "depth.qual", "habitat"), summexpr = c(volume = ~sum(volume)))
habgrids["depth.zone"] = ifelse(habgrids$depth < 2, "surface", "subsurface")
bydepth = summarize_by_strata(habgrids, stratcols = c("date", "id", "code", 
  "depth.zone", "habitat"), summexpr = c(volume = ~sum(volume)))

# stratify habitat by inner vs. out estuary
#lz = c(1700, 5800)
#habgrids["distzone"] = stratify(habgrids$dist, lz)
#bydist = summarize_by_strata(habgrids, stratcols = c("date", "id", "code",
#  "distzone", "habitat"), summexpr = c(volume = ~sum(volume)))

Visualizing habitat

# total habitat
ggplot(overall, aes(x = date, y = volume)) + 
  geom_bar(aes(fill = habitat), stat="identity", position = "stack") +  
  stat_summary(fun.y = sum, geom = "point") + facet_wrap(~ code) + coord_flip()

ggplot(overall, aes(x = code, y = volume, fill = habitat)) + geom_boxplot() + 
  facet_wrap(~habitat, scales = "free", ncol=1)

ggplot(overall, aes(x = volume, fill = habitat)) + facet_grid(code ~ habitat) +
  geom_histogram(binwidth = 1e5, color = "black")

ggplot(overall, aes(x = wse, y = volume, color = habitat)) + 
  geom_smooth(method = "lm") + geom_point(aes(shape = id)) + 
  facet_wrap(~ code, ncol = 1)

ggplot(overall, aes(x = days.since.closure, y = volume, color = habitat)) + 
  geom_smooth() + geom_point() + facet_wrap(~ code, ncol = 1) + 
  facet_wrap(~ habitat, ncol = 1, scales = "free_y")

#ggplot(overall, aes(x = factor(days.since.closure), y = volume, fill = habitat)) + 
#  geom_boxplot() + facet_wrap(~ code, ncol = 1) + facet_wrap(~ habitat, 
#  ncol = 1, scales = "free_y")


ggplot(byelev, aes(x = code, y = volume, fill = habitat)) + geom_boxplot() + 
  facet_grid(elevzone ~ habitat, scales = "free")

ggplot(bypycnocline, aes(x = code, y = volume, fill = habitat)) + geom_boxplot() + 
  facet_grid(pycnozone ~ habitat, scales = "free")




# stratify summary by depth zones
#ggplot(byelev, aes(x = factor(date):factor(id), y = volume, fill = habitat)) + 
#  geom_bar(stat = "identity") + facet_wrap(~ elevzone)

#ggplot(bypycnocline, aes(x = factor(date):factor(id), y = volume, fill = habitat)) + 
#  geom_bar(stat = "identity") + facet_wrap(~ elevzone)


# stratify habitat by inner vs. out estuary
#ggplot(bydist, aes(x = factor(date):factor(id), y = volume, fill = habitat)) + 
#  geom_bar(stat = "identity") + facet_wrap(~ distzone)


mkoohafkan/rremat documentation built on Aug. 6, 2019, 8:59 a.m.