add example script from florian

example/calcul_pressure-lagrangien-2-eulerien-patch-correctif.py

+#!/usr/bin/python
+
+import sys
+import numpy as np
+import netCDF4 as nc
+from numpy import arange, dtype
+import math
+
+#*******************************************************************************
+def main(args):
+
+    """USE
+  {0} mask.nc input.nc output.nc
+
+DESCRIPTION
+  Computes pressure from density and e3t (that changes with SSH).
+
+"""
+    if len(args)<6 or args[1]=='-h' or args[1]=='--help':
+        print( main.func_doc.format(args[0]) )
+        return 0
+
+
+
+#===============================================================================
+# input
+#===============================================================================
+
+    file_pressure = args[1]
+    fichierpression = nc.Dataset(file_pressure)
+
+    pbc_a=fichierpression.variables['Pbc_a'][:,:,:]
+    pbc_G=fichierpression.variables['Pbc_G'][:,:,:]
+    depth=fichierpression.variables['dept'][:]
+    lon=fichierpression.variables['lont'][:]
+    lat=fichierpression.variables['latt'][:]
+
+    file_ssh = args[2]
+    fichierssh = nc.Dataset(file_ssh)
+
+    ssh_a=fichierssh.variables['ssh_a'][:,:]
+    ssh_G=fichierssh.variables['ssh_G'][:,:]
+
+    jj=pbc_a.shape[1]
+    ii=pbc_a.shape[2]
+    kk=pbc_a.shape[0]
+
+    file_mean = args[3]
+    fichierrhop = nc.Dataset(file_mean)
+    rhop_m=fichierrhop.variables['rhop_mean'][0,:,:]
+
+    file_mask = args[5]
+    fichiermask = nc.Dataset(file_mask)
+    gdepw=fichiermask.variables['gdepw'][0,:,:,:]
+    gdept_0_all=fichiermask.variables['gdept_0'][0,74,:,:]
+    gdept_0=fichiermask.variables['gdept_0'][0,0,:,:]
+    tmask = fichiermask.variables['tmask'][0, :, :, :]
+
+
+
+
+#===============================================================================
+# output
+#===============================================================================
+
+    p_file=args[4]
+    print('Creating '+p_file)
+
+    # this will overwrite the file
+    ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC')
+    setattr(ncfile,'history',' '.join(args))
+    ncfile.createDimension('x',ii)
+    ncfile.createDimension('y',jj)
+    ncfile.createDimension('dept',kk)
+
+
+    data = ncfile.createVariable('lont',dtype('float64').char,('y','x'))
+    data[:] = lon
+    setattr(data, 'units',"degrees_north")
+    setattr(data, 'standard_name',"longitude_at_T_location")
+
+    data = ncfile.createVariable('latt',dtype('float64').char,('y','x'))
+    data[:] = lat
+    setattr(data, 'units',"degrees_east")
+    setattr(data, 'standard_name',"latitude_at_T_location")
+
+    data = ncfile.createVariable('dept',dtype('float64').char,('dept'))
+    data[:] = depth
+    setattr(data, 'units',"m")
+    setattr(data, 'positive',"down")
+    setattr(data, 'standard_name',"depth_at_T_location")
+    setattr(data, 'axis',"Z")
+
+    p_corr_a = ncfile.createVariable('P_corr_A',dtype('float32').char,('dept','y','x'))
+    setattr(p_corr_a, 'standard_name',"Amplitude Of Pressure Total")
+    setattr(p_corr_a, 'units',"m-1.kg.s-2")
+    setattr(p_corr_a, 'coordinates', "dept latt lont")
+
+    p_corr_G = ncfile.createVariable('P_corr_G',dtype('float32').char,('dept','y','x'))
+    setattr(p_corr_G, 'standard_name',"Phase Of Pressure Total")
+    setattr(p_corr_G, 'units',"m-1.kg.s-2")
+    setattr(p_corr_G, 'coordinates', "dept latt lont")
+
+#===============================================================================
+# process
+#===============================================================================
+
+    print( 'Dimensions: kk={} jj={} ii={} '.format(kk,jj,ii) )
+
+    p_zr = np.zeros(shape=(kk,jj,ii))
+    p_zi = np.zeros(shape=(kk,jj,ii))
+
+    ssh_zr = np.zeros(shape=(jj,ii))
+    ssh_zi = np.zeros(shape=(jj,ii))
+
+    grav=9.81
+
+
+    for k in range (0,75):
+        for j in range (0,484):
+            for i in range (0,388):
+
+                if math.isnan(pbc_a[k, j, i])==True:
+                    p_zr[k, j, i] == np.NaN;
+                    p_zi[k, j, i] == np.NaN;
+                else:
+                    p_zr[k, j, i]=pbc_a[k, j, i]*math.cos(math.radians(pbc_G[k, j, i]))
+                    p_zi[k, j, i]=pbc_a[k, j, i]*math.sin(math.radians(pbc_G[k, j, i]))
+
+    for j in range (0,484):
+        for i in range (0,388):
+
+            if math.isnan(ssh_a[j, i])==True:
+                ssh_zr[j, i] == np.NaN;
+                ssh_zi[j, i] == np.NaN;
+            else:
+                ssh_zr[j, i]=ssh_a[j, i]*math.cos(math.radians(ssh_G[j, i]))
+                ssh_zi[j, i]=ssh_a[j, i]*math.sin(math.radians(ssh_G[j, i]))
+
+    p_corr_zr = np.zeros(shape=(kk,jj,ii))
+    p_corr_zi = np.zeros(shape=(kk,jj,ii))
+    zz = np.zeros(shape=(jj,ii))
+
+    for k in range (0,75):
+        for j in range(0, 484):
+            for i in range(0, 388):
+
+                zz[j, i] = (0.5 * gdept_0[j,i]) + gdepw[k, j, i]
+                p_corr_zr[k, j, i] = p_zr[k, j, i] + grav * rhop_m[j, i] * ssh_zr[j, i] * ((gdept_0_all[j,i] - zz[j,i]) / gdept_0_all[j,i])
+                p_corr_zi[k, j, i] = p_zi[k, j, i] + grav * rhop_m[j, i] * ssh_zi[j, i] * ((gdept_0_all[j,i] - zz[j,i]) / gdept_0_all[j,i])
+
+
+    p_a = np.zeros(shape=(kk,jj,ii))
+    p_G = np.zeros(shape=(kk,jj,ii))
+
+
+    for k in range (0,75):
+        for j in range (0,484):
+            for i in range (0,388):
+                p_a[k, j, i] = (p_corr_zr[k, j, i] * p_corr_zr[k, j, i] + p_corr_zi[k, j, i] * p_corr_zi[k, j, i])**(0.5)
+                p_G[k, j, i] = math.degrees(math.atan2(p_corr_zi[k, j, i], p_corr_zr[k, j, i]))
+
+    p_a[tmask == 0] = nc.default_fillvals['f4']
+    p_G[tmask == 0] = nc.default_fillvals['f4']
+    p_corr_a[:,:,:] = p_a
+    p_corr_G[:,:,:] = p_G
+
+
+
+
+#===============================================================================
+# clean-up
+#===============================================================================
+
+    fichierpression.close()
+    fichierssh.close()
+    fichierrhop.close()
+    fichiermask.close()
+    ncfile.close()
+
+
+#*******************************************************************************
+if __name__ == '__main__':
+  #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit
+  sys.exit(main(sys.argv))
+

