Recipe: python program reads cubed-sphere data

From GEOS-5
Revision as of 11:29, 6 March 2017 by Weiyuan.jiang (talk | contribs) (Created page with "import sys from netCDF4 import Dataset import numpy as np def get_csdata(filename,var,Time,Level=None): fh = Dataset(filename, mode='r') Times = fh.dimensions['time'...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

import sys from netCDF4 import Dataset import numpy as np

def get_csdata(filename,var,Time,Level=None):

   fh = Dataset(filename, mode='r')
   Times = fh.dimensions['time'].size
   # Levs=fh.dimensions['lev'].size
   nf = fh.dimensions['nf'].size
   Ydim = fh.dimensions['Ydim'].size
   Xdim = fh.dimensions['Xdim'].size
   if Level is None:
       T = fh.variables[var][Time, :, :, :]
   else:
       T = fh.variables[var][Time, Level, :, :, :]
   lon3d = fh.variables['lons'][:, :, :]
   lat3d = fh.variables['lats'][:, :, :]
   fh.close()
   length = nf * Ydim * Xdim
   lat3d = np.reshape(lat3d, (length,))
   lon3d = np.reshape(lon3d, (length,))
   T = np.reshape(T, (length,))
   Tm = T
   # get rid of non values
   if (hasattr(T, 'mask')):
       Tm = T[~T.mask]
       lat3d = lat3d[~T.mask]
       lon3d = lon3d[~T.mask]
   return lon3d, lat3d, Tm  # data without poles
   # use this function if the users like to fill the poles

def get_csdataNS(filename,var,Time,Level=None):

   fh = Dataset(filename, mode='r')
   Times = fh.dimensions['time'].size
   # Levs=fh.dimensions['lev'].size
   nf = fh.dimensions['nf'].size
   Ydim = fh.dimensions['Ydim'].size
   Xdim = fh.dimensions['Xdim'].size
   if Level is None:
       T = fh.variables[var][Time, :, :, :]
   else:
       T = fh.variables[var][Time, Level, :, :, :]
   lon3d = fh.variables['lons'][:, :, :]
   lat3d = fh.variables['lats'][:, :, :]
   fh.close()
   # nearest north pole data
   dataN=np.mean(T[2,(Xdim-1)//2:Xdim//2,(Ydim-1)//2:Ydim//2])
   # nearest south pole data
   dataS=np.mean(T[5,(Xdim-1)//2:Xdim//2,(Ydim-1)//2:Ydim//2])
   length = nf * Ydim * Xdim
   lat3d = np.reshape(lat3d, (length,))
   lon3d = np.reshape(lon3d, (length,))    T = np.reshape(T, (length,))
   Tm = T
   # get rid of non values
   if (hasattr(T, 'mask')):
       Tm = T[~T.mask]
       lat3d = lat3d[~T.mask]
       lon3d = lon3d[~T.mask]
   # add the data on south pole and north pole
   lats=np.copy(lat3d)
   lons=np.copy(lon3d)
   data=np.copy(Tm)
   length = lats.size
   lats.resize(length+4)
   lons.resize(length+4)
   data.resize(length+4)
   lats[length]=90
   lons[length]=0
   data[length]=dataN
   lats[length+1]=90
   lons[length+1]=360
   data[length+1]=dataN
   lats[length+2]=-90
   lons[length+2]=0
   data[length+2]=dataS
   lats[length+3]=-90
   lons[length+3]=360
   data[length+3]=dataS
   return lons,lats,data

</syntaxhighlight>