pyModis/pymodis/qualitymodis.py

#!/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]])
})


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

    def loadData(self):
        """loads the input file to the object"""
        os.path.isfile(self.infile)
        self.ds = gdal.Open(self.infile)

    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

    def setProductGroup(self):
        """read productGroup from Metadata of hdf file"""
        self.productGroup = self.productType[3:5]

    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!")

    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]

    def loadQAArray(self):
        """loads the QA layer to the object"""
        self.qaArray = gdal_array.LoadFile(self.qaLayer)

    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)

    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

    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.")