B
    Qa4                 @   s   yd dl ZW n ek
r(   edY nX d dlZd dlZd dlZd dlmZ d dl	m
Z d dlmZ d dlmZmZmZ dZG dd	 d	ZdS )
    NzCan not import pyresample!)ndimage)RegularGridInterpolator)split_box2sub_boxes)readfileptimeutils0gTXAc               @   s   e Zd ZdZddddejddddddfdd	Zd
d Zdd Zd ddZ	e
d!ddZe
d"ddZdd Zdd Ze
dejddfddZdd ZdejdfddZdS )#resamplea|	  
    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)
    Nnearest   r   
pyresampleTc             C   sX   || _ || _|| _|| _|| _|| _|| _|| _|| _|| _	|	| _
|
| _d| _d| _dS )a  
        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.
        r
   N)lut_filelat_filelon_filesrc_fileSNWE	lalo_stepinterp_method
fill_valuenprocs
max_memorysoftware	print_msgnum_boxradius)selfr   r   r   r   r   r   r   r   r   r   r   r    r   M/home/centos/operations/rsmas_insar/sources/MintPy/mintpy/objects/resample.py__init__S   s    zresample.__init__c             C   sX   | j rt| j | _nd| _| jr2t| j| _nd| _| jdkrT| | j| j| _	dS )zRead metadata

        Note: separate open() from prepare() to handle numpy matrix directly w/o src_file
        by assigning src_meta after open() and before prepare()
        Nr   )
r   r   read_attributelut_metar   src_metar   get_num_boxr   r   )r   r   r   r   open~   s    
zresample.openc             C   s   | j dks| jdkrtdtd| j | jdkrXd| j  krN|   qtdn<| jdkrd| j  krz|   n| 	  | 
| j| j| _dS )z/Prepare aux data before interpolation operationNz-lookup table or source data metadata is None!zresampling software: {}scipyY_FIRSTzXresampling using scipy with lookup table in radar-coord (ISCE / DORIS) is NOT supported!r   )r   r    
ValueErrorprintformatr   keys!prepare_regular_grid_interpolatorprepare_geometry_definition_geo!prepare_geometry_definition_radarget_radius_of_influencer   r   )r   r   r   r   prepare   s    




zresample.preparec       	   
   C   s  | j }tjtjtjtjtjtjg}|jtj	kr@d}|rdt
d n$|j|krdt|rdd}|rdt
d | jdkrt|jdkrt|dd}| j|| j| | j| | j| j|| j| jd}t|jdkrt|dd}n|rt
d	| j t|jdkrt|jd | j| jf|j}tj|jd |d
}x^t|jd D ]L}| j||ddddf | j|dd||ddddf< ||d  q<W |   n| j|| j|dd}|S )a  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
        Fz=input source data is bool type, restrict fill_value to False.r   z@input source data is NOT float, change fill_value from NaN to 0.r      )src_datasrc_defdest_defr   r   r   r   r   zA{} resampling using scipy.interpolate.RegularGridInterpolator ...)maxValuer   NT)r0   r   r   r   r
   )!r   npsingledouble
longdoublecsinglecdoubleclongdoubledtypebool_r&   isnanr   lenshapemoveaxisrun_pyresamplesrc_def_listdest_def_listr   r   r   r   r'   emptylengthwidthr   progressBarrangerun_regular_grid_interpolatorupdateclose)	r   r0   box_indr   r   Zfloat_types	dest_dataprog_barir   r   r   run_resample   sP    	



zresample.run_resample      @c          	      s   d}| rt j| r|dkrt j| d }|dkrvt| d.  fdd  D }tdd |D }W dQ R X n"t	| }t
|d	 t
|d
  }t
t|d d |d  }|S )a  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
        r
   r   )z.h5z.he5rc                s&   g | ]}t  | tjr | jqS r   )
isinstanceh5pyDatasetr?   ).0rO   )fr   r   
<listcomp>  s    z(resample.get_num_box.<locals>.<listcomp>c             S   s   g | ]}t |qS r   )r4   prod)rV   rO   r   r   r   rX     s    NLENGTHWIDTH   i   @)ospathisfilesplitextrT   Filer(   maxr   r   intr4   ceil)r   r   Z	scale_facr   fextZ	ds_shapesZmax_ds_sizeatrr   )rW   r   r!      s    	
zresample.get_num_boxr.   c             C   s>   d|  krd}n(tt| }|tj d t }|| }|S )zQGet radius of influence based on the lookup table resolution in lat/lon directionr$   g     j@g     f@)r(   rb   r4   abspiEARTH_RADIUS)r   r    ratior   Zstep_degZstep_mr   r   r   r,     s    z resample.get_radius_of_influencec       !   	   C   s*  dd }t d| j | jr$| jn| j}| jr6| jn| j}tj|ddd tj	}tj|ddd tj	}|||\}}}d| j
 krt|| }t|| }t|| }	t|| }
| jd	kr|| |jd d
  |
|	 |jd
 d
  f| _n&t| jd d t| jd
 d f| _t d| j | jd	kr|| jd d  || jd d  |	| jd
 d  |
