U
    Qa8                     @   s   d dl Z d dlZd dlZ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mZ dZdZdZd	d
 ZdddZdd Zdd Zdd Zdd Zdd ZdddZedkreej dd  dS )    N)resize)RegularGridInterpolator)
timeseries)ptimereadfile	writefileutilsa  references:
  Yu, C., Li, Z., Penna, N. T., & Crippa, P. (2018). Generic atmospheric correction model for Interferometric
    Synthetic Aperture Radar observations. Journal of Geophysical Research: Solid Earth, 123(10), 9202-9222.
  Yu, C., Li, Z., & Penna, N. T. (2018). Interferometric synthetic aperture radar atmospheric correction
    using a GPS-based iterative tropospheric decomposition model. Remote Sensing of Environment, 204, 109-121.
z--dir ./GACOS
  20060624.ztd
  20060624.ztd.rsc
  20061225.ztd
  20061225.ztd.rsc
  ...
  OR
  20060624.ztd.tif
  20061225.ztd.tif
  ...
zexample:
  tropo_gacos.py -f timeseries.h5 -g inputs/geometryRadar.h5 --dir ./GACOS
  tropo_gacos.py -f geo/geo_timeseries.h5 -g geo/geo_geometryRadar.h5 --dir ./GACOS
c                  C   sr   t jdt jtd t d t d} | jdddddd	 | jd
ddddd	 | jdddddd | jdddd | S )NzBTropospheric correction using GACOS (http://www.gacos.net) delays

)descriptionformatter_classepilogz-fz--filedis_fileTz(timeseries HDF5 file, i.e. timeseries.h5)destrequiredhelpz-gz--geom	geom_filezgeometry file.z--dirz--GACOS-dir	GACOS_dirz./GACOSzAdirectory to downloaded GACOS delays data (default: %(default)s).)r   defaultr   -ocor_dis_filez5Output file name for trospheric corrected timeseries.)r   r   )argparseArgumentParserRawTextHelpFormatter	REFERENCEDIR_DEMOEXAMPLEadd_argument)parser r   H/home/centos/operations/rsmas_insar/sources/MintPy/mintpy/tropo_gacos.pycreate_parser1   s"    r    c                 C   s  t  }|j| d}tj|j|_td|j dD ].}t|| }|r2tj|s2t	d
|q2t|j}t|j}d| krdnd}d| krdnd}|dd	}	|dkr|	d
krd
|	}
|
d7 }
t|
||krDtdd |j|jfD }d}
|
dj
tj|j||d7 }
|
dj
tj|j||d7 }
t|
tjtj|jd|_|js||jdd d |_|S )zCommand line parser.)argsz Use GACOS products at directory:)r   r   zinput file not exist: {}Y_FIRSTgeoradar	PROCESSORisce)gammaroipacz*Radar-coded file from {} is NOT supported!zT
    Try to geocode the time-series and geometry files and re-run with them instead.c                 s   s   | ]}t tj|V  qd S N)lenospathbasename.0ir   r   r   	<genexpr>]   s     z!cmd_line_parse.<locals>.<genexpr>zCInput time-series and geometry file are NOT in the same coordinate!z"
    file {f:<{n}} coordinate: {c})fnczGACOS.h5.r   z	_GACOS.h5)r    
parse_argsr+   r,   abspathr   printvarsisfileFileNotFoundErrorformatr   read_attributer   r   keysget
ValueErrormaxr-   joindirname
tropo_filer   split)iargsr   inpskeyfnameatr1atr2coord1coord2procmsgr3   r   r   r   cmd_line_parseB   s6    

rP   c           	      C   s   t |d t |d  }}t|jdd||fd}t| }t||}tj| |dd }t|||fddddd	}|| }|d
9 }|S )aW  calc single path tropo delay in line-of-sight direction

    Parameters: ztd_file      - str, path of zenith delay file
                atr           - dict, dictionary of attribute for output file
                cos_inc_angle - 2D np.ndarray in float32, cos(inc_angle)
    Returns:    delay         - 2D np.ndarray in float32, LOS delay
    LENGTHWIDTHr   )	pixel_box)box   constantT)ordermodeanti_aliasingpreserve_range)	intut
coordinatebox_pixel2geor   r=   box_geo2pixelreadr   )	ztd_fileatrcos_inc_anglelengthwidthgeo_boxatr_ztdpix_boxdelayr   r   r   get_delay_geol   s    


