## MaxEnt SDM for the whitefin swellshark library(spThin) library(DT) library(terra) library(raster) library(sf) library(corrplot) library(usdm) library(kuenm) library(geosphere) library(plotrix) # 'spThin' Spatial Filter for occurrence records -------------------------- ## Import .csv file WFO(Clean) ### spatial thinning of whitefin swellshark occurrence records for use in SDM: thinned_occurrences <- thin(WFO_clean_, lat.col = "decimalLatitude", long.col = "decimalLongitude", spec.col = "Species", thin.par = 5.5, reps = 100, verbose = TRUE, out.dir = "E:/PeerJ (updated)/(1) Occurence Data/(3) 'sp Thin' Data", out.base = "WFO(spThin)", locs.thinned.list.return = TRUE, write.files = TRUE, max.files = 1, write.log.file = FALSE) # Environmental Data (Current) -------------------------------------------- ### Downloaded 13 environmental layers from Bio-ORACLE v3.0 (https://www.bio-oracle.org/) ### Loading Enviro Data into R Current <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Current Environmental Variables", pattern = "\\.nc$", full.names = TRUE) Current <- rast(Current) ## Checking data plot(Current) ## loading substrate type data into R ## note: substrate type was downloaded from (https://data.csiro.au/collection/csiro:12843) ## as a shapefile and reclassified into a raster in QGIS setwd("E:/PeerJ (updated)/(2) Environmental Data/Raw Substrate Data") Substrate <- raster("Substrate.tif") # Cropping enviro Data to Aus --------------------------------------------- ### Cropping the Bio-ORACLE layers to match the substrate type extent Current <- crop(Current, extent(Substrate)) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to stack all layers together Current <- stack(Current) ## Stacking the substrate type layer with the rest of the Bio-ORACLE layers DataMarea <- stack(Current, Substrate) ###check plot(DataMarea) # Masking enviro data to 'M' area ----------------------------------------- ## loading M area shapefile (created in QGIS) into R setwd("E:/PeerJ (updated)/(3) 'M' area/(3) 'M' area (Dissolved)") Marea <- read_sf("'M' area (Dissolved).shp") ### Masking layers to 'M' Area Marea <- mask(DataMarea, Marea) ## Check plot(Marea) # Checking Multi-collinearity ------------------------------------------------------------------- ## Load spatially thinned occurrence data into R setwd("E:/PeerJ (updated)/(1) Occurence Data/(3) 'sp Thin' Data") WFO <- read.csv("WFO(spThin).csv") coordinates(WFO) <- c("decimalLongitude", "decimalLatitude") ### Extract enviro data for each occurrence record WFOextract <- raster::extract(Marea, WFO) ## Convert to dataframe WFOextract <-as.data.frame(WFOextract) ### Check for multicollinearity VIF <- vifcor(WFOextract, th=0.7, method = "pearson") ## See results VIF ### Exclude to keep only non-correlated variables WFEnvriodata <- exclude(Marea, VIF) ####check plot(WFEnvriodata) ### Write to .asc setwd("E:/PeerJ (updated)/(5) Model/WF/M_variables/Set_1") writeRaster(WFEnvriodata, filename = "Current", format = "ascii", bylayer = T, suffix=names(WFEnvriodata), NAFlag = "-9999", overwrite = T) ## Note we also constructed an initial model using the MaxEnt Java software, and ## environmental factors with a small contribution rate (≤ 1%) were as excluded. # Model Calibration ------------------------------------------------------- ### Script adapted from ENM2020: A Free Online Course and Set of Resources on ## Modeling Species' Niches and Distributions (KUENM section) DOI: https://doi.org/10.17161/bi.v17i.15016 setwd("C:/Model/WF") ## Load spatially thinned occurrence data into R WFO_spThin_ <- read_csv("E:/PeerJ (updated)/(1) Occurence Data/(3) 'sp Thin' Data/WFO(spThin).csv") ## Split occurrence data into 70% for model calibration and 30% internal testing kuenm_occsplit(occ = WFO_spThin_, train.proportion = 0.70, method = "random", save = TRUE, name = "WF") oj <- "WF_joint.csv" ## All 145 occurrence records for the whitefin swellshark otr <- "WF_train.csv" ## 70% occurrence data for model calibration mvars <- "M_variables" ## environmental layers regm <- c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4) ## regularisation multipliers to be evaluated fclas <- c("no.t") ## "no.t" excludes the threshold feature class but keeps all other FCs ## FC combinations (l, q, p, h, lq, lp, lh, qp, qh, ph, lqp, lqh, lph, qph, lqph) mxpath <- "C:/Model" ## the path were the maxent.jar file is bcal <- "batch_cal" candir <- "Candidate_models" ## name of the folder that will contain all calibration model subfolders # calibration ## candidate model creation kuenm_cal(occ.joint = oj , occ.tra = otr , M.var.dir = mvars, batch = bcal , out.dir = candir , max.memory = 1000, reg.mult = regm , f.clas = fclas, maxent.path = mxpath , wait = FALSE, run = TRUE, arg = NULL) ## note wait = FLASE to evaluate candidate models parallelly. All other arguments set as default. ## evaluation of candidate models ote <- "WF_test.csv" ## 30% of the occurrence data for internal testing cresdir <- "Calibration_results" ## name of the folder where evaluation results will be written kuenm_ceval(path = candir , occ.joint = oj , occ.tra = otr, occ.test = ote , batch = bcal , out.eval = cresdir, threshold = 5, rand.percent = 50, iterations = 500, kept = FALSE, selection = "OR_AICc", parallel.proc = FALSE) ## Note: threshold = 5, is the percentage of training data omission error allowed (E) default = 5. ## rand.percnet = 50, is the percentage of data to be used for the bootstraping process when calculating partial ROCs; default = 50 ## iterations = 500, is the number of times that the bootstrap is going to be repeated; default = 500 ## selection = "OR_AICc", is the model selection criterion. Default = "OR_AICc", ## which means that among models that are statistically significant and that present omission rates below the threshold, ## those with delta AICc up to 2 will be selected. # Future variables (SSP1192050)------------------------------------- ### Downloaded future environmental layers from Bio-ORACLE v3.0 (https://www.bio-oracle.org/) ## note: slope remained the same through all scenarios SSP1192050 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP1192050", pattern = "\\.nc$", full.names = TRUE) SSP1192050 <- rast(SSP1192050) ## Crop extent to Aus SSP1192050 <- crop(SSP1192050, extent(Substrate)) ### Mask to 'M' Area setwd("E:/PeerJ/Data 2.0/'M' area/(3) 'M' area (Dissolved)") Marea <- read_sf("'M' area (Dissolved).shp") SSP1192050 <- mask(SSP1192050, Marea) ##check plot(SSP1192050) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP1192050 <- stack(SSP1192050) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP1192050") writeRaster(SSP1192050, filename = "SSP1192050", format = "ascii", bylayer = T, suffix=names(SSP1192050), NAFlag = "-9999", overwrite = T) # Future variables (SSP1192100)------------------------------------- SSP1192100 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP1192100", pattern = "\\.nc$", full.names = TRUE) SSP1192100 <- rast(SSP1192100) ## Crop extent to Aus SSP1192100 <- crop(SSP1192100, extent(Substrate)) ### Mask to 'M' area SSP1192100 <- mask(SSP1192100, Marea) ##check plot(SSP1192100) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP1192100 <- stack(SSP1192100) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP1192100") writeRaster(SSP1192100, filename = "SSP1192100", format = "ascii", bylayer = T, suffix=names(SSP1192100), NAFlag = "-9999", overwrite = T) # Future variables (SSP1262050)------------------------------------- SSP1262050 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP1262050", pattern = "\\.nc$", full.names = TRUE) SSP1262050 <- rast(SSP1262050) ## Crop extent to Aus SSP1262050 <- crop(SSP1262050, extent(Substrate)) ## Mask to 'M' area SSP1262050 <- mask(SSP1262050, Marea) ##check plot(SSP1262050) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP1262050 <- stack(SSP1262050) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP1262050") writeRaster(SSP1262050, filename = "SSP1262050", format = "ascii", bylayer = T, suffix=names(SSP1262050), NAFlag = "-9999", overwrite = T) # Future variables (SSP1262100)------------------------------------- SSP1262100 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP1262100", pattern = "\\.nc$", full.names = TRUE) SSP1262100 <- rast(SSP1262100) ## Crop to extent to Aus SSP1262100 <- crop(SSP1262100, extent(Substrate)) ### Mask to 'M' Area SSP1262100 <- mask(SSP1262100, Marea) #### check plot(SSP1262100) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP1262100 <- stack(SSP1262100) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_13/SSP1262100") writeRaster(SSP1262100, filename = "SSP1262100", format = "ascii", bylayer = T, suffix=names(SSP1262100), NAFlag = "-9999", overwrite = T) # Future variables (SSP2452050)------------------------------------- SSP2452050 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP2452050", pattern = "\\.nc$", full.names = TRUE) SSP2452050 <- rast(SSP2452050) ## Crop to extent to Aus SSP2452050 <- crop(SSP2452050, extent(Substrate)) ### Mask to 'M' Area SSP2452050 <- mask(SSP2452050, Marea) #### check plot(SSP2452050) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP2452050 <- stack(SSP2452050) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP2452050") writeRaster(SSP2452050, filename = "SSP2452050", format = "ascii", bylayer = T, suffix=names(SSP2452050), NAFlag = "-9999", overwrite = T) # Future variables (SSP2452100)------------------------------------- SSP2452100 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP2452100", pattern = "\\.nc$", full.names = TRUE) SSP2452100 <- rast(SSP2452100) ## Crop extent to Aus SSP2452100 <- crop(SSP2452100, extent(Substrate)) ### Mask to 'M' Area SSP2452100 <- mask(SSP2452100, Marea) #### check plot(SSP2452100) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP2452100 <- stack(SSP2452100) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP2452100") writeRaster(SSP2452100, filename = "SSP2452100", format = "ascii", bylayer = T, suffix=names(SSP2452100), NAFlag = "-9999", overwrite = T) # Future variables (SSP3702050)------------------------------------- SSP3702050 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP3702050", pattern = "\\.nc$", full.names = TRUE) SSP3702050 <- rast(SSP3702050) ## Crop extent to Aus SSP3702050 <- crop(SSP3702050, extent(Substrate)) ### Mask to 'M' Area SSP3702050 <- mask(SSP3702050, Marea) #### check plot(SSP3702050) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP3702050 <- stack(SSP3702050) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP3702050") writeRaster(SSP3702050, filename = "SSP3702050", format = "ascii", bylayer = T, suffix=names(SSP3702050), NAFlag = "-9999", overwrite = T) # Future variables (SSP3702100)------------------------------------- SSP3702100 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP3702100", pattern = "\\.nc$", full.names = TRUE) SSP3702100 <- rast(SSP3702100) ## Crop extent to Aus SSP3702100 <- crop(SSP3702100, extent(Substrate)) ### Mask to 'M' Area SSP3702100 <- mask(SSP3702100, Marea) #### check plot(SSP3702100) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP3702100 <- stack(SSP3702100) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP3702100") writeRaster(SSP3702100, filename = "SSP3702100", format = "ascii", bylayer = T, suffix=names(SSP3702100), NAFlag = "-9999", overwrite = T) # Future variables (SSP4602050)------------------------------------- SSP4602050 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP4602050", pattern = "\\.nc$", full.names = TRUE) SSP4602050 <- rast(SSP4602050) ## Crop extent to Aus SSP4602050 <- crop(SSP4602050, extent(Substrate)) ### Mask to 'M' Area SSP4602050 <- mask(SSP4602050, Marea) #### check plot(SSP4602050) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP4602050 <- stack(SSP4602050) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP4602050") writeRaster(SSP4602050, filename = "SSP4602050", format = "ascii", bylayer = T, suffix=names(SSP4602050), NAFlag = "-9999", overwrite = T) # Future variables (SSP4602100)------------------------------------- SSP4602100 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP4602100", pattern = "\\.nc$", full.names = TRUE) SSP4602100 <- rast(SSP4602100) ## Crop extent to Aus SSP4602100 <- crop(SSP4602100, extent(Substrate)) ### Mask to 'M' Area SSP4602100 <- mask(SSP4602100, Marea) #### check plot(SSP4602100) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP4602100 <- stack(SSP4602100) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP4602100") writeRaster(SSP4602100, filename = "SSP4602100", format = "ascii", bylayer = T, suffix=names(SSP4602100), NAFlag = "-9999", overwrite = T) # Future variables (SSP5852050)------------------------------------- SSP5852050 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP5852050", pattern = "\\.nc$", full.names = TRUE) SSP5852050 <- rast(SSP5852050) ## Crop extent to Aus SSP5852050 <- crop(SSP5852050, extent(Substrate)) ### Mask to 'M' Area SSP5852050 <- mask(SSP5852050, Marea) #### check plot(SSP5852050) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP5852050 <- stack(SSP5852050) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP5852050") writeRaster(SSP5852050, filename = "SSP5852050", format = "ascii", bylayer = T, suffix=names(SSP5852050), NAFlag = "-9999", overwrite = T) # Future variables (SSP5852100)------------------------------------- SSP5852100 <- list.files(path = "E:/PeerJ (updated)/(2) Environmental Data/Raw Future Environmental Variables/SSP5852100", pattern = "\\.nc$", full.names = TRUE) SSP5852100 <- rast(SSP5852100) ## Crop extent to Aus SSP5852100 <- crop(SSP5852100, extent(Substrate)) ### Mask to 'M' Area SSP5852100 <- mask(SSP5852100, Marea) #### check plot(SSP5852100) ## Changing Aus from 'Spatraster' to 'Rasterstack' in order to be able to convert to .asc file SSP5852100 <- stack(SSP5852100) ## Convert to .asc file setwd("E:/PeerJ (updated)/Model/WF/G_variables/Set_1/SSP5852100") writeRaster(SSP5852100, filename = "SSP5852100", format = "ascii", bylayer = T, suffix=names(SSP5852100), NAFlag = "-9999", overwrite = T) # Model Projections ------------------------------------------------------- setwd("C:/Model/WF") bfmod <- "batch_model" moddir <- "Final_models" ## name of the output directory to be created and in which all model subdirectories will be created gvars <- "G_variables" ## folder with future environmental variables kuenm_mod(occ.joint = oj, M.var.dir = mvars , out.eval = cresdir, batch = bfmod, rep.type = "Bootstrap", jackknife = TRUE, out.dir = moddir , max.memory = 1000, out.format = "cloglog", project = TRUE, G.var.dir = gvars, ext.type = "ext_clam", write.mess = FALSE, write.clamp = FALSE, maxent.path = mxpath , args = NULL, wait = FALSE, run = TRUE) ## we used the default number of model replicates (default = 10) with the Bootstrap replicate type (rep.type = "Bootstrap") ## jackknife was set to TRUE ## project = TRUE, therefore the model will be projected to scenarios in G.var.dir ## for the whitefin swellshark this was 12 future scenarios under 6 SSPs and two time periods ## ext.type = "ext_clam", is the extrapolation type of projections, we chose extrapolation and clamping (ext_clam) ## out.format = "cloglog", the model output format, we chose cloglog # Evaluation of extrapolation risks (MOP) --------------------------------- ### kuenm_mmop calculates mobility-oriented parity (MOP) layers by comparing environmental values ## between the calibration area and multiple areas or scenarios to which SDMs are transferred outmop <- "MOP_results" ## name of the folder to which MOP results will be written kuenm_mmop(G.var.dir = gvars, M.var.dir = mvars , sets.var = "Set_1", out.mop = outmop , percent = 10, comp.each = 2000, parallel = FALSE, is.swd = FALSE) ## precent = 10, percentage of values sampled from the calibration region to calculate the MOP. ## is.swd = FALSE, whether model calibration and final models were produced using SWD format ## sets.var = "Set_1", variables from G.var.dir and M.var.dir that are going to be compared to create the MOP(s) # Centroids (Range Shift) --------------------------------------------------------------- setwd("E:/PeerJ (updated)/(6) Post Analysis/Centroids") ### Reading centroids shapefiles created in QGIS for current and future scenarios current <- st_read("current.shp") f1192050 <- st_read("1192050.shp") f1262050 <- st_read("1262050.shp") f2452050 <- st_read("2452050.shp") f3702050 <- st_read("3702050.shp") f4602050 <- st_read("4602050.shp") f5852050 <- st_read("5852050.shp") f1192100 <- st_read("1192100.shp") f1262100 <- st_read("1262100.shp") f2452100 <- st_read("2452100.shp") f3702100 <- st_read("3702100.shp") f4602100 <- st_read("4602100.shp") f5852100 <- st_read("5852100.shp") ### Calculating the distance (m) from the current centroid to future centroid distance_f1192050 <- st_distance(current, f1192050) distance_f1262050 <- st_distance(current, f1262050) distance_f2452050 <- st_distance(current, f2452050) distance_f3702050 <- st_distance(current, f3702050) distance_f4602050 <- st_distance(current, f4602050) distance_f5852050 <- st_distance(current, f5852050) distance_f1192100 <- st_distance(current, f1192100) distance_f1262100 <- st_distance(current, f1262100) distance_f2452100 <- st_distance(current, f2452100) distance_f3702100 <- st_distance(current, f3702100) distance_f4602100 <- st_distance(current, f4602100) distance_f5852100 <- st_distance(current, f5852100) # Creating a matrix with coordinates current_x <- st_coordinates(current)[, 1] current_y <- st_coordinates(current)[, 2] f1192050_x <- st_coordinates(f1192050)[, 1] f1192050_y <- st_coordinates(f1192050)[, 2] f1262050_x <- st_coordinates(f1262050)[, 1] f1262050_y <- st_coordinates(f1262050)[, 2] f2452050_x <- st_coordinates(f2452050)[, 1] f2452050_y <- st_coordinates(f2452050)[, 2] f3702050_x <- st_coordinates(f3702050)[, 1] f3702050_y <- st_coordinates(f3702050)[, 2] f4602050_x <- st_coordinates(f4602050)[, 1] f4602050_y <- st_coordinates(f4602050)[, 2] f5852050_x <- st_coordinates(f5852050)[, 1] f5852050_y <- st_coordinates(f5852050)[, 2] f1192100_x <- st_coordinates(f1192100)[, 1] f1192100_y <- st_coordinates(f1192100)[, 2] f1262100_x <- st_coordinates(f1262100)[, 1] f1262100_y <- st_coordinates(f1262100)[, 2] f2452100_x <- st_coordinates(f2452100)[, 1] f2452100_y <- st_coordinates(f2452100)[, 2] f3702100_x <- st_coordinates(f3702100)[, 1] f3702100_y <- st_coordinates(f3702100)[, 2] f4602100_x <- st_coordinates(f4602100)[, 1] f4602100_y <- st_coordinates(f4602100)[, 2] f5852100_x <- st_coordinates(f5852100)[, 1] f5852100_y <- st_coordinates(f5852100)[, 2] # Calculating bearing angle bearing_f1192050 <- bearing(cbind(current_x, current_y), cbind(f1192050_x, f1192050_y)) bearing_f1262050 <- bearing(cbind(current_x, current_y), cbind(f1262050_x, f1262050_y)) bearing_f2452050 <- bearing(cbind(current_x, current_y), cbind(f2452050_x, f2452050_y)) bearing_f3702050 <- bearing(cbind(current_x, current_y), cbind(f3702050_x, f3702050_y)) bearing_f4602050 <- bearing(cbind(current_x, current_y), cbind(f4602050_x, f4602050_y)) bearing_f5852050 <- bearing(cbind(current_x, current_y), cbind(f5852050_x, f5852050_y)) bearing_f1192100 <- bearing(cbind(current_x, current_y), cbind(f1192100_x, f1192100_y)) bearing_f1262100 <- bearing(cbind(current_x, current_y), cbind(f1262100_x, f1262100_y)) bearing_f2452100 <- bearing(cbind(current_x, current_y), cbind(f2452100_x, f2452100_y)) bearing_f3702100 <- bearing(cbind(current_x, current_y), cbind(f3702100_x, f3702100_y)) bearing_f4602100 <- bearing(cbind(current_x, current_y), cbind(f4602100_x, f4602100_y)) bearing_f5852100 <- bearing(cbind(current_x, current_y), cbind(f5852100_x, f5852100_y)) # Converting to compass angle f1192050_compass <- 360 - bearing_f1192050 f1262050_compass <- 360 - bearing_f1262050 f2452050_compass <- 360 - bearing_f2452050 f3702050_compass <- 360 - bearing_f3702050 f4602050_compass <- 360 - bearing_f4602050 f5852050_compass <- 360 - bearing_f5852050 f1192100_compass <- 360 - bearing_f1192100 f1262100_compass <- 360 - bearing_f1262100 f2452100_compass <- 360 - bearing_f2452100 f3702100_compass <- 360 - bearing_f3702100 f4602100_compass <- 360 - bearing_f4602100 f5852100_compass <- 360 - bearing_f5852100 # Converting to radians f1192050_radians <- f1192050_compass * pi / 180 f1262050_radians <- f1262050_compass * pi / 180 f2452050_radians <- f2452050_compass * pi / 180 f3702050_radians <- f3702050_compass * pi / 180 f4602050_radians <- f4602050_compass * pi / 180 f5852050_radians <- f5852050_compass * pi / 180 f1192100_radians <- f1192100_compass * pi / 180 f1262100_radians <- f1262100_compass * pi / 180 f2452100_radians <- f2452100_compass * pi / 180 f3702100_radians <- f3702100_compass * pi / 180 f4602100_radians <- f4602100_compass * pi / 180 f5852100_radians <- f5852100_compass * pi / 180 ## Creating range shift plots for 2050 magnitude <- c(396.4, 305.1, 271.6, 246.4, 175.7, 91.5) angle <- c(4.316, 4.330, 4.336, 4.319, 4.370, 4.400) # Define compass directions directionlabels <- c("N", "NE", "E", "SE", "S", "SW", "W", "NW") # Calculate label positions label.pos <- c(pi/2, pi/4, 0, -pi/4, -pi/2, -3*pi/4, -pi, -5*pi/4) # Radians for compass directions # Create the plot with custom radial limits radial.plot(c(0, magnitude), c(0, angle), lwd = 20, line.col = "NA", # Set to NA for complete transparency labels = directionlabels, label.pos = label.pos, radial.lim = c(0, 400), ticksize = 100, show.grid.label = 1, show.radial.grid = TRUE, grid.unit = " km", grid.col = "black", rad.col = "black") # Define colors for the six arrows colors <- c("#b29269", "#6699cc", "#44aa99", "#882255", "#332288", "#cc79a7") # Adjust colors as needed # Add thicker arrows with adjustments for SE and E directions arrows(x0 = 0, y0 = 0, x1 = magnitude * cos(angle + pi/2), y1 = magnitude * sin(angle + pi/2), length = 0.1, angle = 45, col = colors[1:length(magnitude)], lwd = 4) ## Creating range shift plots for 2100 magnitude21 <- c(1086.7, 944.0, 952.3, 792.2, 327.3, 71.9) angle21 <- c(4.112, 4.203, 4.153, 4.234, 4.332, 4.504) # Define compass directions directionlabels <- c("N", "NE", "E", "SE", "S", "SW", "W", "NW") # Calculate label positions label.pos <- c(pi/2, pi/4, 0, -pi/4, -pi/2, -3*pi/4, -pi, -5*pi/4) # Radians for compass directions # Create the plot with custom radial limits radial.plot(c(0, magnitude21), c(0, angle21), lwd = 4, line.col = "NA", # Set to NA for complete transparency labels = directionlabels, label.pos = label.pos, radial.lim = c(0, 1200), # Overall range from 0 to 1200 ticksize = 100, # Ticks every 200 units show.grid.label = 1, show.radial.grid = TRUE, grid.unit = " km", grid.col = "black", rad.col = "black") # Define colors for the four arrows colors <- c("#b29269", "#882255", "#6699cc", "#44aa99", "#332288", "#cc79a7") # Adjust colors as needed # Add thicker arrows with adjustments for SE and E directions arrows(x0 = 0, y0 = 0, x1 = magnitude21 * cos(angle21 + pi/2), y1 = magnitude21 * sin(angle21 + pi/2), length = 0.1, angle = 45, col = colors[1:length(magnitude21)], lwd = 4)