Global velocities

rm(list=ls())
gc()

##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 528699 28.3    1180377 63.1   660382 35.3
## Vcells 964835  7.4    8388608 64.0  1769918 13.6

# devtools::install_github("GMBA-biodiversity/gmbaR")

list.of.packages <- c("terra","ggplot2","tidyverse","MultiscaleDTM","forcats")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
sapply(list.of.packages, require, character.only = TRUE)

##         terra       ggplot2     tidyverse MultiscaleDTM       forcats 
##          TRUE          TRUE          TRUE          TRUE          TRUE

# set computer
computer = "personal"

if(computer == "muse"){
    setwd("/storage/simple/projects/t_cesab/brunno/Exposure-SDM")
    work_dir <- getwd()
}
if(computer == "matrics"){
    setwd("/users/boliveira/Exposure-SDM")
    work_dir <- getwd()
}
work_dir <- here::here()

source(here::here(work_dir,"R/velocity_functions.R"))

### function to get equal binds from a vector
func_splint <- function(x,interval=4) {
    require(ggplot2)
    is.odd <- function(x) x %% 2 != 0
    
    a <- levels(cut_interval(x,interval-1))
    b <- unlist(strsplit(a,','))
    c <- gsub('[','',b,fixed="TRUE")
    d <- gsub(']','',c,fixed="TRUE")
    e <- gsub('(','',d,fixed="TRUE")
    f <- gsub(')','',e,fixed="TRUE")
    return(as.numeric(c(unique(f))))
}

# Global continental lines
globe_lines <- vect(rnaturalearth::ne_countries(scale = "small"))
globe_lines <- terra::aggregate(globe_lines)

# Mountains
if(!file.exists(here::here("Data/GMBA_Inventory_v2.0_standard.zip"))){
    # download
    download.file("https://data.earthenv.org/mountains/standard/GMBA_Inventory_v2.0_standard.zip",
                  here::here("Data/GMBA_Inventory_v2.0_standard.zip"))
    # unzip
    unzip(here::here("Data/GMBA_Inventory_v2.0_standard.zip"),
          exdir = here::here("Data/GMBA_Inventory_v2.0_standard"))
    mtn <- vect(here::here("Data/GMBA_Inventory_v2.0_standard/GMBA_Inventory_v2.0_standard.shp"))
} else {
    mtn <- vect(here::here("Data/GMBA_Inventory_v2.0_standard/GMBA_Inventory_v2.0_standard.shp"))
}
theandes <- mtn[mtn$Level_02 == "Andes"]

# Mountains slope
terrain_slope <- rast(here::here("Data/elevation_slp_1km.tif"))

1 Load velocity layers

##############################
# Marine
## velocity
mar_vel <- rast(here::here(work_dir,"Data/Velocity_global/Mar_mean_gVel_1960-2009.tif"))
## velocity latitude
mar_vel_lat <- rast(here::here(work_dir,"Data/Velocity_global/Mar_mean_gVelLat_1960-2009.tif"))
## velocity angle
mar_vel_angle <- terra::rast(here::here("Data/Velocity_global/Mar_mean_gVelAngle_1960-2009.tif"))
# mar_vel_angle <- (mar_vel_angle + 180) %% 360

##############################
# Terrestrial 25km
## velocity temperature
ter_vel_mat_25km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_mat_gVel_25km_1960-2009.tif"))
## velocity temperature latitude
ter_vel_mat_lat_25km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_mat_gVelLat_25km_1960-2009.tif"))
## velocity temperature angle
ter_vel_mat_angle_25km <- rast(here::here("Data/Velocity_global/Ter_mat_gVelAngle_25km_1960-2009.tif"))

## velocity precipitation
ter_vel_map_25km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_map_gVel_25km_1960-2009.tif"))
## velocity precipitation latitude
ter_vel_map_lat_25km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_map_gVelLat_25km_1960-2009.tif"))
## velocity precipitation angle
ter_vel_map_angle_25km <- rast(here::here("Data/Velocity_global/Ter_map_gVelAngle_25km_1960-2009.tif"))

## velocity temperature elevation 
ter_vel_mat_ele_1km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_mat_gVelEle_1km_1960-2009.tif"))

##############################
# Terrestrial 1km
## velocity temperature
ter_vel_mat_1km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_mat_gVel_1km_1960-2009.tif"))
## velocity temperature latitude
ter_vel_mat_lat_1km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_mat_gVelLat_1km_1960-2009.tif"))
## velocity temperature angle
ter_vel_mat_angle_1km <- rast(here::here("Data/Velocity_global/Ter_mat_gVelAngle_1km_1960-2009.tif"))

## velocity precipitation
ter_vel_map_1km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_map_gVel_1km_1960-2009.tif"))
## velocity precipitation latitude
ter_vel_map_lat_1km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_map_gVelLat_1km_1960-2009.tif"))
## velocity precipitation angle
ter_vel_map_angle_1km <- rast(here::here("Data/Velocity_global/Ter_map_gVelAngle_1km_1960-2009.tif"))

## velocity precipitation elevation 
ter_vel_map_ele_1km <- rast(here::here(work_dir,"Data/Velocity_global/Ter_map_gVelEle_1km_1960-2009.tif"))

2 Plot Marine velocities

2.1 Unprojected

mar_velrange <- terra::minmax(mar_vel)

the_palette_fc <- leaflet::colorNumeric( 
    palette = "RdBu", 
    domain = c(-max(abs(mar_velrange)),max(abs(mar_velrange))),
    reverse = TRUE)

the_colors <- the_palette_fc(seq(min(mar_velrange), max(mar_velrange), length.out = 15))

plot(mar_vel, col = the_colors, main = "Temperature velocity")

2.2 Latitude

mar_vel_latrange <-  terra::minmax(mar_vel_lat)

the_palette_fc <- leaflet::colorNumeric(
    palette = "RdBu", 
    domain = c(-max(abs(mar_vel_latrange)),max(abs(mar_vel_latrange))),
    reverse = TRUE)

the_colors <- the_palette_fc(seq(min(mar_vel_latrange), max(mar_vel_latrange), length.out = 50))

plot(mar_vel_lat, col = the_colors, main = "Temperature velocity latitude")

2.3 Angle

myramp <- colorRampPalette(colors = c("red","green","blue","yellow","red"))(360)

{
    layout(matrix(c(1,1,1,2), nrow = 1, ncol = 4, byrow=T))
    
    plot(mar_vel_angle, 
         col = myramp)
    
    plotrix::polar.plot(
        start = 90,
        lengths = c(rnorm(360,mean = 1, sd = 0.001)),
        polar.pos = seq(0,360,by=1),
        clockwise = TRUE,
        cex=.01,
        show.grid.labels=0,
        line.col=myramp)
}

2.4 Difference lat vs Unprojected

new_rast <- abs(mar_vel) - abs(mar_vel_lat) 

new_rastrange <-  terra::minmax(new_rast)

the_palette_fc <- leaflet::colorNumeric(
    palette = "RdBu", 
    domain = c(-max(abs(new_rastrange)),max(abs(new_rastrange))),
    reverse = TRUE)

the_colors <- the_palette_fc(seq(min(new_rastrange), max(new_rastrange), length.out = 50))

plot(new_rast, col = the_colors, main = "Difference velocities\nUnprojected - Latitude")

plot(abs(mar_vel_lat),
     abs(mar_vel),
     xlab="Velocity Unprojected",
     ylab = "Velocity latitude (absolut)")

The difference between velocities latitude and unprojected is overall small, suggesting most of the Temperature velocity goes in the trajectory North. Areas highlighted in blue indicate faster velocities in directions away from the North angle (i.e., west or east).

3 Plot Terrestrial velocities

3.1 Temperature

3.1.1 Unprojected

velocity_map(ter_vel_mat_1km, main = "Temperature velocity (1km)")

velocity_map(ter_vel_mat_25km, main = "Temperature velocity (25km)")

hist(ter_vel_mat_1km, main = "Temperature velocity (1km)")

hist(ter_vel_mat_25km, main = "Temperature velocity (25km)")

3.1.2 Latitude

velocity_map(ter_vel_mat_lat_1km, main = "Temperature velocity latitude (1km)")

velocity_map(ter_vel_mat_lat_25km, main = "Temperature velocity latitude (25km)")

hist(ter_vel_mat_lat_1km, main = "Temperature velocity (1km)")

hist(ter_vel_mat_lat_25km, main = "Temperature velocity (25km)")

3.1.3 Angle

angle_map(ter_vel_mat_angle_1km, main = "Temperature velocity angle (1km)")

angle_map(ter_vel_mat_angle_25km, main = "Temperature velocity angle (25km)")

