Skip to content

Commit fed61d4

Browse files
nholschuhnholschuh
committed
Added strainrate calculator
1 parent f28760a commit fed61d4

6 files changed

+285
-13
lines changed

‎add_arrow.py

+45
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
from matplotlib import pyplot as plt
2+
import numpy as np
3+
4+
def add_arrow(lines, position=None, direction='forward', size=15, color='black'):
5+
"""
6+
add an arrow to a line.
7+
8+
line: Line2D object or list of lines
9+
position: x-position of the arrow. If None, mean of xdata is taken
10+
direction: 'left' or 'right'
11+
size: size of the arrow in fontsize points
12+
color: if None, line color is taken.
13+
"""
14+
if isinstance(lines,type([])) == 0:
15+
lines = [lines]
16+
17+
for line in lines:
18+
if color is None:
19+
color = line.get_color()
20+
21+
xdata = line.get_xdata()
22+
ydata = line.get_ydata()
23+
24+
if position is None:
25+
position = xdata.mean()
26+
# find closest index
27+
if direction == 'right':
28+
start_ind = np.argmin(np.absolute(xdata - position))
29+
end_ind = start_ind + 1
30+
elif direction == 'forward':
31+
start_ind = 0;
32+
end_ind = start_ind+1
33+
else:
34+
start_ind = np.argmin(np.absolute(xdata - position))
35+
end_ind = start_ind - 1
36+
37+
#print(start_ind,end_ind)
38+
#print(xdata)
39+
#print(ydata)
40+
line.axes.annotate('',
41+
xytext=(xdata[start_ind], ydata[start_ind]),
42+
xy=(xdata[end_ind], ydata[end_ind]),
43+
arrowprops=dict(arrowstyle="->", color=color),
44+
size=size
45+
)

‎calculate_flowlines.py

