Merge pull request #14 from GeoscienceAustralia/paper_revisions

Paper revisions

misc/monte_carlo_paper_2021/analysis/probabilistic_inundation/plot_depth_1_in_2475.R

@@ -9,7 +9,9 @@ tonga_coast = readOGR('../../elevation/Tonga_coast/Tonga_coast_nearlon180.shp',
     'Tonga_coast_nearlon180')
 
 # Colours and x/y/depth limits
-COLZ = c('white', rev(rainbow(255)[1:200]))
+#COLZ = c('white', rev(rainbow(255)[1:200]))
+library(cptcity)
+COLZ = c('white', cpt('jjg_cbac_seq_cbacYlGnBu09', n=250)[51:250])
 depth_limit = 6
 XLIM = c(184.625, 185.0)
 YLIM = c(-21.29, -21)
@@ -49,7 +51,8 @@ plot_panel<-function(raster_group_matching_strings, raster_group_tiles, xlim=NUL
                 if(is.null(xlim)) xlim=XLIM
                 if(is.null(ylim)) ylim=YLIM
                 plot(r1, asp=1/cos(mean_lat/180*pi), col=COLZ, xlim=xlim, ylim=ylim,
-                     zlim=c(0, depth_limit), maxpixels=Inf, xaxs='i', yaxs='i')
+                     zlim=c(0, depth_limit), maxpixels=Inf, xaxs='i', yaxs='i', cex.axis=1.3, 
+                     cex.lab=1.3, axis.args=list(cex.axis=1.2))
             }else{
                 image(r1, asp=1/cos(mean_lat/180*pi), zlim = c(0, depth_limit), 
                       col=COLZ, add=TRUE, maxpixels=Inf)

misc/monte_carlo_paper_2021/analysis/probabilistic_inundation/plot_depth_vs_exrate_at_parliament.R

@@ -425,13 +425,12 @@ plot_panel<-function(mean_curve_env, plot_title=''){
     plot(peak_depths, mean_curve_env$mean_exrate, 
         t='l', log='y', ylim=c(1e-04, 1e-02), xlim=c(0.1, 10), lwd=3, las=1, 
         xlab="", ylab="",  cex.lab=1.4, cex.axis=1.4, axes=FALSE)
-    axis(side=1, cex.axis=1.5)
-    axis(side=2, at=10**seq(-5, -1), cex.axis=1.5,
-         labels=c('1/100000', '1/10000', '1/1000', '1/100', '1/10'), 
-         las=1)
+    axis(side=1, cex.axis=1.9)
+    axis(side=2, at=10**seq(-5, -1), cex.axis=1.9,
+         labels=c('1/100000', '1/10000', '1/1000', '1/100', '1/10'))
     add_log_axis_ticks(side=2)
-    mtext(side=1, "Tsunami inundation depth (m)", cex=1.7, line=2.5)
-    mtext(side=2, "Mean exceedances per year ", cex=1.7, line=4.5)
+    mtext(side=1, "Tsunami inundation depth (m)", cex=1.9, line=2.5)
+    mtext(side=2, "Mean exceedances per year ", cex=1.9, line=4.5)
     title(main=plot_title, cex.main=2.2)
     points(peak_depths, mean_curve_env$percentile_exrate[3,],  t='l', col='brown', lwd=3)
     points(peak_depths, mean_curve_env$percentile_exrate[2,],  t='l', col='blue', lwd=2)
@@ -443,7 +442,7 @@ plot_panel<-function(mean_curve_env, plot_title=''){
 
     legend('right', c('Mean', 'Median', '16/84 %', '2.5/97.5 %'), 
            col=c('black', 'brown', 'blue', 'green'),
-           lwd=c(3,3,2,2), pch=c(NA, NA, NA, NA), cex=1.8, bty='n')
+           lwd=c(3,3,2,2), pch=c(NA, NA, NA, NA), cex=2.0, bty='n')
 
     }
 
@@ -466,14 +465,13 @@ unsegmented_curve = add_exrate_curves_with_mc_2sd(unsegmented_random_scenarios,
     curve_type='depth', COL='red')
 segmented_curve = add_exrate_curves_with_mc_2sd(segments_union_random_scenarios, 
     curve_type='depth', COL='green')
-axis(side=1, cex.axis=1.5)
-axis(side=2, at=10**seq(-5, -2), cex.axis=1.5,
-     labels=c('1/100000', '1/10000', '1/1000', '1/100'), 
-     las=1)
+axis(side=1, cex.axis=1.8)
+axis(side=2, at=10**seq(-5, -2), cex.axis=1.8,
+     labels=c('1/100000', '1/10000', '1/1000', '1/100'))
 add_log_axis_ticks(side=2)
 abline(v=seq(0,10),  col='darkgrey', lty='dotted')
-mtext('Tsunami inundation depth (m)', side=1, line=2.5, cex=1.7)
-mtext('Mean exceedances per year', side=2, line=4.5, cex=1.7)
+mtext('Tsunami inundation depth (m)', side=1, line=2.5, cex=1.9)
+mtext('Mean exceedances per year', side=2, line=4.5, cex=1.9)
 
 abline(h=1/c(500, 2500, 10000), col='orange')
 

misc/monte_carlo_paper_2021/analysis/probabilistic_inundation/plot_stage_vs_exrate_at_gauge_3458.R

@@ -305,14 +305,14 @@ abline(h=10**seq(-6, 0), col='orange', lty='dotted')
 add_log_axis_ticks(side=2)
 add_log_axis_ticks(side=1)
 legend('bottomleft',  
-       c('High-resolution Monte-Carlo',
+       c('High-res Monte-Carlo',
          'PTHA18: Logic-tree mean',
-         'PTHA18: 95% CI for Monte-Carlo \nexceedance-rate'
+         'PTHA18: 95% CI for Monte-\nCarlo exceedance-rate'
          ),
        col=c('red', 'blue', 'blue'), pch=c(NA, NA, NA), 
        lty=c('solid', 'solid', 'dashed'), lwd = c(3,1,1),
        pt.cex=c(0.3, NA, NA), bg=rgb(1,1,1,alpha=0.0), 
-       box.col=rgb(1,1,1,alpha=0.0), cex=1.2, yjust=1)
+       box.col=rgb(1,1,1,alpha=0.0), cex=1.4, yjust=1)
 title(main = 'A) Monte-Carlo uncertainties', cex.main = 1.8)
 
 # Panel 2

misc/monte_carlo_paper_2021/optimal_sampling/plot_optimal_sampling.R

@@ -7,6 +7,7 @@ sink(file='log_of_plot_optimal_sampling.log')
 
 # Get the PTHA18 access script in its own environment
 ptha18_access_script = '../../../ptha_access/get_PTHA_results.R'
+#ptha18_access_script = '../../../../../AustPTHA/CODE/ptha/ptha_access/get_PTHA_results.R'
 if(!file.exists(ptha18_access_script)) stop('You need to update the path to the get_PTHA_results.R script')
 ptha18 = new.env()
 source(ptha18_access_script, local=ptha18, chdir=TRUE)
@@ -176,7 +177,8 @@ plot_hazard_curve<-function(
                    c('All scenarios in offshore PTHA',
                        paste0('Monte-Carlo estimates (500 only)'),
                        '95% interval (analytical)'), 
-                   lty=c(1, 1, 1), col=c('black', 'grey', 'orange'), lwd=c(2, 2, 2), cex=1.3)
+                   lty=c('solid', 'solid', 'dotdash'), col=c('black', 'grey', 'darkred'), 
+                   lwd=c(2, 2, 2), cex=1.3)
         }
 
     }else{
@@ -189,16 +191,16 @@ plot_hazard_curve<-function(
                      paste0('Monte-Carlo estimates (500 only)'), 
                      'Equivalent Synthetic Catalogue 95% interval'),
                    lty=c('solid', 'solid', 'dashed'), col=c('black', 'grey', 'darkblue'), 
-                   lwd=c(2, 2, 1), cex=1.3, bg=rgb(1,1,1,alpha=0.7), bty='o', box.col=rgb(1,1,1,alpha=0.3))
+                   lwd=c(2, 1, 2), cex=1.3, bg=rgb(1,1,1,alpha=0.7), bty='o', box.col=rgb(1,1,1,alpha=0.3))
         }else{
 
             legend('bottomleft',
                    c('All scenarios in offshore PTHA', 
                      paste0('Monte-Carlo estimates (500 only)'), 
                      '95% interval (analytical)', 
                      'Equivalent Synthetic Catalogue 95% interval'),
-                   lty=c('solid', 'solid', 'solid', 'dashed'), col=c('black', 'grey', 'orange', 'darkblue'), 
-                   lwd=c(2, 2, 2, 1), cex=1.3, bg=rgb(1,1,1,alpha=0.7), bty='o', box.col=rgb(1,1,1,alpha=0.3))
+                   lty=c('solid', 'solid', 'dotdash', 'dashed'), col=c('black', 'grey', 'darkred', 'darkblue'), 
+                   lwd=c(2, 1, 2, 2), cex=1.3, bg=rgb(1,1,1,alpha=0.7), bty='o', box.col=rgb(1,1,1,alpha=0.3))
 
         }
 