angle_map(crop(ter_vel_mat_angle_1km,c(-90,-60,-40,0)), main = "Temperature velocity angle (1km)")

angle_map(crop(ter_vel_mat_angle_25km,c(-90,-60,-40,0)), main = "Temperature velocity angle (25km)")

3.1.4 Difference lat vs Unprojected

abs_ter_vel_mat_1km <- abs(ter_vel_mat_1km)

## |---------|---------|---------|---------|=========================================                                          

abs_ter_vel_mat_lat_1km <- abs(ter_vel_mat_lat_1km) 

## |---------|---------|---------|---------|=========================================                                          

new_rast <- abs_ter_vel_mat_1km - abs(ter_vel_mat_lat_1km) 

## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

new_rastrange <-  terra::minmax(new_rast)

the_palette_fc <- leaflet::colorNumeric(
    palette = "RdBu", 
    domain = c(-max(abs(new_rastrange)),max(abs(new_rastrange))),
    reverse = TRUE)

the_colors <- the_palette_fc(seq(min(new_rastrange), max(new_rastrange), length.out = 50))

plot(new_rast, col = the_colors, main = "Difference velocities\nUnprojected - Latitude (1km)")

plot(abs_ter_vel_mat_1km,
     abs_ter_vel_mat_lat_1km,
     main = "Difference velocities\nUnprojected - Latitude (1km)",
     xlab="Velocity Unprojected",
     ylab = "Velocity latitude (absolut)")

3.2 Precipitation

3.2.1 Unprojected

velocity_map(ter_vel_map_1km, main = "Precipitation velocity (1km)")

velocity_map(ter_vel_map_25km, main = "Precipitation velocity (25km)")

3.2.2 Latitude

velocity_map(ter_vel_map_lat_1km, main = "Precipitation velocity latitude (1km)")

velocity_map(ter_vel_map_lat_25km, main = "Precipitation velocity latitude (25km)")

3.2.3 Angle

angle_map(ter_vel_map_angle_1km, main = "Precipitation velocity angle (1km)")

angle_map(ter_vel_map_angle_25km, main = "Precipitation velocity angle (25km)")

angle_map(crop(ter_vel_map_angle_1km,c(-90,-60,-40,0)), main = "Precipitation velocity angle (1km)")

angle_map(crop(ter_vel_map_angle_25km,c(-90,-60,-40,0)), main = "Precipitation velocity angle (25km)")

3.2.4 Difference lat vs Unprojected

abs_ter_vel_map_1km <- abs(ter_vel_map_1km)

## |---------|---------|---------|---------|=========================================                                          

abs_ter_vel_map_lat_1km <- abs(ter_vel_map_lat_1km) 

## |---------|---------|---------|---------|=========================================                                          

new_rast <- abs_ter_vel_map_1km - abs(ter_vel_map_lat_1km) 

## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

new_rastrange <-  terra::minmax(new_rast)

the_palette_fc <- leaflet::colorNumeric(
    palette = "RdBu", 
    domain = c(-max(abs(new_rastrange)),max(abs(new_rastrange))),
    reverse = TRUE)

the_colors <- the_palette_fc(seq(min(new_rastrange), max(new_rastrange), length.out = 50))

plot(new_rast, col = the_colors, main = "Difference velocities\nUnprojected - Latitude (1km)")

plot(abs_ter_vel_map_1km,
     abs_ter_vel_map_lat_1km,
     main = "Difference velocities\nUnprojected - Latitude (1km)",
     xlab="Velocity Unprojected",
     ylab = "Velocity latitude (absolut)")

3.3 Difference precipitation vs temperature

plot(ter_vel_map_lat_1km,
     ter_vel_map_1km,
     xlab = "Precipitation velocity",
     ylab = "Temperature velocity", 
     cex = 1)

new_rast <- ter_vel_map_1km - ter_vel_map_lat_1km

## |---------|---------|---------|---------|=========================================                                          

velocity_map(new_rast, main = "Difference between velocity in temperature and precipitation (1km)\nTemperature - Precipitation")

plot(ter_vel_map_lat_25km,
     ter_vel_map_25km,
     xlab = "Precipitation velocity",
     ylab = "Temperature velocity", 
     cex = 1)

new_rast <- ter_vel_map_25km - ter_vel_map_lat_25km

velocity_map(new_rast, 
             main = "Difference between velocity in temperature and precipitation (25km)\nTemperature - Precipitation")

Velocities in temperature and precipitation are highly decoupled.

3.4 Elevation

