Commit 318e1f5c authored by Poulet Camille's avatar Poulet Camille
Browse files

Merge branch 'develop' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D into develop

parents 1cd6864c 12b5af12
......@@ -28,11 +28,11 @@ inshore_basin_id,inshore_basin_name,wintering_offshore_basin_id,wintering_offsho
11,Santee,1,Off Florida,4,Bay of Fundy
49,Pee Dee,1,Off Florida,4,Bay of Fundy
8,Cape Fear,1,Off Florida,4,Bay of Fundy
42,New,1,Off Florida,4,Bay of Fundy
39,White Oak,1,Off Florida,4,Bay of Fundy
42,New,2,Mid-Atlantic Bight,4,Bay of Fundy
39,White Oak,2,Mid-Atlantic Bight,4,Bay of Fundy
5,Potomac,2,Mid-Atlantic Bight,4,Bay of Fundy
7,Neuse,1,Off Florida,4,Bay of Fundy
48,Pamlico-Tar,1,Off Florida,4,Bay of Fundy
7,Neuse,2,Mid-Atlantic Bight,4,Bay of Fundy
48,Pamlico-Tar,2,Mid-Atlantic Bight,4,Bay of Fundy
45,Rappahannock,2,Mid-Atlantic Bight,4,Bay of Fundy
12,James,2,Mid-Atlantic Bight,4,Bay of Fundy
46,York,2,Mid-Atlantic Bight,4,Bay of Fundy
......
......@@ -18,6 +18,8 @@ library(dplyr)
library(scales)
library(flextable)
rm(list=ls())
fp <- fp_par(
text.align = "center",
padding.bottom = 20, padding.top = 120,
......
library(marmap)
library(RPostgreSQL)
library(DBI) # database interface
library(Rcpp)
library(rpostgis)
library(sf)
rm(list=ls())
# connexion à la bdd ===============
getPass=function() {
print("password: ")
pass=scan(n=1, what=character(), quiet=TRUE)
cat("\014")
return(pass)
}
password = getPass()
host = "citerne.bordeaux.irstea.priv"
user = "patrick.lambert"
dbname = 'eurodiad'
drv <- dbDriver("PostgreSQL")
conn <- dbConnect(drv, host = host, port = 5432, user = user, password = password, dbname = dbname )
# recupère la bathymetrie de la NOAA
grasp_bathy2 = getNOAA.bathy(lon1 = -90, lon2 = -54, lat1 = 20, lat2 = 65, resolution = 2)
# == cost matrix WITHOUT depth preferundum ======
cost_bathy = grasp_bathy2
# Replace by a null value all cells located at sea
cost_bathy[grasp_bathy2 < 0] = 0
# Replace by a very high value all cells located on land
cost_bathy[grasp_bathy2 >= 0] = 1e12
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
cost_bathy_pref = grasp_bathy2
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
# Plot
plot(cost_bathy_pref, image=TRUE, land=TRUE, lwd=0.1, main ='pref')
# load sea outlet
rm(sea_outlets)
sea_outlets = dbGetQuery(conn, 'SELECT seaoutlet_id, x, y FROM northeastamerica.seaoutlets; ')
# the list is ordered with ordre
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')
# matrice des distance
# reshapePointDistanceMatrix = function(path){
# # calc tthe number of points (need to be confirmed)
# nb = 1
# while (sum(1:(nb - 1)) < length(path)) {
# nb = nb + 1
# }
#
# # calculate the length of each patk in km
# dist = unlist(lapply(path, function(w) st_length(st_sfc(st_multilinestring(list(w)), crs = 4269))))/1000
#
# # fill the square matrix
# pointDistanceMatrix = matrix(data = 0, nrow = nb, ncol = nb)
# k = 0
# for (i in 1:nb) {
# for (j in min((i + 1), nb):nb) {
# if (j != i) {
# k = k + 1
# # cat(i,' ', j,' ', k, '\n')
# pointDistanceMatrix[j, i] = dist[[k]]
# pointDistanceMatrix[i, j] = dist[[k]]
# }
# }
# }
# return(pointDistanceMatrix)
# }
# path= path_pref
# path[[1]] =as.matrix(rbind(sea_outlets[1, c('x', 'y')], path[[1]], sea_outlets[2, c('x', 'y')]))
# path[[1]]
# st_multilinestring(x= path[[1]])
reshapePointDistanceMatrix2 = function(path, sea_outlets){
# calc the number of points (need to be confirmed) --> simplify with nrows(sea_outlets)
nb = 1
while (sum(1:(nb - 1)) < length(path)) {
nb = nb + 1
}
# add departure and final points (from sea_outlets) to the path
k = 0
for (i in 1:nb) {
for (j in min((i + 1), nb):nb) {
if (j != i) {
k = k + 1
#cat(i,' ', j,' ', k, '\n')
path[[k]] =as.matrix(rbind(sea_outlets[i, c('x', 'y')], path[[k]], sea_outlets[j, c('x', 'y')]))
}
}
}
# compute the distance for all the paths in km
dist = unlist(lapply(path, function(w) st_length(st_sfc(st_multilinestring(list(w)), crs = 4326))))/1000
# shape the squared matrix
pointDistanceMatrix = matrix(0, nb, nb)
k = 0
for (i in 1:nb) {
for (j in min((i + 1), nb):nb) {
if (j != i) {
k = k + 1
pointDistanceMatrix[j, i] = dist[k]
pointDistanceMatrix[i, j] = dist[k]
}
}
}
return(pointDistanceMatrix)
}
computeMatrixDistance = function(pointDistanceMatrix, conversionTable, catchmentList = conversionTable$name) {
nbCatchment = length(catchmentList)
distanceMatrix = matrix(0, nbCatchment, nbCatchment, dimnames = list(catchmentList, catchmentList))
for (i in 1:nbCatchment) {
for (j in 1:nbCatchment) {
i_point = which(sea_outlets$seaoutlet_id==conversionTable$seaoutlet_id[conversionTable$name == catchmentList[i]])
j_point = which(sea_outlets$seaoutlet_id==conversionTable$seaoutlet_id[conversionTable$name == catchmentList[j]])
distanceMatrix[i, j] = pointDistanceMatrix[i_point, j_point] +
conversionTable$seapath_length[conversionTable$name == catchmentList[i]] +
conversionTable$seapath_length[conversionTable$name == catchmentList[j]]
}
}
return(distanceMatrix)
}
# --======================================================= distance without pref
#get.depth(cost_bathy_pref, x=sea_outlets$x, y=sea_outlets$y, locator = FALSE)
trans <- trans.mat(cost_bathy)
dist <- lc.dist(trans, sea_outlets[,c('x', 'y')], res = "dist")
path <- lc.dist(trans, sea_outlets[,c('x', 'y')], res = "path")
plot(cost_bathy, image=TRUE, land=TRUE, lwd=0.1, main = 'without pref')
lapply(path, function(w) lines(w, col ='red'))
points( sea_outlets[,c('x', 'y')], col = 'green', pch = 16, cex = .5)
pointDistanceMatrix = reshapePointDistanceMatrix2(path,sea_outlets)
outletDistanceMatrix = computeMatrixDistance(pointDistanceMatrix, conversionTable)
save(pointDistanceMatrix, outletDistanceMatrix, file = 'distanceMatrixNEA17march.Rdata' )
# ==================================== distance with pref
# get.depth(cost_bathy_pref, x=sea_outlets$x, y=sea_outlets$y, locator = FALSE)
trans_pref <- trans.mat(cost_bathy_pref)
dist_pref <- lc.dist(trans_pref, sea_outlets[,c('x', 'y')], res = "dist")
path_pref <- lc.dist(trans_pref, sea_outlets[,c('x', 'y')], res = "path")
plot(cost_bathy_pref, image=TRUE, land=TRUE, lwd=0.1, main = 'with pref')
lapply(path_pref, function(w) lines(w, col ='red'))
points( sea_outlets[,c('x', 'y')], col = 'green', pch = 16, cex = .5)
pointDistanceMatrixWithPref = reshapePointDistanceMatrix2(path_pref,sea_outlets)
outletDistanceMatrixWithPref = computeMatrixDistance(pointDistanceMatrixWithPref, conversionTable)
save(pointDistanceMatrixWithPref, outletDistanceMatrixWithPref, file = 'distanceMatrixNEA17marchWithPref.Rdata')
write.csv(outletDistanceMatrixWithPref, file = '/home/patrick/Documents/workspace/GR3D/data/input/northeastamerica/distanceGridNEA.csv')
# distance between sea outlet and centroid of offshore basins
---
date: "`r Sys.Date()`"
author: "P. Lambert, Camille Poulet"
title: "Distance matrix for NEA"
output:
officedown::rdocx_document:
mapstyles:
Normal: ['First Paragraph']
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, include = FALSE, fig.cap = TRUE)
library(officedown)
library(officer)
library(marmap)
library(DBI)
library(RPostgreSQL)
library(sf)
# see library(dbplyr)
library(dplyr)
library(tidyr)
library(flextable)
```
The computation is based on the depth raster provided by the getNOAA.bathy function with a resolution of 2° x 2°.
```{r}
rm(list = ls())
# need to (re)install both rgdal and raster packages.
# # recupère la bathymetrie de la NOAA
# grasp_bathy2 = getNOAA.bathy(lon1 = -90, lon2 = -54, lat1 = 20, lat2 = 65, resolution = 2, keep = TRUE)
#
# # == cost matrix WITHOUT depth preferundum ======
# cost_bathy = grasp_bathy2
# # Replace by a null value all cells located at sea
# cost_bathy[grasp_bathy2 < 0] = 0
# # Replace by a very high value all cells located on land
# cost_bathy[grasp_bathy2 >= 0] = 1e12
# 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).
#
# cost_bathy_pref = grasp_bathy2
# 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
# # Plot
# 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")
load(file = "bathy_map.RData")
```
```{r}
# connexion à la bdd ===============
# getPass = function() {
# print("password: ")
# pass = scan(n = 1,
# what = character(),
# quiet = TRUE)
# cat("\014")
# return(pass)
# }
# password = getPass()
# host = "citerne.bordeaux.irstea.priv"
# user = "patrick.lambert"
# dbname = 'eurodiad'
# drv <- dbDriver("PostgreSQL")
# conn <-
# dbConnect(
# drv,
# host = host,
# port = 5432,
# user = user,
# password = password,
# dbname = dbname
# )
#
# # load sea outlet
# sea_outlets = dbGetQuery(conn,
# 'SELECT seaoutlet_id, x, y FROM northeastamerica.seaoutlets; ')
#
# # the list is ordered with ordre
# 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'
# )
#
# distanceSeaOutlet_wintering = dbGetQuery(
# conn,
# "select
# seaoutlet_id,
# basin_name,
# 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,
# northeastamerica.offshore_basins ob
# where type = 'WinterBasin';"
# )
#
# inshore_offshore_connections = dbGetQuery(
# conn,
# "select
# *
# from
# northeastamerica.inshore_offshore_connections;"
# )
#
# save(sea_outlets,
# conversionTable,
# distanceSeaOutlet_wintering,
# inshore_offshore_connections,
# file = "outlet.RData")
load("outlet.RData")
```
The distances are calculated as the shortest path avoiding land.
```{r}
computeDistancebetweenTwoPoint = function(data, trans, res = 'dist') {
resultat = sapply(1:nrow(data), function(i) {
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)
}
# distanceToWinteringOffshore <- computeDistancebetweenTwoPoint(data = distanceSeaOutlet_wintering, trans = trans, res ='dist' )
#pathToWinteringOffshore <- computeDistancebetweenTwoPoint(data = distanceSeaOutlet_wintering, trans = trans, res = 'path' )
#save(distanceToWinteringOffshore, pathToWinteringOffshore, file = 'toWinteringOffshore.RData')
load('toWinteringOffshore.RData')
```
```{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.
```{r distanceToWintering, include=TRUE, fig.cap='Distance(in km) to the wintering offshore basins.'}
conversionTable %>% select(name, seaoutlet_id) %>%
inner_join(
distanceToWinteringOffshore %>%
pivot_wider(names_from = basin_name, values_from = distance),
by = c('seaoutlet_id')
) %>% flextable() %>%
autofit()
```
```{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"))
```
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 !
```{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) %>%
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),
by=c('name' = 'inshore_basin_name')) %>%
select(name, closest_wintering_basin, wintering_offshore_basin_name ) %>%
flextable() %>%
autofit()
```
A map to illustrate the paths. Notice that no depth prefredum are considered.
```{r pathTowinteringMap, include=TRUE}
plot(cost_bathy, 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)
```
Markdown is supported
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