Watershed Model Maps

[jdkleiner]

Needed a way to visually inspect hsp2 model runs. This is a quick way to generate a map of model run vs model run to visualize which model segments have been run successfully for a given model_version and runid.

• Goals:
  • To enable similar function to the old standard mapserver version:
    • Example: click here
    • Page with use info: https://github.com/HARPgroup/vahydro/issues/90
  • Identify improvements needed:
    • TBD
  • Where does this cross over with the last WSP mapping round (and current VWP facility): Automate Model Maps hydro-tools#350
  • Docs on using R version to add the ability to manipulate post generation and would be nice
• Script: https://github.com/HARPgroup/vahydro/blob/master/R/modeling/hsp2_validation/Model_Mapper.R
  • Multiple modes of use:
    • Command line with arguments:
      • Use: Rscript.exe Model_Mapper.R basin_prefix model_a runid_a model_b runid_b image_file_basename
      • Ex: Rscript.exe Model_Mapper.R JU "vahydro-1.0" "runid_13" "cbp-6.1" "subsheds" jtest
    • Command Line interactive:
      • Rscript Model_Mapper.R, and just answer the prompts

Example:

• Comparing a new hsp2 run vs. runid_11 is a good way to quickly see which hsp2 segments ran successfully (runid_11 is a good baseline because all segs were run for runid_11 in development of the 2020 State Plan)

source("https://raw.githubusercontent.com/HARPgroup/hydro-tools/master/GIS_functions/model_geoprocessor.R")
export_path <- 'C:/Users/nrf46657/Desktop/GitHub/vahydro/R/modeling/hsp2_validation/'

# specify plot title, rseg where clause, and scenario model_version, runid pairs
plotname <- "potomac"
segswhere <- "hydrocode LIKE '%wshed_P%'"
scenario_a <- c("vahydro-1.0","runid_11")
scenario_b <- c("cbp-6.0","hsp2_2022")

#---------------------------------------------------------------------------
# process geo data
polygons_a <- model_geoprocessor(scenario_a,segswhere)
polygons_b <- model_geoprocessor(scenario_b,segswhere)

# generate & save plot figure
filename <- paste0("Model_Map_",scenario_a[1],scenario_a[2],"_",scenario_b[1],scenario_b[2],"_",plotname,".png")
png(file=paste(export_path,filename,sep=""), width=1500, height=1500)
plot(polygons_a, axes = 1, main=paste0(plotname,"\n",segswhere), cex.main=2, cex.axis=2.5, border = "black")
plot(polygons_b, col = "blue", add = T)
legend("topleft", legend=c(paste0(scenario_a[1],", ",scenario_a[2], " (",length(polygons_a),")"), paste0(scenario_b[1],", ",scenario_b[2], " (",length(polygons_b),")")), fill = c("white","blue"), cex=3.5)
dev.off()

Result: