Skip to content

Commit 96bd321

Browse files
demiangomezdemiangomez
committed
FixPlate with PPP
1 parent 723dc33 commit 96bd321

File tree

3 files changed

+87
-34
lines changed

3 files changed

+87
-34
lines changed

‎com/FixPlate.py‎

Lines changed: 78 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -10,15 +10,15 @@
1010

1111
import argparse
1212
import os
13-
13+
import json
1414
import numpy as np
1515
from pgamit.pyOkada import cosd, sind
1616
from tqdm import tqdm
1717

1818
from pgamit import dbConnection, pyETM
1919
from pgamit.pyLeastSquares import adjust_lsq
2020
from pgamit.Utils import (cart2euler, get_stack_stations, process_stnlist,
21-
stationID)
21+
stationID, file_write)
2222

2323

2424
def build_design(hdata, vdata):
@@ -96,7 +96,9 @@ def main():
9696
metavar='{stack name}',
9797
help='''Stack name to work with. The Euler pole
9898
will be calculated so as to fix the velocities
99-
of the selected sites in this stack.''')
99+
of the selected sites in this stack. To use PPP
100+
solutions, provide any name for this argument and
101+
pass the -ppp switch.''')
100102

101103
parser.add_argument('-include', '--include_stations', nargs='+', type=str,
102104
metavar='{net.stnm}',
@@ -120,6 +122,16 @@ def main():
120122
If not specified, assumed to be the
121123
production directory.''')
122124

