nonergodic-tot_1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b57e5a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#libraries\n",
    "library(pracma)\n",
    "library(ggplot2)\n",
    "library(maps) \n",
    "library(patchwork)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "974501ed",
   "metadata": {},
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in setwd(\"/home/yow004/Documents/GitHub/NonergodicSD\"): cannot change working directory\n",
     "output_type": "error",
     "traceback": [
      "Error in setwd(\"/home/yow004/Documents/GitHub/NonergodicSD\"): cannot change working directory\nTraceback:\n",
      "1. setwd(\"/home/yow004/Documents/GitHub/NonergodicSD\")"
     ]
    }
   ],
   "source": [
    "rm(list = ls())\n",
    "\n",
    "setwd(\"/home/yow004/Documents/GitHub/NonergodicSD\")\n",
    "\n",
    "#user functions\n",
    "source('R_lib/regression/inla/regression_inla_model1_unbounded_hyp.R')\n",
    "\n",
    "\n",
    "#regression name\n",
    "reg_name  <- 'SD_regression_ds1_inla'\n",
    "\n",
    "#input file\n",
    "flatfile_fname <- 'data/compile_all.csv'\n",
    "\n",
    "#output directory\n",
    "dir_out <- 'data/SD_1_ds1_inla/'\n",
    "\n",
    "#INLA regression parameters\n",
    "runinla_flag <- TRUE # TRUE or FALSE\n",
    "res_name     <- 'tot_1'\n",
    "#res_name <- 'tot'\n",
    "#mesh coarseness\n",
    "# #fine\n",
    "# mesh_edge_max     <- 5\n",
    "# mesh_inner_offset <- 15\n",
    "# mesh_outer_offset <- 15\n",
    "# #medium\n",
    "# mesh_edge_max     <- 15\n",
    "# mesh_inner_offset <- 15\n",
    "# mesh_outer_offset <- 50\n",
    "#coarse\n",
    "mesh_edge_max     <- 40\n",
    "mesh_inner_offset <- 40\n",
    "mesh_outer_offset <- 45\n",
    "\n",
    "#approximation options\n",
    "# if flag_gp_approx=TRUE uses int.strategy=\"eb\" and strategy=\"gaussian\"\n",
    "# int.strategy=\"eb\" corresponds to one integration point, and \n",
    "# strategy=\"gaussian\" approximates posteriors as gaussian distributions\n",
    "flag_gp_approx <- TRUE # TRUE or FALSE\n",
    "\n",
    "#number of threads\n",
    "# reduce number of threads if running out of memmory, if not specified\n",
    "# number of CPU threads is used\n",
    "n_threads <- 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee6ae511",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_flatfile <- read.csv(flatfile_fname)\n",
    "df_flatfile[,'tot_trugman']<-df_flatfile[,'PGA']-df_flatfile[,'pPGA']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a1d1576",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "colnames(df_flatfile)[which(names(df_flatfile) == \"eqlat\")] <- \"eqLat\"\n",
    "colnames(df_flatfile)[which(names(df_flatfile) == \"eqlon\")] <- \"eqLon\"\n",
    "colnames(df_flatfile)[which(names(df_flatfile) == \"stalat\")] <- \"staLat\"\n",
    "colnames(df_flatfile)[which(names(df_flatfile) == \"stalon\")] <- \"staLon\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c598612",
   "metadata": {},
   "outputs": [],
   "source": [
    "  # Preprocess Input Data\n",
    "  # ---------------------------\n",
    "  n_data <- nrow(df_flatfile)\n",
    "  #earthquake data\n",
    "  data_eq_all <- df_flatfile[,c('eqid','mag','eqX', 'eqY')]\n",
    "  out_unq  <- UniqueIdxInv(df_flatfile[,'eqid'])\n",
    "  eq_idx   <- out_unq$idx\n",
    "  eq_inv   <- out_unq$inv\n",
    "  data_eq  <- data_eq_all[eq_idx,]\n",
    "  X_eq     <- data_eq[,c(3,4)] #earthquake coordinates\n",
    "  X_eq_all <- data_eq_all[,c(3,4)]\n",
    "  #create earthquake ids for all records (1 to n_eq)\n",
    "  eq_id <- eq_inv\n",
    "  n_eq  <- nrow(data_eq)\n",
    "  \n",
    "  #station data\n",
    "  data_sta_all <- df_flatfile[,c('ssn','Vs30','staX','staY')]\n",
    "  out_unq   <- UniqueIdxInv(df_flatfile[,'ssn'])\n",
    "  sta_idx   <- out_unq$idx\n",
    "  sta_inv   <- out_unq$inv\n",
    "  data_sta  <- data_sta_all[sta_idx,]\n",
    "  X_sta     <- data_sta[,c(3,4)] #station coordinates\n",
    "  X_sta_all <- data_sta_all[,c(3,4)]\n",
    "  #create station indices for all records (1 to n_sta)\n",
    "  sta_id <- sta_inv\n",
    "  n_sta  <- nrow(data_sta)\n",
    "  \n",
    "  #ground-motion observations  \n",
    "  y_data <- df_flatfile[,res_name]\n",
    "  \n",
    "  #UTM zone\n",
    "  utm_zone <- unique(df_flatfile$UTMzone)\n",
    "  utm_no   <- as.numeric(gsub(\"([0-9]+).*$\", \"\\\\1\", utm_zone))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88189cb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_flatfile[,'rsn'] <- 1:n_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6f491d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run INLA, fit model \n",
    "  # ---------------------------\n",
    "  #fixed effects \n",
    "  #---   ---   ---   ---   ---   ---\n",
    "  #prior on the fixed effects\n",
    "  prior_fixed <- list(mean.intercept = 0, prec.intercept = 5,\n",
    "                      mean = (list(intcp=0.0, default=0)),\n",
    "                      prec = (list(intcp=5.0, default=0.01)))\n",
    "  \n",
    "  #covariates\n",
    "  df_inla_covar <- data.frame(intcp = 1, eq = eq_id, sta = sta_id)\n",
    "  \n",
    "  #spatial model\n",
    "  #---   ---   ---   ---   ---   ---\n",
    "  #input arguments\n",
    "  edge_max     <- mesh_edge_max\n",
    "  inner_offset <- mesh_inner_offset\n",
    "  outer_offset <- mesh_outer_offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a111a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#domain mesh\n",
    "mesh <- inla.mesh.2d(loc=rbind(as.matrix(X_eq),as.matrix(X_sta)) ,\n",
    "                    max.edge = c(1,5)*edge_max,\n",
    "                    cutoff = 3, offset = c(inner_offset, outer_offset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1fb4c6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=10, repr.plot.height=10)\n",
    "# Plotting\n",
    "# ---------------------------    \n",
    "#plotting info\n",
    "set1   <- RColorBrewer::brewer.pal(7, \"Set1\") #color map\n",
    "#California\n",
    "map_ca <- subset( map_data(\"state\"), region %in% c(\"california\"))\n",
    "map_ca_utm <- LongLatToUTM(lat=map_ca$lat, lon=map_ca$long, utm_no)\n",
    "map_ca[,c('X','Y')] <- map_ca_utm[,c('X','Y')]/1000\n",
    "#Nevada\n",
    "map_nv <- subset( map_data(\"state\"), region %in% c(\"nevada\"))\n",
    "map_nv_utm <- LongLatToUTM(lat=map_nv$lat, lon=map_nv$long, utm_no)\n",
    "map_nv[,c('X','Y')] <- map_nv_utm[,c('X','Y')]/1000\n",
    "\n",
    "\n",
    "  #Earthquake - Station Mesh\n",
    "  pl_mesh  <- ggplot() + theme_bw() + gg(mesh) +\n",
    "              geom_path(data=map_ca, aes(x=X,y=Y), color='black',size=1) +\n",
    "              geom_point(data=X_eq, aes(x=eqX,y=eqY),size=0.5,color='red') +\n",
    "              geom_point(data=X_sta, aes(x=staX,y=staY),size=2,shape=2,color='dodgerblue') +\n",
    "              scale_size_manual(values=c(2.0,0.5), labels = c('Earthquakes','Stations'), name=element_blank()) + \n",
    "              scale_color_manual(values=c(set1[1],set1[2]), labels = c('Earthquakes','Stations'), name=element_blank()) + \n",
    "              labs(x=\"X (km)\", y=\"Y (km)\") +\n",
    "              xlim(390,800) + ylim(3900,4400) +\n",
    "              theme(plot.title=element_text(size=20), axis.title=element_text(size=20), \n",
    "                    axis.text.y=element_text(size=20), axis.text.x=element_text(size=20),\n",
    "                    legend.key.size = unit(1, 'cm'), legend.text=element_text(size=20), \n",
    "                    legend.position = c(0.20, 0.10))\n",
    "print(pl_mesh)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f542dd59",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1.5\n",
    "out_dir <- dir_out\n",
    "out_fname <- reg_name\n",
    "#prior distributions\n",
    "  #site independent term\n",
    "  prior_omega_1bs <- list(prec = list(prior = \"loggamma\", param = c(0.9, 0.007)))\n",
    "  #spde earthquake prior\n",
    "  spde_eq <- inla.spde2.pcmatern(mesh = mesh, alpha = alpha, # Mesh and smoothness parameter\n",
    "                                 prior.range = c(100, 0.95), # P(range < 100) = 0.95\n",
    "                                 prior.sigma = c(.30, 0.1))  # P(sigma > 0.30) = 0.10\n",
    "  #spde station prior\n",
    "  spde_sta <- inla.spde2.pcmatern(mesh = mesh, alpha = alpha, # Mesh and smoothness parameter\n",
    "                                  prior.range = c(100, 0.95), # P(range < 100) = 0.95\n",
    "                                  prior.sigma = c(.40, 0.1))  # P(sigma > 0.40) = 0.10\n",
    "  \n",
    "  A_eq    <- inla.spde.make.A(mesh, loc = as.matrix(X_eq_all))\n",
    "  idx.eq  <- inla.spde.make.index(\"idx.eq\",spde_eq$n.spde)\n",
    "  A_sta   <- inla.spde.make.A(mesh, loc = as.matrix(X_sta_all))\n",
    "  idx.sta <- inla.spde.make.index(\"idx.sta\",spde_sta$n.spde)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c825b115",
   "metadata": {},
   "outputs": [],
   "source": [
    "#aleatory terms\n",
    "  #---   ---   ---   ---   ---   ---\n",
    "  #prior distributions\n",
    "  prior_phi_0 <- list(prec = list(prior = \"loggamma\", param = c(5.0, 0.5)))\n",
    "  prior_tau_0 <- list(prec = list(prior = \"loggamma\", param = c(4.0, 0.5)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d0a0cf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#inla model\n",
    "  #---   ---   ---   ---   ---   ---\n",
    "  #functional form (with spatial var)\n",
    "  form_inla_spatial <- y ~ 0 + intcp + \n",
    "                           f(eq, model=\"iid\", hyper=prior_tau_0) + f(sta, model=\"iid\", hyper=prior_omega_1bs) +\n",
    "                           f(idx.eq, model = spde_eq) + f(idx.sta, model = spde_sta) \n",
    "  \n",
    "  #build stack\n",
    "  stk_inla_spatial <- inla.stack(data = list(y = y_data),\n",
    "                                 A = list(A_eq, A_sta, 1),\n",
    "                                 effects = list(idx.eq = idx.eq,\n",
    "                                                idx.sta = idx.sta,\n",
    "                                                df_inla_covar),\n",
    "                                 tag = 'model_inla_spatial')\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4312bfe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "  #fit inla model\n",
    "  #---   ---   ---   ---   ---   ---\n",
    "  if(runinla_flag){\n",
    "    #run model (spatial)\n",
    "    if(flag_gp_approx == TRUE){\n",
    "      fit_inla_spatial <- inla(form_inla_spatial, \n",
    "                               data = inla.stack.data(stk_inla_spatial),\n",
    "                               family=\"gaussian\",\n",
    "                               control.family = list(hyper = list(prec = prior_phi_0)),\n",
    "                               control.fixed = prior_fixed,\n",
    "                               control.predictor = list(A = inla.stack.A(stk_inla_spatial)),\n",
    "                               control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),\n",
    "                               control.inla = list(int.strategy='eb', strategy=\"gaussian\"),\n",
    "                               verbose=TRUE,  num.threads=n_threads)\n",
    "    }else{\n",
    "      fit_inla_spatial <- inla(form_inla_spatial, \n",
    "                               data = inla.stack.data(stk_inla_spatial),\n",
    "                               family=\"gaussian\",\n",
    "                               control.family = list(hyper = list(prec = prior_phi_0)),\n",
    "                               control.fixed = prior_fixed,\n",
    "                               control.predictor = list(A = inla.stack.A(stk_inla_spatial)),\n",
    "                               control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),\n",
    "                               verbose=TRUE,  num.threads=n_threads)\n",
    "    }\n",
    "    #save results\n",
    "    dir.create(out_dir, showWarnings=FALSE, recursive=TRUE)\n",
    "    save(fit_inla_spatial, file=file.path(out_dir,paste0(out_fname,'_inla_fit','.Rdata')) )\n",
    "  }else{\n",
    "    #load results\n",
    "    load(file.path(out_dir,paste0(out_fname,'_inla_fit','.Rdata')) )\n",
    "  }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b71c01e",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Post-processing Results\n",
    "  # ---------------------------\n",
    "  #hyper-parameters\n",
    "  hyp_param <- data.frame(matrix(ncol = 6, nrow = 0))\n",
    "  colnames(hyp_param) <- colnames(fit_inla_spatial$summary.hyperpar)\n",
    "\n",
    "  hyp_param['dc_0',]    <- fit_inla_spatial$summary.fixed['intcp',]\n",
    "  #correlation lengths of spatial terms\n",
    "  hyp_param['ell_1e',]  <- fit_inla_spatial$summary.hyperpar['Range for idx.eq',]\n",
    "  hyp_param['ell_1as',] <- fit_inla_spatial$summary.hyperpar['Range for idx.sta',]\n",
    "  #standard deviations of spatial terms\n",
    "  hyp_param['omega_1e',]  <- fit_inla_spatial$summary.hyperpar['Stdev for idx.eq',]\n",
    "  hyp_param['omega_1as',] <- fit_inla_spatial$summary.hyperpar['Stdev for idx.sta',]  \n",
    "  hyp_param['omega_1bs',] <- 1/sqrt(fit_inla_spatial$summary.hyperpar['Precision for sta',] ) \n",
    "  #aleatory terms\n",
    "  hyp_param['phi_0',] <- 1/sqrt( fit_inla_spatial$summary.hyperpar['Precision for the Gaussian observations',] )\n",
    "  hyp_param['tau_0',] <- 1/sqrt( fit_inla_spatial$summary.hyperpar['Precision for eq',] )\n",
    "  #unavailable sd for transformed variables\n",
    "  hyp_param[c('omega_1bs','phi_0','tau_0'),'sd'] <- NA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2495e5ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#projections\n",
    "  prjct_grid_eq  <- inla.mesh.projector(mesh, loc = as.matrix(X_eq))\n",
    "  prjct_grid_sta <- inla.mesh.projector(mesh, loc = as.matrix(X_sta))\n",
    "  \n",
    "  #coefficients    \n",
    "  coeff_1e  <- fit_inla_spatial$summary.random$idx.eq\n",
    "  coeff_1as <- fit_inla_spatial$summary.random$idx.sta\n",
    "  coeff_1bs <- fit_inla_spatial$summary.random$sta\n",
    "  #coeff mean and std\n",
    "  coeff_1e_mu   <- inla.mesh.project(prjct_grid_eq,  coeff_1e$mean)\n",
    "  coeff_1e_sig  <- inla.mesh.project(prjct_grid_eq,  coeff_1e$sd)\n",
    "  coeff_1as_mu  <- inla.mesh.project(prjct_grid_sta, coeff_1as$mean)\n",
    "  coeff_1as_sig <- inla.mesh.project(prjct_grid_sta, coeff_1as$sd)\n",
    "  coeff_1bs_mu  <- coeff_1bs$mean\n",
    "  coeff_1bs_sig <- coeff_1bs$sd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "637eb06f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#mean prediction\n",
    "  y_new_mu <- hyp_param['dc_0','mean'] + coeff_1e_mu[eq_inv] + coeff_1as_mu[sta_inv] + coeff_1bs_mu[sta_inv] \n",
    "  \n",
    "  #residuals\n",
    "  res_tot_mu <- y_data - y_new_mu\n",
    "  res_dB_mu  <- fit_inla_spatial$summary.random$eq$mean[eq_inv]\n",
    "  res_dWS_mu <- res_tot_mu - res_dB_mu\n",
    "  \n",
    "  ## Summarize coefficients and residuals\n",
    "  # ---------------------------\n",
    "  df_flatinfo  <- df_flatfile[,c('rsn','eqid','ssn','eqLat','eqLon','staLat','staLon','eqX','eqY','staX','staY')]\n",
    "  \n",
    "  #summary coefficients\n",
    "  df_coeff <- data.frame(rsn=df_flatinfo$rsn,\n",
    "                         dc_0_mean=hyp_param['dc_0','mean'],\n",
    "                         dc_1e_mean=coeff_1e_mu[eq_inv],  \n",
    "                         dc_1as_mean=coeff_1as_mu[sta_inv],\n",
    "                         dc_1bs_mean=coeff_1bs_mu[sta_inv], \n",
    "                         dc_0_sig=hyp_param['dc_0','sd'], \n",
    "                         dc_1e_sig=coeff_1e_sig[eq_inv], \n",
    "                         dc_1as_sig=coeff_1as_sig[sta_inv], \n",
    "                         dc_1bs_sig=coeff_1bs_sig[sta_inv])\n",
    "  df_coeff <- merge(df_flatinfo, df_coeff, by=c('rsn'))\n",
    "  \n",
    "  #summary predictions and residuals\n",
    "  df_predict_summary <- data.frame(rsn=df_flatinfo$rsn, nerg_mu=y_new_mu, \n",
    "                                   res_tot=res_tot_mu, res_between=res_dB_mu, res_within=res_dWS_mu)\n",
    "  df_predict_summary <- merge(df_flatinfo, df_predict_summary, by=c('rsn'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc6e6a44",
   "metadata": {},
   "outputs": [],
   "source": [
    "  ## Posterior distributions\n",
    "  # ---------------------------\n",
    "  #intercept\n",
    "  post_dc_0 <- as.data.frame(fit_inla_spatial$marginals.fixed$intcp)\n",
    "  #aleatory parameters\n",
    "  post_phi_0 <- as.data.frame(inla.tmarginal(function(x) exp(-x/2), fit_inla_spatial$internal.marginals.hyperpar[['Log precision for the Gaussian observations']]))\n",
    "  post_tau_0 <- as.data.frame(inla.tmarginal(function(x) exp(-x/2), fit_inla_spatial$internal.marginals.hyperpar[['Log precision for eq']]))\n",
    "  #non-ergodic scales\n",
    "  post_omega_1e  <- as.data.frame(inla.tmarginal(function(x) exp( x),   fit_inla_spatial$internal.marginals.hyperpar[['log(Stdev) for idx.eq']]))\n",
    "  post_omega_1as <- as.data.frame(inla.tmarginal(function(x) exp( x),   fit_inla_spatial$internal.marginals.hyperpar[['log(Stdev) for idx.sta']]))\n",
    "  post_omega_1bs <- as.data.frame(inla.tmarginal(function(x) exp(-x/2), fit_inla_spatial$internal.marginals.hyperpar[['Log precision for sta']]))\n",
    "  #correlation length\n",
    "  post_ell_1e   <- as.data.frame(inla.tmarginal(function(x) exp( x), fit_inla_spatial$internal.marginals.hyperpar[['log(Range) for idx.eq']]))\n",
    "  post_ell_1as  <- as.data.frame(inla.tmarginal(function(x) exp( x), fit_inla_spatial$internal.marginals.hyperpar[['log(Range) for idx.sta']]))\n",
    "  \n",
    "  #compute posterior cdfs\n",
    "  post_dc_0$y_int      <- cumtrapz(post_dc_0$x, post_dc_0$y)   / trapz(post_dc_0$x, post_dc_0$y)\n",
    "  post_phi_0$y_int     <- cumtrapz(post_phi_0$x, post_phi_0$y) / trapz(post_phi_0$x, post_phi_0$y)\n",
    "  post_tau_0$y_int     <- cumtrapz(post_tau_0$x, post_tau_0$y) / trapz(post_tau_0$x, post_tau_0$y)\n",
    "  post_omega_1e$y_int  <- cumtrapz(post_omega_1e$x, post_omega_1e$y)   / trapz(post_omega_1e$x, post_omega_1e$y)\n",
    "  post_omega_1as$y_int <- cumtrapz(post_omega_1as$x, post_omega_1as$y) / trapz(post_omega_1as$x, post_omega_1as$y)\n",
    "  post_omega_1bs$y_int <- cumtrapz(post_omega_1bs$x, post_omega_1bs$y) / trapz(post_omega_1bs$x, post_omega_1bs$y)\n",
    "  post_ell_1e$y_int    <- cumtrapz(post_ell_1e$x, post_ell_1e$y)       / trapz(post_ell_1e$x, post_ell_1e$y)\n",
    "  post_ell_1as$y_int   <- cumtrapz(post_ell_1as$x, post_ell_1as$y)     / trapz(post_ell_1as$x, post_ell_1as$y)\n",
    "  \n",
    "  #posterior distributions\n",
    "  #define quantiles\n",
    "  hyp_posterior <- data.frame(quant=seq(0.0,1.0,0.01))\n",
    "  #compute pdf and cdf\n",
    "  if (! all(is.na(post_dc_0$y_int))){\n",
    "    hyp_posterior$dc_0          <- approx(post_dc_0$y_int,      post_dc_0$x,      hyp_posterior$quant)$y\n",
    "    hyp_posterior$dc_0_pdf      <- approx(post_dc_0$y_int,      post_dc_0$y,      hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$dc_0          <- NaN\n",
    "    hyp_posterior$dc_0_pdf      <- NaN\n",
    "  }\n",
    "  if (! all(is.na(post_ell_1e$y_int))){\n",
    "    hyp_posterior$ell_1e        <- approx(post_ell_1e$y_int,    post_ell_1e$x,    hyp_posterior$quant)$y\n",
    "    hyp_posterior$ell_1e_pdf    <- approx(post_ell_1e$y_int,    post_ell_1e$y,    hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$ell_1e        <- NaN\n",
    "    hyp_posterior$ell_1e_pdf    <- NaN\n",
    "  }\n",
    "  if (! all(is.na(post_ell_1as$y_int))){  \n",
    "    hyp_posterior$ell_1as       <- approx(post_ell_1as$y_int,   post_ell_1as$x,   hyp_posterior$quant)$y\n",
    "    hyp_posterior$ell_1as_pdf   <- approx(post_ell_1as$y_int,   post_ell_1as$y,   hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$ell_1as       <- NaN\n",
    "    hyp_posterior$ell_1as_pdf   <- NaN\n",
    "  }\n",
    "  if (! all(is.na(post_omega_1e$y_int))){  \n",
    "    hyp_posterior$omega_1e      <- approx(post_omega_1e$y_int,  post_omega_1e$x,  hyp_posterior$quant)$y\n",
    "    hyp_posterior$omega_1e_pdf  <- approx(post_omega_1e$y_int,  post_omega_1e$y,  hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$omega_1e      <- NaN\n",
    "    hyp_posterior$omega_1e_pdf  <- NaN\n",
    "  }\n",
    "  if (! all(is.na(post_omega_1as$y_int))){  \n",
    "    hyp_posterior$omega_1as     <- approx(post_omega_1as$y_int, post_omega_1as$x, hyp_posterior$quant)$y\n",
    "    hyp_posterior$omega_1as_pdf <- approx(post_omega_1as$y_int, post_omega_1as$y, hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$omega_1as     <- NaN\n",
    "    hyp_posterior$omega_1as_pdf <- NaN\n",
    "  }\n",
    "  if (! all(is.na(post_omega_1bs$y_int))){  \n",
    "    hyp_posterior$omega_1bs     <- approx(post_omega_1bs$y_int, post_omega_1bs$x, hyp_posterior$quant)$y\n",
    "    hyp_posterior$omega_1bs_pdf <- approx(post_omega_1bs$y_int, post_omega_1bs$y, hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$omega_1bs     <- NaN\n",
    "    hyp_posterior$omega_1bs_pdf <- NaN\n",
    "  }\n",
    "  if  (! all(is.na(post_phi_0$y_int))){  \n",
    "    hyp_posterior$phi_0         <- approx(post_phi_0$y_int,     post_phi_0$x,     hyp_posterior$quant)$y\n",
    "    hyp_posterior$phi_0_pdf     <- approx(post_phi_0$y_int,     post_phi_0$y,     hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$phi_0         <- NaN\n",
    "    hyp_posterior$phi_0_pdf     <- NaN\n",
    "  }\n",
    "  if  (! all(is.na(post_tau_0$y_int))){  \n",
    "    hyp_posterior$tau_0         <- approx(post_tau_0$y_int,     post_tau_0$x,     hyp_posterior$quant)$y\n",
    "    hyp_posterior$tau_0_pdf     <- approx(post_tau_0$y_int,     post_tau_0$y,     hyp_posterior$quant)$y\n",
    "  } else {\n",
    "    hyp_posterior$tau_0         <- NaN\n",
    "    hyp_posterior$tau_0_pdf     <- NaN\n",
    "  }\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3f97d62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting\n",
    "  # ---------------------------    \n",
    "  #plotting info\n",
    "  set1   <- RColorBrewer::brewer.pal(7, \"Set1\") #color map\n",
    "  #California\n",
    "  map_ca <- subset( map_data(\"state\"), region %in% c(\"california\"))\n",
    "  map_ca_utm <- LongLatToUTM(lat=map_ca$lat, lon=map_ca$long, utm_no)\n",
    "  map_ca[,c('X','Y')] <- map_ca_utm[,c('X','Y')]/1000\n",
    "  #Nevada\n",
    "  map_nv <- subset( map_data(\"state\"), region %in% c(\"nevada\"))\n",
    "  map_nv_utm <- LongLatToUTM(lat=map_nv$lat, lon=map_nv$long, utm_no)\n",
    "  map_nv[,c('X','Y')] <- map_nv_utm[,c('X','Y')]/1000\n",
    "  \n",
    "  #Earthquake - Station Mesh\n",
    "  pl_mesh  <- ggplot() + theme_bw() + gg(mesh) +\n",
    "              geom_path(data=map_ca, aes(x=X,y=Y), color='black') + geom_path(data=map_nv, aes(x=X,y=Y), color='black')+\n",
    "              geom_point(data=X_eq, aes(x=eqX,y=eqY, size=as.factor('EQ'), color=as.factor('EQ'))) +\n",
    "              geom_point(data=X_sta, aes(x=staX,y=staY, size=as.factor('STA'), color=as.factor('STA'))) +\n",
    "              scale_size_manual(values=c(2.0,0.5), labels = c('Earthquakes','Stations'), name=element_blank()) + \n",
    "              scale_color_manual(values=c(set1[1],set1[2]), labels = c('Earthquakes','Stations'), name=element_blank()) + \n",
    "              labs(x=\"X (km)\", y=\"Y (km)\") +\n",
    "              xlim(390,800) + ylim(3900,4400) +\n",
    "              theme(plot.title=element_text(size=20), axis.title=element_text(size=20), \n",
    "                    axis.text.y=element_text(size=20), axis.text.x=element_text(size=20),\n",
    "                    legend.key.size = unit(1, 'cm'), legend.text=element_text(size=20), \n",
    "                    legend.position = c(0.20, 0.10))\n",
    "\n",
    "  # plot of non-ergodic terms mean and sd of spatially varying event terms\n",
    "  #dc_1e map mean\n",
    "  pl_dc_1e_mu_map <- ggplot() + theme_bw() + xlim(390,800) + ylim(3900,4400)\n",
    "  pl_dc_1e_mu_map <- plot_field(coeff_1e$mean, mesh, xrange=c(-200,800), yrange=c(3400,4750), pl=pl_dc_1e_mu_map) + scico::scale_fill_scico(palette = \"vik\")\n",
    "  pl_dc_1e_mu_map <- pl_dc_1e_mu_map + geom_path(data=map_ca, aes(x=X,y=Y), color='black') + geom_path(data=map_nv, aes(x=X,y=Y), color='black') +\n",
    "                     geom_point(data=X_eq, aes(x=eqX,y=eqY), color=I(\"black\"), size=0.4) +\n",
    "                     labs(x=\"X (km)\", y=\"Y (km)\") + theme(axis.title = element_text(size=20), axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "  #dc_1e map sigma\n",
    "  pl_dc_1e_sd_map <- ggplot() + theme_bw()+ xlim(390,800) + ylim(3900,4400)\n",
    "  pl_dc_1e_sd_map <- plot_field(coeff_1e$sd, mesh, xrange=c(-200,800), yrange=c(3400,4750), pl=pl_dc_1e_sd_map)\n",
    "  pl_dc_1e_sd_map <- pl_dc_1e_sd_map + geom_path(data=map_ca, aes(x=X,y=Y), color='black') + geom_path(data=map_nv, aes(x=X,y=Y), color='black') +\n",
    "                     geom_point(data=X_eq, aes(x=eqX,y=eqY), color=I(\"black\"), size=0.4) +\n",
    "                     labs(x=\"X (km)\", y=\"Y (km)\") + theme(axis.title = element_text(size=20), axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "  #dc_1as map mean\n",
    "  pl_dc_1as_mu_map <- ggplot() + theme_bw()+ xlim(390,800) + ylim(3900,4400)\n",
    "  pl_dc_1as_mu_map <- plot_field(coeff_1as$mean, mesh, xrange=c(-200,800), yrange=c(3400,4750), pl=pl_dc_1as_mu_map) + scico::scale_fill_scico(palette = \"vik\")\n",
    "  pl_dc_1as_mu_map <- pl_dc_1as_mu_map + geom_path(data=map_ca, aes(x=X,y=Y), color='black') + geom_path(data=map_nv, aes(x=X,y=Y), color='black') +\n",
    "                      geom_point(data=X_sta, aes(x=staX,y=staY), color=I(\"black\"), size=0.2) +\n",
    "                      labs(x=\"X (km)\", y=\"Y (km)\") + theme(axis.title = element_text(size=20), axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "  #dc_1as map sigma\n",
    "  pl_dc_1as_sd_map <- ggplot() + theme_bw()+ xlim(390,800) + ylim(3900,4400)\n",
    "  pl_dc_1as_sd_map <- plot_field(coeff_1as$sd, mesh, xrange=c(-200,800), yrange=c(3400,4750), pl=pl_dc_1as_sd_map)\n",
    "  pl_dc_1as_sd_map <- pl_dc_1as_sd_map + geom_path(data=map_ca, aes(x=X,y=Y), color='black') + geom_path(data=map_nv, aes(x=X,y=Y), color='black') +\n",
    "                      geom_point(data=X_sta, aes(x=staX,y=staY), color=I(\"black\"), size=0.2) +\n",
    "                      labs(x=\"X (km)\", y=\"Y (km)\") + theme(axis.title = element_text(size=20), axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "  \n",
    "  #posterior distributions\n",
    "  #dc_0\n",
    "  pl_dc_0_post <- ggplot(post_dc_0, aes(x,y)) + theme_bw() + geom_line() + \n",
    "                  geom_vline(xintercept = hyp_param['dc_0','mean'], colour = \"red\") +\n",
    "                  labs(x = 'dc_0', y = 'posterior', title='Posterior dc_0') + \n",
    "                  theme(plot.title=element_text(size=20), axis.title=element_text(size=20), axis.text.y=element_text(size=20), axis.text.x=element_text(size=20))\n",
    "  #omega_1e\n",
    "  pl_omega_1e_post <- ggplot(post_omega_1e, aes(x,y)) + theme_bw() + geom_line() + \n",
    "                      geom_vline(xintercept = hyp_param['omega_1e','mean'], colour = \"red\") +\n",
    "                      labs(x = 'omega_1e', y = 'posterior', title='Posterior omega_1e') + \n",
    "                      theme(plot.title=element_text(size=20), axis.title=element_text(size=20), axis.text.y=element_text(size=20), axis.text.x=element_text(size=20))\n",
    "  #omega_1as\n",
    "  pl_omega_1as_post <- ggplot(post_omega_1as, aes(x,y)) + theme_bw() + geom_line() + \n",
    "                       geom_vline(xintercept = hyp_param['omega_1as','mean'], colour = \"red\") +\n",
    "                       labs(x = 'omega_1as', y = 'posterior', title='Posterior omega_1as') + \n",
    "                       theme(plot.title=element_text(size=20), axis.title=element_text(size=20), axis.text.y=element_text(size=20), axis.text.x=element_text(size=20))\n",
    "  #omega_1bs\n",
    "  pl_omega_1bs_post <- ggplot(post_omega_1bs, aes(x,y)) + theme_bw() + geom_line() + \n",
    "                       geom_vline(xintercept = hyp_param['omega_1bs','mean'], colour = \"red\") +\n",
    "                       labs(x = 'omega_1bs', y = 'posterior', title='Posterior omega_1bs') + \n",
    "                       theme(plot.title=element_text(size=20), axis.title=element_text(size=20), axis.text.y=element_text(size=20), axis.text.x=element_text(size=20))\n",
    "  #ell_1e\n",
    "  pl_ell_1e_post <- ggplot(post_ell_1e, aes(x,y)) + theme_bw() + geom_line() + \n",
    "                    geom_vline(xintercept = hyp_param['ell_1e','mean'], colour = \"red\") +\n",
    "                    labs(x = 'ell_1e', y = 'posterior', title='Posterior ell_1e') + \n",
    "                    theme(plot.title=element_text(size=20), axis.title=element_text(size=20), axis.text.y=element_text(size=20), axis.text.x=element_text(size=20))\n",
    "  #ell_1as\n",
    "  pl_ell_1as_post <- ggplot(post_ell_1as, aes(x,y)) + theme_bw() + geom_line() + \n",
    "                     geom_vline(xintercept = hyp_param['ell_1as','mean'], colour = \"red\") +\n",
    "                     labs(x = 'ell_1as', y = 'posterior', title='Posterior ell_1as') + \n",
    "                     theme(plot.title=element_text(size=20), axis.title=element_text(size=20), axis.text.y=element_text(size=20), axis.text.x=element_text(size=20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e8043d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=24, repr.plot.height=8)\n",
    "pl_mesh + pl_dc_1e_mu_map + pl_dc_1as_mu_map "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62615a83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write Out\n",
    "  # ---------------------------   \n",
    "  fig_dir <- file.path(out_dir ,'figures')\n",
    "  #create output directories\n",
    "  dir.create(out_dir, showWarnings = FALSE)\n",
    "  dir.create(fig_dir, showWarnings = FALSE)\n",
    "  #data files\n",
    "  # ---   ---   ---   ---   ---\n",
    "  write.csv(as.data.frame(t(hyp_param)), file=file.path(out_dir,paste0(out_fname,'_inla_hyperparameters','.csv')) )\n",
    "  write.csv(df_predict_summary,          file=file.path(out_dir,paste0(out_fname,'_inla_residuals',       '.csv')), row.names = FALSE )\n",
    "  write.csv(df_coeff,  \t          file=file.path(out_dir,paste0(out_fname,'_inla_coefficients',    '.csv')), row.names = FALSE )\n",
    "  write.csv(hyp_posterior,  \t          file=file.path(out_dir,paste0(out_fname,'_inla_hyperposterior',  '.csv')), row.names = FALSE )\n",
    "  \n",
    "  #figures\n",
    "  # ---   ---   ---   ---   ---\n",
    "  #mesh\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_mesh','.png')),  plot=pl_mesh,  device='png')\n",
    "  #spatial distribution of coefficients\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_map_dc_1e_mu','.png')),  plot=pl_dc_1e_mu_map,  device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_map_dc_1e_sd','.png')),  plot=pl_dc_1e_sd_map,  device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_map_dc_1as_mu','.png')), plot=pl_dc_1as_mu_map, device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_map_dc_1as_sd','.png')), plot=pl_dc_1as_sd_map, device='png')\n",
    "  #posterior distribution\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_post_dc_0','.png')),      plot=pl_dc_0_post,      device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_post_omega_1e','.png')),  plot=pl_omega_1e_post,  device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_post_omega_1as','.png')), plot=pl_omega_1as_post, device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_post_omega_1bs','.png')), plot=pl_omega_1bs_post, device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_post_ell_1e','.png')),    plot=pl_ell_1e_post,    device='png')\n",
    "  ggsave(file.path(fig_dir,paste0(out_fname,'_post_ell_1as','.png')),   plot=pl_ell_1as_post,   device='png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e940ce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load non-ergodic coefficients and residuals\n",
    "fname_hparam  <- file.path(dir_out, paste0(reg_name,'_inla_hyperparameters.csv'))\n",
    "fname_gmotion <- file.path(dir_out, paste0(reg_name,'_inla_residuals.csv'))\n",
    "fname_coeff   <- file.path(dir_out, paste0(reg_name,'_inla_coefficients.csv'))\n",
    "df_hparam  <- read.csv(fname_hparam)\n",
    "df_gmotion <- read.csv(fname_gmotion)\n",
    "df_coeff   <- read.csv(fname_coeff)\n",
    "\n",
    "#merge gm-flatfile with non-ergodic coeffs and res \n",
    "df_gmotion <- merge(df_flatfile,df_gmotion)\n",
    "df_coeff   <- merge(df_flatfile,df_coeff)\n",
    "\n",
    "#print mean and standard deviations of aleatory terms\n",
    "fprintf('Between-event mean:\\t %.3f\\n', mean(df_gmotion$res_between))\n",
    "fprintf('Within-event mean:\\t %.3f\\n',  mean(df_gmotion$res_within))\n",
    "fprintf('Between-event std:\\t %.3f\\n',  sd(df_gmotion$res_between))\n",
    "fprintf('Within-event std:\\t %.3f\\n',   sd(df_gmotion$res_within))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ed10453",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#between-event residuals \n",
    "pl_dB_mag <- ggplot(df_gmotion, aes(mag, res_between)) + theme_bw() + ylim(-3,3) + \n",
    "              geom_point() + geom_hline(yintercept=0,linetype=\"dashed\") +\n",
    "              labs(x='magnitude', y='between-event residuals') + \n",
    "              theme(axis.title = element_text(size=20), \n",
    "                    axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "\n",
    "#within-event residuals \n",
    "pl_dW_mag <- ggplot(df_gmotion, aes(mag, res_within)) + theme_bw() +ylim(-3,3) + \n",
    "              geom_point() + geom_hline(yintercept=0,linetype=\"dashed\") +\n",
    "              labs(x='magnitude', y='within-event residuals') + \n",
    "              theme(axis.title = element_text(size=20), \n",
    "                    axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "\n",
    "pl_dW_Rrup <- ggplot(df_gmotion, aes(R, res_within)) + theme_bw() + ylim(-3,3) + \n",
    "               geom_point() + geom_hline(yintercept=0,linetype=\"dashed\") +\n",
    "               labs(x='rupture distance', y='within-event residuals') + \n",
    "               theme(axis.title = element_text(size=20), \n",
    "                     axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "\n",
    "pl_dB_mag + pl_dW_mag + pl_dW_Rrup\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfaf91ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#spatially varying earthquake adjustment\n",
    "pl_dc1e_mu_mag <- ggplot(df_coeff, aes(mag, dc_1e_mean)) + theme_bw() +\n",
    "                    geom_point() + geom_hline(yintercept=0,linetype=\"dashed\") +\n",
    "                    labs(x='magnitude', y='mu(delta c_1e)') + \n",
    "                    theme(axis.title = element_text(size=20), \n",
    "                          axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "\n",
    "#spatially varying site adjustment\n",
    "pl_dc1as_mu_vs30 <- ggplot(df_coeff, aes(log(Vs30), dc_1as_mean)) + theme_bw() +\n",
    "                      geom_point() + geom_hline(yintercept=0,linetype=\"dashed\") +\n",
    "                      labs(x='log V_S30', y='mu(delta c_1as)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "\n",
    "#spatially independent site adjustment\n",
    "pl_dc1bs_mu_vs30 <- ggplot(df_coeff, aes(log(Vs30), dc_1bs_mean)) + theme_bw() +\n",
    "                      geom_point() + geom_hline(yintercept=0,linetype=\"dashed\") +\n",
    "                      labs(x='log V_S30', y='mu(delta c_1bs)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20))\n",
    "\n",
    "pl_dc1e_mu_mag+pl_dc1as_mu_vs30+pl_dc1bs_mu_vs30\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67e44c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#UTM zone\n",
    "utm_zone <- unique(df_flatfile$UTMzone)\n",
    "utm_no   <- as.numeric(gsub(\"([0-9]+).*$\", \"\\\\1\", utm_zone))\n",
    "#California\n",
    "map_ca     <- subset( map_data(\"state\"), region %in% c(\"california\"))\n",
    "map_ca_utm <- LongLatToUTM(lat=map_ca$lat, lon=map_ca$long, utm_no)\n",
    "map_ca[,c('X','Y')] <- map_ca_utm[,c('X','Y')]/1000\n",
    "#Nevada\n",
    "map_nv     <- subset( map_data(\"state\"), region %in% c(\"nevada\"))\n",
    "map_nv_utm <- LongLatToUTM(lat=map_nv$lat, lon=map_nv$long, utm_no)\n",
    "map_nv[,c('X','Y')] <- map_nv_utm[,c('X','Y')]/1000\n",
    "\n",
    "#base map\n",
    "pl_map <- ggplot() + theme_bw() + coord_fixed(ratio = 1) + xlim(390,800) + ylim(3900,4400)+ \n",
    "            geom_path(data=map_ca, aes(x=X,y=Y), color='black') + \n",
    "            geom_path(data=map_nv, aes(x=X,y=Y), color='black')\n",
    "\n",
    "#spatially varying earthquake adjustment\n",
    "pl_dc1e_mu_map <- pl_map + geom_point(data=df_coeff, aes(eqX, eqY, color=dc_1e_mean)) +\n",
    "                      scale_color_gradientn(colours = rainbow(10)) +\n",
    "                      labs(x='X (km)', y='Y (km)', color='mu(delta c_1,E)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20),\n",
    "                            legend.text=element_text(size=20), legend.title=element_text(size=20))\n",
    "\n",
    "pl_dc1e_sd_map <- pl_map + geom_point(data=df_coeff, aes(eqX, eqY, color=dc_1e_sig)) +\n",
    "                      scale_color_gradientn(colours = rev(rainbow(10))) +\n",
    "                      labs(x='X (km)', y='Y (km)', color='sd(delta c_1,E)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20),\n",
    "                            legend.text=element_text(size=20), legend.title=element_text(size=20))\n",
    "\n",
    "#spatially varying site adjustment\n",
    "pl_dc1as_mu_map <- pl_map + geom_point(data=df_coeff, aes(staX, staY, color=dc_1as_mean)) +\n",
    "                      scale_color_gradientn(colours = rainbow(10)) +\n",
    "                      labs(x='X (km)', y='Y (km)', color='mu(delta c_1a,S)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20),\n",
    "                            legend.text=element_text(size=20), legend.title=element_text(size=20))\n",
    "\n",
    "pl_dc1as_sd_map <- pl_map + geom_point(data=df_coeff, aes(staX, staY, color=dc_1as_sig)) +\n",
    "                      scale_color_gradientn(colours = rev(rainbow(10))) +\n",
    "                      labs(x='X (km)', y='Y (km)', color='sd(delta c_1a,S)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20),\n",
    "                            legend.text=element_text(size=20), legend.title=element_text(size=20))\n",
    "\n",
    "#spatially independent site adjustment\n",
    "pl_dc1bs_mu_map <- pl_map + geom_point(data=df_coeff, aes(staX, staY, color=dc_1bs_mean)) +\n",
    "                      scale_color_gradientn(colours = rainbow(10)) +\n",
    "                      labs(x='X (km)', y='Y (km)', color='mu(delta c_1b,S)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20),\n",
    "                            legend.text=element_text(size=20), legend.title=element_text(size=20))\n",
    "\n",
    "pl_dc1bs_sd_map <- pl_map + geom_point(data=df_coeff, aes(staX, staY, color=dc_1bs_sig)) +\n",
    "                      scale_color_gradientn(colours = rev(rainbow(10))) +\n",
    "                      labs(x='X (km)', y='Y (km)', color='sd(delta c_1b,S)') + \n",
    "                      theme(axis.title = element_text(size=20), \n",
    "                            axis.text.y = element_text(size=20), axis.text.x = element_text(size=20),\n",
    "                            legend.text=element_text(size=20), legend.title=element_text(size=20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39399cf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "options(repr.plot.width=15, repr.plot.height=15)\n",
    "pl_dc1e_mu_map + pl_dc1e_sd_map + pl_dc1as_mu_map+pl_dc1as_sd_map + pl_dc1bs_mu_map+ pl_dc1bs_sd_map +  plot_layout(ncol = 2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5f38f24",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b73c90d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c12b62ba",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7abf2ec2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8665b77f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}