example/calcul_pressure-lagrangien.py

+#!/usr/bin/python
+
+import sys
+import numpy as np
+import netCDF4 as nc
+from numpy import arange, dtype
+
+#*******************************************************************************
+def main(args):
+
+    """USE
+  {0} mask.nc input.nc output.nc
+
+DESCRIPTION
+  Computes pressure from density and e3t (that changes with SSH).
+
+"""
+    if len(args)<4 or args[1]=='-h' or args[1]=='--help':
+        print( main.func_doc.format(args[0]) )
+        return 0
+
+# ===============================================================================
+# mask file
+# ===============================================================================
+
+    fichiermask=args[1]
+
+    mask_file=nc.Dataset(fichiermask)
+    depth=mask_file.variables['gdept_1d'][:]
+    tmask=mask_file.variables['tmask'][0,:,:,:]
+    mask_file.close()
+
+#===============================================================================
+# input
+#===============================================================================
+
+    fichierjulien=args[2]
+    lefichier = nc.Dataset(fichierjulien)
+
+    e3tv=lefichier.variables['e3t']
+    ji=e3tv.shape[3]
+    jj=e3tv.shape[2]
+    jk=e3tv.shape[1]
+    jl=e3tv.shape[0]
+
+    e3wv=lefichier.variables['e3w']
+    rhopv=lefichier.variables['rhop']
+    itv=lefichier.variables['time_counter']
+
+    lon=lefichier.variables['nav_lon_grid_T'][:]
+    lat=lefichier.variables['nav_lat_grid_T'][:]
+
+#===============================================================================
+# output
+#===============================================================================
+
+    p_file=args[3]
+    print('Creating '+p_file)
+
+    # this will overwrite the file
+    ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC')
+    setattr(ncfile,'history',' '.join(args))
+    ncfile.createDimension('x',ji)
+    ncfile.createDimension('y',jj)
+    ncfile.createDimension('dept',jk)
+    ncfile.createDimension('time') # unlimited
+
+    data = ncfile.createVariable('lont',dtype('float64').char,('y','x'))
+    data[:] = lon
+    setattr(data, 'units',"degrees_north")
+    setattr(data, 'standard_name',"longitude_at_T_location")
+
+    data = ncfile.createVariable('latt',dtype('float64').char,('y','x'))
+    data[:] = lat
+    setattr(data, 'units',"degrees_east")
+    setattr(data, 'standard_name',"latitude_at_T_location")
+
+    data = ncfile.createVariable('dept',dtype('float64').char,('dept'))
+    data[:] = depth
+    setattr(data, 'units',"m")
+    setattr(data, 'positive',"down")
+    setattr(data, 'standard_name',"depth_at_T_location")
+    setattr(data, 'axis',"Z")
+
+    otv = ncfile.createVariable('time',dtype('float64').char,('time'))
+    setattr(otv, 'units',itv.units)
+    setattr(otv, 'standard_name',"time")
+    setattr(otv, 'axis',"T")
+
+    pv = ncfile.createVariable('Pbc',dtype('float32').char,('time','dept','y','x'))
+    setattr(pv, 'standard_name',"Amplitude Of Pressure Total")
+    setattr(pv, 'units',"m-1.kg.s-2")
+    setattr(pv, 'coordinates', "time dept latt lont")
+
+#===============================================================================
+# process
+#===============================================================================
+
+    print( 'Dimensions: jl={} jk={} jj={} ji={}'.format(jl,jk,jj,ji) )
+
+    p_bc=np.zeros(shape=(jk,jj,ji))
+
+    grav=9.81
+
+    for l in range(0,jl):
+
+        time=itv[l]
+        print( 'Time frame {}/{}'.format(l,jl) )
+
+        e3t=e3tv[l, 0, :, :]
+        e3w=e3wv[l, :, :, :]
+        rhop=rhopv[l, :, :, :]
+
+        for k in range(0,jk):
+            if k==0:
+                p0=0.5*grav*e3t*rhop[k,:,:]
+            else:
+                p0+=0.5*grav*e3w[k,:,:]*(rhop[k,:,:]+rhop[k-1,:,:])
+
+            p0*=tmask[k,:,:]
+            p_bc[k,:,:]=p0
+
+        p_bc[tmask==0]=nc.default_fillvals['f4']
+        pv[l,:,:,:] = p_bc
+        otv[l] = time
+
+#===============================================================================
+# clean-up
+#===============================================================================
+
+    lefichier.close()
+    ncfile.close()
+
+
+#*******************************************************************************
+if __name__ == '__main__':
+  #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit
+  sys.exit(main(sys.argv))

