import os import netCDF4 import csv import re import textwrap from datetime import datetime from numpy import any as numpy_any from numpy import arange as numpy_arange from numpy import arccos as numpy_arccos from numpy import arcsin as numpy_arcsin from numpy import array as numpy_array from numpy import clip as numpy_clip from numpy import column_stack as numpy_column_stack from numpy import cos as numpy_cos from numpy import deg2rad as numpy_deg2rad from numpy import dtype as numpy_dtype from numpy import empty as numpy_empty from numpy import mean as numpy_mean from numpy import pi as numpy_pi from numpy import rad2deg as numpy_rad2deg from numpy import resize as numpy_resize from numpy import result_type as numpy_result_type from numpy import sin as numpy_sin from numpy import transpose as numpy_transpose from numpy import where as numpy_where from netCDF4 import date2num as netCDF4_date2num from .. import __version__, __Conventions__, __file__ from ..domain import Domain from ..coordinatereference import CoordinateReference from ..field import Field, FieldList from ..cellmethods import CellMethods from ..cfdatetime import Datetime from ..coordinate import DimensionCoordinate, AuxiliaryCoordinate from ..functions import RTOL, ATOL, equals from ..units import Units from ..functions import (open_files_threshold_exceeded, close_one_file, abspath) from ..data.data import Data, Partition, PartitionMatrix from .filearray import UMFileArray from .functions import _open_um_file from .umread.umfile import UMFileException # -------------------------------------------------------------------- # Constants # -------------------------------------------------------------------- _pi_over_180 = numpy_pi/180.0 # PP missing data indicator _pp_rmdi = -1.0e+30 # Reference surface pressure in Pascals _pstar = 1.0e5 # -------------------------------------------------------------------- # Characters used in decoding LBEXP into a runid # -------------------------------------------------------------------- _characters = ('a','b','c','d','e','f','g','h','i','j','k','l','m', 'n','o','p','q','r','s','t','u','v','w','x','y','z', '0','1','2','3','4','5','6','7','8','9') _n_characters = len(_characters) # -------------------------------------------------------------------- # Number matching regular expression # -------------------------------------------------------------------- _number_regex = '([-+]?\d*\.?\d+(e[-+]?\d+)?)' # -------------------------------------------------------------------- # Date-time object that copes with non-standard calendars # -------------------------------------------------------------------- netCDF4_netcdftime_datetime = netCDF4.netcdftime.datetime # -------------------------------------------------------------------- # Caches for derived values # -------------------------------------------------------------------- _cache_latlon = {} _cached_runid = {} _cached_latlon = {} _cached_time = {} _cached_ctime = {} _cached_size_1_height_coordinate = {} _cached_z_coordinate = {} _cached_date2num = {} _cached_model_level_number_coordinate = {} _Units = { None : Units(), '' : Units(''), '1' : Units('1'), 'Pa' : Units('Pa'), 'm' : Units('m'), 'hPa' : Units('hPa'), 'K' : Units('K'), 'degrees' : Units('degrees'), 'degrees_east' : Units('degrees_east'), 'degrees_north' : Units('degrees_north'), 'days' : Units('days'), 'gregorian 1752-9-13': Units('days since 1752-9-13', 'gregorian'), '365_day 1752-9-13' : Units('days since 1752-9-13', '365_day'), '360_day 0-1-1' : Units('days since 0-1-1', '360_day'), } # -------------------------------------------------------------------- # Names of PP integer and real header items # -------------------------------------------------------------------- _header_names = ('LBYR', 'LBMON', 'LBDAT', 'LBHR', 'LBMIN', 'LBDAY', 'LBYRD', 'LBMOND', 'LBDATD', 'LBHRD', 'LBMIND', 'LBDAYD', 'LBTIM', 'LBFT', 'LBLREC', 'LBCODE', 'LBHEM', 'LBROW', 'LBNPT', 'LBEXT', 'LBPACK', 'LBREL', 'LBFC', 'LBCFC', 'LBPROC', 'LBVC', 'LBRVC', 'LBEXP', 'LBEGIN', 'LBNREC', 'LBPROJ', 'LBTYP', 'LBLEV', 'LBRSVD1', 'LBRSVD2', 'LBRSVD3', 'LBRSVD4', 'LBSRCE', 'LBUSER1', 'LBUSER2', 'LBUSER3', 'LBUSER4', 'LBUSER5', 'LBUSER6', 'LBUSER7', 'BRSVD1', 'BRSVD2', 'BRSVD3', 'BRSVD4', 'BDATUM', 'BACC', 'BLEV', 'BRLEV', 'BHLEV', 'BHRLEV', 'BPLAT', 'BPLON', 'BGOR', 'BZY', 'BDY', 'BZX', 'BDX', 'BMDI', 'BMKS') # -------------------------------------------------------------------- # Positions of PP header items in their arrays # -------------------------------------------------------------------- (lbyr, lbmon, lbdat, lbhr, lbmin, lbday, lbyrd, lbmond, lbdatd, lbhrd, lbmind, lbdayd, lbtim, lbft, lblrec, lbcode, lbhem, lbrow, lbnpt, lbext, lbpack, lbrel, lbfc, lbcfc, lbproc, lbvc, lbrvc, lbexp, lbegin, lbnrec, lbproj, lbtyp, lblev, lbrsvd1, lbrsvd2, lbrsvd3, lbrsvd4, lbsrce, lbuser1, lbuser2, lbuser3, lbuser4, lbuser5, lbuser6, lbuser7, ) = range(45) (brsvd1, brsvd2, brsvd3, brsvd4, bdatum, bacc, blev, brlev, bhlev, bhrlev, bplat, bplon, bgor, bzy, bdy, bzx, bdx, bmdi, bmks, ) = range(19) # -------------------------------------------------------------------- # Assign CF standard name attributes to PP axis codes. (The full list # of field code keys may be found at # http://cms.ncas.ac.uk/html_umdocs/wave/@header.) # -------------------------------------------------------------------- _coord_standard_name = { 0 : None, # Sigma (or eta, for hybrid coordinate data). 1 : 'air_pressure', # Pressure (mb). 2 : 'altitude', # Height above sea level (km). 3 : 'atmosphere_hybrid_sigma_pressure_coordinate', # Eta (U.M. hybrid coordinates) only. 4 : 'depth', # Depth below sea level (m) 5 : 'model_level_number', # Model level. 6 : 'air_potential_temperature', # Theta 7 : 'atmosphere_sigma_coordinate', # Sigma only. 8 : None, # Sigma-theta 10 : 'latitude', # Latitude (degrees N). 11 : 'longitude', # Longitude (degrees E). 13 : 'region', # Site number (set of parallel rows or columns e.g.Time series) 14 : 'atmosphere_hybrid_height_coordinate', 15 : 'height', 20 : 'time', # Time (days) (Gregorian calendar (not 360 day year)) 21 : 'time', # Time (months) 22 : 'time', # Time (years) 23 : 'time', # Time (model days with 360 day model calendar) 40 : None, # pseudolevel 99 : None, # Other -10 : 'grid_latitude', # Rotated latitude (degrees). -11 : 'grid_longitude', # Rotated longitude (degrees). -20 : 'radiation_wavelength', } # -------------------------------------------------------------------- # Assign CF long names to PP axis codes. # -------------------------------------------------------------------- _coord_long_name = {} # -------------------------------------------------------------------- # Assign CF units strings to PP axis codes. # -------------------------------------------------------------------- #_coord_units = { _axiscode_to_units = { 0 : '1', # Sigma (or eta, for hybrid coordinate data) 1 : 'hPa', # air_pressure 2 : 'm', # altitude 3 : '1', # atmosphere_hybrid_sigma_pressure_coordinate 4 : 'm', # depth 5 : '1', # model_level_number 6 : 'K', # air_potential_temperature 7 : '1', # atmosphere_sigma_coordinate 10 : 'degrees_north', # latitude 11 : 'degrees_east', # longitude 13 : '', # region 14 : '1', # atmosphere_hybrid_height_coordinate 15 : 'm', # height 20 : 'days', # time (gregorian) 23 : 'days', # time (360_day) 40 : '1', # pseudolevel -10 : 'degrees', # rotated latitude (not an official axis code) -11 : 'degrees', # rotated longitude (not an official axis code) } _axiscode_to_Units = { 0 : _Units['1'], # Sigma (or eta, for hybrid coordinate data) 1 : _Units['hPa'], # air_pressure 2 : _Units['m'], # altitude 3 : _Units['1'], # atmosphere_hybrid_sigma_pressure_coordinate 4 : _Units['m'], # depth 5 : _Units['1'], # model_level_number 6 : _Units['K'], # air_potential_temperature 7 : _Units['1'], # atmosphere_sigma_coordinate 10 : _Units['degrees_north'], # latitude 11 : _Units['degrees_east'], # longitude 13 : _Units[''], # region 14 : _Units['1'], # atmosphere_hybrid_height_coordinate 15 : _Units['m'], # height 20 : _Units['days'], # time (gregorian) 23 : _Units['days'], # time (360_day) 40 : _Units['1'], # pseudolevel -10 : _Units['degrees'], # rotated latitude (not an official axis code) -11 : _Units['degrees'], # rotated longitude (not an official axis code) } # -------------------------------------------------------------------- # Assign CF axis attributes to PP axis codes. # -------------------------------------------------------------------- _coord_axis = { 1 : 'Z', # air_pressure 2 : 'Z', # altitude 3 : 'Z', # atmosphere_hybrid_sigma_pressure_coordinate 4 : 'Z', # depth 5 : 'Z', # model_level_number 6 : 'Z', # air_potential_temperature 7 : 'Z', # atmosphere_sigma_coordinate 10 : 'Y', # latitude 11 : 'X', # longitude 13 : None, # region 14 : 'Z', # atmosphere_hybrid_height_coordinate 15 : 'Z', # height 20 : 'T', # time (gregorian) 23 : 'T', # time (360_day) 40 : None, # pseudolevel -10 : 'Y', # rotated latitude (not an official axis code) -11 : 'X', # rotated longitude (not an official axis code) } # -------------------------------------------------------------------- # Assign CF positive attributes to PP axis codes. # -------------------------------------------------------------------- _coord_positive = { 1 : 'down', # air_pressure 2 : 'up', # altitude 3 : 'down', # atmosphere_hybrid_sigma_pressure_coordinate 4 : 'down', # depth 5 : None, # model_level_number 6 : 'up', # air_potential_temperature 7 : 'down', # atmosphere_sigma_coordinate 10 : None, # latitude 11 : None, # longitude 13 : None, # region 14 : 'up', # atmosphere_hybrid_height_coordinate 15 : 'up', # height 20 : None, # time (gregorian) 23 : None, # time (360_day) 40 : None, # pseudolevel -10 : None, # rotated latitude (not an official axis code) -11 : None, # rotated longitude (not an official axis code) } # -------------------------------------------------------------------- # Translate LBVC codes to PP axis codes. (The full list of field code # keys may be found at # http://cms.ncas.ac.uk/html_umdocs/wave/@fcodes.) # -------------------------------------------------------------------- _lbvc_to_axiscode = { 1 : 2, # altitude (Height) 2 : 4, # depth (Depth) 3 : None, # (Geopotential (= g*height)) 4 : None, # (ICAO height) 6 : 5, # model_level_number 7 : None, # (Exner pressure) 8 : 1, # air_pressure (Pressure) 9 : 3, # atmosphere_hybrid_sigma_pressure_coordinate (Hybrid pressure) 10 : 7, # atmosphere_sigma_coordinate (Sigma (= p/surface p)) ## dch check 16 : None, # (Temperature T) 19 : 6, # air_potential_temperature (Potential temperature) 27 : None, # (Atmospheric) density 28 : None, # (d(p*)/dt . p* = surface pressure) 44 : None, # (Time in seconds) 65 : 14, # atmosphere_hybrid_height_coordinate (Hybrid height) 129 : None, # Surface 176 : 10, # latitude (Latitude) 177 : 11, # longitude (Longitude) } # -------------------------------------------------------------------- # Names of PP extra data codes # -------------------------------------------------------------------- _extra_data_name = { 1 : 'x', 2 : 'y', 3 : 'y_domain_lower_bound', 4 : 'x_domain_lower_bound', 5 : 'y_domain_upper_bound', 6 : 'x_domain_upper_bound', 7 : 'z_domain_lower_bound', 8 : 'x_domain_upper_bound', 9 : 'title', 10 : 'domain_title', 11 : 'x_lower_bound', 12 : 'x_upper_bound', 13 : 'y_lower_bound', 14 : 'y_upper_bound', } # -------------------------------------------------------------------- # Model identifier codes. These are the the last four digits of # LBSRCE. # -------------------------------------------------------------------- _lbsrce_model_codes = {1111 : 'UM'} # -------------------------------------------------------------------- # LBCODE values for unrotated latitude longitude grids # -------------------------------------------------------------------- _true_latitude_longitude_lbcodes = set((1, 2)) # -------------------------------------------------------------------- # LBCODE values for rotated latitude longitude grids # -------------------------------------------------------------------- _rotated_latitude_longitude_lbcodes = set((101, 102, 111)) _axis = {'t' : 'dim0', 'z' : 'dim1', 'y' : 'dim2', 'x' : 'dim3', 'r' : 'dim4', 'p' : 'dim5', 'area': None, } class UMField(object): ''' ''' _debug = False def __init__(self, var, fmt, byte_ordering, word_size, um_version, set_standard_name, height_at_top_of_model, **kwargs): ''' **Initialization** :Parameters: var : cf.um.umread.umfile.Var byte_ordering : str 'little_endian' or 'big_endian' word_size : int Word size in bytes (4 or 8). fmt : str 'PP' or 'FF' um_version : number set_standard_name : bool If True then set the standard_name CF property. height_at_top_of_model : float, optional **kwargs : *optional* Keyword arguments specifying CF properties for the UM field. ''' self._nonzero = False self.fmt = fmt self.height_at_top_of_model = height_at_top_of_model self.byte_ordering = byte_ordering self.word_size = word_size self.atol = ATOL() self.domain = Domain() cf_properties = {} attributes = {} self.fields = [] filename = abspath(var.file.path) self.filename = filename groups = var.group_records_by_extra_data() n_groups = len(groups) if n_groups == 1: # There is one group of records groups_nz = [var.nz] groups_nt = [var.nt] elif n_groups > 1: # There are multiple groups of records, distinguished by # different extra data. groups_nz = [] groups_nt = [] groups2 = [] for group in groups: group_size = len(group) if group_size == 1: # There is only one record in this group split_group= False nz = 1 elif group_size > 1: # There are multiple records in this group # Find the lengths of runs of identical times times = [(self.header_vtime(rec), self.header_dtime(rec)) for rec in group] lengths = [len(tuple(g)) for k, g in itertools.groupby(times)] if len(set(lengths)) == 1: # Each run of identical times has the same # length, so it is possible that this group # forms a variable of nz x nt records. split_group = False nz = lengths.pop() z0 = [self.z for rec in group[:nz]] for i in range(nz, group_size, nz): z1 = [self.header_z(rec) for rec in group[i:i+nz]] if z1 != z0: split_group = True break else: # Different runs of identical times have # different lengths, so it is not possible for # this group to form a variable of nz x nt # records. split_group = True nz = 1 #--- End: if if split_group: # This group doesn't form a complete nz x nt # matrix, so split it up into 1 x 1 groups. groups2.extend([[rec] for rec in group]) groups_nz.extend([1] * group_size) groups_nt.extend([1] * group_size) else: # This group forms a complete nz x nt matrix, so # it may be considered as a variable in its own # right and doesn't need to be split up. groups2.append(group) groups_nz.append(nz) groups_nt.append(group_size/nz) #--- End: for groups = groups2 #--- End: if rec0 = groups[0][0] int_hdr = rec0.int_hdr self.int_hdr_dtype = int_hdr.dtype int_hdr = int_hdr.tolist() real_hdr = rec0.real_hdr.tolist() self.int_hdr = int_hdr self.real_hdr = real_hdr # ------------------------------------------------------------ # Set some metadata quantities which are guaranteed to be the # same for all records in a variable # ------------------------------------------------------------ LBNPT = int_hdr[lbnpt] LBROW = int_hdr[lbrow] LBTIM = int_hdr[lbtim] LBCODE = int_hdr[lbcode] LBPROC = int_hdr[lbproc] LBVC = int_hdr[lbvc] LBUSER5 = int_hdr[lbuser5] BPLAT = real_hdr[bplat] BPLON = real_hdr[bplon] BDX = real_hdr[bdx] BDY = real_hdr[bdy] self.lbnpt = LBNPT self.lbrow = LBROW self.lbtim = LBTIM self.lbproc = LBPROC self.lbvc = LBVC self.bplat = BPLAT self.bplon = BPLON self.bdx = BDX self.bdy = BDY # ------------------------------------------------------------ # Set some derived metadata quantities which are (as good as) # guaranteed to be the same for all records in a variable # ------------------------------------------------------------ self.lbtim_ia, ib = divmod(LBTIM, 100) self.lbtim_ib, ic = divmod(ib, 10) if ic == 1: calendar = 'gregorian' elif ic == 4: calendar = '365_day' else: calendar = '360_day' self.calendar = calendar self.reference_time_Units() header_um_version, source = divmod(int_hdr[lbsrce], 10000) if header_um_version > 0 and int(um_version) == um_version: #len(um_version) <= 3: # header_um_version = str(header_um_version) model_um_version = header_um_version self.um_version = header_um_version else: model_um_version = None self.um_version = um_version # Set source source = _lbsrce_model_codes.setdefault(source, None) if source is not None and model_um_version is not None: source += ' vn%s' % model_um_version if source: cf_properties['source'] = source # ------------------------------------------------------------ # Set the T, Z, Y and X axis codes. These are guaranteed to be # the same for all records in a variable. # ------------------------------------------------------------ if LBCODE == 1 or LBCODE == 2: # 1 = Unrotated regular lat/long grid # 2 = Regular lat/lon grid boxes (grid points are box # centres) ix = 11 iy = 10 elif LBCODE == 101 or LBCODE == 102: # 101 = Rotated regular lat/long grid # 102 = Rotated regular lat/lon grid boxes (grid points # are box centres) ix = -11 # rotated longitude (not an official axis code) iy = -10 # rotated latitude (not an official axis code) elif LBCODE >= 10000: # Cross section ix, iy = divmod(divmod(LBCODE, 10000)[1], 100) else: ix = None iy = None iz = _lbvc_to_axiscode.setdefault(LBVC, None) # Set it from the calendar type if iy in (20, 23) or ix in (20, 23): # Time is dealt with by x or y it = None elif calendar == 'gregorian': it = 20 else: it = 23 self.ix = ix self.iy = iy self.iz = iz self.it = it self.cf_info = {} # Set a identifying name based on the submodel and STASHcode # (or field code). stash = int_hdr[lbuser4] submodel = int_hdr[lbuser7] self.stash = stash # The STASH code has been set in the PP header, so try to find # its standard_name from the conversion table stash_records = _stash2standard_name.get((submodel, stash), None) um_Units = None long_name = None um_condition = None if stash_records: um_version = self.um_version for (long_name, units, valid_from, valid_to, standard_name, cf_info, um_condition) in stash_records: # Check that conditions are met if not self.test_um_version(valid_from, valid_to, um_version): continue if um_condition: if not self.test_um_condition(um_condition, LBCODE, BPLAT, BPLON): continue # Still here? Then we have our standard_name, etc. if standard_name: if set_standard_name: cf_properties['standard_name'] = standard_name else: attributes['_standard_name'] = standard_name cf_properties['long_name'] = long_name.rstrip() um_Units = _Units.get(units, None) if um_Units is None: um_Units = Units(units) _Units[units] = um_Units self.um_Units = um_Units self.cf_info = cf_info break #--- End: for #--- End: if if stash: section, item = divmod(stash, 1000) um_stash_source = 'm%02ds%02di%03d' % (submodel, section, item) cf_properties['um_stash_source'] = um_stash_source identity = 'UM_%s_vn%s' % (um_stash_source, self.um_version) else: identity = 'UM_%d_fc%d_vn%s' % (submodel, int_hdr[lbfc], self.um_version) if um_Units is None: self.um_Units = _Units[None] if um_condition: identity += '_%s' % um_condition if long_name is None: cf_properties['long_name'] = identity for recs, nz, nt in zip(groups, groups_nz, groups_nt): self.recs = recs self.nz = nz self.nt = nt self.z_recs = recs[:nz] self.t_recs = recs[::nz] LBUSER5 = recs[0].int_hdr.item(lbuser5,) self.cell_method_axis_name = {'area': 'area'} self.down_axes = set() self.z_axis = 'z' # ------------------------------------------------------------ # Get the extra data for this group # ------------------------------------------------------------ extra = recs[0].get_extra_data() self.extra = extra # ------------------------------------------------------------ # Set some derived metadata quantities # ------------------------------------------------------------ if self._debug: print self.__dict__ self.printfdr() # ------------------------------------------------------------ # Create the 'T' dimension coordinate # ------------------------------------------------------------ axiscode = it if axiscode is not None: c = self.time_coordinate(axiscode) # ------------------------------------------------------------ # Create the 'Z' dimension coordinate # ------------------------------------------------------------ axiscode = iz if axiscode is not None: # Get 'Z' coordinate from LBVC if axiscode == 3: c = self.atmosphere_hybrid_sigma_pressure_coordinate(axiscode) elif axiscode == 2 and 'height' in self.cf_info: # Create the height coordinate from the information # given in the STASH to standard_name conversion table height, units = self.cf_info['height'] c = self.size_1_height_coordinate(axiscode, height, units) elif axiscode == 14: c = self.atmosphere_hybrid_height_coordinate(axiscode) else: c = self.z_coordinate(axiscode) # Create a model_level_number auxiliary coordinate LBLEV = int_hdr[lblev] if LBVC in (2, 9, 65) or LBLEV in (7777, 8888): # CHECK! self.LBLEV = LBLEV c = self.model_level_number_coordinate(aux=bool(c)) #--- End: if # ------------------------------------------------------------ # Create the 'Y' dimension coordinate # ------------------------------------------------------------ axiscode = iy yc = None if axiscode is not None: if axiscode in (20, 23): # 'Y' axis is time-since-reference-date if extra.get('y', None) is not None: c = self.time_coordinate_from_extra_data(axiscode, 'y') else: LBUSER3 = int_hdr[lbuser3] if LBUSER3 == LBROW: self.lbuser3 = LBUSER3 c = self.time_coordinate_from_um_timeseries(axiscode, 'y') else: yc = self.xy_coordinate(axiscode, 'y') #--- End: if # ------------------------------------------------------------ # Create the 'X' dimension coordinate # ------------------------------------------------------------ axiscode = ix xc = None if axiscode is not None: if axiscode in (20, 23): # X axis is time since reference date if extra.get('x', None) is not None: c = self.time_coordinate_from_extra_data(axiscode, 'x') else: LBUSER3 = int_hdr[lbuser3] if LBUSER3 == LBNPT: self.lbuser3 = LBUSER3 c = self.time_coordinate_from_um_timeseries(axiscode, 'x') else: xc = self.xy_coordinate(axiscode, 'x') #--- End: if # -10: rotated latitude (not an official axis code) # -11: rotated longitude (not an official axis code) if (iy, ix) == (-10, -11) or (iy, ix) == (-11, -10): # ---------------------------------------------------- # Create a ROTATED_LATITUDE_LONGITUDE coordinate # reference # ---------------------------------------------------- transform = CoordinateReference( name='rotated_latitude_longitude', grid_north_pole_latitude=BPLAT, grid_north_pole_longitude=BPLON, coords=(_axis['y'], _axis['x'])) # -------------------------------------------------------- # Create UNROTATED, 2-D LATITUDE and LONGITUDE auxiliary # coordinates # -------------------------------------------------------- self.latitude_longitude_2d_aux_coordinates(yc, xc, transform) # Insert the coordinate reference into the domain self.domain.insert_ref(transform, copy=False) #--- End: if # ------------------------------------------------------------ # Create a RADIATION WAVELENGTH dimension coordinate # ------------------------------------------------------------ try: rwl, rwl_units = self.cf_info['below'] except (KeyError, TypeError): pass else: c = self.radiation_wavelength_coordinate(rwl, rwl_units) # Set LBUSER5 to zero so that it is not confused for a # pseudolevel LBUSER5 = 0 #--- End: try # ------------------------------------------------------------ # Create a PSEUDOLEVEL dimension coordinate. This must be # done *after* the possible creation of a radiation # wavelength dimension coordinate. # ------------------------------------------------------------ if LBUSER5 != 0: self.pseudolevel_coordinate(LBUSER5) attributes['int_hdr'] = int_hdr[:] attributes['real_hdr'] = real_hdr[:] attributes['file'] = filename attributes['id'] = identity cf_properties['Conventions'] = __Conventions__ cf_properties['runid'] = self.decode_lbexp() cf_properties['lbproc'] = str(LBPROC) cf_properties['lbtim'] = str(LBTIM) cf_properties['stash_code'] = str(stash) # ------------------------------------------------------------ # Create cell methods # ------------------------------------------------------------ cell_methods = self.create_cell_methods() if cell_methods is not None: cf_properties['cell_methods'] = cell_methods # ------------------------------------------------------------ # Set the data and extra data # ------------------------------------------------------------ data = self.create_data() cf_properties['_FillValue'] = data.fill_value # ------------------------------------------------------------ # Create the field # ------------------------------------------------------------ # Add kwargs to the CF properties cf_properties.update(kwargs) field = Field(domain=self.domain, data=self.data, axes=self.data_axes, properties=cf_properties, attributes=attributes, copy=False) # Check for decreasing axes that aren't decreasing down_axes = self.down_axes if down_axes: field.flip(down_axes, i=True) # Force cyclic X axis for paritcular values of LBHEM if int_hdr[lbhem] in (0, 1, 2, 4): field.cyclic('X', period=360) self.fields.append(field) #--- End: for self._nonzero = True #--- End: def def __nonzero__(self): ''' x.__nonzero__() <==> bool(x) ''' return self._nonzero #--- End: if def __repr__(self): ''' x.__repr__() <==> repr(x) ''' return self.fdr() #--- End: def def __str__(self): ''' x.__str__() <==> str(x) ''' out = [self.fdr()] attrs = ('endian', 'reftime', 'vtime', 'dtime', 'um_version', 'source', 'it', 'iz', 'ix', 'iy', 'site_time_cross_section', 'timeseries', 'file') for attr in attrs: out.append('%s=%s' % (attr, getattr(self, attr, None))) out.append('') return '\n'.join(out) #--- End: def def atmosphere_hybrid_height_coordinate(self, axiscode): ''' **From appendix A of UMDP F3** From UM Version 5.2, the method of defining the model levels in PP headers was revised. At vn5.0 and 5.1, eta values were used in the PP headers to specify the levels of model data, which was of limited use when plotting data on model levels. From 5.2, the PP headers were redefined to give information on the height of the level. Given a 2D orography field, the height field for a given level can then be derived. The height coordinates for PP-output are defined as: Z(i,j,k)=Zsea(k)+C(k)*orography(i,j) where Zsea(k) and C(k) are height based hybrid coefficients. Zsea(k) = eta_value(k)*Height_at_top_of_model C(k)=[1-eta_value(k)/eta_value(first_constant_rho_level)]**2 forlevels less than or equal to first_constant_rho_level C(k)=0.0 for levels greater than first_constant_rho_level where eta_value(k) is the eta_value for theta or rho level k. The eta_value is a terrain-following height coordinate; full details are given in UMDP15, Appendix B. The PP headers store Zsea and C as follows :- * 46 = bulev = brsvd1 = Zsea of upper layer boundary * 52 = blev = Zsea of level * 53 = brlev = Zsea of lower layer boundary * 47 = bhulev = brsvd2 = C of upper layer boundary * 54 = bhlev = C of level * 55 = bhrlev = C of lower layer boundary :Parameters: axiscode : int :Returns: out : cf.DimensionCoordinate or None ''' domain = self.domain zdim = _axis['z'] # Insert new Z axis domain.insert_axis(self.nz, key=zdim) # "a" auxiliary coordinate array = numpy_array([rec.real_hdr[blev] for rec in self.z_recs], # Zsea dtype=float) bounds0 = numpy_array([rec.real_hdr[brlev] for rec in self.z_recs], #Zsea lower dtype=float) bounds1 = numpy_array([rec.real_hdr[brsvd1] for rec in self.z_recs], #Zsea upper dtype=float) bounds = numpy_column_stack((bounds0, bounds1)) ac = AuxiliaryCoordinate() ac = self.coord_data(ac, array, bounds, units=_Units['m']) ac.id = 'UM_atmosphere_hybrid_height_coordinate_a' ac.long_name = 'height based hybrid coeffient a' key_a = domain.insert_aux(ac, axes=[zdim], copy=False) # atmosphere_hybrid_height_coordinate dimension coordinate TOA_height = bounds1.max() if TOA_height <= 0: TOA_height = self.height_at_top_of_model if not TOA_height: dc = None else: array = array / TOA_height bounds = bounds / TOA_height dc = DimensionCoordinate() dc = self.coord_data(dc, array, bounds, units=_Units['']) dc.standard_name = 'atmosphere_hybrid_height_coordinate' dc = self.coord_axis(dc, axiscode) dc = self.coord_positive(dc, axiscode, zdim) domain.insert_dim(dc, key=zdim, copy=False) #--- End: if # "b" auxiliary coordinate array = numpy_array([rec.real_hdr[bhlev] for rec in self.z_recs], dtype=float) bounds0 = numpy_array([rec.real_hdr[bhrlev] for rec in self.z_recs], dtype=float) bounds1 = numpy_array([rec.real_hdr[brsvd2] for rec in self.z_recs], dtype=float) bounds = numpy_column_stack((bounds0, bounds1)) ac = AuxiliaryCoordinate() ac = self.coord_data(ac, array, bounds, units=_Units['1']) ac.id = 'UM_atmosphere_hybrid_height_coordinate_b' ac.long_name = 'height based hybrid coeffient b' key_b = domain.insert_aux(ac, axes=[zdim], copy=False) if bool(dc): self.cell_method_axis_name['z'] = dc.identity() # atmosphere_hybrid_height_coordinate coordinate reference ref = CoordinateReference( name='atmosphere_hybrid_height_coordinate', a=key_a, b=key_b, coords=(key_a, key_b), coord_terms=('a', 'b')) self.domain.insert_ref(ref, copy=False) #--- End: if return dc #--- End: def def depth_coordinate(self, axiscode): ''' :Parameters: axiscode : int :Returns: out : cf.DimensionCoordinate or None ''' dc = self.model_level_number_coordinate(aux=False) domain = self.domain zdim = _axis['z'] array = numpy_array([rec.real_hdr[blev] for rec in self.z_recs], dtype=float) bounds0 = numpy_array([rec.real_hdr[brlev] for rec in self.z_recs], dtype=float) bounds1 = numpy_array([rec.real_hdr[brsvd1] for rec in self.z_recs], dtype=float) bounds = numpy_column_stack((bounds0, bounds1)) ac = AuxiliaryCoordinate() ac = self.coord_data(ac, array, bounds, units=_Units['m']) ac.id = 'UM_atmosphere_hybrid_height_coordinate_ak' ac.long_name = 'atmosphere_hybrid_height_coordinate_ak' domain.insert_aux(ac, axes=[zdim], copy=False) array = numpy_array([rec.real_hdr[bhlev] for rec in self.z_recs], dtype=float) bounds0 = numpy_array([rec.real_hdr[bhrlev] for rec in self.z_recs], dtype=float) bounds1 = numpy_array([rec.real_hdr[brsvd2] for rec in self.z_recs], dtype=float) bounds = numpy_column_stack((bounds0, bounds1)) ac = AuxiliaryCoordinate() ac = self.coord_data(ac, array, bounds, units=_Units['1']) ac.id = 'UM_atmosphere_hybrid_height_coordinate_bk' ac.long_name = 'atmosphere_hybrid_height_coordinate_bk' domain.insert_aux(ac, axes=[zdim], copy=False) if dc: self.cell_method_axis_name['z'] = dc.identity() return dc #--- End: def def atmosphere_hybrid_sigma_pressure_coordinate(self, axiscode): ''' atmosphere_hybrid_sigma_pressure_coordinate when not an array axis :Parameters: axiscode : int :Returns: out : cf.DimensionCoordinate ''' # 46 BULEV Upper layer boundary or BRSVD(1) # # 47 BHULEV Upper layer boundary or BRSVD(2) # # For hybrid levels: # - BULEV is B-value at half-level above. # - BHULEV is A-value at half-level above. # # For hybrid height levels (vn5.2-, Smooth heights) # - BULEV is Zsea of upper layer boundary # * If rho level: Zsea for theta level above # * If theta level: Zsea for rho level above # - BHLEV is C of upper layer boundary # * If rho level: C for theta level above # * If theta level: C for rho level above array = [] bounds = [] ak_array = [] ak_bounds = [] bk_array = [] bk_bounds = [] for rec in self.z_recs: BLEV, BRLEV, BHLEV, BHRLEV, BULEV, BHULEV = self.header_bz(rec) array.append(BLEV + BHLEV/_pstar) bounds.append([BRLEV + BHRLEV/_pstar, BULEV + BHULEV/_pstar]) ak_array.append(BHLEV) ak_bounds.append((BHRLEV, BHULEV)) bk_array.append(BLEV) bk_bounds.append((BRLEV , BULEV)) #--- End: for array = numpy_array(array , dtype=float) bounds = numpy_array(bounds , dtype=float) ak_array = numpy_array(ak_array , dtype=float) ak_bounds = numpy_array(ak_bounds, dtype=float) bk_array = numpy_array(bk_array , dtype=float) bk_bounds = numpy_array(bk_bounds, dtype=float) domain = self.domain zdim = _axis['z'] dc = DimensionCoordinate() dc = self.coord_data( dc, array, bounds, units=_axiscode_to_Units.setdefault(axiscode, None)) dc = self.coord_positive(dc, axiscode, zdim) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) domain.insert_dim(dc, key=zdim, copy=False) ac = AuxiliaryCoordinate() ac = self.coord_data(ac, ak_array, ak_bounds, units=_Units['Pa']) ac.id = 'UM_atmosphere_hybrid_sigma_pressure_coordinate_ak' ac.long_name = 'atmosphere_hybrid_sigma_pressure_coordinate_ak' domain.insert_aux(ac, axes=[zdim], copy=False) ac = AuxiliaryCoordinate() ac = self.coord_data(ac, bk_array, bk_bounds, units=_Units['1']) domain.insert_aux(ac, axes=[zdim], copy=False) ac.id = 'UM_atmosphere_hybrid_sigma_pressure_coordinate_bk' ac.long_name = 'atmosphere_hybrid_sigma_pressure_coordinate_bk' self.cell_method_axis_name['z'] = dc.identity() return dc #--- End: def def create_cell_methods(self): '''Create the cell methods ''' cell_methods = [] LBPROC = self.lbproc LBTIM_IB = self.lbtim_ib tmean_proc = 0 if LBTIM_IB in (2, 3) and LBPROC in (128, 192, 2176, 4224, 8320): tmean_proc = 128 LBPROC -= 128 # ------------------------------------------------------------ # Area cell methods # ------------------------------------------------------------ # -10: rotated latitude (not an official axis code) # -11: rotated longitude (not an official axis code) if self.ix in (10, 11, 12, -10, -11) and self.iy in (10, 11, 12, -10, -11): cf_info = self.cf_info if 'where' in cf_info: cell_methods.append('area: mean') cell_methods.append(cf_info['where']) if 'over' in cf_info: cell_methods.append(cf_info['over']) #--- End: if if LBPROC == 64: cell_methods.append('x: mean') # dch : do special zonal mean as as in pp_cfwrite # ------------------------------------------------------------ # Vertical cell methods # ------------------------------------------------------------ if LBPROC == 2048: cell_methods.append('z: mean') # ------------------------------------------------------------ # Time cell methods # ------------------------------------------------------------ if LBTIM_IB == 0 or LBTIM_IB == 1: cell_methods.append('t: point') elif LBPROC == 4096: cell_methods.append('t: minimum') elif LBPROC == 8192: cell_methods.append('t: maximum') if tmean_proc == 128: if LBTIM_IB == 2: cell_methods.append('t: mean') elif LBTIM_IB == 3: cell_methods.append('t: mean within years') cell_methods.append('t: mean over years') #--- End: if if not cell_methods: return None cell_methods = CellMethods(' '.join(cell_methods)) cell_method_axis_name = self.cell_method_axis_name for c in cell_methods: names0 = c.names[0] c.axes = [_axis[name] for name in names0] c.names = [cell_method_axis_name[name] for name in names0] #--- End: for return cell_methods #--- End: def def coord_axis(self, c, axiscode): axis = _coord_axis.setdefault(axiscode, None) if axis is not None: c.axis = axis return c #--- End: def def coord_data(self, c, array=None, bounds=None, units=None, fill_value=None, climatology=False): ''' Set the data array of a coordinate construct. :Parameters: c : cf.DimensionCoordinate or cf.AuxiliaryCoordinate data : array-like, optional The data array. bounds : array-like, optional The Cell bounds for the data array. units : cf.Units, optional The units of the data array. fill_value : optional climatology : bool, optional Whether or not the coordinate construct is a time climatology. By default it is not. :Returns: out : cf.Coordinate ''' if array is not None: array = Data(array, units=units, fill_value=fill_value) if bounds is not None: bounds = Data(bounds, units=units, fill_value=fill_value) if climatology: c.climatology = True #--- End: if c.insert_data(array, bounds=bounds, copy=False) return c #--- End: def def coord_names(self, c, axiscode): ''' ''' standard_name = _coord_standard_name.setdefault(axiscode, None) if standard_name is not None: c.setprop('standard_name', standard_name) c.ncvar = standard_name else: long_name = _coord_long_name.setdefault(axiscode, None) if long_name is not None: c.long_name = long_name return c #--- End: def def coord_positive(self, c, axiscode, dim): positive = _coord_positive.setdefault(axiscode, None) if positive is not None: c.positive = positive if positive == 'down' and axiscode != 4: self.down_axes.add(dim) #--- End: def return c #--- End: def def ctime(self, rec): ''' ''' reftime = self.refUnits LBVTIME = tuple(self.header_vtime(rec)) LBDTIME = tuple(self.header_dtime(rec)) key = (LBVTIME, LBDTIME, self.refunits, self.calendar) ctime = _cached_ctime.get(key, None) if ctime is None: # LTIME = list(LBDTIME) # LTIME[0] = LBVTIME[0] ctime = Datetime(*LBDTIME) ctime.year = LBVTIME[0] if ctime < Datetime(*LBVTIME): ctime.year += 1 ctime = Data(ctime, reftime).array.item() _cached_ctime[key] = ctime #--- End: if return ctime #--- End: def def header_vtime(self, rec): ''' Return the list [LBYR, LBMON, LBDAT, LBHR, LBMIN] for the given record. :Parameters: rec : :Returns: out : list :Examples: >>> u.header_vtime(rec) [1991, 1, 1, 0, 0] ''' return rec.int_hdr[lbyr:lbmin+1] #--- End: def def header_dtime(self, rec): ''' Return the list [LBYRD, LBMOND, LBDATD, LBHRD, LBMIND] for the given record. :Parameters: rec : :Returns: out : list :Examples: >>> u.header_dtime(rec) [1991, 2, 1, 0, 0] ''' return rec.int_hdr[lbyrd:lbmind+1] #--- End: def def header_bz(self, rec): ''' Return the list [BLEV, BRLEV, BHLEV, BHRLEV, BULEV, BHULEV] for the given record. :Parameters: rec : :Returns: out : list :Examples: >>> u.header_bz(rec) ''' real_hdr = rec.real_hdr return (real_hdr[blev:bhrlev+1].tolist() + # BLEV, BRLEV, BHLEV, BHRLEV real_hdr[brsvd1:brsvd2+1].tolist()) # BULEV, BHULEV #--- End: def def header_lz(self, rec): ''' Return the list [LBLEV, LBUSER5] for the given record. :Parameters: rec : :Returns: out : list :Examples: >>> u.header_lz(rec) ''' int_hdr = rec.int_hdr return [int_hdr.item(lblev,), int_hdr.item(lbuser5,)] #--- End: def def header_z(self, rec): ''' Return the list [LBLEV, LBUSER5, BLEV, BRLEV, BHLEV, BHRLEV, BULEV, BHULEV] for the given record. :Parameters: rec : :Returns: out : list :Examples: >>> u.header_z(rec) ''' # ------------------------------------------------------------ # These header items are used by the compare_levels function # in compare.c # ------------------------------------------------------------ return self.header_lz + self.header_bz #--- End: def def create_data(self): ''' Sets the `!data` and `!data_axes` attributes. :Returns: None ''' LBROW = self.lbrow LBNPT = self.lbnpt yx_shape = (LBROW, LBNPT) yx_size = LBROW * LBNPT nz = self.nz nt = self.nt recs = self.recs units = self.um_Units data_type_in_file = self.data_type_in_file filename = self.filename data_axes = [_axis['y'], _axis['x']] if len(recs) == 1: # -------------------------------------------------------- # 0-d partition matrix # -------------------------------------------------------- rec = recs[0] data = Data(UMFileArray(file=filename, ndim=2, shape=yx_shape, size=yx_size, dtype=data_type_in_file(rec), header_offset=rec.hdr_offset, data_offset=rec.data_offset, disk_length=rec.disk_length), units=units, fill_value=rec.real_hdr.item(bmdi,)) else: # -------------------------------------------------------- # 1-d or 2-d partition matrix # -------------------------------------------------------- file_data_types = set() word_sizes = set() # Find the partition matrix shape pmshape = [n for n in (nt, nz) if n > 1] pmndim = len(pmshape) partitions = [] empty_list = [] partitions_append = partitions.append zero_to_LBROW = (0, LBROW) zero_to_LBNPT = (0, LBNPT) if pmndim == 1: # ---------------------------------------------------- # 1-d partition matrix # ---------------------------------------------------- data_ndim = 3 if nz > 1: pmaxes = [_axis[self.z_axis]] data_shape = (nz, LBROW, LBNPT) data_size = nz * yx_size else: pmaxes = [_axis['t']] data_shape = (nt, LBROW, LBNPT) data_size = nt * yx_size partition_shape = [1, LBROW, LBNPT] for i, rec in enumerate(recs): # Find the data type of the array in the file file_data_type = data_type_in_file(rec) file_data_types.add(file_data_type) subarray = UMFileArray(file=filename, ndim=2, shape=yx_shape, size=yx_size, dtype=file_data_type, header_offset=rec.hdr_offset, data_offset=rec.data_offset, disk_length=rec.disk_length) partitions_append(Partition( subarray = subarray, location = [(i, i+1), zero_to_LBROW, zero_to_LBNPT], shape = partition_shape, axes = data_axes, flip = empty_list, part = empty_list, Units = units)) #--- End: for # Populate the 1-d partition matrix matrix = numpy_array(partitions, dtype=object) else: # ---------------------------------------------------- # 2-d partition matrix # ---------------------------------------------------- pmaxes = [_axis['t'], _axis[self.z_axis]] data_shape = (nt, nz, LBROW, LBNPT) data_size = nt * nz * yx_size data_ndim = 4 partition_shape = [1, 1, LBROW, LBNPT] for i, rec in enumerate(recs): # Find T and Z axis indices t, z = divmod(i, nz) # Find the data type of the array in the file file_data_type = data_type_in_file(rec) file_data_types.add(file_data_type) subarray = UMFileArray(file=filename, ndim=2, shape=yx_shape, size=yx_size, dtype=file_data_type, header_offset=rec.hdr_offset, data_offset=rec.data_offset, disk_length=rec.disk_length) partitions_append(Partition( subarray=subarray, location=[(t, t+1), (z, z+1), zero_to_LBROW, zero_to_LBNPT], shape=partition_shape, axes=data_axes, flip=empty_list, part=empty_list, Units=units)) #--- End: for # Populate the 2-d partition matrix matrix = numpy_array(partitions, dtype=object) matrix.resize(pmshape) #--- End: if data_axes = pmaxes + data_axes # Set the data array data = Data(units=units, fill_value=recs[0].real_hdr.item(bmdi,)) data._axes = data_axes data._shape = data_shape data._ndim = data_ndim data._size = data_size data.partitions = PartitionMatrix(matrix, pmaxes) data.dtype = numpy_result_type(*file_data_types) #--- End: if self.data = data self.data_axes = data_axes return data #---End: def def decode_lbexp(self): '''Decode the integer value of LBEXP in the PP header into a runid. If this value has already been decoded, then it will be returned from the cache, otherwise the value will be decoded and then added to the cache. :Returns: out : str A string derived from LBEXP. If LBEXP is a negative integer then that number is returned as a string. :Examples: >>> self.decode_lbexp() 'aaa5u' >>> self.decode_lbexp() '-34' ''' LBEXP = self.int_hdr[lbexp] runid = _cached_runid.get(LBEXP, None) if runid is not None: # Return a cached decoding of this LBEXP return runid if LBEXP < 0: runid = str(LBEXP) else: # Convert LBEXP to a binary string, filled out to 30 bits with # zeros bits = bin(LBEXP) bits = bits.lstrip('0b').zfill(30) # Step through 6 bits at a time, converting each 6 bit chunk into # a decimal integer, which is used as an index to the characters # lookup list. runid = [] for i in xrange(0,30,6): index = int(bits[i:i+6], 2) if index < _n_characters: runid.append(_characters[index]) #--- End: for runid = ''.join(runid) #--- End: def # Enter this runid into the cache _cached_runid[LBEXP] = runid # Return the runid return runid #--- End: def def dtime(self, rec): ''' ''' reftime = self.refUnits units = self.refunits calendar = self.calendar LBDTIME = tuple(self.header_dtime(rec)) key = (LBDTIME, units, calendar) time = _cached_date2num.get(key, None) if time is None: # It is important to use the same time_units as vtime if self.calendar == 'gregorian': time = netCDF4_date2num( datetime(*LBDTIME), units, calendar) else: time = netCDF4_date2num( netCDF4_netcdftime_datetime(*LBDTIME), units, calendar) _cached_date2num[key] = time #--- End: if return time #--- End: def def fdr(self): '''Return a the contents of PP field headers as strings. This is a bit like printfdr in the UKMO IDL PP library. :Returns: out : list ''' out2 = [] for i, rec in enumerate(self.recs): out = ['Field %d:' % i] x = ['%s::%s' % (name, value) for name, value in zip(_header_names, self.int_hdr + self.real_hdr)] x = textwrap.fill(' '.join(x), width=79) out.append(x.replace('::', ': ')) if self.extra: out.append('EXTRA DATA:') for key in sorted(self.extra): out.append('%s: %s' % (key, str(self.extra[key]))) #--- End: if out.append('file: '+self.filename) out.append('format, byte order, word size: %s, %s, %d' % (self.fmt, self.byte_ordering, self.word_size)) out.append('') out2.append('\n'.join(out)) #--- End: for return out2 #--- End: def def latitude_longitude_2d_aux_coordinates(self, yc, xc, transform): ''' ''' BDX = self.bdx BDY = self.bdy LBNPT = self.lbnpt LBROW = self.lbrow BPLAT = self.bplat BPLON = self.bplon # Create the unrotated latitude and longitude arrays if we # couldn't find them in the cache cache_key = (LBNPT, LBROW, BDX, BDY, BPLAT, BPLON) lat, lon = _cache_latlon.get(cache_key, (None, None)) if lat is None: lat, lon = self.unrotated_latlon(yc.varray, xc.varray, BPLAT, BPLON) atol = self.atol if abs(BDX) >= atol and abs(BDY) >= atol: _cache_latlon[cache_key] = (lat, lon) #--- End: if # if xc.hasbounds and yc.hasbounds: # cache_key = ('bounds',) + cache_key # lat_bounds, lon_bounds = _cache_latlon.get(cache_key, (None, None)) # print lat_bounds # if lat_bounds is None: # print 'CALC BOUNDS' # xb = numpy_empty(xc.size+1) # yb = numpy_empty(yc.size+1) # xb[:-1] = xc.bounds.subspace[ :, 0].squeeze(1, i=True).array # xb[-1 ] = xc.bounds.subspace[-1, 1].squeeze(1, i=True).array # yb[:-1] = yc.bounds.subspace[ :, 0].squeeze(1, i=True).array # yb[-1 ] = yc.bounds.subspace[-1, 1].squeeze(1, i=True).array # # lat_bounds, lon_bounds = self.unrotated_latlon(yb, xb, BPLAT, BPLON) # # print lat_bounds # print lat_bounds.shape # yyy = numpy_empty(lat.shape + (4,)) # # # # print lat.shape, yyy.shape # # atol = self.atol # if abs(BDX) >= atol and abs(BDY) >= atol: # _cache_latlon[cache_key] = (lat, lon) #--- End: if axes = [_axis['y'], _axis['x']] for axiscode, array in zip((10, 11), (lat, lon)): ac = AuxiliaryCoordinate() ac = self.coord_data(ac, array, units=_axiscode_to_Units.setdefault(axiscode, None)) ac = self.coord_names(ac, axiscode) key = self.domain.insert_aux(ac, axes=axes, copy=False) transform.coords.add(key) #--- End: for #--- End: def def model_level_number_coordinate(self, aux=False): '''model_level_number dimension or auxiliary coordinate :Parameters: aux : bool :Returns: out : cf.AuxiliaryCoordinate or cf.DimensionCoordinate or None ''' array = tuple([rec.int_hdr.item(lblev,) for rec in self.z_recs]) key = array c = _cached_model_level_number_coordinate.get(key, None) if c is not None: if aux: self.domain.insert_aux(c, axes=[_axis['z']], copy=True) else: self.domain.insert_dim(c, key=_axis['z'], copy=True) self.cell_method_axis_name['z'] = c.identity() else: array = numpy_array(array, dtype=self.int_hdr_dtype) if array.min() < 0: return array = numpy_where(array==9999, 0, array) axiscode = 5 if aux: c = AuxiliaryCoordinate() c = self.coord_data(c, array, units=Units('1')) c = self.coord_names(c, axiscode) self.domain.insert_aux(c, axes=[_axis['z']], copy=False) else: c = DimensionCoordinate() c = self.coord_data(c, array, units=Units('1')) c = self.coord_names(c, axiscode) c = self.coord_axis(c, axiscode) self.domain.insert_dim(c, key=_axis['z'], copy=False) self.cell_method_axis_name['z'] = c.identity() #--- End: if _cached_model_level_number_coordinate[key] = c #--- End: if return c #--- End: def def data_type_in_file(self, rec): '''Return the data type of the data array. :Parameters: rec : umfile.Rec :Returns: out : numpy.dtype :Examples: ''' # Find the data type if rec.int_hdr.item(lbuser2,) == 3: # Boolean return numpy_dtype(bool) else: # Int or float return rec.get_type_and_num_words()[0] # rec_file = rec.file ## data_type = rec_file.c_interface.get_type_and_length( # data_type = rec_file.c_interface.get_type_and_num_words(rec.int_hdr)[0] # if data_type == 'int': # # Integer # data_type = 'int%d' % (rec_file.word_size * 8) # else: # # Float # data_type = 'float%d' % (rec_file.word_size * 8) # #--- End: if # # return numpy_dtype(data_type) #--- End: def def printfdr(self): '''Print out the contents of PP field headers. This is a bit like printfdr in the UKMO IDL PP library. :Examples: >>> u.printfdr() ''' for header in self.fdr(): print header #--- End: def def pseudolevel_coordinate(self, LBUSER5): ''' ''' if self.nz == 1: array = numpy_array((LBUSER5,), dtype=self.int_hdr_dtype) else: # 'Z' aggregation has been done along the pseudolevel axis array = numpy_array([rec.int_hdr.item(lbuser5,) for rec in self.z_recs], dtype=self.int_hdr_dtype) self.z_axis = 'p' #--- End: if axiscode = 40 dc = DimensionCoordinate() dc = self.coord_data( dc, array, units=_axiscode_to_Units.setdefault(axiscode, None)) dc.long_name = 'pseudolevel' # for PP stash_code %d' % self.stash dc.id = 'UM_pseudolevel' self.domain.insert_dim(dc, key=_axis['p'], copy=False) self.cell_method_axis_name['p'] = dc.identity() return dc #--- End: def def radiation_wavelength_coordinate(self, rwl, rwl_units): ''' ''' array = numpy_array((rwl,), dtype=float) bounds = numpy_array(((0.0, rwl)), dtype=float) units = _Units.get(rwl_units, None) if units is None: units = Units(rwl_units) _Units[rwl_units] = units axiscode = -20 dc = DimensionCoordinate() dc = self.coord_data(dc, array, bounds, units=units) dc = self.coords_names(dc, axiscode) self.domain.insert_dim(dc, key=_axis['r'], copy=False) self.cell_method_axis_name['r'] = dc.identity() return dc #--- End: def def reference_time_Units(self): ''' ''' time_units = 'days since %d-1-1' % self.int_hdr[lbyr] calendar = self.calendar key = time_units+' calendar='+calendar units = _Units.get(key, None) if units is None: units = Units(time_units, calendar) _Units[key] = units #--- End: if self.refUnits = units self.refunits = time_units return units #--- End: def def size_1_height_coordinate(self, axiscode, height, units): # Create the height coordinate from the information given in the # STASH to standard_name conversion table key = (axiscode, height, units) dc = _cached_size_1_height_coordinate.get(key, None) zdim = _axis['z'] if dc is not None: copy = True else: height_units = _Units.get(units, None) if height_units is None: height_units = Units(units) _Units[units] = height_units array = numpy_array((height,), dtype=float) dc = DimensionCoordinate() dc = self.coord_data(dc, array, units=height_units) dc = self.coord_positive(dc, axiscode, zdim) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) _cached_size_1_height_coordinate[key] = dc copy = False #--- End: def self.domain.insert_dim(dc, key=zdim, copy=copy) self.cell_method_axis_name['z'] = dc.identity() return dc #--- End: def def test_um_condition(self, um_condition, LBCODE, BPLAT, BPLON): '''Return True if a field satisfies the condition specified for a STASH code to standard name conversion. :Parameters: um_condition : str LBCODE : int BPLAT : float BPLON : float :Returns: out : bool True if a field satisfies the condition specified, False otherwise. :Examples: >>> ok = u.test_um_condition('true_latitude_longitude', ...) ''' if um_condition == 'true_latitude_longitude': if LBCODE in _true_latitude_longitude_lbcodes: return True # Check pole location in case of incorrect LBCODE atol = self.atol if (abs(BPLAT-90.0) <= atol + RTOL()*90.0 and abs(BPLON) <= atol): return True elif um_condition == 'rotated_latitude_longitude': if LBCODE in _rotated_latitude_longitude_lbcodes: return True # Check pole location in case of incorrect LBCODE atol = self.atol if not (abs(BPLAT-90.0) <= atol + RTOL()*90.0 and abs(BPLON) <= atol): return True else: raise ValueError( "Unknown UM condition in STASH code conversion table: '%s'" % um_condition) # Still here? Then the condition has not been satisfied. return #--- End: def def test_um_version(self, valid_from, valid_to, um_version): '''Return True if the UM version applicable to tghis field is within the given range. If possible, the UM version is derived from the PP header and stored in the metadata object. Otherwise it is taken from the *um_version* parameter. :Parameters: valid_from : int, float or None valid_to : int, float or None um_version : int or float :Returns: out : bool True if the UM version applicable to this field is within the range, False otherwise. :Examples: >>> ok = u.test_um_version(401, 505, 1001) >>> ok = u.test_um_version(401, None, 606.3) >>> ok = u.test_um_version(None, 405, 401) ''' if valid_to is None: if valid_from <= um_version: return True elif valid_from is None: if um_version <= valid_to: return True elif valid_from <= um_version <= valid_to: return True return False # if valid_from is '': # valid_from = None # # if valid_to is '': # if valid_from <= um_version: # return True # elif valid_from <= um_version <= valid_to: # return True # # return False #--- End: def def time_coordinate(self, axiscode): ''' Return the T dimension coordinate ''' recs = self.t_recs vtimes = numpy_array([self.vtime(rec) for rec in recs], dtype=float) dtimes = numpy_array([self.dtime(rec) for rec in recs], dtype=float) IB = self.lbtim_ib if IB <= 1 or vtimes.item(0,) >= dtimes.item(0,): array = vtimes bounds = None climatology = False elif IB == 3: # The field is a time mean from T1 to T2 for each year # from LBYR to LBYRD ctimes = numpy_array([self.ctime(rec) for rec in recs]) array = 0.5*(vtimes + ctimes) bounds = numpy_column_stack((vtimes, dtimes)) climatology = True else: array = 0.5*(vtimes + dtimes) bounds = numpy_column_stack((vtimes, dtimes)) climatology = False #--- End: if dc = DimensionCoordinate() dc = self.coord_data(dc, array, bounds, units=self.refUnits, climatology=climatology) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) self.domain.insert_dim(dc, key=_axis['t'], copy=False) self.cell_method_axis_name['t'] = dc.identity() return dc #--- End: def def time_coordinate_from_extra_data(self, axiscode, axis): ''' ''' extra = self.extra array = extra[axis] bounds = extra.get(axis+'_bounds', None) calendar = self.calendar if calendar == '360_day': units = _Units['360_day 0-1-1'] elif calendar == 'gregorian': units = _Units['gregorian 1752-9-13'] elif calendar == '365_day': units = _Units['365_day 1752-9-13'] else: units = None dc = DimensionCoordinate() dc = self.coord_data(dc, array, bounds, units=units) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) self.domain.insert_dim(dc, key=_axis[axis], copy=False) self.cell_method_axis_name[axis] = dc.identity() self.cell_method_axis_name['t'] = self.cell_method_axis_name[axis] return dc #--- End: def def time_coordinate_from_um_timeseries(self, axiscode, axis): # This PP/FF field is a timeseries. The validity time is # taken to be the time for the first sample, the data time # for the last sample, with the others evenly between. rec = self.recs[0] vtime = self.vtime(rec) dtime = self.dtime(rec) size = self.lbuser3 - 1.0 delta = (dtime - vtime)/size array = numpy_arange(vtime, vtime+delta*size, size, dtype=float) dc = DimensionCoordinate() dc = self.coord_data(dc, array, units=units) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) self.domain.insert_dim(dc, key=_axis[axis], copy=False) self.cell_method_axis_name['t'] = dc.identity() return dc #--- End: def def vtime(self, rec): ''' ''' reftime = self.refUnits units = self.refunits calendar = self.calendar # LBVTIME = tuple(rec.int_hdr[lbyr: lbmin+1]) LBVTIME = tuple(self.header_vtime(rec)) key = (LBVTIME, units, calendar) time = _cached_date2num.get(key, None) if time is None: # It is important to use the same time_units as dtime if self.calendar == 'gregorian': time = netCDF4_date2num( datetime(*LBVTIME), units, calendar) else: time = netCDF4_date2num( netCDF4_netcdftime_datetime(*LBVTIME), units, calendar) _cached_date2num[key] = time #--- End: if return time #--- End: def def dddd(self): for axis_code, extra_type in zip((11 , 10 ), ('x', 'y')): coord_type = extra_type + '_domain_bounds' if coord_type in p.extra: p.extra[coord_type] # Create, from extra data, an auxiliary coordinate should # with 1) data and bounds, if the upper and lower be # bounds have no missing values; or 2) data but no the # bounds, if the upper bound has missing values axis # but the lower bound does not. # which file_position = ppfile.tell() # has bounds = p.extra[coord_type][...] # axis_code # Reset the file pointer after reading the extra # 13 # data into a numpy array ppfile.seek(file_position, os.SEEK_SET) data = None if numpy_any(bounds[..., 1] == _pp_rmdi): # dch also test in bmdi? if not numpy_any(bounds[...,0] == _pp_rmdi): # dch also test in bmdi? data = bounds[...,0] bounds = None else: data = numpy_mean(bounds, axis=1) if (data, bounds) != (None, None): aux = 'aux%(auxN)d' % locals() auxN += 1 # Increment auxiliary number coord = _create_Coordinate(domain, aux, axis_code, p=p, array = data, aux=True, bounds_array = bounds, pubattr = {'axis': None}, dimensions = [xdim]) # DCH # xdim? # should # be # the # axis # which # has # axis_code # 13 #--- End: if else: coord_type = '%s_domain_lower_bound' % extra_type if coord_type in p.extra: # Create, from extra data, an auxiliary # coordinate with data but no bounds, if the # data noes not contain any missing values file_position = ppfile.tell() data = p.extra[coord_type][...] # Reset the file pointer after reading the # extra data into a numpy array ppfile.seek(file_position, os.SEEK_SET) if not numpy_any(data == _pp_rmdi): # dch also test in bmdi? aux = 'aux%(auxN)d' % locals() auxN += 1 # Increment auxiliary number coord = _create_Coordinate(domain, aux, axis_code, p=p, aux=True, array=numpy_array(data), pubattr={'axis': None}, dimensions=[xdim])# DCH xdim? #--- End: if #--- End: for #--- End: if # -------------------- def unrotated_latlon(self, rotated_lat, rotated_lon, pole_lat, pole_lon): ''' Create 2-d arrays of unrotated latitudes and longitudes. :Parameters: rotated_lat, rotated_lon, pole_lat, pole_lon ''' # Make sure rotated_lon and pole_lon is in [0, 360) pole_lon = pole_lon % 360.0 # Convert everything to radians pole_lon *= _pi_over_180 pole_lat *= _pi_over_180 cos_pole_lat = numpy_cos(pole_lat) sin_pole_lat = numpy_sin(pole_lat) # Create appropriate copies of the input rotated arrays rot_lon = rotated_lon.copy() rot_lat = rotated_lat.view() # Make sure rotated longitudes are between -180 and 180 rot_lon %= 360.0 rot_lon = numpy_where(rot_lon < 180.0, rot_lon, rot_lon-360) # Create 2-d arrays of rotated latitudes and longitudes in radians nlat = rot_lat.size nlon = rot_lon.size rot_lon = numpy_resize(numpy_deg2rad(rot_lon), (nlat, nlon)) rot_lat = numpy_resize(numpy_deg2rad(rot_lat), (nlon, nlat)) rot_lat = numpy_transpose(rot_lat, axes=(1,0)) # Find unrotated latitudes CPART = numpy_cos(rot_lon) * numpy_cos(rot_lat) sin_rot_lat = numpy_sin(rot_lat) x = cos_pole_lat * CPART + sin_pole_lat * sin_rot_lat x = numpy_clip(x, -1.0, 1.0) unrotated_lat = numpy_arcsin(x) # Find unrotated longitudes x = -cos_pole_lat*sin_rot_lat + sin_pole_lat*CPART x /= numpy_cos(unrotated_lat) # dch /0 or overflow here? surely # lat could be ~+-pi/2? if so, # does x ~ cos(lat)? x = numpy_clip(x, -1.0, 1.0) unrotated_lon = -numpy_arccos(x) unrotated_lon = numpy_where(rot_lon > 0.0, -unrotated_lon, unrotated_lon) if pole_lon >= self.atol: SOCK = pole_lon - numpy_pi else: SOCK = 0 unrotated_lon += SOCK # Convert unrotated latitudes and longitudes to degrees unrotated_lat = numpy_rad2deg(unrotated_lat) unrotated_lon = numpy_rad2deg(unrotated_lon) # Return unrotated latitudes and longitudes return unrotated_lat, unrotated_lon #--- End: def def xy_coordinate(self, axiscode, axis): ''' Create an X or Y dimension coordinate from header entries or extra data. :Parameters: axiscode : int axis : str 'x' or 'y' :Returns: out : cf.DimensionCoordinate ''' if axis == 'y': delta = self.bdy origin = self.real_hdr[bzy] size = self.lbrow else: delta = self.bdx origin = self.real_hdr[bzx] size = self.lbnpt if abs(delta) > self.atol: # Create regular coordinates from header items if axiscode == 11 or axiscode == -11: origin -= divmod(origin + delta*size, 360.0)[0] * 360 while origin + delta*size > 360.0: origin -= 360.0 while origin + delta*size < -360.0: origin += 360.0 #--- End: if array = numpy_arange(origin+delta, origin+delta*(size+0.5), delta, dtype=float) # Create the coordinate bounds if axiscode in (13, 31, 40, 99): # The following axiscodes do not have bounds: # 13 = Site number (set of parallel rows or columns # e.g.Time series) # 31 = Logarithm to base 10 of pressure in mb # 40 = Pseudolevel # 99 = Other bounds = None else: delta_by_2 = 0.5 * delta bounds = numpy_empty((size, 2), dtype=float) bounds[:, 0] = array - delta_by_2 bounds[:, 1] = array + delta_by_2 else: # Create coordinate from extra data array = self.extra.get(axis, None) bounds = self.extra.get(axis+'_bounds', None) #--- End: if dc = DimensionCoordinate() dc = self.coord_data( dc, array, bounds, units=_axiscode_to_Units.setdefault(axiscode, None)) dc = self.coord_positive(dc, axiscode, _axis[axis]) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) self.domain.insert_dim(dc, key=_axis[axis], copy=False) self.cell_method_axis_name[axis] = dc.identity() return dc #--- End: def def z_coordinate(self, axiscode): '''Create a Z dimension coordinate from BLEV :Parameters: axiscode : int :Returns: out : cf.DimensionCoordinate ''' z_recs = self.z_recs array = tuple([rec.real_hdr.item(blev,) for rec in z_recs]) bounds0 = tuple([rec.real_hdr[brlev] for rec in z_recs]) # lower level boundary bounds1 = tuple([rec.real_hdr[brsvd1] for rec in z_recs]) # bulev if _coord_positive.get(axiscode, None) == 'down': bounds0, bounds1 = bounds1, bounds0 # key = (axiscode, array, bounds0, bounds1) # dc = _cached_z_coordinate.get(key, None) # if dc is not None: # copy = True # else: copy = False array = numpy_array(array, dtype=float) bounds0 = numpy_array(bounds0, dtype=float) bounds1 = numpy_array(bounds1, dtype=float) bounds = numpy_column_stack((bounds0, bounds1)) if (bounds0 == bounds1).all(): bounds = None else: bounds = numpy_column_stack((bounds0, bounds1)) dc = DimensionCoordinate() dc = self.coord_data( dc, array, bounds=bounds, units=_axiscode_to_Units.setdefault(axiscode, None)) dc = self.coord_positive(dc, axiscode, _axis['z']) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) # _cached_z_coordinate[key] = dc # #--- End: if self.domain.insert_dim(dc, key=_axis['z'], copy=copy) self.cell_method_axis_name['z'] = dc.identity() return dc #--- End: def def z_reference_coordinate(self, axiscode): ''' ''' array = numpy_array([rec.real_hdr.item(brlev,) for rec in self.z_recs], dtype=float) LBVC = self.lbvc key = (axiscode, LBVC, array) dc = _cached_z_reference_coordinate.get(key, None) if dc is not None: copy = True else: if not 128 <= LBVC <= 139: bounds = [] for rec in self.z_recs: BRLEV = rec.real_hdr.item(brlev,) BRSVD1 = rec.real_hdr.item(brsvd1,) if abs(BRSVD1-BRLEV) >= ATOL: bounds = None break bounds.append((BRLEV, BRSVD1)) #--- End: for else: bounds = None if bounds: bounds = numpy_array((bounds,), dtype=float) dc = DimensionCoordinate() dc = self.coord_data( dc, array, bounds, units=_axiscode_to_Units.setdefault(axiscode, None)) dc = self.coord_axis(dc, axiscode) dc = self.coord_names(dc, axiscode) if not dc.get('positive', True): # ppp dc.flip(i=True) _cached_z_reference_coordinate[key] = dc copy = False #--- End: def self.domain.insert_dim(dc, key=_axis['z'], copy=copy) return dc #--- End: def #--- End: class _stash2standard_name = {} def load_stash2standard_name(table=None, delimiter='!'): '''Load a STASH to standard name conversion table. :Parameters: table : str, optional Use the conversion table at this file location. By default the table will be looked for at ``os.path.join(os.path.dirname(cf.__file__),'etc/STASH_to_CF.txt')`` delimiter : str, optional The delimiter of the table columns. By default, ``!`` is taken as the delimiter. :Returns: None *Examples:* >>> load_stash2standard_name() >>> load_stash2standard_name('my_table.txt') >>> load_stash2standard_name('my_table2.txt', ',') ''' # 0 Model # 1 STASH code # 2 STASH name # 3 units # 4 valid from UM vn # 5 valid to UM vn # 6 standard_name # 7 CF extra info # 8 PP extra info if table is None: # Use default conversion table package_path = os.path.dirname(__file__) table = os.path.join(package_path, 'etc/STASH_to_CF.txt') #--- End: if lines = csv.reader(open(table, 'r'), delimiter=delimiter, skipinitialspace=True) raw_list = [] [raw_list.append(line) for line in lines] # Get rid of comments for line in raw_list[:]: if line[0].startswith('#'): raw_list.pop(0) continue break #--- End: for # Convert to a dictionary which is keyed by (submodel, STASHcode) # tuples (model, stash, name, units, valid_from, valid_to, standard_name, cf, pp) = range(9) stash2sn = {} for x in raw_list: key = (int(x[model]), int(x[stash])) if not x[units]: x[units] = None try: cf_info = {} if x[cf]: for d in x[7].split(): if d.startswith('height='): cf_info['height'] = re.split(_number_regex, d, re.IGNORECASE)[1:4:2] if cf_info['height'] == '': cf_info['height'][1] = '1' if d.startswith('below_'): cf_info['below'] = re.split(_number_regex, d, re.IGNORECASE)[1:4:2] if cf_info['below'] == '': cf_info['below'][1] = '1' if d.startswith('where_'): cf_info['where'] = d.replace('where_', 'where ', 1) if d.startswith('over_'): cf_info['over'] = d.replace('over_', 'over ', 1) x[cf] = cf_info except IndexError: pass try: x[valid_from] = float(x[valid_from]) except ValueError: x[valid_from] = None try: x[valid_to] = float(x[valid_to]) except ValueError: x[valid_to] = None x[pp] = x[pp].rstrip() line = (x[name:],) if key in stash2sn: stash2sn[key] += line else: stash2sn[key] = line #--- End: for _stash2standard_name.clear() _stash2standard_name.update(stash2sn) # return stash2sn #--- End: def # --------------------------------------------------------------------- # Create the STASH code to standard_name conversion dictionary # --------------------------------------------------------------------- #_stash2standard_name = load_stash2standard_name() load_stash2standard_name() def read(filename, um_version=405, verbose=False, aggregate=True, byte_ordering=None, word_size=None, set_standard_name=True, height_at_top_of_model=None): '''Read fields from a PP file or UM fields file. The file may be big or little endian, 32 or 64 bit :Parameters: filename : file or str A string giving the file name, or an open file object, from which to read fields. um_version : number, optional The Unified Model (UM) version to be used when decoding the PP header. Valid versions are, for example, ``402`` (v4.2), ``606.3`` (v6.6.3) and ``1001`` (v10.1). The default version is ``405`` (v4.5). The version is ignored if it can be inferred from the PP headers, which will generally be the case for files created at versions 5.3 and later. Note that the PP header can not encode tertiary version elements (such as the ``3`` in ``606.3``), so it may be necessary to provide a UM version in such cases. verbose : bool, optional set_standard_name : bool, optional :Returns: out : FieldList The fields in the file. :Examples: >>> f = read('file.pp') >>> f = read('*/file[0-9].pp', um_version=708) ''' history = 'Converted from UM by cf-python v%s' % __version__ f = _open_um_file(filename) # print 'um.read.read.py um_version=', repr(um_version) um = [UMField(var, f.format, f.byte_ordering, f.word_size, um_version, set_standard_name, history=history, height_at_top_of_model=height_at_top_of_model) for var in f.vars] # # Clear the cache of unrotated latitude and longitude arrays # _cache_latlon.clear() return FieldList([field for x in um for field in x.fields if field]) #--- End: def def _atmosphere_hybrid_sigma_pressure_coordinate(f): # atmosphere_hybrid_sigma_pressure_coordinate real_hdr = f.real_hdr BLEV, BRLEV, BHLEV, BHRLEV = real_hdr[blev:bhrlev+1] BRSVD1, BRSVD2 = real_hdr[brsvd1:brsvd2+1] if 'z' not in pmaxes: indices = [0] * pmndim indices[pmaxes.index('z')] = slice(1, None, None) for rec in var.recs[tuple(indices)]: BLEV, BRLEV, BHLEV, BHRLEV = rec.real_hdr[blev:bhrlev+1] BRSVD1, BRSVD2 = rec.real_hdr[brsvd1:brsvd2+1] array.append(BLEV + BHLEV/_pstar) bounds.append([BRLEV + BHRLEV/_pstar, BRSVD1 + BRSVD2/_pstar]) ak_array.append(BHLEV) ak_bounds.append([BHRLEV, BRSVD2]) bk_array.append(BLEV) bk_bounds.append([BRLEV , BRSVD1]) array = numpy_array(array , dtype=float) bounds = numpy_array(bounds , dtype=float) ak_array = numpy_array(ak_array , dtype=float) ak_bounds = numpy_array(ak_bounds, dtype=float) bk_array = numpy_array(bk_array , dtype=float) bk_bounds = numpy_array(bk_bounds, dtype=float) domain = f.field.domain coord = _create_Coordinate(domain, vdim, axis_code=axis_code, p=p, array=array, bounds_array=bounds, dimensions=[vdim]) coord = _create_Coordinate( domain, aux, axis_code=None, p=p, pubattr={'standard_name': 'atmosphere_hybrid_sigma_pressure_coordinate_ak'}, units=_units['Pa'], array=ak_array, bounds_array=ak_array, dimensions=[vdim], aux=True) coord = _create_Coordinate( domain, aux, axis_code=None, p=p, pubattr={'standard_name': 'atmosphere_hybrid_sigma_pressure_coordinate_bk'}, units=_units['1'], array=bk_array, bounds_array=bk_array, dimensions=[vdim], aux=True) else: cache_key = ('atmosphere_hybrid_sigma_pressure_coordinate', BLEV, BHLEV, BRLEV, BHRLEV, BRSVD1, BRSVD2) cache_key_ak = ('atmosphere_hybrid_sigma_pressure_coordinate_ak', BLEV, BHLEV, BRLEV, BHRLEV, BRSVD1, BRSVD2) cache_key_bk = ('atmosphere_hybrid_sigma_pressure_coordinate_bk', BLEV, BHLEV, BRLEV, BHRLEV, BRSVD1, BRSVD2) if cache_key in _cached_coordinate: c = _cached_coordinate[cache_key] domain.insert_dim(c, key=vdim, copy=True) c = _cached_coordinate[cache_key_ak] domain.insert_aux(c, axes=[vdim], copy=True) c = _cached_coordinate[cache_key_bk] domain.insert_aux(c, axes=[vdim], copy=True) else: array = numpy_array((BLEV + BHLEV/_pstar,), dtype=float) bounds = numpy_array(((BRLEV + BHRLEV/_pstar, BRSVD1 + BRSVD2/_pstar)), dtype=float) ak_array = numpy_array((BHLEV,) , dtype=float) ak_bounds = numpy_array(((BHRLEV, BRSVD2)), dtype=float) bk_array = numpy_array((BLEV,) , dtype=float) bk_bounds = numpy_array(((BRLEV, BRSVD1)), dtype=float) coord = _create_Coordinate(domain, vdim, axis_code=axis_code, p=p, array=array, bounds_array=bounds, dimensions=[vdim], cache_key=cache_key) coord = _create_Coordinate( domain, aux, axis_code=None, p=p, pubattr={'standard_name': 'atmosphere_hybrid_sigma_pressure_coordinate_ak'}, units=_units['Pa'], array=ak_array, bounds_array=ak_array, dimensions=[vdim], aux=True, cache_key=cache_key_ak) coord = _create_Coordinate( domain, aux, axis_code=None, p=p, pubattr={'standard_name': 'atmosphere_hybrid_sigma_pressure_coordinate_bk'}, units=_units['1'], array=bk_array, bounds_array=bk_array, dimensions=[vdim], aux=True, cache_key=cache_key_bk) #--- End: def def is_um_file(filename): '''Return True if a file is a PP file or UM fields file. Note that the file type is determined by inspecting the file's contents and any file suffix is not not considered. :Parameters: filename : str :Returns: out : bool :Examples: >>> is_um_file('myfile.pp') True >>> is_um_file('myfile.nc') False >>> is_um_file('myfile.pdf') False >>> is_um_file('myfile.txt') False ''' try: f = _open_um_file(filename) except: return False try: f.close_fd() except: pass return True #--- End: def ''' Problems: Z and P coordinates /home/david/data/pp/aaaao/aaaaoa.pmh8dec.03328.pp /net/jasmin/chestnut/data-24/david/testpp/026000000000c.fc0607.000128.0000.00.04.0260.0020.1491.12.01.00.00.pp skipping variable stash code=0, 0, 0 because: grid code not supported umfile: error condition detected in routine list_copy_to_ptr_array umfile: error condition detected in routine process_vars umfile: error condition detected in routine file_parse OK 2015-04-01 /net/jasmin/chestnut/data-24/david/testpp/026000000000c.fc0619.000128.0000.00.04.0260.0020.1491.12.01.00.00.pp skipping variable stash code=0, 0, 0 because: grid code not supported umfile: error condition detected in routine list_copy_to_ptr_array umfile: error condition detected in routine process_vars umfile: error condition detected in routine file_parse OK 2015-04-01 /net/jasmin/chestnut/data-24/david/testpp/lbcode_10423.pp skipping variable stash code=0, 0, 0 because: grid code not supported umfile: error condition detected in routine list_copy_to_ptr_array umfile: error condition detected in routine process_vars umfile: error condition detected in routine file_parse OK 2015-04-01 /net/jasmin/chestnut/data-24/david/testpp/lbcode_11323.pp skipping variable stash code=0, 0, 0 because: grid code not supported umfile: error condition detected in routine list_copy_to_ptr_array umfile: error condition detected in routine process_vars umfile: error condition detected in routine file_parse OK 2015-04-01 EXTRA_DATA: /net/jasmin/chestnut/data-24/david/testpp/ajnjgo.pmm1feb.pp SLOW: (Not any more! 2015-04-01) /net/jasmin/chestnut/data-24/david/testpp/xgdria.pdk949a.pp /net/jasmin/chestnut/data-24/david/testpp/xhbmaa.pm27sep.pp RUN LENGTH ENCODED dump (not fields file) /home/david/data/um/xhlska.dak69h0 Field 115 (stash code 9) dch@eslogin008:/nerc/n02/n02/dch> ff2pp xgvwko.piw96b0 xgvwko.piw96b0.pp file xgvwko.piw96b0 is a byte swapped 64 bit ieee um file '''