In mkoohafkan/rremat: Russian River Estuary data collected by the Largier Lab

Analyze event

# plot setup
baseplot = function(p) 
  ggplot(p, aes(x = dist, y = elev)) + 
  xlim(min(thesegrids$dist), max(thesegrids$dist)) +
  ylim(min(thesegrids$elev), max(thesegrids$elev)) +
  ylab("elevation above NAVD29 (m)") +
  xlab("Distance from river mouth (m)") 

lyrs = list(
  ta.qual = geom_raster(aes(fill = ta.qual)),
  sa.qual = geom_raster(aes(fill = sa.qual)),
  oa.qual = geom_raster(aes(fill = oa.qual)),
  habitat = geom_raster(aes(fill = habitat))
)

legend.settings = theme(legend.position = "right", 
  legend.text = element_text(size = 6), 
  legend.title = element_text(size = 6))

habfill = list(
  habitat = scale_fill_manual("habitat\nquality\n(overall)", 
    values = c(
      "unsuitable" = "#000000",
      "optimal" = "#1f78b4",
      "Sub-optimal (T)" = "#a6cee3",
      "stressful (T)" = "#ffff99",
      "sub-optimal (DO)" = "#33a02c",
      "sub-optimal (S)" = "#6a3d9a",
      "sub-optimal (T, DO)" = "#fb9a99",
      "sub-optimal (T, S)" = "#fdbf6f",
      "sub-optimal (DO, S)" = "#b2df8a",
      "stressful sub-optimal (T, DO)" = "#e31a1c",
      "stressful sub-optimal (T, S)" = "#ff7f00",
      "sub-optimal (T, DO, S)" = "#cab2d6",
      "stressful sub-optimal (T, DO, S)" = "#b15928"    
    )),
  ta.qual = scale_fill_manual("habitat\nquality\n(temperature)", 
    values = setNames(brewer.pal(4, "RdYlGn"), c("unsuitable", 
    "stressful", "suitable", "optimal"))),
  sa.qual = scale_fill_manual("habitat\nquality\n(salinity)", 
    values = setNames(brewer.pal(4, "RdYlGn"), c("ion.exporting", 
    "induced.ion.exporting", "ion.neutral", "ion.importing"))),
  oa.qual = scale_fill_manual("habitat\nquality\n(dissolved oxygen)", 
    values = setNames(brewer.pal(4, "RdYlGn"), c("unsuitable", "limited", 
    "suitable", "no.impairment")))
)

habcolor = list(
  habitat = scale_color_manual("habitat\nquality\n(overall)", 
        values = c(
      "unsuitable" = "#000000",
      "optimal" = "#1f78b4",
      "Sub-optimal (T)" = "#a6cee3",
      "stressful (T)" = "#ffff99",
      "sub-optimal (DO)" = "#33a02c",
      "sub-optimal (S)" = "#6a3d9a",
      "sub-optimal (T, DO)" = "#fb9a99",
      "sub-optimal (T, S)" = "#fdbf6f",
      "sub-optimal (DO, S)" = "#b2df8a",
      "stressful sub-optimal (T, DO)" = "#e31a1c",
      "stressful sub-optimal (T, S)" = "#ff7f00",
      "sub-optimal (T, DO, S)" = "#cab2d6",
      "stressful sub-optimal (T, DO, S)" = "#b15928"    
    )),
  ta.qual = scale_color_manual("habitat\nquality\n(temperature)", 
    values = setNames(brewer.pal(4, "RdYlGn"), c("unsuitable", 
    "stressful", "suitable", "optimal"))),
  sa.qual = scale_color_manual("habitat\nquality\n(salinity)", 
    values = setNames(brewer.pal(4, "RdYlGn"), c("ion.exporting", 
    "induced.ion.exporting", "ion.neutral", "ion.importing"))),
  oa.qual = scale_color_manual("habitat\nquality\n(dissolved oxygen)", 
    values = setNames(brewer.pal(4, "RdYlGn"), c("unsuitable", "limited", 
    "suitable", "no.impairment")))
)