example/calcul_pressure_eulerien.py

+#!/usr/bin/python
+
+import sys
+import numpy as np
+import netCDF4 as nc
+from numpy import arange, dtype
+
+#*******************************************************************************
+def main(args):
+
+    """USE
+  {0} mask.nc input.nc output.nc
+
+DESCRIPTION
+  Computes pressure from density and e3t (that changes with SSH).
+
+"""
+    if len(args)<4 or args[1]=='-h' or args[1]=='--help':
+        print( main.func_doc.format(args[0]) )
+        return 0
+
+# ===============================================================================
+# mask file
+# ===============================================================================
+
+    fichiermask=args[1]
+
+    mask_file=nc.Dataset(fichiermask)
+    depth=mask_file.variables['gdept_1d'][:]
+    tmask=mask_file.variables['tmask'][0,:,:,:]
+    e3t = mask_file.variables['e3t'][:,:,:,:]
+    e3w = mask_file.variables['e3w'][:,:,:,:]
+    mask_file.close()
+
+#===============================================================================
+# input
+#===============================================================================
+
+    fichierjulien=args[2]
+    lefichier = nc.Dataset(fichierjulien)
+
+    rhopv = lefichier.variables['rhop']
+    ji=rhopv.shape[3]
+    jj=rhopv.shape[2]
+    jk=rhopv.shape[1]
+    jl=rhopv.shape[0]
+
+    ssh1 = lefichier.variables['ssh']
+
+
+
+    itv=lefichier.variables['time_counter']
+
+    lon=lefichier.variables['nav_lon_grid_T'][:]
+    lat=lefichier.variables['nav_lat_grid_T'][:]
+
+#===============================================================================
+# output
+#===============================================================================
+
+    p_file=args[3]
+    print('Creating '+p_file)
+
+    # this will overwrite the file
+    ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC')
+    setattr(ncfile,'history',' '.join(args))
+    ncfile.createDimension('x',ji)
+    ncfile.createDimension('y',jj)
+    ncfile.createDimension('dept',jk)
+    ncfile.createDimension('time') # unlimited
+
+    data = ncfile.createVariable('lont',dtype('float64').char,('y','x'))
+    data[:] = lon
+    setattr(data, 'units',"degrees_north")
+    setattr(data, 'standard_name',"longitude_at_T_location")
+
+    data = ncfile.createVariable('latt',dtype('float64').char,('y','x'))
+    data[:] = lat
+    setattr(data, 'units',"degrees_east")
+    setattr(data, 'standard_name',"latitude_at_T_location")
+
+    data = ncfile.createVariable('dept',dtype('float64').char,('dept'))
+    data[:] = depth
+    setattr(data, 'units',"m")
+    setattr(data, 'positive',"down")
+    setattr(data, 'standard_name',"depth_at_T_location")
+    setattr(data, 'axis',"Z")
+
+    otv = ncfile.createVariable('time',dtype('float64').char,('time'))
+    setattr(otv, 'units',itv.units)
+    setattr(otv, 'standard_name',"time")
+    setattr(otv, 'axis',"T")
+
+    pv = ncfile.createVariable('Ptot',dtype('float32').char,('time','dept','y','x'))
+    setattr(pv, 'standard_name',"Amplitude Of Pressure Total")
+    setattr(pv, 'units',"m-1.kg.s-2")
+    setattr(pv, 'coordinates', "time dept latt lont")
+
+#===============================================================================
+# process
+#===============================================================================
+
+    print( 'Dimensions: jl={} jk={} jj={} ji={}'.format(jl,jk,jj,ji) )
+
+    p_tot=np.zeros(shape=(jk,jj,ji))
+
+    grav=9.81
+
+    for l in range(0,jl):
+
+        time=itv[l]
+        print( 'Time frame {}/{}'.format(l,jl) )
+
+        e3t1=e3t[0, 0, :, :]
+        e3w1=e3w[0, :, :, :]
+        rhop=rhopv[l, :, :, :]
+        ssh=ssh1[l, :, :]
+
+        for k in range(0,jk):
+            if k==0:
+                p0=0.5*grav*e3t1*rhop[k,:,:]+(rhop[k,:,:]*grav*ssh[:,:])
+            else:
+                p0+=0.5*grav*e3w1[k,:,:]*(rhop[k,:,:]+rhop[k-1,:,:])
+
+            p0*=tmask[k,:,:]
+            p_tot[k,:,:]=p0
+
+        p_tot[tmask==0]=nc.default_fillvals['f4']
+        pv[l,:,:,:] = p_tot
+        otv[l] = time
+
+#===============================================================================
+# clean-up
+#===============================================================================
+
+    lefichier.close()
+    ncfile.close()
+
+
+#*******************************************************************************
+if __name__ == '__main__':
+  #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit
+  sys.exit(main(sys.argv))

