source: cows/trunk/cows/service/imps/data_reader_geoplot_backend/data_readers/modis_file_reader.py @ 5684

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

Fixed a problem in the masking of the landtype_common variable from the modis_file_reader. The index 0 data was bing masked along with the cases where all the indicies are 0.0.

Also made some small improvements to the geoplot layers.

Line 
1import cdms2 as cdms
2import numpy as N
3
4import logging
5
6from cows.service.imps.data_reader_geoplot_backend.data_readers.mosid_config import config
7
8log = logging.getLogger(__name__)
9
10class ModisFileReader(object):
11   
12    def __init__(self, fileoruri):
13       
14        if fileoruri != config['modis_variable']:
15            raise Exception("Unknown fileoruri value, fileoruri=%s, modis_var_name=%s" % (fileoruri, modis_var_name))
16       
17        self.fileoruri = fileoruri
18        fin = cdms.open(config['modis_file'])
19        self.var = fin(config['modis_variable'])
20        fin.close()
21       
22    def getWMSLayerInfo(self):
23       
24        name = 'landtype_common'
25        title = 'Common Land Type'
26        abstract = 'The most common land type in the area'
27        dimensions = {}
28        units = 'land usage index'
29        #!TODO: check the CRS
30        crss = ['CRS:84']
31        bb = self._getBBox()
32 
33        yield (name, title, abstract, dimensions, units, crss, bb) 
34       
35        for i in range(16):
36            name = "landtype_%s" % (i,)
37            title = 'Land Type index %s'  % (i,) 
38            abstract = 'The proportion of the area of land type index %s'  % (i,)
39           
40            yield (name, title, abstract, dimensions, units, crss, bb)
41           
42    def _getBBox(self):
43        try: 
44            latBounds = self.var.getLatitude().getBounds()
45            lonBounds = self.var.getLongitude().getBounds()
46           
47            bb = [lonBounds.min(), latBounds.min(), lonBounds.max(), latBounds.max()]
48        except:
49            log.exception("Exception occurred retriving bounds:")
50           
51            #default to global
52            bb=[-180,-90,180,90]
53        return bb       
54
55    def getNetcdfVar(self, name,  dimValues):
56        "Opens up the csml and retrieves the variable described by the dimensions"
57       
58        convertedDimVals = self._convertDimValues(dimValues)
59       
60        if name == 'landtype_common':
61            var = self._getLandtypeCommon(convertedDimVals)
62       
63        elif name.find('landtype_') == 0:
64           
65           
66            indexStr = name[len('landtype_'):]
67           
68            try:
69                index = int(indexStr)
70            except:
71                raise Exception("Error converting %s to an int, name = %s" % (indexStr, name))
72           
73            var = self._getVarNumber(convertedDimVals, index)
74                 
75       
76        else:
77            raise Exception("Unknown layer name %s" % (name,))
78       
79     
80        latMin = var.getLatitude().getBounds()[0].min()
81        latMax =  var.getLatitude().getBounds()[-1].max()
82        lonMin = var.getLongitude().getBounds()[0].min()
83        lonMax =  var.getLongitude().getBounds()[-1].max()
84        nLat = len(var.getLatitude())
85        nLon = len(var.getLongitude())
86       
87        log.debug("var: latMin = %s, latMax = %s, lonMin = %s, lonMax = %s, nLat = %s, nLon = %s" % (latMin, latMax, lonMin, lonMax, nLat, nLon))
88        log.debug("var: getAxisIds = %s, shape=%s" % (var.getAxisIds(), var[:].shape))   
89               
90               
91        return var
92   
93    def _convertDimValues(self, dimValues):
94        """
95        Converts the string dimension values to floats (except for time values)
96        """
97        convertedVals = {}
98       
99        for dimval in dimValues:
100            if dimval != 'time':
101                convertedVals[dimval]=float(dimValues[dimval])
102            else:
103                #remove any trailing Zs from time string
104                if dimValues[dimval] [-1:] in ['Z', 'z']:
105                    convertedVals[dimval]=dimValues[dimval][:-1]
106                else:
107                    convertedVals[dimval] = dimValues[dimval] 
108
109        return convertedVals
110   
111    def _getLandtypeCommon(self, convertedDimValues):
112        # the problem is we need to know when all the land types are 0
113        # and put that as a mask, otherwise we need to use the index of the
114        # maxiumum
115       
116        (nx, ny, nz) = self.var.shape
117        log.debug("nx = %s, ny=%s, nz=%s" % (nx, ny, nz))
118       
119        #assuming the axis are lat, lon index, sum over the index axis
120        sumArr = self.var[:].sum(axis=2)
121       
122        #create a mask where this sum is 0.0
123        mask = sumArr == 0.0
124       
125        #get the index of the maximum for the actual data ( should be 0-15)
126        arr = self.var[:].argmax(2)
127       
128        common_var = cdms.createVariable(arr, axes=self.var.getAxisList()[0:2],
129                                         mask=mask)
130       
131       
132        return common_var
133       
134    def _getVarNumber(self, convertedDimValues, varIndex):
135        var_n = self.var(index=varIndex, squeeze=True)
136        var_n.set_fill_value(0.0)
137        return var_n
138       
139    @staticmethod
140    def isDataPresent(fileoruri):
141       
142        if config['modis_variable'] == None:
143            return False
144        else:
145            return fileoruri == config['modis_variable']
146
147
148if __name__ == '__main__':
149   
150    import logging
151    logging.basicConfig(level=logging.DEBUG, format='%(name)-10s %(asctime)s ln:%(lineno)-3s %(levelname)-8s\n %(message)s\n')
152   
153    config['modis_variable'] = 'landtype'
154    config['modis_file'] = '/data/pnorton/qesdi_data/modis/modis_data/MCD12Q1_10min/MODIS_classification_IGBP.nc'
155   
156   
157    reader = ModisFileReader('landtype')
158   
159    var = reader._getLandtypeCommon({})
160   
161    print var.id
Note: See TracBrowser for help on using the repository browser.