bring my plot routines to the actual state, using the latest data lo…

…ading routin
FESOM · Oct 6, 2020 · 7a0b034 · 7a0b034
1 parent 6ffab28
commit 7a0b034
Show file tree

Hide file tree

Showing 15 changed files with 2,336 additions and 5,843 deletions.
diff --git a/view_pscholz/main_plot_horiz.ipynb b/view_pscholz/main_plot_horiz.ipynb
diff --git a/view_pscholz/main_plot_horiz_climanom.ipynb b/view_pscholz/main_plot_horiz_climanom.ipynb
diff --git a/view_pscholz/main_plot_horiz_climanom_3x2lev.ipynb b/view_pscholz/main_plot_horiz_climanom_3x2lev.ipynb
diff --git a/view_pscholz/main_plot_index.ipynb b/view_pscholz/main_plot_index.ipynb
diff --git a/view_pscholz/main_plot_line.ipynb b/view_pscholz/main_plot_line.ipynb
diff --git a/view_pscholz/main_plot_xmoc.ipynb b/view_pscholz/main_plot_xmoc.ipynb
diff --git a/view_pscholz/set_inputarray.py b/view_pscholz/set_inputarray.py
@@ -42,6 +42,8 @@ def set_inputarray():
 	inputarray['proj'     ] 	   = 'cyl' # 'ortho', 'cyl', 'npstere'  
 	inputarray['proj_lon' ] 	   = 0 #only for ortho
 	inputarray['proj_lat' ] 	   = 75 #only for ortho
+
+	inputarray['use_cavity' ] 	   = False #only for ortho
 	#___________________________________________________________________________
 	return inputarray
 
diff --git a/view_pscholz/sub_climatology.py b/view_pscholz/sub_climatology.py
@@ -11,6 +11,7 @@
 from mpl_toolkits.basemap import Basemap
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from colormap_c2c import colormap_c2c
+from sub_fesom_plot  import fesom_plot_lmask
 
 #+_____________________________________________________________________________+
 #|                                                                             |
@@ -190,7 +191,7 @@ def clim_vinterp(data_in,levels):
 #___DO VERTICAL INTERPOLATE AVERAGE OVER CERTAIN LAYERS_________________________
 #
 #_______________________________________________________________________________
-def clim_plot_anom(clim,figsize=[],do_subplot=[]):
+def clim_plot_anom(mesh, clim,figsize=[],do_subplot=[]):
     from set_inputarray import inputarray
 
     if len(figsize)==0 : figsize=[12,8]
@@ -278,7 +279,7 @@ def clim_plot_anom(clim,figsize=[],do_subplot=[]):
     #| adjustable colormap                                                     |
     #+_________________________________________________________________________+
     cnumb = 20; # minimum number of colors
-    if len(clim.cnumb)!=0: cnumb=clim.cnumb[0]
+    if not clim.cnumb == False: cnumb=clim.cnumb
     cmax  = np.nanmax(clim.anom[idx_box])
     cmin  = np.nanmin(clim.anom[idx_box])
     cref  = 0
@@ -340,8 +341,10 @@ def clim_plot_anom(clim,figsize=[],do_subplot=[]):
                         dashes=[1,1e-10],
                         fontsize=fsize)
     map.drawmapboundary(fill_color='0.9',linewidth=1.0)
-    map.drawcoastlines(color='k',linewidth=0.5)
-    map.fillcontinents(color='0.6')
+
+    #map.drawcoastlines(color='k',linewidth=0.5)
+    #map.fillcontinents(color='0.6')
+    fesom_plot_lmask(map,mesh,ax,'0.6')
 
     #___________________________________________________________________________
     # arange colorbar position and labels

diff --git a/view_pscholz/sub_fesom_data.py b/view_pscholz/sub_fesom_data.py
diff --git a/view_pscholz/sub_fesom_data_netcdf4.py b/view_pscholz/sub_fesom_data_netcdf4.py
@@ -42,6 +42,7 @@ class fesom_data(object):
     rescale                     = []
     which_mean                  = 'monthly'
     anom                        = False