125+
parser.add_argument('-ppp', '--ppp_solutions', action='store_true',
126+
default=False,
127+
help='''Use PPP solutions instead of GAMIT. The
128+
input stack name will be ignored.''')
129+
130+
parser.add_argument('-json', '--save_json', action='store_true',
131+
default=False,
132+
help='''Save json files for the plotted ETMs.
133+
Needs -plot to work.''')
134+
123135
parser.add_argument('-save', '--save_stack', type=str,
124136
metavar='{new stack name}',
125137
help='''Save the time series in the plate-fixed
@@ -168,13 +180,24 @@ def main():
168180
station = stn['StationCode']
169181
network = stn['NetworkCode']
170182

171-
rs = cnn.query_float(f'''SELECT etms.*, stations.lat, stations.lon
172-
FROM etms INNER JOIN stations
173-
USING ("NetworkCode", "StationCode")
174-
WHERE ("NetworkCode", "StationCode", "stack",
175-
"object") =
176-
(\'{network}\', \'{station}\', \'{stack}\',
177-
\'polynomial\')''', as_dict=True)
183+
if not args.ppp_solutions:
184+
# use a GAMIT stack
185+
rs = cnn.query_float(f'''SELECT etms.*, stations.lat, stations.lon
186+
FROM etms INNER JOIN stations
187+
USING ("NetworkCode", "StationCode")
188+
WHERE ("NetworkCode", "StationCode", "stack",
189+
"object") =
190+
(\'{network}\', \'{station}\', \'{stack}\',
191+
\'polynomial\')''', as_dict=True)
192+
else:
193+
# use PPP solutions
194+
rs = cnn.query_float(f'''SELECT etms.*, stations.lat, stations.lon
195+
FROM etms INNER JOIN stations
196+
USING ("NetworkCode", "StationCode")
197+
WHERE ("NetworkCode", "StationCode", "soln",
198+
"object") =
199+
(\'{network}\', \'{station}\', \'ppp\',
200+
\'polynomial\')''', as_dict=True)
178201

179202
if len(rs):
180203
params = np.array(rs[0]['params'])
@@ -193,13 +216,24 @@ def main():
193216
station = stn['StationCode']
194217
network = stn['NetworkCode']
195218

196-
rs = cnn.query_float(f'''SELECT etms.*, stations.lat, stations.lon
197-
FROM etms INNER JOIN stations
198-
USING ("NetworkCode", "StationCode")
199-
WHERE ("NetworkCode", "StationCode",
200-
"stack", "object") =
201-
(\'{network}\', \'{station}\', \'{stack}\',
202-
\'polynomial\')''', as_dict=True)
219+
if not args.ppp_solutions:
220+
# use a GAMIT stack
221+
rs = cnn.query_float(f'''SELECT etms.*, stations.lat, stations.lon
222+
FROM etms INNER JOIN stations
223+
USING ("NetworkCode", "StationCode")
224+
WHERE ("NetworkCode", "StationCode",
225+
"stack", "object") =
226+
(\'{network}\', \'{station}\', \'{stack}\',
227+
\'polynomial\')''', as_dict=True)
228+
else:
229+
# use PPP solutions
230+
rs = cnn.query_float(f'''SELECT etms.*, stations.lat, stations.lon
231+
FROM etms INNER JOIN stations
232+
USING ("NetworkCode", "StationCode")
233+
WHERE ("NetworkCode", "StationCode",
234+
"soln", "object") =
235+
(\'{network}\', \'{station}\', \'ppp\',
236+
\'polynomial\')''', as_dict=True)
203237

204238
if len(rs):
205239
params = np.array(rs[0]['params'])
@@ -270,13 +304,11 @@ def main():
270304
tqdm.write(' -- Obs. nok : %i' % index[~index].shape[0])
271305
tqdm.write(' -- wrms : %.3f mm/yr' % (factor[0, 0] * 1000.))
272306
tqdm.write(' ==== EULER POLE SUMMARY ====')
273-
tqdm.write(''' -- XYZ (mas/yr mas/yr mas/yr) : %8.4f \xB1 %.3f %9.4f
274-
\xB1 %.3f %7.4f \xB1 %.3f'''
307+
tqdm.write(''' -- XYZ (mas/yr mas/yr mas/yr) : %8.4f \xB1 %.3f %9.4f \xB1 %.3f %7.4f \xB1 %.3f'''
275308
% (C[0, 0]*k, sigma[0, 0]*k,
276309
C[1, 0]*k, sigma[0, 1]*k,
277310
C[2, 0]*k, sigma[0, 2]*k))
278-
tqdm.write(''' -- llr (deg deg deg/Myr) : %8.4f \xB1 %.3f %9.4f
279-
\xB1 %.3f %7.4f \xB1 %.3f'''
311+
tqdm.write(''' -- llr (deg deg deg/Myr) : %8.4f \xB1 %.3f %9.4f \xB1 %.3f %7.4f \xB1 %.3f'''
280312
% (lat, np.rad2deg(np.sqrt(cov_lla[1, 1])),
281313
lon, np.rad2deg(np.sqrt(cov_lla[2, 2])),
282314
rot, np.rad2deg(np.sqrt(cov_lla[0, 0])) * 1e-9 * 1e6))
@@ -287,15 +319,24 @@ def main():
287319
v = np.zeros((3, 1))
288320
v[0:3 if len(vref) > 0 else 2] = A @ C
289321
model = pyETM.Model(pyETM.Model.VEL, velocity=v, fit=True)
290-
etm = pyETM.GamitETM(cnn, stn['NetworkCode'], stn['StationCode'],
291-
stack_name=stack, models=[model],
292-
plot_remove_jumps=True)
322+
323+
if not args.ppp_solutions:
324+
etm = pyETM.GamitETM(cnn, stn['NetworkCode'], stn['StationCode'],
325+
stack_name=stack, models=[model],
326+
plot_remove_jumps=True)
327+
else:
328+
etm = pyETM.PPPETM(cnn, stn['NetworkCode'], stn['StationCode'],
329+
models=[model], plot_remove_jumps=True)
293330

294331
xfile = os.path.join(args.directory, '%s.%s_%s'
295332
% (etm.NetworkCode,
296333
etm.StationCode, 'plate-fixed'))
297334
etm.plot(xfile + '.png', plot_missing=False)
298335

336+
if args.save_json:
337+
obj = etm.todictionary(time_series=True, model=True)
338+
file_write(xfile + '.json', json.dumps(obj, indent=4, sort_keys=False))
339+
299340
if save_stack:
300341
stations = get_stack_stations(cnn, args.stack_name[0])
301342

@@ -313,12 +354,19 @@ def main():
313354
v[0:3 if len(vref) > 0 else 2] = A @ C
314355
model = pyETM.Model(pyETM.Model.VEL, velocity=v, fit=True)
315356
try:
316-
etm = pyETM.GamitETM(cnn,
317-
stn['NetworkCode'],
318-
stn['StationCode'],
319-
stack_name=stack,
320-
models=[model],
321-
plot_remove_jumps=True)
357+
if not args.ppp_solutions:
358+
etm = pyETM.GamitETM(cnn,
359+
stn['NetworkCode'],
360+
stn['StationCode'],
361+
stack_name=stack,
362+
models=[model],
363+
plot_remove_jumps=True)
364+
else:
365+
etm = pyETM.PPPETM(cnn,
366+
stn['NetworkCode'],
367+
stn['StationCode'],
368+
models=[model],
369+
plot_remove_jumps=True)
322370

323371
tqdm.write(' -- Saving station %s to stack %s (%i/%i)'
324372
% (stationID(stn), save_stack,

‎pgamit/pyETM.py‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -285,9 +285,11 @@ def __str__(self):
285285
class pyETMException_NoDesignMatrix(pyETMException):
286286
pass
287287

288+
288289
class pyETMException_Model(pyETMException):
289290
pass
290291

292+
291293
def distance(lon1, lat1, lon2, lat2):
292294
"""
293295
Calculate the great circle distance between two points
@@ -376,6 +378,7 @@ def __init__(self, cnn, NetworkCode, StationCode):
376378

377379
self.type = 'ppp'
378380
self.stack_name = 'ppp'
381+
self.project = 'from_ppp'
379382

380383
# get the station from the stations table
381384
stn = cnn.query('SELECT * FROM stations WHERE "NetworkCode" = \'%s\' AND "StationCode" = \'%s\''
@@ -3203,7 +3206,7 @@ def push_params(self, cnn, params=None, reset_polynomial=False, reset_periodic=F
32033206
then year is needed). It can be = None or not passed
32043207
{'object': 'periodic', 'frequencies': list[float]}
32053208
frequencies: a list of integers of the frequencies to fit. This value must be
3206-
passed as 1/days. If no periodic terms are requested, then pass [].
3209+
passed as days. If no periodic terms are requested, then pass [].
32073210
{'object': 'jump', 'Year': int, 'DOY': int, 'action': str, 'jump_type': int, 'relaxation': list}
32083211
Year: sets the year part of the discontinuity date. Mandatory.
32093212
DOY: sets the day of year part of the discontinuity date. Mandatory
@@ -3242,9 +3245,11 @@ def push_params(self, cnn, params=None, reset_polynomial=False, reset_periodic=F
32423245
reset_periodic = True
32433246

32443247
# check that frequencies are not negative
3245-
for f in params['frequencies']:
3248+
for i, f in enumerate(params['frequencies']):
32463249
if f <= 0:
32473250
raise pyETMException('Cannot insert negative frequencies for periodic components')
3251+
else:
3252+
params['frequencies'][i] = 1/f
32483253

32493254
# add the fields for station and network
32503255
params['NetworkCode'] = self.NetworkCode

‎pgamit/pyOkada.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -185,12 +185,12 @@ def __init__(self, cnn, lat, lon, sdate, edate):
185185
# seismic score came back > 0, add jump
186186
self.table.append([j['mag'], Date(datetime=j['date']), j['lon'], j['lat'],
187187
etm.CO_SEISMIC_JUMP_DECAY,
188-
j['id'] + ': M%.1f' % j['mag'] + ' ' + j['location']])
188+
j['id'] + ': M%.1f' % j['mag'] + ' ' + j['location'] + ' -> %.0f km' % dist])
189189
elif p_score > 0:
190190
# seismic score came back == 0, but post-seismic score > 0 add jump
191191
self.table.append([j['mag'], Date(datetime=j['date']), j['lon'], j['lat'],
192192
etm.CO_SEISMIC_DECAY,
193-
j['id'] + ': M%.1f' % j['mag'] + ' ' + j['location']])
193+
j['id'] + ': M%.1f' % j['mag'] + ' ' + j['location'] + ' -> %.0f km' % dist])
194194

195195

196196
class Score(object):

0 commit comments

Comments
 (0)