rk   c           
      C   s~   t | \}}t|}tj|dd\}}t|}| | f}t||dddd}||}	|	|j	}	|	| }	|	d9 }	|	S )a`  calc single path tropo delay in line-of-sight direction

    Parameters: ztd_file      - str, path of zenith delay file
                cos_inc_angle - 2D np.ndarray in (len, wid) in float32, cos(inc_angle)
                pts_new       - 2D np.ndarray in (len*wid, 2) in float32
    Returns:    delay         - 2D np.ndarray in float32, LOS delay
    rU   )	dimensionnearestFr   )methodbounds_error
fill_valuer[   )
r   ra   npflipudr]   get_lat_lonflattenRGIreshapeshape)
rb   rd   pts_newZ	delay_ztdrh   latslonsZpts_ztdZRGI_funcrj   r   r   r   get_delay_radar   s"    	

r{   c                    s  t |}|d }|dkr(t| }n6|dkrPt |d }t|d}ntd|g }t	j
t|t	jd}	t|D ]\\}
 fdd	d
D }dd	 |D }t|dkr||d  q~td d|	|
< q~t	|	dkrt	||	  }dd fdd}||| |dkr&dS d|d< d|d< dD ]}|| kr:|| q:t|d t|d  }}t|}t	j|t	jd}|j|f|gt	j|||fdgd}tj| ||d td| t j|ddd }t	|t	j d }d | krd}n6td! t||\}}t	 |!d"d#|!d"d#f}tj"|d$}t#|D ]}
||
 ||
 }d | krt$|||}nt%|||}|
|
d# d|d|g}tj&| |d|dd% |j'|
d# t(j)*|d& qP|+  | S )'z2calculate delay time-series and write to HDF5 file	FILE_TYPEr   .unwDATE12-z'un-supported displacement file type: {})dtypec              	      s"   g | ]}t j d |qS )z{}{})r+   r,   rB   r<   )r/   fext)r   date_strr   r   
<listcomp>   s     z.calculate_delay_timeseries.<locals>.<listcomp>)z.ztdz.ztd.tifc                 S   s   g | ]}t j|r|qS r   )r+   r,   exists)r/   r2   r   r   r   r      s      r   z7WARNING: NO ztd file found for {}! Continue without it.Fc                 S   s   t | }|d |d fS )NrQ   rR   )r   r=   )rI   rc   r   r   r   get_dataset_size   s    
z4calculate_delay_timeseries.<locals>.get_dataset_sizec              	      s   t d t d  d}tj | dddkr<d}t d nt d d	d
 | D } |ksxt fdd|D rd}t d| nZt d t d@}t|d dd d d d f dkrd}t d nt d W 5 Q R X t d| |S )Nzupdate mode: ONzoutput file: {}skipF)out_filein_file	print_msgrunzF1) output file either do NOT exist or is NOT newer than all ZTD files.z61) output file exists and is newer than all ZTD files.c                 S   s    g | ]}t td |d qS )z\d{8}r   )strrefindallr.   r   r   r   r      s     zCcalculate_delay_timeseries.<locals>.run_or_skip.<locals>.<listcomp>c                 3   s   | ]}|t   kV  qd S r)   )r   get_date_listr.   rD   r   r   r1      s     zBcalculate_delay_timeseries.<locals>.run_or_skip.<locals>.<genexpr>zc2) output file does NOT have the same len/wid as the geometry file {} or does NOT contain all dateszO2) output file has the same len/wid as the geometry file and contains all datesrr   r[   r   z$3) output file is NOT fully written.z 3) output file is fully written.zrun or skip: {})	r8   r<   r]   run_or_skipanyh5pyFilerq   all)	ztd_filesrD   r   flag	date_listr2   )r   r   r   r      s.    
$
z/calculate_delay_timeseries.<locals>.run_or_skipr   NmUNIT)REF_DATEREF_XREF_YREF_LATREF_LONrQ   rR   )dater   )metadataz!read incidenceAngle from file: {}incidenceAngledatasetNameg     f@r"   z&get pixel coordinates in geometry filer[   rU   )maxValue)datar   blockr   )suffix),r   r=   r   r   r   yyyymmddrE   r@   r<   rq   onesr*   bool_	enumerateappendr8   r   arraytolistr>   popr\   string_r   float32r   layout_hdf5ra   cospir]   rs   hstackrv   progressBarrangerk   r{   write_hdf5_blockupdater+   r,   r-   close)rD   r   r   r   rc   ftyper   date12r   r   r0   fnamesr   rH   re   rf   num_datedatesds_name_dict	inc_anglerd   rx   ry   rz   prog_barrb   rj   r   r   )r   r   r   r   calculate_delay_timeseries   s|    

#
r   c                 C   sH   t d t d ddlm} | |d|dg}t dd| || |S )	NO
------------------------------------------------------------------------------z=correcting relative delay for input time-series using diff.pyr   )diffr   z--forcezdiff.py )r8   mintpyr   rB   main)r   rD   r   r   rF   r   r   r   correct_timeseriesE  s    
r   c                 C   s   t d t d t d|  tj| dd\}}t|d d\}}t d||| tj|d	|dd
 }|tj|d	|dd
 8 }|dtj t	|d  9 }t d| t
|| || |S )Nr   z1correcting relative delay for input interferogramzread data from {}phaser   r~   r   z)calc tropospheric delay for {}-{} from {}ztimeseries-{}r   g      
WAVELENGTHzwrite corrected data to {})r8   r<   r   ra   r   r   rE   rq   r   floatr   write)r   rD   r   r   rc   date1date2Ztropor   r   r   correct_single_ifgramR  s    r   c                 C   s~   t | }t|j|j|j|jd t|jd }|dkrNt|j|j|j	d n,|dkrlt
|j|j|j	d ntd| d S )N)rD   r   r   r   r|   r   )r   rD   r   r}   zJinput file {} is not timeseries nor .unw, correction is not supported yet.)rP   r   rD   r   r   r   r   r=   r   r   r   r8   r<   )rF   rG   r   r   r   r   r   e  s&    r   __main__rU   )N)N)!r+   sysr   r   r   numpyrq   skimage.transformr   scipy.interpolater   ru   mintpy.objectsr   mintpy.utilsr   r   r   r   r]   r   r   r   r    rP   rk   r{   r   r   r   r   __name__argvr   r   r   r   <module>	   s.   
*$% 