| jd
 d  f| _tt| jd | jd
  | jd  | _tt| jd | jd  | jd
  | _| jd
 | jd | j  | jd
 | jd | jd | jd
 | j  f| _t d| j t d| j| j g | _g | _g | _g | _tdd| j| jf| jddd| _xt| jD ]\}}| jd
krt d|d
 | j| |d |d
  }|d |d  }| jd
 | jd |d
 d   }| jd
 | jd |d d   }| jd | jd
 |d d   }| jd | jd
 |d d   }tj|||d |||d f \}}|| |
|	  }|| ||  }||d k rt d tjj ||||dd}t!|\}}t"|t"|t#|t#|f}ndd|jd
 |jd f}tj$j%||d
 |d |d |d f ||d
 |d |d |d f d}tj$j&||d}| j'| | j'| | j'| qW n|t(| j
d t(| j
d g| _t d| j t(| j
d }t(| j
d  }| jsddt| j
d! t| j
d" f}ntt| jd | | jd
  tt| jd
 | | jd  tt| jd | | jd
  tt| jd | | jd  f}|| jd |d   || jd |d
   || jd
 |d   || jd
 |d   f| _t d#| j |j\| _| _|d |d
  }|d |d  }| jd
 | jd |d
 d   }| jd
 | jd |d d   }| jd | jd
 |d d   }	| jd | jd
 |d d   }
tj|||d |	|
|d f \}} |g| _tj$j&| |dg| _dd| j| jfg| _tj$j%||dg| _d
| _d	S )$zKGet src_def and dest_def for lookup table in radar-coord (from ISCE, DORIS)c             S   s2  t | dk|dk}t |t j}x| |gD ]}t j|| dd\}}xt |t |d kr
tj|dd}t	
||kd }t ||t | kd }	|d |d  }
||	d  |
d	  }||	d
 d  |
d	  }|t ||k||k9 }t j|| dd\}}qHW q,W d| |dk< d||dk< | ||fS )zmask pixels with abnormal values (0, etc.)
            This is found on sentinelStack multiple swath lookup table file.
            g        
   )binsg333333?r.   )cutoffr   r
   g       @r/   g     V@)r4   multiplyarrayr<   	histogramrb   sumutmedian_abs_deviation_thresholdr   labelwhereargmax)latlon	zero_maskmaskdataZ	bin_valueZbin_edgeZbin_value_thresZ	bin_labelidxZbin_stepZd_minZd_maxr   r   r   mark_lat_lon_anomoly"  s     zHresample.prepare_geometry_definition_radar.<locals>.mark_lat_lon_anomolyz4read latitude / longitude from lookup table file: {}latitude)datasetNamer   	longituder$   Nr
   g      g      ?z-output pixel size in (lat, lon) in degree: {}g       @r.      z0output area extent in (S, N, W, E) in degree: {}z)output file row / column number: ({}, {})yT)box	num_split	dimensionr   z)preparing geometry for dest_box {}/{}: {}g      ?y              ?z0searching relevant box covering the current SNWEi  )radius_of_influence)lonslatsY_STEPX_STEPz,input pixel size in (lat, lon) in degree: {}X_FIRSTr[   rZ   z/input area extent in (S, N, W, E) in degree: {}))r&   r'   r   r   r   r   readastyper4   float32r    r(   nanmaxnanminr   r?   rg   r   rc   rintrE   rF   src_box_listrB   dest_box_listrC   r   r   	enumeratemgridprZdata_reduceZ get_valid_index_from_lonlat_gridru   minrb   geometrySwathDefinitionGridDefinitionappendfloat)!r   r}   r   r   Zlut_latZlut_lonrz   Zsrc_lat0Zsrc_lat1Zsrc_lon0Zsrc_lon1rO   dest_boxlat_numlon_numlat0lat1lon0lon1dest_latdest_lonZsrc_areaZ	dest_areaflagidx_rowidx_colsrc_boxr1   r2   src_lensrc_widsrc_latsrc_lonr   r   r   r+     s    "**
    "   $    z*resample.prepare_geometry_definition_radarc          	   C   s  dd }d| j  krt| jd t| jd f| _td| j t| jd }t| jd }| jsddt| jd	 t| jd
 f}ntt	
| jd | | jd  tt	
| jd | | jd  tt	
| jd | | jd  tt	
| jd | | jd  f}t|d dt|d dt|d t| jd	 t|d t| jd
 f}|| jd |d   || jd |d   || jd |d   || jd |d   f| _td| j |d |d  | _|d |d  | _t| j d
 }t| j d	 }t	jd|d |d d|d |d f \}}tj| jd|dd }	tj| jd|dd }
