pythonFLAC/34.phase_interpolate.py at main · chingchenn/pythonFLAC · GitHub

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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 13 14:15:55 2022

@author: chingchen
"""

#-------------------------------------------------------------------
import sys, os
import numpy as np

import flac
import flac_interpolate as fi
import flac_gravity3 as fg
import matplotlib
import matplotlib.pyplot as plt

###############################################

#model = sys.argv[1]
# frame = int(sys.argv[2])
model = 'Nazca_aa06'
model = 'Nazca_v2_01'
#model = 'Ref_Cocos'
# frame = 10

def make_phase_interploate(model,frame):
    path='/home/jiching/geoflac/'
    path = '/Users/chingchen/Desktop/model/'


    os.chdir(path+model)

    fl = flac.Flac()
    savepath='/scratch2/jiching/data/'
    savepath = '/Users/chingchen/Desktop/data/'
    figpath='/scratch2/jiching/figure/'
    figpath = '/Users/chingchen/Desktop/figure/'
    def find_trench_index(z):
        '''Returns the i index of trench location.'''
        zz = z[:,0]
        # the highest point defines the forearc
        imax = zz.argmax()
        # the trench is the lowest point west of forearc
        i = zz[:imax].argmin()
        return i


    def interpolate_phase(frame, xtrench):
        # domain bounds in km
        fi.xmin = xtrench + left
        fi.xmax = xtrench + right
        fi.zmin = down
        fi.zmax = up

        # resolution in km
        fi.dx = dx
        fi.dz = dz

        xx, zz, ph = fi.interpolate(frame, 'phase')
        return xx, zz, ph
    # domain bounds
    left = -200
    right = 1000
    up = 10
    down = -200
    dx = 1.5
    dz = 0.6
    x, z = fl.read_mesh(frame)

    itrench = find_trench_index(z)
    xtrench = x[itrench,0]

    # get interpolated phase either from previous run or from original data
    phasefile = 'intp3-phase.%d' % frame
    xx, zz, ph = interpolate_phase(frame, xtrench)
    f = open(phasefile, 'w')
    f.write('%d %d\n' % xx.shape)
    flac.printing(xx, zz, ph, stream=f)
    f.close()

    # get topography and gravity at uniform spacing
    px, topo, topomod, gravity = fg.compute_gravity2(frame)
    # convert to km and mGal before saving
    px *= 1e-3
    topo *= 1e-3
    topomod *= 1e-3
    gravity *= 1e5
    gfile = 'topo-grav.%d' % frame
    f = open(gfile, 'w')
    flac.printing(px, topo, gravity, topomod, stream=f)
    f.close()


    ###############
    model = os.path.split(os.getcwd())[-1]
    phgrd = model+'_phase_%d.grd' % frame

    xmin = xtrench + left
    xmax = xtrench + right
    zmin = down
    zmax = up
    aspect_ratio = float(up - down) / (right - left)
    width = 6.5
    height = width * aspect_ratio
    shiftz = height + 0.3

    # height of gravity plot
    height2 = 1.0

    # gravity grid
    gravgridsize = 50
    gmin = int(gravity.min() / gravgridsize - 1) * gravgridsize
    gmax = int(gravity.max() / gravgridsize + 1) * gravgridsize
    gravann = max(abs(gmin), abs(gmax))
    # topography grid
    topogridsize = 2
    tpmin = int(topo.min() / topogridsize - 1) * topogridsize
    tpmax = int(topo.max() / topogridsize + 1) * topogridsize
    topoann = max(abs(tpmin), abs(tpmax))
    # interval of temperature contours
    cint = 200

    if not os.path.exists(phgrd):
        cmd = '''tail -n +2 %(phasefile)s | gmt xyz2grd -G%(phgrd)s -I%(dx)f/%(dz)f -R%(xmin)f/%(xmax)f/%(zmin)f/%(zmax)f
        mv  %(phgrd)s %(savepath)s
        ''' % locals()
        print(cmd)
        os.system(cmd)


frame_list=[30,60,120,140]
# end = fl.nrec
for frame in range(2,221,2):
#for frame in frame_list:
    make_phase_interploate(model,frame)
    print('---------------end of '+str(frame)+'---------------')