xan/awips2
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
awips2/ncep/gov.noaa.nws.ncep.edex.plugin.ncgrib/NcgribDecoder.py
Richard Peter 8c9e5f214b Issue #2298: Remove setPluginName python calls. 
Former-commit-id: 6b93f9ad06 [formerly b666f99d3dc79b90525177072e2448077c39545a]
Former-commit-id: 799fe37bca
2013-09-12 12:49:23 -05:00

1602 lines
No EOL
64 KiB
Python
Raw Blame History

##
# 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 gov.noaa.nws.ncep.common.dataplugin.ncgrib import NcgribRecord
from gov.noaa.nws.ncep.common.dataplugin.ncgrib import NcgribModel
from gov.noaa.nws.ncep.common.dataplugin.ncgrib import NcgribParameter
from gov.noaa.nws.ncep.common.dataplugin.ncgrib import Ncgrib1Parameter
from gov.noaa.nws.ncep.common.dataplugin.ncgrib.util import Ncgrib1ParameterLookup
from gov.noaa.nws.ncep.common.dataplugin.ncgrib.spatial.projections import LambertConformalNcgridCoverage
from gov.noaa.nws.ncep.common.dataplugin.ncgrib.spatial.projections import LatLonNcgridCoverage
from gov.noaa.nws.ncep.common.dataplugin.ncgrib.spatial.projections import MercatorNcgridCoverage
from gov.noaa.nws.ncep.common.dataplugin.ncgrib.spatial.projections import PolarStereoNcgridCoverage
from com.raytheon.uf.common.dataplugin.level import Level
from com.raytheon.uf.common.dataplugin.level import LevelFactory
from gov.noaa.nws.ncep.edex.plugin.ncgrib.spatial import NcgribSpatialCache
from gov.noaa.nws.ncep.edex.plugin.ncgrib.util import NcgribModelCache
from gov.noaa.nws.ncep.edex.util.ncgrib import NcgribTableLookup
from gov.noaa.nws.ncep.edex.util.ncgrib import NcgribModelLookup
#from gov.noaa.nws.ncep.common.dataplugin.ncgrib.util import NcgribModelLookup
from gov.noaa.nws.ncep.edex.util.grib2vars import Grib2VarsTableLookup
from gov.noaa.nws.ncep.edex.util.grib2vcrd import Grib2VcrdTableLookup
from gov.noaa.nws.ncep.edex.plugin.ncgrib import Ncgrib1Decoder
from gov.noaa.nws.ncep.edex.util.ncgrib import NcgribParamTranslator
# 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 ncgrib 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.
# 10/13/10 #276 llin modified for NC GRIB.
# 06/28/11 xguo Added codes to get correct VAR/VCD
# 07/12/11 xguo Changed Grib1Decoder() to Ncgrib1Decoder()
# 07/26/11 xguo Added codes to handle Derived Ensemble Data
# 09/08/11 xguo Added one column grib data
# 09/21/11 xguo Check derived ensemble data
# 10/04/11 xguo Remove '_' from model name
# 11/02/11 xguo Added codes to decode firewx
# 11/08/11 xguo Adjusted glevel1/glevel2 for PRES/PDLY/POTV
# 11/22/11 xguo Updated Level infor in model
# Sep 04, 2013 2298 rjpeter Removed setPluginName call
#
class NcgribDecoder():
##
# Initializes the ncgrib decoder
#
# @param text: Unused
# @param filePath: The file to decode
##
def __init__(self, text=None, filePath=None):
# Assign public file name
self.fileName = filePath
self.dicipline = -1
self.inputFile = None
self.abbr = None
self.derived = None
self.VCD = -1
self.count = 1
self.AddColumn = -1
##
# Decodes the ncgrib file
#
# @return: List of decoded NcgribRecords
# @rtype: List
##
def decode(self):
# The NcgribRecords to be returned back to Java
records = []
filePointer = 0;
version = -1;
if os.path.exists(self.fileName):
try:
version = grib2.checkVersion(self.fileName)
except:
LogStream.logProblem("Error opening file [", self.fileName, "]: ", sys.exc_info()[1])
return records
else:
LogStream.logProblem("The file does not exist: [", self.fileName, "]")
return records
decodeFile = None
if version == 1:
grib1Decoder = Ncgrib1Decoder()
return grib1Decoder.decode(self.fileName)
else:
decodeFile = self.fileName
fileList = decodeFile.split("/")
self.inputFile = fileList[len(fileList) - 1]
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:
numFieldsSave = numFields
metadata = metadataResults['metadata']
record = self._getData(gribFile, metadata, recordIndex, fieldIndex)
if record != None:
self._addRecord(records, record)
if self.abbr == "uW" or self.abbr == "uWmean" or self.abbr == "uWsprd" or self.abbr == "uWprob":
#Extract a second field if it exists
metadataResults = grib2.getMetadata(gribFile, recordIndex, fieldIndex + 1, 0)
numFields = metadataResults['numFields']
fieldIndex = fieldIndex + 1
if numFields != - 2:
metadata = metadataResults['metadata']
record = self._getData(gribFile, metadata, recordIndex, fieldIndex)
if record != None:
self._addRecord(records, record)
numFields = numFieldsSave
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 NcgribRecord object
# @rtype: NcgribRecord
##
def _getData(self, fptr, metadata, recordIndex, fieldIndex):
interval = -9999
# get discipline
self.dicipline = metadata[19]
# Extracts data from ncgrib 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']
self.VCD = self.pds[9]
lenPds = len(self.pds)
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)
nx = gdsSectionValues['coverage'].getNx().intValue()
ny = gdsSectionValues['coverage'].getNy().intValue()
if self.AddColumn == 1:
nx = nx - 1
# Correct the data according to the scan mode found in the gds section.
scanMode = gdsSectionValues['scanMode']
if scanMode is not None:
# Flip the grid vertically
if scanMode == 64:
if not thinnedGrid:
numpyDataArray = numpy.resize(data, (ny, nx))
numpyDataArray = numpy.flipud(numpyDataArray)
# Flip the grid horizontally
elif scanMode == 128:
if not thinnedGrid:
numpyDataArray = numpy.resize(data, (ny, nx))
numpyDataArray = numpy.fliplr(numpyDataArray)
# Flip vertically and horizontally (rotate 180 degrees)
elif scanMode == 192:
if not thinnedGrid:
numpyDataArray = numpy.resize(data, (ny, nx))
numpyDataArray = numpy.rot90(numpyDataArray, 2)
# No modification necessary
else:
if not thinnedGrid:
numpyDataArray = data
else:
if not thinnedGrid:
numpyDataArray = data
# check sub gridding
#
# Synchronizes model information with the database.
#
pdsSectionValues['model'].setGridid(gdsSectionValues['coverage'].getName())
self._createModelName(pdsSectionValues['model'], self.inputFile)
modelName = pdsSectionValues['model'].getModelName()
if modelName == "ghmNest" or modelName == "ghm6th" or \
modelName == "hwrfNest" or modelName == "gfs" or \
modelName == "fireWxNest" or modelName == "hwrf":
pdsSectionValues['model'].generateId(self.inputFile)
else:
pdsSectionValues['model'].generateId()
spatialCache = NcgribSpatialCache.getInstance()
subCoverage = spatialCache.getSubGridCoverage(modelName)
if subCoverage is not None:
subGrid = spatialCache.getSubGrid(modelName)
# resize the data array
numpyDataArray = numpy.resize(numpyDataArray, (ny, nx))
startx = subGrid.getStartX()
starty = subGrid.getStartY()
nx = subGrid.getNX()
ny = subGrid.getNY()
numpyDataArray = numpyDataArray[starty:starty + ny,startx:startx + nx]
# update the number of points
metadata[4] = nx * ny
# set the new coverage
gdsSectionValues['coverage'] = subCoverage
numpyDataArray = numpy.resize(numpyDataArray, (1, metadata[4]))
pdsSectionValues['model'].setLocation(gdsSectionValues['coverage'])
if pdsSectionValues['model'].getParameterName() == MISSING:
model = pdsSectionValues['model']
else:
model = NcgribModelCache.getInstance().getModel(pdsSectionValues['model'])
# Construct the DataTime object
if pdsSectionValues['endTime'] is None:
dataTime = DataTime(idSectionValues['refTime'], pdsSectionValues['forecastTime'])
duration = dataTime.getValidPeriod().getDuration()
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)
duration = dataTime.getValidPeriod().getDuration()
newAbbr = NcgribParamTranslator.getInstance().translateParameter(2,pdsSectionValues['model'],dataTime)
if newAbbr is None:
if pdsSectionValues['model'].getParameterName() != MISSING and duration > 0:
abbrStr = pdsSectionValues['model'].getParameterAbbreviation() + str(duration / 3600000) + "hr"
# pdsSectionValues['model'].setParameterAbbreviation(abbrStr)
model.setParameterAbbreviation(abbrStr)
else:
# pdsSectionValues['model'].setParameterAbbreviation(newAbbr)
model.setParameterAbbreviation(newAbbr)
# self.abbr = pdsSectionValues['model'].getParameterAbbreviation()
self.abbr = model.getParameterAbbreviation()
# if pdsSectionValues['model'].getParameterName() == MISSING:
# model = pdsSectionValues['model']
# else:
# model = NcgribModelCache.getInstance().getModel(pdsSectionValues['model'])
# Construct the ncgribRecord
record = NcgribRecord()
record.setDataTime(dataTime)
record.setMasterTableVersion(idSectionValues['masterTableVersion'])
record.setLocalTableVersion(idSectionValues['localTableVersion'])
record.setRefTimeSignificance(idSectionValues['sigRefTime'])
#record.setProductionStatus(idSectionValues['productionStatus'])
record.setProcessedDataType(idSectionValues['typeProcessedData'])
record.setLocalSectionUsed(localSectionValues is not None)
record.setLocalSection(localSectionValues)
record.setHybridGrid(hybridCoordList is not None)
record.setHybridCoordList(hybridCoordList)
record.setThinnedGrid(False)
record.setThinnedPts(thinnedPts)
record.setModelName(model.getModelName())
record.setFileName(self.inputFile)
tokens = self.inputFile.split(".")
if len(tokens) >= 3 and tokens[2] == "firewxnest":
record.setEventName(tokens[2]);
else :
record.setEventName(tokens[0])
record.setModelInfo(model)
record.setResCompFlags(Integer(gdsSectionValues['resCompFlags']))
discipline = self.dicipline;
record.setDiscipline(discipline)
record.setCategory(int(pdsSectionValues['category']))
record.setParameterId(int(pdsSectionValues['parameterId']))
pdsTemplateNumber = metadata[10]
# pdt = metadata[10]
record.setPdt(pdsTemplateNumber)
# In order match press unit (mb) in GEMPACK, do this scale,
# expect to remove the scale in near future
category = record.getCategory()
parameterId = record.getParameterId()
discipline = record.getDiscipline()
# pdt = record.getPdt()
g2scale = Grib2VarsTableLookup.getG2varsScale(discipline, category, parameterId, pdsTemplateNumber)
if g2scale != 0:
scale = pow(10,g2scale )
missscale = -999999 * scale
for i in range(0,len(numpyDataArray)):
numpyDataArray[i] = numpyDataArray[i] * scale
numpyDataArray = numpy.where(numpyDataArray == missscale, -999999, numpyDataArray)
if self.AddColumn == 1:
self.AddColumn = -1
numpyDataArray1 = None
nx1 = nx + 1
metadata[4] = nx1 * ny
numpyDataArray1 = numpy.zeros((ny, nx1), numpy.float32)
numpyDataArray = numpy.resize(numpyDataArray, (ny, nx))
for i in range(ny):
for j in range(nx):
numpyDataArray1[i,j] = numpyDataArray[i,j]
numpyDataArray1[i,nx] = numpyDataArray[i,0]
numpyDataArray1 = numpy.resize(numpyDataArray1, (1, metadata[4]))
record.setMessageData(numpyDataArray1)
else:
record.setMessageData(numpyDataArray)
if pdsTemplateNumber == 8 :
if lenPds > 27 :
interval = self.pds[26]
record.setInterval(interval)
record.setProcessType(int(pdsSectionValues['processedDataType']))
record.setVcrdId1(int(pdsSectionValues['verticalCoordinate1']))
record.setVcrdId2(int(pdsSectionValues['verticalCoordinate2']))
record.setDecodedLevel1( pdsSectionValues['level1'] )
record.setGridVersion(2)
if pdsSectionValues['level2'] is not None:
record.setDecodedLevel2( pdsSectionValues['level2'] )
# Special case handling for ffg grids
if(model.getCenterid() == 9):
from com.raytheon.edex.plugin import PluginFactory
record.constructDataURI()
ncgribDao = PluginFactory.getInstance().getPluginDao("ncgrib")
result = ncgribDao.executeNativeSql("select max(gridVersion) from awips.grib where datauri like '" +
record.getDataURI() + "%'")
resultCount = result.getResultCount()
if(resultCount == 1 and result.getRowColumnValue(0, 0) is not None):
newVersion = result.getRowColumnValue(0, 0).intValue() + 1
record.setGridVersion(newVersion)
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:
return numpy.resize(numpy.frombuffer(localSectionData, numpy.float32), (1, len(localSectionData)))
# 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 = {}
model = NcgribModel()
endTime = None
forecastTime = 0
duration = 0
centerID = idSection[0]
subcenterID = idSection[1]
pdsTemplateNumber = metadata[10]
model.setPdsTemplate(pdsTemplateNumber)
# default to UNKNOWN
model.setLevel(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 = NcgribTableLookup.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 = NcgribTableLookup.getInstance().getTableValue(centerID, subcenterID, GENPROCESS_TABLE+"center"+str(centerID), pdsTemplate[4])
levelName = None;
levelUnit = None;
ncgribLevel = NcgribTableLookup.getInstance().getTableValue(centerID, subcenterID, LEVELS_TABLE, pdsTemplate[9])
if ncgribLevel is not None:
levelName = ncgribLevel.getAbbreviation();
levelUnit = ncgribLevel.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 levelName=='SFC' and levelOneValue != float(0):
levelOneValue=float(0)
# 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))
# Special case handling for specific PDS Templates
if pdsTemplateNumber == 1 or pdsTemplateNumber == 11:
pdst15 = pdsTemplate[15]
model.setTypeEnsemble(Integer(pdst15))
# Use the following codes with correct grib headers
# print "Type of Ensemble Forecast: ", pdst15, " perturbation number: ", pdsTemplate[16]
#if ( pdst15 == 0 or self.fileName.find ('ctl1') != -1 ):
# self.derived = 'ctl1'
#elif ( pdst15 == 1 or self.fileName.find ('ctl2') != -1 ):
# self.derived = 'ctl2'
#elif pdst15 == 2:
# self.derived = 'n'
# to do
#elif pdst15 == 3:
# self.derived = 'p'
# To do
# set perturbation number when code value = 2/3/192 to avoid 0
if pdst15 == 2 or pdst15 == 3 or pdst15 == 192:
model.setPerturbationNumber(str(pdsTemplate[16]))
model.setNumForecasts(Integer(pdsTemplate[17]))
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 == 0 or derivedForecast == 1 or derivedForecast == 6):
# parameterAbbreviation= parameterAbbreviation+"mean"
self.derived = 'mean'
elif (derivedForecast == 2 or derivedForecast == 3 ):
# parameterAbbreviation= parameterAbbreviation+"sprd"
self.derived = 'sprd'
elif (derivedForecast >= 193 and derivedForecast <= 195 ):
# parameterAbbreviation= parameterAbbreviation+"prob"
self.derived = 'prob'
# Use the following codes with correct grib headers
# if (derivedForecast >= 192 and derivedForecast <= 195 ):
#if ( derivedForecast == 193 or self.fileName.find ('10p') != -1):
# self.derived = '10p'
#elif ( derivedForecast == 194 or self.fileName.find ('50p') != -1):
# self.derived = '50p'
#elif ( derivedForecast == 195 or self.fileName.find ('90p') != -1):
# self.derived = '90p'
#elif ( derivedForecast == 192 or self.fileName.find ('mode') != -1):
# self.derived = 'mode'
model.setTypeEnsemble(Integer(pdsTemplate[15]))
model.setNumForecasts(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:
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]
# Constructing the NcgribModel object
model.setCenterid(centerID)
model.setSubcenterid(subcenterID)
model.setBackGenprocess(pdsTemplate[3])
model.setGenprocess(pdsTemplate[4])
model.setParameterName(parameterName)
model.setParameterAbbreviation(parameterAbbreviation)
model.setParameterUnit(parameterUnit)
tokens = self.inputFile.split(".")
if len(tokens) >= 3 and tokens[2] == "firewxnest":
model.setEventName(tokens[2])
else:
model.setEventName(tokens[0])
# Constructing the Level object
level = LevelFactory.getInstance().getLevel(levelName, levelOneValue, levelTwoValue, levelUnit)
model.setLevel(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:
model.setParameterName("Unknown")
model.setParameterAbbreviation("Unknown")
model.setParameterUnit("Unknown")
# endtime needs to be used to calculate forecastTime and forecastTime should be used for startTime of interval
pdsFields['forecastTime'] = forecastTime
pdsFields['endTime'] = endTime
pdsFields['model'] = model
pdsFields['level1'] = levelOneValue #pdsTemplate[11]
pdsFields['level2'] = levelTwoValue #pdsTemplate[14]
pdsFields['category'] = pdsTemplate[0]
pdsFields['parameterId'] = pdsTemplate[1]
pdsFields['processedDataType'] = pdsTemplate[2]
pdsFields['verticalCoordinate1'] = pdsTemplate[9]
pdsFields['verticalCoordinate2'] = pdsTemplate[12]
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 = LatLonNcgridCoverage()
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])
Lon1 = lo1 + lo2
Lon2 = lo2 + dx
if Lon1 == 360.0 or Lon2 == 360.0:
nx = nx + 1
lo2 = lo2 + dx
self.AddColumn = 1
coverage.setSpacingUnit(DEFAULT_SPACING_UNIT2)
coverage.setNx(Integer(nx))
coverage.setNy(Integer(ny))
coverage.setLa1(la1)
coverage.setLo1(lo1)
coverage.setLa2(la2)
coverage.setLo2(lo2)
coverage.setDx(dx)
coverage.setDy(dy)
coverage.setId(coverage.hashCode())
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 = MercatorNcgridCoverage()
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]
Lon1 = lo1 + lo2
Lon2 = lo2 + dx
if Lon1 == 360.0 or Lon2 == 360.0:
nx = nx + 1
lo2 = lo2 + dx
self.AddColumn = 1
coverage.setSpacingUnit(DEFAULT_SPACING_UNIT)
coverage.setMajorAxis(majorAxis)
coverage.setMinorAxis(minorAxis)
coverage.setNx(Integer(nx))
coverage.setNy(Integer(ny))
coverage.setLa1(la1)
coverage.setLo1(lo1)
coverage.setLatin(latin)
coverage.setLa2(la2)
coverage.setLo2(lo2)
coverage.setDx(dx)
coverage.setDy(dy)
coverage.setId(coverage.hashCode())
coverage = self._getGrid(coverage)
# Polar Stereographic projection
elif gdsTemplateNumber == 20:
coverage = PolarStereoNcgridCoverage()
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])
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.setLa1(la1)
coverage.setLo1(lo1)
coverage.setLov(lov)
coverage.setDx(dx)
coverage.setDy(dy)
coverage.setId(coverage.hashCode())
coverage = self._getGrid(coverage)
# Lambert Conformal projection
elif gdsTemplateNumber == 30:
coverage = LambertConformalNcgridCoverage()
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.setLa1(la1)
coverage.setLo1(lo1)
coverage.setLov(lov)
coverage.setDx(dx)
coverage.setDy(dy)
coverage.setLatin1(latin1)
coverage.setLatin2(latin2)
coverage.setId(coverage.hashCode())
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 str(gdsTemplateNumber) == "32768":
coverage = LatLonNcgridCoverage()
majorAxis, minorAxis = self._getEarthShape(gdsTemplate)
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])
Lon1 = lo1 + lo2
Lon2 = lo2 + dx
if Lon1 == 360.0 or Lon2 == 360.0:
nx = nx + 1
lo2 = lo2 + dx
self.AddColumn = 1
coverage.setSpacingUnit(DEFAULT_SPACING_UNIT2)
coverage.setNx(Integer(nx))
coverage.setNy(Integer(ny))
coverage.setLa1(la1)
coverage.setLo1(lo1)
coverage.setLa2(la2)
coverage.setLo2(lo2)
coverage.setDx(dx)
coverage.setDy(dy)
coverage.setId(coverage.hashCode())
coverage = self._getGrid(coverage)
# 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):
cache = NcgribSpatialCache.getInstance()
# Check the cache first
grid = cache.getGrid(temp)
# If not found, create a new GribCoverage and store in the cache
if grid is None:
cache.putGrid(temp, True)
grid = cache.getGrid(temp.getId())
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)
##
# 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 = float(gdsTemplate[2])
majorAxis = float(gdsTemplate[2])
# 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:
if gdsTemplate[3] > 0:
minorAxis = gdsTemplate[4] * gdsTemplate[3] * 1000
else:
minorAxis = gdsTemplate[4] * 1000
if gdsTemplate[5] > 0:
majorAxis = gdsTemplate[6] * gdsTemplate[5] * 1000
else:
majorAxis = gdsTemplate[6] * 1000
# 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:
if gdsTemplate[3] > 0:
minorAxis = gdsTemplate[4] * gdsTemplate[3]
else:
minorAxis = gdsTemplate[4]
if gdsTemplate[5] > 0:
majorAxis = gdsTemplate[6] * gdsTemplate[5]
else:
majorAxis = gdsTemplate[6]
# 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 majorAxis, 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 _createModelName(self, model, fileName):
center = model.getCenterid()
subcenter = model.getSubcenterid()
gridid = model.getGridid()
process = model.getGenprocess()
gridModel = NcgribModelLookup.getInstance().getModel(center, subcenter, gridid, process, fileName, model)
if gridModel is None:
name = "NewGrid:" + str(center) + ":" + str(subcenter) + ":" + str(process) + ":" + gridid
tokens = fileName.split(".")
hurricane = tokens[0]
basin = hurricane[-1]
trackno = hurricane[-3:-1]
if trackno.isdigit() or tokens[2] =="firewxnest" or tokens[2] == "hysplit":
if trackno.isdigit():
basins = "lewcs"
if basin in basins:
#name = "GHM:" + str(center) + ":" + str(subcenter) + ":" + str(process) + ":" + gridid
#hurricaneName = hurricane[:len(hurricane)-3]
name = "ghm"
if tokens[2] == "gribn3":
name = "ghmNest"
elif tokens[2] == "grib6th":
name = "ghm6th"
elif tokens[2] == "hwrfprs_n":
name = "hwrfNest"
elif tokens[2] == "hwrfprs_p":
name = "hwrf"
elif tokens[2] == "firewxnest":
name = "fireWxNest"
else:
name = "hysplit"
NcgribModelLookup.getInstance().setModel(center, subcenter, gridid, process, name)
gridModel = NcgribModelLookup.getInstance().getModel(center, subcenter, gridid, process, fileName,model)
#name = gridModel.getName()
else:
name = gridModel.getName()
model.setModelName(name)
##
# Add record to record list if its parameters are in
# NCEP VARS and VCRD control files.
#
# @param records: the record list.
# @param record: the record is going to be checked.
##
def _addRecord(self, records, record):
category = record.getCategory()
discipline = record.getDiscipline()
pdt = record.getPdt()
parameterId = record.getParameterId()
vcrd1 = record.getVcrdId1()
vcrd2 = record.getVcrdId2()
g2varsId = -1
g2vcrdId = -1
vcord=''
scale=''
parm=''
accumHour=0
charHour=''
#calculate accumulative hour
interval = record.getInterval()
if pdt == 8 :
accumHour = interval
if accumHour < 12 :
charHour = '0' + str(accumHour)
else :
charHour = str(accumHour)
g2vcrdId = Grib2VcrdTableLookup.getG2vcrdId(vcrd1, vcrd2)
if str(g2vcrdId).isdigit() :
# get and set vcord and scale
vcord = Grib2VcrdTableLookup.getGnamByVcrdId(vcrd1, vcrd2)
record.setVcord(vcord)
scale = float( Grib2VcrdTableLookup.getScaleByVcrdId(vcrd1, vcrd2) )
scale = pow ( 10, scale)
#record.setScale(scale)
glevel1 = record.getDecodedLevel1()
if (glevel1 != 0 and scale != 1 ) :
record.setGlevel1(int(round(glevel1 * scale)))
else :
record.setGlevel1(int(round(glevel1)))
if vcrd2 == 255 :
record.setGlevel2(-9999)
else :
glevel2 = record.getDecodedLevel2()
if (glevel2 != 0 and scale != 1 ) :
record.setGlevel2(int(round(glevel2 * scale)))
else :
record.setGlevel2( int(round(glevel2) ))
levelName = record.getModelInfo().getLevelName()
levelOneValue = 1.0*record.getGlevel1()
levelTwoValue = 1.0*record.getGlevel2()
levelUnit = record.getModelInfo().getLevelUnit()
level = LevelFactory.getInstance().getLevel(levelName, levelOneValue, levelTwoValue, levelUnit)
record.getModelInfo().setLevel(level)
else :
record.setGlevel1( int(record.getDecodedLevel1) )
record.setGlevel2( int(record.getDecodedLevel2) )
LogStream.logEvent("No vertical coordinate ID for first fixed surfaced level" + str(vcrd1) + "], and second fixed surfaced [" +
str(vcrd2) + "]");
if str(g2vcrdId).isdigit() :
g2varsId = Grib2VarsTableLookup.getG2varsId(discipline, category, parameterId, pdt)
# print ("5=", discipline, category, parameterId, pdt,g2varsId, parm)
if g2varsId > 0 :
parm = Grib2VarsTableLookup.getVarGnam(discipline, category, parameterId, pdt)
modelName = record.getModelName()
#if self.derived == 'mean':
# parm = parm + "ENMW"
# if modelName.find('Mean') == -1:
# modelName = modelName +"Mean"
#if self.derived == 'sprd':
# parm = parm + "ENSA"
# if modelName.find('Spread') == -1 :
# modelName = modelName +"Spread"
#elif self.derived == 'prob':
# if modelName.find('PROB') == -1 :
# modelName = modelName +"PROB"
#elif self.derived == 'mode':
#parm = parm + "ENMO"
# if modelName.find('MODE') == -1 :
# modelName = modelName +"MODE"
#elif self.derived == '10p':
# parm = parm + "PROB"
# if modelName.find ('10P') == -1:
# modelName = modelName +"10P"
#elif self.derived == '50p':
# if modelName.find ('50P') == -1:
# modelName = modelName + "50P"
#elif self.derived == '90p':
# if modelName.find ('90P') == -1:
# modelName = modelName + "90P"
#elif self.derived == 'ctl1' or self.derived == 'ctl2':
# if modelName.find ('CTL') == -1:
# modelName = modelName + self.derived
record.setModelName(modelName)
record.getModelInfo().setModelName(modelName)
if pdt == 8 :
record.setParm(parm.replace("--", charHour))
else :
record.setParm(parm)
record.setProcessType(self.count)
if g2varsId > 0 :
# the parameter ID is in the table so append the record
records.append(record)
else :
if discipline != 255:
LogStream.logEvent("No variable ID for discipline:", discipline, "category =",category, "parameterId =", parameterId, "pdt=", pdt);
self.derived = None
self.count = self.count + 1
Reference in a new issue View git blame Copy permalink