diff --git a/exploration/NEA_calibration_offline/distance_inshore_wintering/distanceMatrixNEA.Rmd b/exploration/NEA_calibration_offline/distance_inshore_wintering/distanceMatrixNEA.Rmd index 97bc6bad4ec50cca811c635f4bbc18b1086eacee..bbb76c599f1153864ac083aef3e3755f41f72ac4 100644 --- a/exploration/NEA_calibration_offline/distance_inshore_wintering/distanceMatrixNEA.Rmd +++ b/exploration/NEA_calibration_offline/distance_inshore_wintering/distanceMatrixNEA.Rmd @@ -25,6 +25,8 @@ library(dplyr) library(tidyr) library(flextable) + +library(units) ``` The computation is based on the depth raster provided by the getNOAA.bathy function with a resolution of 2° x 2°. @@ -45,17 +47,17 @@ rm(list = ls()) # plot(cost_bathy, image = TRUE, land = TRUE, lwd = 0.1) # # # == cost matrix WITH depth preferundum ====== -# # cost_bathy_pref[cost_bathy_pref < 0] = 0 -# bath_inf = -75 -# bath_sup = -25 - -# informations relatives aux jeunes stades/adultes en mer dans la section D du rapport ASMFC sur l'alose américaine (voir PJ - page 33 du rapport ou 52 du PDF). -# Cette partie est divisée en plusieurs rubriques, qui regroupe les informations sur les préférédums environnementaux (i.e. une section dédiées aux températures, l'autre à la salinité, ect.). -# La partie concernant les profondeurs figure p35. -# Tu as un tableau récapitulatif qui reprend l'ensemble des facteurs influençant la distribution en mer, avec les ranges spécifiques à la p38 du rapport (ou 57 du PDF). +# Greene, Karen E., Jennifer L. Zimmerman, R. Wilson Laney, et Jessie C. Thomas-Blate. « Atlantic coast +# diadromous fish habitat: a review of utilization, threats, recommendations for conservation, and research +# needs ». Atlantic States Marine Fisheries Commission Habitat Management Series 464 (2009): 276. +# see table page 39 + +# bath_inf = -100 +# bath_sup = -50 # # cost_bathy_pref = grasp_bathy2 -# cost_bathy_pref[grasp_bathy2 <= 0 ] = 10 +# cost_bathy_pref[grasp_bathy2 >= 0] = 1e12 +# cost_bathy_pref[grasp_bathy2 < 0 ] = 10 # cost_bathy_pref[grasp_bathy2 >= bath_inf & grasp_bathy2 <= bath_sup] = 0 # # Replace by a very high value all cells located on land # cost_bathy_pref[grasp_bathy2 > 0] = 1e12 @@ -63,6 +65,7 @@ rm(list = ls()) # plot(cost_bathy_pref, image=TRUE, land=TRUE, lwd=0.1, main ='pref') # # trans <- trans.mat(cost_bathy) + # trans_pref <- trans.mat(cost_bathy_pref) # # save # save(list = ls(all.names = TRUE), file = "bathy_map.RData") @@ -72,7 +75,7 @@ load(file = "bathy_map.RData") ``` ```{r} -# connexion à la bdd =============== +# # connexion à la bdd =============== # getPass = function() { # print("password: ") # pass = scan(n = 1, @@ -98,99 +101,145 @@ load(file = "bathy_map.RData") # # # load sea outlet # sea_outlets = dbGetQuery(conn, -# 'SELECT seaoutlet_id, x, y FROM northeastamerica.seaoutlets; ') +# 'select +# distinct s.seaoutlet2_id, +# s.x, +# s.y +# from +# northeastamerica.seaoutlets s +# inner join northeastamerica.basin_outlet bo +# using (seaoutlet2_id) +# inner join northeastamerica.v_river_basin_to_gr3d vrbtgd +# using(basin_id) +# order by +# y;') # -# # the list is ordered with ordre +# # the list is ordered with sea outlet latitude # conversionTable = dbGetQuery( # conn, -# 'SELECT gid, -# B.name, -# st_length(seapath::geography)/1000 as seapath_length, -# seaoutlet_id, ordre -# FROM northeastamerica.basin_outlet B -# INNER JOIN northeastamerica.seaoutlets S on (st_endpoint(B.seapath) = S.seaoutlet_geom) -# inner join northeastamerica.v_nea_basins using (gid) -# ORDER BY ordre' +# 'select +# basin_id, +# bo.basin_name, +# st_length(seapath::geography)/1000 as seapath_length, +# seaoutlet2_id +# from +# northeastamerica.basin_outlet bo +# inner join northeastamerica.v_river_basin_to_gr3d vrbtgd +# using(basin_id) +# order by +# st_y(seaoutlet_geom)' # ) # -# distanceSeaOutlet_wintering = dbGetQuery( +# combinationsSeaOutlet_wintering = dbGetQuery( # conn, -# "select -# seaoutlet_id, +# "with part1 as +# (select +# distinct s.* +# from +# northeastamerica.seaoutlets s +# inner join northeastamerica.basin_outlet bo +# using (seaoutlet2_id) +# inner join northeastamerica.v_river_basin_to_gr3d vrbtgd +# using (basin_id)) +# select +# seaoutlet2_id, # basin_name, +# function, # x as x_seaoutlet, # y as y_seaoutlet, # st_x(centroid_geom) as x_winter_offshore, # st_y(centroid_geom) as y_winter_offshore # from -# northeastamerica.seaoutlets, +# part1, # northeastamerica.offshore_basins ob -# where type = 'WinterBasin';" -# ) +# where function = 'WINTERING'; ") # +# # # inshore_offshore_connections = dbGetQuery( # conn, # "select -# * +# ioc.* # from -# northeastamerica.inshore_offshore_connections;" +# northeastamerica.inshore_offshore_connections ioc +# inner join northeastamerica.v_river_basin_to_gr3d vrbtgd on +# ioc.inshore_basin_id = vrbtgd.basin_id;" # ) # +# +# # save(sea_outlets, # conversionTable, -# distanceSeaOutlet_wintering, +# combinationsSeaOutlet_wintering, # inshore_offshore_connections, # file = "outlet.RData") + load("outlet.RData") ``` -The distances are calculated as the shortest path avoiding land. +# Distance between inshore and wintering offshore basins + +The distances are calculated as the shortest path avoiding land. In that case the distance computed by the lc.dict take into account the weights used to define the preferenda. To have the actual distance is obtainad with path. ```{r} -computeDistancebetweenTwoPoint = function(data, trans, res = 'dist') { +computePathBetweenTwoPoints = function(data, trans, crs = 4326) { resultat = sapply(1:nrow(data), function(i) { - lc.dist( + + slicedData <- data %>% slice(i) + path = lc.dist( trans, - data %>% slice(i) %>% select(x = x_seaoutlet, y = y_seaoutlet) %>% - bind_rows(data %>% slice(i) %>% select(x = x_winter_offshore , y = y_winter_offshore)) , - res = res - ) - }) - if (res == 'dist') - output <- data %>% select(seaoutlet_id , basin_name) %>% mutate(distance = resultat %>% unlist(use.names = FALSE)) - else - output <- data %>% select(seaoutlet_id , basin_name) %>% mutate(path = resultat) - return(output) + slicedData %>% select(x = x_departure, y = y_departure) %>% + bind_rows(slicedData %>% select(x = x_arrival , y = y_arrival)) , + res = 'path') + + + # add departure and arrival points + extendedPath <- slicedData %>% select(x = x_departure, y = y_departure) %>% + bind_rows(data.frame(x = path[[1]][,1], y = path[[1]][,2])) %>% + bind_rows(slicedData %>% select(x = x_arrival, y = y_arrival)) + class(list(extendedPath)) + + # distance in km + distance = st_length(st_sfc(st_multilinestring(list(as.matrix(extendedPath))), crs = 4326)) + units(distance) <- make_units(km) + + return(list(data = slicedData %>% mutate(distance = distance), path = extendedPath)) + }, simplify = FALSE) + + return(resultat) } -# distanceToWinteringOffshore <- computeDistancebetweenTwoPoint(data = distanceSeaOutlet_wintering, trans = trans, res ='dist' ) -#pathToWinteringOffshore <- computeDistancebetweenTwoPoint(data = distanceSeaOutlet_wintering, trans = trans, res = 'path' ) +getDistanceResult = function(dataResult){ + do.call(rbind.data.frame, lapply(seq_along(dataResult), function(x, data) {data[[x]]$data}, data = dataResult)) +} + + +pathToWinteringOffshore <- computePathBetweenTwoPoints(data = combinationsSeaOutlet_wintering %>% + rename(x_departure = x_seaoutlet, + y_departure = y_seaoutlet, + x_arrival = x_winter_offshore, + y_arrival = y_winter_offshore), + trans = trans_pref) -#save(distanceToWinteringOffshore, pathToWinteringOffshore, file = 'toWinteringOffshore.RData') +save(pathToWinteringOffshore, file = 'toWinteringOffshore.RData') load('toWinteringOffshore.RData') +distanceToWinteringOffshore <- getDistanceResult(pathToWinteringOffshore) ``` -```{r} -# verification des distances -# pathToWinteringOffshore %>% mutate(distance = -# lapply(pathToWinteringOffshore[,'path'], function(w) st_length(st_sfc(st_multilinestring(list(w)), crs = 4326))) %>% unlist()/1000) %>% -# inner_join(distanceToWinteringOffshore, by =c('seaoutlet_id', 'basin_name' )) - -``` -Distances between sea outlets (which are connected to river basin) and the 3 wintering offshore basins are then calculated. +Distances between sea outlets (which are connected to river basin) and the 3 wintering offshore basins are then calculated considering a depth preferendum between `r -bath_inf` and `r -bath_sup` m. ```{r distanceToWintering, include=TRUE, fig.cap='Distance(in km) to the wintering offshore basins.'} -conversionTable %>% select(name, seaoutlet_id) %>% + +conversionTable %>% select(basin_name, seaoutlet2_id) %>% inner_join( - distanceToWinteringOffshore %>% + distanceToWinteringOffshore%>% pivot_wider(names_from = basin_name, values_from = distance), - by = c('seaoutlet_id') + by = c('seaoutlet2_id') ) %>% flextable() %>% autofit() ``` @@ -198,16 +247,16 @@ conversionTable %>% select(name, seaoutlet_id) %>% ```{r} shortestPathToWinteringOffshore <- pathToWinteringOffshore %>% inner_join(distanceToWinteringOffshore %>% - group_by(seaoutlet_id) %>% - arrange(distance) %>% - mutate(rank = order(distance)) %>% - filter(rank == 1), by = c("seaoutlet_id", "basin_name")) + group_by(seaoutlet2_id) %>% + arrange(distance) %>% + mutate(rank = order(distance)) %>% + filter(rank == 1), by = c("seaoutlet_id", "basin_name")) ``` -By selecting the shortest distance, it is possible to verify the consistency of GR3D connections with this criterion. Presently there are some discrepancies near Cape Fear ! +By selecting the shortest distance, GR3D connections is consistent with this criterion. ```{r winterinbBasinconsistency, include=TRUE, tab.cap="Consistency between GR3D wintering basin and selection of wintering basin based on minimal distance."} -conversionTable %>% select(name, seaoutlet_id) %>% +conversionTable %>% select(basin_name, seaoutlet2_id) %>% inner_join(shortestPathToWinteringOffshore %>% select(seaoutlet_id, closest_wintering_basin = basin_name), by = c('seaoutlet_id')) %>% inner_join(inshore_offshore_connections %>% select(inshore_basin_name, wintering_offshore_basin_name), @@ -217,13 +266,57 @@ conversionTable %>% select(name, seaoutlet_id) %>% autofit() ``` -A map to illustrate the paths. Notice that no depth prefredum are considered. +A map to illustrate the paths. ```{r pathTowinteringMap, include=TRUE} -plot(cost_bathy, image=TRUE, land=TRUE, lwd=0.1, main = 'without pref') +plot(cost_bathy_pref, image=TRUE, land=TRUE, lwd=0.1, main = 'without pref') rien = lapply(shortestPathToWinteringOffshore[,'path'], function(w) lines(w, col = 'red')) points(sea_outlets[,c('x', 'y')], col = 'green', pch = 16, cex = .5) ``` +# Distances between river basins + +```{r} +# dataForPathBetweenSeaOutlets <- expand_grid(departure_row_nb = 1:nrow(sea_outlets), arrival_row_nb= 1:nrow(sea_outlets)) %>% +# filter(departure_row_nb < arrival_row_nb) %>% +# inner_join(sea_outlets %>% +# mutate(row_nb = row_number()) %>% +# select(row_nb, departure_seaoutlet2_id = seaoutlet2_id, x_departure = x, y_departure = y), +# by = c('departure_row_nb' = 'row_nb')) %>% +# inner_join(sea_outlets %>% +# mutate(row_nb = row_number()) %>% +# select(row_nb, arrival_seaoutlet2_id = seaoutlet2_id, x_arrival = x, y_arrival = y ), +# by = c('arrival_row_nb' = 'row_nb') ) %>% +# select(departure_seaoutlet2_id, arrival_seaoutlet2_id, x_departure, y_departure, arrival_seaoutlet2_id, x_arrival, y_arrival) +# +# pathBetweenSeaOutlets = pathBetweenSeaOulets +# pathBetweenSeaOutlets = computePathBetweenTwoPoints(dataForPathBetweenSeaOutlets, trans_pref) +# +# save(pathBetweenSeaOutlets, file = 'betweenSeaOutletsWithPref.RData') +load('betweenSeaOutletsWithPref.RData') +distanceBetweenSeaOulets <- getDistanceResult(pathBetweenSeaOulets) + +``` + + +```{r} +plot(cost_bathy_pref, image=TRUE, land=TRUE, lwd=0.1, main = 'with pref') +lapply(pathBetweenSeaOulets, function(w) lines(w, col ='red')) +points( sea_outlets[,c('x', 'y')], col = 'green', pch = 16, cex = .5) +``` + +```{r} +expand_grid(departure_id = conversionTable$basin_id, arrival_id = conversionTable$basin_id) %>% + inner_join(conversionTable %>% select(basin_id, departure_seaoutlet2_id = seaoutlet2_id), + by = c('departure_id' = 'basin_id')) %>% + inner_join(conversionTable %>% select(basin_id, arrival_seaoutlet2_id = seaoutlet2_id), + by = c('arrival_id' = 'basin_id')) + + + + + + +```