ras2fim debugging

Step 1) Verify the model

One of the first things we should do is verify that this model doesn’t exhibit any of the hallmarks of models which RAS2FIM was not developed for: A RAS2FIM compliant HEC-RAS model. One of the most direct ways is to open it up using the HEC-RAS GUI, but if it can run through RRASSLER that’s a pretty good indicator that it’s compliant. Let’s test this process with a model!

library(data.table)
library(RRASSLER)

xs_to_test = 3

# Path to raw model
source_model_path <- file.path(fs::path_package("", package = "RRASSLER"),"extdata","sample_ras","FEMA-R6-BLE-sample-dataset","12090301","12090301_models","Model","Alum Creek-Colorado River","ALUM 114","ALUM 114.prj",fsep = .Platform$file.sep)
# source_model_path <- file.path("~/Dropbox/root/projects/RRASSLER/RRASSLER/inst/extdata/sample_ras/FEMA-R6-BLE-sample-dataset/12090301/12090301_models/Model/Alum Creek-Colorado River/ALUM 114/ALUM 114.prj")

current_model_projection = "EPSG:2277"
current_model_units = "English Units"

# Path to ras catalog
# ras_dbase <- file.path(fs::path_package("", package = "RRASSLER"),"extdata","sample_output","ras_catalog",fsep = .Platform$file.sep)
# ras_dbase <- file.path("~/Dropbox/root/projects/RRASSLER/RRASSLER/inst/extdata/sample_output/ras_cat")
# ras_dbase <- file.path("C:/Users/stick/root/Dropbox/root/projects/RRASSLER/RRASSLER/inst/extdata/sample_output/ras_catalog")

## Pick a dem
# dap_catalog <- jsonlite::read_json('https://raw.githubusercontent.com/mikejohnson51/climateR-catalogs/refs/heads/main/docs/catalog.json', simplifyVector = TRUE)
# dem <- climateR::dap(URL = dem_catalog$URL, AOI = AOI::aoi_get(AOI_lines %>% sf::st_buffer(50)))
# dap_catalog <- jsonlite::stream_in(url('https://raw.githubusercontent.com/mikejohnson51/climateR-catalogs/refs/heads/main/docs/catalog.json'))
# dem_catalog = dplyr::filter(dap_catalog,id == "USGS 3DEP",asset=="10m CONUS DEM")

# Precompute and clean
## Guess at values
prj_files <- list.files(dirname(source_model_path),pattern = utils::glob2rx(glue::glue("{tools::file_path_sans_ext(basename(source_model_path))}.prj$")),full.names = TRUE,ignore.case = TRUE,recursive = TRUE)
for (potential_file in prj_files) {
  # potential_file <- prj_files[1]
  file_text <- read.delim(potential_file, header = FALSE)
  
  if (any(c('PROJCS', 'GEOGCS', 'DATUM', 'PROJECTION') == file_text)) {
    current_model_projection = potential_file
  } else if (grepl("SI Units", file_text, fixed = TRUE)) {
    current_last_modified = as.integer(as.POSIXct(file.info(potential_file)$mtime))
    current_model_units = "SI Units"
  } else if (grepl("English Units", file_text, fixed = TRUE)) {
    current_last_modified = as.integer(as.POSIXct(file.info(potential_file)$mtime))
    current_model_units = "English Units"
  }
}

root_path <- file.path(dirname(source_model_path),glue::glue("{tools::file_path_sans_ext(basename(source_model_path))}.g01"),fsep = .Platform$file.sep)
database_approach <- RRASSLER::parse_model_to_xyz(geom_path = root_path,units = current_model_units,proj_string = current_model_projection,quiet = FALSE,is_verbose = TRUE)
g_approach <- RRASSLER::process_ras_g_to_xyz(geom_path = root_path,units = current_model_units,proj_string = current_model_projection,vdat = FALSE,quiet = FALSE)
ghdf_approach <- RRASSLER::process_ras_hdf_to_xyz(geom_path = paste0(root_path,".hdf"),units = current_model_units,proj_string = current_model_projection,vdat = FALSE,quiet = FALSE)

AOI_lines <- sfheaders::sf_linestring(obj = database_approach[[1]], x = "x", y = "y", linestring_id = "xid", keep = FALSE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) 

desired_xs <- 3

