filter surface field

.gitignore

@@ -1,4 +1,5 @@
 *.DS_Store
+paropy/scripts/SHT_example.py
 
 # VS Code
 **/.vscode

paropy/plot_utils.py

@@ -7,6 +7,7 @@
 """
 
 import numpy as np
+import scipy.special as sp
 import math
 from scipy.integrate import cumtrapz
 
@@ -25,19 +26,13 @@ def rad_to_deg(phi,theta):
     for val in theta:
         theta_deg[i]=math.degrees(val)-90
         i=i+1    
-    # phi, theta = np.meshgrid(phi_degree, theta_degree) # indexing='xy')
 
     return (phi_deg, theta_deg)
 
-def get_Z_lim(Z):
-    '''Choose Z limit for plot to the nearest sig. fig. modulo 5'''
-    Z_lim = np.max([np.abs(Z.min()),np.abs(Z.max())])
-    index = np.ceil(-np.log10(Z_lim))
-    modulo = Z_lim % (5*10**(-index))
-    if modulo!=Z_lim:
-        Z_lim = Z_lim - (Z_lim % (5*10**(-index)))
-    else:
-        Z_lim = np.round(Z_lim, index)
+def get_Z_lim(Z,dp=1):
+    '''Choose Z limit for plot to dp decimal places'''
+    Z_lim = np.max(np.abs(Z))
+    Z_lim = np.round(Z_lim, dp)
 
     return Z_lim        
 
@@ -122,3 +117,4 @@ def polar_minimum_latitude(theta,Br):
     pm_lat_south = 90 - theta[idx_south]*180/np.pi
 
     return pm_lat_north, pm_lat_south
+

paropy/routines.py

@@ -6,6 +6,7 @@
 @author: wongj
 """
 import numpy as np
+import scipy.special as sp
 import math
 import h5py
 

paropy/scripts/filter_surface_field.py

@@ -25,15 +25,16 @@
 # directory = '/Users/wongj/Desktop/data/{}'.format(run_ID)
 
 timestamp = '16.84707134'
+l_max = 133 # max. spherical harmonic degree from simulation
+l_trunc = 14 # SH degree truncation
 
 fig_aspect = 1  # figure aspect ratio
 n_levels = 61  # no. of contour levels
 
 saveOn = 1  # save figures?
-saveDir = '/home/wongj/Work/figures/surface'  # path to save files
-# saveDir = '/Users/wongj/Documents/isterre/parody/figures/surface'
-#%%----------------------------------------------------------------------------
-# Load data
+saveDir = '/home/wongj/Work/figures/filter_surface_field'  # path to save files
+# saveDir = '/Users/wongj/Documents/isterre/parody/figures/filter_surface_field'
+#%% Load data
 St_file = 'St={}.{}'.format(timestamp, run_ID)
 filename = '{}/{}'.format(directory, St_file)
 
@@ -42,13 +43,21 @@
  nr, ntheta, nphi, azsym, radius, theta, phi, Vt, Vp,
  Br, dtBr) = surfaceload(filename)
 
-#%%----------------------------------------------------------------------------
-# Plot
+#%% SH transform
+m_max = l_trunc
+sh = shtns.sht(l_trunc, m_max)
+nlat, nlon = sh.set_grid(nphi=nphi, nlat = ntheta)
+vr = Br.T.astype('float64') # NOTE: array has to be dtype='float64' and not 'float32'
+
+coeff = sh.analys(vr) # spatial to spectral
+out = sh.synth(coeff) # spectral to spatial
+
+#%% Plot
 w, h = plt.figaspect(fig_aspect)
 fig, ax = plt.subplots(1, 1, figsize=(1.5*w, h),
                        subplot_kw={'projection': ccrs.Mollweide()})
 X, Y = rad_to_deg(phi, theta)
-Z = Br.T
+Z = out
 Z_lim = get_Z_lim(Z)
 levels = np.linspace(-Z_lim, Z_lim, n_levels)
 c = ax.contourf(X, Y, Z, levels, transform=ccrs.PlateCarree(), cmap='PuOr_r',
@@ -61,3 +70,12 @@
 cbar.ax.tick_params(length=6)
 ax.gridlines()
 ax.set_global()
+
+if saveOn == 1:
+    if not os.path.exists('{}/{}'.format(saveDir, run_ID)):
+        os.makedirs('{}/{}'.format(saveDir, run_ID))
+    fig.savefig('{}/{}/{}.png'.format(saveDir, run_ID, timestamp),
+                format='png', dpi=200, bbox_inches='tight')
+    fig.savefig('{}/{}/{}.pdf'.format(saveDir, run_ID, timestamp),
+                format='pdf', dpi=200, bbox_inches='tight')
+    print('Figures saved as {}/{}/{}.*'.format(saveDir, run_ID, timestamp))

paropy/scripts/surface_timeavg.py

@@ -23,7 +23,7 @@
 from paropy.plot_utils import rad_to_deg, get_Z_lim
 from paropy.routines import surface_timeavg
 
-matplotlib.use('Agg')  # backend for no display
+# matplotlib.use('Agg')  # backend for no display
 #%%--------------------------------------------------------------------------%%
 # INPUT PARAMETERS
 #----------------------------------------------------------------------------%%