plot_title = function(plotvar){
  if(plotvar == "habitat"){
    pt = "overall habitat quality"
  } else if(plotvar == "ta.qual"){
    pt = "temperature habitat quality"
  } else if(plotvar == "sa.qual"){
    pt = "salinity habitat quality"
  } else if(plotvar == "oa.qual"){
    pt = "dissolved oxygen habitat quality"
  } else {
    stop("value of 'plotvar' not recognized")
  }
  pt
}


################
# grids
############# 

plotvars = c("habitat", "ta.qual", "sa.qual", "oa.qual")
for(plotvar in plotvars){
  pt = plot_title(plotvar)
  for(d in unique(paste(thesegrids$date))){
    p = baseplot(filter(thesegrids, date == d))  + 
      lyrs[[plotvar]] + habfill[[plotvar]] + legend.settings +
      ggtitle(paste(pt, d)) 
    ggsave(paste0("C:/Users/Michael/Desktop/john_plots/", plotvar, "_", d, ".png"),
      p, width = 11, height = 8.5)
  }
}

##########
# overall
##########

plotvars = c("habitat", "ta.qual", "sa.qual", "oa.qual")
for(plotvar in plotvars){
  pt = plot_title(plotvar)
  hablevels = unique(thesegrids[[plotvar]])
  overall = summarize_by_strata(thesegrids, stratcols = c("date", "id",
    plotvar, "code", "days.since.closure"), summexpr = c(volume = ~sum(volume)))
  overall.spread = overall %>% spread_(plotvar, "volume", fill = 0) 
  overall.spread["total.volume"] = rowSums(overall.spread[hablevels])
  overall.gathered = overall.spread %>% gather_(plotvar, "volume", 
    gather_cols = hablevels)
  overall.gathered["volume.frac"] = overall.gathered$volume/
                                      overall.gathered$total.volume
  overall.gathered = left_join(overall.gathered, 
    unique(thesegrids[c("date", "id", "wse")]), by = c("date", "id"))

  if(closureplot){ # closed
    p = ggplot(overall.gathered, aes_string(x = "date", y = "volume", 
      fill = plotvar)) + 
      geom_area(position = "stack") +
      geom_vline(aes(xintercept = as.numeric(date)), size = 1, linetype = "dashed") + 
      habfill[[plotvar]] + legend.settings + xlab("") +
      scale_y_continuous(name="Volume (m3)", labels = comma)
  } 
  else { # open
    p = ggplot(overall.gathered, aes_string(x = "factor(date)", y = "volume", 
      fill = plotvar)) + geom_bar(stat = "identity", position = "stack", color = "black") + 
      habfill[[plotvar]] + legend.settings + xlab("") +
      scale_y_continuous(name="Volume (m3)", labels = comma)  
  }
  ggsave(paste0("C:/Users/Michael/Desktop/john_plots/", "overall_", plotvar, 
    ".png"), p, width = 10, height = 8)
}

#####################################
# stratify by depth
#####################################

thesegrids["depthzone"] = factor(ifelse(thesegrids$depth < 2,
  "surface (< 2m below surface)", "subsurface (> 2m below surface)"))

plotvars = c("habitat", "ta.qual", "sa.qual", "oa.qual")
for(plotvar in plotvars){
  pt = plot_title(plotvar)
  hablevels = unique(thesegrids[[plotvar]])
  bydepth = summarize_by_strata(thesegrids, 
    stratcols = c("date", "id", "depthzone", plotvar), 
    summexpr = c(volume = ~sum(volume)))
  bydepth.spread = bydepth %>% spread_(plotvar, "volume", fill = 0) 
  bydepth.spread["total.volume"] = rowSums(bydepth.spread[hablevels])
  for(h in hablevels)
    bydepth.spread[paste0("frac.", h)] = bydepth.spread[[h]]/
                                           bydepth.spread$total.volume 
  bydepth.spread = left_join(bydepth.spread, 
    unique(thesegrids[c("date", "wse")]), "date")
  bydepth.gathered = bydepth.spread %>% gather_(plotvar, "volume", 
    gather_cols = hablevels)
  bydepth.gathered["volume.frac"] = bydepth.gathered$volume/
                                      bydepth.gathered$total.volume

  if(closureplot){ # closed
    p = ggplot(bydepth.gathered, aes_string(x = "date", y = "volume", fill = plotvar)) +
      geom_area(position = "stack") +
      geom_vline(aes(xintercept = as.numeric(date)), size = 1, linetype = "dashed") + 
      habfill[[plotvar]] + legend.settings + xlab("") +
      scale_y_continuous(name="Volume (m3)", labels = comma) +
      facet_wrap(~depthzone) 
  } 
  else { # open
    p = ggplot(bydepth.gathered, aes_string(x = "factor(date)", y = "volume", fill = plotvar)) + 
      geom_bar(stat = "identity", position = "stack", color = "black") + 
      habfill[[plotvar]] + legend.settings + xlab("") +
      scale_y_continuous(name="Volume (m3)", labels = comma) +
      facet_wrap(~depthzone)
  }
  ggsave(paste0("C:/Users/Michael/Desktop/john_plots/", "bydepth_", plotvar, 
    ".png"), p, width = 10, height = 8)
}



