Commit 21bbbfe0 authored by patrick.lambert's avatar patrick.lambert
Browse files

compute distances with minimum cost

Showing with 157 additions and 64 deletions
+157 -64
......@@ -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'))
```
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment