stack.py

"""Classes for HDF5/MintPy file operations."""
############################################################
# Program is part of MintPy                                #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
# Author: Zhang Yunjun, Heresh Fattahi, 2018               #
############################################################
# Recommend import:
#   from mintpy.objects import timeseries, ifgramStack, geometry


import datetime as dt
import itertools
import os
import time

import h5py
import numpy as np
from scipy import ndimage

from mintpy.utils import ptime

##------------------ Global Variables ---------------------##

DATA_TYPE_DICT = {
    'bool': np.bool_, 'byte': np.bool_, 'flag': np.bool_,
    'int': np.int16, 'int16': np.int16, 'short': np.int16, 'int32': np.int32,
    'int64': np.int64, 'long': np.int64,
    'float': np.float32, 'float32': np.float32,
    'float_': np.float64, 'float64': np.float64,
    'complex': np.complex64, 'complex64': np.complex64, 'cpx_float32': np.complex64,
    'cfloat': np.complex64, 'cfloat32': np.complex64,
    'complex128': np.complex128, 'complex_': np.complex128, 'cpx_float64': np.complex128,
}

TIMESERIES_KEY_NAMES = ['timeseries', 'HDFEOS', 'giantTimeseries']

TIMESERIES_DSET_NAMES = [
    'timeseries',
    'raw',
    'troposphericDelay',
    'topographicResidual',
    'ramp',
    'displacement',
]

GEOMETRY_DSET_NAMES = [
    # coordinates
    'height',
    'latitude',
    'longitude',
    'rangeCoord',
    'azimuthCoord',
    # others
    'incidenceAngle',
    'azimuthAngle',
    'slantRangeDistance',
    'shadowMask',
    'waterMask',
    'commonMask',
    'bperp',
]

IFGRAM_DSET_NAMES = [
    # interferogram
    'unwrapPhase',
    'unwrapPhase_bridging_phaseClosure',
    'unwrapPhase_bridging',
    'unwrapPhase_phaseClosure',
    'coherence',
    'connectComponent',
    'wrapPhase',
    'magnitude',
    # offset
    'azimuthOffset',
    'azimuthOffsetStd',
    'rangeOffset',
    'rangeOffsetStd',
    'offsetSNR',
    'refPhase',
]

DSET_UNIT_DICT = {
    # interferogram
    'unwrapPhase'      : 'radian',
    'coherence'        : '1',
    'connectComponent' : '1',
    'wrapPhase'        : 'radian',
    'magnitude'        : '1',

    # offset
    'azimuthOffset'    : 'pixel',
    'azimuthOffsetStd' : 'pixel',
    'rangeOffset'      : 'pixel',
    'rangeOffsetStd'   : 'pixel',
    'offsetSNR'        : '1',

    # geometry
    'height'             : 'm',
    'latitude'           : 'degree',
    'longitude'          : 'degree',
    'rangeCoord'         : '1',
    'azimuthCoord'       : '1',
    'incidenceAngle'     : 'degree',
    'azimuthAngle'       : 'degree',
    'slantRangeDistance' : 'm',
    'shadowMask'         : '1',
    'waterMask'          : '1',
    'commonMask'         : '1',
    'bperp'              : 'm',

    # time-series
    'timeseries'          : 'm',
    'raw'                 : 'm',
    'troposphericDelay'   : 'm',
    'topographicResidual' : 'm',
    'ramp'                : 'm',
    'displacement'        : 'm',

    # others
    'temporalCoherence' : '1',
    'velocity'          : 'm/year',
    'acceleration'      : 'm/year^2',
    'mask'              : '1',

    # giant
    'giantTimeseries' : 'mm',
    'recons'          : 'mm',
    'rawts'           : 'mm',
    'sar_aps'         : 'mm',
    'igram_aps'       : 'mm',
    'figram'          : 'mm',
    'igram'           : 'mm',
    'cmask'           : '1',
    'ifgcnt'          : '1',

    # binary
    'unw'       : 'radian',
    'int'       : 'radian',
    'flat'      : 'radian',
    'cor'       : '1',
    'dem'       : 'm',
    'hgt'       : 'm',
    'hgt_sim'   : 'm',
    'intensity' : '1',
}



