--- title: "Estimate a small network of GNSS station" output: rmarkdown::html_vignette bibliography: REFERENCES.bib vignette: > %\VignetteIndexEntry{Estimate a small network of GNSS station} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r 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 ) # knitr::opts_chunk$set(eval = FALSE) ``` # Load packages ```{r, warning=FALSE, message=FALSE} library(gmwmx2) library(dplyr) library(raster) # library(rnaturalearth) library(geodata) library(shape) library(tibble) library(tidygeocoder) library(sf) ``` # Define some functions for plotting ```{r} # 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 ```{r} # Estimate little network in France all_station <- download_all_stations_ngl() # download selected stations selected_station <- c("BSCN", "CERN", "SCDA", "GLRA", "STPS") selected_station <- c("BSCN") # selected_station = c("BSCN") 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", stochastic_model = "wn + pl") fit_E <- gmwmx2(station_data, n_seasonal = 2, component = "E", stochastic_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 and create raster for plot ```{r} # 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 ) ``` # Plot estimated N-E velocities and associated uncertainty ```{r} # sett 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" ) ```