3.4.1 Temperature

velocity_map(ter_vel_mat_ele_1km, main = "Temperature velocity elevation (m/year)")

tmp_r <- terra::mask(terra::crop(ter_vel_mat_ele_1km,theandes),theandes)
tmp_r2 <- terra::mask(terra::crop(terrain_slope,theandes),theandes)

{
    par(mfrow=c(1,2))
    velocity_map(tmp_r, 
                 main = "Temperature velocity elevation (m/year)")
    plot(tmp_r2,
         main = "Terrain slope",
         col = rev(terrain.colors(50)))
    par(mfrow=c(1,1))
}

The maximum velocity in elevation is 1m/year.

plot(tmp_r2, tmp_r,
     xlab="Terrain slope",
     ylab="Temperature velocity elevation (m/year)",
     cex=1)

# mapview::mapview(tmp_r)

Velocity in elevation tends to be higher when the terrain is flat.

3.4.2 Precipitation

velocity_map(ter_vel_map_ele_1km, main = "Precipitation velocity elevation (m/year)")

Precipitation velocities are very high in mountains!

tmp_r <- terra::mask(terra::crop(ter_vel_map_ele_1km,theandes),theandes)
tmp_r2 <- terra::mask(terra::crop(terrain_slope,theandes),theandes)

{
    par(mfrow=c(1,2))
    velocity_map(tmp_r, 
                 main = "Precipitation velocity elevation (m/year)")
    plot(tmp_r2,
         main = "Terrain slope",
         col = rev(terrain.colors(50)))
    par(mfrow=c(1,1))
}

plot(tmp_r2, tmp_r,
     xlab="Terrain slope",
     ylab="Precipitation velocity elevation (m/year)",
     cex=1)

# mapview::mapview(tmp_r)

3.4.3 Difference precipitation vs temperature

new_rast <- ter_vel_mat_ele_1km - ter_vel_map_ele_1km

## |---------|---------|---------|---------|=========================================                                          

velocity_map(new_rast, main = "Difference between velocity in temperature and precipitation\nTemperature - Precipitation")

velocity_map(terra::mask(terra::crop(new_rast,theandes),theandes), 
             main = "The Andes")

4 Mountains

Extract velocities in major mountain systems

data_global <- vel_elev_mtn[which(vel_elev_mtn$mtn=="Global"),]

vel_elev_mtn %>%
    dplyr::mutate(mtn = fct_reorder(mtn, vel_mat),
                  color = ifelse(mtn=="Global","a","b")) %>%
    ggplot(aes(x = vel_mat, y = mtn, color=color))+
    ylab("")+
    xlab("Temperature velocity")+
    geom_point()+
    guides(color="none")+
    scale_y_discrete(position = "right")+
    geom_errorbar(aes(xmin=vel_mat-vel_mat_sd, xmax=vel_mat+vel_mat_sd), 
                  width=.2, 
                  position=position_dodge(0.05))+
    theme_bw()

vel_elev_mtn %>%
    dplyr::mutate(mtn = fct_reorder(mtn, vel_mat),
                  color = ifelse(mtn=="Global","a","b")) %>%
    ggplot(aes(x = vel_map, y = mtn, color=color))+
    ylab("")+
    xlab("Precipitation velocity")+
    geom_point()+
    guides(color="none")+
    scale_y_discrete(position = "right")+
    geom_errorbar(aes(xmin=vel_map-vel_map_sd, xmax=vel_map+vel_map_sd), 
                  width=.2, 
                  position=position_dodge(0.05))+
    theme_bw()

vel_elev_mtn %>%
    dplyr::mutate(mtn = fct_reorder(mtn, vel_mat),
                  Mountain = mtn) %>%
    ggplot(aes(x = abs(vel_mat), y = abs(vel_map), color=Mountain))+
    geom_smooth(method = "lm", color = "black")+
    geom_point()+
    xlab("Temperature velocity")+
    ylab("Precipitation velocity")+
    guides(color="none")+
    geom_errorbar(aes(xmin=abs(vel_mat)-vel_mat_sd, xmax=abs(vel_mat)+vel_mat_sd),
                  width=.2, alpha=.5,
                  position=position_dodge(0.05))+
    geom_errorbar(aes(ymin=abs(vel_map)-vel_map_sd, ymax=abs(vel_map)+vel_map_sd),
                  width=.2, alpha=.5,
                  position=position_dodge(0.05))+
    theme_classic()