+2-2
Original file line numberDiff line numberDiff line change
@@ -66,8 +66,8 @@ def calculate_flowlines(input_xr,seed_points,uv_varnames=['u','v'],xy_varnames=[
6666
xs = temp_xs
6767
ys = temp_ys
6868
else:
69-
xs = np.empty([2,1])
70-
ys = np.empty([2,1])
69+
xs = np.empty([0,len(seed_points)])
70+
ys = np.empty([0,len(seed_points)])
7171

7272

7373
#################### Here is the backward calculation

‎find_cresisfiles.py

+6-1
Original file line numberDiff line numberDiff line change
@@ -58,14 +58,19 @@ def find_cresisfiles(y,m=0,d=0,seg=0,frame=0):
5858

5959
processing_types = sorted(glob.glob(root_dir+season['season']+'/*/'))
6060
search_types = ['standard','music','surf','DEM']
61+
6162
dir_names = [[],[],[],[]]
6263
found_files = [[],[],[],[]]
6364

6465
for ind0,ptype in enumerate(search_types):
6566
type_fdrs,type_fdrs_ind = ndh.str_compare(processing_types,ptype)
6667

6768
for ind1,type_fdr in enumerate(type_fdrs):
68-
file_opts = sorted(glob.glob(type_fdr+dayseg_str+'/'+filestr+'.mat'))
69+
if ind0 < 3:
70+
file_opts = sorted(glob.glob(type_fdr+dayseg_str+'/'+filestr+'.mat'))
71+
else:
72+
file_opts = sorted(glob.glob(type_fdr+dayseg_str+'/'+'_'.join(filestr.split('_')[1:])+'_bottom.mat'))
73+
6974
for ind2, file_select in enumerate(file_opts):
7075
found_files[ind0].append(file_select)
7176
temp_dir_name = file_select.split('/')

‎process_Music_pickedpdf.py

+21-10
Original file line numberDiff line numberDiff line change
@@ -17,20 +17,27 @@
1717
import NDH_Tools as ndh
1818

1919

20-
def process_Music_pickedpdf(fn,data_dir,surf_load,music_load,surf_save,only_edgetrims=0):
20+
def process_Music_pickedpdf(fn,data_dir,surf_load,music_load,surf_save,only_edgetrims=0,remove_framedir=1):
2121
"""
2222
% (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
2323
%
24-
% This function does the standard load, transformation, and plotting
25-
% that is common in the CReSIS radar analysis workflow
24+
% This function takes a picked PDF for a MUSIC image and digitizes
25+
% the picked layers and associated edgetrims
2626
%
2727
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2828
% The inputs are:
29-
%
29+
% fn -- The name of the picked pdf file to process
30+
% data_dir -- The directory that contains the season information for the file
31+
% surf_load -- The name (eg. CSARP_surf_ndh) that loads in the original surface files
32+
% music_load -- The name (eg. CSARP_music3D_ndh) that loads in the original music files
33+
% surf_save -- The name of the directory you want to save the new surf files into (CSARP_surf_iPad)
34+
% only_edgetrims -- [0] In case you only want to extract the edge-trim values
3035
%
3136
%%%%%%%%%%%%%%%
3237
% The outputs are:
33-
%
38+
% This will save a new surf file into the "surf_save" directory. Because of some funny ways
39+
% Python saves .mat files, you need to run a complementary matlab function on your picked files afterward
40+
% "XXXXXXXXXXX"
3441
%
3542
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3643
"""
@@ -56,7 +63,7 @@ def process_Music_pickedpdf(fn,data_dir,surf_load,music_load,surf_save,only_edge
5663
surf_fn = data_dir+surf_load+'/'+day_seg+'/'+local_fn+'.mat'
5764
music_fn = data_dir+music_load+'/'+day_seg+'/'+local_fn+'.mat'
5865
surf_data = ndh.loadmat(surf_fn)
59-
times = ndh.loadmat(music_fn,['Time'])['Time']
66+
times = ndh.loadmat(music_fn,['Time'])['Time'][0]
6067
surf_pick = np.ones(surf_data['surf']['y'][3].shape)*np.NaN
6168

6269
######## These are the properties of the original file that need to be provided to the image processing
@@ -79,7 +86,10 @@ def process_Music_pickedpdf(fn,data_dir,surf_load,music_load,surf_save,only_edge
7986
########## The following converts a pdf to multiple images
8087
print('Starting the pdf deconstruction for: '+local_fn)
8188

82-
comb_deconstruct_dir = '/'.join(fn.split('/')[0:-1])
89+
comb_deconstruct_dir = '/'.join(fn.split('/')[0:-1])+'/temp_frame_deconstruction/'
90+
if os.path.isdir(comb_deconstruct_dir) == 0:
91+
os.makedirs(comb_deconstruct_dir)
92+
8393
os_cmd = 'convert -quality 20 -density 144 %s %s/%s' % (fn,comb_deconstruct_dir,'Frame_%03d.png')
8494
os.system(os_cmd)
8595
frame_list = sorted(glob.glob(comb_deconstruct_dir+'/*.png'))
@@ -179,9 +189,10 @@ def process_Music_pickedpdf(fn,data_dir,surf_load,music_load,surf_save,only_edge
179189
##########################################################################################################
180190
# Part 7 #################################################################################################
181191
########## Here we clean up the temporary directory and save the output
182-
os_cmd = 'rm -r %s/Frame*' % (comb_deconstruct_dir)
183-
print('When ready, run: '+os_cmd)
184-
print(' ')
192+
if remove_framedir != 1:
193+
os_cmd = 'rm -r %s' % (comb_deconstruct_dir[:-1])
194+
print('When ready, run: '+os_cmd)
195+
print(' ')
185196

186197
if len(error_frames) > 0:
187198
print('To test the pixelcoords function, run:')

‎rotate_to_axis.py

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
import numpy as np
2+
3+
def rotate_to_axis(x,y):
4+
"""
5+
% (C) Nick Holschuh - Amherst College - 2022 ([email protected])
6+
% This function rounds a value to a given spacing / order of magnitude
7+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8+
% The inputs are as follows:
9+
%
10+
% data_to_round -- the array that you want to collapse onto the new spacing
11+
% spacing -- the gap between sequential values
12+
% shift -- the offset from 0 to round to
13+
%
14+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15+
% The outputs are as follows:
16+
%
17+
% new_data -- the values of the original array, rounded to their target
18+
% data_inds -- the number of steps away from the lowest value, using the given spacing
19+
%
20+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21+
"""
22+
23+
24+
xmin = np.min(x)
25+
xmax = np.max(x)
26+
ymin = np.min(y)
27+
ymax = np.max(y)
28+
29+
dx = xmax-xmin
30+
dy = ymax-ymin
31+
theta = np.arctan(dy/dx)
32+
33+
unrotate_mat = np.array([[np.cos(-theta),-np.sin(-theta)],[np.sin(-theta),np.cos(-theta)]])
34+
new_xy = np.matmul(np.stack([x,y]).T,unrotate_mat)
35+
new_x = new_xy[:,0]
36+
new_y = new_xy[:,1]
37+
38+
return new_x,new_y

‎strainrate_calculator.py

+173
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
def strainrate_calculator(x_axis,y_axis,u,v,strain_selections,debug_flag=0):
2+
'''
3+
% (C) Nick Holschuh - Penn State University - 2016 ([email protected])
4+
% This calculates the maximum longitudinal strain ("along flow"), the
5+
% minimum longitudinal strain ("across flow"), and the maximum shear, using
6+
% properties of the eigenvectors and eigenvalues of DV tensor. This
7+
% solution is presented in (Hackl et al. 2009), although it is foundational
8+
% theory in continuum mechanics.
9+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10+
% The inputs are as follows:
11+
%
12+
% x_axis - The xaxis values that are associated with the velocity grids
13+
% y_axis - The yaxis values that are associated with the velocity grids
14+
% u - The velocity component in the orientation of the x axis
15+
% v - The velocity component in the orientation of the y axis
16+
% strain_selections - This input allows you to select which of the output
17+
% products you would like to retain. It should be a matrix with 0s and 1s
18+
% in the following positions, to indicate which values you are most
19+
% interested in.
20+
%%%%%%%
21+
% 1 - Max Longitudinal Strain Rate
22+
% 2 - Min Longitudinal Strain Rate
23+
% 3 - Max Shear Strain Rate
24+
% 4 - Max Longitudinal Orientation
25+
% 5 - Min Longitudinal Orientation
26+
% 6 - Rotation Matrix
27+
% 7 - Vertical Strain Rate (assuming incompressibility)
28+
%
29+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30+
% The Output is as follows:
31+
% sr - a Cell array, containing 6 entries.
32+
% 1 - The principle longitudinal strain rate (scaler)
33+
% 2 - The secondary longitudinal strian rate (scaler)
34+
% 3 - The maximum shear strain rate (scaler)
35+
% 4 - The orientation of maximum longitudinal strain (vector)
36+
% 5 - The orientation of minimum longitudinal strain (vector)
37+
% 6 - The rotation matrix (matrix, with positions corresponding to:
38+
% [ 1 2;
39+
% 3 4 ]
40+
% 7 - The vertical strain rate value
41+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42+
%%
43+
'''
44+
45+
import numpy as np
46+
47+
##################################################################
48+
##################################### Variable Initialization
49+
while len(strain_selections) < 7:
50+
strain_selections.append(0)
51+
52+
# Initialize the principle eigenvector matrix (1), if desired
53+
if strain_selections[0] == 1:
54+
e1 = np.zeros(u.shape + (1,))
55+
else:
56+
e1 = 0
57+
58+
# Initialize the principle eigenvector matrix (2), if desired
59+
if strain_selections[1] == 1:
60+
e2 = np.zeros(u.shape + (1,))
61+
else:
62+
e2 = 0
63+
64+
# Initialize the principle eigenvector matrix (3), if desired
65+
if strain_selections[2] == 1:
66+
ss = np.zeros(u.shape + (1,))
67+
else:
68+
ss = 0
69+
70+
# Initialize the principle eigenvector matrix (4), if desired
71+
if strain_selections[3] == 1:
72+
ev1 = np.zeros(u.shape + (2,))
73+
else:
74+
ev1 = 0
75+
76+
# Initialize the secondary eigenvector matrix (5), if desired
77+
if strain_selections[4] == 1:
78+
ev2 = np.zeros(u.shape + (2,))
79+
else:
80+
ev2 = 0
81+
82+
# Initialize the rotation matrix (6), if desired
83+
if strain_selections[5] == 1:
84+
rm = np.zeros(u.shape + (4,))
85+
else:
86+
rm = 0
87+
88+
# Initialize the vertical strain rate matrix (7), if desired
89+
if strain_selections[6] == 1:
90+
vs = np.zeros(u.shape + (1,))
91+
else:
92+
vs = 0
93+
##################################################################
94+
95+
# Generate the change in velocity components
96+
du_y, du_x = np.gradient(u)
97+
dv_y, dv_x = np.gradient(v)
98+
99+
dx = x_axis[1] - x_axis[0]
100+
dy = y_axis[1] - y_axis[0]
101+
102+
dudx = du_x / dx
103+
dudy = du_y / dy
104+
dvdx = dv_x / dx
105+
dvdy = dv_y / dy
106+
107+
# Fill in the components of the strain rate tensor
108+
shear1 = 0.5 * (dudy + dvdx)
109+
rotation1 = 0.5 * (dudy - dvdx)
110+
rotation2 = 0.5 * (dvdx - dudy)
111+
112+
zeromat = np.zeros_like(shear1)
113+
114+
# Initialize arrays to store results
115+
e1 = np.zeros_like(zeromat)
116+
e2 = np.zeros_like(zeromat)
117+
ss = np.zeros_like(zeromat)
118+
ev1 = np.zeros((len(zeromat), len(zeromat[0]), 2))
119+
ev2 = np.zeros((len(zeromat), len(zeromat[0]), 2))
120+
rm = np.zeros((len(zeromat), len(zeromat[0]), 4))
121+
vs = np.zeros_like(zeromat)
122+
123+
for i in range(len(zeromat)):
124+
for j in range(len(zeromat[0])):
125+
# Compute eigenvalues and eigenvectors
126+
calc_mat = [[dudx[i, j], 0.5 * (dudy[i, j] + dvdx[i, j])],
127+
[0.5 * (dudy[i, j] + dvdx[i, j]), dvdy[i, j]]]
128+
129+
#print(calc_mat)
130+
if strain_selections[3] == 1 or strain_selections[4] == 1:
131+
e_val, e_vec = np.linalg.eig(calc_mat)
132+
else:
133+
e_val = np.linalg.eigvals([[dudx[i, j], 0.5 * (dudy[i, j] + dvdx[i, j])],
134+
[0.5 * (dudy[i, j] + dvdx[i, j]), dvdy[i, j]]])
135+
136+
#print(e_val)
137+
#print('-----------------')
138+
139+
# Store results based on strain_selections
140+
if strain_selections[0] == 1:
141+
e1[i, j] = e_val[0]
142+
if strain_selections[1] == 1:
143+
e2[i, j] = e_val[1]
144+
if strain_selections[2] == 1:
145+
ss[i, j] = (np.max(e_val) - np.min(e_val)) / 2
146+
if strain_selections[3] == 1:
147+
ev1[i, j] = e_vec[:, 0]
148+
if strain_selections[4] == 1:
149+
ev2[i, j] = e_vec[:, 1]
150+
if strain_selections[5] == 1:
151+
rm[i, j] = [0, rotation1[i, j], rotation2[i, j], 0]
152+
if strain_selections[6] == 1:
153+
vs[i, j] = -e_val[0] - e_val[1]
154+
155+
# Debugging
156+
if debug_flag == 1:
157+
test_vs = -dudx - dudy
158+
import matplotlib.pyplot as plt
159+
plt.imshow(vs - test_vs)
160+
plt.clim(-0.01, 0.01)
161+
plt.colorbar()
162+
163+
# Store results
164+
sr = [e1, e2, ss, ev1, ev2, rm, vs]
165+
sr_meta = ['Max Longitudinal Strain Rate',
166+
'Min Longitudinal Strain Rate',
167+
'Max Shear Strain Rate',
168+
'Max Longitudinal Orientation',
169+
'Min Longitudinal Orientation',
170+
'Rotation Matrix',
171+
'Vertical Strain Rate']
172+
173+
return sr,sr_meta

0 commit comments

Comments
 (0)