Recipe: python program reads cubed-sphere data

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