## # This software was developed and / or modified by Raytheon Company, # pursuant to Contract DG133W-05-CQ-1067 with the US Government. # # U.S. EXPORT CONTROLLED TECHNICAL DATA # This software product contains export-restricted data whose # export/transfer/disclosure is restricted by U.S. law. Dissemination # to non-U.S. persons whether in the United States or abroad requires # an export license or other authorization. # # Contractor Name: Raytheon Company # Contractor Address: 6825 Pine Street, Suite 340 # Mail Stop B8 # Omaha, NE 68106 # 402.291.0100 # # See the AWIPS II Master Rights File ("Master Rights File.pdf") for # further licensing information. ## import grib2 import numpy from math import pow import time, os, sys, math import LogStream import tempfile from matplotlib.mlab import griddata from java.lang import Float from java.lang import Double from java.lang import Integer from java.util import Calendar from java.util import Date from java.util import GregorianCalendar from javax.measure.unit import SI from javax.measure.unit import NonSI from javax.measure.unit import Unit from com.raytheon.uf.common.time import DataTime from com.raytheon.uf.common.time import TimeRange from com.raytheon.uf.common.geospatial import MapUtil from com.raytheon.uf.common.serialization import SerializationUtil from com.raytheon.uf.common.dataplugin.grid import GridRecord from com.raytheon.uf.common.gridcoverage import LambertConformalGridCoverage from com.raytheon.uf.common.gridcoverage import LatLonGridCoverage from com.raytheon.uf.common.gridcoverage import MercatorGridCoverage from com.raytheon.uf.common.gridcoverage import PolarStereoGridCoverage from com.raytheon.uf.common.gridcoverage.lookup import GridCoverageLookup from com.raytheon.uf.common.gridcoverage import Corner from com.raytheon.edex.plugin.grib.util import GribModelLookup from com.raytheon.uf.common.dataplugin.level.mapping import LevelMapper from com.raytheon.uf.common.dataplugin.level import Level from com.raytheon.uf.common.dataplugin.level import LevelFactory from com.raytheon.edex.plugin.grib.spatial import GribSpatialCache from com.raytheon.edex.util.grib import GribTableLookup from com.raytheon.edex.util import Util from com.raytheon.edex.plugin.grib import Grib1Decoder from com.raytheon.edex.util.grib import GribParamTranslator from com.raytheon.uf.common.parameter import Parameter; from com.raytheon.uf.common.parameter.mapping import ParameterMapper; PLUGIN_NAME = "grid" # Static values for accessing parameter lookup tables PARAMETER_TABLE = "4.2" GENPROCESS_TABLE = "A" LEVELS_TABLE = "4.5" DOT = "." DASH = "-" SPACE = " " MISSING = "Missing" # Static values for converting forecast times to seconds SECONDS_PER_MINUTE = 60 SECONDS_PER_HOUR = 3600 SECONDS_PER_DAY = 86400 # Assumes 31 days in 1 month SECONDS_PER_MONTH = 2678400 #Assumes 365 days in 1 year SECONDS_PER_YEAR = 977616000 # Default values for earth shape MAJOR_AXIS_DEFAULT = 6378160 MINOR_AXIS_DEFAULT = 6356775 # Default values for dx/dy spacing of grids DEFAULT_SPACING_UNIT = "km" DEFAULT_SPACING_UNIT2 = "degree" # Quasi-regular values (grids 37-44) THINNED_GRID_PTS = 73 THINNED_GRID_MIDPOINT = 37 THINNED_GRID_SPACING = 1.25 THINNED_GRID_SIZE = 3447 THINNED_GRID_REMAPPED_SIZE = 5329 THINNED_XI_LINSPACE = numpy.linspace(0, THINNED_GRID_PTS - 1, THINNED_GRID_PTS) THINNED_YI_LINSPACE = numpy.linspace(0, THINNED_GRID_PTS - 1, THINNED_GRID_PTS) # Map of latitudes (north and south) to number of points on a quasi-regular (thinned) grid THINNED_GRID_PT_MAP = {0:73, 1.25:73, 2.50:73, 3.75:73, 5.0:73, 6.25:73, 7.50:73, 8.75:73, 10.0:72, 11.25:72, 12.50:72, 13.75:71, 15.0:71, 16.25:71, 17.50:70, 18.75:70, 20.0:69, 21.25:69, 22.50: 68, 23.75:67, 25.00:67, 26.25:66, 27.50:65, 28.75:65, 30.0:64, 31.25:63, 32.50:62, 33.75:61, 35.00:60, 36.25:60, 37.50:59, 38.75:58, 40.00:57, 41.25:56, 42.5:55, 43.75:54, 45.00:52, 46.25:51, 47.50:50, 48.75:49, 50.00:48, 51.25:47, 52.50:45, 53.75:44, 55.00:43, 56.25:42, 57.50:40, 58.75:39, 60.00:38, 61.25:36, 62.50:35, 63.75:33, 65.00:32, 66.25:30, 67.50:29, 68.75:28, 70.00:26, 71.25:25, 72.50:23, 73.75:22, 75.00:20, 76.25:19, 77.50:17, 78.75:16, 80.00:14, 81.25:12, 82.50:11, 83.75:9, 85.00:8, 86.25:6, 87.50:5, 88.75:3, 90.00:2} THINNED_GRID_VALUES = THINNED_GRID_PT_MAP.values() # # Python implementation of the grib decoder. This decoder uses the python ctypes # library to access the NCEP grib decoder for extracting data # # # SOFTWARE HISTORY # # Date Ticket# Engineer Description # ------------ ---------- ----------- -------------------------- # 04/7/09 #1994 bphillip Initial Creation. # Mar 25, 2013 1821 bsteffen Reshape grib data arrays in # place to improve performance. # class GribDecoder(): ## # Initializes the grib decoder # # @param text: Unused # @param filePath: The file to decode ## def __init__(self, text=None, filePath=None): # Assign public file name self.fileName = filePath ## # Decodes the grib file # # @return: List of decoded GribRecords # @rtype: List ## def decode(self): # The GribRecords to be returned back to Java records = [] #tokens = self.fileName.rsplit("/") #if tokens[len(tokens) - 1].startswith("H"): # return records filePointer = 0; version = -1; decodeFile = None if version == 1: grib1Decoder = Grib1Decoder() return grib1Decoder.decode(self.fileName) else: decodeFile = self.fileName if decodeFile == None: LogStream.logProblem("Could not get final filename to decode: [", self.fileName, "]") return records gribFile = open(decodeFile, "rb") # Define some basic navigation variable for extracting grib records recordIndex = 0 fieldIndex = 0 numFields = 1 try: # Iterate over and decode each record in the file while numFields != -1: while fieldIndex < numFields: # Extract the metadata to the metadata array metadataResults = grib2.getMetadata(gribFile, recordIndex, fieldIndex + 1, 0) numFields = metadataResults['numFields'] fieldIndex = fieldIndex + 1 if numFields != -1: metadata = metadataResults['metadata'] record = self._getData(gribFile, metadata, recordIndex, fieldIndex) if record != None: records.append(record) recordIndex = recordIndex + 1 fieldIndex = 0 except: LogStream.logProblem("Error processing file [", self.fileName, "]: ", LogStream.exc()) finally: gribFile.close() return records ## # Decodes a single record contained in the grib file # # @param fptr: The C file pointer to the file # @param metadata: The extracted metadata # @param recordIndex: The index of the record being decoded in the file # @param fieldIndex: The index of the field of the record in the file # @return: Decoded GribRecord object # @rtype: GribRecord ## def _getData(self, fptr, metadata, recordIndex, fieldIndex): # Extracts data from grib record via C call to getData dataResults = grib2.getData(fptr, recordIndex, fieldIndex) data = dataResults['data'] localSectionValues = None bitMap = None # Extracts data from the ID section idSectionValues = self._decodeIdSection(dataResults['idSection']) self.id = dataResults['idSection'] # Extracts data from the Local section if 'localSection' in dataResults: localSectionValues = self._decodeLocalSection(dataResults['localSection']) # Extracts data from the gds template gdsSectionValues = self._decodeGdsSection(metadata, dataResults['gdsTemplate']) self.gds = dataResults['gdsTemplate'] # Extracts data from the pds template pdsSectionValues = self._decodePdsSection(metadata, dataResults['idSection'], dataResults['pdsTemplate']) self.pds = dataResults['pdsTemplate'] if 'bitmap' in dataResults: bitMap = dataResults['bitmap'] # Construct the DataTime object if pdsSectionValues['endTime'] is None: dataTime = DataTime(idSectionValues['refTime'], pdsSectionValues['forecastTime']) else: # endTime defines forecast time based on the difference to refTime since forecastTime is the start of the valid period timeRange = TimeRange(idSectionValues['refTime'].getTimeInMillis() + (pdsSectionValues['forecastTime'] * 1000), pdsSectionValues['endTime'].getTimeInMillis()) forecastTime = int(float(pdsSectionValues['endTime'].getTimeInMillis() - idSectionValues['refTime'].getTimeInMillis()) / 1000) dataTime = DataTime(idSectionValues['refTime'], forecastTime, timeRange) hybridCoordList = None if 'coordList' in dataResults: hybridCoordList = numpy.resize(coordList, (1, coordList.size)) numpyDataArray = None thinnedPts = None thinnedGrid = gdsSectionValues['thinned'] # Special case for thinned grids. # Map the thinned grid on to a square lat/lon grid if thinnedGrid: optValues = dataResults['listOps'] optList = numpy.zeros(len(optValues), numpy.int32) for i in range(0, len(optValues)): optList[i] = optValues[i] dataArray = numpy.zeros(len(data), numpy.float32) for i in range(0, len(data)): dataArray[i] = data[i] # Temporary place holder pending Numpy update numpyDataArray = numpy.zeros((THINNED_GRID_PTS, THINNED_GRID_PTS), numpy.float32) # The list of points per parallel for thinned grids thinnedPts = numpy.resize(numpy.frombuffer(optList, numpy.int32)[::-1], (1, len(optList))) # Temporary arrays to hold the (grid, not lat/lon) coordinates of the data x = numpy.zeros(THINNED_GRID_SIZE) y = numpy.zeros(THINNED_GRID_SIZE) # The index in the original data dataIndex = 0 for row in range(THINNED_GRID_PTS): pts = optList[row] rowSpace = numpy.linspace(0, THINNED_GRID_PTS, pts) for curridx in range(pts): x[dataIndex] = rowSpace[curridx] y[dataIndex] = row dataIndex = dataIndex + 1 # grid the data. numpyDataArray = griddata(x, y, dataArray, THINNED_XI_LINSPACE, THINNED_YI_LINSPACE).astype(numpy.float32) # Apply the bitmap if one is provided and set masked values to missing value if metadata[18] == 0: data = numpy.where(bitMap == 0, -999999, data) # Check for fill value provided if complex packing is used drs = dataResults['drsTemplate'] if metadata[14] == 2 or metadata[14] == 3: primaryFill = Float.intBitsToFloat(drs[7]) secondaryFill = Float.intBitsToFloat(drs[8]) if drs[6] == 1: data = numpy.where(data == primaryFill, -999999, data) elif drs[6] == 2: data = numpy.where(data == primaryFill, -999999, data) data = numpy.where(data == secondaryFill, -999999, data) nx = gdsSectionValues['coverage'].getNx().intValue() ny = gdsSectionValues['coverage'].getNy().intValue() # Correct the data according to the scan mode found in the gds section. scanMode = gdsSectionValues['scanMode'] if scanMode is not None: if not thinnedGrid: numpyDataArray = numpy.reshape(data, (ny, nx)) # Check if rows are scanned in opposite direction. If so, we need to flip them around if scanMode & 16 == 16: # Check if x or y points are scanned consecutively if scanMode & 32 == 32: # y points are scanned consecutively i = 0 while i < numpyDataArray.shape[1]: theColumn = numpy.zeros(numpyDataArray.shape[0]) for j in range(0, numpyDataArray.shape[0]): theColumn[j] = numpyDataArray[j][i] for j in range(0, numpyDataArray.shape[0]): numpyDataArray[j][i] = theColumn[numpyDataArray.shape[0] - j - 1] i = i + 2 else: # x points are scanned consecutively i = 1 while i < numpyDataArray.shape[0]: theRow = numpy.array(numpyDataArray[i], copy=True) numpyDataArray[i] = theRow[::-1] i = i + 2 # Check y direction scan mode if scanMode & 64 == 64: numpyDataArray = numpy.flipud(numpyDataArray) # Check x direction scan mode if scanMode & 128 == 128: numpyDataArray = numpy.fliplr(numpyDataArray) else: if not thinnedGrid: numpyDataArray = data origCoverage = gdsSectionValues['coverage'] # check sub gridding modelName = self._createModelName(pdsSectionValues, origCoverage) spatialCache = GribSpatialCache.getInstance() gridCoverage = gdsSectionValues['coverage'] subCoverage = spatialCache.getSubGridCoverage(modelName, gridCoverage) if subCoverage is not None: subGrid = spatialCache.getSubGrid(modelName, gridCoverage) # resize the data array numpyDataArray = numpy.reshape(numpyDataArray, (ny, nx)) startx = subGrid.getUpperLeftX() starty = subGrid.getUpperLeftY() subnx = subGrid.getNX() subny = subGrid.getNY() endY = starty + subny endX = startx + subnx # handle world wide grid wrap if (endX > nx): subGridDataArray = numpy.zeros((subny, subnx), numpy.float32) midx = nx - startx subGridDataArray[0:subny, 0:midx] = numpyDataArray[starty:endY, startx:nx] subGridDataArray[0:subny, midx:subnx] = Util.GRID_FILL_VALUE numpyDataArray = subGridDataArray else: numpyDataArray = numpyDataArray[starty:endY, startx:endX] # update the number of points nx = subnx ny = subny metadata[4] = nx * ny # set the new coverage gdsSectionValues['coverage'] = subCoverage numpyDataArray = numpy.reshape(numpyDataArray, (1, metadata[4])) newAbbr = GribParamTranslator.getInstance().translateParameter(2, pdsSectionValues['parameterAbbreviation'], pdsSectionValues['centerid'], pdsSectionValues['subcenterid'], pdsSectionValues['genprocess'], dataTime, gridCoverage) if newAbbr is None: if pdsSectionValues['parameterName'] != MISSING and dataTime.getValidPeriod().getDuration() > 0: pdsSectionValues['parameterAbbreviation'] = pdsSectionValues['parameterAbbreviation'] + str(dataTime.getValidPeriod().getDuration() / 3600000) + "hr" else: pdsSectionValues['parameterAbbreviation'] = newAbbr pdsSectionValues['parameterAbbreviation'] = pdsSectionValues['parameterAbbreviation'].replace('_', '-') # Construct the GribRecord record = GridRecord() record.setPluginName(PLUGIN_NAME) record.setDataTime(dataTime) record.setMessageData(numpyDataArray) record.setLocation(gdsSectionValues['coverage']) record.setLevel(pdsSectionValues['level']) record.setDatasetId(modelName) record.addExtraAttribute("centerid", Integer(pdsSectionValues['centerid'])) record.addExtraAttribute("subcenterid", Integer(pdsSectionValues['subcenterid'])) record.addExtraAttribute("genprocess", Integer(pdsSectionValues['genprocess'])) record.addExtraAttribute("backGenprocess", Integer(pdsSectionValues['backGenprocess'])) record.addExtraAttribute("pdsTemplate", Integer(pdsSectionValues['pdsTemplateNumber'])) record.addExtraAttribute("gridid", origCoverage.getName()) if "numForecasts" in pdsSectionValues: record.addExtraAttribute("numForecasts", pdsSectionValues['numForecasts']) record.setEnsembleId(pdsSectionValues['ensembleId']) param = Parameter(pdsSectionValues['parameterAbbreviation'], pdsSectionValues['parameterName'], pdsSectionValues['parameterUnit']) GribParamTranslator.getInstance().getParameterNameAlias(modelName, param) record.setParameter(param) # record.setResCompFlags(Integer(gdsSectionValues['resCompFlags'])) #check if forecast used flag needs to be removed self._checkForecastFlag(pdsSectionValues, origCoverage, record.getDataTime()) return record ## # Decodes the values from the id section into a dictionary # @param idSectionData: The values of the ID section of the grib file # @return: A dictionary containing the values of the ID section # @rtype: dictionary ## def _decodeIdSection(self, idSectionData): # Map to hold the values idSection = {} # GRIB master tables version number (currently 2) (see table 1.0) idSection['masterTableVersion'] = idSectionData[2] # Version number of GRIB local tables used to augment Master Table (see Table 1.1) idSection['localTableVersion'] = idSectionData[3] # Significance of reference time (See table 1.2) idSection['sigRefTime'] = idSectionData[4] # The reference time as a java.util.GregorianCalendar object idSection['refTime'] = GregorianCalendar(idSectionData[5], idSectionData[6] - 1, idSectionData[7], idSectionData[8], idSectionData[9], idSectionData[10]) # Production Status of Processed Data in the GRIB message (see table 1.3) idSection['productionStatus'] = idSectionData[11] # Type of processed data in this GRIB message (See table 1.4) idSection['typeProcessedData'] = idSectionData[12] return idSection ## # Extracts the local section into a numpy array # @param localSectionData: the values of the local section of the grib file # @return: The local section as a numpy array if present, else None is returned # @rtype: numpy array else None if local section not present ## def _decodeLocalSection(self, localSectionData): # Extract the local section and resize into a numpy array if len(localSectionData) > 0: localData = numpy.zeros(len(localSectionData),numpy.int32) for i in range(0,len(localSectionData)): localData[i] = localSectionData[i] return localData # Return None if local section is not present return None ## # Decodes the values in the PDS template # # @param metadata: The metadata information # @param idSection: The ID section values # @param pdsTemplate: The PDS template values # @return: Dictionary of PDS information # @rtype: Dictionary ## def _decodePdsSection(self, metadata, idSection, pdsTemplate): # Dictionary to hold information extracted from PDS template pdsFields = {} endTime = None forecastTime = 0 duration = 0 centerID = idSection[0] subcenterID = idSection[1] pdsTemplateNumber = metadata[10] # Default to null pdsFields['ensembleId'] = None pdsFields['pdsTemplateNumber'] = pdsTemplateNumber # default to UNKNOWN pdsFields['level'] = LevelFactory.getInstance().getLevel(LevelFactory.UNKNOWN_LEVEL, float(0)); # Templates 0-11 are ordered the same for the most part and can therefore be processed the same # Exception cases are handled accordingly if pdsTemplateNumber <= 12: # Get the basic level and parameter information if (pdsTemplate[0] == 255): parameterName = MISSING parameterAbbreviation = MISSING parameterUnit = MISSING else: metadata19 = metadata[19] pds0 = pdsTemplate[0] tableName = PARAMETER_TABLE + DOT + str(metadata19) + DOT + str(pds0) parameter = GribTableLookup.getInstance().getTableValue(centerID, subcenterID, tableName, pdsTemplate[1]) if parameter is not None: parameterName = parameter.getName() if parameter.getD2dAbbrev() is not None: parameterAbbreviation = parameter.getD2dAbbrev() else: parameterAbbreviation = parameter.getAbbreviation() parameterUnit = parameter.getUnit() else: LogStream.logEvent("No parameter information for center[" + str(centerID) + "], subcenter[" + str(subcenterID) + "], tableName[" + tableName + "], parameter value[" + str(pdsTemplate[1]) + "]"); parameterName = MISSING parameterAbbreviation = MISSING parameterUnit = MISSING genprocess = GribTableLookup.getInstance().getTableValue(centerID, subcenterID, GENPROCESS_TABLE+"center"+str(centerID), pdsTemplate[4]) levelName = None; levelUnit = None; gribLevel = GribTableLookup.getInstance().getTableValue(centerID, subcenterID, LEVELS_TABLE, pdsTemplate[9]) if gribLevel is not None: levelName = gribLevel.getAbbreviation(); levelUnit = gribLevel.getUnit() else: LogStream.logEvent("No level information for center[" + str(centerID) + "], subcenter[" + str(subcenterID) + "], tableName[" + LEVELS_TABLE + "], level value[" + str(pdsTemplate[9]) + "]"); if levelName is None or len(levelName) == 0: levelName = LevelFactory.UNKNOWN_LEVEL # Convert the forecast time to seconds forecastTime = self._convertToSeconds(pdsTemplate[8], pdsTemplate[7]) # Scale the level one value if necessary if pdsTemplate[10] == 0 or pdsTemplate[11] == 0: levelOneValue = float(pdsTemplate[11]) else: levelOneValue = float(pdsTemplate[11] * pow(10, pdsTemplate[10] * - 1)) levelTwoValue = levelOneValue # If second level is present, scale if necessary if pdsTemplate[12] == 255: levelTwoValue = Level.getInvalidLevelValue() elif pdsTemplate[12] == 1: levelTwoValue = Level.getInvalidLevelValue() else: if pdsTemplate[13] == 0 or pdsTemplate[14] == 0: levelTwoValue = float(pdsTemplate[14]) else: levelTwoValue = float(pdsTemplate[14] * pow(10, pdsTemplate[13] * - 1)) if levelName=='SFC' and levelOneValue != float(0): levelOneValue=float(0) if levelName=='EATM': levelOneValue=float(0) levelTwoValue=float(Level.getInvalidLevelValue()) # Special case handling for specific PDS Templates if pdsTemplateNumber == 1 or pdsTemplateNumber == 11: typeEnsemble = Integer(pdsTemplate[15]).intValue() perturbationNumber = Integer(pdsTemplate[16]).intValue() pdsFields['numForecasts'] = Integer(pdsTemplate[17]) if(typeEnsemble == 0): pdsFields['ensembleId'] = "ctlh" + str(perturbationNumber); elif(typeEnsemble == 1): pdsFields['ensembleId'] = "ctll" + str(perturbationNumber); elif(typeEnsemble == 2): pdsFields['ensembleId'] = "n" + str(perturbationNumber); elif(typeEnsemble == 3): pdsFields['ensembleId'] = "p" + str(perturbationNumber); else: pdsFields['ensembleId'] = str(typeEnsemble) + "." + str(perturbationNumber); if pdsTemplateNumber == 11: endTime = GregorianCalendar(pdsTemplate[18], pdsTemplate[19] - 1, pdsTemplate[20], pdsTemplate[21], pdsTemplate[22], pdsTemplate[23]) numTimeRanges = pdsTemplate[24] numMissingValues = pdsTemplate[25] statisticalProcess = pdsTemplate[26] elif pdsTemplateNumber == 2 or pdsTemplateNumber == 12: derivedForecast = pdsTemplate[15] if (derivedForecast == 1 or derivedForecast == 0 ): parameterAbbreviation= parameterAbbreviation+"mean" elif (derivedForecast == 2 or derivedForecast == 3 or derivedForecast == 4 ): parameterAbbreviation= parameterAbbreviation+"sprd" pdsFields['typeEnsemble'] = Integer(pdsTemplate[15]) pdsFields['numForecasts'] = Integer(pdsTemplate[16]) if(pdsTemplateNumber == 12): endTime = GregorianCalendar(pdsTemplate[17], pdsTemplate[18] - 1, pdsTemplate[19], pdsTemplate[20], pdsTemplate[21], pdsTemplate[22]) numTimeRanges = pdsTemplate[23] numMissingValues = pdsTemplate[24] statisticalProcess = pdsTemplate[25] elif pdsTemplateNumber == 5 or pdsTemplateNumber == 9: parameterUnit = "%" probabilityNumber = pdsTemplate[15] forecastProbabilities = pdsTemplate[16] probabilityType = pdsTemplate[17] scaleFactorLL = pdsTemplate[18] scaledValueLL = pdsTemplate[19] scaleFactorUL = pdsTemplate[20] scaledValueUL = pdsTemplate[21] if(pdsTemplateNumber == 9): endTime = GregorianCalendar(pdsTemplate[22], pdsTemplate[23] - 1, pdsTemplate[24], pdsTemplate[25], pdsTemplate[26], pdsTemplate[27]) numTimeRanges = pdsTemplate[28] numMissingValues = pdsTemplate[29] statisticalProcess = pdsTemplate[30] if(probabilityType == 1 or probabilityType ==2): if(scaleFactorUL == 0): parameterAbbreviation = parameterAbbreviation+"_"+str(scaledValueUL) else: parameterAbbreviation = parameterAbbreviation+"_"+str(scaledValueUL)+"E"+str(scaleFactorUL) elif(probabilityType == 0): if(scaleFactorLL == 0): parameterAbbreviation = parameterAbbreviation+"_"+str(scaledValueLL) else: parameterAbbreviation = parameterAbbreviation+"_"+str(scaledValueLL)+"E"+str(scaleFactorLL) elif pdsTemplateNumber == 8: endTime = GregorianCalendar(pdsTemplate[15], pdsTemplate[16] - 1, pdsTemplate[17], pdsTemplate[18], pdsTemplate[19], pdsTemplate[20]) numTimeRanges = pdsTemplate[21] numMissingValues = pdsTemplate[22] statisticalProcess = pdsTemplate[23] elif pdsTemplateNumber == 10: endTime = GregorianCalendar(pdsTemplate[16], pdsTemplate[17] - 1, pdsTemplate[18], pdsTemplate[19], pdsTemplate[20], pdsTemplate[21]) numTimeRanges = pdsTemplate[22] numMissingValues = pdsTemplate[23] statisticalProcess = pdsTemplate[24] if(pdsTemplate[2] == 6 or pdsTemplate[2] == 7): parameterAbbreviation = parameterAbbreviation+"erranl" parameterAbbreviation = ParameterMapper.getInstance().lookupBaseName(parameterAbbreviation, "grib"); # Constructing the GribModel object pdsFields['centerid'] = centerID pdsFields['subcenterid'] = subcenterID pdsFields['backGenprocess'] = pdsTemplate[3] pdsFields['genprocess'] = pdsTemplate[4] pdsFields['parameterName'] = parameterName pdsFields['parameterAbbreviation'] = parameterAbbreviation pdsFields['parameterUnit'] = parameterUnit # Constructing the Level object level = LevelMapper.getInstance().lookupLevel(levelName, 'grib', levelOneValue, levelTwoValue, levelUnit) pdsFields['level'] = level # Derived forecasts based on all ensemble members at a horizontal # level or in a horizontal layer, in a continuous or non-continuous # time interval. #elif pdsTemplateNumber == 12: # pass # Derived forecasts based on a cluster of ensemble members over a # rectangular area at a horizontal level or in a horizontal layer, # in a continuous or non-continuous time interval. elif pdsTemplateNumber == 13: pass # Derived forecasts based on a cluster of ensemble members over a # circular area at a horizontal level or in a horizontal layer, in # a continuous or non-continuous time interval. elif pdsTemplateNumber == 14: pass # Radar Product elif pdsTemplateNumber == 20: pass # Satellite Product Template # NOTE:This template is deprecated. Template 31 should be used instead. elif pdsTemplateNumber == 30: pass # Satellite Product Template elif pdsTemplateNumber == 31: pass # CCITT IA5 character string elif pdsTemplateNumber == 254: pass # Cross-section of analysis and forecast at a point in time. elif pdsTemplateNumber == 1000: pass # Cross-section of averaged or otherwise statistically processed analysis or forecast over a range of time. elif pdsTemplateNumber == 1001: pass # Cross-section of analysis and forecast, averaged or otherwise statistically-processed over latitude or longitude. elif pdsTemplateNumber == 1002: pass # Hovmoller-type grid with no averaging or other statistical processing elif pdsTemplateNumber == 1100: pass # Reserved or Missing else: pass #Temporary fix to prevent invalid values getting persisted #to the database until the grib decoder is fully implemented if pdsTemplateNumber >= 13: pdsFields['parameterName'] ="Unknown" pdsFields['parameterAbbreviation'] ="Unknown" pdsFields['parameterUnit'] ="Unknown" # endtime needs to be used to calculate forecastTime and forecastTime should be used for startTime of interval pdsFields['forecastTime'] = forecastTime pdsFields['endTime'] = endTime return pdsFields ## # Decodes spatial information from the GDS template # @param metadata: The metadata information # @param gdsTemplate: The GDS Template values # @return: Dictionary of GDS information # @rtype: Dictionary ## def _decodeGdsSection(self, metadata, gdsTemplate): # Dictionary to hold information extracted from PDS template gdsFields = {} coverage = None scanMode = None resCompFlags = None thinned = False gdsTemplateNumber = metadata[7] # Latitude/Longitude projection if gdsTemplateNumber == 0: coverage = LatLonGridCoverage() majorAxis, minorAxis = self._getEarthShape(gdsTemplate) # la1 = self._correctLat(self._divideBy10e6(gdsTemplate[11])) # lo1 = self._correctLon(self._divideBy10e6(gdsTemplate[12])) # la2 = self._correctLat(self._divideBy10e6(gdsTemplate[14])) # lo2 = self._correctLon(self._divideBy10e6(gdsTemplate[15])) la1 = self._divideBy10e6(gdsTemplate[11]) lo1 = self._divideBy10e6(gdsTemplate[12]) la2 = self._divideBy10e6(gdsTemplate[14]) lo2 = self._divideBy10e6(gdsTemplate[15]) scanMode = gdsTemplate[18] resCompFlags = gdsTemplate[13] # Check for quasi-regular grid if metadata[5] > 0: # Quasi-regular grid detected thinned = True nx = THINNED_GRID_PTS ny = THINNED_GRID_PTS dx = THINNED_GRID_SPACING dy = THINNED_GRID_SPACING metadata[4] = THINNED_GRID_REMAPPED_SIZE else: # Not a quasi-regular grid nx = gdsTemplate[7] ny = gdsTemplate[8] dx = self._divideBy10e6(gdsTemplate[16]) dy = self._divideBy10e6(gdsTemplate[17]) # According to the grib2 spec 65.535 is completely valid, however it # is impossible to define anything larger than a 5x2 grid with this # spacing so we assume it is invalid and try to calculate a better # value. 65.535 was chosen because it is the value encoded in the # GFS161 model and it is completely wrong. This value is probably # an artifact of converting from grib1 to grib2 since in grib1 this # value would be encoded as an unsigned short with all bits as 1 # which is a special value in grib1, but in grib2 its just wrong if dx >= 65.535: dx = abs(lo1-lo2)/nx if dy >= 65.535: dy = abs(la1-la2)/ny coverage.setSpacingUnit(DEFAULT_SPACING_UNIT2) coverage.setNx(Integer(nx)) coverage.setNy(Integer(ny)) coverage.setLa1(la1) coverage.setLo1(lo1) coverage.setDx(dx) coverage.setDy(dy) corner = GribSpatialCache.determineFirstGridPointCorner(scanMode) coverage.setFirstGridPointCorner(corner) coverage = self._getGrid(coverage) # Rotated Latitude/Longitude projection elif gdsTemplateNumber == 1: pass # Stretched Latitude/Longitude projection elif gdsTemplateNumber == 2: pass # Rotated and Stretched Latitude/Longitude projection elif gdsTemplateNumber == 3: pass # Mercator projection elif gdsTemplateNumber == 10: coverage = MercatorGridCoverage() majorAxis, minorAxis = self._getEarthShape(gdsTemplate) nx = gdsTemplate[7] ny = gdsTemplate[8] la1 = self._correctLat(self._divideBy10e6(gdsTemplate[9])) lo1 = self._correctLon(self._divideBy10e6(gdsTemplate[10])) latin = self._correctLat(self._divideBy10e6(gdsTemplate[12])) la2 = self._correctLat(self._divideBy10e6(gdsTemplate[13])) lo2 = self._correctLon(self._divideBy10e6(gdsTemplate[14])) dx = self._divideBy10e6(gdsTemplate[17]) dy = self._divideBy10e6(gdsTemplate[18]) scanMode = gdsTemplate[15] resCompFlags = gdsTemplate[11] coverage.setSpacingUnit(DEFAULT_SPACING_UNIT) coverage.setMajorAxis(majorAxis) coverage.setMinorAxis(minorAxis) coverage.setNx(Integer(nx)) coverage.setNy(Integer(ny)) coverage.setLatin(latin) coverage.setLa1(la1) coverage.setLo1(lo1) coverage.setDx(dx) coverage.setDy(dy) corner = GribSpatialCache.determineFirstGridPointCorner(scanMode) coverage.setFirstGridPointCorner(corner) coverage = self._getGrid(coverage) # Polar Stereographic projection elif gdsTemplateNumber == 20: coverage = PolarStereoGridCoverage() majorAxis, minorAxis = self._getEarthShape(gdsTemplate) nx = gdsTemplate[7] ny = gdsTemplate[8] la1 = self._correctLat(self._divideBy10e6(gdsTemplate[9])) lo1 = self._correctLon(self._divideBy10e6(gdsTemplate[10])) lov = self._correctLon(self._divideBy10e6(gdsTemplate[13])) lad = self._correctLat(self._divideBy10e6(gdsTemplate[12])) dx = self._divideBy10e6(gdsTemplate[14]) dy = self._divideBy10e6(gdsTemplate[15]) scanMode = gdsTemplate[17] resCompFlags = gdsTemplate[11] coverage.setSpacingUnit(DEFAULT_SPACING_UNIT) coverage.setMajorAxis(majorAxis) coverage.setMinorAxis(minorAxis) coverage.setNx(Integer(nx)) coverage.setNy(Integer(ny)) coverage.setLov(lov) coverage.setLad(lad) coverage.setLa1(la1) coverage.setLo1(lo1) coverage.setDx(dx) coverage.setDy(dy) corner = GribSpatialCache.determineFirstGridPointCorner(scanMode) coverage.setFirstGridPointCorner(corner) coverage = self._getGrid(coverage) # Lambert Conformal projection elif gdsTemplateNumber == 30: coverage = LambertConformalGridCoverage() majorAxis, minorAxis = self._getEarthShape(gdsTemplate) nx = gdsTemplate[7] ny = gdsTemplate[8] la1 = self._correctLat(self._divideBy10e6(gdsTemplate[9])) lo1 = self._correctLon(self._divideBy10e6(gdsTemplate[10])) lov = self._correctLon(self._divideBy10e6(gdsTemplate[13])) dx = self._divideBy10e6(gdsTemplate[14]) dy = self._divideBy10e6(gdsTemplate[15]) latin1 = self._correctLat(self._divideBy10e6(gdsTemplate[18])) latin2 = self._correctLat(self._divideBy10e6(gdsTemplate[19])) scanMode = gdsTemplate[17] resCompFlags = gdsTemplate[11] coverage.setSpacingUnit(DEFAULT_SPACING_UNIT) coverage.setMajorAxis(majorAxis) coverage.setMinorAxis(minorAxis) coverage.setNx(Integer(nx)) coverage.setNy(Integer(ny)) coverage.setLov(lov) coverage.setLa1(la1) coverage.setLo1(lo1) coverage.setDx(dx) coverage.setDy(dy) coverage.setLatin1(latin1) coverage.setLatin2(latin2) corner = GribSpatialCache.determineFirstGridPointCorner(scanMode) coverage.setFirstGridPointCorner(corner) coverage = self._getGrid(coverage) # Albers Equal Area projection elif gdsTemplate == 31: pass # Gaussian Latitude/Longitude projection elif gdsTemplate == 40: pass # Rotated Gaussian Latitude/Longitude projection elif gdsTemplate == 41: pass # Stretched Gaussian Latitude/Longitude projection elif gdsTemplate == 42: pass # Rotated and Stretched Gaussian Latitude/Longitude projection elif gdsTemplate == 43: pass # Spherical Harmonic Coefficients elif gdsTemplate == 50: pass # Rotated Spherical Harmonic Coefficients elif gdsTemplate == 51: pass # Stretched Spherical Harmonic Coefficients elif gdsTemplate == 52: pass # Rotated and Stretched Spherical Harmonic Coefficients elif gdsTemplate == 53: pass # Space View Perspective or Orthographic elif gdsTemplate == 90: pass # Triangular Grid based on Icosahedron elif gdsTemplate == 100: pass # Equatorial Azimuthal Equidistance projection elif gdsTemplate == 110: pass # Azimuth-Range projection elif gdsTemplate == 120: pass # Curvilinear Orthogonal projection elif gdsTemplate == 204: pass # Cross Section Grid with Points Equally spaced on the horizontal elif gdsTemplate == 1000: pass # Hovmoller Diagram with Points Equally spaced on the horizontal elif gdsTemplate == 1100: pass # Time Section grid elif gdsTemplate == 1200: pass # Rotated Latitude/Longitude (Arakawa Staggered E-Grid) elif gdsTemplate == 32768: pass # Missing elif gdsTemplate == 65535: pass gdsFields['scanMode'] = scanMode gdsFields['coverage'] = coverage gdsFields['thinned'] = thinned gdsFields['resCompFlags'] = resCompFlags return gdsFields ## # Gets a grid from the cache. If not found, one is created and stored to the cache # # @param temp: A GridCoverage object withough geometry or crs information populated # @return: A GribCoverage object # @rtype: GribCoverage ## def _getGrid(self, temp): # Check the cache first grid = GribSpatialCache.getInstance().getGrid(temp) # If not found, create a new GribCoverage and store in the cache if grid is None: grid = GridCoverageLookup.getInstance().getCoverage(temp, True) return grid ## # Divides a number by 1000 # # @param number: A number to be divided by 1000 # @return: The provided number divided by 1000 # @rtype: float ## def _divideBy10e3(self, number): return float(float(number) / 1000) ## # Divides a number by 1000000 # # @param number: A number to be divided by 1000000 # @return: The provided number divided by 1000000 # @rtype: float ## def _divideBy10e6(self, number): return float(float(number) / 1000000) ## # Convert a scaledValue and scaleFactor to the unscaled value # # @param scaledValue: The scaled value # @param scaleFactor: The scale factor # @return: The unscaled value # @rtype: float ## def _convertScaledValue(self, scaledValue, scaleFactor): return float(scaledValue) / 10**scaleFactor ## # Corrects a longitude to fall within the geotools required bounds of -180 and 180 # # @param lon: The longitude to be corrected # @return: The corrected longitude # @rtype: float ## def _correctLon(self, lon): if lon < 0: lon = lon % 360 else: lon = lon % 360 if lon > 180: lon = (180 - lon % 180) * - 1 elif lon < - 180: lon = (180 - (- lon % 180)) return lon ## # Corrects a latitude to fall within the geotools required bounds of -90 and 90 # # @param lat: The latitude to be corrected # @return: The corrected latitude # @rtype: float ## def _correctLat(self, lat): if lat < 0: lat = lat % -180 else: lat = lat % 180 if lat > 90: lat = 90 - lat % 90 elif lat < - 90: lat = (90 - (- lat % 90)) * - 1 return lat ## # Gets the shape of the earth based on Table 3.2 # # @param gdsTemplate:The gdsTemplate values # @return: The minor and major axis sizes of the earth # @rtype: long, long ## def _getEarthShape(self, gdsTemplate): # Shape of the earth which keys into Table 3.2 number = gdsTemplate[0] # # Determine the shape of Earth based on Table 3.2 # # Earth assumed spherical with radius = 6,367,470.0 m if number == 0: minorAxis = 6367470.0 majorAxis = 6367470.0 # Earth assumed spherical with radius specified (in m) by data producer elif number == 1: minorAxis = self._convertScaledValue(gdsTemplate[2], gdsTemplate[1]) majorAxis = minorAxis if majorAxis < 6000000.0 or minorAxis < 6000000.0: LogStream.logEvent("Invalid earth shape majorAxis,minorAxis = " + str(majorAxis) + "," + str(minorAxis) + " defaulting to 6367470.0,6367470.0") minorAxis = majorAxis = 6367470.0 # Earth assumed oblate spheriod with size as determined by IAU in 1965 # (major axis = 6,378,160.0 m, minor axis = 6,356,775.0 m, f = 1/297.0) elif number == 2: minorAxis = 6356775.0 majorAxis = 6378160.0 # Earth assumed oblate spheriod with major and minor axes specified (in km) by data producer elif number == 3: minorAxis = self._convertScaledValue(gdsTemplate[4], gdsTemplate[3]) * 1000 if minorAxis < 6000000.0: LogStream.logEvent("Invalid earth shape minorAxis = " + str(minorAxis) + " defaulting to " + MINOR_AXIS_DEFAULT) minorAxis = MINOR_AXIS_DEFAULT majorAxis = self._convertScaledValue(gdsTemplate[6], gdsTemplate[5]) * 1000 if majorAxis < 6000000.0: LogStream.logEvent("Invalid earth shape majorAxis = " + str(majorAxis) + " defaulting to " + MAJOR_AXIS_DEFAULT) majorAxis = MAJOR_AXIS_DEFAULT # Earth assumed oblate spheriod as defined in IAG-GRS80 model # (major axis = 6,378,137.0 m, minor axis = 6,356,752.314 m, f = 1/298.257222101) elif number == 4: minorAxis = 6356752.314 majorAxis = 6378137.0 # Earth assumed represented by WGS84 (as used by ICAO since 1998) elif number == 5: minorAxis = 6356752.314245 majorAxis = 6378137.0 # Earth assumed spherical with radius = 6,371,229.0 m elif number == 6: minorAxis = 6371229.0 majorAxis = 6371229.0 # Earth assumed oblate spheroid with major and minor axes specified (in m) by data producer elif number == 7: minorAxis = self._convertScaledValue(gdsTemplate[4], gdsTemplate[3]) if minorAxis < 6000000.0: LogStream.logEvent("Invalid earth shape minorAxis = " + str(minorAxis) + " defaulting to " + MINOR_AXIS_DEFAULT) minorAxis = MINOR_AXIS_DEFAULT majorAxis = self._convertScaledValue(gdsTemplate[6], gdsTemplate[5]) if majorAxis < 6000000.0: LogStream.logEvent("Invalid earth shape majorAxis = " + str(majorAxis) + " defaulting to " + MAJOR_AXIS_DEFAULT) majorAxis = MAJOR_AXIS_DEFAULT # Earth model assumed spherical with radius 6,371,200 m, # but the horizontal datum of the resulting Latitude/Longitude field is # the WGS84 reference frame elif number == 8: minorAxis = 6371200.0 majorAxis = 6371200.0 else: minorAxis = MINOR_AXIS_DEFAULT majorAxis = MAJOR_AXIS_DEFAULT return float(majorAxis), float(minorAxis) ## # Converts a value in the specified unit (according to table 4.4) to seconds # # @param value: The value to convert to seconds # @param fromUnit: The value from Table 4.4 to convert from # @return: The number of seconds of the provided value # @rtype: long ## def _convertToSeconds(self, value, fromUnit): retVal = value # Convert from minutes if fromUnit == 0: retVal = value * SECONDS_PER_MINUTE # Convert from hours elif fromUnit == 1: retVal = value * SECONDS_PER_HOUR # Convert from days elif fromUnit == 2: retVal = value * SECONDS_PER_DAY # Convert from months elif fromUnit == 3: retVal = value * SECONDS_PER_MONTH # Convert from years elif fromUnit == 4: retVal = value * SECONDS_PER_YEAR # Convert from decades elif fromUnit == 5: retVal = value * 10 * SECONDS_PER_YEAR # Convert from Normal (30 years) elif fromUnit == 6: retVal = value * 30 * SECONDS_PER_YEAR # Convert from centuries elif fromUnit == 7: retVal = value * 100 * SECONDS_PER_YEAR # Convert from 3 hours elif fromUnit == 10: retVal = value * 3 * SECONDS_PER_HOUR # Convert from 6 hours elif fromUnit == 11: retVal = value * 6 * SECONDS_PER_HOUR # Convert from 12 horus elif fromUnit == 12: retVal = value * 12 * SECONDS_PER_HOUR return retVal def _getGridModel(self, pdsSectionValues, grid): center = pdsSectionValues['centerid'] subcenter = pdsSectionValues['subcenterid'] process = pdsSectionValues['genprocess'] gridModel = GribModelLookup.getInstance().getModel(center, subcenter, grid, process) return gridModel def _createModelName(self, pdsSectionValues, grid): center = pdsSectionValues['centerid'] subcenter = pdsSectionValues['subcenterid'] process = pdsSectionValues['genprocess'] return GribModelLookup.getInstance().getModelName(center, subcenter, grid, process) def _checkForecastFlag(self, pdsSectionValues, grid, dataTime): gridModel = self._getGridModel(pdsSectionValues, grid) if gridModel is None: return else: if gridModel.getAnalysisOnly(): dataTime.getUtilityFlags().remove(FLAG.FCST_USED)