@@ -213,8 +215,8 @@ plot_hazard_curve<-function(
 
         # Add lines to the plot -- because it is log-log we should not use zeros -- instead use a very small
         # threshold
-        points(peak_stage_vals, pmax(equivalent_synthetic_lower, 1e-100), t='l', col='darkblue', lty='dashed')
-        points(peak_stage_vals, pmax(equivalent_synthetic_upper, 1e-100), t='l', col='darkblue', lty='dashed')
+        points(peak_stage_vals, pmax(equivalent_synthetic_lower, 1e-100), t='l', col='darkblue', lty='dashed', lwd=2)
+        points(peak_stage_vals, pmax(equivalent_synthetic_upper, 1e-100), t='l', col='darkblue', lty='dashed', lwd=2)
 
     }
 
@@ -246,8 +248,8 @@ plot_hazard_curve<-function(
             lower_CI[i] = tmp[1] + qnorm(0.025)*sqrt(tmp[2])
             upper_CI[i] = tmp[1] + qnorm(0.975)*sqrt(tmp[2])
         }
-        points(peak_stage_vals, lower_CI, t='l', col='orange', lty='solid')
-        points(peak_stage_vals, upper_CI, t='l', col='orange', lty='solid')
+        points(peak_stage_vals, lower_CI, t='l', col='darkred', lty='dotdash', lwd=2)
+        points(peak_stage_vals, upper_CI, t='l', col='darkred', lty='dotdash', lwd=2)
     }
 
     #