########################
#stratify by transect locations
##########################

hablevels = unique(thesegrids$habitat)
a = unique(ctd$dist)
acuts = a[-length(a)] + diff(a)/2
thesegrids["transectid"] = stratify(thesegrids$dist, acuts) 
levels(thesegrids$transectid) <- a
thesegrids$transectid = as.numeric(as.character(thesegrids$transectid))

plotvars = c("habitat", "ta.qual", "sa.qual", "oa.qual")
for(plotvar in plotvars){
  pt = plot_title(plotvar)
  hablevels = unique(thesegrids[[plotvar]])
  bytid = summarize_by_strata(thesegrids, 
    stratcols = c("date", "id", "transectid", plotvar), 
    summexpr = c(volume = ~sum(volume)))
  names(bytid)[names(bytid) == "transectid"] = "dist"
  bytid.spread = bytid %>% spread_(plotvar, "volume", fill = 0) 
  bytid.spread["total.volume"] = rowSums(bytid.spread[hablevels])
  for(h in hablevels)
    bytid.spread[paste0("frac.", h)] = bytid.spread[[h]]/bytid.spread$total.volume 
  bytid.spread = left_join(bytid.spread, unique(thesegrids[c("date", "wse")]), 
    "date")
  bytid.gathered = bytid.spread %>% gather_(plotvar, "volume", 
    gather_cols = hablevels)
  bytid.gathered["volume.frac"] = bytid.gathered$volume/bytid.gathered$total.volume

  p = ggplot(bytid.gathered, aes_string(x = "date", y = "volume", 
    color = plotvar)) + geom_line(size = 1) + geom_point(size = 3) + 
    geom_line(aes(y = total.volume), color = "black", linetype = "dashed", size = 1) + 
    habcolor[[plotvar]] + legend.settings + xlab("") +
    scale_y_continuous(name="Volume (m3)", labels = comma) +
    facet_wrap(~dist) + theme(legend.position = "bottom")
  ggsave(paste0("C:/Users/Michael/Desktop/john_plots/", "bytransect_", plotvar, 
    ".png"), p, width = 10, height = 8)
}

##########################
# mapping by transect id
#########################

library(arcpyr)
arcpy.initialize()
sa.initialize()
for(thisdate in unique(paste(bytid.spread$date))){
  # write grid for this date
  thisgrid = filter(bytid.spread, date == thisdate)
  thisgrid$dist = paste((thisgrid$dist))
  thisid = as.character(unique(thisgrid$id))
  suffix = paste0("_", gsub("-", "", thisdate), thisid, "_tid")
  thiswse = unique(thisgrid$wse)
  write.csv(thisgrid, file = "C:/GIS workspace/RRE/rrzone_bytid.txt", row.names=FALSE)  
  # copy grid to gdb
  TableToTable_conversion("C:/GIS workspace/RRE/rrzone_bytid.txt", 
    "C:/GIS workspace/RRE/habitat.gdb", "rrzone_hab")
  # define wse extent
  RasterCalculator(paste0('rastermask = Con(bathy <= ', thiswse, ', 1)'), 
    list(bathy = "C:/GIS workspace/RRE/habitat.gdb/bathymetry_NGVD_meters"),
    list(rastermask = "C:/GIS workspace/RRE/habitat.gdb/rastermask"))
  RasterToPolygon_conversion("C:/GIS workspace/RRE/habitat.gdb/rastermask", 
    "C:/GIS workspace/RRE/habitat.gdb/polymask", "SIMPLIFY", "VALUE")
  # clip zones to wse extent
  Clip_analysis("C:/GIS workspace/RRE/habitat.gdb/markerzones_tid_dissolve", 
    "C:/GIS workspace/RRE/habitat.gdb/polymask", 
    paste0("C:/GIS workspace/RRE/habitat.gdb/rrzone", suffix))
  # join habitat data to clipped zones
  JoinField_management(paste0("C:/GIS workspace/RRE/habitat.gdb/rrzone", suffix), 
    "transectid", "C:/GIS workspace/RRE/habitat.gdb/rrzone_hab", "dist")
  # clean up
  Delete_management("C:/GIS workspace/RRE/habitat.gdb/rastermask")
  Delete_management("C:/GIS workspace/RRE/habitat.gdb/polymask")
  Delete_management("C:/GIS workspace/RRE/habitat.gdb/rrzone_hab")
  file.remove("C:/GIS workspace/RRE/rrzone_bytid.txt")
}

