source: TI02-CSML/trunk/csml/csmllibs/csmlcrs.py @ 1979

import string
import cdms # just while testing

class CRSystem(object):
    def __init__(self, srsName, axes):
        self.srsName=srsName
        self.axes=axes
        self.srsDimension=len(self.axes)
        self.axisLabels=''
        for axis in self.axes:
            self.axisLabels=self.axisLabels + axis + ' '
        self.timeAxis=None
        self.lonAxis=None
        self.latAxis=None

class CRSCatalogue(object):
    def __init__(self):
        # dictionary to hold known CRSystems:
        self.systems={}
        
        # define lon lat pressure time CRS:        
        crs=CRSystem(srsName='ndg:crs:xypt', axes =['Lon', 'Lat','Pressure','Time'])
        crs.lonAxis=0
        crs.latAxis=1
        crs.timeAxis=3
        self.systems['ndg:crs:xypt']=crs
        
        # define lon lat height time CRS:        
        crs=CRSystem(srsName='ndg:crs:xyht', axes =['Lon', 'Lat','Height','Time'])
        crs.lonAxis=0
        crs.latAxis=1
        crs.timeAxis=3
        self.systems['ndg:crs:xyht']=crs
        
        # define lon lat time CRS:
        crs=CRSystem(srsName='ndg:crs:xyt', axes =['Lon', 'Lat','Time'])
        self.systems['ndg:crs:xyt']=crs
        
        #define unknown 1D CRS:
        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
        self.systems['ndg:crs:unknown']=crs
        
        #define unknown 3D CRS:
        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
        self.systems['ndg:crs:unknown2d']=crs
        
        #define unknown 3D CRS:
        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
        self.systems['ndg:crs:unknown3d']=crs
        
        #define unknown 4D CRS:
        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
        self.systems['ndg:crs:unknown4d']=crs
        
            
    def getCRS(self, srsName):
        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object 
        try:
            return self.systems[srsName]
        except:
            return self.systems['ndg:crs:unknown']
    
    def determineCRS(self, axes, units):
        '''given any list of axis names and a list of units for these axes attempt to determine the CRS and return the CRSystem object and an axis order list.
        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
        returns self.systems['ndg:crs:xypt'], [3,2,1,0]
        Assumes the units are ordered to correspond with the axes!
        '''
        axisorder=[]
        if len(axes)==2:
            crs=self.systems['ndg:crs:unknown2d']
        if len(axes)==3:
            crs=self.systems['ndg:crs:unknown3d']
        elif len(axes)==4:
            crs=self.systems['ndg:crs:unknown4d']
        else :
            crs=self.systems['ndg:crs:unknown']
        crsMap=[]
        #this can be extended to accomodate more units
        for axis in axes:
            unit=units[axes.index(axis)]
            unittype='unknown'
            if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
                unittype='Time'
            elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc']:
                unittype='Time'
            elif string.lower(unit) in ['mbar', 'pa','level']:
                unittype='Pressure'
            elif string.lower(unit) in ['m', 'km']:
                unittype='Height'
            elif string.lower(unit) in ['degrees_north']:
                unittype='Lat'
            elif string.lower(unit) in ['degrees_east']:
                unittype='Lon'
            crsMap.append(unittype)
        
        #now try and match up crsMap with known crsystems.        
        #axisorder stores the order of the supplied axes against the known crs.
        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
        #axisorder would be [2, 0, 1]  
        
        for system in self.systems:
            if len(self.systems[system].axes) != len(crsMap):
                #wrong dimensionality, skip this crs
                continue
            match=0
            axisorder=[]
            for ax in self.systems[system].axes:
                if ax in crsMap:
                    axisorder.append(crsMap.index(ax))
                else:
                    match = 1
                    break
            if match == 0: 
                crs=self.systems[system]
                break
                        
        return crs, axisorder
        
def main():
    cat=CRSCatalogue()
  

    #test every crs in this netcdf file
    crsdictionary={}
    f=cdms.open('/home/dom/Desktop/SVN/trunk/csml/testfiles/gridseries/xaaqda@pxs19c1.nc')
    for var in f.variables:
        v=f.variables[var]
        axes=v.getAxisIds()
        units =[]
        for axis in axes:
            ax=f.getAxis(axis)
            units.append(ax.attributes['units'])
        crsdictionary[str(axes)]=units
    n=2
    #test each definition in dictionary
    for key in crsdictionary:
        n=n+1
        print '\n TEST %d'%n
        print '*************************************'
       #make list from string.
        axs=[]
        for item in key.split():
            axs.append(item)
        uns=crsdictionary[key]
        print 'AXES: %s'%axs
        print 'UNITS: %s'%uns
        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
        print 'CRS: %s'%crs.srsName
        print 'ORDER:  %s'%axisorder
        #print 'axis order =%s'%axisorder
        labels=''
        
        print axes
        # ['t', 'ht', 'latitude_1', 'longitude_1']
        print crs.axes
         # ['Lon', 'Lat', 'Pressure', 'Time']
        print axisorder
        # [3, 2, 1, 0]
        
        for axis in axes:
            labels=labels + axis + ' '
        
        print '\n Pseudo CSML:'
        print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
        print '<axisLabels>%s</axisLabels>'%labels
        for axis in axes:
            print '<GridOrdinate>'
 
            unit=units[axes.index(axis)]
            print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
            print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
            print '</GridOrdinate>'
        print '</GridSeriesDomain>'
        
        print '**************************************'
   #test getting CRS by name:
    #print 'getting ndg:crs:xypt by name'
    #crs=cat.getCRS('ndg:crs:xypt')
    #print 'returns %s'%crs.srsName
       
    
if __name__=="__main__":
    main()