## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( fig.width = 10, # Set default plot width (adjust as needed) fig.height = 8, # Set default plot height (adjust as needed) fig.align = "center" # Center align all plots ) # so as not to execute this vignette during package development knitr::opts_chunk$set(eval = FALSE) ## ----warning=FALSE, message=FALSE--------------------------------------------- # library(gmwmx2) # library(dplyr) # library(raster) # # library(rnaturalearth) # library(geodata) # library(shape) # library(tibble) # library(tidygeocoder) # library(sf) ## ----------------------------------------------------------------------------- # # define some function for plotting # ellipse <- function(hlaxa = 1, hlaxb = 1, theta = 0, xc = 0, yc = 0, newplot = F, npoints = 100, fill = F, fillColor = "black", ...) { # a <- seq(0, 2 * pi, length = npoints + 1) # x <- hlaxa * cos(a) # y <- hlaxb * sin(a) # alpha <- angle(x, y) # rad <- sqrt(x^2 + y^2) # xp <- rad * cos(alpha + theta) + xc # yp <- rad * sin(alpha + theta) + yc # if (newplot) { # plot(xp, yp, type = "l", ...) # } else { # lines(xp, yp, ...) # if (fill == T) { # polygon(xp, yp, border = F, col = fillColor) # } # } # # invisible() # } # # # angle <- function(x, y) { # angle2 <- function(xy) { # x <- xy[1] # y <- xy[2] # if (x > 0) { # atan(y / x) # } else { # if (x < 0 & y != 0) { # atan(y / x) + sign(y) * pi # } else { # if (x < 0 & y == 0) { # pi # } else { # if (y != 0) { # (sign(y) * pi) / 2 # } else { # NA # } # } # } # } # } # apply(cbind(x, y), 1, angle2) # } # # # define function to make color transparent # make_transparent <- function(colors, alpha = 0.5) { # # Ensure alpha is between 0 and 1 # if (alpha < 0 || alpha > 1) { # stop("Alpha value must be between 0 and 1") # } # # # Convert colors to RGB and add alpha # transparent_colors <- sapply(colors, function(col) { # rgb_val <- col2rgb(col) / 255 # rgb(rgb_val[1], rgb_val[2], rgb_val[3], alpha = alpha) # }) # # return(transparent_colors) # } ## ----------------------------------------------------------------------------- # # Estimate a small network in France # all_station <- download_all_stations_ngl() # # # download selected stations # selected_station <- c("BSCN", "CERN", "SCDA", "GLRA", "STPS") # df_network <- all_station %>% filter(station_name %in% selected_station) # df_network # # df_estimated_velocities <- data.frame(matrix(NA, nrow = dim(df_network)[1], ncol = 6)) # for (station_index in seq_along(df_network$station_name)) { # station_name <- df_network$station_name[station_index] # # extract station # station_data <- download_station_ngl(station_name = station_name) # fit_N <- gmwmx2(station_data, n_seasonal = 2, component = "N", model = wn() + pl()) # fit_E <- gmwmx2(station_data, n_seasonal = 2, component = "E", model = wn() + pl()) # df_estimated_velocities[station_index, 1] <- station_name # df_estimated_velocities[station_index, 2:6] <- c(fit_N$beta_hat[2], fit_N$std_beta_hat[2], fit_E$beta_hat[2], fit_E$std_beta_hat[2], dim(fit_N$design_matrix_X)[1]) # # cat(paste0("Processing station ", station_name, " ", station_index, "/", length(df_network$station_name), "\n")) # } # # colnames(df_estimated_velocities) <- c("station_name", "estimated_trend_N", "std_estimated_trend_N", "estimated_trend_E", "std_estimated_trend_E", "time_series_length") # # df_estimated_velocities$estimated_trend_N_scaled <- df_estimated_velocities$estimated_trend_N * 365.25 # df_estimated_velocities$std_estimated_trend_N_scaled <- df_estimated_velocities$std_estimated_trend_N * 365.25 # df_estimated_velocities$estimated_trend_E_scaled <- df_estimated_velocities$estimated_trend_E * 365.25 # df_estimated_velocities$std_estimated_trend_E_scaled <- df_estimated_velocities$std_estimated_trend_E * 365.25 # # # transform longitude # df_network$longitude2 <- ifelse(df_network$longitude < -180, # df_network$longitude + 360, # df_network$longitude # ) # # # merge with location # df_estimated_velocities_and_location <- dplyr::left_join(df_estimated_velocities, df_network, by = "station_name") # # # print estimated North and East velocities # head(df_estimated_velocities) # # knitr::kable(df_estimated_velocities) ## ----------------------------------------------------------------------------- # # load elevation data # tmp <- tempdir() # elevation_data_swiss <- geodata::elevation_30s(country = "Switzerland", path = tmp) # elevation_data_france <- geodata::elevation_30s(country = "France", path = tmp) # elevation_data_italy <- geodata::elevation_30s(country = "Italy", path = tmp) # # # # # # load raster # elevation_raster_swiss <- raster(elevation_data_swiss) # elevation_raster_france <- raster(elevation_data_france) # elevation_raster_italy <- raster(elevation_data_italy) # # # create combined raster # combined_raster <- merge( # elevation_raster_swiss, elevation_raster_france, # elevation_raster_italy # ) ## ----------------------------------------------------------------------------- # # set plot limits # xlims <- c(2, 7) # ylims <- c(44, 48) # # # Custom color scale # custom_colors <- c( # "#33660059", "#33CB6659", "#BAE39159", "#FEDBB859", "#F2C98859", # "#E5B75859", "#D8A52759", "#A7991F59", "#A38F1959", "#A1851359", "#9E7B0D59", "#9B710759", # "#98660059", "#A1595959", "#B1767659", "#B6929259", "#C1AFAF59", "#CBCBCB59", "#E4E4E459", "#FEFEFE59" # ) # # # Define breaks for the color scale based on the raster values # raster_values <- values(combined_raster) # min_val <- min(raster_values, na.rm = TRUE) # max_val <- max(raster_values, na.rm = TRUE) # # # Generate breaks based on the range of the raster data # num_colors <- length(custom_colors) # breaks <- seq(min_val, max_val, length.out = num_colors + 1) # # # # plot # plot(NA, # xlim = xlims, ylim = ylims, las = 1, # ylab = "", xlab = "", # xaxt = "n", yaxt = "n" # ) # # axis(side = 1, at = seq(0, 12, by = 2), labels = (paste0(seq(0, 12, by = 2), " E"))) # axis(side = 2, at = seq(40, 50, by = 2), labels = (paste0(seq(40, 50, by = 2), " N")), las = 1) # # # Plot the elevation data # raster::plot(combined_raster, # col = custom_colors, # # breaks=breaks, # add = TRUE, # Add to the existing plot # legend = FALSE # ) # Disable default legend # # # add axis # for (i in seq(-150, 150, by = 2)) { # abline(v = i, col = "grey80", lty = 5) # } # for (i in seq(-90, 90, by = 2)) { # abline(h = i, col = "grey80", lty = 5) # } # # # add points for station data # points(df_network$longitude2, df_network$latitude, pch = 16) # # # # set param for graph # shift <- 0 # scale_arrow <- 20 # arrow_width <- .1 # arrow_lwd <- 2 # my_col <- c("#e96bff") # scale_ellipses <- 3500 # my_col_trans <- make_transparent(my_col, alpha = .3) # # for (i in seq(nrow(df_estimated_velocities_and_location))) { # ellipse( # hlaxa = as.numeric(df_estimated_velocities_and_location[i, "std_estimated_trend_E_scaled"]) * scale_ellipses, # hlaxb = as.numeric(df_estimated_velocities_and_location[i, "std_estimated_trend_N_scaled"]) * scale_ellipses, # theta = 0, # xc = as.numeric(df_estimated_velocities_and_location[i, "longitude2"]) + as.numeric(df_estimated_velocities_and_location[i, "estimated_trend_E_scaled"]) * scale_arrow, # yc = as.numeric(df_estimated_velocities_and_location[i, "latitude"]) + as.numeric(df_estimated_velocities_and_location[i, "estimated_trend_N_scaled"]) * scale_arrow, # fill = T, # fillColor = my_col_trans[1], # lw = 1, # col = my_col_trans[1] # ) # # x0 <- as.numeric(df_estimated_velocities_and_location[i, "longitude2"]) # y0 <- as.numeric(df_estimated_velocities_and_location[i, "latitude"]) # x1 <- as.numeric(df_estimated_velocities_and_location[i, "longitude2"] + df_estimated_velocities_and_location[i, "estimated_trend_E_scaled"] * scale_arrow) # y1 <- as.numeric(df_estimated_velocities_and_location[i, "latitude"] + df_estimated_velocities_and_location[i, "estimated_trend_N_scaled"] * scale_arrow) # # shape::Arrows( # x0 = x0, y0 = y0, x1 = x1, y1 = y1, # col = my_col, # arr.type = "triangle", # arr.length = .10, # arr.width = arrow_width, # lwd = arrow_lwd # ) # } # # # add # text( # x = df_estimated_velocities_and_location$longitude2, y = df_estimated_velocities_and_location$latitude, # labels = df_estimated_velocities_and_location$station_name, pos = 3, cex = 0.8, col = "black" # ) # # # define cities # df_city <- tibble(address = c("Genève", "Clermont-Ferrand", "Dijon")) # # # # Load geolocalisation of cities # df_geo <- df_city %>% # geocode_combine( # queries = list( # list(method = "osm") # ), # global_params = list(address = "address"), # cascade = FALSE # ) # # df_city_2 <- cbind(df_city, df_geo) # df_city_2$City <- df_city_2$address # # # Add city to map # points(x = df_city_2$lon, y = df_city_2$lat, pch = 15, col = "black") # # cex_size_city <- .7 # for (i in seq(dim(df_city_2)[1])) { # text( # x = df_city_2$lon[i], y = df_city_2$lat[i], # labels = df_city_2$City[i], # adj = c(0, 2), col = "black", # cex = cex_size_city # ) # } # # # # add a legend # twenty_mm_per_year <- .02 # twenty_mm_per_year_mm_per_year_scaled <- twenty_mm_per_year * scale_arrow # x_start <- xlims[2] - 5 # y <- ylims[1] + .3 # segments(x0 = x_start, y0 = y, x1 = x_start + twenty_mm_per_year_mm_per_year_scaled, y1 = y) # delta_y <- .1 # segments(x0 = x_start, x1 = x_start, y0 = y + delta_y, y1 = y - delta_y) # segments(x0 = x_start + twenty_mm_per_year_mm_per_year_scaled, x1 = x_start + twenty_mm_per_year_mm_per_year_scaled, y0 = y + delta_y, y1 = y - delta_y) # txt_size <- .8 # text( # x = x_start + twenty_mm_per_year_mm_per_year_scaled / 2, # y = y + .1, # pos = 3, cex = txt_size, # labels = "20 mm/year" # )