resample.py

"""Class wrapped around pyresample/scipy for resampling."""
############################################################
# 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

with warnings.catch_warnings():
    warnings.simplefilter('ignore', UserWarning)
    import pyresample as pr

from mintpy.constants import EARTH_RADIUS
from mintpy.objects.cluster import split_box2sub_boxes
from mintpy.utils import ptime, readfile, utils0 as ut


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(f'resampling software: {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
        """
        vprint = print if print_msg else lambda *args, **kwargs: None

        # 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
            vprint('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
            vprint('input source data is NOT float, change fill_value from NaN to 0.')

        # run resampling
        kwargs = dict(
            interp_method=self.interp_method,
            fill_value=fill_value,
            print_msg=print_msg,
        )

        if self.software == 'pyresample':
            # move 1st/time dimension to the last
            # so that rows/cols axis are the first, 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,
                nprocs=self.nprocs,
                **kwargs,
            )

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

        else:
            vprint(f'{self.interp_method} resampling using scipy.interpolate.RegularGridInterpolator ...')
            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]):
                    prog_bar.update(i+1)
                    dest_data[i, :, :] = self.run_regular_grid_interpolator(src_data=src_data[i, :, :], **kwargs)
                prog_bar.close()
            else:
                dest_data = self.run_regular_grid_interpolator(src_data=src_data, **kwargs)

        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 split
        """
        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
            atr = readfile.read_attribute(src_file)
            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:
                max_ds_size = int(atr['LENGTH']) * int(atr['WIDTH'])

            # scale the size for complex arrays
            if atr.get('DATA_TYPE', 'float32').startswith('complex'):
                max_ds_size *= 6 if atr['DATA_TYPE'].endswith('128') else 3

            # calc num_box
            num_box = int(np.ceil((max_ds_size * 4 * 6) / (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


    def get_lat_lon_step(self, src_lat0, src_lat1, src_lon0, src_lon1):
        """Get/check the lat/lon step size for geocoding.

        Approach 1: ensure the same pixel area before / after resampling
            Treat the pixel in radar coordinates as an rotated rectangle. Use the bounding
                box of the rotated rectangle for the ratio between lat and lon steps. Then
                scale the lat and lon step size to ensure the same area between the pixels
                in radar and geo coordinates.
            This is recommended, but it requires metadata on the pixel size info.
            Link: https://math.stackexchange.com/questions/4001034

        Approach 2: ensure the same matrix shape before / after resampling
            This is the backup approach if approach 1 does not work.

        Parameters: src_lat0/1   - float, max/min latitude  in degree of the input data file
                    src_lon0/1   - float, min/max longitude in degree of the input data file
        Returns:    lat/lon_step - float, output step size in latitude/longitude in degree
        """

        if self.lalo_step is not None:
            # use input lat/lon step
            lat_step, lon_step = self.lalo_step

        else:
            # calculate default lat/lon step
            # approach 1: ensure the same pixel area
            print('calculate output pixel size using approach 1 '
                  '(same pixel area before/after resampling)')

            # check required metadata
            meta = {**self.lut_meta, **self.src_meta}
            key_list = [
                'AZIMUTH_PIXEL_SIZE',
                'HEADING',
                'HEIGHT',
                'RANGE_PIXEL_SIZE',
                'STARTING_RANGE',
                'WIDTH',
            ]
            key_not_found = [x for x in key_list if x not in meta.keys()]
            if len(key_not_found) == 0:

                # azimuth angle (rotation angle) in radian
                az_angle = np.deg2rad(abs(ut.heading2azimuth_angle(float(meta['HEADING']))))

                # radar pixel size in meter
                az_step = ut.azimuth_ground_resolution(meta)
                rg_step = ut.range_ground_resolution(meta)

                # geo pixel size in meter
                x_step = rg_step * abs(np.cos(az_angle)) + az_step * abs(np.sin(az_angle))
                y_step = rg_step * abs(np.sin(az_angle)) + az_step * abs(np.cos(az_angle))
                scale_factor = np.sqrt((rg_step * az_step) / (x_step * y_step))
                x_step *= scale_factor
                y_step *= scale_factor

                # geo pixel size in degree
                lat_c = (src_lat0 + src_lat1) / 2.
                lon_step = np.rad2deg(x_step / (EARTH_RADIUS * np.cos(np.deg2rad(lat_c))))
                lat_step = np.rad2deg(y_step / EARTH_RADIUS) * -1.

            else:
                # approach 2: ensure the same matrix shape
                msg = 'WARNING: NOT all required metadata are found. '
                msg += f'Missing metadata: {key_not_found}. Switch to approach 2.\n'
                msg += 'calculate output pixel size using approach 2 '
                msg += '(same matrix shape before/after resampling)'
                print(msg)

                # ensure the same matrix shape before / after geocoding
                # if not enough metadata found for the above
                lut_len = int(self.lut_meta['LENGTH'])
                lut_wid = int(self.lut_meta['WIDTH'])
                lat_step = (src_lat1 - src_lat0) / (lut_len - 1)
                lon_step = (src_lon1 - src_lon0) / (lut_wid - 1)

        # ensure lat/lon step sign
        lat_step = abs(lat_step) * -1.
        lon_step = abs(lon_step)

        return lat_step, lon_step


    ##------------------------------ 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 find_valid_lat_lon(lat, lon):
            """mask pixels with abnormal values (0, etc.)
            This is found on sentinelStack multiple swath lookup table file.

            Parameters: lat/lon - 2D np.ndarray in float32, latitude/longitude in degrees.
            Returns:    lat/lon - 2D np.ndarray in float32, latitude/longitude in degrees.
                        mask    - 2D np.ndarray in bool, 1/0 for valid/invalid pixels.
            """
            # 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 continuous 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(f'read latitude / longitude from lookup table file: {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 = find_valid_lat_lon(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)
            self.lalo_step = self.get_lat_lon_step(src_lat0, src_lat1, src_lon0, src_lon1)
            print(f'output pixel size in (lat, lon) in degree: {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(f'output area extent in (S, N, W, E) in degree: {self.SNWE}')
            print(f'output file row / column number: ({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)
            # and update num_box based on the actual dest_box_list
            self.dest_box_list, self.num_box = 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):
                msg = f'[{i+1}/{self.num_box}] preparing geometry for dest_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 = abs(src_lat1 - src_lat0) * abs(src_lon1 - src_lon0)
                dest_area = abs(lat1 - lat0) * abs(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)
                    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))
                    msg += f' --> reduced src_box: {src_box}'

                else:
                    src_box = (0, 0, lut_lat.shape[1], lut_lat.shape[0])
                    msg += f' --> full src_box: {src_box}'
                print(msg)

                # 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(f'input pixel size in (lat, lon) in degree: {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(f'input area extent in (S, N, W, E) in degree: {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(f'output pixel size in (lat, lon) in degree: {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(f'output area extent in (S, N, W, E) in degree: {self.SNWE}')

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

            # 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 = f'{interp_method} resampling with pyresample.kd_tree '
                msg += f'using {nprocs} CPU cores in {num_segment} segments ...'
                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 = f'{interp_method} resampling with pyresample.bilinear '
                msg += f'using {nprocs} CPU cores ...'
                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 RegularGridInterpolator 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
        interp_func = RegularGridInterpolator(
            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] = interp_func(self.dest_pts)

        return geo_data