################################ timeseries class begin ################################
class timeseries:
    """
    Time-series object for displacement of a set of SAR images from the same platform and track.
    It contains three datasets in root level: date, bperp and timeseries.

    File structure: https://mintpy.readthedocs.io/en/latest/api/data_structure/#timeseries
    """

    def __init__(self, file=None):
        self.file = file
        self.name = 'timeseries'

    def close(self, print_msg=True):
        try:
            self.f.close()
            if print_msg:
                print(f'close timeseries file: {os.path.basename(self.file)}')
        except:
            pass
        return None

    def open(self, print_msg=True):
        if print_msg:
            print(f'open {self.name} file: {os.path.basename(self.file)}')
        self.get_metadata()
        self.get_size()
        self.get_date_list()
        self.numPixel = self.length * self.width

        with h5py.File(self.file, 'r') as f:
            try:
                self.pbase = f['bperp'][:]
                self.pbase -= self.pbase[self.refIndex]
            except:
                self.pbase = None

        # time info
        self.dateFormat = ptime.get_date_str_format(self.dateList[0])
        self.times = np.array([dt.datetime.strptime(i, self.dateFormat) for i in self.dateList])
        # add hh/mm/ss info to the datetime objects
        if 'T' not in self.dateFormat or all(i.hour==0 and i.minute==0 for i in self.times):
            if 'CENTER_LINE_UTC' in self.metadata.keys():
                utc_sec = float(self.metadata['CENTER_LINE_UTC'])
                self.times = np.array([i + dt.timedelta(seconds=utc_sec) for i in self.times])
        self.tbase = np.array([(i.days + i.seconds / (24 * 60 * 60))
                               for i in (self.times - self.times[self.refIndex])],
                              dtype=np.float32)

        # list of float for year, 2014.95
        self.yearList = [i.year + (i.timetuple().tm_yday-1)/365.25 for i in self.times]
        self.sliceList = [f'{self.name}-{i}' for i in self.dateList]
        return None

    def get_metadata(self):
        with h5py.File(self.file, 'r') as f:
            self.metadata = dict(f.attrs)
            dates = f['date'][:]
        for key, value in self.metadata.items():
            try:
                self.metadata[key] = value.decode('utf8')
            except:
                self.metadata[key] = value

        # ref_date/index
        dateList = [i.decode('utf8') for i in dates]
        if 'REF_DATE' not in self.metadata.keys():
            self.metadata['REF_DATE'] = dateList[0]
        self.refIndex = dateList.index(self.metadata['REF_DATE'])
        self.metadata['START_DATE'] = dateList[0]
        self.metadata['END_DATE'] = dateList[-1]
        return self.metadata

    def get_size(self):
        with h5py.File(self.file, 'r') as f:
            self.numDate, self.length, self.width = f[self.name].shape[-3:]
        return self.numDate, self.length, self.width

    def get_date_list(self):
        with h5py.File(self.file, 'r') as f:
            self.dateList = [i.decode('utf8') for i in f['date'][:]]
        return self.dateList

    def read(self, datasetName=None, box=None, squeeze=True, print_msg=True):
        """Read dataset from timeseries file
        Parameters: self : timeseries object
                    datasetName : (list of) string in YYYYMMDD format
                    box : tuple of 4 int, indicating x0,y0,x1,y1 of range
        Returns:    data : 2D or 3D dataset
        Examples:   from mintpy.objects import timeseries
                    tsobj = timeseries('timeseries_ERA5_demErr.h5')
                    data = tsobj.read(datasetName='20161020')
                    data = tsobj.read(datasetName='20161020', box=(100,300,500,800))
                    data = tsobj.read(datasetName=['20161020','20161026','20161101'])
                    data = tsobj.read(box=(100,300,500,800))
        """
        if print_msg:
            print(f'reading {self.name} data from file: {self.file} ...')
        self.open(print_msg=False)

        # convert input datasetName into list of dates
        if not datasetName or datasetName == 'timeseries':
            datasetName = []
        elif isinstance(datasetName, str):
            datasetName = [datasetName]
        datasetName = [i.replace('timeseries', '').replace('-', '') for i in datasetName]

        with h5py.File(self.file, 'r') as f:
            ds = f[self.name]
            if isinstance(ds, h5py.Group):  # support for old mintpy files
                ds = ds[self.name]

            # Get dateFlag - mark in time/1st dimension
            dateFlag = np.zeros((self.numDate), dtype=np.bool_)
            if not datasetName:
                dateFlag[:] = True
            else:
                for e in datasetName:
                    dateFlag[self.dateList.index(e)] = True

            # Get Index in space/2_3 dimension
            if box is None:
                box = [0, 0, self.width, self.length]
            xsize = box[2] - box[0]
            ysize = box[3] - box[1]

            # read
            num_slice = np.sum(dateFlag)
            inds = np.where(dateFlag)[0].tolist()

            if num_slice / dateFlag.size < 0.05:
                # single indexing if only a small fraction is read
                data = np.zeros((num_slice, ysize, xsize), dtype=ds.dtype)
                for i, ind in enumerate(inds):
                    data[i] = ds[ind,
                                 box[1]:box[3],
                                 box[0]:box[2]]
            else:
                data = ds[:,
                          box[1]:box[3],
                          box[0]:box[2]][dateFlag]

            if squeeze and any(i == 1 for i in data.shape):
                data = np.squeeze(data)
        return data


    def write2hdf5(self, data, outFile=None, dates=None, bperp=None, metadata=None, refFile=None, compression=None):
        """
        Parameters: data  : 3D array of float32
                    dates : 1D array/list of string in YYYYMMDD format
                    bperp : 1D array/list of float32 (optional)
                    metadata : dict
                    outFile : string
                    refFile : string
                    compression : string or None
        Returns: outFile : string
        Examples:
            from mintpy.objects import timeseries

            ##Generate a new timeseries file
            tsobj = timeseries('timeseries.h5')
            timeseries.write(data, dates=dateList, bperp=bperp, metadata=atr)

            ##Generate a timeseries with same attributes and same date/bperp info
            tsobj = timeseries('timeseries_demErr.h5')
            timeseries.write(data, refFile='timeseries.h5')
        """
        if not outFile:
            outFile = self.file
        if refFile:
            refobj = timeseries(refFile)
            refobj.open(print_msg=False)
            if metadata is None:
                metadata = refobj.metadata
            if dates is None:
                dates = refobj.dateList
            if bperp is None:
                bperp = refobj.pbase
            # get ref file compression type if input compression is None
            if compression is None:
                with h5py.File(refFile, 'r') as rf:
                    compression = rf[TIMESERIES_DSET_NAMES[0]].compression
            refobj.close(print_msg=False)
        data = np.array(data, dtype=np.float32)
        dates = np.array(dates, dtype=np.bytes_)
        bperp = np.array(bperp, dtype=np.float32)
        metadata = dict(metadata)
        metadata['FILE_TYPE'] = self.name

        # directory
        outDir = os.path.dirname(os.path.abspath(outFile))
        if not os.path.isdir(outDir):
            os.makedirs(outDir)
            print(f'create directory: {outDir}')

        # 3D dataset - timeseries
        print(f'create timeseries HDF5 file: {outFile} with w mode')
        with h5py.File(outFile, 'w') as f:
            print(('create dataset /timeseries of {t:<10} in size of {s} '
                   'with compression={c}').format(t=str(data.dtype),
                                                  s=data.shape,
                                                  c=compression))
            f.create_dataset('timeseries',
                             data=data,
                             chunks=True,
                             compression=compression)

            # 1D dataset - date / bperp
            print(f'create dataset /dates      of {str(dates.dtype):<10} in size of {dates.shape}')
            f.create_dataset('date', data=dates)

            if bperp.shape != ():
                print(f'create dataset /bperp      of {str(bperp.dtype):<10} in size of {bperp.shape}')
                f.create_dataset('bperp', data=bperp)

            # Attributes
            for key, value in metadata.items():
                f.attrs[key] = str(value)
        print(f'finished writing to {outFile}')
        return outFile

    def timeseries_std(self, maskFile=None, outFile=None):
        """Calculate the standard deviation (STD) for acquisition of time-series,
           output result to a text file.
        """
        # Calculate STD
        data = self.read()
        if maskFile:
            mask = singleDataset(maskFile).read()
            print('read mask from file: '+maskFile)
            data[:, mask == 0] = np.nan
        self.std = np.nanstd(data, axis=(1, 2))

        # Write text file
        header = 'Standard Deviation in space for each acquisition of time-series\n'
        header += f'Timeseries file: {self.file}\n'
        header += f'Mask file: {maskFile}\n'
        header += 'Date\t\tSTD (m)'
        if not outFile:
            outFile = os.path.join(os.path.dirname(os.path.abspath(self.file)),
                                   f'std_{os.path.splitext(os.path.basename(self.file))[0]}.txt')
        np.savetxt(outFile, np.hstack((np.array(self.dateList).reshape(-1, 1), self.std.reshape(-1, 1))),
                   fmt='%s', delimiter='\t', header=header)
        print(f'save timeseries STD to text file: {outFile}')
        return outFile

    def timeseries_rms(self, maskFile=None, outFile=None):
        """Calculate the Root Mean Square for each acquisition of time-series
            and output result to a text file.
        """
        # Get date list
        date_list = self.get_date_list()
        num_date  = len(date_list)

        # Get mask
        if maskFile and os.path.isfile(maskFile):
            print('read mask from file: '+maskFile)
            mask = singleDataset(maskFile).read()

        # Calculate RMS one date at a time
        self.rms = np.zeros(num_date) * np.nan
        print(f'reading {self.name} data from file: {self.file} ...')
        prog_bar = ptime.progressBar(maxValue=num_date)
        for i in range(num_date):
            data = self.read(datasetName=f'{date_list[i]}', print_msg=False)
            if maskFile and os.path.isfile(maskFile):
                data[mask == 0] = np.nan
            self.rms[i] = np.sqrt(np.nanmean(np.square(data), axis=(0, 1)))
            prog_bar.update(i+1, suffix=f'{i+1}/{num_date}')
        prog_bar.close()

        # Write text file
        header = 'Root Mean Square in space for each acquisition of time-series\n'
        header += f'Timeseries file: {self.file}\n'
        header += f'Mask file: {maskFile}\n'
        header += 'Date\t\tRMS (m)'
        if not outFile:
            outFile = os.path.join(os.path.dirname(os.path.abspath(self.file)),
                                   f'rms_{os.path.splitext(os.path.basename(self.file))[0]}.txt')
        np.savetxt(outFile, np.hstack((np.array(self.dateList).reshape(-1, 1), self.rms.reshape(-1, 1))),
                   fmt='%s', delimiter='\t', header=header)
        print(f'save timeseries RMS to text file: {outFile}')
        return outFile

    def spatial_average(self, maskFile=None, box=None, reverseMask=False, threshold=None):
        self.open(print_msg=False)
        data = self.read(box=box)
        if maskFile and os.path.isfile(maskFile):
            print('read mask from file: '+maskFile)
            mask = singleDataset(maskFile).read(box=box)
            data[:, mask == int(reverseMask)] = np.nan

        # calculate area ratio if threshold is specified
        # percentage of pixels with value above the threshold
        if threshold is not None:
            data[data > threshold] = 1
            data[data <= threshold] = 0

        dmean = np.nanmean(data, axis=(1, 2))
        return dmean, self.dateList

    def temporal_average(self):
        print(f'calculating the temporal average of timeseries file: {self.file}')
        self.open(print_msg=False)
        data = self.read(squeeze=False)
        dmean = np.nanmean(data, axis=0)
        return dmean


    def temporal_derivative(self, out_file):
        print(f'calculating the temporal derivative of timeseries file: {self.file}')

        # read
        print('reading timeseries data')
        self.open(print_msg=False)
        ts_data = self.read(print_msg=False)

        # calculate
        print('calculate the 1st derivative of timeseries data')
        ts_data_1d = np.zeros(ts_data.shape, np.float32)
        ts_data_1d[1:, :, :] = np.diff(ts_data, n=1, axis=0)

        # write
        if not out_file:
            fbase = os.path.splitext(self.file)[0]
            out_file = f'{fbase}_1stDiff.h5'
        self.write2hdf5(ts_data_1d, outFile=out_file, refFile=self.file)

        return out_file


    def temporal_filter(self, time_win=1.0, filter_type='guassian', out_file=None):
        """Filter the time-series in time with a moving window.

        Parameters: time_win    - float, sigma of Gaussian distribution in time in months (30.4 days)
                    filter_type - str, filter type: gaussian, median
                    out_file    - str, output file name
        Returns:    out_file    - str, output file name
        """
        # use integer type if possible for shorter default output file name
        time_win = int(time_win) if float(time_win).is_integer() else time_win

        # output file
        if not out_file:
            fbase = os.path.splitext(self.file)[0]
            out_file = f'{fbase}_temp{filter_type.capitalize()}{time_win}.h5'
        print(f'output file: {out_file}')

        # read
        self.open()
        ts_data = self.read().reshape(self.numDate, -1)

        print('-'*50)
        ts_data_filt = np.zeros(ts_data.shape, np.float32)

        if filter_type == 'gaussian':
            print(f'temporal filtering via a Gaussian window of {time_win} months')
            tbase = self.tbase.reshape(-1, 1) / 365.25 * 12  # months (30.4 days)
            prog_bar = ptime.progressBar(maxValue=self.numDate)
            for i in range(self.numDate):

                # calc weight from Gaussian (normal) distribution in time
                tbase_diff = tbase[i] - tbase
                weight = np.exp(-0.5 * (tbase_diff**2) / (time_win**2))
                weight /= np.sum(weight)

                # smooth the current acquisition via Gaussian weighting
                ts_data_filt[i, :] = np.sum(ts_data * weight, axis=0)

                prog_bar.update(i+1, suffix=self.dateList[i])
            prog_bar.close()

        elif filter_type == 'median':
            print(f'temporal filtering via scipy.ndimage.median_filter of {time_win} acquisitions')
            ts_data_filt = ndimage.median_filter(
                ts_data,
                size=(time_win, 1),
                mode='nearest',
            )

        else:
            raise ValueError(f'un-supported temporal filter: {filter_type}!')
        del ts_data

        # prepare for writing: temporal referencing + reshape
        ts_data_filt -= ts_data_filt[self.refIndex, :]
        ts_data_filt = np.reshape(ts_data_filt, (self.numDate, self.length, self.width))

        # write
        self.write2hdf5(ts_data_filt, outFile=out_file, refFile=self.file)

        return


    def save2bl_list_file(self, out_file='bl_list.txt'):
        """Generate bl_list.txt file from timeseries h5 file."""
        self.open(print_msg=False)
        date6_list = [i[2:8] for i in self.dateList]
        pbase_list = self.pbase.tolist()
        print(f'write baseline list info to file: {out_file}')
        with open(out_file, 'w') as f:
            for d, pbase in zip(date6_list, pbase_list):
                f.write(f'{d}\t{pbase}\n')
        return out_file

