Source code for pymodis.qualitymodis

#!/usr/bin/env python
#  class to convert/process modis data
#
#  (c) Copyright Ingmar Nitze 2013
#  Authors: Ingmar Nitze, Luca Delucchi
#  Email: initze at ucc dot ie
#  Email: luca dot delucchi at iasma dot it
#
##################################################################
#
#  This MODIS Python class is licensed under the terms of GNU GPL 2.
#  This program is free software; you can redistribute it and/or
#  modify it under the terms of the GNU General Public License as
#  published by the Free Software Foundation; either version 2 of
#  the License, or (at your option) any later version.
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#  See the GNU General Public License for more details.
#
##################################################################
"""A class for the extraction and transformation of MODIS quality layers to
specific information

Classes:

* :class:`QualityModis`

"""

# python 2 and 3 compatibility
from __future__ import print_function
from builtins import dict

import os
try:
    import numpy as np
except ImportError:
    raise ImportError('Numpy library not found, please install it')

try:
    import osgeo.gdal as gdal
    import osgeo.gdal_array as gdal_array
except ImportError:
    try:
        import gdal
        import gdal_array
    except ImportError:
        raise ImportError('Python GDAL library not found, please install '
                          'python-gdal')


VALIDTYPES = dict({'13': list(map(str, list(range(1, 10)))), '11': list(map(str, list(range(1, 6))))})

PRODUCTPROPS = dict({
    'MOD13Q1': ([2], ['QAGrp1']),
    'MYD13Q1': ([2], ['QAGrp1']),
    'MOD13A1': ([2], ['QAGrp1']),
    'MYD13A1': ([2], ['QAGrp1']),
    'MOD13A2': ([2], ['QAGrp1']),
    'MYD13A2': ([2], ['QAGrp1']),
    'MOD13A3': ([2], ['QAGrp1']),
    'MYD13A3': ([2], ['QAGrp1']),
    'MOD13C1': ([2], ['QAGrp1']),
    'MYD13C1': ([2], ['QAGrp1']),
    'MOD13C2': ([2], ['QAGrp1']),
    'MYD13C2': ([2], ['QAGrp1']),
    'MOD11A1': ([1, 5], ['QAGrp2', 'QAGrp2']),
    'MYD11A1': ([1, 5], ['QAGrp2', 'QAGrp2']),
    'MOD11A2': ([1, 5], ['QAGrp4', 'QAGrp4']),
    'MYD11A2': ([1, 5], ['QAGrp4', 'QAGrp4']),
    'MOD11B1': ([1, 5, -2], ['QAGrp2', 'QAGrp2', 'QAGrp3']),
    'MYD11B1': ([1, 5, -2], ['QAGrp2', 'QAGrp2', 'QAGrp3']),
    'MOD11C1': ([1, 5, -2], ['QAGrp2', 'QAGrp2', 'QAGrp3']),
    'MYD11C1': ([1, 5, -2], ['QAGrp2', 'QAGrp2', 'QAGrp3']),
    'MOD11C2': ([1, 6], ['QAGrp2', 'QAGrp2']),
    'MYD11C2': ([1, 6], ['QAGrp2', 'QAGrp2']),
    'MOD11C3': ([1, 6], ['QAGrp2', 'QAGrp2']),
    'MYD11C3': ([1, 6], ['QAGrp2', 'QAGrp2'])
})


QAindices = dict({
    'QAGrp1': (16, [[-2, None], [-6, -2], [-8, -6], [-9, -8],
               [-10, -9], [-11, -10], [-14, -11], [-15, -14],
               [-16, -15]]),
    'QAGrp2': (7, [[-2, None], [-3, -2], [-4, -3], [-6, -4],
               [-8, -6]]),
    'QAGrp3': (7, [[-3, None], [-6, -3], [-7, -6]]),
    'QAGrp4': (8, [[-2, None], [-4, -2], [-6, -4], [-8, -6]])
})


[docs]class QualityModis(): """A Class for the extraction and transformation of MODIS quality layers to specific information :param str infile: the full path to the hdf file :param str outfile: the full path to the parameter file """ def __init__(self, infile, outfile, qType=None, qLayer=None, pType=None): """Function to initialize the object""" self.infile = infile self.outfile = outfile self.qType = qType self.qLayer = qLayer self.qaGroup = None self.pType = pType
[docs] def loadData(self): """loads the input file to the object""" os.path.isfile(self.infile) self.ds = gdal.Open(self.infile)
[docs] def setProductType(self): """read productType from Metadata of hdf file""" if self.pType == None: self.productType = self.ds.GetMetadata()['SHORTNAME'] else: self.productType = self.pType
[docs] def setProductGroup(self): """read productGroup from Metadata of hdf file""" self.productGroup = self.productType[3:5]
[docs] def setQAGroup(self): """set QA dataset group type""" if self.productType in list(PRODUCTPROPS.keys()): self.qaGroup = PRODUCTPROPS[self.productType][1][int(self.qLayer)-1] else: print("Product version is currently not supported!")
[docs] def setQALayer(self): """function sets the input path of the designated QA layer""" self.qaLayer = self.ds.GetSubDatasets()[PRODUCTPROPS[self.productType][0][int(self.qLayer)-1]][0]
[docs] def loadQAArray(self): """loads the QA layer to the object""" self.qaArray = gdal_array.LoadFile(self.qaLayer)
[docs] def qualityConvert(self, modisQaValue): """converts encoded Bit-Field values to designated QA information""" startindex = QAindices[self.qaGroup][1][int(self.qType)-1][0] endindex = QAindices[self.qaGroup][1][int(self.qType)-1][1] return int(np.binary_repr(modisQaValue, QAindices[self.qaGroup][0])[startindex: endindex], 2)
[docs] def exportData(self): """writes calculated QA values to physical .tif file""" qaDS = gdal.Open(self.qaLayer) dr = gdal.GetDriverByName('GTiff') outds = dr.Create(self.outfile, self.ncols, self.nrows, 1, gdal.GDT_Byte) outds.SetProjection(qaDS.GetProjection()) outds.SetGeoTransform(qaDS.GetGeoTransform()) outds.GetRasterBand(1).WriteArray(self.qaOut) outds = None qaDS = None
[docs] def run(self): """Function defines the entire process""" self.loadData() self.setProductType() self.setProductGroup() #self.setDSversion() self.setQAGroup() self.setQALayer() self.loadQAArray() self.nrows, self.ncols = self.qaArray.shape print("Conversion started !") self.qaOut = np.zeros_like(self.qaArray, dtype=np.int8) if self.productGroup in ['11', '13'] and self.qType in VALIDTYPES[self.productGroup] and self.qaGroup != None: for val in np.unique(self.qaArray): ind = np.where(self.qaArray == val) self.qaOut[ind] = self.qualityConvert(self.qaArray[ind][0]) self.exportData() print("Export finished!") else: print("This MODIS type is currently not supported.")