import jarray #import Helper from util import Helper from visad import * from visad.util import DataUtility from java.io import FileInputStream from java.util import Properties from visad.python.JPythonMethods import * from java.lang import Float, Math, System, Object from java.net import URL from visad.data.units import * false = 0 true = 1 #--------------------------------------------------------------------------------- class AbstractNameAssociationTable: def __init__(self, filename): self.filename = filename #- this table's filename fis = FileInputStream("./multispectral/properties/"+filename) p = Properties() p.load(fis) self.p = p self.instrumentNameAttribute = p.getProperty("instrumentNameAttribute") self.instrumentName = p.getProperty("instrumentName") self.observation = p.getProperty("observation") self.observationChannels = p.getProperty("observationChannels") self.channelIntervals = p.getProperty("obsChannelInterval") self.channelIndex = p.getProperty("channelIndex") self.obsElement = p.getProperty("obsElement") self.obsLine = p.getProperty("obsLine") self.obsLatitude = p.getProperty("obsLatitude") self.obsLongitude = p.getProperty("obsLongitude") self.segmentBounds = p.getProperty("segmentBounds") self.obsFOVangle = p.getProperty("obsFOVangle") self.obsTimeOfDay = p.getProperty("obsTimeOfDay") self.geoLonElement = p.getProperty("geoLonElement") self.geoLonLine = p.getProperty("geoLonLine") self.geoLatElement = p.getProperty("geoLatElement") self.geoLatLine = p.getProperty("geoLatLine") #- verify this name translation table with file opened by reader # def verifyTable(self, reader): names = [] a = self.observation.split(",") for i in xrange(len(a)): names.append(a[i]) a = self.channelIndex.split(",") for i in xrange(len(a)): names.append(a[i]) a = self.obsElement.split(",") for i in xrange(len(a)): names.append(a[i]) a = self.obsLine.split(",") for i in xrange(len(a)): names.append(a[i]) names.append(self.obsLatitude) names.append(self.obsLongitude) for i in xrange(len(names)): if names[i] != None: ok = reader.verifyName(names[i]) if (ok == false): return false return true #--------------------------------------------------------------------------- #- Find a local/remote reader for this ; return an # instanceof MultiSpectralData obj which has methods for accessing # spectral/images in terms of VisAD data objects, and dataset query. # Note: MultiSpectralData is generic and can be subclassed for # issues related to specific instrument data. #- channelPropertyFile applies mainly to AIRS, but other instruments # may have the spectral files defined outside the primary data file. def load(filename, channelPropertyFile=None): from os import listdir files = listdir("./multispectral/properties") #- translation tables tablesList = [] for x in files: table = AbstractNameAssociationTable(x) tablesList.append(table) import file from file import * readers = file.__all__ reader = None readerTable = None for x in readers: try: reader = eval(x+"."+x+"(filename)") except: pass else: reader_name = x for y in tablesList: """ #- old code, files should have an attribute for the instrument type if instrName == None: if reader.getInstrumentName(y.instrumentNameAttribute) == y.instrumentName: readerTable = y else: # just get the table that matches the supplied name if instrName == y.instrumentName: readerTable = y """ if y.verifyTable(reader) == true: readerTable = y #break if reader == None: raise Exception, "No local/remote read support for: "+filename if readerTable == None: print "no translation table found, using default" readerTable = tablesList[files.index("default.txt")] instrName = readerTable.instrumentName #-- reader needs access to name association table reader.readerTable = readerTable #-- Use adapter class. This forces specific reader to implement # methods required by MultiSpectralData which will be # referenced by applications. from file.MultiSpectralFile import MultiSpectralFile reader = MultiSpectralFile(reader) #-- most instruments can use the default, MultiSpectralData, but MODIS # needs its own. if instrName == "MODIS": return ModisData(filename, reader, readerTable) #- extends MultiSpectralData elif instrName == "AIRS": return AirsData(filename, reader, readerTable, channelPropertyFile) #- extends MultiSpectralData elif instrName == "MAS": return MasData(filename, reader, readerTable) elif instrName == "MSG": return MsgData(filename, reader, readerTable) else: return MultiSpectralData(filename, reader, readerTable) #--------------------------------------------------------------------------------------------- class MultiSpectralData: def __init__(self, filename, reader, readerTable, channelPropertyFile=None, master=None): if master == None: self.filename = filename self.subset = None self.reader = reader self.geo_reader = reader #- geolocation may exist in separate file, eg. MODIS self.readerTable = readerTable self.channelPropertyFile = channelPropertyFile self.instrumentName = readerTable.instrumentName self.observationName = (readerTable.observation).split(",") self.channelIndexName = (readerTable.channelIndex).split(",") self.numObsDataSets = len(self.observationName) self.numChannelsPerDataSet = [] self.obsChannelsName = readerTable.observationChannels if self.obsChannelsName != None: self.obsChannelsName = (self.obsChannelsName).split(",") self.obsChannelIntervalsName = readerTable.channelIntervals if self.obsChannelIntervalsName != None: self.obsChannelIntervalsName = (self.obsChannelIntervalsName).split(",") self.obsElementName = (readerTable.obsElement).split(",") self.obsLineName = (readerTable.obsLine).split(",") self.geoLonElementName = readerTable.geoLonElement self.geoLatElementName = readerTable.geoLatElement self.geoLonLineName = readerTable.geoLonLine self.geoLatLineName = readerTable.geoLatLine self.obsLatitudeName = readerTable.obsLatitude self.obsLongitudeName = readerTable.obsLongitude self.segmentBounds = readerTable.segmentBounds self.obsTimeOfDayName = readerTable.obsTimeOfDay self.obsFOVangleName = readerTable.obsFOVangle self.makeRealTypes() self.makeSpectrumDomainType() self.numObsElements = self.reader.getNumObsElements(self.obsElementName[0]) self.numObsLines = self.reader.getNumObsLines(self.obsLineName[0]) self.numObsChannels = 0 for i in xrange(self.numObsDataSets): self.numChannelsPerDataSet.append(self.reader.getNumObsChannels(self.channelIndexName[i])) self.numObsChannels += self.numChannelsPerDataSet[i] if self.geoLonElementName != None and self.geoLonLineName != None: self.numGeoElements = self.reader.getNumGeoElements(self.geoLonElementName) self.numGeoLines = self.reader.getNumGeoLines(self.geoLonLineName) else: self.numGeoElements = self.numObsElements self.numGeoLines = self.numObsLines if self.numObsChannels <= 75: self.broadBand = true else: self.broadBand = false #-- for MODIS internal L1B low-res nav self.ratio = self.numObsLines/self.numGeoLines #--sub sampling factor, really a work around for the Java-NetCDF which # doesnt currently have stride/skip reading self.subFactor = 1 if self.numObsElements > 1000: #self.subFactor = 4 self.subFactor = 1 self.makeTrackSegments() #- get channels arrays self.band_numbers = None self.observationChannels = None self.observationChannels = self.getObservationChannels() self.observationChannelIntervals = None self.observationChannelIntervals = self.getObservationChannelIntervals() #- init geo range self.obsLonMinMax = None self.obsLatMinMax = None self.currentSegmentIndex = 0 self.currentSpectrum = None self.currentImage = None self.currentImageCS = None self.observationField = None self.observationDomain = None #---------------------------------------------------------------------------- #self.obsDimOrder = self.reader.getObsDimOrder() self.obsDimOrder = self.reader.getDimOrder(self.observationName[0],[self.channelIndexName[0],self.obsElementName[0],self.obsLineName[0]]) #-default(assumption) 0: channel, 1: element, 2: line; rightmost varies fastest self.xlength = self.numObsLines self.ylength = self.numObsElements self.spatialTupleArray = [self.obsLineType, self.obsElementType] self.invert = false #if (self.obsDimOrder.index(1) > self.obsDimOrder.index(2)): if (self.obsDimOrder[1] > self.obsDimOrder[2]): self.xlength = self.numObsElements self.ylength = self.numObsLines self.spatialTupleArray[0] = self.obsElementType self.spatialTupleArray[1] = self.obsLineType self.invert = true if self.geoLatElementName != None and self.geoLatLineName != None: dim_list = [self.geoLatElementName, self.geoLatLineName] else: dim_list = [self.obsElementName[0], self.obsLineName[0]] self.geoDimOrder = self.reader.getDimOrder(self.obsLatitudeName, dim_list) #-------------------------------------------------------------------------- self.obsIsRadiance = true #- false indicates brightness temperature self.native_param = "radiance" self.paramName = self.native_param self.makeRangeTypes() self.makeSpectrumRangeType(self.paramName) self.makeImageRangeType(self.paramName) self.geoOffset = 0 else: self.copy(master) #----------------------------------------------------------------------------- def copy(self, master): self.filename = master.filename self.reader = master.reader self.geo_reader = master.geo_reader self.readerTable = master.readerTable self.channelPropertyFile = master.channelPropertyFile self.instrumentName = master.instrumentName self.observationName = master.observationName self.obsChannelsName = master.obsChannelsName self.channelIndexName = master.channelIndexName self.obsElementName = master.obsElementName self.obsLineName = master.obsLineName self.geoLonElementName = master.geoLonElementName self.geoLatElementName = master.geoLatElementName self.geoLonLineName = master.geoLonLineName self.geoLatLineName = master.geoLatLineName self.obsLatitudeName = master.obsLatitudeName self.obsLongitudeName = master.obsLongitudeName self.segmentBounds = master.segmentBounds self.obsTimeOfDayName = master.obsTimeOfDayName self.obsFOVangleName = master.obsFOVangleName self.obsChannelIntervalsName = master.obsChannelIntervalsName self.observationType = master.observationType self.obsChannelsIndexType = master.obsChannelsIndexType self.obsElementType = master.obsElementType self.obsLineType = master.obsLineType self.obsLongitudeType = master.obsLongitudeType self.obsLatitudeType = master.obsLatitudeType self.obsChannelsType = master.obsChannelsType self.numObsElements = master.numObsElements self.numObsLines = master.numObsLines self.numObsChannels = master.numObsChannels self.numGeoElements = master.numGeoElements self.numGeoLines = master.numGeoLines self.ratio = master.ratio self.subFactor = master.subFactor self.trackSegments = master.trackSegments self.numberSegments = master.numberSegments self.observationChannels = master.observationChannels self.observationChannelIntervals = master.observationChannelIntervals self.obsLonMinMax = master.obsLonMinMax self.obsLatMinMax = master.obsLatMinMax self.currentSegmentIndex = master.currentSegmentIndex self.currentSpectrum = master.currentSpectrum self.currentImage = master.currentImage self.currentImageCS = master.currentImageCS self.observationField = None self.obsDimOrder = master.obsDimOrder self.xlength = master.xlength self.ylength = master.ylength self.invert = master.invert self.spatialTupleArray = master.spatialTupleArray self.geoDimOrder = master.geoDimOrder self.fileIndexes = master.fileIndexes self.obsIsRadiance = master.obsIsRadiance self.native_param = master.native_param self.paramName = master.paramName self.numObsDataSets = master.numObsDataSets self.numChannelsPerDataSet = master.numChannelsPerDataSet self.imageRangeType = master.imageRangeType self.spectrumRangeType = master.spectrumRangeType self.broadBand = master.broadBand self.observationDomain = master.observationDomain self.geoOffset = master.geoOffset def clone(self): return MultiSpectralData(self.filename,self.reader,self.readerTable,master=self) #--- VisAD RealTypes def makeRealTypes(self): name = self.readerTable.p.getProperty(self.observationName[0]) if name == None: name = self.observationName[0] self.observationType = getRealType(name) name = self.readerTable.p.getProperty(self.channelIndexName[0]) if name == None: name = self.channelIndexName[0] self.obsChannelsIndexType = getRealType(name) name = self.readerTable.p.getProperty(self.obsElementName[0]) if name == None: name = self.obsElementName[0] self.obsElementType = getRealType(name) name = self.readerTable.p.getProperty(self.obsLineName[0]) if name == None: name = self.obsLineName[0] self.obsLineType = getRealType(name) self.obsLongitudeType = RealType.Longitude self.obsLatitudeType = RealType.Latitude def makeSpectrumDomainType(self): if self.obsChannelsName != None: centimeter = ScaledUnit(0.01, CommonUnit.meter, "cm") inv_cm = centimeter.pow(-1) inv_cm = ScaledUnit(1, inv_cm, "cm^-1") self.obsChannelsType = getRealType("wavenumber", inv_cm) else: self.obsChannelsType = self.obsChannelsIndexType def makeRangeTypes(self): self.rangeTypes = [] self.rangeTypes.append(RealType.getRealType("radiance")) self.rangeTypes.append(RealType.getRealType("brightnessTemp")) self.spectrumRangeTypeDict = {"radiance" : RealType.getRealType("radiance"), "brightnessTemp" : RealType.getRealType("brightnessTemp")} self.imageRangeTypeDict = self.spectrumRangeTypeDict def makeSpectrumRangeType(self, paramName): self.spectrumRangeType = self.spectrumRangeTypeDict[paramName] def makeImageRangeType(self, paramName): self.imageRangeType = self.imageRangeTypeDict[paramName] def makeTrackSegments(self): if self.segmentBounds == None: self.trackSegments = None else: self.trackSegments = self.reader.getTrackSegments(self.segmentBounds) if (self.trackSegments != None): """ any_bad = false for i in xrange(len(self.trackSegments)): if self.trackSegments[i] > self.numObsLines: any_bad = true if any_bad == true: self.numberSegments = 1 self.trackSegments = [1, self.numObsLines] else: """ self.numberSegments = len(self.trackSegments)/2 else: #-- not defined, so set to 1 self.numberSegments = 1 #self.trackSegments = [0, self.numObsLines-1] #- one's based here self.trackSegments = [1, self.numObsLines] def getObservationChannels(self): self.fileIndexes = None if self.observationChannels == None: if self.obsChannelsName != None: self.observationChannels = [] for i in xrange(self.numObsDataSets): self.observationChannels += list(self.reader.getValues(self.obsChannelsName[i])) else: self.observationChannels = [i for i in xrange(self.numObsChannels)] return self.observationChannels def getObservationChannelIntervals(self): if self.observationChannelIntervals == None: self.observationChannelIntervals = self.reader.getValues(self.obsChannelIntervalsName) return self.observationChannelIntervals def getObservations(self, line, element, copy=false): if (self.invert == true): tmp = element element = line line = tmp line, element = self.missingObs.dataToSegment(line,element) if (self.subset != None): line *= self.step element *= self.step if (self.invert == true): line += self.subset.start_y element += self.subset.start_x else: line += self.subset.start_x element += self.subset.start_y else: line *= self.step element *= self.step ll = int(line + (self.trackSegments[self.currentSegmentIndex*2]-1)) #convert one to zero based #-- file access vals = jarray.zeros(self.numObsChannels, 'f') cnt = 0 for i in xrange(self.numObsDataSets): start = [0,element+self.geoOffset,ll+self.geoOffset] count = [self.numChannelsPerDataSet[i],1,1] nstart = [None, None, None] ncount = [None, None, None] for j in xrange(len(start)): nstart[self.obsDimOrder[j]] = start[j] ncount[self.obsDimOrder[j]] = count[j] vv = self.reader.getValues(self.observationName[i], nstart, ncount) for k in xrange(len(vv)): vals[cnt] = float(vv[k]) cnt += 1 vals = self.getSpectrumRange(vals, self.paramName) #-- make the VisAD data object if self.observationField == None or copy == true: self.observationDomain = Gridded1DSet(self.obsChannelsType, [self.observationChannels], self.numObsChannels) if self.broadBand == false: self.observationField = FlatField(FunctionType(self.obsChannelsType, self.spectrumRangeType), self.observationDomain) self.observationField.setSamples([vals]) else: if (isinstance(self.spectrumRangeType, RealTupleType)): fields = [] for i in xrange(self.spectrumRangeType.getDimension()): range_type = RealTupleType([self.obsChannelsType]+[self.spectrumRangeType.getComponent(i)]) ff = FlatField(FunctionType(RealType.Generic, range_type), Integer1DSet(self.numObsChannels)) ff.setSamples([self.observationChannels, vals[i]]) fields.append(ff) self.observationField = Tuple(fields) else: range_type = RealTupleType(self.obsChannelsType, self.spectrumRangeType) self.observationField = FlatField(FunctionType(RealType.Generic, range_type ), Integer1DSet(self.numObsChannels)) self.observationField.setSamples([self.observationChannels, vals]) else: if self.broadBand == false: self.observationField.setSamples([vals]) else: if (isinstance(self.spectrumRangeType, RealTupleType)): for i in xrange(self.spectrumRangeType.getDimension()): self.observationField.getComponent(i).setSamples([self.observationChannels, vals[i]]) else: self.observationField.setSamples([self.observationChannels, vals]) return self.observationField def getImageSegment(self, channelIndex, segmentIndex=0, subset=None): if channelIndex < 0: return None dataset_index = self.channelIndexToDataSetIndex(channelIndex) channelIndex = self.channelIndexToFileIndex(dataset_index, channelIndex) self.currentSegmentIndex = segmentIndex segStartLine = int(self.trackSegments[self.currentSegmentIndex*2]-1) #convert one to zero based segEndLine = int(self.trackSegments[self.currentSegmentIndex*2+1]-1) if subset != None: #-- subset logic self.subset = subset if subset.end_x == None or subset.end_y == None: self.startElem = 0 self.nElements = int((self.numObsElements)/subset.step) self.startLine = segStartLine self.nLines = int((int(segEndLine-self.startLine+1))/subset.step) else: t_line = subset.start_x t_elem = subset.start_y n_lines = int((subset.end_x - subset.start_x)/subset.step) + 1 n_elems = int((subset.end_y - subset.start_y)/subset.step) + 1 if (self.invert == true): t_elem = subset.start_x t_line = subset.start_y n_elems = int((subset.end_x - subset.start_x)/subset.step) + 1 n_lines = int((subset.end_y - subset.start_y)/subset.step) + 1 self.startElem = t_elem self.startLine = segStartLine + t_line self.nElements = n_elems self.nLines = n_lines self.step = self.subset.step if float(self.ratio)/float(self.step) <= 1: self.geo_step = int(1.0/(float(self.ratio)/float(self.step))) nGeoLines = self.nLines nGeoElements = self.nElements else: self.geo_step = 1 nGeoLines = int(float(self.nLines)/float(self.ratio)/float(self.step)) nGeoElements = int(float(self.nElements)/float(self.ratio)/float(self.step)) else: self.subset = None self.startElem = 0 self.nElements = self.numObsElements self.startLine = segStartLine self.nLines = int(segEndLine-self.startLine+1) self.step = 1 if float(self.ratio)/float(self.step) <= 1: self.geo_step = int(1.0/(float(self.ratio)/float(self.step))) nGeoLines = self.nLines nGeoElements = self.nElements else: self.geo_step = 1 nGeoLines = int(float(self.nLines)/float(self.ratio)/float(self.step)) nGeoElements = int(float(self.nElements)/float(self.ratio)/float(self.step)) start = [channelIndex,self.startElem+self.geoOffset,self.startLine+self.geoOffset] count = [1,self.nElements,self.nLines] stride = [1, self.step, self.step] nstart = [None, None, None] ncount = [None, None, None] nstride = [None, None, None] for i in xrange(len(start)): nstart[self.obsDimOrder[i]] = start[i] ncount[self.obsDimOrder[i]] = count[i] nstride[self.obsDimOrder[i]] = stride[i] if self.step > 1: #- this hack will go away when Java-Netcdf implements stride reading self.z = self.reader.getValues(self.observationName[dataset_index], nstart, ncount, nstride) else: self.z = self.reader.getValues(self.observationName[dataset_index], nstart, ncount, None) start = [self.startElem/self.ratio, self.startLine/self.ratio] count = [nGeoElements, nGeoLines] stride = [self.geo_step, self.geo_step] nstart = [] ncount = [] nstride = [] for i in xrange(len(start)): nstart.append(start[self.geoDimOrder[i]]) ncount.append(count[self.geoDimOrder[i]]) nstride.append(stride[self.geoDimOrder[i]]) if self.step > 1: self.lons = self.geo_reader.getValues(self.obsLongitudeName, nstart, ncount, nstride) self.lats = self.geo_reader.getValues(self.obsLatitudeName, nstart, ncount, nstride) else: self.lons = self.geo_reader.getValues(self.obsLongitudeName, nstart, ncount, None) self.lats = self.geo_reader.getValues(self.obsLatitudeName, nstart, ncount, None) if self.invert == true: self.ylength = self.nLines self.xlength = self.nElements self.gylength = nGeoLines self.gxlength = nGeoElements else: self.xlength = self.nLines self.ylength = self.nElements self.gxlength = nGeoLines self.gylength = nGeoElements self.missingObs = self.initMissing(self.lons, self.lats, self.z) self.xlength, self.ylength = self.missingObs.getLengths() self.z = self.missingObs.removeMissing(self.z) self.lons = self.missingObs.removeMissing(self.lons) self.lats = self.missingObs.removeMissing(self.lats) if self.ratio == 1: self.gxlength = self.xlength self.gylength = self.ylength try: self.gs = Gridded2DSet(RealTupleType(self.obsLatitudeType,self.obsLongitudeType), [self.lats, self.lons], self.gxlength,self.gylength,None,None,None,0,0) self.cs = GridCoordinateSystem(self.gs) except: self.cs = None print "*** Error Making Image CoordinateSystem ***" self.imageDomType = RealTupleType(self.spatialTupleArray[0], self.spatialTupleArray[1], self.cs,None) self.imageDomain = makeDomain(self.imageDomType, 0,((self.xlength)/(self.xlength/self.gxlength))-1,self.xlength, 0,((self.ylength)/(self.ylength/self.gylength))-1,self.ylength) self.subsetDomain = makeDomain(self.imageDomType, 0, (self.xlength-1)*self.step, self.xlength, 0, (self.ylength-1)*self.step, self.ylength) self.imageFunctionType = FunctionType(self.imageDomType, self.imageRangeType) self.imageField = FlatField(self.imageFunctionType, self.imageDomain) self.z = self.getImageRange(self.z,self.paramName,dataset_index,channelIndex) indexes = jarray.zeros(2, 'i') minmax = Helper.minmax(self.z, indexes) minloc = self.cs.toReference(self.imageDomain.indexToValue([indexes[0]])) maxloc = self.cs.toReference(self.imageDomain.indexToValue([indexes[1]])) self.obsMin = [minloc[0][0], minloc[1][0], minmax[0]] self.obsMax = [maxloc[0][0], maxloc[1][0], minmax[1]] self.imageField.setSamples([self.z], false) return self.imageField def getSameImageSegment(self, channelIndex): if channelIndex >= 0: dataset_index = self.channelIndexToDataSetIndex(channelIndex) channelIndex = self.channelIndexToFileIndex(dataset_index, channelIndex) start = [channelIndex,self.startElem+self.geoOffset,self.startLine+self.geoOffset] count = [1,self.nElements,self.nLines] stride = [1, self.step, self.step] """ nstart = [] ncount = [] nstride = [] for i in xrange(len(start)): nstart.append(start[self.obsDimOrder[i]]) ncount.append(count[self.obsDimOrder[i]]) nstride.append(stride[self.obsDimOrder[i]]) """ nstart = [None, None, None] ncount = [None, None, None] nstride = [None, None, None] for j in xrange(len(start)): nstart[self.obsDimOrder[j]] = start[j] ncount[self.obsDimOrder[j]] = count[j] nstride[self.obsDimOrder[j]] = stride[j] if self.step > 1: self.z = self.reader.getValues(self.observationName[dataset_index],nstart,ncount,nstride) else: self.z = self.reader.getValues(self.observationName[dataset_index],nstart,ncount) self.z = self.missingObs.removeMissing(self.z) self.z = self.getImageRange(self.z, self.paramName,dataset_index, channelIndex) indexes = jarray.zeros(2, 'i') minmax = Helper.minmax(self.z, indexes) minloc = self.cs.toReference(self.imageDomain.indexToValue([indexes[0]])) maxloc = self.cs.toReference(self.imageDomain.indexToValue([indexes[1]])) self.obsMin = [minloc[0][0], minloc[1][0], minmax[0]] self.obsMax = [maxloc[0][0], maxloc[1][0], minmax[1]] if (self.subFactor > 1): tmpField = FlatField(self.imageFunctionType,self.imageDomain) tmpField.setSamples([self.z]) tmpDomain = makeDomain(self.imageDomType,0,self.xlength-1,self.xlength/self.subFactor, 0,self.ylength-1,self.ylength/self.subFactor) tmpField = tmpField.resample(tmpDomain) self.imageField.setSamples(Helper.getValues(tmpField,false), false) else: self.imageField.setSamples([self.z], false) def getSameImageSegmentCopy(self, channelIndex): if channelIndex >= 0: dataset_index = self.channelIndexToDataSetIndex(channelIndex) channelIndex = self.channelIndexToFileIndex(dataset_index, channelIndex) start = [channelIndex,self.startElem+self.geoOffset,self.startLine+self.geoOffset] count = [1,self.nElements,self.nLines] stride = [1,self.step,self.step] """ nstart = [] ncount = [] nstride = [] for i in xrange(len(start)): nstart.append(start[self.obsDimOrder[i]]) ncount.append(count[self.obsDimOrder[i]]) nstride.append(stride[self.obsDimOrder[i]]) """ nstart = [None, None, None] ncount = [None, None, None] nstride = [None, None, None] for j in xrange(len(start)): nstart[self.obsDimOrder[j]] = start[j] ncount[self.obsDimOrder[j]] = count[j] nstride[self.obsDimOrder[j]] = stride[j] if self.step > 1: self.z = self.reader.getValues(self.observationName[dataset_index],nstart,ncount,nstride) else: self.z = self.reader.getValues(self.observationName[dataset_index],nstart,ncount) self.z = self.missingObs.removeMissing(self.z) new_imageField = FlatField(self.imageFunctionType, self.imageDomain) self.z = self.getImageRange(self.z, self.paramName, dataset_index, channelIndex) indexes = jarray.zeros(2, 'i') minmax = Helper.minmax(self.z, indexes) minloc = self.cs.toReference(self.imageDomain.indexToValue([indexes[0]])) maxloc = self.cs.toReference(self.imageDomain.indexToValue([indexes[1]])) self.obsMin = [minloc[0][0], minloc[1][0], minmax[0]] self.obsMax = [maxloc[0][0], maxloc[1][0], minmax[1]] new_imageField.setSamples([self.z], false) if self.subFactor > 1: tmpDomain = makeDomain(self.imageDomType, 0,self.xlength-1,self.xlength/self.subFactor, 0,self.ylength-1,self.ylength/self.subFactor) new_imageField = new_imageField.resample(tmpDomain) return new_imageField def changeParameter(self, paramName, channelIndex, line, element): self.paramName = paramName self.spectrumRangeType = self.spectrumRangeTypeDict[paramName] self.imageRangeType = self.imageRangeTypeDict[paramName] self.imageFunctionType = FunctionType(self.imageDomType, self.imageRangeType) self.imageField = self.getSameImageSegmentCopy(channelIndex) self.observationField = self.getObservations(line, element, true) return self.imageField, self.observationField def getGeoDomain(self): if self.obsLonMinMax == None: start = [0,0] count = [self.numGeoElements, self.numGeoLines] nstart = [] ncount = [] for i in xrange(len(start)): nstart.append(start[self.geoDimOrder[i]]) ncount.append(count[self.geoDimOrder[i]]) #if self.numGeoElements > 512: if self.numGeoElements > 5000: ncount[0] = ncount[0]/10 ncount[1] = ncount[1]/10 nstride = [10,10] else: nstride = None nstride = None lons = self.geo_reader.getValues(self.obsLongitudeName, nstart, ncount, nstride) lats = self.geo_reader.getValues(self.obsLatitudeName, nstart, ncount, nstride) self.obsLonMinMax = Helper.minmax(lons, -180, 180) self.obsLatMinMax = Helper.minmax(lats, -90, 90) return self.obsLatMinMax, self.obsLonMinMax def getTimeObs(self, line, element): if (self.invert == true): tmp = element element = line line = tmp line, element = self.missingObs.dataToSegment(line, element) if (self.subset != None): if (self.invert == true): line += self.subset.start_y element += self.subset.start_x else: line += self.subset.start_x element += self.subset.start_y ll = int(line + (self.trackSegments[self.currentSegmentIndex*2]-1)) return self.reader.getValue(self.obsTimeOfDayName, [element, ll]) def getFOVangleObs(self, line, element): if (self.invert == true): tmp = element element = line line = tmp line, element = self.missingObs.dataToSegment(line,element) if (self.subset != None): if (self.invert == true): line += self.subset.start_y element += self.subset.start_x else: line += self.subset.start_x element += self.subset.start_y ll = int(line + (self.trackSegments[self.currentSegmentIndex*2]-1)) return self.reader.getValue(self.obsFOVangleName, [element, ll]) def getImageValue(self, lat, lon): val = self.imageField.evaluate(RealTuple(self.imageField.getType().getDomain().getCoordinateSystem().getReference(), [lat, lon]), Data.NEAREST_NEIGHBOR, Data.NO_ERRORS) return val def getRealTypes(self): #return [self.obsElementType, self.obsLineType, None, self.obsChannelsType, self.observationType] return [self.obsElementType, self.obsLineType, None, self.obsChannelsType, self.imageRangeType] def getGeoRealTypes(self): return [self.obsLatitudeType, self.obsLongitudeType] def getObservationChannelRange(self): if self.broadBand == false: #return (self.observationChannels[0], self.observationChannels[self.numObsChannels-1]) return (min(self.observationChannels), max(self.observationChannels)) else: #- return band numbers, ie MODIS or GOES return (min(self.observationChannels), max(self.observationChannels)) #return (self.band_numbers[0], self.band_numbers[self.numObsChannels-1]) def getObsValidRange(self): #- valid is probably a misnomer here; mostly for setRange if self.paramName == "radiance": return [0, 150] elif self.paramName == "brightnessTemp": return [180, 320] def getNumSegments(self): return self.numberSegments def initMissing(self, longitudes, latitudes, values): class Missing: #- this class must implement these methods def __init__(self, longitudes, latitudes, values, xlength, ylength, invert): self.longitudes = longitudes self.latitudes = latitudes self.xlength = xlength self.ylength = ylength self.invert = invert def dataToSegment(self, line, element): return line, element def removeMissing(self, values): return values def getLengths(self): return self.xlength, self.ylength m = Missing(longitudes, latitudes, values, self.xlength, self.ylength, self.invert) return m def channelIndexToFileIndex(self, dataset_index, channel_index): return channel_index def channelIndexToDataSetIndex(self, index): return 0 def getSpectrumRange(self, values, paramName): if self.native_param == paramName: return values elif paramName == "brightnessTemp": return self.radianceToBrightnessTempSpectrum(values) elif paramName == "radiance": return self.brightnessTempToRadiance(values,channelIndex) def getImageRange(self, values, paramName, datasetIndex, channelIndex): if self.native_param == paramName: self.imageRangeUnit = ScaledUnit(1, Parser.parse("mW/steradian/m^2/cm^-1"), "mW/ster/m^2/cm^-1") return values elif paramName == "brightnessTemp": self.imageRangeUnit = Parser.parse("K") return self.radianceToBrightnessTemp(values,channelIndex) elif paramName == "radiance": self.imageRangeUnit = ScaledUnit(1, Parser.parse("mW/steradian/m^2/cm^-1"), "mW/ster/m^2/cm^-1") return self.brightnessTempToRadiance(values,channelIndex) def radianceToBrightnessTemp(self, values, channelIndex): # Converts radiances [mW/ster/m2/cm^-1] to BT [K] # Input: nu array of wavenmbers [cm^-1] # B radiances [mW/ster/m2/cm^-1] # Output: bt brightness temperature in [K] # Paolo Antonelli # Wed Feb 25 16:43:05 CST 1998 c1=1.191066E-5; #- mW/m2/ster/cm^-4 c2=1.438833; #- K*cm wavenumber = self.observationChannels[channelIndex] nu = wavenumber nvals = len(values) new_values = jarray.zeros(nvals, 'f') for i in xrange(nvals): B = values[i] bt=c2*nu/(Math.log(((c1*nu*nu*nu)/B)+1.0)); new_values[i] = bt return new_values def radianceToBrightnessTempSpectrum(self, values): # Converts radiances [mW/ster/m2/cm^-1] to BT [K] # Input: nu array of wavenmbers [cm^-1] # B radiances [mW/ster/m2/cm^-1] # Output: bt brightness temperature in [K] # Paolo Antonelli # Wed Feb 25 16:43:05 CST 1998 c1=1.191066E-5; #- mW/m2/ster/cm^-4 c2=1.438833; #- K*cm nvals = len(values) new_values = jarray.zeros(nvals, 'f') for i in xrange(nvals): nu = self.observationChannels[i] B = values[i] bt=c2*nu/(Math.log(((c1*nu*nu*nu)/B)+1.0)); new_values[i] = bt return new_values def brightnessTempToRadiance(self, values, wavenumber): # compute radiance given wavenumbers and brightness temperature. # # Inputs: # freq wavenumbers (Nx1) in 1/cm # bt brightnes temperature (Nx1) in Kelvin # Outputs: # radiance Planck radiances (Nx1) in mW/(m2 sr. 1/cm) # # see also RADTOT.M, RAD2BT.M, TTORAD.M # # DCT 11/11-99 # # fundamental constants: # (Cohen, E.R. and B.N. Taylor, The 1986 CODATA recommended values # of the fundamental physical constants, Journal of Research of # the National Bureau of Standard, 92(2), March-April 1987.) h = 6.6260755E-34; # Planck constant in Js c = 2.99792458E8; # photon speed in m/s k = 1.380658E-23; # Boltzmann constant in J/K c1 = 2*h*c*c*1e8; c2 = h*c/k*1e2; freq = wavenumber nvals = len(values) for i in xrange(nvals): bt = values[i] radiance = 1e3 * c1*(freq*freq*freq)/(Math.exp((c2*freq)/bt)-1); new_values[i] = radiance return new_values #------------------------------------------------------------------------------------------------ #- eventually break these out into separate modules class ModisData(MultiSpectralData): def __init__(self, filename, reader, readerTable,master=None): MultiSpectralData.__init__(self, filename, reader, readerTable, master) #-look for MODIS geo file in same directory as L1B if master == None: if filename.find("1000m") != -1: geo_filename = filename.replace("1000m", "geo") elif filename.find("MOD02") != -1: geo_filename = filename.replace("MOD02", "MOD03") elif filename.find("MYD02") != -1: geo_filename = filename.replace("MYD02", "MYD03") reader_name = reader.reader.__class__ reader_name = reader_name.__name__ if ((filename.find("MYD") != -1) or (filename.find("a1.") != -1)): self.modis_platform = "Aqua" else: self.modis_platform = "Terra" import file from file import * #try: """ self.geo_reader = eval(reader_name+"."+reader_name+"(geo_filename)") dim_list = ["mframes", "nscans*10"] self.geoDimOrder = self.geo_reader.getDimOrder(self.obsLatitudeName, dim_list) if self.geoDimOrder == []: dim_list = ["mframes:MODIS_Swath_Type_GEO", "nscans*10:MODIS_Swath_Type_GEO"] self.geoDimOrder = self.geo_reader.getDimOrder(self.obsLatitudeName, dim_list) """ #except: #self.geo_reader = self.reader self.getScaleOffset() self.geoOffset = 2 else: self.geo_reader = master.geo_reader self.geoDimOrder = master.geoDimOrder self.modis_platform = master.modis_platform self.radScale = master.radScale self.radOffset = master.radOffset self.reflScale = master.reflScale self.reflOffset = master.reflOffset self.geoOffset = master.geoOffset def getObservationChannels(self): channels = [] self.band_numbers = [] ff = open('./ancillary/modis/default_cntr_wn_bw.txt', 'r') lines = ff.readlines() nlines = len(lines) self.bandDict = {} cnt = 0 for i in xrange(nlines): if lines[i].startswith('!'): pass else: fields = lines[i].split() band_number = fields[0] center_wn = fields[1] low_wn = fields[2] hi_wn = fields[3] center_wn = float(center_wn) if band_number == "13lo": center_wn -= center_wn/1000 if band_number == "13hi": center_wn += center_wn/1000 if band_number == "14lo": center_wn -= center_wn/1000 if band_number == "14hi": center_wn += center_wn/1000 self.band_numbers.append(band_number) channels.append(center_wn) self.bandDict.update({band_number : [cnt,center_wn]}) self.bandDict.update({center_wn : band_number}) cnt += 1 self.fileIndexes = QuickSort.sort(channels) new_channels = jarray.zeros(len(channels), 'f') for i in xrange(len(new_channels)): new_channels[i] = channels[self.fileIndexes[i]] return new_channels #- override def makeTrackSegments(self): if self.numObsLines < 2600: self.numberSegments = 1 self.trackSegments = [1, self.numObsLines] return """ self.trackSegments = [] segLen = 2000 start = 401 #-track segment array is one's based self.numberSegments = int((self.numObsLines-100)/segLen) for i in xrange(self.numberSegments): start = start + segLen*i self.trackSegments.append(start) self.trackSegments.append(start + (segLen-1)) """ self.numberSegments = 1 self.trackSegments = [401, self.numObsLines-100] def getScaleOffset(self): self.radScale = [] self.radOffset = [] self.reflScale = [] self.reflOffset = [] for i in xrange(self.numObsDataSets): radScale, radOffset = self.reader.getObsScaleOffset(self.observationName[i], "radiance_scales", "radiance_offsets") reflScale, reflOffset = self.reader.getObsScaleOffset(self.observationName[i], "reflectance_scales", "reflectance_offsets") if radScale == None or radOffset == None: self.radScale.append(None) self.radOffset.append(None) else: scales = jarray.zeros(len(radScale), 'f') offsets = jarray.zeros(len(radOffset), 'f') for i in xrange(len(scales)): scales[i] = radScale[i] offsets[i] = radOffset[i] self.radScale.append(scales) self.radOffset.append(offsets) if reflScale == None or reflOffset == None: self.reflScale.append(None) self.reflOffset.append(None) else: scales = jarray.zeros(len(reflScale), 'f') offsets = jarray.zeros(len(reflOffset), 'f') for i in xrange(len(scales)): scales[i] = reflScale[i] offsets[i] = reflOffset[i] self.reflScale.append(scales) self.reflOffset.append(offsets) def getSpectrumRange(self, values, paramName): if paramName == "radiance": rad_scales = jarray.zeros(self.numObsChannels, 'f') rad_offsets = jarray.zeros(self.numObsChannels, 'f') cnt = 0 for i in xrange(self.numObsDataSets): System.arraycopy(self.radScale[i], 0, rad_scales, cnt, self.numChannelsPerDataSet[i]) System.arraycopy(self.radOffset[i], 0, rad_offsets, cnt, self.numChannelsPerDataSet[i]) cnt += self.numChannelsPerDataSet[i] vals1_20 = jarray.zeros(len(values), 'f') vals21_36 = jarray.zeros(len(values), 'f') System.arraycopy(values, 0, vals1_20, 0, len(values)) System.arraycopy(values, 0, vals21_36, 0, len(values)) vals1_20 = Helper.unpack(vals1_20, rad_scales, rad_offsets) vals21_36 = Helper.unpack(vals21_36, rad_scales, rad_offsets) for i in xrange(len(values)): if i <= 21: vals1_20[i] = vals1_20[i] vals21_36[i] = Float.NaN else: vals1_20[i] = Float.NaN vals21_36[i] = vals21_36[i] new_vals1_20 = jarray.zeros(len(values), 'f') new_vals21_36 = jarray.zeros(len(values), 'f') for i in xrange(len(values)): new_vals1_20[i] = vals1_20[self.fileIndexes[i]] new_vals21_36[i] = vals21_36[self.fileIndexes[i]] return [new_vals1_20,new_vals21_36] #- Brightness Temperature for emissive bands, Reflectance for reflective elif paramName == "brightnessTemp": rad_scales = jarray.zeros(self.numObsChannels, 'f') rad_offsets = jarray.zeros(self.numObsChannels, 'f') ref_scales = jarray.zeros(self.numObsChannels, 'f') ref_offsets = jarray.zeros(self.numObsChannels, 'f') cnt1 = 0 cnt2 = 0 for i in xrange(self.numObsDataSets): if self.reflScale[i] != None: System.arraycopy(self.reflScale[i], 0, ref_scales, cnt1, self.numChannelsPerDataSet[i]) System.arraycopy(self.reflOffset[i], 0, ref_offsets, cnt1, self.numChannelsPerDataSet[i]) cnt1 += self.numChannelsPerDataSet[i] System.arraycopy(self.radScale[i], 0, rad_scales, cnt2, self.numChannelsPerDataSet[i]) System.arraycopy(self.radOffset[i], 0, rad_offsets, cnt2, self.numChannelsPerDataSet[i]) cnt2 += self.numChannelsPerDataSet[i] vals1_20 = jarray.zeros(len(values), 'f') vals21_36 = jarray.zeros(len(values), 'f') System.arraycopy(values, 0, vals1_20, 0, len(values)) System.arraycopy(values, 0, vals21_36, 0, len(values)) vals1_20 = Helper.unpack(vals1_20, ref_scales, ref_offsets) vals21_36 = Helper.unpack(vals21_36, rad_scales, rad_offsets) for i in xrange(len(values)): if i <= 21: vals1_20[i] = vals1_20[i] vals21_36[i] = Float.NaN else: vals1_20[i] = Float.NaN vals21_36[i] = vals21_36[i] for i in xrange(16): vals21_36[22+i] = (self.radianceToBrightnessTemp([vals21_36[22+i]], 22+i))[0] new_vals1_20 = jarray.zeros(len(values), 'f') new_vals21_36 = jarray.zeros(len(values), 'f') for i in xrange(len(values)): new_vals1_20[i] = vals1_20[self.fileIndexes[i]] new_vals21_36[i] = vals21_36[self.fileIndexes[i]] return [new_vals1_20,new_vals21_36] def getImageRange(self, values, paramName, dataset_index, index): if paramName == "radiance": unit = Parser.parse("W/m^2/micrometer/steradian") self.imageRangeUnit = ScaledUnit(1, unit, "W/m^2/\xb5m/ster") return Helper.unpack(values, self.radScale[dataset_index][index], self.radOffset[dataset_index][index]) elif paramName == "brightnessTemp": self.imageRangeUnit = Parser.parse("K") if dataset_index < 3: return Helper.unpack(values, self.reflScale[dataset_index][index], self.reflOffset[dataset_index][index]) else: vals = Helper.unpack(values, self.radScale[dataset_index][index], self.radOffset[dataset_index][index]) return self.radianceToBrightnessTemp(vals, (index + 22)) def radianceToBrightnessTempSpec(self, values): new_values = jarray.zeros(len(values), 'd') for i in xrange(len(values)): new_values[i] = (self.radianceToBrightnessTemp([values[i]], self.observationChannels[i]))[0] return new_values def radianceToBrightnessTemp(self, values, channelIndex): return Helper.modis_radiance_to_brightnessTemp(self.modis_platform, int(self.band_numbers[channelIndex]), values) #-- should use units and look in configuration file? def getObsValidRange(self): if self.paramName == "radiance": return [[0,500],[0, 10]] elif self.paramName == "brightnessTemp": return [[0,1],[100, 320]] def channelIndexToDataSetIndex(self, channelIndex): idx = self.fileIndexes[channelIndex] #- get the unsorted index if idx <= 1: return 0 elif idx <= 6: return 1 elif idx <= 21: return 2 else: return 3 def clone(self): return ModisData(self.filename,self.reader,self.readerTable,master=self) def channelIndexToFileIndex(self, datasetIndex, channelIndex): idx = self.fileIndexes[channelIndex] if datasetIndex == 0: return idx else: a = 0 for i in xrange(datasetIndex): a += self.numChannelsPerDataSet[i] return idx - a def makeRangeTypes(self): self.spectrumRangeTypeDict = { "radiance": RealTupleType(RealType.getRealType("Radiance_1-19_26"), RealType.getRealType("Radiance_20-36")), "brightnessTemp": RealTupleType(RealType.getRealType("Reflectance"), RealType.getRealType("BrightnessTemp"))} self.imageRangeTypeDict = { "radiance": RealType.getRealType("Radiance"), "brightnessTemp": RealType.getRealType("Refl/BT")} self.rangeTypes = [] self.rangeTypes.append(RealType.getRealType("Radiance_1-20")) self.rangeTypes.append(RealType.getRealType("Radiance_21-36")) self.rangeTypes.append(RealType.getRealType("Reflectance")) self.rangeTypes.append(RealType.getRealType("BrightnessTemp")) def makeSpectrumDomainType(self): micrometer = ScaledUnit(0.000001, CommonUnit.meter, "\xb5m") self.obsChannelsType = RealType.getRealType("wavelength", micrometer) def makeSpectrumRangeType(self, paramName): if paramName == "radiance": self.spectrumRangeType = RealTupleType(self.rangeTypes[0], self.rangeTypes[1]) elif paramName == "brightnessTemp": self.spectrumRangeType = RealTupleType(self.rangeTypes[2], self.rangeTypes[3]) def makeImageRangeType(self, paramName): if paramName == "radiance": self.imageRangeType = RealType.getRealType("radiance") elif paramName == "brightnessTemp": self.imageRangeType = RealType.getRealType("brightnessTemp") def bandNameToValue(self, text): return self.observationChannels[list(self.fileIndexes).index((self.bandDict[text])[0])] def bandValueToName(self, value): idx = self.observationDomain.valueToIndex([[value]])[0] return self.band_numbers[self.fileIndexes[idx]] #------------------------------------------------------------------------------------------- class MsgData(MultiSpectralData): def __init__(self, filename, reader, readerTable,master=None): MultiSpectralData.__init__(self, filename, reader, readerTable, master) def getObservationChannels(self): self.fileIndexes = None self.observationChannels = self.reader.getValues(self.obsChannelsName[0]) self.bandDict = {} self.band_numbers = ['1','2','3','4','5','6','7','8','9','10','11'] self.fileIndexes = QuickSort.sort(self.observationChannels) cnt = 0 for i in xrange(len(self.observationChannels)): band_number = self.band_numbers[i] center_wn = self.observationChannels[i] self.bandDict.update({band_number : [cnt,center_wn]}) self.bandDict.update({center_wn : band_number}) cnt += 1 return self.observationChannels def channelIndexToDataSetIndex(self, channelIndex): if channelIndex <= 2: return 0 else: return 1 def channelIndexToFileIndex(self, datasetIndex, channelIndex): if datasetIndex == 0: return channelIndex else: a = 0 for i in xrange(datasetIndex): a += self.numChannelsPerDataSet[i] return channelIndex - a def makeRangeTypes(self): self.spectrumRangeTypeDict = { "radiance": RealTupleType(RealType.getRealType("Radiance_1-3"), RealType.getRealType("Radiance_4-11")), "brightnessTemp": RealTupleType(RealType.getRealType("Reflectance"), RealType.getRealType("BrightnessTemp"))} self.imageRangeTypeDict = { "radiance": RealType.getRealType("Radiance"), "brightnessTemp": RealType.getRealType("Refl/BT")} self.rangeTypes = [] self.rangeTypes.append(RealType.getRealType("Radiance_1-3")) self.rangeTypes.append(RealType.getRealType("Radiance_4-11")) self.rangeTypes.append(RealType.getRealType("Reflectance")) self.rangeTypes.append(RealType.getRealType("BrightnessTemp")) def makeSpectrumDomainType(self): micrometer = ScaledUnit(0.000001, CommonUnit.meter, "\xb5m") self.obsChannelsType = RealType.getRealType("wavelength", micrometer) def makeSpectrumRangeType(self, paramName): if paramName == "radiance": self.spectrumRangeType = RealTupleType(self.rangeTypes[0], self.rangeTypes[1]) elif paramName == "brightnessTemp": self.spectrumRangeType = RealTupleType(self.rangeTypes[2], self.rangeTypes[3]) def makeImageRangeType(self, paramName): if paramName == "radiance": self.imageRangeType = RealType.getRealType("radiance") elif paramName == "brightnessTemp": self.imageRangeType = RealType.getRealType("brightnessTemp") def getSpectrumRange(self,values, paramName): vals1_3 = jarray.zeros(len(values), 'f') vals4_11 = jarray.zeros(len(values), 'f') for i in xrange(len(values)): vals1_3[i] = Float.NaN vals4_11[i] = Float.NaN if paramName == "radiance": for i in xrange(len(values)): if i <= 2: vals1_3[i] = values[i] else: vals4_11[i] = values[i] elif paramName == "brightnessTemp": tmp = Helper.msg_radiance_to_reflectance([0,1,2], values[0:3]) for i in xrange(len(tmp)): vals1_3[i] = tmp[i] tmp = Helper.msg_radiance_to_brightnessTemperature([0,1,2,3,4,5,6,7], values[4:12]) for i in xrange(len(tmp)): vals4_11[i+3] = tmp[i] return [vals1_3,vals4_11] def getImageRange(self, values, paramName, dataset_index, index): if paramName == "radiance": unit = Parser.parse("mW/m^2/micrometer/steradian") self.imageRangeUnit = ScaledUnit(1, unit, "mW/m^2/\xb5m/ster") return values elif paramName == "brightnessTemp": self.imageRangeUnit = Parser.parse("K") if dataset_index < 1: return self.radianceToReflectance(values, index) else: return self.radianceToBrightnessTemp(values, index) def radianceToReflectance(self, values, index): return Helper.msg_radiance_to_reflectance(index, values) def radianceToBrightnessTemp(self, values, channelIndex): return Helper.msg_radiance_to_brightnessTemperature(channelIndex, values) def bandNameToValue(self, text): return (self.bandDict[text])[1] def bandValueToName(self, value): idx = self.observationDomain.valueToIndex([[value]])[0] return self.band_numbers[idx] def getObsValidRange(self): if self.paramName == "radiance": return [[0,10],[0, 100]] elif self.paramName == "brightnessTemp": return [[0,1],[100, 320]] def clone(self): return MsgData(self.filename,self.reader,self.readerTable,master=self) #------------------------------------------------------------------------------------------------ #-- break this out into separate module class AirsData(MultiSpectralData): def __init__(self,filename,reader,readerTable,channelPropertyFile,master=None): MultiSpectralData.__init__(self,filename,reader,readerTable,channelPropertyFile,master) if master != None: self.rad_quality = master.rad_quality #- override def getObservationChannels(self): #if self.channelPropertyFile == None: # return MultiSpectralData.getObservationChannels(self) #ff = open('./ancillary/airs/L2.chan_prop.2003.01.10.v6.6.8.anc', 'r') ff = open('./ancillary/airs/L2.chan_prop.2003.11.19.v6.6.9.anc', 'r') lines = ff.readlines() nlines = len(lines) channels = [] good_indexes = [] self.rad_quality = [] cnt = 0 for i in xrange(nlines): if lines[i].startswith('!'): pass else: fields = lines[i].split() srf_centroid_freq = fields[1] rq = fields[11] l2_ignore = fields[12] channels.append(float(srf_centroid_freq)) good_indexes.append(cnt) self.rad_quality.append(int(l2_ignore)) cnt += 1 self.fileIndexes = QuickSort.sort(channels) new_channels = jarray.zeros(len(channels), 'f') for i in xrange(len(new_channels)): new_channels[i] = channels[self.fileIndexes[i]] return new_channels def makeSpectrumDomainType(self): centimeter = ScaledUnit(0.01, CommonUnit.meter, "cm") inv_cm = centimeter.pow(-1) inv_cm = ScaledUnit(1, inv_cm, "cm^-1") self.obsChannelsType = getRealType("wavenumber", inv_cm) def getSpectrumRange(self, values, paramName): #first, apply QC flags from channel property file nn = len(values) new_values = jarray.zeros(len(values), 'f') for i in xrange(nn): val = values[self.fileIndexes[i]] rq = self.rad_quality[self.fileIndexes[i]] if rq == 0: new_values[i] = val else: new_values[i] = Float.NaN if self.native_param == paramName: return new_values elif paramName == "brightnessTemp": return self.radianceToBrightnessTempSpectrum(new_values) elif paramName == "radiance": return self.brightnessTempToRadiance(new_values,channelIndex) #-override def channelIndexToFileIndex(self, dataset_index, idx): return self.fileIndexes[idx] #- override def getObsValidRange(self): if self.paramName == "radiance": return [0, 150] elif self.paramName == "brightnessTemp": return [180, 320] #- override def initMissing(self, longitudes, latitudes, values): class Missing: def __init__(self, longitudes, latitudes, values, xlength, ylength, invert): self.longitudes = longitudes self.latitudes = latitudes self.good_lines = [] nlines = xlength nelems = ylength self.invert = invert if invert == true: nlines = ylength nelems = xlength for i in xrange(nlines): cnt = 0 for j in xrange(nelems): if invert == true: ii = j + i*nelems else: ii = i + j*nlines lat = latitudes[ii] lon = longitudes[ii] val = values[ii] if ((lon <= 180 and lon >= -180) and (lat <= 90 and lat >= -90) and val > -9000): cnt += 1 if cnt == nelems: self.good_lines.append(i) self.n_good_lines = len(self.good_lines) self.nelems = nelems self.xlength = self.n_good_lines self.ylength = self.nelems if invert == true: self.ylength = self.n_good_lines self.xlength = self.nelems def dataToSegment(self, line, element): return self.good_lines[line], element def removeMissing(self, values): new_values = jarray.zeros(self.n_good_lines*self.nelems, 'f') for i in xrange(self.n_good_lines): for j in xrange(self.nelems): if self.invert == true: ii = j + i*self.nelems iii = j + self.good_lines[i]*self.nelems else: ii = i + j*self.n_good_lines iii = self.good_lines[i] + j*self.nlines new_values[ii] = values[iii] return new_values def getLengths(self): return self.xlength, self.ylength m = Missing(longitudes, latitudes, values, self.xlength, self.ylength, self.invert) return m def clone(self): return AirsData(self.filename,self.reader,self.readerTable,self.channelPropertyFile,master=self) #- Modis Airborne Simulator class MasData(MultiSpectralData): def __init__(self, filename, reader, readerTable,master=None): MultiSpectralData.__init__(self, filename, reader, readerTable, master) #-look for MODIS geo file in same directory as L1B if master == None: self.getScaleOffset() else: self.radScale = master.radScale def getScaleOffset(self): self.radScale, dum = self.reader.getObsScaleOffset(self.observationName[0], "scale_factor", "offset") def getImageRange(self, values, paramName, dataset_index, index): if paramName == "radiance": unit = Parser.parse("W/m^2/steradian/micrometer") self.imageRangeUnit = ScaledUnit(1, unit, "W/m^2/ster/\xb5m") return Helper.unpack(values, self.radScale[index], 0) elif paramName == "brightnessTemp": self.imageRangeUnit = Parser.parse("K") vals = Helper.unpack(values, self.radScale[index], 0) return self.radianceToBrightnessTemp(vals, index) def clone(self): return MasData(self.filename,self.reader,self.readerTable,master=self)