example/diff_eulerien-lagrangien.py

+#!/usr/bin/python
+
+import sys
+import numpy as np
+import netCDF4 as nc
+from numpy import arange, dtype
+import math
+
+#*******************************************************************************
+def main(args):
+
+    """USE
+  {0} mask.nc input.nc output.nc
+
+DESCRIPTION
+  Computes pressure from density and e3t (that changes with SSH).
+
+"""
+    if len(args)<4 or args[1]=='-h' or args[1]=='--help':
+        print( main.func_doc.format(args[0]) )
+        return 0
+
+
+
+#===============================================================================
+# input
+#===============================================================================
+
+    file_pressure_eulerien = args[1]
+    fichierpression1 = nc.Dataset(file_pressure_eulerien)
+
+    ptot_a=fichierpression1.variables['Ptot_a'][:,:,:]
+    ptot_G=fichierpression1.variables['Ptot_G'][:,:,:]
+    depth=fichierpression1.variables['dept'][:]
+    lon=fichierpression1.variables['lont'][:]
+    lat=fichierpression1.variables['latt'][:]
+
+    jj=ptot_a.shape[1]
+    ii=ptot_a.shape[2]
+    kk=ptot_a.shape[0]
+
+
+    file_pressure_lagpatch = args[2]
+    fichierpression2 = nc.Dataset(file_pressure_lagpatch)
+    P_corr_A=fichierpression2.variables['P_corr_A'][:,:,:]
+    P_corr_G=fichierpression2.variables['P_corr_G'][:,:,:]
+
+
+
+#===============================================================================
+# output
+#===============================================================================
+
+    p_file=args[3]
+    print('Creating '+p_file)
+
+    # this will overwrite the file
+    ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC')
+    setattr(ncfile,'history',' '.join(args))
+    ncfile.createDimension('x',ii)
+    ncfile.createDimension('y',jj)
+    ncfile.createDimension('dept',kk)
+
+
+    data = ncfile.createVariable('lont',dtype('float64').char,('y','x'))
+    data[:] = lon
+    setattr(data, 'units',"degrees_north")
+    setattr(data, 'standard_name',"longitude_at_T_location")
+
+    data = ncfile.createVariable('latt',dtype('float64').char,('y','x'))
+    data[:] = lat
+    setattr(data, 'units',"degrees_east")
+    setattr(data, 'standard_name',"latitude_at_T_location")
+
+    data = ncfile.createVariable('dept',dtype('float64').char,('dept'))
+    data[:] = depth
+    setattr(data, 'units',"m")
+    setattr(data, 'positive',"down")
+    setattr(data, 'standard_name',"depth_at_T_location")
+    setattr(data, 'axis',"Z")
+
+    p_diff_a = ncfile.createVariable('P_corr_A',dtype('float32').char,('dept','y','x'))
+    setattr(p_diff_a, 'standard_name',"Amplitude Of Pressure Total")
+    setattr(p_diff_a, 'units',"m-1.kg.s-2")
+    setattr(p_diff_a, 'coordinates', "dept latt lont")
+
+    p_diff_G = ncfile.createVariable('P_corr_G',dtype('float32').char,('dept','y','x'))
+    setattr(p_diff_G, 'standard_name',"Phase Of Pressure Total")
+    setattr(p_diff_G, 'units',"m-1.kg.s-2")
+    setattr(p_diff_G, 'coordinates', "dept latt lont")
+
+#===============================================================================
+# process
+#===============================================================================
+
+    p_diff_a[:,:,:] = ptot_a[:,:,:] - P_corr_A[:,:,:]
+    p_diff_G[:,:,:] = ptot_G[:,:,:] - P_corr_G[:,:,:]
+
+
+#===============================================================================
+# clean-up
+#===============================================================================
+
+    fichierpression1.close()
+    fichierpression2.close()
+    ncfile.close()
+
+
+#*******************************************************************************
+if __name__ == '__main__':
+  #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit
+  sys.exit(main(sys.argv))
+