Changeset 3221


Ignore:
Timestamp:
15/01/08 15:21:07 (12 years ago)
Author:
domlowe
Message:

rewriting of code that determines CRS system

Location:
TI02-CSML/trunk/csml/csmllibs
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • TI02-CSML/trunk/csml/csmllibs/csmlbuilder.py

    r3144 r3221  
    4040        self.createFeatures() 
    4141        self.insertXlinks() 
    42         self.insertBoundingBoxes() 
     42        try: 
     43            self.insertBoundingBoxes() 
     44        except: 
     45            print 'warning, scanner could not calculate bounding boxes' 
    4346        self.saveFile() 
    4447        self.printToScreen() 
  • TI02-CSML/trunk/csml/csmllibs/csmlcrs.py

    r3144 r3221  
    22import sys      
    33import csml.csmllibs.csmlextra 
     4from copy import copy 
    45 
    56class CRSystem(object): 
     
    9495        crs.latAxis=None 
    9596        crs.timeAxis=1 
     97         
     98        # define british national grid plus time CRS: 
     99        crs=CRSystem(srsName='ndg:crs:bngt', axes =['eastings','northings','time']) 
     100        self.systems['ndg:crs:bngt']=crs 
     101        crs.lonAxis=0 
     102        crs.latAxis=1 
     103        crs.timeAxis=2 
     104         
    96105         
    97106        #define rotated_latitude_longitude CRS (test) 
     
    227236        return crs, axisorder 
    228237 
    229     def determineCRS(self, axes=None, units=None, knownCRSAxes=None,stdNames=None): 
    230         '''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. 
    231         e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east']) 
    232         returns self.systems['ndg:crs:xypt'], [3,2,1,0] 
    233         Assumes the units are ordered to correspond with the axes 
    234         '''      
    235         axisorder=[] 
    236         crs =None 
    237         crsMap=None 
    238          
    239         #stdNames is a list, but may be a list of Nones, in which case set = None. 
    240         if type(stdNames) is list: 
    241             namesok=False 
    242             for item in stdNames: 
    243                 if item is not None: 
    244                     namesok=True 
    245             if namesok is False: 
    246                 stdNames=None 
    247          
    248                         
    249                  
    250         if knownCRSAxes is not None: 
    251             for item in knownCRSAxes: 
    252                 if item is not None: 
    253                     crsMap=knownCRSAxes 
    254                     break          
     238     
     239    def _compareRealAxisNames(self, axes): 
     240        ''' will find complete exact matches only''' 
     241        crs = None 
     242        axisorder={} 
     243        unsortedaxes=copy(axes) 
     244        for system in self.systems: 
     245            crsaxes=self.systems[system].axes 
     246            unsortedcrsaxes=copy(crsaxes) 
     247            crsaxes.sort() 
     248            axes.sort() 
     249            if crsaxes==axes: 
     250                #found exact match  - stop checking 
     251                crs=self.systems[system] 
     252                axisorder=self._getAxisOrder(unsortedaxes,unsortedcrsaxes) 
     253                break 
     254        return crs, axisorder 
     255         
     256    def _getAxisOrder(self, axes, crsaxes): 
     257        #given a list of axis names and a list of comparable crsaxes returns an axis order list e.g. [3,2,1,0] where '3' means the first item in axes corresponds to the item at index 3 in the crsaxes list. 
     258        axisorder={} 
     259        for ax in axes: 
     260                if ax in crsaxes: 
     261                    axisorder[ax]=crsaxes.index(ax) 
     262        return axisorder 
     263     
     264    def _compareNamesAndUnits(self, axes, units, crsMap, stdNames=None): 
     265        #build a map of axes and units - use std names to help! 
     266        crs = None 
     267        axisorder={} 
    255268        if not crsMap: 
    256             #build a map of axes and units - use std names to help! 
    257             crsMap=[] 
    258269            for axis in axes: 
    259270                unittype = None 
     
    264275                    unit=units[axes.index(axis)] 
    265276                    unittype=self.getUnitType(unit)           
    266                 crsMap.append(unittype)               
     277                crsMap.append(unittype)  
     278             
    267279        match=0 
    268280        #now try and match up crsMap with known crsystems.         
     
    284296            if match == 1:  
    285297                crs=self.systems[system]               
    286                 break 
    287          
    288         if match==0:  #if couldn't find a match - generate a pseudo crs       
    289             crs, axisorder=self.__generatePseudoCRS(crsMap, axes, units) 
    290          
     298                break     
     299        return crs, axisorder, crsMap 
     300     
     301    def _addUnitsToCRSDefinition(self, crs, units, axisorder): 
     302        '''if no units are defined in the crs, set the units to be those of the original data''' 
    291303        if crs.units is None: 
    292             #if no units are defined, set the units to be those of the original data 
    293304            if units is not None: 
    294305                orderedunits=[] 
     
    296307                    orderedunits.append(units[axisorder[a]]) 
    297308                    crs.units=csml.csmllibs.csmlextra.stringify(orderedunits) 
     309        return crs 
     310                 
     311    def determineCRS(self, axes=None, units=None, knownCRSAxes=None,stdNames=None): 
     312        '''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. 
     313        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east']) 
     314        returns self.systems['ndg:crs:xypt'], {'latitude': 2, 'pressure': 1, 'longitude': 3, 'time': 0} 
     315        Assumes the units are ordered to correspond with the axes 
     316        '''  
     317        axisorder={} 
     318        crs =None 
     319        #perform various tests to try and determine which crs this is. 
     320         
     321        #test: compare real names of axes with names of axes in coordinate reference systems 
     322        #looks for exact matches only.  
     323        if axes is not None: 
     324            crs, axisorder = self._compareRealAxisNames(copy(axes)) 
     325         
     326        #set up crsMap - may be empty 
     327        crsMap=[] 
     328        if knownCRSAxes is not None: 
     329            for item in knownCRSAxes: 
     330                if item is not None: 
     331                    crsMap=knownCRSAxes 
     332                    break  
     333         
     334         
     335        if crs is None: 
     336            #try and match based on units, names and stdNames 
     337            crs,axisorder, crsMap=self._compareNamesAndUnits(axes, units, crsMap, stdNames=None) 
     338        if crs is None: 
     339            #try to generate a pseudo crs 
     340            crs, axisorder=self.__generatePseudoCRS(crsMap, axes, units) 
     341             
     342        #Finally, if no units are defined in the crs, set the units to be those of the original data 
     343        if crs is not None: 
     344            crs=self._addUnitsToCRSDefinition(crs, units, axisorder) 
    298345        return crs, axisorder 
    299      
     346       
    300347    def _identifyAxisByStdName(self,  stdName): 
    301348        '''given a CF standard name return the axis name''' 
     
    315362    crsdictionary={} 
    316363    f=cdms.open('/home/dom/svn/trunk/csml/testfiles/gridseries/xaaqda@pxs09c1.nc') 
     364    #f=cdms.open('/home/dom/EN/air_frost_1961-1990.nc') 
    317365    for var in f.variables: 
    318366        v=f.variables[var] 
     
    340388        print '----------------------------' 
    341389        crs, axisorder= cat.determineCRS(axes=axs, units=uns) 
    342         #print 'CRS: %s'%crs.srsName 
    343         print 'ORDER:  %s'%axisorder 
    344         print 'LABELS: %s'%crs.axisLabels 
    345         print 'CRS AXES: %s'%crs.axes 
     390        print crs, axisorder 
     391        if crs is not None: 
     392            print 'CRS found: %s'%crs.srsName 
     393            print 'ORDER:  %s'%axisorder 
     394            print 'LABELS: %s'%crs.axisLabels 
     395            print 'CRS AXES: %s'%crs.axes 
     396        else: 
     397            print 'CRS not found' 
    346398        labels='' 
    347399         
    348400         
    349         for axis in axes: 
    350             labels=labels + axis + ' ' 
    351          
    352         #print '\n Pseudo CSML:' 
    353         #print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels) 
    354         #print '<axisLabels>%s</axisLabels>'%labels 
    355401        #for axis in axes: 
    356             #print '<GridOrdinate>' 
    357             #unit=units[axes.index(axis)] 
    358             #print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]] 
    359             #print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis 
    360             #print '</GridOrdinate>' 
    361         #print '</GridSeriesDomain>' 
    362          
    363         #print '**************************************' 
    364    #test getting CRS by name: 
     402            #labels=labels + axis + ' ' 
     403         
     404        ##print '\n Pseudo CSML:' 
     405        ##print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels) 
     406        ##print '<axisLabels>%s</axisLabels>'%labels 
     407        ##for axis in axes: 
     408            ##print '<GridOrdinate>' 
     409            ##unit=units[axes.index(axis)] 
     410            ##print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]] 
     411            ##print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis 
     412            ##print '</GridOrdinate>' 
     413        ##print '</GridSeriesDomain>' 
     414         
     415        ##print '**************************************' 
     416   ##test getting CRS by name: 
    365417    #print 'getting ndg:crs:xypt by name' 
    366418    #crs=cat.getCRS('ndg:crs:xypt') 
Note: See TracChangeset for help on using the changeset viewer.