@@ -265,24 +267,24 @@ plot_hazard_curve<-function(
     # Add normal distribution
     xs_local = seq(min(exrate_ts_store$mean), max(exrate_ts_store$mean), len=500)
     points(xs_local, dnorm(xs_local, mean=mean_analytical, sd=sqrt(var_analytical)), t='l',
-           col='orange', lwd=2)
+           col='darkred', lty='dotdash', lwd=3)
 
     if(!add_hardcoded_normal_distribution_to_second_panel){
         legend('topright', 
                paste0("Normal distribution (analytical\n",
                       "mean and variance)"),
-               lwd=2, col='orange', lty='solid', pch=NA, cex=1.25, bty='n')
+               lwd=3, col='darkred', lty='dotdash', pch=NA, cex=1.25, bty='n')
     }else{
         # Useful when we want to compare the spread of Monte-Carlo results with other results
         x_vals = seq(hist_xlim[1], hist_xlim[2], len=201)
         y_vals = dnorm(x_vals, mean=0.00122335093286598, sd=sqrt(0.0000000304813319339911))
-        points(x_vals, y_vals, t='l', col='darkred', lwd=2, lty='dashed')
+        points(x_vals, y_vals, t='l', col='skyblue', lwd=3, lty='dotted')
 
         legend('topright', c(
                paste0("Normal distribution (analytical\n",
                       "mean and variance)"),
                paste0("Stratified-sampling (Figure 2)")),
-               lwd=c(2,2), col=c('orange', 'darkred'), lty=c('solid', 'dashed'), 
+               lwd=c(3,3), col=c('darkred', 'skyblue'), lty=c('dotdash', 'dotted'), 
                pch=c(NA,NA), cex=1.25, bty='n')
     }
     dev.off()