################################ timeseries class end ##################################




################################# geometry class begin #################################
class geometry:
    """ Geometry object.

    File structure: https://mintpy.readthedocs.io/en/latest/api/data_structure/#geometry
    """

    def __init__(self, file=None):
        self.file = file
        self.name = 'geometry'

    def close(self, print_msg=True):
        try:
            self.f.close()
            if print_msg:
                print(f'close geometry file: {os.path.basename(self.file)}')
        except:
            pass

    def open(self, print_msg=True):
        if print_msg:
            print(f'open {self.name} file: {os.path.basename(self.file)}')
        self.get_metadata()
        self.get_size()
        self.numPixel = self.length * self.width
        self.geocoded = False
        if 'Y_FIRST' in self.metadata.keys():
            self.geocoded = True

        with h5py.File(self.file, 'r') as f:
            self.datasetNames = [i for i in f.keys() if isinstance(f[i], h5py.Dataset)]
            self.sliceList = list(self.datasetNames)
            if 'bperp' in f.keys():
                self.dateList = [i.decode('utf8') for i in f['date'][:]]
                self.numDate = len(self.dateList)
                # Update bperp datasetNames
                try:
                    self.sliceList.remove('bperp')
                except:
                    pass
                self.sliceList += ['bperp-'+d for d in self.dateList]
            else:
                self.dateList = None

    def get_size(self):
        with h5py.File(self.file, 'r') as f:
            dsName = [i for i in f.keys() if i in GEOMETRY_DSET_NAMES][0]
            dsShape = f[dsName].shape
            if len(dsShape) == 3:
                self.length, self.width = dsShape[1:3]
            else:
                self.length, self.width = dsShape
        return self.length, self.width

    def get_metadata(self):
        with h5py.File(self.file, 'r') as f:
            self.metadata = dict(f.attrs)
        for key, value in self.metadata.items():
            try:
                self.metadata[key] = value.decode('utf8')
            except:
                self.metadata[key] = value
        return self.metadata

    def read(self, datasetName=GEOMETRY_DSET_NAMES[0], box=None, print_msg=True):
        """Read 2D / 3D dataset with bounding box in space
        Parameters: datasetName : (list of) string, to point to specific 2D dataset, e.g.:
                        height
                        incidenceAngle
                        bperp
                        ...
                        bperp-20161020
                        bperp-20161026
                        bperp-...
                    box : tuple of 4 int, for (x0,y0,x1,y1)
                    print_msg : bool
        Returns: data : 2D or 3D array
        Example:
            obj = geometry('./inputs/geometryRadar.h5')
            obj.read(datasetName='height')
            obj.read(datasetName='incidenceAngle')
            obj.read(datasetName='bperp')
            obj.read(datasetName='bperp-20161020')
            obj.read(datasetName=['bperp-20161020',
                                  'bperp-20161026'])
        """
        self.open(print_msg=False)
        if box is None:
            box = (0, 0, self.width, self.length)

        if datasetName is None:
            datasetName = GEOMETRY_DSET_NAMES[0]
        elif isinstance(datasetName, str):
            datasetName = [datasetName]

        with h5py.File(self.file, 'r') as f:
            familyName = datasetName[0].split('-')[0]
            ds = f[familyName]
            if print_msg:
                print(f'reading {familyName:<15} data from file: {self.file} ...')

            if len(ds.shape) == 1:
                data = ds[:]
            elif len(ds.shape) == 2:
                data = ds[box[1]:box[3], box[0]:box[2]]
            else:
                # get dateFlag - mark in time/1st dimension
                dateFlag = np.zeros((ds.shape[0]), dtype=np.bool_)
                datasetName = [i.replace(familyName, '').replace('-', '') for i in datasetName]
                if any(not i for i in datasetName):
                    dateFlag[:] = True
                else:
                    for e in datasetName:
                        dateFlag[self.dateList.index(e)] = True

                # read
                data = ds[:,
                          box[1]:box[3],
                          box[0]:box[2]][dateFlag]

                if any(i == 1 for i in data.shape):
                    data = np.squeeze(data)
        return data
