resample.py

############################################################
# Program is part of MintPy                                #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
# Author: Zhang Yunjun, 2018                               #
############################################################
# Recommend import:
#     from mintpy.objects.resample import resample


import os
import warnings
import h5py
import numpy as np
from scipy import ndimage
from scipy.interpolate import RegularGridInterpolator as RGI
with warnings.catch_warnings():
    warnings.simplefilter('ignore', UserWarning)
    import pyresample as pr
from mintpy.objects.cluster import split_box2sub_boxes
from mintpy.utils import readfile, ptime, utils0 as ut

EARTH_RADIUS = 6378122.65   # m


class resample:
    """
    Geometry Definition objects for geocoding using:
    1) pyresample (http://pyresample.readthedocs.org)
    2) scipy.interpolate.RegularGridInterpolator:
       (https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html)

    Example:
        from mintpy.objects.resample import resample
        from mintpy.utils import readfile, attribute as attr

        ##### opt 1 - entire matrix (by not changing max_memory=0)
        res_obj = resample(lut_file='./inputs/geometryRadar.h5', src_file='velocity.h5')
        # OR use ISCE-2 lookup table files instead of MintPy geometry file in HDF5 format
        res_obj = resample(lut_file='../merged/geom_reference/lat.rdr', src_file='velocity.h5',
                           lat_file='../merged/geom_reference/lat.rdr',
                           lon_file='../merged/geom_reference/lon.rdr')

        res_obj.open()
        res_obj.prepare()

        # read data / attribute
        rdr_data, atr = readfile.read(src_file)
        # resample data
        box = res_obj.src_box_list[0]
        geo_data = res_obj.run_resample(src_data=rdr_data[box[1]:box[3], box[0]:box[2]])
        # update attribute
        atr = attr.update_attribute4radar2geo(atr, res_obj=res_obj)

        ##### opt 2 - block-by-block IO (by setting max_memory=4)
        res_obj = resample(lut_file='./inputs/geometryRadar.h5', src_file='timeseries.h5', max_memory=4)
        res_obj.open()
        res_obj.prepare()

        # prepare output: metadata and initial file
        atr = readfile.read_attribute(src_file)
        atr = attr.update_attribute4radar2geo(atr, res_obj=res_obj)
        writefile.layout_hdf5(outfile, metadata=atr, ref_file=src_file)

        # block-by-block IO
        for i in range(res_obj.num_bix):
            src_box = res_obj.src_box_list[i]
            dest_box = res_obj.dest_box_list[i]
            # read
            data = readfile.read(src_file, box=src_box)[0]
            # resample
            data = res_obj.run_resample(src_data=rdr_data, box_ind=i)
            # write
            block = [0, data.shape[0],
                     dest_box[1], dest_box[3],
                     dest_box[0], dest_box[2]]
            writefile.write_hdf5_block(outfile,
                                       data=data,
                                       datasetName='timeseries',
                                       block=block)
    """

    def __init__(self, lut_file, src_file=None, SNWE=None, lalo_step=None, interp_method='nearest', fill_value=np.nan,
                 nprocs=1, max_memory=0, software='pyresample', print_msg=True, lat_file=None, lon_file=None):
        """
        Parameters: lut_file      - str, path of lookup table file, containing datasets:
                                    latitude / longitude      for lut_file in radar-coord
                                    azimuthCoord / rangeCoord for lut_file in geo-coord
                    src_file      - str, path of data file to be resampled.
                    SNWE          - tuple of 4 float, indicating south/north/west/east
                                    coordinates at pixel outer boundary (consistent with Y/X_FIRST; not pixel center)
                    lalo_step     - list of 2 float, output step in lat/lon direction in degree / meter
                                    input lalo_step is used ONLY IF 1) radar2geo = True AND 2) lut_file is in radar-coord
                    interp_method - str, interpolation / resampling method, nearest / linear
                    fill_value    - number, fill value for extrapolation pixels
                    nprocs        - int, number of processors used in parallel
                    max_memory    - float, maximum memory to use
                                    set to 0 or negative value to disable block-by-block IO (default)
                    software      - str, interpolation software, pyresample / scipy
                    lat/lon_file  - str, path of the ISCE-2 lat/lon.rdr file
                                    To geocode file with ISCE-2 lookup table files directly,
                                    without using/loading geometry files in HDF5/MintPy format.
        """
        # input variables
        self.lut_file = lut_file
        # use isce lat/lon.rdr file, as an alternative to lut_file with mintpy geometry HDF5 file
        self.lat_file = lat_file
        self.lon_file = lon_file

        self.src_file = src_file
        self.SNWE = SNWE
        self.lalo_step = lalo_step
        self.interp_method = interp_method
        self.fill_value = fill_value

        self.nprocs = nprocs
        self.max_memory = max_memory
        self.software = software
        self.print_msg = print_msg

        # initial variables
        self.num_box = 1
        self.radius = None


    def open(self):
        """Read metadata

        Note: separate open() from prepare() to handle numpy matrix directly w/o src_file
        by assigning src_meta after open() and before prepare()
        """
        # read metadata: lookup table and source data
        if self.lut_file:
            self.lut_meta = readfile.read_attribute(self.lut_file)
        else:
            self.lut_meta = None

        if self.src_file:
            self.src_meta = readfile.read_attribute(self.src_file)
        else:
            self.src_meta = None

        # get num_box
        if self.software == 'pyresample':
            self.num_box = self.get_num_box(self.src_file, self.max_memory)


    def prepare(self):
        """Prepare aux data before interpolation operation"""

        # check metadata: lookup table and source data
        if self.lut_meta is None or self.src_meta is None:
            raise ValueError('lookup table or source data metadata is None!')

        # prepare geometry for resampling
        print('resampling software: {}'.format(self.software))
        if self.software == 'scipy':
            if 'Y_FIRST' in self.lut_meta.keys():
                # gamma / roipac
                self.prepare_regular_grid_interpolator()
            else:
                raise ValueError('resampling using scipy with lookup table in radar-coord (ISCE / DORIS) is NOT supported!')

        elif self.software == 'pyresample':
            if 'Y_FIRST' in self.lut_meta.keys():
                # gamma / roipac
                self.prepare_geometry_definition_geo()

            else:
                # isce / doris
                self.prepare_geometry_definition_radar()

            # aux info
            self.radius = self.get_radius_of_influence(self.lalo_step, self.src_meta)


    def run_resample(self, src_data, box_ind=0, print_msg=True):
        """Run interpolation operation for input 2D/3D data
        Parameters: src_data   - 2D/3D np.array, source data to be resampled
                    box_ind    - int, index of the current box of interest
                                 for multiple boxes with pyresample only
                    print_msg  - bool
        Returns:    dest_data  - 2D/3D np.array, resampled data
        """
        # adjust fill_value for each source data / block
        fill_value = self.fill_value
        float_types = [np.single, np.double, np.longdouble, np.csingle, np.cdouble, np.clongdouble]
        if src_data.dtype == np.bool_:
            fill_value = False
            if print_msg:
                print('input source data is bool type, restrict fill_value to False.')

        elif src_data.dtype not in float_types and np.isnan(fill_value):
            fill_value = 0
            if print_msg:
                print('input source data is NOT float, change fill_value from NaN to 0.')

        ## pyresample
        if self.software == 'pyresample':
            # move 1st/time dimension to the last
            # so that rows/cols axis are the frist, as required by pyresample
            if len(src_data.shape) == 3:
                src_data = np.moveaxis(src_data, 0, -1)

            # resample source data into target data
            dest_data = self.run_pyresample(src_data=src_data,
                                            src_def=self.src_def_list[box_ind],
                                            dest_def=self.dest_def_list[box_ind],
                                            radius=self.radius,
                                            interp_method=self.interp_method,
                                            fill_value=fill_value,
                                            nprocs=self.nprocs,
                                            print_msg=self.print_msg)

            # move 1st/time dimension back
            if len(dest_data.shape) == 3:
                dest_data = np.moveaxis(dest_data, -1, 0)

        ## scipy
        else:
            if print_msg:
                print('{} resampling using scipy.interpolate.RegularGridInterpolator ...'.format(self.interp_method))
            if len(src_data.shape) == 3:
                dest_data = np.empty((src_data.shape[0], self.length, self.width), src_data.dtype)
                prog_bar = ptime.progressBar(maxValue=src_data.shape[0], print_msg=print_msg)
                for i in range(src_data.shape[0]):
                    dest_data[i, :, :] = self.run_regular_grid_interpolator(src_data=src_data[i, :, :],
                                                                            interp_method=self.interp_method,
                                                                            fill_value=fill_value,
                                                                            print_msg=True)
                    prog_bar.update(i+1)
                prog_bar.close()
            else:
                dest_data = self.run_regular_grid_interpolator(src_data=src_data,
                                                               interp_method=self.interp_method,
                                                               fill_value=fill_value,
                                                               print_msg=True)
        return dest_data


    @staticmethod
    def get_num_box(src_file=None, max_memory=0, scale_fac=3.0):
        """Get the number of boxes to split to not exceed the max memory
        Parameters: src_file   - str, path of source data file (largest one if multiple)
                    max_memory - float, memory size in GB
                    scale_fac  - float, scale factor from data size to memory used
                                 empirically estimated.
        Returns:    num_box    - int, number of boxes to be splitted
        """
        num_box = 1

        # auto split into list of boxes if max_memory > 0
        if src_file and os.path.isfile(src_file) and max_memory > 0:
            # get max dataset shape
            fext = os.path.splitext(src_file)[1]
            if fext in ['.h5', '.he5']:
                with h5py.File(src_file, 'r') as f:
                    ds_shapes = [f[i].shape for i in f.keys()
                                 if isinstance(f[i], h5py.Dataset)]
                    max_ds_size = max([np.prod(i) for i in ds_shapes])
            else:
                atr = readfile.read_attribute(src_file)
                max_ds_size = int(atr['LENGTH']) * int(atr['WIDTH'])

            # calc num_box
            num_box = int(np.ceil((max_ds_size * 4 * 4) / (max_memory * 1024**3)))

        return num_box


    @staticmethod
    def get_radius_of_influence(lalo_step, src_meta, ratio=3):
        """Get radius of influence based on the lookup table resolution in lat/lon direction"""
        if 'Y_FIRST' in src_meta.keys():
            # geo2radar
            radius = 100e3
        else:
            # radar2geo
            step_deg = max(np.abs(lalo_step))
            step_m = step_deg * np.pi / 180.0 * EARTH_RADIUS
            radius = step_m * ratio
        return radius


    ##------------------------------ resample using pyresample -------------------------------------##

    def prepare_geometry_definition_radar(self):
        """Get src_def and dest_def for lookup table in radar-coord (from ISCE, DORIS)"""

        def mark_lat_lon_anomoly(lat, lon):
            """mask pixels with abnormal values (0, etc.)
            This is found on sentinelStack multiple swath lookup table file.
            """
            # ignore pixels with zero value
            zero_mask = np.multiply(lat != 0., lon != 0.)

            # ignore anomaly non-zero values
            # by get the most common data range (d_min, d_max) based on histogram
            mask = np.array(zero_mask, np.bool_)
            for data in [lat, lon]:
                bin_value, bin_edge = np.histogram(data[mask], bins=10)
                # if there is anomaly, histogram won't be evenly distributed
                while np.max(bin_value) > np.sum(zero_mask) * 0.3:
                    # find the continous bins where the largest bin is --> normal data range
                    bin_value_thres = ut.median_abs_deviation_threshold(bin_value, cutoff=3)
                    bin_label = ndimage.label(bin_value > bin_value_thres)[0]
                    idx = np.where(bin_label == bin_label[np.argmax(bin_value)])[0]
                    # convert to min/max data value
                    bin_step = bin_edge[1] - bin_edge[0]
                    d_min = bin_edge[idx[0]] - bin_step / 2.
                    d_max = bin_edge[idx[-1]+1] + bin_step / 2.
                    mask *= np.multiply(data >= d_min, data <= d_max)
                    bin_value, bin_edge = np.histogram(data[mask], bins=10)

            # set invalid pixels to fixed values
            lat[mask == 0] = 90.
            lon[mask == 0] = 0.
            return lat, lon, mask


        # read lookup table: lat/lon at pixel center
        # src  for radar2geo
        # dest for geo2radar
        print('read latitude / longitude from lookup table file: {}'.format(self.lut_file))
        lat_file = self.lat_file if self.lat_file else self.lut_file
        lon_file = self.lon_file if self.lon_file else self.lut_file
        lut_lat = readfile.read(lat_file, datasetName='latitude')[0].astype(np.float32)
        lut_lon = readfile.read(lon_file, datasetName='longitude')[0].astype(np.float32)
        lut_lat, lut_lon, mask = mark_lat_lon_anomoly(lut_lat, lut_lon)

        # radar2geo (with block-by-block support)
        if 'Y_FIRST' not in self.src_meta.keys():

            # src_lat/lon0/1
            src_lat0 = np.nanmax(lut_lat[mask])
            src_lat1 = np.nanmin(lut_lat[mask])
            src_lon0 = np.nanmin(lut_lon[mask])
            src_lon1 = np.nanmax(lut_lon[mask])

            # parameter 1 - lalo_step (output grid)
            if self.lalo_step is None:
                try:
                    # ensure the same pixel area before / after geocoding
                    merged_meta = {**self.lut_meta, **self.src_meta}
                    lat_c = (src_lat0 + src_lat1) / 2.
                    lat_step, lon_step = ut.auto_lat_lon_step_size(merged_meta, lat_c)

                except KeyError:
                    # ensure the same matrix shape before / after geocoding
                    # if not enough metadata found for the above
                    lat_step = (src_lat1 - src_lat0) / (lut_lat.shape[0] - 1)
                    lon_step = (src_lon1 - src_lon0) / (lut_lat.shape[1] - 1)
                self.lalo_step = (abs(lat_step) * -1., abs(lon_step))

            else:
                # ensure lat/lon step sign
                self.lalo_step = (abs(self.lalo_step[0]) * -1.,
                                  abs(self.lalo_step[1]) * 1.)
            print('output pixel size in (lat, lon) in degree: {}'.format(self.lalo_step))

            # parameter 2 / 3 - SNWE (at pixel outer boundary; output grid) / length & width
            if self.SNWE is None:
                self.SNWE = (src_lat1 + self.lalo_step[0] / 2.0,
                             src_lat0 - self.lalo_step[0] / 2.0,
                             src_lon0 - self.lalo_step[1] / 2.0,
                             src_lon1 + self.lalo_step[1] / 2.0)
            self.length = int(np.rint((self.SNWE[0] - self.SNWE[1]) / self.lalo_step[0]))
            self.width  = int(np.rint((self.SNWE[3] - self.SNWE[2]) / self.lalo_step[1]))
            # adjust SNWE ending coordinate (S, E) for precise alignment
            self.SNWE = (self.SNWE[1] + self.lalo_step[0] * self.length,
                         self.SNWE[1],
                         self.SNWE[2],
                         self.SNWE[2] + self.lalo_step[1] * self.width)
            print('output area extent in (S, N, W, E) in degree: {}'.format(self.SNWE))
            print('output file row / column number: ({}, {})'.format(self.length, self.width))

            # parameter 4 - list of boxes & geometry definitions
            self.src_box_list = []
            self.src_def_list = []
            self.dest_box_list = []
            self.dest_def_list = []

            # split dest_box (in grid)
            self.dest_box_list = split_box2sub_boxes(box=(0, 0, self.width, self.length),
                                                     num_split=self.num_box,
                                                     dimension='y',
                                                     print_msg=True)

            # dest_box --> src_box / src_def / dest_def
            for i, dest_box in enumerate(self.dest_box_list):
                if self.num_box > 1:
                    print('preparing geometry for dest_box {}/{}: {}'.format(i+1, self.num_box, dest_box))

                # dest_lat/lon at pixel center
                lat_num = dest_box[3] - dest_box[1]
                lon_num = dest_box[2] - dest_box[0]
                lat0 = self.SNWE[1] + self.lalo_step[0] * (dest_box[1] + 0.5)
                lat1 = self.SNWE[1] + self.lalo_step[0] * (dest_box[3] - 0.5)
                lon0 = self.SNWE[2] + self.lalo_step[1] * (dest_box[0] + 0.5)
                lon1 = self.SNWE[2] + self.lalo_step[1] * (dest_box[2] - 0.5)
                dest_lat, dest_lon = np.mgrid[lat0:lat1:lat_num*1j,
                                              lon0:lon1:lon_num*1j]

                # src_box
                src_area = (src_lat1 - src_lat0) * (src_lon1 - src_lon0)
                dest_area = (lat1 - lat0) * (lon1 - lon0)
                if dest_area < src_area * 0.5:
                    # reduction of swath data
                    # https://pyresample.readthedocs.io/en/latest/data_reduce.html
                    # get src_box (in swath) from lat/lon (from dest_box in grid)
                    print('searching relevant box covering the current SNWE')
                    flag = pr.data_reduce.get_valid_index_from_lonlat_grid(dest_lon,
                                                                           dest_lat,
                                                                           lut_lon,
                                                                           lut_lat,
                                                                           radius_of_influence=3000)
                    idx_row, idx_col = np.where(flag)
                    src_box = (np.min(idx_col), np.min(idx_row),
                               np.max(idx_col), np.max(idx_row))
                else:
                    src_box = (0, 0, lut_lat.shape[1], lut_lat.shape[0])

                # geometry definition
                src_def = pr.geometry.SwathDefinition(lons=lut_lon[src_box[1]:src_box[3],
                                                                   src_box[0]:src_box[2]],
                                                      lats=lut_lat[src_box[1]:src_box[3],
                                                                   src_box[0]:src_box[2]])
                dest_def = pr.geometry.GridDefinition(lons=dest_lon, lats=dest_lat)

                self.src_box_list.append(src_box)
                self.src_def_list.append(src_def)
                self.dest_def_list.append(dest_def)


        # geo2radar (WITHOUT block-by-block support)
        else:
            # parameter 1 - lalo_step (input grid)
            self.lalo_step = [float(self.src_meta['Y_STEP']),
                              float(self.src_meta['X_STEP'])]
            print('input pixel size in (lat, lon) in degree: {}'.format(self.lalo_step))

            # parameter 2 - SNWE (input grid)
            lat0 = float(self.src_meta['Y_FIRST'])
            lon0 = float(self.src_meta['X_FIRST'])
            if not self.SNWE:
                # default SNWE --> src_box
                src_box = (0, 0, int(self.src_meta['WIDTH']), int(self.src_meta['LENGTH']))
            else:
                # custom input SNWE --> src_box
                # to align SNWE to precisely to source file in geo-coord
                src_box = (int(np.rint((self.SNWE[2] - lon0) / self.lalo_step[1])),  # x0 - W
                           int(np.rint((self.SNWE[1] - lat0) / self.lalo_step[0])),  # y0 - N
                           int(np.rint((self.SNWE[3] - lon0) / self.lalo_step[1])),  # x1 - E
                           int(np.rint((self.SNWE[0] - lat0) / self.lalo_step[0])))  # y1 - S
            # src_box --> SNWE
            self.SNWE = (lat0 + self.lalo_step[0] * src_box[3],  # S - y1
                         lat0 + self.lalo_step[0] * src_box[1],  # N - y0
                         lon0 + self.lalo_step[1] * src_box[0],  # W - x0
                         lon0 + self.lalo_step[1] * src_box[2])  # E - x1
            print('input area extent in (S, N, W, E) in degree: {}'.format(self.SNWE))

            # parameter 3 - length / width (output grid)
            self.length, self.width = lut_lat.shape

            # src_lat/lon (at pixel center)
            src_len = src_box[3] - src_box[1]
            src_wid = src_box[2] - src_box[0]
            src_lat0 = self.SNWE[1] + self.lalo_step[0] * (src_box[1] + 0.5)
            src_lat1 = self.SNWE[1] + self.lalo_step[0] * (src_box[3] - 0.5)
            src_lon0 = self.SNWE[2] + self.lalo_step[1] * (src_box[0] + 0.5)
            src_lon1 = self.SNWE[2] + self.lalo_step[1] * (src_box[2] - 0.5)
            src_lat, src_lon = np.mgrid[src_lat0:src_lat1:src_len*1j,
                                        src_lon0:src_lon1:src_wid*1j]

            # parameter 4 - list of boxes & geometry definitions
            self.src_box_list = [src_box]
            self.src_def_list = [pr.geometry.GridDefinition(lons=src_lon, lats=src_lat)]
            self.dest_box_list = [(0, 0, self.width, self.length)]
            self.dest_def_list = [pr.geometry.SwathDefinition(lons=lut_lon, lats=lut_lat)]
            self.num_box = 1

        return


    def prepare_geometry_definition_geo(self):
        """Get src_def and dest_def for lookup table from Gamma and ROI_PAC."""

        def project_yx2lalo(yy, xx, SNWE, laloScale):
            """scale/project coordinates in pixel number into lat/lon
            based on bbox and step
            """
            lats = SNWE[1] + yy * laloScale[0]
            lons = SNWE[2] + xx * laloScale[1]
            lats[np.isnan(lats)] = 90.
            lats[np.isnan(lons)] = 90.
            lons[lats == 90.] = 0
            return lats, lons

        # radar2geo
        # block-by-block support is NOT implemented because:
        # writefile.write_hdf5_block requires split in dest_box for writing
        # while data_reduce.get_valid_index_from_lonlat_grid() requires split in grid (src_box)
        #   resulting in uneven / overlapped flag matrix in dest_box
        if 'Y_FIRST' not in self.src_meta.keys():

            # parameter 1 - lalo_step (output grid)
            # custom input lalo_step is not supported
            # however, it could be added by rounding the input step to the nearest
            # multiple of Y/X_STEP
            self.lalo_step = (float(self.lut_meta['Y_STEP']),
                              float(self.lut_meta['X_STEP']))
            print('output pixel size in (lat, lon) in degree: {}'.format(self.lalo_step))

            # parameter 2 - SNWE (output grid)
            lat0 = float(self.lut_meta['Y_FIRST'])
            lon0 = float(self.lut_meta['X_FIRST'])
            if not self.SNWE:
                # default SNWE --> dest_box
                dest_box = (0, 0, int(self.lut_meta['WIDTH']), int(self.lut_meta['LENGTH']))
            else:
                # custom input SNWE --> dest_box
                # to align SNWE to precisely to lookup table file in geo-coord
                dest_box = (int(np.rint((self.SNWE[2] - lon0) / self.lalo_step[1])),  # x0 - W
                            int(np.rint((self.SNWE[1] - lat0) / self.lalo_step[0])),  # y0 - N
                            int(np.rint((self.SNWE[3] - lon0) / self.lalo_step[1])),  # x1 - E
                            int(np.rint((self.SNWE[0] - lat0) / self.lalo_step[0])))  # y1 - S
                # check coverage box
                dest_box = (max(dest_box[0], 0),
                            max(dest_box[1], 0),
                            min(dest_box[2], int(self.lut_meta['WIDTH'])),
                            min(dest_box[3], int(self.lut_meta['LENGTH'])))
            # dest_box --> SNWE
            self.SNWE = (lat0 + self.lalo_step[0] * dest_box[3],  # S - y1
                         lat0 + self.lalo_step[0] * dest_box[1],  # N - y0
                         lon0 + self.lalo_step[1] * dest_box[0],  # W - x0
                         lon0 + self.lalo_step[1] * dest_box[2])  # E - x1
            print('output area extent in (S, N, W, E) in degree: {}'.format(self.SNWE))

            # parameter 3 - length / width (output grid)
            self.length = dest_box[3] - dest_box[1]
            self.width = dest_box[2] - dest_box[0]

            # parameter 4 - list of boxes & geometry definitions

            # src_y/x at pixel center
            src_len = int(self.src_meta['LENGTH'])
            src_wid = int(self.src_meta['WIDTH'])
            src_y, src_x = np.mgrid[0.5:src_len-0.5:src_len*1j,
                                    0.5:src_wid-0.5:src_wid*1j]

            # dest_y/x at pixel center
            dest_y = readfile.read(self.lut_file, datasetName='azimuthCoord', box=dest_box)[0]
            dest_x = readfile.read(self.lut_file, datasetName='rangeCoord', box=dest_box)[0]
            if 'SUBSET_XMIN' in self.src_meta.keys():
                print('input data file was cropped before.')
                dest_y[dest_y != 0.] -= float(self.src_meta['SUBSET_YMIN'])
                dest_x[dest_x != 0.] -= float(self.src_meta['SUBSET_XMIN'])

            # src/desc_y/x --> src/dest_lat/lon
            # scale y/x to the range of lat/lon to be able to use pyresample.geometryDefinition
            commMask = np.multiply(np.multiply(dest_y > 0, dest_y < src_len),
                                   np.multiply(dest_x > 0, dest_x < src_wid))
            idx_row, idx_col = np.where(commMask)
            commSNWE = (self.SNWE[1] + self.lalo_step[0] * np.max(idx_row),
                        self.SNWE[1] + self.lalo_step[0] * np.min(idx_row),
                        self.SNWE[2] + self.lalo_step[1] * np.min(idx_col),
                        self.SNWE[2] + self.lalo_step[1] * np.max(idx_col))
            laloScale = ((commSNWE[0] - commSNWE[1]) / src_len,
                         (commSNWE[3] - commSNWE[2]) / src_wid)

            src_lat, src_lon = project_yx2lalo(src_y, src_x, commSNWE, laloScale)
            dest_y[dest_y == 0.] = np.nan
            dest_x[dest_x == 0.] = np.nan
            dest_lat, dest_lon = project_yx2lalo(dest_y, dest_x, commSNWE, laloScale)

            # src_def and dest_def
            self.src_box_list = [(0, 0, src_wid, src_len)]
            self.src_def_list = [pr.geometry.GridDefinition(lons=src_lon, lats=src_lat)]
            self.dest_box_list = [dest_box]
            self.dest_def_list = [pr.geometry.SwathDefinition(lons=dest_lon, lats=dest_lat)]
            self.num_box = 1


        # geo2radar [not finished]
        else:
            raise ValueError('geo2radar with lookup table in geo-coord it NOT supported yet!')

            # provide dest_lat/dest_lon to get dest_x, dest_y
            dest_y_lt = readfile.read(self.lut_file, datasetName='azimuthCoord', box=dest_box)[0]
            dest_x_lt = readfile.read(self.lut_file, datasetName='rangeCoord', box=dest_box)[0]

            # please check this is right, I am not sure what self is. Just copy the above code
            lat0 = float(self.lut_meta['Y_FIRST'])
            lon0 = float(self.lut_meta['X_FIRST'])
            lat_step = float(self.lut_meta['Y_STEP'])
            lon_step = float(self.lut_meta['X_STEP'])

            yy = int((dest_lat - lat0)/lat_step)
            xx = int((dest_lon - lon0)/lon_step)

            row, col = yy.shape
            yy = yy.flatten()
            xx = xx.flatten()

            dest_y = dest_y_lt[yy,xx] # rows in radar coord
            dest_x = dest_x_lt[yy,xx] # column in radar coord

            dest_y = dest_y.reshape(row,col)
            dest_x = dest_x.reshape(row,col)
            # please Yunjun to finish the left part to update self

            # src_y/x
            #raise NotImplementedError('Not implemented yet for GAMMA and ROIPAC products')
        return


    @staticmethod
    def run_pyresample(src_data, src_def, dest_def, radius, interp_method='nearest', fill_value=np.nan,
                       nprocs=1, print_msg=True):
        """Resample input src_data into dest_data

        Parameters: src_data      - 2/3D np.ndarray, source data to be resampled
                    src_def       - pyresample geometry definition for source data
                    dest_def      - pyresample geometry definition for destination data
                    radius        - float, radius of influence
                    interp_method - str, interpolation method, nearest / linear
                    fill_value    - number, fill value for extrapolation
                    nprocs        - int, number of processors
        Returns:    dest_data     - 2/3D np.ndarray, source data to be resampled
        Example:    dest_data = reObj.run_pyresample(src_data, src_def, dest_def, radius,
                                                     interp_method=inps.interpMethod,
                                                     fill_value=inps.fillValue)
        """
        # get number of segments
        num_segment = int(np.ceil(dest_def.size / 1e6 + 0.5))

        # resample
        if interp_method.startswith('near'):
            if print_msg:
                msg = '{} resampling with pyresample.kd_tree '.format(interp_method)
                msg += 'using {} CPU cores in {} segments ...'.format(nprocs, num_segment)
                print(msg)
            dest_data = pr.kd_tree.resample_nearest(src_def,
                                                    src_data,
                                                    dest_def,
                                                    radius_of_influence=radius,
                                                    fill_value=fill_value,
                                                    nprocs=nprocs,
                                                    segments=num_segment,
                                                    epsilon=0.5)

        elif interp_method.endswith('linear'):
            if print_msg:
                msg = '{} resampling with pyresample.bilinear '.format(interp_method)
                msg += 'using {} CPU cores ...'.format(nprocs)
                print(msg)
            dest_data = pr.bilinear.resample_bilinear(src_data,
                                                      src_def,
                                                      dest_def,
                                                      radius=radius,
                                                      fill_value=fill_value,
                                                      neighbours=32,
                                                      nprocs=nprocs,
                                                      segments=num_segment,
                                                      epsilon=0)

        # for debug
        debug_mode = False
        if debug_mode:
            import matplotlib.pyplot as plt
            fig, ((ax11, ax12, ax13), (ax21, ax22, ax23)) = plt.subplots(nrows=2, ncols=3, figsize=(12, 8))
            dest_lats = np.array(dest_def.lats);    dest_lats[dest_lats == 90.] = np.nan
            dest_lons = np.array(dest_def.lons);    dest_lons[dest_lons == 0.] = np.nan
            im = ax11.imshow(src_def.lats);         fig.colorbar(im, ax=ax11);   ax11.set_title('src_lats')
            im = ax12.imshow(src_def.lons);         fig.colorbar(im, ax=ax12);   ax12.set_title('src_lons')
            im = ax13.imshow(src_data);             fig.colorbar(im, ax=ax13);   ax13.set_title('src_data')
            im = ax21.imshow(dest_lats);            fig.colorbar(im, ax=ax21);   ax21.set_title('dest_lats')
            im = ax22.imshow(dest_lons);            fig.colorbar(im, ax=ax22);   ax22.set_title('dest_lons')
            im = ax23.imshow(dest_data);            fig.colorbar(im, ax=ax23);   ax23.set_title('dest_data')
            plt.show()

        return dest_data


    ##------------------------------ resample using scipy.interpolate ------------------------------##

    def prepare_regular_grid_interpolator(self):
        """Prepare aux data for RGI module"""

        # radar2geo
        if 'Y_FIRST' not in self.src_meta.keys():
            # source points in regular grid
            src_len = int(self.src_meta['LENGTH'])
            src_wid = int(self.src_meta['WIDTH'])
            self.src_pts = (np.arange(src_len) + 0.5,
                            np.arange(src_wid) + 0.5)

            # destination points
            dest_y = readfile.read(self.lut_file, datasetName='azimuthCoord')[0]
            dest_x = readfile.read(self.lut_file, datasetName='rangeCoord')[0]
            if 'SUBSET_XMIN' in self.src_meta.keys():
                print('input data file was cropped before.')
                dest_y[dest_y != 0.] -= float(self.src_meta['SUBSET_YMIN'])
                dest_x[dest_x != 0.] -= float(self.src_meta['SUBSET_XMIN'])
            self.interp_mask = np.multiply(np.multiply(dest_y > 0, dest_y < src_len),
                                           np.multiply(dest_x > 0, dest_x < src_wid))
            self.dest_pts = np.hstack((dest_y[self.interp_mask].reshape(-1, 1),
                                       dest_x[self.interp_mask].reshape(-1, 1)))

            # parameter 1/2/3 - lalo_step / SNWE / length & width
            lat_num = int(self.lut_meta['LENGTH'])
            lon_num = int(self.lut_meta['WIDTH'])
            lat_step = float(self.lut_meta['Y_STEP'])
            lon_step = float(self.lut_meta['X_STEP'])
            lat0 = float(self.lut_meta['Y_FIRST'])
            lon0 = float(self.lut_meta['X_FIRST'])
            self.lalo_step = (lat_step, lon_step)
            self.SNWE = (lat0 + lat_step * lat_num,
                         lat0,
                         lon0,
                         lon0 + lon_step * lon_num)
            self.length = lat_num
            self.width = lon_num

            # parameter 4 - list of boxes
            self.src_box_list = [(0, 0, src_wid, src_len)]
            self.dest_box_list = [(0, 0, lon_num, lat_num)]

            # WITHOUT block-by-block support
            # however, date-by-date resampling is already used
            # one could add date-by-date IO in geocoded.py to memory efficiency


        # geo2radar
        else:
            raise ValueError('geo2radar with lookup table in geo-coord it NOT supported yet!')


    def run_regular_grid_interpolator(self, src_data, interp_method='nearest', fill_value=np.nan, print_msg=True):
        """Interpolate 2D matrix"""
        # prepare interpolation function
        rgi_func = RGI(self.src_pts,
                       src_data,
                       method=interp_method,
                       bounds_error=False,
                       fill_value=fill_value)

        # prepare output matrix
        geo_data = np.empty((self.length, self.width), src_data.dtype)
        geo_data.fill(fill_value)

        # interpolate output matrix
        geo_data[self.interp_mask] = rgi_func(self.dest_pts)
        return geo_data