""" This function accepts the fields needed by WRF and more specifically: based on the GFS Vtable. Assumptions: - all the 3D fields need to be on isobaric levels when they are passed to the encoding function - every 3D field must be supplied with a corresponding surface field - the COARDS ordering of dimensions is assumed i.e. time, level, lat, lon Usage: """ ############################################################################# # TODO # - change the 'data' structure to a dictionary creating appropriate labels # from the dates ############################################################################# import os import sys # the syntax to import nio depends on its version. We try all options from the # newest to the oldest try: import Nio as nio except: try: import PyNGL.Nio as nio except: import PyNGL_numpy.Nio as nio # The following is only needed for hn3.its.monash.edu.au # Unfortunately we must get rid of spurious entries in the sys.path that can # lead to the loading of conflicting packages for dir_idx in range(len(sys.path)-1,-1,-1): if 'python2.4' in sys.path[dir_idx]: sys.path.pop(dir_idx) import numpy as n import pylab as p from matplotlib.numerix import ma import grib2 from string import zfill import pywrf.viz.utils as vu from pywrf.misc.met_utils import * #data_directory = '/Users/val/data/laps_nc_files/press/first_try/' #sfc_data_directory = '/Users/val/data/laps_surface_data/' #data_directory = '/mnt/mac/Users/val/data/laps_nc_files/press/first_try/' #sfc_data_directory = '/mnt/mac/Users/val/data/laps_surface_data/' #data_directory = '/media/DATA/wrf/canberra_nc_files/press/' #sfc_data_directory = '/media/DATA/wrf/canberra_nc_files/sfc/' #data_directory = '/media/DATA/wrf/canberra/press/' #sfc_data_directory = '/media/DATA/wrf/canberra/sfc/' #data_directory = '/Volumes/COMMON/wrf/canberra/press/' #sfc_data_directory = '/Volumes/COMMON/wrf/canberra/sfc/' #data_directory = '/Volumes/DATA/wrf/canberra/press/' #sfc_data_directory = '/Volumes/DATA/wrf/canberra/sfc/' #data_directory = '/Users/val/Desktop/workspace/plotting/canberra/press/' #sfc_data_directory = '/Users/val/Desktop/workspace/plotting/canberra/sfc/' #data_directory = \ # '/nfs/1/home/vbisign/wrf/wps_data/pygrib2_test/canberra/press/analyses/' #sfc_data_directory = \ # '/nfs/1/home/vbisign/wrf/wps_data/pygrib2_test/canberra/sfc/' data_directory = \ '/nfs/1/home/vbisign/wrf/wps_data/run_003/press/' sfc_data_directory = \ #sfc_data_directory = '/Users/val/Desktop/workspace/plotting/canberra/sfc/' # It is assumed that most people will use a single nc file as the source # of their fields. Neverthless, the field nc_file_name is available for those # that need to use multiple data sources. # the code is split in two parts... the first part gathers all the # information that is specific to your data source... stuff like the location # of the netcdf files, the name of the variables of interest in your files, # and any special treatment of the fields like changing units, scaling the # data or anything that is needed to get the data in the form required by the # encoding function. # This translates to simply having to generate a fairly pedantic dictionary # that spells out what needs to be done with the fields during the grib # encoding # the second part is the encoding function itself in which it is assumed the # data is passed in a dictionary containing the data as numpy arrays in the # right units # Unfortunately it proves hard to both generalize the data structure and # maintain flexibility... hence at the price of making the code longer I am # going to process one field at the time so that it is clear to anyone what # the steps that are (may be) involved are... # I still believe it is a good idea to collect all the information needed in # one big dictionary for ease of retrieval ############################################################################# # Preliminaries # This is the data structure that I will pass to the encoder data = [] # to generate a grib2 file using Jeffrey Whitaker's grib2 class, it is necessary # to provide the constructor with two pieces of information: # - the discipline of the data to be encoded # - the identification section # this means the class will need to be reconstructed every time the discipline # changes def missing(octects): # To say that something is not available in a grib field they use a # sequence of 1s as long as the space allocated (number of octects). # When this sequence is read as an integer, it corresponds to the maximum # representable number given the octects*8 bits of memory, hence: # PS: refer to the grib documentation to know how many octects are # allocated to any particular field return 2**(8*octects) - 1 def prepare_sfc_data(f, fv): # work out which file to pick based on the base date of the first file # this is based on the (fragile) assumption that nobody changed the file # names from the standard laps... # this was needed for the first batch of data... valid_date = fv['valid_date'].get_value()[0] valid_time = fv['valid_time'].get_value()[0] # the following is old stuff if 0: # the data in each file 'starts' at 12:00pm (UTC) if valid_time < 1200: valid_date -= 1 #sfc_file = 'laps' + str(valid_date)[4:] + '_surface.nc' print sfc_data_directory + sfc_file # this is new and (hopefully) better if 1: # there is one sfc data file for each of the analysis containing the # full 73 hours of diagnostic fields # in the simplified case that we do not span analysis -> we start again # with a different file at every analysis: base_date = fv['base_date'].get_value() base_time = fv['base_time'].get_value() sfc_file = 'regprog.laps_pt375.' \ + str(base_date)[4:] + '.' + str(base_time).zfill(4)[:2] \ + '.slvfld.n3.nc' # if we use data from a single forecast we just name the appropriate file # sfc_file = 'laps0117_surface.nc' sfc_f = nio.open_file(sfc_data_directory + sfc_file) sfc_fv = sfc_f.variables # I need this to work out which time to pick in the file... sfc_valid_time = sfc_fv['valid_time'].get_value() time_idx = sfc_valid_time.tolist().index(valid_time) return sfc_file, sfc_f, sfc_fv, time_idx ############################################################################# ############################################################################# # this variables remain unchanged across the different times # BTW these (with the exception of the significance_ref_time) should not # affect WRF # Id of originating centre (Common Code Table C-1) -> 1 for Melbourne # Id of originating centre (Common Code Table C-1) -> 7 for NCEP centre = 7 # Id of orginating sub-centre (local table) -> 0 for no subcentre # Id of orginating sub-centre (local table) -> 3 for NCEP Central Operations subcentre = 0 # GRIB Master Tables Version Number (Code Table 1.0) master_table = 1 # GRIB Local Tables Version Number (Code Table 1.1) -> 0 for not used # GRIB Local Tables Version Number (Code Table 1.1) -> 1 for let's see what happens local_table = 1 # Significance of Reference Time (Code Table 1.2) -> 1 for start of forecast significance_ref_time = 1 # Production status of data (Code Table 1.3) -> 2 for research product production_status = 2 # idsect[12]=Type of processed data (Code Table 1.4) -> 1 -> forecast product type_processed_data = 1 ############################################################################# # the first dimension of data will be the times to be processed # you will have to modify this to suit your data # this could mean just having a time index to get different times from the # same file or just pick different files for different times. files = os.listdir(data_directory) #for file in files[:2]: for file in files: f = nio.open_file(data_directory + file) fv = f.variables base_date = str(fv['base_date'].get_value()) year = int(base_date[0:4]) month = int(base_date[4:6]) day = int(base_date[6:8]) # padding needed to use a string method on an integer base_time = zfill(str(fv['base_time'].get_value()),4) hour = int(base_time[0:2]) minute = int(base_time[2:]) second = 0 # I did not generate this in the preamble above because I needed to wait # for the base date and time info from the netcdf files... # NB the order of the variables in the list MATTERS! idsect = [ centre, subcentre, master_table, local_table, significance_ref_time, year, month, day, hour, minute, second, production_status, type_processed_data ] print file # gdsinfo - Sequence containing information needed for the grid definition # section. gdsinfo = n.zeros((5,), dtype=int) # gdsinfo[0] = Source of grid definition (see Code Table 3.0) -> 0 -> will be # specified in octects 13-14 gdsinfo[0] = 0 # gdsinfo[1] = Number of grid points in the defined grid. # TODO make it more general #lat = fv['lat'].get_value() # the idx was chosen to stop the data at 55S and avoid the rubbish far south lat_cheat_idx = 30 # 1 for the full range lon_cheat_idx = 1 lat = fv['lat'].get_value()[lat_cheat_idx:] lon = fv['lon'].get_value()[:-lon_cheat_idx] gdsinfo[1] = len(lat)*len(lon) # gdsinfo[2] = Number of octets needed for each additional grid points defn. # Used to define number of points in each row for non-reg grids # (=0 for regular grid). gdsinfo[2] = 0 # gdsinfo[3] = Interp. of list for optional points defn (Code Table 3.11) # 0 -> There is no appended list gdsinfo[3] = 0 # gdsinfo[4] = Grid Definition Template Number (Code Table 3.1) # 0 -> Latitude/Longitude (or equidistant cylindrical, or Plate Carree) gdsinfo[4] = 0 # gdtmpl - Contains the data values for the specified Grid Definition # Template ( NN=gds[4] ). Each element of this integer array contains # an entry (in the order specified) of Grid Definition Template 3.NN gdtmpl = n.zeros((19,), dtype=n.int64) # In the following we refer to the lat/lon template (grid template 3.0) # Oct: 15 Shape of the Earth (See Table 3.2) # 0 -> spherical with radius = 6,367,470.0 m # TODO double check the radius they used in Tan's code gdtmpl[0] = 0 # Oct: 16 Scale Factor of radius of spherical Earth # TODO hopefully (given the last option) the next 5 are ignored... # if 1 octet is the basic unit then a value of 255 should be equivalent # to all bits set to 1 -> this could be broken if more than one octect # is assigned -> maybe I should use the maximum represantable number # tailored to the number of octects? # are there any rules that make the encoder ignore automatically certain # options e.g. the oblate spheroid one when a spherical shape is assumed gdtmpl[1] = missing(1) # Oct: 17-20 Scale value of radius of spherical Earth # NB they (grib2 std) interprets all bits set to 1 as a missing value # hence the value is base*(word_width*words_assigned) - 1 # base is 2 since we are working in binary # word_width is 8 since they use octects # words_assigned is the number of words reserved to represent the value # -1 gives us the greatest representable number since we start from 0 gdtmpl[2] = missing(1) # Oct: 21 Scale factor of major axis of oblate spheroid Earth gdtmpl[3] = missing(1) # Oct: 22-25 Scaled value of major axis of oblate spheroid Earth gdtmpl[4] = missing(4) # Oct: 26 Scale factor of minor axis of oblate spheroid Earth gdtmpl[5] = missing(1) # Oct: 27-30 Scaled value of minor axis of oblate spheroid Earth gdtmpl[6] = missing(4) # Oct: 31-34 Number of points along a parallel gdtmpl[7] = len(lon) # Oct: 35-38 Number of points along a meridian gdtmpl[8] = len(lat) # Oct: 39-42 Basic angle of the initial production domain (see Note 1) gdtmpl[9] = missing(4) # Oct: 43-46 Subdivisions of basic angle used to define extreme longitudes # and latitudes, and direction increments (see Note 1) gdtmpl[10] = missing(4) # Oct: 47-50 La1 - Latitude of first grid point (see Note 1) #gdtmpl.append(in_vars['lat'].get_value()[-1]*1e6) gdtmpl[11] = lat[-1]*1e6 #gdtmpl[11] = lat[0]*1e6 # Oct: 51-54 Lo1 - Longitude of first grid point (see Note 1) #gdtmpl.append(in_vars['lon'].get_value()[0]*1e6) gdtmpl[12] = lon[0]*1e6 # Oct: 55 Resolution and component flags (see Table 3.3) # This should set all the bits to 0 meaning # bit value meaning # 1-2 ? Reserved # 3 0 i direction increments not given # 4 0 j direction increments not given # 5 0 Resolved u and v components of vector quantities # relative to easterly and northerly directions # 6-8 0 Reserved - set to zero #gdtmpl[13] = 0 #gdtmpl[13] = 12 gdtmpl[13] = 48 # Oct: 56-59 La2 - Latitude of last grid point (see Note1) #gdtmpl.append(in_vars['lat'].get_value()[0]*1e6) gdtmpl[14] = lat[0]*1e6 #gdtmpl[14] = lat[-1]*1e6 # Oct: 60-63 Lo2 - Longitude of last grid point (see Note 1) #gdtmpl.append(in_vars['lon'].get_value()[-1]*1e6) gdtmpl[15] = lon[-1]*1e6 # Oct: 64-67 i (west-east) Direction increment (see Note1) #gdtmpl[16] = missing(4) gdtmpl[16] = 0.375*1e6 # Oct: 68-71 j (south-north) Direction increment (see Note1) #gdtmpl[17] = missing(4) gdtmpl[17] = 0.375*1e6 #gdtmpl[17] = -0.375*1e6 #gdtmpl[17] = 33143 # Oct: 72 Scanning mode (flags see Table 3.4) # TODO double check if this is legit considering the flags # This should set all the bits to 0 meaning # bit value meaning # 1 0 Points in the first row or column scan in +i (+x) # direction # 2 0 Points in the first row or column scan in -i (-x) # direction # 3 0 Adjacent pints in i(x) direction are consecutive # 4 0 All rows scan in the same direction # 5-8 0 Reserved (set to zero) # NB I suspect the 3 bit might be important for the 'majoring' of the array? # # 2^7 2^6 2^5 2^4 2^3 2^2 2^1 2^0 # # 1 2 3 4 5 6 7 8 gdtmpl[18] = 0 # 0 0 0 0 0 0 0 0 #gdtmpl[18] = 16 # 0 0 0 1 0 0 0 0 #gdtmpl[18] = 32 # 0 0 1 0 0 0 0 0 #gdtmpl[18] = 64 # 0 1 0 0 0 0 0 0 #gdtmpl[18] = 128 # 1 0 0 0 0 0 0 0 #gdtmpl[18] = 192 # 1 1 0 0 0 0 0 0 #gdtmpl[18] = 160 # 1 0 1 0 0 0 0 0 #gdtmpl[18] = 144 # 1 0 0 1 0 0 0 0 #gdtmpl[18] = 96 # 0 1 1 0 0 0 0 0 #gdtmpl[18] = 80 # 0 1 0 1 0 0 0 0 #gdtmpl[18] = 48 # 0 0 1 1 0 0 0 0 #gdtmpl[18] = 224 # 1 1 1 0 0 0 0 0 #gdtmpl[18] = 208 # 1 1 0 1 0 0 0 0 #gdtmpl[18] = 176 # 1 0 1 1 0 0 0 0 #gdtmpl[18] = 112 # 0 1 1 1 0 0 0 0 #gdtmpl[18] = 240 # 1 1 1 1 0 0 0 0 data.append( { 'idsect' : idsect, 'gdsinfo' : gdsinfo, 'gdtmpl' : gdtmpl, 'discipline': { # the first number in the list is the code for the discipline 'atmos': { 'code':0, 'fields':{'2d':{}, '3d':{}} }, 'ocean': { 'code':10, 'fields':{'2d':{}, '3d':{}} }, 'land': { 'code':2, 'fields':{'2d':{}, '3d':{}} } } } ) # this is the part where we extract the data from the nc files and put # into the structure we will pass later to the encoder # let's start with the atmospheric products specifically the surface # fields # each variables comes with all the metadata necessary to do its grib # encoding using the addfield method of the encoder class i.e. pdtnum, # pdtmpl, drtnum, drtmpl, field # -1 stands for the last one (time) added dummy = data[-1]['discipline']['atmos']['fields']['2d'] # Let's prepare the surface pressure data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 3 -> mass pdtmpl[0] = 3 # Oct 11: Parameter number (see Code table 4.2). 0 -> pressure pdtmpl[1] = 0 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> ground or water surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # drtnum - Data Representation Template Number (see Code Table 5.0). # 0 -> grid point data - simple packing drtnum = 0 # drtmpl - Sequence with the data values for the specified Data Representation # Template (N=drtnum). Each element of this integer array contains an # entry (in the order specified) of Data Representation Template 5.N Note # that some values in this template (eg. reference values, number of # bits, etc...) may be changed by the data packing algorithms. Use this to # specify scaling factors and order of spatial differencing, if desired. drtmpl = n.zeros((5,),dtype=int) # NB to make sense of the following remember (from the wmo std): # The original data value Y (in the units of code table 4.2) can be recovered # with the formula: # Y * 10^D= R + (X1+X2) * 2^E # For simple packing and all spectral data # E = Binary scale factor, # D = Decimal scale factor # R = Reference value of the whole field, # X1 = 0, # X2 = Scaled (encoded) value. # Oct: 12-15. Reference value (R) (IEEE 32-bit floating-point value) drtmpl[0] = 0 # Oct: 16-17. Binary scale factor (E) drtmpl[1] = 0 # Oct: 18-19. Decimal scale factor (D) drtmpl[2] = 0 # Oct: 20. Number of bits used for each packed value for simple packing, or drtmpl[3] = 24 # for each group reference value for complex packing or spatial differencing # Oct: 21. Type of original field values (see Code Table 5.1) # 0 -> floating point drtmpl[4] = 0 cheat = fv['sfc_pres'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] if 0: dummy['sfc_pres'] = { 'long_name' : 'surface pressure', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['sfc_pres'].get_value()[0] #'field' : fv['sfc_pres'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } # Let's prepare the mean sea level pressure data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 3 -> mass pdtmpl[0] = 3 # Oct 11: Parameter number (see Code table 4.2). 1 -> MSLP pdtmpl[1] = 1 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 101 -> sea level pdtmpl[9] = 101 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['mslp'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] # to go from mb to Pa #cheat *= 100. cheat = cheat * 100. field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['mslp'] = { 'long_name' : 'mean sea level pressure', # in my data it comes in mb and I want in Pascals 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['mslp'].get_value()[0]*100. #'field' : fv['mslp'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx]*100. 'field' : field } # Let's prepare the 2m temperature data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> temperature pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 0 -> temperature pdtmpl[1] = 0 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> ground or water surface pdtmpl[9] = 103 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> 2m above ground pdtmpl[11] = 2 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['sfc_temp'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['sfc_temp'] = { 'long_name' : 'surface temperature', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['sfc_temp'].get_value()[0] #'field' : fv['sfc_temp'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } # Let's prepare the 2m temperature data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> temperature pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 0 -> temperature pdtmpl[1] = 0 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> 2m above ground pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['skin_temp'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['skin_temp'] = { 'long_name' : 'skin temperature', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['skin_temp'].get_value()[0] #'field' : fv['skin_temp'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } # the next fields will come from a different file... sfc_file, sfc_f, sfc_fv, time_idx = prepare_sfc_data(f,fv) # DEBUG #p.figure() #p.contour(dummy['skin_temp']['field']) #p.figure() #p.contour(sfc_fv['skin_temp'].get_value()[time_idx]) # Let's prepare the 2m RH data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 1 -> moisture pdtmpl[0] = 1 # Oct 11: Parameter number (see Code table 4.2). 1 -> RH pdtmpl[1] = 1 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 103 -> fixed height above ground pdtmpl[9] = 103 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> 2m above ground pdtmpl[11] = 2 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. # The surface files do not come with mixing ratio, but we can calculate # that from the dew point temperature... #dew_point = sfc_fv['scrn_dewpt'].get_value()[time_idx] dew_point = sfc_fv['scrn_dewpt'].get_value()[time_idx,lat_cheat_idx:,:-lon_cheat_idx] vapour_pressure = goff_gratch(dew_point) #sfc_temp = sfc_fv['scrn_temp'].get_value()[time_idx] sfc_temp = sfc_fv['scrn_temp'].get_value()[time_idx,lat_cheat_idx:,:-lon_cheat_idx] sat_vapour_pressure = goff_gratch(sfc_temp) #sat_vapour_pressure = ma.masked_where(sat_vapour_pressure == 0., # sat_vapour_pressure) sfc_rh = vapour_pressure / sat_vapour_pressure * 100 del sfc_temp, vapour_pressure, dew_point, sat_vapour_pressure cheat = sfc_rh field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['sfc_rh'] = { 'long_name' : '2m relative humidity', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : sfc_rh.filled(fill_value=0.) #'field' : sfc_rh 'field' : field } # Let's prepare the zonal 10m wind data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 2 -> momentum pdtmpl[0] = 2 # Oct 11: Parameter number (see Code table 4.2). 2 -> u-component of wind pdtmpl[1] = 2 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 103 -> fixed height above ground pdtmpl[9] = 103 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> 10m above ground pdtmpl[11] = 10 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = sfc_fv['q10m_zonal_wnd'].get_value()[time_idx,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['sfc_zonal_wnd'] = { 'long_name' : '10m zonal wind', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : sfc_fv['q10m_zonal_wnd'].get_value()[time_idx] #'field' : sfc_fv['q10m_zonal_wnd'].get_value()[time_idx,lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } # Let's prepare the meridional 10m wind data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 2 -> momentum pdtmpl[0] = 2 # Oct 11: Parameter number (see Code table 4.2). 3 -> v-component of wind pdtmpl[1] = 3 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 103 -> fixed height above ground pdtmpl[9] = 103 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> 10m above ground pdtmpl[11] = 10 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = sfc_fv['q10m_merid_wnd'].get_value()[time_idx,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['sfc_merid_wnd'] = { 'long_name' : '10m meridional wind', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : sfc_fv['q10m_merid_wnd'].get_value()[time_idx] #'field' : sfc_fv['q10m_merid_wnd'].get_value()[time_idx,lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } if 1: # Let's prepare the topography data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 3 -> mass pdtmpl[0] = 3 # Oct 11: Parameter number (see Code table 4.2). 5 -> geopotential height pdtmpl[1] = 5 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['topog'].get_value()[lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['topog'] = { 'long_name' : 'topography', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['topog'].get_value() #'field' : fv['topog'].get_value()[lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } if 1: # Let's prepare the water equivalent accumulated snow depth # data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 1 -> moisture pdtmpl[0] = 1 # Oct 11: Parameter number (see Code table 4.2). 13 -> # water equivalent of accumulated snow depth pdtmpl[1] = 13 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. #dummy_field = n.zeros((len(lon),len(lat)),dtype=float) dummy_field = n.zeros((len(lat),len(lon)),dtype=float) dummy['weasd'] = { 'long_name' : 'water equivalent of accumulated snow depth', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, 'field' : dummy_field } ####################################################################### # Finished with the 2d fields # Let's move on to the 3d ones ####################################################################### dummy = data[-1]['discipline']['atmos']['fields']['3d'] # Let's prepare the air temperature data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> temperature pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 0 -> temperature pdtmpl[1] = 0 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> isobaric level pdtmpl[9] = 100 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['air_temp'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['air_temp'] = { 'long_name' : '3D air temperature', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['air_temp'].get_value()[0], #'field' : fv['air_temp'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx], 'field' : field, #'field' : '', 'lvl' : fv['lvl'].get_value() } # Let's prepare the 3d relative humidity data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 1 -> moisture pdtmpl[0] = 1 # Oct 11: Parameter number (see Code table 4.2). 1 -> RH pdtmpl[1] = 1 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> isobaric level pdtmpl[9] = 100 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. #mix_rto = fv['mix_rto'].get_value()[0] mix_rto = fv['mix_rto'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] spec_hum = n.zeros(mix_rto.shape,dtype=float) spec_hum = mix_rto / (mix_rto + 1) #temp = fv['air_temp'].get_value()[0] temp = fv['air_temp'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] #print temp.min(), temp.max() #print temp.min()-273.16, temp.max()-273.16 sat_vapour_pres = goff_gratch(temp) #print sat_vapour_pres.min(), sat_vapour_pres.max() pres = fv['lvl'].get_value() sat_spec_hum = n.zeros(mix_rto.shape,dtype=float) for lvl_idx in range(len(pres)): sat_spec_hum[lvl_idx] = 0.622 * sat_vapour_pres[lvl_idx] \ / pres[lvl_idx] #sat_spec_hum = ma.masked_where(sat_spec_hum <= 1.e-9, #sat_spec_hum = ma.masked_where(sat_spec_hum <= 1.e-6, #sat_spec_hum) rh = n.zeros(mix_rto.shape,dtype=float) rh = spec_hum / sat_spec_hum * 100 #sys.exit() del mix_rto, spec_hum, temp, sat_vapour_pres, pres, sat_spec_hum #cheat = rh #field = n.zeros(cheat.shape) #for lat_idx in range(cheat.shape[0]): # field[-(lat_idx+1)] = cheat[lat_idx] cheat = rh field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['rh'] = { 'long_name' : '3D relative humidity', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : rh.filled(fill_value=0.), #'field' : rh, 'field' : field, 'lvl' : fv['lvl'].get_value() } #sys.exit() # Let's prepare the 3d zonal wind data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 2 -> momentum pdtmpl[0] = 2 # Oct 11: Parameter number (see Code table 4.2). 2 -> u-component of the # wind pdtmpl[1] = 2 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> isobaric level pdtmpl[9] = 100 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['zonal_wnd'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['zonal_wnd'] = { 'long_name' : '3D zonal wind', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['zonal_wnd'].get_value()[0], #'field' : fv['zonal_wnd'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx], 'field' : field, #'field' : '', 'lvl' : fv['lvl'].get_value() } # Let's prepare the 3d zonal wind data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 2 -> momentum pdtmpl[0] = 2 # Oct 11: Parameter number (see Code table 4.2). 2 -> v-component of the # wind pdtmpl[1] = 3 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> isobaric level pdtmpl[9] = 100 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat =fv['merid_wnd'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['merid_wnd'] = { 'long_name' : '3D meridional wind', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['merid_wnd'].get_value()[0], #'field' : fv['merid_wnd'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx], 'field' : field, #'field' : '', 'lvl' : fv['lvl'].get_value() } # Let's prepare the 3d geopotential height data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 3 -> mass pdtmpl[0] = 3 # Oct 11: Parameter number (see Code table 4.2). 5 -> geopotential height pdtmpl[1] = 5 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> isobaric level pdtmpl[9] = 100 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['geop_ht'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['geop_ht'] = { 'long_name' : 'geopotential height', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['geop_ht'].get_value()[0], #'field' : fv['geop_ht'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx], 'field' : field, #'field' : '', 'lvl' : fv['lvl'].get_value() } ####################################################################### # Finished with the atmospheric discipline # Let's move on to the oceanographic one ####################################################################### if 0: # Preliminaries: we need to unpack the land_sea_ice mask from the laps # files. #sea_lnd_ice_mask = fv['sea_lnd_ice_mask'].get_value()[0] sea_lnd_ice_mask = fv['sea_lnd_ice_mask'].get_value()[0,lat_cheat_idx:,:-lon_cheat_idx] shape = sea_lnd_ice_mask.shape land_cover = n.zeros(shape, 'i') ice_cover = n.zeros(shape, 'i') # I assume (know) that the mask has 3 dimensions of which the first is time # which I am going to ignore for now # I am sure this could be done more elegantly with a numpy in-built # function, but I am not sure about which one would it be for i in range(shape[0]): for j in range(shape[1]): if sea_lnd_ice_mask[i,j] == 0: land_cover[i,j] = 0 ice_cover[i,j] = 0 elif sea_lnd_ice_mask[i,j] == 1: land_cover[i,j] = 1 ice_cover[i,j] = 0 elif sea_lnd_ice_mask[i,j] == 2: # the BoM only considers ice over water land_cover[i,j] = 0 ice_cover[i,j] = 1 else: print 'Illegal value in the sea_land_ice mask!' sys.exit() dummy = data[-1]['discipline']['ocean']['fields']['2d'] # Let's prepare the ice cover data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 2 -> ice pdtmpl[0] = 2 # Oct 11: Parameter number (see Code table 4.2). 0 -> ice-cover # 1=ice, 0=no-ice (over water) pdtmpl[1] = 0 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = ice_cover field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['ice_cover'] = { 'long_name' : 'ice cover', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : ice_cover 'field' : field } ####################################################################### # Finished with the oceanographic discipline # Let's move on to the land one ####################################################################### dummy = data[-1]['discipline']['land']['fields']['2d'] if 0: # Let's prepare the land cover data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> vegetation/biomass pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 0 -> land-cover # (1=land, 0=ocean) pdtmpl[1] = 0 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. # I seem to remember there was a problem with the land cover that I could # fix using the decimal scale factor... let's try land_cover_decimal_factor = 1 # Oct: 18-19. Decimal scale factor (D) #drtmpl[2] = land_cover_decimal_factor #print 'land_cover drtmpl', drtmpl cheat = land_cover * 10**land_cover_decimal_factor field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['land_cover'] = { 'long_name' : 'land cover', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : n.array([0,0,land_cover_decimal_factor,24,0]), #'drtmpl' : drtmpl, #'field' : land_cover * 10**land_cover_decimal_factor 'field' : field } # Now put it back to what it was # Oct: 18-19. Decimal scale factor (D) #drtmpl[2] = 0 if 1: # Let's prepare the topography data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> vegetation/biomass pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 7 -> model terrain height pdtmpl[1] = 7 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 1 -> surface pdtmpl[9] = 1 # Oct 24: Scale factor of first fixed surface pdtmpl[10] = 0 # Oct 25-28: Scaled value of first fixed surface pdtmpl[11] = 0 # Oct 29: Type of second fixed surface (see Code table 4.5) pdtmpl[12] = missing(1) # Oct 30: Scale factor of second fixed surface pdtmpl[13] = missing(1) # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. cheat = fv['topog'].get_value()[lat_cheat_idx:,:-lon_cheat_idx] field = n.zeros(cheat.shape) for lat_idx in range(cheat.shape[0]): field[-(lat_idx+1)] = cheat[lat_idx] dummy['topog'] = { 'long_name' : 'topography', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #'field' : fv['topog'].get_value() #'field' : fv['topog'].get_value()[lat_cheat_idx:,:-lon_cheat_idx] 'field' : field } dummy = data[-1]['discipline']['land']['fields']['3d'] if 1: # Let's prepare the 3d soil temperature data and metadata # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> vegetation/biomass pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 2 -> soil temperature pdtmpl[1] = 2 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 106 -> depth below land surface pdtmpl[9] = 106 # Oct 24: Scale factor of first fixed surface # 3 -> preserve mm accuracy # 2 -> preserve cm accuracy soil_lvl_scaling = 2 pdtmpl[10] = soil_lvl_scaling # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) # 106 -> depth below land surface pdtmpl[12] = 106 # Oct 30: Scale factor of second fixed surface # 3 -> preserve mm accuracy pdtmpl[13] = soil_lvl_scaling # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. # Let's interpolate linearly to the Noah levels # NB for the deepest level we actually have an extrapolation laps_soil_lvl = fv['soil_lvl'].get_value() noah_soil_lvl = [0.10, 0.40, 1.0, 2.0] #laps_soil_temp = fv['soil_temp'].get_value()[0] laps_soil_temp = fv['soil_temp'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] shape = laps_soil_temp.shape noah_soil_temp = n.zeros(shape, dtype=float) max_soil_lvl = 4 for soil_idx in range(max_soil_lvl): if soil_idx < max_soil_lvl - 1: x_a = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx] x_b = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx + 1] y_a = laps_soil_temp[soil_idx] y_b = laps_soil_temp[soil_idx + 1] else: # this accomplishes the extrapolation x_a = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx - 1] x_b = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx] y_a = laps_soil_temp[soil_idx - 1] y_b = laps_soil_temp[soil_idx] #print x_a, '\n\n', x_b, '\n\n', y_a, '\n\n', y_b, '\n\n' #print soil_idx, x_a.shape, x_b.shape, y_a.shape, y_b.shape x = n.ones((shape[1],shape[2])) * noah_soil_lvl[soil_idx] noah_soil_temp[soil_idx] = vu.lin_interp(x_a, x_b, y_a, y_b, x) cheat = noah_soil_temp field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['soil_temp'] = { 'long_name' : '3D volumetric soil temperature', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #TODO: make sure the scaling is right #'field' : noah_soil_temp, 'field' : field, 'lvl' : noah_soil_lvl } if 1: # Let's prepare the 3d volumetric soil moisture # Product Definition Template Number (see Code Table 4.0) # 0 -> Analysis or forecast at a horizontal level or in a # horizontal layer at a point in time. (see Template 4.0) pdtnum = 0 # array of zeros (int) of size 15 pdtmpl = n.zeros((15,),dtype=int) # Follows the definition of template 4.0 # Oct 10: Parameter category (see Code table 4.1). 0 -> vegetation/biomass pdtmpl[0] = 0 # Oct 11: Parameter number (see Code table 4.2). 9 -> volumetric soil # moisture # pdtmpl[1] = 9 # I don't like this choice but maybe it is going to work more easily with # WRF and other NCAR utilities like cnvgrib pdtmpl[1] = 192 # Oct 12: Type of generating process (see Code table 4.3). 2 -> forecast pdtmpl[2] = 2 # Oct 13: Background generating process identifier # (defined by originating centre) pdtmpl[3] = missing(1) # Oct 14: Analysis or forecast generating process identified # (defined by originating centre) pdtmpl[4] = missing(1) # Oct 15-16: Hours of observational data cutoff after reference time pdtmpl[5] = missing(2) # Oct 17: Minutes of observational data cutoff after reference time pdtmpl[6] = missing(1) # Oct 18: Indicator of unit of time range (see Code table 4.4) # 1 -> hours # NB WPS expects this to be hours! pdtmpl[7] = 1 # Oct 19-22: Forecast time in units defined by octet 18 pdtmpl[8] = fv['forc_hrs'].get_value()[0] # Oct 23: Type of first fixed surface (see Code table 4.5) # 106 -> depth below land surface pdtmpl[9] = 106 # Oct 24: Scale factor of first fixed surface # 2 -> preserve cm accuracy pdtmpl[10] = 2 # Oct 25-28: Scaled value of first fixed surface # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding pdtmpl[11] = missing(3) # Oct 29: Type of second fixed surface (see Code table 4.5) # 106 -> depth below land surface pdtmpl[12] = 106 # Oct 30: Scale factor of second fixed surface # 2 -> preserve cm accuracy pdtmpl[13] = 2 # Oct 31-34: Scaled value of second fixed surfaces pdtmpl[14] = missing(3) # -> missing as there is will be defined later when the data is sliced in # horizontal slabs for its encoding # NB I am using the same drtnum and drtmpl for all data so I want # replicate it in the following fields if you need to change this on a # per-variable basis this is the place to redefine it. If you choose to do # so then make sure you define one for each variable that follows or be # aware that the last one defined will apply to all that follows. # in the ECMWF soil scheme the data is saved scaled by the depth of the # soil layer (or so I understand...) # TODO work out this 'not writable' rubbish #print laps_soil_mois.flags['WRITEABLE'] #laps_soil_mois.flags['WRITEABLE'] = True #laps_soil_mois = fv['soil_mois'].get_value()[0] #laps_soil_mois = fv['soil_mois'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] tmp = fv['soil_mois'].get_value()[0,:,lat_cheat_idx:,:-lon_cheat_idx] shape = tmp.shape laps_soil_mois = n.zeros(shape) laps_soil_lvl = fv['soil_lvl'].get_value() for soil_idx in range(len(laps_soil_lvl)): #laps_soil_mois[soil_idx] /= laps_soil_lvl[soil_idx] laps_soil_mois[soil_idx] = tmp[soil_idx] / laps_soil_lvl[soil_idx] # Let's interpolate linearly to the Noah levels # NB for the deepest level we actually have an extrapolation noah_soil_lvl = [0.10, 0.40, 1.0, 2.0] noah_soil_mois = n.zeros(laps_soil_mois.shape, dtype=float) max_soil_lvl = 4 for soil_idx in range(max_soil_lvl): if soil_idx < max_soil_lvl - 1: x_a = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx] x_b = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx + 1] y_a = laps_soil_mois[soil_idx] y_b = laps_soil_mois[soil_idx + 1] else: # this accomplishes the extrapolation x_a = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx - 1] x_b = n.ones((shape[1],shape[2])) * laps_soil_lvl[soil_idx] y_a = laps_soil_mois[soil_idx - 1] y_b = laps_soil_mois[soil_idx] x = n.ones((shape[1],shape[2])) * noah_soil_lvl[soil_idx] noah_soil_mois[soil_idx] = vu.lin_interp(x_a, x_b, y_a, y_b, x) cheat = noah_soil_mois field = n.zeros(cheat.shape) for lvl_idx in range(cheat.shape[0]): for lat_idx in range(cheat.shape[1]): field[lvl_idx,-(lat_idx+1)] = cheat[lvl_idx,lat_idx] dummy['soil_mois'] = { 'long_name' : '3D volumetric soil moisture', 'pdtnum' : pdtnum, 'pdtmpl' : pdtmpl, 'drtnum' : drtnum, 'drtmpl' : drtmpl, #TODO: make sure the scaling is right #'field' : noah_soil_mois, 'field' : field, 'lvl' : noah_soil_lvl } f.close() sfc_f.close() #sys.exit() #p.show() # Multiple times can be included in a grib files as sequential messages so the # outermost 'dimension' of the data structure should be time # As I am going to use the GFS grib files from ncep as a template for my own # files, I will include only one time in each of the files. # It should be easy to wrap to modify the following code to concatenate # multiple times def do_the_encoding(data): dummy_idx = 0 for time_slice in data: g2e = grib2.grib2.Grib2Encode( time_slice['discipline']['atmos']['code'], time_slice['idsect']) g2e.addgrid(time_slice['gdsinfo'],time_slice['gdtmpl']) dummy = time_slice['discipline']['atmos']['fields']['2d'] for var_key in dummy.keys(): print 'Processing ' + dummy[var_key]['long_name'] g2e.addfield( dummy[var_key]['pdtnum'], dummy[var_key]['pdtmpl'], dummy[var_key]['drtnum'], dummy[var_key]['drtmpl'], dummy[var_key]['field'] ) dummy = time_slice['discipline']['atmos']['fields']['3d'] for var_key in dummy.keys(): print 'Processing ' + dummy[var_key]['long_name'] for lvl_idx in range(len(dummy[var_key]['lvl'])): print '\tlevel ' + str(lvl_idx) # '*100' because it must be in Pa dummy[var_key]['pdtmpl'][11] = \ dummy[var_key]['lvl'][lvl_idx]*100 g2e.addfield( dummy[var_key]['pdtnum'], dummy[var_key]['pdtmpl'], dummy[var_key]['drtnum'], dummy[var_key]['drtmpl'], dummy[var_key]['field'][lvl_idx] ) g2e.end() file_name = 'laps_' + str(time_slice['idsect'][5]) + '_' \ + str(zfill(time_slice['idsect'][6],2)) + '_' \ + str(zfill(time_slice['idsect'][7],2)) + '_' \ + str(zfill(time_slice['idsect'][8],2)) + 'Z_' \ + str(zfill(dummy[var_key]['pdtmpl'][8],2)) + '_hrs_fcst.grb2' #f = open('dummy' + str(dummy_idx) + '.grb2', 'w') f = open(file_name, 'w') f.write(g2e.msg) f.close() g2e = grib2.grib2.Grib2Encode( time_slice['discipline']['ocean']['code'], time_slice['idsect']) g2e.addgrid(time_slice['gdsinfo'],time_slice['gdtmpl']) dummy = time_slice['discipline']['ocean']['fields']['2d'] for var_key in dummy.keys(): print 'Processing ' + dummy[var_key]['long_name'] g2e.addfield( dummy[var_key]['pdtnum'], dummy[var_key]['pdtmpl'], dummy[var_key]['drtnum'], dummy[var_key]['drtmpl'], dummy[var_key]['field'] ) if 0: dummy = time_slice['discipline']['ocean']['fields']['3d'] for var_key in dummy.keys(): print 'Processing ' + dummy[var_key]['long_name'] for lvl_idx in range(len(dummy[var_key]['lvl'])): print '\tlevel ' + str(lvl_idx) # '*100' because it must be in Pa dummy[var_key]['pdtmpl'][11] = \ dummy[var_key]['lvl'][lvl_idx]*100 g2e.addfield( dummy[var_key]['pdtnum'], dummy[var_key]['pdtmpl'], dummy[var_key]['drtnum'], dummy[var_key]['drtmpl'], dummy[var_key]['field'][lvl_idx] ) g2e.end() #f = open('dummy' + str(dummy_idx) + '.grb2', 'a') f = open(file_name, 'a') f.write(g2e.msg) f.close() g2e = grib2.grib2.Grib2Encode( time_slice['discipline']['land']['code'], time_slice['idsect']) g2e.addgrid(time_slice['gdsinfo'],time_slice['gdtmpl']) dummy = time_slice['discipline']['land']['fields']['2d'] for var_key in dummy.keys(): print 'Processing ' + dummy[var_key]['long_name'] g2e.addfield( dummy[var_key]['pdtnum'], dummy[var_key]['pdtmpl'], dummy[var_key]['drtnum'], dummy[var_key]['drtmpl'], dummy[var_key]['field'] ) dummy = time_slice['discipline']['land']['fields']['3d'] for var_key in dummy.keys(): print 'Processing ' + dummy[var_key]['long_name'] for lvl_idx in range(len(dummy[var_key]['lvl'])): print '\tlevel ' + str(lvl_idx) if lvl_idx == 0: dummy[var_key]['pdtmpl'][11] = 0 else: dummy[var_key]['pdtmpl'][11] = \ round(dummy[var_key]['lvl'][lvl_idx-1]* 10**soil_lvl_scaling) dummy[var_key]['pdtmpl'][14] = \ round(dummy[var_key]['lvl'][lvl_idx]* 10**soil_lvl_scaling) g2e.addfield( dummy[var_key]['pdtnum'], dummy[var_key]['pdtmpl'], dummy[var_key]['drtnum'], dummy[var_key]['drtmpl'], dummy[var_key]['field'][lvl_idx] ) g2e.end() #f = open('dummy' + str(dummy_idx) + '.grb2', 'a') f = open(file_name, 'a') f.write(g2e.msg) f.close() #os.system('cnvgrib -g21 ' + file_name + ' ' + file_name[:-4] + 'grb') #os.system('/Users/val/src/cnvgrib-1.1.3/cnvgrib -g21 ' + file_name + ' ' + file_name[:-4] + 'grb') os.system('cnvgrib -g21 ' + file_name + ' ' + file_name[:-4] + 'grb') os.system('rm ' + file_name) dummy_idx += 1