@@ -317,6 +319,7 @@ plot_hazard_curve('stratified', fig_title = 'Exceedance_rate_stratified_target_p
 plot_hazard_curve('stratified_importance', fig_title = 'Exceedance_rate_stratified_importance_target_point.png',
     add_hardcoded_normal_distribution_to_second_panel=TRUE)
 
+#stop()
 #
 # Optimal sampling
 #
@@ -437,6 +440,9 @@ plot_sampling_effort<-function(
     MIN_BAR_HT = 0*(const_samples > 0) # To show a bar can try using a small non-zero value
     BAR_H_OFFSET = 1.5
 
+    library(cptcity)
+    COLZ = cpt(pal='cb_seq_YlOrRd_06', n=6)[3:6]
+
     png('Optimal_sampling_effort.png', width=9, height=6, units='in', res=300)
     par(mar=c(4,4.5,2,1))
     par(mfrow=c(2,1))
@@ -449,7 +455,7 @@ plot_sampling_effort<-function(
 
             plot(optimal_samples[[i]]$Mw + (i-BAR_H_OFFSET)/BAR_GROUP,
                  pmax(optimal_samples[[i]]$Nsamples, MIN_BAR_HT),
-                 t='h', lwd=BAR_LWD, lend=1, col=i, ylim=PLOT_YLIM,
+                 t='h', lwd=BAR_LWD, lend=1, col=COLZ[i], ylim=PLOT_YLIM,
                  xlab='', ylab='Optimal # Samples', las=1,
                  main='Stratified sampling', cex.main=1.8, cex.lab=1.5)
 
@@ -460,12 +466,12 @@ plot_sampling_effort<-function(
 
             points(optimal_samples[[i]]$Mw + (i-BAR_H_OFFSET)/BAR_GROUP,
                    pmax(optimal_samples[[i]]$Nsamples, MIN_BAR_HT),
-                   t='h', lwd=BAR_LWD, lend=1, col=i)
+                   t='h', lwd=BAR_LWD, lend=1, col=COLZ[i])
         }
     }
     points(optimal_samples[[1]]$Mw + (0-BAR_H_OFFSET)/BAR_GROUP,
            pmax(const_samples, MIN_BAR_HT),
-           t='h', lwd=BAR_LWD, lend=1, col='purple')
+           t='h', lwd=BAR_LWD, lend=1, col='black')
     grid(col='orange')
 
     # Report the variance reduction that would be obtained by optimising the
@@ -477,7 +483,7 @@ plot_sampling_effort<-function(
            c('Equal in all bins    (VR = 1.00)', 
               paste0('Threshold = ', threshold_stages, 
               ' m   (VR = ', round(legend_VR_local_optimised ,2),')')), 
-           fill=c('purple', 1:length(threshold_stages)))
+           fill=c('black', COLZ))
 
     # Plot with importance sampling
     for(i in 1:length(threshold_stages)){
@@ -488,7 +494,7 @@ plot_sampling_effort<-function(
 
             plot(optimal_samples_IS[[i]]$Mw + (i-2.5)/BAR_GROUP,
                  pmax(optimal_samples_IS[[i]]$Nsamples, MIN_BAR_HT),
-                 t='h', lwd=BAR_LWD, lend=1, col=i, ylim=PLOT_YLIM,
+                 t='h', lwd=BAR_LWD, lend=1, col=COLZ[i], ylim=PLOT_YLIM,
                  xlab='', ylab='Optimal # Samples', las=1,
                  main='Stratified/importance-sampling', cex.main=1.8, cex.lab=1.5)
 
@@ -499,12 +505,12 @@ plot_sampling_effort<-function(
 
             points(optimal_samples_IS[[i]]$Mw + (i-2.5)/BAR_GROUP,
                    pmax(optimal_samples_IS[[i]]$Nsamples, MIN_BAR_HT),
-                   t='h', lwd=BAR_LWD, lend=1, col=i)
+                   t='h', lwd=BAR_LWD, lend=1, col=COLZ[i])
         }
     }
     grid(col='orange')
     points(optimal_samples_IS[[1]]$Mw + (0-2.5)/BAR_GROUP, pmax(const_samples, MIN_BAR_HT),
-           t='h', lwd=BAR_LWD, lend=1, col='purple')
+           t='h', lwd=BAR_LWD, lend=1, col='black')
 
     # Report the variance reduction that would be obtained by optimising the
     # number of samples for the specific threshold (independent of the