#######################
# mapping grids
######################3

hablevels = unique(thesegrids$habitat)
bydist = summarize_by_strata(thesegrids, 
  stratcols = c("date", "id", "dist", "habitat", 
  "days.since.closure"), summexpr = c(volume = ~sum(volume)))
bydist.spread = bydist %>% spread(habitat, volume, fill = 0) 
bydist.spread["total.volume"] = rowSums(bydist.spread[hablevels])
for(h in hablevels)
  bydist.spread[paste0("frac.", h)] = bydist.spread[[h]]/bydist.spread$total.volume 
bydist.spread = left_join(bydist.spread, unique(thesegrids[c("date", "wse")]), 
  "date")

#bydist.gathered = bydist.spread %>% gather_("habitat", "volume", hablevels)
#bydist.gathered["volume.frac"] = bydist.gathered$volume/bydist.gathered$total.volume
#bydist.gathered = left_join(bydist.gathered, 
#  unique(thesegrids[c("date", "dist", "wse")]), c("dist","date"))

library(arcpyr)
for(thisdate in unique(paste(bydist$date))){
  # write grid for this date
  thisgrid = filter(bydist.spread, date == thisdate)
  thisgrid$dist = factor(thisgrid$dist)
  thisid = as.character(unique(thisgrid$id))
  suffix = paste0("_", gsub("-", "", thisdate), thisid)
  thiswse = unique(thisgrid$wse)
  write.csv(thisgrid, file = "C:/GIS workspace/RRE/rrzone_hab.txt")  
  # copy grid to gdb
  TableToTable_conversion("C:/GIS workspace/RRE/rrzone_hab.txt", 
    "C:/GIS workspace/RRE/habitat.gdb", "rrzone_hab")
  # define wse extent
  RasterCalculator(paste0('rastermask = Con(bathy <= ', thiswse, ', 1)'), 
    list(bathy = "C:/GIS workspace/RRE/habitat.gdb/bathymetry_NGVD_meters"),
    list(rastermask = "C:/GIS workspace/RRE/habitat.gdb/rastermask"))
  RasterToPolygon_conversion("C:/GIS workspace/RRE/habitat.gdb/rastermask", 
    "C:/GIS workspace/RRE/habitat.gdb/polymask", "SIMPLIFY", "VALUE")
  # clip zones to wse extent
  Clip_analysis("C:/GIS workspace/RRE/habitat.gdb/markerzones_poly_dissolve", 
    "C:/GIS workspace/RRE/habitat.gdb/polymask", 
    paste0("C:/GIS workspace/RRE/habitat.gdb/rrzone", suffix))
  # join habitat data to clipped zones
  JoinField_management(paste0("C:/GIS workspace/RRE/habitat.gdb/rrzone", suffix), 
    "gridcode", "C:/GIS workspace/RRE/habitat.gdb/rrzone_hab", "dist")
  # clean up
  Delete_management("C:/GIS workspace/RRE/habitat.gdb/rastermask")
  Delete_management("C:/GIS workspace/RRE/habitat.gdb/polymask")
  Delete_management("C:/GIS workspace/RRE/habitat.gdb/rrzone_hab")
  file.remove("C:/GIS workspace/RRE/rrzone_hab.txt")
}