-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpub_demo_svl_2.0.r
245 lines (217 loc) · 11.6 KB
/
pub_demo_svl_2.0.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
# This file includes a demonstration of how to use the method "svl" to identify the source
# vent location of tephra fall deposits based on thickness or maximum clast size
# measurements.
# Author: Qingyuan Yang, Marcus Bursik, and E. Bruce Pitman.
# GPL: License. Use at your own risk.
###
# The work was supported by National Science Foundation Hazard SEES grant number 1521855
# to G. Valentine, M. Bursik, E.B. Pitman and A.K. Patra, and National Science Foundation
# DMS grant number 1621853 to A.K. Patra, M. Bursik, and E.B. Pitman.
# We appreciate your comments, suggestions, and feedback.
# Please feel free to contact us through vhub or email: qyang5@buffalo.edu
########################################## Demonstration ################################
# In this demonstration, we use a thickness dataset of the North Mono Bed 2
# (Sieh and Bursik, 1986) as an example. North Mono Bed 2 was erupted from Upper Dome, which is a
# part of the Mono-Inyo Craters (eastern central California).
# It is noted that the method can also be applied to maximum clast size measurements.
# Set working directory.
setwd("/the/directory/where/you/put/your/data/")
#---
# The input dataset should have three columns:
# The first two columns are the coordinates (x or E and y or N) of sample sites in the UTM system.
# The third column should include the corresponding thickness or maximum clast size
# measurements in millimeters, with a minimum value of 1 mm.
# The input dataset should be a ".csv" file, and separated by ",".
# The head or column names of the three columns should be "x", "y", and "zr".
# In the case of North Mono Bed 2, the x and y coordinates are in UTM Zone 11.
# Read dataset from the working directory.
td = read.table("nmb2.csv", header=TRUE, sep = ",")
#---
# Check data structure (optional command).
str(td)
# Return:
#'data.frame': 114 obs. of 3 variables:
#$ x : num 320364 328807 327410 329283 327802 ...
#$ y : num 4194955 4192135 4193500 4194283 4194929 ...
#$ zr: num 0 0 0 0 0 0 0 0 0 0 ...
# NOTE: if the data include observations of zero thickness; it is necessary to exclude them for the
# method to work (optional command).
#
#
#---
# Exclude zero-thickness observations.
td = subset(td, zr!=0)
#---
# Check data structure in a different way (optional command).
# Command below shows the first six rows of the dataset.
head(td)
# Return:
# x y zr
#12 326341.2 4199522 22
#15 327685.3 4203216 10
#26 320851.0 4197091 38
#27 321276.5 4197466 82
#28 320946.3 4197606 96
#29 321606.7 4197460 78
# Note: the 1st column with # is just sample numbers; these are not
# part of the dataset. The data are not in order (1, 2, 3,...)
# because zero-thickness observations have been excluded.
#---
# Define the initial guess of source vent location.
# These coordinates should be within the area where the vent is likely
# to be.
# We recommend users to try different coordinates within that area to
# run the method,
# and then check if the results converge to the same point.
sv_assumed = cbind(x = 322227, y = 4181034)
# In the example case, assuming that we known the vent is from the Mono-Inyo Craters,
# the values of the initial guess are the coordinates of Obsidian Dome, which is > 14 km
# south of Upper Dome, the true vent location of North Mono Bed 2.
#---
# Source functions of the method "svl".
# Note: the semi-empirical model proposed by Yang and Bursik (2016)
# is used here.
source("/where/you/put/your/source/code/pub_svl_2.0_exponential.r")
#---
# Run the method "svl".
# Brief description of the method:
# Starting with an initial guess, "sv_assumed", on vent location, the
# value is updated in each iteration (loop).
# Start loop: <----------------|
# The method proposes different possible vent locations around "sv_assumed" |
# that are at a certain distance (distance: h = 1000 m in this case) from it, |
# and compares if they are closer or farther to the true vent compared |
# with "sv_assumed". |
# |
# If "sv_assumed" is closer to the true vent: |
# we shrink the search radius (controlled by r * h, 0.7 * 1000 = 700 here) |
# to improve the resolution, |
# and do the same thing with the smaller search radius; --------------->----/
# If one of the nearby points is closer to the true vent: |
# the guess on vent location, "sv_assumed", |
# is updated/replaced (see function "compare" for more details), |
# and with the same search radius, the method does |
# the same thing with the new guess on vent location.------------------>----/
#
# In the actual implementation of the method, the search radius is updated
# in a slightly more efficient manner (see detailed description in our work, whose
# link will be provided soon).
# This does not affect the notation in this file.
#
# As described above, after each iteration,
# the search radius gets either smaller or unchanged.
# When the search radius is sufficiently small, the vent is found.
# End loop.
# The number of iterations is controlled by "runs".
#---
# How to assign values to the input parameters:
# sv = sv_assumed, the initial guess on the source vent location;
# td: the input dataset;
# runs: the number of iterations;
# h: the initial search radius (in meters);
# r: the shrink ratio in the search radius.
# Suggested value: 0.7
# nlps, h_ratio: these are related to estimating the wind direction in each
# iteration. It is sufficient to keep them fixed as 20 and 1, respectively.
# numb: the number of rows of data points within "td" that will be used as input.
# To use the complete dataset, set it to "1:nrow(td)".
#---
#*** Values of "runs" and "h" need to be selected with care. ***
# Here we use the example of the North Mono Bed 2 to show how the values are determined.
#
# Assuming that we know the tephra was erupted from the Mono-Inyo Craters,
# which cover a range of >15 km in the y direction, and about 5 km in the x-direction.
#
# By specifying the initial search radius to be 1000 m (h = 1000), we ensure that it takes at most
# ~15 steps/iterations before the method shrinks the search radius.
# By the time it shrinks the search radius, it can be assumed approximately that the true vent
# is within ~1500 m of the current point.
#
# Consider the worst-case scenario, namely the search radius keeps shrinking. This means that
# the searching resolution is improving in each iteration at the lowest rate.
# In this case, the initial search radius keeps shrinking at a rate of 0.7 (r = 0.7).
# For the resolution to be at ~10 m scale, the method requires 13 more iterations: 1000*(0.7)^13 = 9.61.
# Based on this, it is ensured that runs = 40 (greater than 15+13)
# iterations are sufficient.
#
# Based on the argument above, we recommend users to follow this strategy to assign values of
# "h" and "runs":
# h = (the maximum height or width of the region of interest)/10
# For "runs", if we want to reach the solution at a ~10 m scale, then
# we need to find out the value "x" in the brackets of the function: h*0.7^(x) = ~10.
# And then set "runs = 10 + x + E", where E could be an integer ranging from 5 to 20, a fudge-factor.
# Greater value in E means more additional iterations,
# which is likely to yield a more accurate estimate.
result_linear = gd_simplified(sv = sv_assumed, td, runs = 40, h = 1000, r = 0.7, nlps = 20, h_ratio = 1, numb = 1:nrow(td))
# Check the output.
result_linear
# Return:
# x y output_ang h (Intercept) dist dd ssr rsquare
#322224.4 4197601 172.3507 0.06281024 2.231973 -0.0003055393 -0.0002597007 2.877154 0.9419849
# Interpretation:
# x and y: estimated vent coordinates of the North Mono Bed 2. This location is within the extent of
# its true source, Upper Dome.
# output_ang: estimated wind direction (from north clockwise). In this case, it is the UPWIND direction.
# See note below on how to tell if the estimated wind direction
# is downwind or upwind.
# h: the length of the search radius (m) from the last iteration. If this value is small
# enough (e.g., <10), then the method converged.
# If this value is comparable to the initial search radius (say 490 = 1000*0.7*0.7),
# then users should try to increase the number of iterations,
# namely increase the value "runs" (e.g., set it to be 60),
# and run the method again.
# If the method is run with more iterations,
# but the resultant "h" is still comparable
# to the initial search radius, that means the method fails
# to converge.
# Then users should try to change the initial
# guess on vent location (sv_assumed), and run the method.
# (intercept): fitted coefficients for the semi-empirical model.
# dist: same as above.
# dd: same as above.
# ssr: sum of squared residuals from the fitting given estimated vent location and
# wind direction.
# rsquare: r-squared value for the final estimate given estimated vent location and wind direction.
#
# Given sparse data, it is important to check if the fitted coefficients are physical.
# For the exponential model, this can be done by checking if dist+dd<0. If their sum dist+dd > 0 or dist > 0, this means that
# along the dispersal direction, the thickness or maximum clast size increases with distance, which is generally unphysical.
#------------------------------------------
# Use the semi-empirical (power-law) model proposed by Gonzalez-Mellado and De la Cruz-Reyna (2010).
# Source the functions of the method "svl".
source("/where/you/put/your/source/code/pub_svl_2.0_power-law.r")
#---
# Run the method "svl" with identical inputs.
result_power = gd_simplified(sv = sv_assumed, td, runs = 40, h = 1000, r = 0.7, nlps = 20, h_ratio = 1, numb = 1:nrow(td))
# Check the result.
result_power
# Return:
# x y output_ang h (Intercept) dif log(dist) ssr rsquare
# 322859.7 4198078 337.9795 0.5338782 5.194675 -0.0001700573 -0.4008759 2.084604 0.9579659
# Interpretation:
# See lines 166-190.
# The estimated vent location is within the area of Upper Dome, the true vent of
# this tephra deposit.
# ### the output_ang in this case is the DOWNwind direction.
# See note below on how to tell if the estimated wind direction is downwind or upwind.
# If log(dist) > 0 or (Intercept) < 0, then unphysical prediction occurs.
######
#---
# Note on how to tell if the estimated wind direction (output_ang) is downwind or upwind:
# Due to the design of "svl", the output "output_ang" could be either the upwind or downwind direction.
# This is related to the gradient descent method adopted in "svl".
# To tell if it is upwind or downwind,
# users need to examine the relationship between the fitted coefficients:
# If "svl" is combined with the exponential method of Yang and Bursik (2016):
# if dd < 0, then "output_ang" is the upwind direction;
# if dd > 0, then "output_ang" is the downwind direction.
# If "svl" is combined with the power-law method of Gonzalez-Mellado and De la Cruz-Reyna (2010):
# if dif < 0, then "output_ang" is the downwind direction;
# if dif > 0, then "output_ang" is the upwind direction.
# References:
# Gonzalez-Mellado, A. O., and S. De la Cruz-Reyna. "A simple semi-empirical approach to model thickness of
# ash-deposits for different #eruption scenarios." Natural Hazards and Earth System Sciences 10.11 (2010): 2241.
# Sieh K, Bursik M. Most recent eruption of the Mono Craters, eastern central California. Journal of Geophysical
# Research: Solid Earth, 1986, 91(B12): 12539-12571.
# Yang Q, Bursik M. A new interpolation method to model thickness, isopachs, extent, and volume of tephra fall
# deposits. Bulletin of Volcanology, 2016, 78(10): 68.