+    use_cavity                  = False
 
     #___INIT DATA OBJECT________________________________________________________
     #
@@ -54,6 +55,7 @@ def __init__(self,inputarray):
             self.proj_lon    = inputarray['proj_lon']
             self.proj_lat    = inputarray['proj_lat']
             self.which_plot  = inputarray['which_plot']
+            self.use_cavity  = inputarray['use_cavity']
 
 #___LOAD FESOM2.0 DATA AS TIME AVERAGED HORIZONTAL SLICE________________________
 #
@@ -81,10 +83,22 @@ def fesom_load_data_horiz_netcdf4(mesh,data,             \
         return data
     #_______________________________________________________________________________
     # plot triangle area interpolated to node
+    elif data.var=='nodearea':
+        if len(mesh.nodes_2d_area)==0: mesh.fesom_calc_nodearea()
+        data.value 	= mesh.nodes_2d_area
+        data.sname,data.lname, data.unit, data.cmap= 'triarea', 'Area', 'km^2', 'rygbw'
+        return data
+    #_______________________________________________________________________________
+    # plot triangle area 
     elif data.var=='triarea':
         if len(mesh.nodes_2d_area)==0: mesh.fesom_calc_triarea()
-        data.value 	= mesh.nodes_2d_area
-        data.sname,data.lname, data.unit, data.cmap= 'triarea', 'Area', 'km^2', 'cmocean.cm.balance'
+        data.value 	= np.concatenate((mesh.elem0_2d_area,mesh.elem0_2d_area[mesh.pbndtri_2d_i]))
+        data.sname,data.lname, data.unit, data.cmap= 'triarea', 'Area', 'km^2', 'rygbw'
+        return data
+
+    elif data.var=='cavity_depth':
+        data.value 	= -mesh.nodes_2d_cg
+        data.sname, data.lname, data.unit, data.cmap = 'depth', 'Cavity Depth', 'm', 'wbgyr'
         return data
 
     #___________________________________________________________________________
@@ -108,6 +122,11 @@ def fesom_load_data_horiz_netcdf4(mesh,data,             \
     # compute dimension contained in file
     nti,nsi,ndi = do_filedims(fname_data, nyi, do_output)
 
+    # fix ndi issue for blowup file he can determine automatical if its a 2d or 3d 
+    # variable 
+    if which_files=='blowup_oce' and any(x in data.var for x in ['eta_n','d_eta','ssh_rhs','ssh_rhs_old','zbar_n_bot','zbar_e_bot','bottom_node_thickness','bottom_node_thickness','heat_flux','water_flux']):
+        ndi=0
+
     #___________________________________________________________________________
     # get attributes of variable 
     data        = do_fileattr(fname_data, data, var_list)
@@ -216,10 +235,18 @@ def do_multiyear_fname_list(data, which_files, do_output):
                 var_list = ['temp', 'salt',[]]
                 fname_list[0].append(data.path+'/'+do_fname_mask(which_files,var_list[0],data.runid,str(ayi[yi])))
                 fname_list[1].append(data.path+'/'+do_fname_mask(which_files,var_list[1],data.runid,str(ayi[yi])))
-            elif any(x in data.var for x in ['pgf_x','pgf_y','pgf_xy','pgf']):
+            elif any(x in data.var for x in ['pgf_x','pgf_y','pgf_xy']):
                 var_list = ['pgf_x', 'pgf_y',[]]
                 fname_list[0].append(data.path+'/'+do_fname_mask(which_files,var_list[0],data.runid,str(ayi[yi])))
-                fname_list[1].append(data.path+'/'+do_fname_mask(which_files,var_list[1],data.runid,str(ayi[yi])))    
+                fname_list[1].append(data.path+'/'+do_fname_mask(which_files,var_list[1],data.runid,str(ayi[yi])))
+            elif any(x in data.var for x in ['pgfa_x','pgfa_y','pgfa_xy']):
+                var_list = ['pgfa_x', 'pgfa_y',[]]
+                fname_list[0].append(data.path+'/'+do_fname_mask(which_files,var_list[0],data.runid,str(ayi[yi])))
+                fname_list[1].append(data.path+'/'+do_fname_mask(which_files,var_list[1],data.runid,str(ayi[yi])))
+            elif any(x in data.var for x in ['pgfb_x','pgfb_y','pgfb_xy']):
+                var_list = ['pgfb_x', 'pgfb_y',[]]
+                fname_list[0].append(data.path+'/'+do_fname_mask(which_files,var_list[0],data.runid,str(ayi[yi])))
+                fname_list[1].append(data.path+'/'+do_fname_mask(which_files,var_list[1],data.runid,str(ayi[yi])))        
             elif any(x in data.var for x in ['uv','u','v']):
                 var_list = ['u', 'v',[]]
                 fname_list[0].append(data.path+'/'+do_fname_mask(which_files,var_list[0],data.runid,str(ayi[yi])))
@@ -490,13 +517,28 @@ def do_zinterp(mesh, idata, idepth, ndi, nsi, sel_levidx,do_output):
             weight[zi,0], weight[zi,1] = sel_levidx.index(auxidx), sel_levidx.index(auxidx+1)
             if ndi!=mesh.nlev: 
                 depth = mesh.zmid
+                #_______________________________________________________________
                 if   nsi==mesh.n2dn: isvalid = mesh.nodes_2d_izg[:nsi]-1>=auxidx+1
                 elif nsi==mesh.n2de: isvalid = mesh.elem0_2d_iz[:nsi]-1>=auxidx+1
+                #_______________________________________________________________
+                # case of cavity
+                if (mesh.use_cavity):
+                    if   nsi==mesh.n2dn: isvalid = np.logical_and(mesh.nodes_2d_icg[:nsi]<=auxidx,isvalid)
+                    elif nsi==mesh.n2de: isvalid = np.logical_and(mesh.elem0_2d_ic[:nsi] <=auxidx,isvalid)
+                #_______________________________________________________________
                 valid_lay[isvalid,zi] = valid_lay[isvalid,zi] + 1.0
             else             :
                 depth = mesh.zlev
+                #_______________________________________________________________
                 if   nsi==mesh.n2dn: isvalid = mesh.nodes_2d_izg[:nsi]>=auxidx+1
                 elif nsi==mesh.n2de: isvalid = mesh.elem0_2d_iz[:nsi]>=auxidx+1
+                #_______________________________________________________________
+                # case of cavity
+                if (mesh.use_cavity):
+                    if   nsi==mesh.n2dn: isvalid = np.logical_and(mesh.nodes_2d_icg[:nsi]<=auxidx,isvalid)
+                    elif nsi==mesh.n2de: isvalid = np.logical_and(mesh.elem0_2d_ic[:nsi] <=auxidx,isvalid)
+
+                #_______________________________________________________________
                 valid_lay[isvalid,zi] = valid_lay[isvalid,zi] + 1.0
             weight[zi,2] = (idepth[zi]+depth[auxidx])/(depth[auxidx]-depth[auxidx+1])    
             div[isvalid] = div[isvalid]+1.0
@@ -545,7 +587,41 @@ def do_zinterp(mesh, idata, idepth, ndi, nsi, sel_levidx,do_output):
                 idata[it,isvalid]                 = idata[it,isvalid]/div[isvalid]
                 # set nodes where are no valid layers for interpolation to nan 
                 idata[it,np.logical_not(isvalid)] = np.nan
-
+
+    # if data.depth is empty --> still fill the bottom topography with nan's
+    else:
+        if idata.ndim==3:
+            valid_lay= np.zeros((nsi,len(idepth))) 
+
+            ndimax = mesh.nodes_2d_izg.max()-1
+            if ndi==mesh.nlev: ndimax = mesh.nodes_2d_izg.max()
+
+            for zi in range(0,ndi):
+                if ndi!=mesh.nlev: 
+                    #_______________________________________________________________
+                    if   nsi==mesh.n2dn: isvalid = mesh.nodes_2d_izg[:nsi]-1>=zi+1
+                    elif nsi==mesh.n2de: isvalid = mesh.elem0_2d_iz[:nsi]-1>=zi+1
+                    #_______________________________________________________________
+                    # case of cavity
+                    if (mesh.use_cavity):
+                        if   nsi==mesh.n2dn: isvalid = np.logical_and(mesh.nodes_2d_icg[:nsi]<=zi,isvalid)
+                        elif nsi==mesh.n2de: isvalid = np.logical_and(mesh.elem0_2d_ic[:nsi] <=zi,isvalid)
+                    #_______________________________________________________________
+                else             :
+                    depth = mesh.zlev
+                    #_______________________________________________________________
+                    if   nsi==mesh.n2dn: isvalid = mesh.nodes_2d_izg[:nsi]>=zi+1
+                    elif nsi==mesh.n2de: isvalid = mesh.elem0_2d_iz[:nsi]>=zi+1
+                    #_______________________________________________________________
+                    # case of cavity
+                    if (mesh.use_cavity):
+                        if   nsi==mesh.n2dn: isvalid = np.logical_and(mesh.nodes_2d_icg[:nsi]<=zi,isvalid)
+                        elif nsi==mesh.n2de: isvalid = np.logical_and(mesh.elem0_2d_ic[:nsi] <=zi,isvalid)
+
+                    #_______________________________________________________________
+
+                idata[:,np.logical_not(isvalid),zi] = np.nan
+
     return(idata)
 
 #___DO NECCESSARY POSTPROCESSING________________________________________________
@@ -751,10 +827,22 @@ def do_select_timeidx(data ,nti, nyi, nmi, do_loadloop, do_output):
         if   nti==1:
             sel_timeidx = [0]
             if nmi!=12: print(' --> WARNING --> no monthly information in file --> load annual instead')
+            data.time = np.arange(data.year[0],data.year[1]+1,1)
 
         # file contains monthly data
         elif nti==12:                   
             sel_timeidx = [x-1 for x in data.month]
+            if data.month.size==12:
+                aux_time=list()
+                for year in range(box.year[0],box.year[1]+1):
+                    aux_time.extend([(x-1)/12 + year for x in box.month])
+                data.time=np.array(aux_time)
+            else:
+                aux_mon = np.arange(0,len(box.month),1)
+                aux_time=list()
+                for year in range(box.year[0],box.year[1]+1):
+                    aux_time.extend([x/len(box.month) + year for x in aux_mon])
+                data.time=np.array(aux_time)
 
         # file contains 5 daily data    
         elif nti==73:                   
@@ -890,6 +978,7 @@ def fesom_data_copy(data):
     copy.anom                           = data.anom  
     copy.value2                         = data.value2
     copy.which_mean                     = data.which_mean
+    copy.use_cavity                     = data.use_cavity
 
     #___________________________________________________________________________
     return(copy)
diff --git a/view_pscholz/sub_fesom_mesh.py b/view_pscholz/sub_fesom_mesh.py
@@ -54,13 +54,21 @@ class fesom_mesh:
     nodes_2d_zg, nodes_2d_izg = [], []
     nodes_2d_ig               = []
 
+    # cavity info 
+    use_cavity                = False
+    nodes_2d_cg, nodes_2d_icg = [], []
+    nodes_2d_c,  nodes_2d_ic  = [], []
+    elem0_2d_ic               = []
+
+    # rotated coordiantes
     nodes_2d_xr, nodes_2d_yr  = [], []
 
     pbndn_2d_i                = []
 
     #____mesh elem info__________________________
     elem0_2d_i, elem_2d_i     = [], []
     elem0_2d_iz,pbndtri_2d_i  = [], []
+
     abnormtri_2d_i            = []
 
     #____fesom lsmask polygon____________________
@@ -96,6 +104,10 @@ def __init__(self,inputarray):
         # load mesh from file
         self.fesom_load_mesh()
 
+        if (inputarray['use_cavity']==True):
+            self.use_cavity = True
+            self.fesom_load_cavity()
+
         #_______________________________________________________________________
         # rotate mesh from rot to geo coordinates
         if (inputarray['mesh_rotate'     ]==True):
@@ -106,7 +118,9 @@ def __init__(self,inputarray):
             self.nodes_2d_zg  = self.nodes_2d_z
             self.nodes_2d_izg = self.nodes_2d_iz
             self.nodes_2d_ig  = self.nodes_2d_i
-
+            if (inputarray['use_cavity']==True):
+                self.nodes_2d_cg  = self.nodes_2d_c
+                self.nodes_2d_icg = self.nodes_2d_ic
         #_______________________________________________________________________
         # change grid focus from -180...180 to e.g. 0...360
         if (inputarray['mesh_focus']!=0):
@@ -177,6 +191,29 @@ def fesom_load_mesh(self):
         self.elem0_2d_iz= file_content.values.astype('uint16') - 1
         self.elem0_2d_iz= self.elem0_2d_iz.squeeze()
 
+    #+___LOAD FESOM CAVITY INFORMATION FROM FILE_______________________________+
+    #|                                                                         |
+    #+_________________________________________________________________________+
+    def fesom_load_cavity(self):
+
+        print(' --> read cavity files')
+
+        #____load number of levels at each node_________________________________
+        print('     > cavity_nlvls.out') 
+        file_content    = pa.read_csv(self.path+'cavity_nlvls.out', delim_whitespace=True, skiprows=0, \
+                                       names=['numb_of_lev'])
+        self.nodes_2d_ic= file_content.values.astype('uint16') - 1
+        self.nodes_2d_ic= self.nodes_2d_ic.squeeze()
+        self.nodes_2d_c = np.float32(self.zlev[self.nodes_2d_ic])
+
+        #____load number of levels at each elem_________________________________
+        print('     > cavity_elvls.out') 
+        file_content    = pa.read_csv(self.path+'cavity_elvls.out', delim_whitespace=True, skiprows=0, \
+                                       names=['numb_of_lev'])
+        self.elem0_2d_ic= file_content.values.astype('uint16') - 1
+        self.elem0_2d_ic= self.elem0_2d_ic.squeeze()
+
+
     #+___ROTATE GRID ROT-->GEO_________________________________________________+
     #| input : r2g (default) change coordinate from rotated-->geo              |
     #|         g2r           change coordinate from geo-->rotated              |
@@ -282,6 +319,10 @@ def fesom_grid_rot_r2g(self,str_mode='r2g'):
             self.nodes_2d_zg = self.nodes_2d_z
             self.nodes_2d_izg = self.nodes_2d_iz
             self.nodes_2d_ig = self.nodes_2d_i
+            if (self.use_cavity == True):
+                self.nodes_2d_cg = self.nodes_2d_c
+                self.nodes_2d_icg = self.nodes_2d_ic
+
             if (str_mode == 'focus'):
                 self.fesom_remove_pbnd()
         elif (str_mode == 'g2r'):
@@ -345,6 +386,10 @@ def fesom_remove_pbnd(self):
         self.nodes_2d_zg  = np.concatenate((self.nodes_2d_z,self.nodes_2d_z[self.pbndn_2d_i  ]))
         self.nodes_2d_ig  = np.concatenate((self.nodes_2d_i,self.nodes_2d_i[self.pbndn_2d_i  ]))
         self.nodes_2d_izg = np.concatenate((self.nodes_2d_iz,self.nodes_2d_iz[self.pbndn_2d_i]))
+        if (self.use_cavity==True):
+            self.nodes_2d_cg  = np.concatenate((self.nodes_2d_c,self.nodes_2d_c[self.pbndn_2d_i  ]))
+            self.nodes_2d_icg = np.concatenate((self.nodes_2d_ic,self.nodes_2d_ic[self.pbndn_2d_i]))
+
         self.n2dna=self.n2dn+self.pbndn_2d_i.size # new (augmented) n2d
 
         #____loop over all periodic bnd. segments___________________________________