source: cows/trunk/cows/service/imps/data_reader_geoplot_backend/data_readers/csml_data_reader.py @ 5891

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/cows/trunk/cows/service/imps/data_reader_geoplot_backend/data_readers/csml_data_reader.py@5891
Revision 5891, 7.5 KB checked in by pnorton, 11 years ago (diff)

Fixed a small problem with the code that masks NAN values.

Line 
1
2import os
3
4import csml
5import cdms2 as cdms
6from copy import copy
7
8from cows.service.imps.csmlbackend.config import config
9from cows.service.imps.csmlbackend.csmlcommon import CSMLConnector
10from cows.service.imps.csmlbackend.wms_csmllayer import CSMLwmsDimension
11import numpy
12import logging
13
14log = logging.getLogger(__name__)
15
16#
17# Global CSML connector instance
18#
19globalCSMLConnector = CSMLConnector()
20
21
22class CSMLDataReader(object):
23   
24    def __init__(self, fileoruri):
25        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue()
26        self.connector = globalCSMLConnector
27        self.fileoruri = fileoruri
28        self.ds = self.connector.getCsmlDoc(fileoruri)
29        self.varcache = {} 
30       
31
32   
33    def getWMSLayerInfo(self):
34       
35        for feature in csml.csmllibs.csmlextra.listify(self.ds.featureCollection.featureMembers):
36            log.debug("feature.id = %s" % (feature.id,))
37            title, abstract, dimensions, units, crss=self._getWMSInfo(feature)
38            bb = self._getBBox(feature)
39            name = feature.id
40           
41            yield (name, title, abstract, dimensions, units, crss, bb) 
42
43   
44    def _getFeatureTitleAndAbstract(self, feature):
45        ''' given a csml feature, return basic info about the layer/feature/coverage
46        @return:   title, abstract'''
47
48        try:
49            title=feature.name.CONTENT
50        except:
51            title=''
52           
53        try:
54            abstract=feature.description.CONTENT
55        except:
56            abstract=title
57       
58        return title, abstract   
59
60    def _getBBox(self, feature):
61        try: 
62            bb = feature.getCSMLBoundingBox().getBox()
63        except:
64            #default to global
65            bb=[-180,-90,180,90]
66        return bb       
67
68    def _getWMSInfo(self, feature):
69        ''' given a csml feature, return info about the layer/feature
70        @return:   title, abstract, dimensions, units, crss '''
71
72        title, abstract = self._getFeatureTitleAndAbstract(feature)
73        units=feature.getDomainUnits() 
74        dimensions={}
75        tmpunits=copy(units)
76        tmpunits.reverse()
77        domain = feature.getDomain()
78       
79        for dim in feature.getAxisLabels():
80            nextdim=CSMLwmsDimension(domain, dim, tmpunits.pop())
81           
82            if dim not in ['latitude', 'longitude']:
83                dimensions[dim]=nextdim
84               
85        crs=feature.getNativeCRS()
86        crss=[self._crscat.getCRS(crs).twoD]
87       
88        if 'EPSG:4326' in crss:
89            crss.append('CRS:84')
90            crss.append('WGS84')
91                   
92        #the units to return are the units of measure.
93        try:
94            units=feature.value.rangeSet.valueArray.valueComponent.uom
95        except:
96            units='unknown units'
97           
98        return title, abstract, dimensions, units, crss
99
100    def _getFeature(self, id):
101        for feature in csml.csmllibs.csmlextra.listify(self.ds.featureCollection.featureMembers):
102            if feature.id == id:
103                return feature
104           
105        raise Exception("Feature with id %s not found" % (id,))
106 
107
108    def getNetcdfVar(self, featureId,  dimValues):
109        "Opens up the csml and retrieves the variable described by the dimensions"
110       
111
112        log.debug("featureId = %s, dimValues = %s" % (featureId, dimValues))
113       
114        dimList = list(dimValues.items())
115        dimList.sort()
116       
117        cacheKey = "%s:%s" % (featureId, dimList)
118
119        if cacheKey not in self.varcache:
120       
121            feature = self._getFeature(featureId)
122           
123            convertedDimVals = self._convertDimValues(dimValues)
124                 
125            if type(feature) == csml.parser.GridSeriesFeature:
126               
127                randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
128               
129                log.debug("getting csml feature")
130               
131                result= feature.subsetToGridSeries(config['tmpdir'], 
132                                            ncname=randomname, **convertedDimVals)
133             
134                #for now have to read netcdf back from
135                #disk (limitiation of CSML api)
136                netcdf=cdms.open(result[1])
137                                               
138                #and then delete the temporary file
139                os.system('rm %s'%result[1])
140               
141                log.debug("removed temp file %s" % (result[1],))
142               
143            else:
144                raise NotImplementedError
145           
146            variable =  netcdf(feature.id, squeeze=1)
147           
148            #try to set any NAN variable to masked variables
149            try:
150                #replace any NaN's with masked values
151                are_nan = numpy.isnan(variable)
152
153                if are_nan.any():
154                    # if the mask is just a single value we need to expand it
155                    if variable.mask.shape != variable.shape:
156                       
157                        if variable.mask:
158                            variable.mask = numpy.ones(variable.shape)
159                        else:
160                            variable.mask = numpy.zeros(variable.shape)
161                       
162                    variable[are_nan] = variable.getMissing()
163                    variable.mask[are_nan] = True
164           
165            except:
166                log.exception("Exception occurred while trying to fix NAN numbers in variable.")
167                raise
168           
169            self.varcache[cacheKey] = variable
170        else:
171            variable = self.varcache[cacheKey]
172
173        return variable
174   
175       
176    def _convertDimValues(self, dimValues):
177        """
178        Converts the string dimension values to floats (except for time values)
179        """
180        convertedVals = {}
181       
182        for dimval in dimValues:
183            if dimval != 'time':
184                convertedVals[dimval]=float(dimValues[dimval])
185            else:
186                #remove any trailing Zs from time string
187                if dimValues[dimval] [-1:] in ['Z', 'z']:
188                    convertedVals[dimval]=dimValues[dimval][:-1]
189                else:
190                    convertedVals[dimval] = dimValues[dimval] 
191
192        return convertedVals
193   
194   
195    def getConfigAxisXMLFile(self):
196       
197        xmlPath = None
198        for m in self._getMetadataElements():
199            log.debug("m.text = %s" % (m.text,))
200            metadataValue = m.text.strip()
201            if metadataValue.find('AxisConfigXML') == 0:
202                xmlPath = metadataValue.split('=')[1]
203       
204        log.debug("xmlPath = %s" % (xmlPath,))
205        return xmlPath
206   
207    def _getMetadataElements(self):
208       
209        featureCollectionElt = None
210        for c in self.ds.elem.getchildren(): 
211            if c.tag.find('CSMLFeatureCollection') > -1:
212                featureCollectionElt = c
213                break
214       
215        metadataElements = []
216        if featureCollectionElt != None:
217            for c in featureCollectionElt.getchildren():
218                if c.tag.find("metaDataProperty") > -1:
219                    metadataElements.append(c)   
220
221        return metadataElements
222       
223    @staticmethod
224    def isDataPresent(fileoruri):
225        global globalCSMLConnector
226       
227        log.debug("globalCSMLConnector.list() = %s" % ([x for x in globalCSMLConnector.list()],))
228        for file in globalCSMLConnector.list():
229            log.debug("file = %s" % (file,))
230            if file == fileoruri:
231                return True
232           
233        return False
234   
235   
Note: See TracBrowser for help on using the repository browser.