# +-======-+ # Copyright (c) 2003-2007 United States Government as represented by # the Admistrator of the National Aeronautics and Space Administration. # All Rights Reserved. # # THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, # REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN # COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS # REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). # THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN # INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR # REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, # DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED # HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE # RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. # # Government Agency: National Aeronautics and Space Administration # Government Agency Original Software Designation: GSC-15354-1 # Government Agency Original Software Title: GEOS-5 GCM Modeling Software # User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov # Government Agency Point of Contact for Original Software: # Dale Hithon, SRA Assistant, (301) 286-2691 # # +-======-+ """ Very simple object-oriented Python interface to the GFIO library. It can be used to read and write COARDS compatible files. """ from numpy import linspace, ones, zeros from GFIO_ import * class GFIO(object): def __init__(self,filename=None,mode='r'): """ Open a file GFIO file for reading or writting. All relevant metadata regarding variables and coordinates are retrieved and stored in the object. Examine the contents of __dict__ to see what is available. """ if filename is None: self.filename = None self.fid = None return self.filename = filename # Open an existing file # --------------------- if mode == 'r': self.fid, rc = gfioopen(filename,1) else: self.fid, rc = gfioopen(filename,0) if rc: raise IOError, "cannot open GFIO file " + filename # Fetch dimensions # ---------------- self.im,self.jm,self.km,self.lm,self.nvars,self.ngatts,rc = gfiodiminquire(self.fid) if rc: raise IOError, "cannot get dimensions for GFIO file " + filename # Fetch metadata # -------------- km_ = max(self.km,1) self.title, self.source, self.contact, self.amiss, \ self.lon, self.lat, self.levs, self.levunits, \ self.yyyymmdd, self.hhmmss, self.timinc, \ self.vname, self.vtitle, self.vunits, self.kmvar, \ self.valid_range, self.packing_range,rc \ = gfioinquire(self.fid,self.im,self.jm,km_,self.lm,self.nvars) if rc==3: raise IOError, "cannot get COARDS metadata for GFIO file " \ + filename + ", rc = " + str(rc) # Create variable property dictionaries # -------------------------------------- self.vname = self.vname.split(':')[:-1] self.vtitle = self.vtitle.split(':')[:-1] self.vunits = self.vunits.split(':')[:-1] # Hash with number of levels by variable name # ------------------------------------------- self.kmvar_name = {} for i in range(self.nvars): self.kmvar_name[self.vname[i]] = self.kmvar[i] #--- def read(self,vname,nymd=None,nhms=None,kbeg=None,kount=None,squeeze=True): """ Reads a variable at a given time/date, or the first time on file if nymd/nhms is not specified. By default it returns all vertical levels, but a range can be specified with the first vertical index (kbeg) and the number of vertical levels to be read (kount). Example: ------- from gfio import GFIO f = GFIO('test_py.nc') # open existing file u = f.read('u') # read variable u """ if nymd is None: nymd, nhms, dt, rc = gfiogetbegdatetime(self.fid) if rc: raise IOError, "cannot get initial date/time for GFIO file "+self.filename if kbeg is None: try: kmv = self.kmvar_name[vname] except: raise IOError, "variable <"+vname+"> not present in GFIO file" if kmv == 0: kbeg = 0 # 2D kount = 1 else: kbeg = 1 # 3D kount = self.km var, rc = gfiogetvar(self.fid,vname,nymd,nhms,self.im,self.jm,kbeg,kount) if rc: raise IOError, "cannot read <"+vname+"> from GFIO file "+self.filename if squeeze: var = var.squeeze() return var #--- def create(self, filename, vname, nymd_beg, nhms_beg, lon=None, lat=None, refine=None, res=None, vtitle=None, vunits=None, timinc=60000, kmvar=None, levs=[1000.,], levunits='hPa', title='Produced with PyGFIO', source='NASA/GSFC/GMAO', contact='unknown', amiss=1.e+20, valid_range=None, packing_range=None, prec=0): """ Creates a GFIO dataset.The following parameters are mandatory: filename --- output file name vname --- list of variables to write to file nymd_beg --- first date on file, e.g., 20080630 nhms_beg --- first time on file, e.g., 12000000 For the *horizontal* coordinates, either lon, lat --- list of longitudes/latitudes, in degrees or refine --- refinement for the standard GEOS-5 4x5 A-Grid. For example, refine=1 produces a 4 x 5 grid refine=2 produces a 2 x2.50 grid refine=4 produces a 1 x1,25 grid refine=8 produces a 0.50x0.625 grid refine=16 produces a 0.25x0.3125 grid or even res --- single letter denoting GEOS-5 resolution, res='a' produces a 4 x 5 grid res='b' produces a 2 x2.50 grid res='c' produces a 1 x1,25 grid res='d' produces a 0.50x0.625 grid res='e' produces a 0.25x0.3125 grid For 3D datasets, you must specify: levs --- vertical levels levunits --- units of the vertical levels, per COARDS conventions For datasets that will have more than 1 time step in it, you must specify timinc --- time step, in HHMMSS format (integer) It is also recommended that you specify vtitle --- list of variable long names (default=vname) vunits --- list of variable units (default='unknown') Optional parameters: amiss --- missing values valid_range --- array of shape (2,nvars) Variable valid range; method write() will return an error if a value is outside of this range. IMPORTANT: If packing is not desired for a given variable, YOU MUST set both components of valid_range to amiss; this is the default. packing_range --- array of shape (2,nvars) Packing range to be used for 16-bit packing of each variable. IMPORTANT: If packing is not desired for a given variable, YOU MUST set both components of packing_range to amiss. NOTE: * The packing algorithm sets all values outside the packing range to missing. * The larger the packing range, the greater the loss of precision. prec --- Desired precision of data: 0 = 32 bit 1 = 64 bit Consult the GFIO documentation for additional information. Example: ------- from gfio import GFIO filename = 'test_py.nc' vname = ['u','v'] vtitle = ['zonal wind','meridional wind'] vunits = ['m/s','m/s'] nymd = 20080630 nhms = 120000 f = GFIO() f.create(filename, vname, nymd, nhms, res='c', vtitle=vtitle, vunits=vunits) u = zeros((f.im,f.jm)) v = ones((f.im,f.jm)) f.write('u',nymd,nhms,u) f.write('v',nymd,nhms,v) f.close() """ self.filename = filename nvars = len(vname) # Defaults # -------- if lon is None or lat is None: if res is not None: if res=='a': refine = 1 if res=='b': refine = 2 if res=='c': refine = 4 if res=='d': refine = 8 if res=='e': refine = 16 if refine is None: raise IOError, 'must specify either lat/lon, res or refine' dx = 5. / refine dy = 4. / refine im = int(360. / dx) jm = int(180. / dy + 1) lon = linspace(-180.,180.,im,endpoint=False) lat = linspace(-90.,90.,jm) if vtitle is None: vtitle = vname if vunits is None: vunits = nvars * ['unknown',] if valid_range is None: valid_range = amiss * ones((2,nvars)) if packing_range is None: packing_range = amiss * ones((2,nvars)) if kmvar is None: kmvar = zeros(nvars) # Do not know howw to pass string array to f2py # --------------------------------------------- vname_ = ":".join(vname) + ':' vtitle_ = ":".join(vtitle) + ':' vunits_ = ":".join(vunits) + ':' self.fid, rc = gfiocreate(filename,title,source,contact,amiss,\ lon,lat,levs,levunits,\ nymd_beg,nhms_beg,timinc,\ vname_,vtitle_,vunits_,kmvar,\ valid_range,packing_range,prec) if rc: raise IOError, "cannot create "+self.filename+", rc=%d"%rc # Save relevant metadata # ---------------------- self.im, self.jm = (lon.size,lat.size) self.lon, self.lat = (lon,lat) self.levs, self.levunits = (levs,levunits) self.nvars = nvars self.vname, self.vtitle, self.vunits = (vname,vtitle,vunits) self.title, self.source, self.contact = (title,source,contact) self.kmvar = kmvar self.amiss = amiss self.timinc = timinc self.valid_range = valid_range self.packing_range = packing_range self.nymd_beg, self.nhms_beg = (nymd_beg,nhms_beg) # Hash with number of levels by variable name # ------------------------------------------- self.kmvar_name = {} for i in range(self.nvars): self.kmvar_name[self.vname[i]] = self.kmvar[i] #--- def write(self,vname,nymd,nhms,var,kbeg=None,kount=None): """ Writes a variable to file at a given date/time. See method create() for an example. """ # 2D variables # ------------ if self.kmvar_name[vname] == 0: kbeg = 0 kount = 1 # 3D variables # ------------ else: if kbeg is None: kbeg = 1 if kount is None: kount = self.kmvar_name[vname] rc = gfioputvar(self.fid,vname,nymd,nhms,kbeg,var) if rc: raise IOError, "cannot write <%s> to file %s, rc=%d "%\ (vname,self.filename,rc) #--- def interp(self,vname,lon,lat,nymd=None,nhms=None,kbeg=None,kount=None,squeeze=True): """ Reads a variable at a given time/date, or the first time on file if nymd/nhms is not specified, and interpolates it to the observation locations given by the list of (lon,lat) on input. By default it returns all vertical levels of the variable, but a range can be specified with the first vertical index (kbeg) and the number of vertical levels to be read (kount). """ # First read the gridded field # ---------------------------- var = self.read(vname,nymd,nhms,kbeg,kount,squeeze=False) # Check consistency of longitudes # ------------------------------- if any(lon<0): if any(self.lon>180.): raise ValueError, \ "inconsistent longitudes, obs in [-180,180] but grid in [0,360]" elif any(lon>180): if any(self.lon<0.): raise ValueError, \ "inconsistent longitudes, obs in [0,360] but grid in [-180,180]" # Next interpolate (assumes global, zonally periodic grids) # --------------------------------------------------------- ivar = gfiointerpxy(lon,lat,self.lon[0],var) if squeeze: ivar = ivar.squeeze() return ivar #--- def close(self): """Closes a GFIO file.""" if self.fid is not None: rc = gfioclose(self.fid) #--- def __del__(self): if self.fid is not None: rc = gfioclose(self.fid) #........................................................................... if __name__ == "__main__": # f = GFIO('sample.nc') filename = 'test_py.nc' vname = ['u','v'] vtitle = ['zonal wind','meridional wind'] vunits = ['m/s','m/s'] nymd = 20080630 nhms = 120000 f = GFIO() f.create(filename, vname, nymd, nhms, res='c', vtitle=vtitle, vunits=vunits) u = zeros((f.im,f.jm)) v = ones((f.im,f.jm)) f.write('u',nymd,nhms,u) f.write('v',nymd,nhms,v) f.close()