# Step 1 plots
# A plot
g_approach_streamline_p <- sf::st_transform(g_approach[[3]],sf::st_crs("EPSG:4326"))
g_approach_pts <- g_approach[[1]]
g_approach_pts_p <- sfheaders::sf_point(obj = g_approach_pts, x = "x", y = "y", z = "z", keep = TRUE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% sf::st_transform(sf::st_crs("EPSG:4326"))
g_approach_transects <- sfheaders::sf_linestring(obj = g_approach_pts, x = "x", y = "y", linestring_id = "xid", keep = FALSE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) 
g_approach_transects_p <- sf::st_transform(g_approach_transects,sf::st_crs("EPSG:4326"))
  
bm_fullmod <- basemaps::basemap_terra(ext = sf::st_bbox(sf::st_buffer(g_approach[[3]],200)),map_service = "esri", map_type = "world_imagery", verbose = FALSE) %>% as("SpatRaster")
bm_fullmod_light <- bm_fullmod + 20

# A plot
g_approach_streamline_p <- sf::st_transform(g_approach[[3]],sf::st_crs("EPSG:4326"))
g_approach_pts <- g_approach[[1]]
g_approach_pts_p <- sfheaders::sf_point(obj = g_approach_pts, x = "x", y = "y", z = "z", keep = TRUE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% sf::st_transform(sf::st_crs("EPSG:4326"))
g_approach_transects <- sfheaders::sf_linestring(obj = g_approach_pts, x = "x", y = "y", linestring_id = "xid", keep = FALSE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) 
g_approach_transects_p <- sf::st_transform(g_approach_transects,sf::st_crs("EPSG:4326"))

## Geo
full_g_geo <- ggplot2::ggplot() +
  tidyterra::geom_spatraster_rgb(data = bm_fullmod_light) +
  ggplot2::geom_sf(data = g_approach_streamline_p ,colour = "#ffffff", linewidth = 2.5, arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  ggplot2::geom_sf(data = g_approach_streamline_p, color = "#3182bd", arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  ggplot2::geom_sf(data = g_approach_transects_p[g_approach_transects_p$xid==desired_xs,] ,colour = "#ffffff", linewidth = 2.5, arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  ggplot2::geom_sf(data = g_approach_transects_p, color = "red", arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm")))  +
  cowplot::theme_half_open() +
  ggplot2::coord_sf(expand = FALSE) +
  ggplot2::scale_x_continuous(breaks = seq(from = sf::st_bbox(g_approach_transects_p)$xmin, to = sf::st_bbox(g_approach_transects_p)$xmax, by = (sf::st_bbox(g_approach_transects_p)$xmax-sf::st_bbox(g_approach_transects_p)$xmin)/2))

## Single XS
g_target_xs <- g_approach_transects_p[g_approach_transects_p$xid==desired_xs,]
g_target_xs_pts <- g_approach_pts[g_approach_pts$xid==g_target_xs$xid,]
g_target_xs_pts_p <- sfheaders::sf_point(obj = g_target_xs_pts, x = "x", y = "y", z = "z", keep = TRUE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% sf::st_transform(sf::st_crs("EPSG:4326"))
bm <- basemaps::basemap_terra(ext = sf::st_bbox(sf::st_buffer(g_target_xs,200)),map_service = "esri", map_type = "world_imagery", verbose = FALSE) %>% as("SpatRaster") # mapview::mapview(sf::st_buffer(g_target_xs,200))
g_target_xs_pts_p$ras_z <- g_target_xs_pts$z


# B plot
ghdf_approach_streamline_p <- sf::st_transform(ghdf_approach[[3]],sf::st_crs("EPSG:4326"))
ghdf_approach_pts <- ghdf_approach[[1]]
ghdf_approach_pts_p <- sfheaders::sf_point(obj = ghdf_approach_pts, x = "x", y = "y", z = "z", keep = TRUE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% sf::st_transform(sf::st_crs("EPSG:4326"))
ghdf_approach_transects <- sfheaders::sf_linestring(obj = ghdf_approach_pts, x = "x", y = "y", linestring_id = "xid", keep = FALSE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) 
ghdf_approach_transects_p <- sf::st_transform(ghdf_approach_transects,sf::st_crs("EPSG:4326"))
bm <- basemaps::basemap_terra(ext = sf::st_bbox(sf::st_buffer(ghdf_approach[[3]],200)),map_service = "esri", map_type = "world_imagery", verbose = FALSE) %>% as("SpatRaster")

## Geo
full_ghdf_geo <- ggplot2::ggplot() +
  tidyterra::geom_spatraster_rgb(data = bm_fullmod_light) +
  ggplot2::geom_sf(data = ghdf_approach_streamline_p ,colour = "#ffffff", linewidth = 2.5, arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  ggplot2::geom_sf(data = ghdf_approach_streamline_p, color = "#3182bd", arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  ggplot2::geom_sf(data = ghdf_approach_transects_p[ghdf_approach_transects_p$xid==desired_xs,] ,colour = "#ffffff", linewidth = 2.5, arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  ggplot2::geom_sf(data = ghdf_approach_transects_p, color = "red", arrow = ggplot2::arrow(angle = 45, ends = "last",type = "open", length = ggplot2::unit(0.5, "cm"))) +
  cowplot::theme_half_open() + 
  ggplot2::coord_sf(expand = FALSE) +
  ggplot2::scale_x_continuous(breaks = seq(from = sf::st_bbox(ghdf_approach_transects_p)$xmin, to = sf::st_bbox(ghdf_approach_transects_p)$xmax, by = (sf::st_bbox(ghdf_approach_transects_p)$xmax-sf::st_bbox(ghdf_approach_transects_p)$xmin)/2))


## Single XS
ghdf_target_xs <- ghdf_approach_transects_p[ghdf_approach_transects_p$xid==desired_xs,]
ghdf_target_xs_pts <- ghdf_approach_pts[ghdf_approach_pts$xid==ghdf_target_xs$xid,]
ghdf_target_xs_pts_p <- sfheaders::sf_point(obj = ghdf_target_xs_pts, x = "x", y = "y", z = "z", keep = TRUE) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% sf::st_transform(sf::st_crs("EPSG:4326"))
bm <- basemaps::basemap_terra(ext = sf::st_bbox(sf::st_buffer(ghdf_target_xs,200)),map_service = "esri", map_type = "world_imagery", verbose = FALSE) %>% as("SpatRaster") # mapview::mapview(sf::st_buffer(g_target_xs,200))
ghdf_target_xs_pts_p$ras_z <- ghdf_target_xs_pts$z

# Combined
full_extract_comp <- ggplot2::ggplot() +
  ggplot2::geom_point(data = g_target_xs_pts_p, ggplot2::aes(x = xid_d, y = ras_z, color = 'g file')) +
  ggplot2::geom_point(data = ghdf_target_xs_pts_p, ggplot2::aes(x = xid_d, y = ras_z, color = 'g.hdf file')) +
  ggplot2::geom_line(data = g_target_xs_pts_p, ggplot2::aes(x = xid_d, y = ras_z, color = 'g file')) +
  ggplot2::geom_line(data = ghdf_target_xs_pts_p, ggplot2::aes(x = xid_d, y = ras_z,color = 'g.hdf file')) +
  ggplot2::xlab("Station (meters)") + ggplot2::ylab("Elevation (meters)") +
  ggplot2::labs(title = "RRASSLER read cross sections", subtitle = glue::glue("Cross section {desired_xs} from the two extractions"), family = "serif") +
  cowplot::theme_half_open() +
  cowplot::background_grid() +
  ggplot2::geom_text(data = g_target_xs_pts_p, x=max(g_target_xs_pts_p$xid_d)/2, y=max(g_target_xs_pts_p$ras_z), 
                     label = glue::glue("Stats"), 
                     family = "sans serif", size=6, color="#000000") +
  ggplot2::geom_text(data = g_target_xs_pts_p, x=max(g_target_xs_pts_p$xid_d)/2, y=max(g_target_xs_pts_p$ras_z)-1, 
                     label = glue::glue("g file: n points:{nrow(g_approach_pts_p)} | sum of length of xs: {sum(sf::st_length(g_approach_transects_p))}"), 
                     family = "sans serif", size=4,color=scico::scico(6, palette = 'bamako')[3]) +
  ggplot2::geom_text(data = g_target_xs_pts_p, x=max(g_target_xs_pts_p$xid_d)/2, y=max(g_target_xs_pts_p$ras_z)-2, 
                     label = glue::glue("g.hdf file: n points:{nrow(ghdf_approach_pts_p)} | sum of length of xs: {sum(sf::st_length(ghdf_approach_transects_p))}"), 
                     family = "sans serif", size=4,color=scico::scico(6, palette = 'bilbao')[5]) +
  ggplot2::scale_color_manual(name='Extracted data',breaks=c('g file', 'g.hdf file'),values=c('g file'=scico::scico(6, palette = 'bamako')[3], 'g.hdf file'= scico::scico(6, palette = 'bilbao')[5])) +
  ggplot2::theme(legend.position = c(.92, .15),legend.justification = c("right", "top"),legend.box.just = "right",legend.margin = ggplot2::margin(6, 6, 6, 6),legend.background = ggplot2::element_rect(fill = "lightgrey"))

Step 2) Verify the conflation data

Lets also verify that we have a domain to conflate to:

# 1. Ensure the hull is in the same CRS as the file (usually EPSG:4326 or 5070 for NWM)
# Assuming nwm_flows.gpkg is in EPSG:5070 as per your conus_crs variable:
ras_hull_5070 <- sf::st_transform(ras_hull, 5070)

# 2. Convert to WKT string
hull_wkt <- sf::st_as_text(sf::st_geometry(ras_hull_5070))

# 3. Read only the intersecting features
nwm_subset <- sf::st_read(
  file.path(ras2fim_national_datasets, "nwm_flows.gpkg"),
  wkt_filter = hull_wkt,
  quiet = !is_verbose
)

ras2fim debugging

Step 1) Verify the model

Step 2) Verify the conflation data

Test V1

Test V2