import numpy as np import string import qg_vars import pygrib from ncepgrib2 import Grib2Encode import re, sys import time,subprocess from os import system """ QGtxt2GRIB.py Author: Greg Blumberg (OU/CIMMS; wblumberg@ou.edu) Creation Date: 5/13/2016 This script converts the flat text files output by Eric Thaler's QG solver directly into individual GRIB files. The code in this script is heavily derived from Paul Wolyn's original Python scripts that would convert the QG netCDF files into GRIB files. In order to run this script, you'll need to use these command line arguments: python QGtxt2GRIB.py - the model name in all capitol letters (e.g. GFS, RAP) - the model initalization time (YYYYMMDD_HHMM) - the forecast hour of the text files (090) - the path to the directory where the flat files exist - the output directory to put the GRIB files. Example: $ python QGtxt2GRIB.py GFS 20160504_1200 090 ./ ./ Requires: - PyGRIB - Numpy - qg_vars Interdispersed across this document are comments that refer to the GRIB2 Documentation http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc.shtml A lot of numbers and arrays are used that refer to things in the GRIB Tables. """ fillVal = -99999.0 CONUS211List = ["GFS_QG","NAM_QG","RAP_QG","ECM_QG","DGX_QG","GEM_QG","SRE_QG"] def createGRIBFile(model,parm,initialTime,forecastHoursList,levelList,inData,scaleFactor=1.0): """ createGRIBFile Inputs ------ model (mainly to get mapping information for gribFile) parm - Used to get get correct information for grib file initialTime - used to get initial time in gribfile forecastHoursList - List of forecast hours in seconds levelList - list for data levels indata - 2D numpy array with data to plot scalefactor is used to scale the value (if needed) Returns ------- Nothing """ #Call this function to load parmInfoDict. I use a function call to keep the code cleaner. # Defines things like the units in a dictionary. print "Getting parmInfoDict..." parmInfoDict = getparmInfo() # If the parameter is not defined in this script, return nothing. if parm not in parmInfoDict.keys(): print ("parm: "+parm+" not defined in getparmInfo") return print "Getting productID..." productID = getProductID() productIDKeys = productID.keys() # If the model type we ran the QG code is not listed, exit the program. if model not in productIDKeys: print ("Must add model: "+model+ \ " to productID dictionary in getproductID") return pdtmpl = {} # Get information for various levels in the QG data netCDF file by calling a function levelDict = definelevels(); levelDictKeys = levelDict.keys() # Create the ID section (tells the GRIB file something about the time of the data) # Also known as Section 1 of the GRIB file idsect = createidsect(initialTime) # Set the discipline of the data that is going in the GRIB (see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table0-0.shtml) discipline = 0 # This function retruns lists with the mapping information for the GRIB file print "Getting the mapping information..." # Get the mapping information if the model output is on the CONUS 211 grid # The mapping info needs to be changed for other grids. # This fills out the information for Section 3 of the GRIB file (Section 2 is local information or something like that that we don't need) if model in CONUS211List: gdsinfo, gdtmpl = getMappingInfoCONUS211() # Can add other elif's for other mappings else: print ("No mapping information for: "+model) return # Make first call to get metadata about the fields the GRIB will hold. # For each level and valid time in the below loops, make changes to specific elements # of these arrays to modify the metadata (e.g., pressure level, units, etc.) about the field # This fills out Sections 4 & 5 of the GRIB file. pdtnum, pdtmpl, drtnum, drtmpl = getFieldInfo() pdtmpl[0] = parmInfoDict[parm][0] pdtmpl[1] = parmInfoDict[parm][1] drtmpl[2] = parmInfoDict[parm][2] pdtmpl[4] = productID[model] # Loop over all the forecast times and levels, creating GRIB files for each. for i in range(len(forecastHoursList)): # Loop over the forecast hours for j in range(len(levelList)): # Loop over the vertical levels if 'MB ' + str(levelList[j]) not in levelDictKeys: # If you don't want to make a GRIB file for this level, don't. print ("No level information for: "+levelList[j]) continue # Create a grib file for each time and level. Provide information # on the discipline, the ID, and the Grid Definition Template. grib = openGribFile(discipline, idsect,gdsinfo,gdtmpl) # Alter some values in pdtmpl for this level # These changes are defined in def definelevels() for k in range(len(levelDict['MB ' + str(levelList[j])])): # Add in more information about the units and pressure levels to the Product Definition Template pdtmpl[k+9] = levelDict['MB ' + str(levelList[j])][k] # Define valid time of the grid (Need to update the Product Definition Template) pdtmpl[8] = int(forecastHoursList[i])*3600 # Dimensions of inData are [forecastHour][level][x][y] packData = inData[i][j][:][:] # Sending a masked array into addfield allows for the # GRIB encoder to account for missing data, so use Numpy to make a masked array # Missing values in netcdf are 1.e+37 # Create a masked array of missing values missingArray = np.zeros_like(packData,dtype = bool) missingArray = np.where(packData > 1.e+30, True,False) # Adjust grid based on scalar and do not adjust missing values packData = np.where(packData > 1.e+30, packData,packData*scaleFactor) # Creating the masked array encodeArray = np.ma.array(packData,mask = missingArray) # Add the field to the grib file. There can be multiple # fields in a grib2 file as long as they all have the same # grid mapping information described in Section 3. This command takes in the Product Definition Template information # the Data Representation Template information and the array to be encoded into the GRIB file. grib.addfield(pdtnum,pdtmpl,drtnum,drtmpl,encodeArray) # Add end of grib2 file grib.end() # Create output GRIB filename #levelForGrib = re.sub('\s+','',str('MB ' + levelList[j])) levelForGrib = str(levelList[j]) fhr = str(forecastHoursList[i]) outGribDir = outputDir outGribFile = "%s/LDAD-GRIB-%s-%s-%s-%s-%s.grb2"%\ (outGribDir,model,parm,initialTime,levelForGrib,fhr) # Write the GRIB message to a file. writeAndProcessGrib(outGribFile,grib.msg) def openGribFile(discipline, idsect,gdsinfo,gdtmpl): """ This creates an instance of the grib file. The contents of the grib file is written to grib.msg. The grib object is returned so the user can add the field (the data) later. Adds information to the GRIB file about the grid. Inputs ------ discipline - see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table0-0.shtml idsect - see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect1.shtml gdsinfo - see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect3.shtml gdtmpl - see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-1.shtml """ grib = Grib2Encode(discipline, idsect) # This adds the grid (mapping) information to the grib2 file grib.addgrid(gdsinfo,gdtmpl) return grib def writeAndProcessGrib(outGribFile,gribContents): """ Write the GRIB file. """ print ("Writing: %s"%(outGribFile)) gribOut = open(outGribFile,"wb") gribOut.write(gribContents) gribOut.close() def createidsect(initialTime): """ Input is initial time from the filename. It has the form: YYYYMMDD_HHMM Creates an array for Section 1 (the identification section): See: http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect1.shtml Returns ------- idsect (list) """ year = int(initialTime[0:4]) month = int(initialTime[4:6]) day = int(initialTime[6:8]) hour = int(initialTime[9:11]) minute = int(initialTime[11:13]) idsect = [] # 1- Originating center = 9 (other NWS) idsect.append(9) # 2- Subcenter WFO BOU = 29 idsect.append(29) # 3- Grib2 version (try 14 from Nov 5, 2014) idsect.append(14) # 4- local tables used (Try 1) - Start time of model/forecast idsect.append(1) # 5- Reference Time = Try start of forecast (1) idsect.append(1) # 6-12 - Reference time 5 = year (4 digits) 6= month, 7=day, 8 = hour, 9=minute, 10=seconds idsect.append(year) idsect.append(month) idsect.append(day) idsect.append(hour) idsect.append(minute) idsect.append(0) # 13- Production status (0 = operational) idsect.append(0) # 14- Type of processed data (1 = forecast) idsect.append(1) return idsect def getMappingInfoCONUS211(): """ Returns two lists which are used by the grib2 file for mapping Grid definition for CONUS211 (Lambert Conformal) See Section 3 for more info: http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect3.shtml Returns ------- gdsinfo (list) - Grid Definition Section (see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect3.shtml) gdtmpl (list) - Grid Definition Template (see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-1.shtml) """ # Define the gdsinfo array gdsinfo = [] # 0- Grid definition givien in detailed information =0, predefined grid =1 (wish I could get this to work) gdsinfo.append(0) # 1- Number of grid points (Nx*ny) # for CONUS 211, Nx = 93, Ny = 65 if (model == "RAP_QG") : # The RAP QG solutions are output on a different grid for some reason. gdsinfo.append(76*57) else : gdsinfo.append(93*65) # 2- Not sure what this is, but it =0 for a regular grid gdsinfo.append(0) # 3- Not sure what this is, but it is 0 if no additional list gdsinfo.append(0) # 4- grid definition template. Lambert conformal is 30 gdsinfo.append(30) # This list contains grid data information for the table/projection specified on gdsinfo[4] # See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-1.shtml gdtmpl = [] # 0- How shape of the earth is defined http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-2.shtml gdtmpl.append(6) # Earth assumed spherical with raduis of 6371229m # 1- Scale factor for radius of the earth? Assume o since spherical? gdtmpl.append(0) # 2- Radius of the earth in meters. (Got value from NDFD infomation) gdtmpl.append(0) # 3-6- Since earth is spherical, next 4 entries are zero gdtmpl.append(0) gdtmpl.append(0) gdtmpl.append(0) gdtmpl.append(0) # 7- Nx if (model == "RAP_QG"): # Custom results for the RAP QG solution. gdtmpl.append(76) else: gdtmpl.append(93) # 8- Ny if (model == "RAP_QG"): gdtmpl.append(57) else: gdtmpl.append(65) # 9-10- Latitude and longitude of first grid point (must be integers *10E_6? if (model == "RAP_QG"): gdtmpl.append(16280000) gdtmpl.append(233862000) # must be between 0 and 360? else: gdtmpl.append(12190000) gdtmpl.append(226541000) # must be between 0 and 360? # 11- Resolutuion and flag components in 8 digit bindary? # digit 1 = 0, 2= 0, 3= 1, 4= 1, 5= 1, 6=0, 7=0, 8=0 ) # This gives 00110000 = 32+16+8 = 56? (set dogot 5 to 1 for winds (vectors))? gdtmpl.append(56) # Binary 0011100 # 12-13- Latitude and longitude where delx and dely are specified *1.e+6 (From netcdf file) gdtmpl.append(25000060) gdtmpl.append(265000000) # -100.5548 # gdtmpl.append(40606060) # gdtmpl.append(259445200) # -100.5548 # 14-15- The delx and dely at the lat/lon specified previously (in meters*1E+3) (from netcdf file) gdtmpl.append(81270500) gdtmpl.append(81270500) # gdtmpl.append(78041740) # gdtmpl.append(78042930) # 16- Projection Center flag. 8 digit binary (and it is all zeroes) gdtmpl.append(0) # 17- Scanning mode, This is another 8 digit binary. From the local WRF, it is 01000000 = 64 gdtmpl.append(64) # 18- Lat 1 in lat*1E+6 gdtmpl.append(25000000) # 19- Lat 2 in lat*1E+6 gdtmpl.append(25000000) # 20-21- Lat and Lon of South Pole Projection. Since not in projection set to 0 gdtmpl.append(0) gdtmpl.append(0) return gdsinfo, gdtmpl def getFieldInfo(): """ This routine defines most of the grib2 information for the QGgh field QGgh means the QG Geopotential Height field. See Sections 4 and 5 for Tables 4 and 5: http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect4.shtml http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect5.shtml Returns ------- pdtnum (int) - Product Definition Template Number (see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-0.shtml) pdtmpl (list) - Product Definition Template drtnum (int) - Data Representation Template Number (see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table5-0.shtml) drtmpl (list) - Data Representation Template """ pdtnum = 0 # Since this code works on NWP forecast data, use Table 4.0 pdtmpl = [] # See : http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp4-0.shtml # 0 - Parameter catagory - Try mass (3) for geopotential height pdtmpl.append(3) # 1 - Under mass, geopotential height is: 5 pdtmpl.append(5) # 2 - Generating process- 2= forecast pdtmpl.append(2) # 3 - Background generating process? Try 0 pdtmpl.append(0) # 4 - Forecast generating process ID. Try 255 for missing pdtmpl.append(0) # 5&6 - Hours and minutes of obseravtional cutoff- Not applicable to this code pdtmpl.append(0) pdtmpl.append(0) # 7 - Indicator of time unit range (11= 6 hours, 10 = 3 hours, 13 = seconds) pdtmpl.append(13) # 8 - forccast time in units defined by previous entry. pdtmpl.append(2) # 9 - Fixed surface type; 100 means it's isobaric in Pa # see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-5.shtml pdtmpl.append(100) # 10 - Scale factor of first surface. Try 0 pdtmpl.append(0) # 11 - Value of first surface. 500mb = 50000 Pa pdtmpl.append(50000) # 12 - Second Fixed surface type. If none enter missing (255) pdtmpl.append(255) # 13 - Scale factor of first surface. Try 0 pdtmpl.append(0) # 14 - None so try 0 pdtmpl.append(0) # Data representation number 0= grid point data drtnum = 0 drtmpl = [] # 0-3- Defined later in packing?, set to zero # These will change with packing routine. drtmpl.append(0) drtmpl.append(0) drtmpl.append(0) drtmpl.append(0) # 4- data type 0=floating point drtmpl.append(0) return pdtnum, pdtmpl, drtnum, drtmpl def definelevels(): """ Returns a dictionary which gives level information based in the lvlel string in the netcdf file These are elements 9-14 in drtmpl (see the append statements in getFieldInfo()) These describe the different pressure levels that are exported by the QG code. # I think the first element of the array comes from: http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-5.shtml """ levelDict = {} levelDict["MB 1000"] = [100, 0, 100000, 255, 0, 0] levelDict["MB 950"] = [100, 0, 95000, 255, 0, 0] levelDict["MB 900"] = [100, 0, 90000, 255, 0, 0] levelDict["MB 850"] = [100, 0, 85000, 255, 0, 0] levelDict["MB 800"] = [100, 0, 80000, 255, 0, 0] levelDict["MB 750"] = [100, 0, 75000, 255, 0, 0] levelDict["MB 700"] = [100, 0, 70000, 255, 0, 0] levelDict["MB 650"] = [100, 0, 65000, 255, 0, 0] levelDict["MB 600"] = [100, 0, 60000, 255, 0, 0] levelDict["MB 550"] = [100, 0, 55000, 255, 0, 0] levelDict["MB 500"] = [100, 0, 50000, 255, 0, 0] levelDict["MB 450"] = [100, 0, 45000, 255, 0, 0] levelDict["MB 400"] = [100, 0, 40000, 255, 0, 0] levelDict["MB 350"] = [100, 0, 35000, 255, 0, 0] levelDict["MB 300"] = [100, 0, 30000, 255, 0, 0] levelDict["MB 250"] = [100, 0, 25000, 255, 0, 0] levelDict["MB 200"] = [100, 0, 20000, 255, 0, 0] levelDict["MB 150"] = [100, 0, 15000, 255, 0, 0] levelDict["MB 100"] = [100, 0, 10000, 255, 0, 0] levelDict["SFC"] = [ 1, 0, 0 , 255, 0, 0] return levelDict def getparmInfo(): """ Returns a dictionary with parameter info and scaling factor for each parameter. Parameters are different variables and the units that will be contained in the GRIB file. The dictionary contains arrays whose elements are: Element 0 = pdtmpl[0] (see parameter category http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-1.shtml) Element 1 = pdtmpl[1] (see parameter number (requires category) http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-2.shtml) Element 2 can be considered number of decimal places to include when saving the data to the grib2 file. A value of 4 will save the data to 4 decimal places (ie. 0.0001) Since omeg is close to zero, use 4 decimal places. """ parmInfoDict = {} parmInfoDict["QGgh"] = [3, 5, 1] parmInfoDict["omeg"] = [2, 8, 4] # in Pa/s, netcdf is in mb/hr parmInfoDict["QGtv"] = [0, 0, 2] #parmInfoDict["QGrunp"] = [1, 8, 3] # in meters parmInfoDict["QGpcpn"] = [1, 7, 8] # in meters/second parmInfoDict["QGcond"] = [1, 217, 12] parmInfoDict["QGaccp"] = [1, 8, 3] # in meters parmInfoDict["QGug"] = [2 , 2, 4] parmInfoDict["QGvg"] = [2 , 3, 4] parmInfoDict["psig0"] = [2, 4, 4] parmInfoDict["omeg0"] = [2, 8, 4] parmInfoDict["gpht0"] = [3,2,4] parmInfoDict["gpht"] = [3,2,6] parmInfoDict["p"] = [3, 5, 1] return parmInfoDict # This defines the model and product ID. def getProductID(): """ Returns a dictionary with the product ID. I think this is an ID that defines where the data came from. This modifies pdtmpl[4] for each model. In a sense, we are creating a product ID table for center 9 (rest of NWS) subcenter 29 )WFO BOU) The key must match a value of the variable "model" in "def netcdfToGrib2" """ productID = {} productID["DGX_QG"] = 0 productID["ECM_QG"] = 1 productID["GEM_QG"] = 2 productID["GFS_QG"] = 3 productID["NAM_QG"] = 4 productID["RAP_QG"] = 5 productID["SRE_QG"] = 6 return productID # Returns a 3-d cube of data 2-d * time def readSfcGrid(varName, model, timeStr, path, fcstHour): cube = [] fileName = varName + "_" + model + "_" + timeStr + "_" + fcstHour + "_srfc" + ".dat" fileName = path + fileName grid = readFile(fileName) cube.append(grid) cube = np.array(cube) return cube # Returns 4-d cube of data 2-d * height * time def read4DCube(varName, model, timeStr, path, fcstHour): # read each level, time and return the 4-d cube bigCube = [] cube = [] for level in range(1000, 50, -50): levelStr = str(level) while len(levelStr) < 4 : levelStr = "0" + levelStr fileName = varName + "_" + model + "_" + timeStr + "_" + fcstHour + "_" + levelStr + ".dat" fileName = path + fileName grid = readFile(fileName) cube.append(grid) cube = np.array(cube) bigCube.append(cube) return np.array(bigCube) def readFile(filename): """ Reads a flat file from the QG solver into an array. """ return np.genfromtxt(filename, dtype=float) # Print out usage statement. print 'python QGtxt2GRIB.py ' print 'EX: python QGtxt2GRIB.py GFS 20160504_1200 090 ./ ./' pcpn_only = 'n' Model=sys.argv[1] # model name runTime = initialTime = sys.argv[2]#'20160504_1200' dirName = sys.argv[4] # input directory outputDir = sys.argv[5] fcstHour= sys.argv[3] levelList = np.arange(1000, 50, -50) forecastHoursList = [fcstHour] data = {} # Read in the forecast parameters and make GRIBs for j in range(0,len(qg_vars.ForecastParms)) : print "Reading file",qg_vars.ForecastParms[j][0],"to populate netCDF variable",qg_vars.ForecastParms[j][1],qg_vars.ForecastParms[j][2] d = read4DCube(qg_vars.ForecastParms[j][0], Model, runTime, dirName, fcstHour) d = np.where(np.equal(d, fillVal), 1E37, d) print "\tShape of the array:",d.shape #store4DCube(myFile,qg_vars.ForecastParms[j][1],d) data[qg_vars.ForecastParms[j][1]] = d # Read in the forecast surface parameters and make GRIBs for j in range(0,len(qg_vars.SurfaceForecastParms)) : if (Model == "GEN") and ((qg_vars.SurfaceForecastParms[j][0] == "cond") or (qg_vars.SurfaceForecastParms[j][0] == "accp") or (qg_vars.SurfaceForecastParms[j][0] == "runp")) : continue print "Reading file",qg_vars.SurfaceForecastParms[j][0],"to populate netCDF variable",qg_vars.SurfaceForecastParms[j][1],qg_vars.SurfaceForecastParms[j][2] d = readSfcGrid(qg_vars.SurfaceForecastParms[j][0], Model, runTime, dirName, fcstHour) d = np.where(np.equal(d,fillVal), 1E37, d) print "\tShape of the array:",d.shape data[qg_vars.SurfaceForecastParms[j][1]] = d #storeSfcGrids(myFile,SurfaceForecastParms[j][1],d) # Read in supplimental parameters and make the GRIBs for j in range(0,len(qg_vars.SupplementalParms)) : print "Reading file",qg_vars.SupplementalParms[j][0],"to populate netCDF variable",qg_vars.SupplementalParms[j][1],qg_vars.SupplementalParms[j][2] d = read4DCube(qg_vars.SupplementalParms[j][0], Model, runTime, dirName, fcstHour) d = np.where(np.equal(d,fillVal), 1E37, d) #store4DCube(myFile,qg_vars.SupplementalParms[j][1],d) print "\tShape of the array:",d.shape data[qg_vars.SupplementalParms[j][1]] = d write_yn = 'n' if (write_yn == "y") : for j in range(0,len(qg_vars.EducationalParms)) : print "Reading file",qg_vars.EducationalParms[j][0],"to populate netCDF variable",qg_vars.EducationalParms[j][1],qg_vars.EducationalParms[j][2] d = read4DCube(qg_vars.EducationalParms[j][0], Model, runTime, dirName, fcstHour) d = np.where(np.equal(d, fillVal), 1E37, d) #store4DCube(myFile,EducationalParms[j][1],d) for j in range(0,len(qg_vars.SurfaceEducationalParms)) : print "Reading file",qg_vars.SurfaceEducationalParms[j][0],"to populate netCDF variable",qg_vars.SurfaceEducationalParms[j][1],qg_vars.SurfaceEducationalParms[j][2] d = readSfcGrid(qg_vars.SurfaceEducationalParms[j][0], Model, runTime, dirName, fcstHour) d = np.where(np.equal(d,fillVal), 1E37, d) # storeSfcGrids(myFile,SurfaceEducationalParms[j][1],d) else : print "Reading file accp to populate netCDF variable QGaccp : Accumulated QG precipitation" d = readSfcGrid("accp", Model, runTime, dirName, fcstHour) d = np.where(np.equal(d,fillVal), 1E37, d) #storeSfcGrids(myFile,"QGaccp",d) # List parmList that tells the program which GRIBs to create parmList = ["psig0","omeg0","gpht0","QGgh","QGtv","omeg","gpht","QGpcpn","QGcond"]####,"QGaccp" Model = model = Model.upper() + '_QG' # Loop over the variables and create GRIB files for each one for parm in parmList: print "Making GRIBs for the parameter:", parm if parm in ['QGrunp', 'QGpcpn', 'QGaccp']: scaleFactor = 1000. else: scaleFactor = 1.0 forecastHourList = [fcstHour] inData = data[parm] # Make a GRIB! createGRIBFile(Model, parm, initialTime, forecastHoursList, levelList, inData, scaleFactor=scaleFactor)