################################# geometry class end ###################################



################################# ifgramStack class begin ##############################
class ifgramStack:
    """ Interferograms Stack object.

    File structure: https://mintpy.readthedocs.io/en/latest/api/data_structure/#ifgramstack
    """

    def __init__(self, file=None):
        self.file = file
        self.name = 'ifgramStack'

    def close(self, print_msg=True):
        try:
            self.f.close()
            if print_msg:
                print(f'close {self.name} file: {os.path.basename(self.file)}')
        except:
            pass

    def open(self, print_msg=True):
        """
        Time format/rules:
            All datetime.datetime objects named with time
            All string in YYYYMMDD        named with date (following roipac)
        """
        if print_msg:
            print(f'open {self.name} file: {os.path.basename(self.file)}')
        self.get_metadata()
        self.get_size()
        self.read_datetimes()
        self.numPixel = self.length * self.width

        # time info
        self.date12List = [f'{i}_{j}' for i, j in zip(self.mDates, self.sDates)]
        self.tbaseIfgram = np.array([i.days + i.seconds / (24 * 60 * 60)
                                     for i in (self.sTimes - self.mTimes)],
                                    dtype=np.float32)

        with h5py.File(self.file, 'r') as f:
            self.dropIfgram = f['dropIfgram'][:]
            self.pbaseIfgram = f['bperp'][:]

            # get existed datasetNames in the order of IFGRAM_DSET_NAMES
            dsNames = [i for i in f.keys()
                       if (isinstance(f[i], h5py.Dataset)
                           and f[i].shape[-2:] == (self.length, self.width))]
            self.datasetNames = [i for i in IFGRAM_DSET_NAMES if i in dsNames]
            self.datasetNames += [i for i in dsNames if i not in IFGRAM_DSET_NAMES]

        # Get sliceList for self.read()
        self.sliceList = []
        for dsName in self.datasetNames:
            self.sliceList += [f'{dsName}-{i}' for i in self.date12List]

        # Time in timeseries domain
        self.dateList = self.get_date_list(dropIfgram=False)
        self.numDate = len(self.dateList)

        # Reference pixel
        try:
            self.refY = int(self.metadata['REF_Y'])
            self.refX = int(self.metadata['REF_X'])
        except:
            self.refY = None
            self.refX = None
        try:
            self.refLat = float(self.metadata['REF_LAT'])
            self.refLon = float(self.metadata['REF_LON'])
        except:
            self.refLat = None
            self.refLon = None

    def get_metadata(self):
        # read metadata from root level
        with h5py.File(self.file, 'r') as f:
            self.metadata = dict(f.attrs)
            dates = f['date'][:].flatten()

        # decode metadata
        for key, value in self.metadata.items():
            try:
                self.metadata[key] = value.decode('utf8')
            except:
                self.metadata[key] = value

        # START/END_DATE
        dateList = sorted(i.decode('utf8') for i in dates)
        self.metadata['START_DATE'] = dateList[0]
        self.metadata['END_DATE'] = dateList[-1]
        return self.metadata

    def get_size(self, dropIfgram=False, datasetName=None):
        with h5py.File(self.file, 'r') as f:
            # get default datasetName
            if datasetName is None:
                datasetName = [i for i in ['unwrapPhase', 'rangeOffset', 'azimuthOffset'] if i in f.keys()][0]

            # get 3D size
            self.numIfgram, self.length, self.width = f[datasetName].shape

            # update 1st dimension size
            if dropIfgram:
                self.numIfgram = np.sum(f['dropIfgram'][:])
        return self.numIfgram, self.length, self.width

    def read_datetimes(self):
        """Read date1/2 into array of datetime.datetime objects"""
        with h5py.File(self.file, 'r') as f:
            dates = f['date'][:]

        # grab the date string format
        self.dateFormat = ptime.get_date_str_format(dates[0, 0])

        # convert date from str to datetime.datetime objects
        self.mDates = np.array([i.decode('utf8') for i in dates[:, 0]])
        self.sDates = np.array([i.decode('utf8') for i in dates[:, 1]])
        self.mTimes = np.array([dt.datetime.strptime(i, self.dateFormat) for i in self.mDates])
        self.sTimes = np.array([dt.datetime.strptime(i, self.dateFormat) for i in self.sDates])

    def read(self, datasetName='unwrapPhase', box=None, print_msg=True, dropIfgram=False):
        """Read 3D dataset with bounding box in space
        Parameters: datasetName : string, to point to specific 2D dataset, e.g.:
                        unwrapPhase
                        coherence
                        connectComponent
                        ...
                        unwrapPhase-20161020_20161026
                        unwrapPhase-...
                        coherence-20161020_20161026
                        ...
                        ['unwrapPhase-20161020_20161026',
                         'unwrapPhase-20161020_20161101',
                         ...]
                    box : tuple of 4 int, for (x0,y0,x1,y1)
                    print_msg : bool
        Returns: data : 2D or 3D array
        Example:
            obj = ifgramStack('./inputs/ifgramStack.h5')
            obj.read(datasetName='unwrapPhase')
            obj.read(datasetName='coherence')
            obj.read(datasetName='unwrapPhase-20161020_20161026')
            obj.read(datasetName=['unwrapPhase-20161020_20161026',
                                  'unwrapPhase-20161020_20161101'])
        """
        self.get_size(dropIfgram=False)
        date12List = self.get_date12_list(dropIfgram=False)

        # convert input datasetName into list
        if datasetName is None:
            datasetName = ['unwrapPhase']
        elif isinstance(datasetName, str):
            datasetName = [datasetName]

        with h5py.File(self.file, 'r') as f:
            familyName = datasetName[0].split('-')[0]
            ds = f[familyName]
            if print_msg:
                print(f'reading {familyName} data from file: {self.file} ...')

            # get dateFlag - mark in time/1st dimension
            dateFlag = np.zeros((self.numIfgram), dtype=np.bool_)
            datasetName = [i.replace(familyName, '').replace('-', '') for i in datasetName]
            if any(not i for i in datasetName):
                if dropIfgram:
                    dateFlag = f['dropIfgram'][:]
                else:
                    dateFlag[:] = True
            else:
                for e in datasetName:
                    dateFlag[date12List.index(e)] = True

            # get index in space/2-3 dimension
            if box is None:
                box = (0, 0, self.width, self.length)

            # read
            data = ds[:,
                      box[1]:box[3],
                      box[0]:box[2]][dateFlag]

            if any(i == 1 for i in data.shape):
                data = np.squeeze(data)
        return data

    def spatial_average(self, datasetName='coherence', maskFile=None, box=None, useMedian=False,
                        reverseMask=False, threshold=None):
        """ Calculate the spatial average."""
        if datasetName is None:
            datasetName = 'coherence'

        if useMedian:
            print(f'calculating spatial median of {datasetName} in file {self.file} ...')
        else:
            print(f'calculating spatial mean of {datasetName} in file {self.file} ...')

        # read mask
        if maskFile and os.path.isfile(maskFile):
            print('read mask from file: '+maskFile)
            mask = singleDataset(maskFile).read(box=box)
        else:
            maskFile = None

        # calculation
        with h5py.File(self.file, 'r') as f:
            dset = f[datasetName]
            numIfgram = dset.shape[0]
            dmean = np.zeros((numIfgram), dtype=np.float32)

            prog_bar = ptime.progressBar(maxValue=numIfgram)
            for i in range(numIfgram):
                prog_bar.update(i+1, suffix=f'{i+1}/{numIfgram}')

                # read
                data = dset[i, box[1]:box[3], box[0]:box[2]]
                if maskFile:
                    data[mask == int(reverseMask)] = np.nan

                # ignore ZERO value for coherence
                if datasetName == 'coherence':
                    data[data == 0] = np.nan

                # calculate area ratio if threshold is specified
                # percentage of pixels with value above the threshold
                if threshold is not None:
                    data[data > threshold] = 1
                    data[data <= threshold] = 0

                if useMedian:
                    dmean[i] = np.nanmedian(data)
                else:
                    dmean[i] = np.nanmean(data)
            prog_bar.close()
        return dmean, self.date12List

    # Functions considering dropIfgram value
    def get_date12_list(self, dropIfgram=True):
        with h5py.File(self.file, 'r') as f:
            dates = f['date'][:]
            if dropIfgram:
                dates = dates[f['dropIfgram'][:], :]
        mDates = np.array([i.decode('utf8') for i in dates[:, 0]])
        sDates = np.array([i.decode('utf8') for i in dates[:, 1]])
        date12List = [f'{i}_{j}' for i, j in zip(mDates, sDates)]
        return date12List

    def get_drop_date12_list(self):
        with h5py.File(self.file, 'r') as f:
            dates = f['date'][:]
            dates = dates[~f['dropIfgram'][:], :]
        mDates = np.array([i.decode('utf8') for i in dates[:, 0]])
        sDates = np.array([i.decode('utf8') for i in dates[:, 1]])
        date12List = [f'{i}_{j}' for i, j in zip(mDates, sDates)]
        return date12List

    def get_date_list(self, dropIfgram=False):
        with h5py.File(self.file, 'r') as f:
            dates = f['date'][:]
            if dropIfgram:
                dates = dates[f['dropIfgram'][:], :]
        mDates = [i.decode('utf8') for i in dates[:, 0]]
        sDates = [i.decode('utf8') for i in dates[:, 1]]
        dateList = sorted(list(set(mDates + sDates)))
        return dateList

    def get_reference_phase(self, unwDatasetName='unwrapPhase', skip_reference=False, dropIfgram=False):
        """Get reference value
        Parameters: unwDatasetName : string, unwrapPhase, or unwrapPhase_unwCor
                    skip_reference : bool, skip reference value (for simulation only)
                    dropIfgram : bool, skip ifgrams marked as dropped or not
        Returns:    ref_phase : 1D np.array in size of (num_ifgram,) in float32
        """
        self.open(print_msg=False)
        if skip_reference:
            ref_phase = np.zeros(self.get_size(dropIfgram=dropIfgram)[0], np.float32)
            print('skip checking reference pixel info - This is for offset and testing ONLY.')
        elif 'REF_Y' not in self.metadata.keys():
            raise ValueError('No REF_X/Y found!\nrun reference_point.py to select reference pixel.')
        else:
            print(f'reference pixel in y/x: ({self.refY}, {self.refX}) from dataset: {unwDatasetName}')
            ref_phase = self.read(datasetName=unwDatasetName,
                                  box=(self.refX, self.refY, self.refX+1, self.refY+1),
                                  dropIfgram=dropIfgram,
                                  print_msg=False)
        return ref_phase

    def nonzero_mask(self, datasetName=None, print_msg=True, dropIfgram=True):
        """Return the common mask of pixels with non-zero value in dataset of all ifgrams.
           Ignoring dropped ifgrams
        """
        self.open(print_msg=False)
        with h5py.File(self.file, 'r') as f:
            if datasetName is None:
                datasetName = [i for i in ['connectComponent', 'unwrapPhase']
                               if i in f.keys()][0]
            print(f'calculate the common mask of pixels with non-zero {datasetName} value')

            dset = f[datasetName]
            mask = np.ones(dset.shape[1:3], dtype=np.bool_)
            dropIfgramFlag = np.ones(dset.shape[0], dtype=np.bool_)
            if dropIfgram:
                dropIfgramFlag = self.dropIfgram
            num2read = np.sum(dropIfgramFlag)
            idx2read = np.where(dropIfgramFlag)[0]

            # Loop to save memory usage
            prog_bar = ptime.progressBar(maxValue=num2read)
            for i in range(num2read):
                prog_bar.update(i+1, suffix=f'{i+1}/{num2read}')
                data = dset[idx2read[i], :, :]
                mask[data == 0.] = 0
                mask[np.isnan(data)] = 0
            prog_bar.close()
        return mask

    def temporal_average(self, datasetName='coherence', dropIfgram=True, max_memory=4):
        self.open(print_msg=False)
        if datasetName is None:
            datasetName = 'coherence'
        print(f'calculate the temporal average of {datasetName} in file {self.file} ...')

        # index of pairs to read
        ifgram_flag = np.ones(self.numIfgram, dtype=np.bool_)
        if dropIfgram:
            ifgram_flag = self.dropIfgram
            if np.all(ifgram_flag == 0.):
                raise Exception('ALL interferograms are marked as dropped, '
                                'can not calculate temporal average.')

        # temporal baseline for phase
        # with unit of years (float64 for very short tbase of UAVSAR data)
        if 'unwrapPhase' in datasetName:
            phase2range = -1 * float(self.metadata['WAVELENGTH']) / (4.0 * np.pi)
            tbase = np.array(self.tbaseIfgram, dtype=np.float64) / 365.25
            tbase = tbase[ifgram_flag]

        with h5py.File(self.file, 'r') as f:
            dset = f[datasetName]

            # reference value for phase
            ref_val = None
            if ('unwrapPhase' in datasetName
                   and self.refY is not None and 0 <= self.refY <= self.width
                   and self.refX is not None and 0 <= self.refX <= self.length):
                ref_val = dset[:, self.refY, self.refX][ifgram_flag]

            # get step size and number
            ds_size = np.sum(ifgram_flag, dtype=np.int64) * self.length * self.width * 4
            num_step = int(np.ceil(ds_size * 3 / (max_memory * 1024**3)))
            row_step = int(np.rint(self.length / num_step / 10) * 10)
            num_step = int(np.ceil(self.length / row_step))

            # calculate lines by lines
            dmean = np.zeros(dset.shape[1:3], dtype=np.float32)
            prog_bar = ptime.progressBar(maxValue=num_step)
            for i in range(num_step):
                r0 = i * row_step
                r1 = min(r0 + row_step, self.length)
                prog_bar.update(i+1, suffix=f'lines {r1}/{self.length}')

                # read
                data = dset[:, r0:r1, :][ifgram_flag]

                # referencing / normalizing for phase
                if 'unwrapPhase' in datasetName:
                    # spatial referencing
                    if ref_val is not None:
                        data -= np.tile(ref_val.reshape(-1, 1, 1), (1, data.shape[1], data.shape[2]))
                    # phase to phase velocity
                    for j in range(data.shape[0]):
                        data[j,:,:] *= (phase2range / tbase[j])

                # use nanmean to better handle NaN values
                dmean[r0:r1, :] = np.nanmean(data, axis=0)
            prog_bar.close()
        return dmean

    def get_max_connection_number(self):
        date12_list = self.get_date12_list()
        A = self.get_design_matrix4timeseries(date12_list, refDate=0)[0]
        num_conn = np.zeros(A.shape[0], dtype=np.int16)
        for i in range(A.shape[0]):
            Ai = A[i, :]
            num_conn[i] = np.where(Ai == 1)[0] - np.where(Ai == -1)[0]
        return np.max(num_conn)


    def split2boxes(self, max_memory=4, dim0_size=None, print_msg=True):
        """Split into chunks in rows to reduce memory usage.

        Parameters: max_memory - float, max memory to use in GB
                    dim0_size  - the 1st dimension size of all used datasets
                                 e.g., dim0_size = num_pair * 2 + num_date
                    print_msg  - bool
        Returns:    box_list   - list of tuple of 4 int
                    num_box    - int, number of boxes
        """
        self.open(print_msg=False)
        length = self.length
        width = self.width

        # dimension in time: phase/offset, weight, timeseries, etc.
        if not dim0_size:
            # for time series estimation
            dim0_size = self.numIfgram * 2 + self.numDate
        ds_size = dim0_size * length * width * 4

        num_box = int(np.ceil(ds_size * 1.5 / (max_memory * 1024**3)))
        y_step = int(np.ceil((length / num_box) / 10) * 10)
        num_box = int(np.ceil(length / y_step))
        if print_msg and num_box > 1:
            print('maximum memory size: %.1E GB' % max_memory)
            print('split %d lines into %d patches for processing' % (length, num_box))
            print('    with each patch up to %d lines' % y_step)

        # y_step / num_box --> box_list
        box_list = []
        for i in range(num_box):
            y0 = i * y_step
            y1 = min([length, y0 + y_step])
            box = (0, y0, width, y1)
            box_list.append(box)

        return box_list, num_box


    # Functions for closure phase bias
    def get_closure_phase_index(self, conn, dropIfgram=True):
        """Get the indices of interferograms that forms the given connection level closure loop.

        Parameters: conn       - int, connection level
                    dropIfgram - bool, exclude the dropped interferograms.
        Returns:    cp_idx     - 2D np.ndarray in int16 in size of (num_cp, conn + 1)
                                 Each row for the indices of interferograms for one closure loop.
                                 num_cp <= num_date - conn
        """
        date12_list = self.get_date12_list(dropIfgram=False)
        date_list = self.get_date_list(dropIfgram=dropIfgram)
        num_date = len(date_list)

        # get the closure index
        cp_idx = []
        for i in range(num_date - conn):
            # compose the connection-n pairs
            cp_date12_list = []
            for j in range(conn):
                cp_date12_list.append(f'{date_list[i+j]}_{date_list[i+j+1]}')
            cp_date12_list.append(f'{date_list[i]}_{date_list[i+conn]}')

            # add to cp_idx, ONLY IF all pairs exist for this closure loop
            if all(x in date12_list for x in cp_date12_list):
                cp_idx.append([date12_list.index(x) for x in cp_date12_list])

        # list(list) to 2D array
        cp_idx = np.array(cp_idx, dtype=np.int16)
        cp_idx = np.unique(cp_idx, axis=0)

        return cp_idx


    def get_sequential_closure_phase(self, box, conn, post_proc=None):
        """Computes wrapped sequential closure phases for a given connection level.

        Reference: Equation (21) in Zheng et al. (2022, TGRS)
        For conn = 5, seq_closure_phase = p12 + p23 + p34 + p45 + p56 - p16.

        Parameters: box       - tuple of 4 int, bounding box in (x0, y0, x1, y1)
                    conn      - int, connection level of the closure phase
                    post_proc - str, post processing of the closure phase:
                                None - 3D array in float32, seq closure phase
                                sum  - 2D array in complex64, sum  in time of the complex seq closure phase
                                mean - 2D array in complex64, mean in time of the complex seq closure phase
        Returns:    cp_w      - 3D np.ndarray in float32 in size of (num_cp, box_len, box_wid)
                                wrapped sequential  closure phase for the given connection level.
                    sum_cp    - None or 2D np.ndarray in complex64 in size of (box_len, box_width)
                                wrapped average seq closure phase for the given connection level,
                                controlled by post_proc.
                    num_cp    - int, number of  seq closure phase for the given connection level.
        """
        # basic info
        num_date = len(self.get_date_list(dropIfgram=True))
        box_wid = box[2] - box[0]
        box_len = box[3] - box[1]

        ## get the closure index
        cp_idx = self.get_closure_phase_index(conn=conn, dropIfgram=True)
        num_cp = cp_idx.shape[0]
        print(f'number of closure measurements expected: {num_date - conn}')
        print(f'number of closure measurements found   : {num_cp}')

        if not post_proc:
            if num_cp < num_date - conn:
                msg = f'num_cp ({num_cp}) < num_date - conn ({num_date - conn})'
                msg += ' --> some interferograms are missing!'
                raise Exception(msg)
        else:
            if num_cp < 1:
                raise Exception(f"No triplets found at connection level: {conn}!")

        ## read data
        # Note by Yujie Zheng on Dec 2, 2022: I removed the lines for spatial referencing,
        # because no moisture change should have a naturally zero closure phase. On the other
        # hand, choosing a reference point in the wrong place can lead to misinterpretations.
        # Tests in both Central Valley and Bristol Dry lakes confirms this.
        ds_name = 'wrapPhase' if 'wrapPhase' in self.datasetNames else 'unwrapPhase'
        print(f'reading {ds_name} to compute closure phases')
        phase = self.read(datasetName=ds_name, box=box, print_msg=False)

        # apply spatial referencing (for ARIA only)
        # to avoid the abnormal result as shown in https://github.com/insarlab/MintPy/pull/1063
        # This may be caused by the phase stitching during product preparation via ARIA-tools,
        # which could have broken the temporal consistency of the native unwrapped phase.
        processor = self.metadata.get('mintpy.load.processor', 'isce')
        if ds_name == 'unwrapPhase' and processor in ['aria']:
            print(f'apply spatial referencing to {processor} products')
            ref_phase = self.get_reference_phase(dropIfgram=False)
            for i in range(phase.shape[0]):
                mask = phase[i] != 0.
                phase[i][mask] -= ref_phase[i]

        ## calculate the 3D complex seq closure phase
        cp_w = np.zeros((num_cp, box_len, box_wid), dtype=np.complex64)
        for i in range(num_cp):

            # calculate (un)wrapped closure phase
            idx_plus, idx_minor = cp_idx[i, :-1], cp_idx[i, -1]
            cp0_w = np.sum(phase[idx_plus], axis=0) - phase[idx_minor]

            # re-wrapping to ensure the wrapped closure phase
            cp_w[i] = np.exp(1j * cp0_w)

        ## post-processing
        if not post_proc:
            sum_cp = None

        elif post_proc == 'sum':
            sum_cp = np.sum(cp_w, axis=0)

        elif post_proc == 'mean':
            sum_cp = np.mean(cp_w, axis=0)

        else:
            raise ValueError(f'un-recognized post_proc={post_proc}! Available choices: sum, mean.')

        return np.angle(cp_w), sum_cp, num_cp


    # Functions for unwrapping error correction
    @staticmethod
    def get_design_matrix4triplet(date12_list):
        """Generate the design matrix of ifgram triangle for unwrap error correction using phase closure

        Parameters: date12_list : list of string in YYYYMMDD_YYYYMMDD format
        Returns:    C : 2D np.array in size of (num_tri, num_ifgram) consisting 0, 1, -1
                        for 3 SAR acquisition in t1, t2 and t3 in time order,
                        ifg1 for (t1, t2) with 1
                        ifg2 for (t1, t3) with -1
                        ifg3 for (t2, t3) with 1
        Examples:   obj = ifgramStack('./inputs/ifgramStack.h5')
                    date12_list = obj.get_date12_list(dropIfgram=True)
                    C = ifgramStack.get_design_matrix4triplet(date12_list)
        """
        # Date info
        date12_list = list(date12_list)
        # Use tuples of the dates, which are hashable for the lookup dict
        date12_tuples = [tuple(d.split('_')) for d in date12_list]

        # Create an inverse map from tuple(date1, date1) -> index in ifg list
        ifg_to_idx = {ifg: idx for idx, ifg in enumerate(date12_tuples)}

        # Get the unique SAR dates present in the interferogram list
        date_list = sorted(set(itertools.chain.from_iterable(date12_tuples)))

        # Start with all possible triplets, narrow down based on ifgs present
        closure_list = itertools.combinations(date_list, 3)

        M = len(date12_tuples)  # Number of igrams, number of rows
        C_list = []
        for date1, date2, date3 in closure_list:
            ifg12 = (date1, date2)
            ifg23 = (date2, date3)
            ifg13 = (date1, date3)
            # Check if any ifg is not available in the current triple. Skip if so
            try:
                idx12 = ifg_to_idx[ifg12]
                idx23 = ifg_to_idx[ifg23]
                idx13 = ifg_to_idx[ifg13]
            except KeyError:
                continue

            # Add the +/-1 row of the matrix to our list
            row = np.zeros(M, dtype=np.int8)
            row[idx12] = 1
            row[idx23] = 1
            row[idx13] = -1
            C_list.append(row)

        if len(C_list) == 0:
            print(f'\nWARNING: No triangles found from input date12_list:\n{date12_list}!\n')
            return None

        return np.stack(C_list).astype(np.float32)


    # Functions for network inversion / time series estimation
    @staticmethod
    def get_design_matrix4timeseries(date12_list, refDate=None):
        """Return design matrix of the input ifgramStack for timeseries estimation

        Parameters: date12_list - list of string in YYYYMMDD_YYYYMMDD format
                    refDate     - str, date in YYYYMMDD format
                                  set to None for the 1st date
                                  set to 'no' to disable reference date
        Returns:    A - 2D array of float32 in size of (num_ifgram, num_date-1)
                    B - 2D array of float32 in size of (num_ifgram, num_date-1)
        Examples:   obj = ifgramStack('./inputs/ifgramStack.h5')
                    A, B = obj.get_design_matrix4timeseries(obj.get_date12_list(dropIfgram=True))
                    A = ifgramStack.get_design_matrix4timeseries(date12_list, refDate='20101022')[0]
                    A = ifgramStack.get_design_matrix4timeseries(date12_list, refDate=0)[0] #do not omit the 1st column
        """
        # Date info
        date12_list = list(date12_list)
        date1s = [i.split('_')[0] for i in date12_list]
        date2s = [i.split('_')[1] for i in date12_list]
        date_list = sorted(list(set(date1s + date2s)))
        num_ifgram = len(date12_list)
        num_date = len(date_list)

        # tbase in the unit of years
        date_format = ptime.get_date_str_format(date_list[0])
        dates = np.array([dt.datetime.strptime(i, date_format) for i in date_list])
        tbase = [i.days + i.seconds / (24 * 60 * 60) for i in (dates - dates[0])]
        tbase = np.array(tbase, dtype=np.float32) / 365.25

        # calculate design matrix
        # A for minimizing the residual of phase
        # B for minimizing the residual of phase velocity
        A = np.zeros((num_ifgram, num_date), np.float32)
        B = np.zeros((num_ifgram, num_date), np.float32)
        for i in range(num_ifgram):
            ind1, ind2 = (date_list.index(d) for d in date12_list[i].split('_'))
            A[i, ind1] = -1
            A[i, ind2] = 1
            # support date12_list with the first date NOT being the earlier date
            if ind1 < ind2:
                B[i, ind1:ind2] = tbase[ind1 + 1:ind2 + 1] - tbase[ind1:ind2]
            else:
                B[i, ind2:ind1] = tbase[ind2:ind1] - tbase[ind2 + 1:ind1 + 1]

        # Remove reference date as it can not be resolved
        if refDate != 'no':
            # default refDate
            if refDate is None:
                # for single   reference network, use the same reference date
                # for multiple reference network, use the first date
                if len(set(date1s)) == 1:
                    refDate = date1s[0]
                else:
                    refDate = date_list[0]

            # apply refDate
            if refDate:
                ind_r = date_list.index(refDate)
                A = np.hstack((A[:, 0:ind_r], A[:, (ind_r+1):]))
                B = B[:, :-1]

        return A, B


    def get_perp_baseline_timeseries(self, dropIfgram=True):
        """Get spatial perpendicular baseline in timeseries from ifgramStack, ignoring dropped ifgrams"""
        # read pbase of interferograms
        with h5py.File(self.file, 'r') as f:
            pbaseIfgram = f['bperp'][:]
            if dropIfgram:
                pbaseIfgram = pbaseIfgram[f['dropIfgram'][:]]

        # estimate pbase of time-series
        date12List = self.get_date12_list(dropIfgram=dropIfgram)
        A = self.get_design_matrix4timeseries(date12List)[0]
        pbaseTimeseries = np.zeros(A.shape[1]+1, dtype=np.float32)
        pbaseTimeseries[1:] = np.linalg.lstsq(A, pbaseIfgram, rcond=None)[0]
        return pbaseTimeseries

    def update_drop_ifgram(self, date12List2Drop):
        """Update dropIfgram dataset based on input date12List2Drop"""
        if date12List2Drop is None:
            return
        date12ListAll = self.get_date12_list(dropIfgram=False)

        # check date12 list to drop: old vs new
        date12ListKeptOld = self.get_date12_list(dropIfgram=True)
        date12List2DropOld = sorted(list(set(date12ListAll) - set(date12ListKeptOld)))
        if date12List2DropOld == date12List2Drop:
            print('The same date12List2Drop / dropIfgram is already marked in the file, skip updating dropIfgram.')
            return

        with h5py.File(self.file, 'r+') as f:
            print(f'open file {self.file} with r+ mode')
            print('update HDF5 dataset "/dropIfgram".')
            f['dropIfgram'][:] = np.array([i not in date12List2Drop for i in date12ListAll], dtype=np.bool_)

            # update MODIFICATION_TIME for all datasets in IFGRAM_DSET_NAMES
            for dsName in IFGRAM_DSET_NAMES:
                if dsName in f.keys():
                    print(f'update MODIFICATION_TIME in HDF5 dataset "/{dsName}"')
                    f[dsName].attrs['MODIFICATION_TIME'] = str(time.time())
                    time.sleep(1)   #to distinguish the modification time of input files