@@ -515,7 +521,7 @@ plot_sampling_effort<-function(
            c('Equal in all bins    (VR = 1.00)', 
               paste0('Threshold = ', threshold_stages, 
                      ' m   (VR = ', round(legend_VR_local_optimised,2),')')), 
-           fill=c('purple', 1:length(threshold_stages)))
+           fill=c('black', COLZ))
     dev.off()
 
 
@@ -531,7 +537,7 @@ plot_sampling_effort<-function(
 
     abline(h=const_samples[1], col='purple', lty='dashed')
     points(optimal_samples[[1]]$Mw + BAR_EPS, mean_optimal_IS_including_constant, ylim=c(0, 150),
-        t='h', lwd=BAR_LWD*1.3, lend=1, col='darkblue')
+        t='h', lwd=BAR_LWD*1.3, lend=1, col='skyblue4')
 
     title(main='Selected non-uniform sampling effort and extra variance-reduction', cex.main=1.5)
 
@@ -540,7 +546,7 @@ plot_sampling_effort<-function(
     legend_VR_chosen = unlist(lapply(optimal_samples,
         function(x) get_variances(x, mean_optimal_including_constant)$constant_on_chosen))
 
-    white_t = rgb(1,1,1,alpha=0.3)
+    white_t = rgb(1,1,1,alpha=0.0)
     legend('topleft', 
            paste0('Threshold = ', rep(threshold_stages, 1), 
                   ' (VR = ', format(round(legend_VR_chosen, 2)), ')'),
@@ -550,8 +556,8 @@ plot_sampling_effort<-function(
            paste0('Threshold = ', rep(threshold_stages, 1), 
                   ' (VR = ', format(round(legend_VR_chosen_IS, 2)), ')'),
            title = 'Stratified/importance-sampling', 
-           text.col='darkblue',
-           title.col='darkblue',
+           text.col='skyblue4',
+           title.col='skyblue4',
            box.col=white_t, cex=1.1, bg=white_t)
 
     text(7.5, 53, 'Uniform sampling', col='purple', cex=1.4)
@@ -756,6 +762,7 @@ dev.off()
 #
 
 ptha18_source_rate_env = new.env()
+#source('../../../../../AustPTHA/CODE/ptha/ptha_access/get_detailed_PTHA18_source_zone_info.R',
 source('../../../ptha_access/get_detailed_PTHA18_source_zone_info.R',
        local=ptha18_source_rate_env, chdir=TRUE)
 
@@ -1038,11 +1045,11 @@ for(nm_i in names_plot_order){
 
     YLIM = max(c(max(ds_stratified_importance$y), max(ds_stratified$y)))
     plot(ds_stratified$x, ds_stratified$y, t='l', col='black', ylim=c(0, YLIM),
-        xlab='', ylab='Density', cex.lab=1.2, cex.axis=1.2,
+        xlab='', ylab='Density', cex.lab=1.4, cex.axis=1.4,
         main=all_titles[[nm_i]], cex.main=1.8, lwd=2)
     points(ds_stratified_importance, col='red', t='l', lty='dashed', lwd=2)
     grid(col='lightblue', lty='dotted')
-    mtext(side=1, bquote(paste('Monte-Carlo exceedance-rate @ ', Q^T, '=2 m')), cex=0.8, line=3)
+    mtext(side=1, bquote(paste('Monte-Carlo exceedance-rate (', Q^T, '=2m)')), cex=0.95, line=3)
     text(max(ds_stratified$x)*0.97, YLIM*0.7,
          paste0('VR = ', signif(VR_analytical, 3), '\n',
                 '     ', '(', signif(ideal_VR_analytical[[nm_i]], 3), ')'),