d| j  krtd |	|	dk  t| j d 8  < |
|
dk  t| j d 8  < t	t	|	dk|	|k t	|
dk|
|k }t	|\}}| jd | jd t	|  | jd | jd t	|  | jd | jd t	|  | jd | jd t	|  f}|d |d  | |d |d  | f}|||||\}}t	j|	|	dk< t	j|
|
dk< ||	|
||\}}dd||fg| _tjj||dg| _|g| _tjj||dg| _d| _ntdtj| jd|dd }tj| jd|dd }t| jd }t| jd }t| jd }t| jd }t|| | }t|| | }|j \}}|! }|! }|||f }	|||f }
|	"||}	|
"||}
dS )zAGet src_def and dest_def for lookup table from Gamma and ROI_PAC.c             S   sX   |d | |d   }|d ||d   }d|t |< d|t |< d||dk< ||fS )zfscale/project coordinates in pixel number into lat/lon
            based on bbox and step
            r
   r   r   g     V@)r4   r=   )yyxxr   	laloScaler   r   r   r   r   project_yx2lalo  s    zAresample.prepare_geometry_definition_geo.<locals>.project_yx2lalor$   r   r   z-output pixel size in (lat, lon) in degree: {}r   r   r[   rZ   r   r
   r.   z0output area extent in (S, N, W, E) in degree: {}g      ?y              ?azimuthCoord)r   r   
rangeCoordSUBSET_XMINz#input data file was cropped before.g        SUBSET_YMIN)r   r   z>geo2radar with lookup table in geo-coord it NOT supported yet!N)#r    r(   r   r   r   r&   r'   r   rc   r4   r   rb   r   rE   rF   r   r   r   r   rn   ru   nanr   r   r   r   rB   r   r   rC   r   r%   r?   flattenreshape)r   r   r   r   r   r   r   Zsrc_yZsrc_xdest_ydest_xZcommMaskr   r   ZcommSNWEr   r   r   r   r   Z	dest_y_ltZ	dest_x_ltlat_steplon_stepr   r   rowcolr   r   r   r*     s    "   $ 
z(resample.prepare_geometry_definition_geoc             C   s  t t|jd d }|drf|rHd|}	|	d||7 }	t|	 tjj	|| |||||dd}
nL|
dr|rd|}	|	d	|7 }	t|	 tjj| ||||d
||dd	}
d}|rddlm} |jdddd\}\\}}}\}}}t|j}tj||dk< t|j}tj||dk< ||j}|j||d |d ||j}|j||d |d || }|j||d |d ||}|j||d |d ||}|j||d |d ||
}|j||d |d |  |
S )a  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)
        g    .Ag      ?Znearz&{} resampling with pyresample.kd_tree z%using {} CPU cores in {} segments ...)r   r   r   segmentsepsilonlinearz'{} resampling with pyresample.bilinear zusing {} CPU cores ...    r   )r   r   Z
neighboursr   r   r   FNr   r.   )      )nrowsncolsfigsizeg     V@g        )axZsrc_latsZsrc_lonsr0   	dest_lats	dest_lonsrM   )rc   r4   rd   size
startswithr'   r&   r   Zkd_treeZresample_nearestendswithbilinearZresample_bilinearmatplotlib.pyplotpyplotsubplotsro   r   r   r   imshowcolorbar	set_titleshow)r0   r1   r2   r   r   r   r   r   Znum_segmentmsgrM   
debug_modepltfigZax11Zax12Zax13Zax21Zax22Zax23r   r   imr   r   r   rA   ^  sp    




$    
  

  

  

  

  
zresample.run_pyresamplec             C   s  d| j  krt| j d }t| j d }t|d t|d f| _tj| jddd }tj| jddd }d	| j  krt	d
 ||dk  t
| j d 8  < ||dk  t
| j d	 8  < tt|dk||k t|dk||k | _t|| j dd|| j ddf| _t| jd }t| jd }t
| jd }t
| jd }t
| jd }	t
| jd }
||f| _|	||  |	|
|
||  f| _|| _|| _dd||fg| _dd||fg| _ntddS )zPrepare aux data for RGI moduler$   rZ   r[   g      ?r   )r   r   r   r   z#input data file was cropped before.g        r   r/   r
   r   r   r   z>geo2radar with lookup table in geo-coord it NOT supported yet!N)r    r(   rc   r4   arangesrc_ptsr   r   r   r&   r   rn   interp_maskhstackr   dest_ptsr   r   r   rE   rF   r   r   r%   )r   r   r   r   r   r   r   r   r   r   r   r   r   r   r)     s>    

	z*resample.prepare_regular_grid_interpolatorc             C   sH   t | j||d|d}t| j| jf|j}|| || j|| j	< |S )zInterpolate 2D matrixF)methodbounds_errorr   )
RGIr   r4   rD   rE   rF   r;   fillr   r   )r   r0   r   r   r   Zrgi_funcZgeo_datar   r   r   rI     s    
z&resample.run_regular_grid_interpolator)r   T)Nr   rQ   )r.   )__name__
__module____qualname____doc__r4   r   r   r"   r-   rP   staticmethodr!   r,   r+   r*   rA   r)   rI   r   r   r   r   r      s&   7*
@ < D4r   )r   r   ImportErrorr]   rT   numpyr4   r#   r   scipy.interpolater   r   mintpy.objects.clusterr   mintpy.utilsr   r   r   rr   ri   r   r   r   r   r   <module>
   s   