################################# ifgramStack class end ################################



########################################################################################
class singleDataset:
    def __init__(self, file=None):
        self.file = file

    def read(self, box=None):
        with h5py.File(self.file, 'r') as f:
            dsName = list(f.keys())[0]
            data = f[dsName][:]

        if box is not None:
            data = data[box[1]:box[3],
                        box[0]:box[2]]
        return data


################################# HDF-EOS5 class begin #################################
class HDFEOS:
    """
    Time-series object in HDF-EOS5 format for Univ of Miami's InSAR Time-series Web Viewer
        Link: http://insarmaps.miami.edu
    It contains a "timeseries" group and three datasets: date, bperp and timeseries.

    File structure: https://mintpy.readthedocs.io/en/latest/hdfeos5/#file_structure
    """

    def __init__(self, file=None):
        self.file = file
        self.name = 'HDFEOS'
        self.datasetGroupNameDict = {'displacement'       : 'observation',
                                     'raw'                : 'observation',
                                     'troposphericDelay'  : 'observation',
                                     'topographicResidual': 'observation',
                                     'ramp'               : 'observation',
                                     'date'               : 'observation',
                                     'temporalCoherence'  : 'quality',
                                     'mask'               : 'quality',
                                     'coherence'          : 'quality',
                                     'variance'           : 'quality',
                                     'uncertainty'        : 'quality',
                                     'height'             : 'geometry',
                                     'incidenceAngle'     : 'geometry',
                                     'slantRangeDistance' : 'geometry',
                                     'azimuthAngle'       : 'geometry',
                                     'shadowMask'         : 'geometry',
                                     'waterMask'          : 'geometry',
                                     'bperp'              : 'geometry'
                                    }

    def close(self, print_msg=True):
        try:
            self.f.close()
            if print_msg:
                print(f'close timeseries file: {os.path.basename(self.file)}')
        except:
            pass

    def open(self, print_msg=True):
        if print_msg:
            print(f'open {self.name} file: {os.path.basename(self.file)}')
        self.get_metadata()
        self.length = int(self.metadata['LENGTH'])
        self.width = int(self.metadata['WIDTH'])

        self.sliceList = []
        with h5py.File(self.file, 'r') as f:
            gname = 'HDFEOS/GRIDS/timeseries/observation'
            g = f[gname]
            self.dateList = [i.decode('utf8') for i in g['date'][:]]
            self.pbase = g['bperp'][:]
            self.numDate = len(self.dateList)

            # get slice list
            self.sliceList += [f'{gname}/displacement-{i}' for i in self.dateList]
            for gname in ['HDFEOS/GRIDS/timeseries/quality',
                          'HDFEOS/GRIDS/timeseries/geometry']:
                g = f[gname]
                for key in g.keys():
                    if isinstance(g[key], h5py.Dataset) and len(g[key].shape) == 2:
                        self.sliceList.append(f'{gname}/{key}')

    def get_metadata(self):
        with h5py.File(self.file, 'r') as f:
            self.metadata = dict(f.attrs)
            dates = f['HDFEOS/GRIDS/timeseries/observation/date'][:]
        for key, value in self.metadata.items():
            try:
                self.metadata[key] = value.decode('utf8')
            except:
                self.metadata[key] = value
        self.metadata['FILE_TYPE'] = self.name

        # ref_date/index
        dateList = [i.decode('utf8') for i in dates]
        if 'REF_DATE' not in self.metadata.keys():
            self.metadata['REF_DATE'] = dateList[0]
        self.refIndex = dateList.index(self.metadata['REF_DATE'])
        return self.metadata

    def get_date_list(self):
        with h5py.File(self.file, 'r') as f:
            g = f['HDFEOS/GRIDS/timeseries/observation']
            self.dateList = [i.decode('utf8') for i in g['date'][:]]
        return self.dateList

    def read(self, datasetName=None, box=None, print_msg=True):
        """Read dataset from HDF-EOS5 file
        Parameters: self : HDFEOS object
                    datasetName : (list of) str
                    box : tuple of 4 int, for (x0, y0, x1, y1)
                    print_msg : bool
        Returns:    data: 2D or 3D array
        Example:    obj = HDFEOS('S1_IW1_128_0593_0597_20141213_20171221.he5')
                    obj.read(datasetName='displacement')
                    obj.read(datasetName='displacement-20150915')
                    obj.read(datasetName=['displacement-20150915',
                                          'displacement-20150921'])
                    obj.read(datasetName='incidenceAngle')
        """
        self.open(print_msg=False)
        if box is None:
            box = [0, 0, self.width, self.length]

        if datasetName is None:
            datasetName = [TIMESERIES_DSET_NAMES[-1]]
        elif isinstance(datasetName, str):
            datasetName = [datasetName]

        with h5py.File(self.file, 'r') as f:
            familyName = datasetName[0].split('-')[0]
            groupName = self.datasetGroupNameDict[familyName]
            ds = f[f'HDFEOS/GRIDS/timeseries/{groupName}/{familyName}']
            if print_msg:
                print(f'reading {familyName} data from file: {self.file} ...')

            if len(ds.shape) == 1:
                data = ds[:]
            elif len(ds.shape) == 2:
                data = ds[box[1]:box[3], box[0]:box[2]]
            else:
                # get dateFlag - mark in time/1st dimension
                dateFlag = np.zeros((ds.shape[0]), dtype=np.bool_)
                datasetName = [i.replace(familyName, '').replace('-', '') for i in datasetName]
                if any(not i for i in datasetName):
                    dateFlag[:] = True
                else:
                    for e in datasetName:
                        dateFlag[self.dateList.index(e)] = True

                # read
                data = ds[:,
                          box[1]:box[3],
                          box[0]:box[2]][dateFlag]

                # squeeze/shrink dimension whenever it is possible
                if any(i == 1 for i in data.shape):
                    data = np.squeeze(data)
        return data
################################# HDF-EOS5 class end ###################################