B
    Qa                 @   s&  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mZmZmZ d dlmZmZmZmZmZmZ ejZdZddd	d
dgZe	dZdZdZdZdd Zd(ddZ dd Z!d)ddZ"dd Z#dd Z$dd Z%d*ddZ&d d! Z'd+d"d#Z(d,d$d%Z)e*d&kr"e)ej+d'd  dS )-    N)linalg)get_template_content)
timeseriesgiantTimeseriesHDFEOScluster)	arg_groupptime	time_funcreadfile	writefileutilszmintpy.velocity.	startDateendDateexcludeDate	bootstrapbootstrapCountvelocitya  references:
  Fattahi, H., and F. Amelung (2015), InSAR bias and uncertainty due to the systematic and stochastic
  tropospheric delay, Journal of Geophysical Research: Solid Earth, 120(12), 8758-8773, doi:10.1002/2015JB012419.

  Efron, B., and R. Tibshirani (1986), Bootstrap methods for standard errors, confidence intervals,
  and other measures of statistical accuracy, Statistical science, 54-75, doi:10.1214/ss/1177013815.
a  example:
  timeseries2velocity.py  timeseries_ERA5_demErr.h5
  timeseries2velocity.py  timeseries_ERA5_demErr_ramp.h5  -t KyushuT73F2980_2990AlosD.template
  timeseries2velocity.py  timeseries.h5  --start-date 20080201  --end-date 20100508
  timeseries2velocity.py  timeseries.h5  --exclude-date exclude_date.txt

  timeseries2velocity.py  LS-PARAMS.h5
  timeseries2velocity.py  NSBAS-PARAMS.h5
  timeseries2velocity.py  TS-PARAMS.h5

  # bootstrapping for STD calculation
  timeseries2velocity.py timeseries_ERA5_demErr.h5 --bootstrap

  # complex time functions
  timeseries2velocity.py timeseries_ERA5_ramp_demErr.h5 --poly 3 --period 1 0.5 --step 20170910
  timeseries2velocity.py timeseries_ERA5_demErr.h5      --poly 1 --exp 20170910 90
  timeseries2velocity.py timeseries_ERA5_demErr.h5      --poly 1 --log 20170910 60.4
  timeseries2velocity.py timeseries_ERA5_demErr.h5      --poly 1 --log 20170910 60.4 200 --log 20171026 200.7
z-exclude_date.txt:
20040502
20060708
20090103
c              C   sh  t jdt jtd t d t d} | jddd | jddd	d
d | jdddd | jddddd | jddddd | jdddtddd | jddd d!d" | d#}|jd$d%d&d'd |jd(d)d*d+d |jd,d-d.d/g d0t	 d1 | d2d3}|jd4d5d6dd7d |jd8d9d:td;d<d= t
| } | d>d?}|jd@dAdBddCd |jdDdEdFdGdHdI t
| } | S )JNz4Estimate velocity / time functions from time-series.
)descriptionformatter_classepilogtimeseries_filez(Time series file for velocity inversion.)helpz
--templatez-ttemplate_fileztemplate file with options)destr   z--ts-std-filets_std_filez2Time-series STD file for velocity STD calculation.z-oz--outputoutfilezoutput file namez--updateupdate_mode
store_truezEnable update mode, and skip estimation if:
1) output velocity file already exists, readable and newer than input file
2) all configuration parameters are the same.)r   actionr   z--ref-yxref_yx)YX   z)Change referene point Y X for estimation.)r   metavartypenargsr   z
--ref-dateref_dateDATEz%Change reference date for estimation.)r   r%   r   zdates of interestz--start-datez-sr   z"start date for velocity estimationz
--end-datez-er   z end date for velocity estimationz	--excludez--exr   +zsdate(s) not included in velocity estimation, i.e.:
--exclude 20040502 20060708 20090103
--exclude exclude_date.txt
)r   r'   defaultr   Zbootstrappingz3estimating the mean / STD of the velocity estimatorz--bootstrapz--bootstrappingr   zLEnable bootstrapping to estimate the mean and STD of the velocity estimator.z--bcz--bootstrap-countr   i  z>number of iterations for bootstrapping (default: %(default)s).)r   r&   r+   r   zResidual filez4Save residual displacement time-series to HDF5 file.z
--save-resz--save_residualsave_resz8Save the residual displacement time-series to HDF5 file.z
--res-filez--residual-fileres_fileztimeseriesResidual.h5zJOutput file name for the residual time-series file (default: %(default)s).)r   r+   r   )argparseArgumentParserRawTextHelpFormatterTEMPLATE	REFERENCEEXAMPLEadd_argumentintadd_argument_groupDROP_DATE_TXTr   add_timefunc_argumentadd_memory_argument)parserdater   resid r=   P/home/centos/operations/rsmas_insar/sources/MintPy/mintpy/timeseries2velocity.pycreate_parserJ   sF    






r?   c             C   s   t  }|j| d}t|jd |_|jdkr>td|j|jrd|j	dkrdd|_t
d t
d |jrt
d	 |jdks|js|js|js|jrtd
|jrt|j|}t|}|S )zCommand line parser.)args	FILE_TYPE)r   r   r   z!input file is {}, NOT timeseries!   FzIbootstrap-count should be larger than 1, otherwise it does not make sensez/turn OFF bootstrapping and continue without it.zbootstrapping is turned ON.zMbootstrapping currently support polynomial ONLY and ONLY with the order of 1!)r?   
parse_argsr   read_attributer   key	Exceptionformatr   r   print
polynomialperiodicstepexplog
ValueErrorr   read_template2inpsinit_exp_log_dicts)iargsr:   inpsr=   r=   r>   cmd_line_parse   s"    
"rS   c             C   sN  t  | _| jrx| jD ]}|d |dd  }}t|dkrt|dkrdt|t | j|< qd	|}|d7 }|d7 }|d7 }t
|qt
d		|qW t  | _| jrJx| jD ]}|d |dd  }}t|dkr8t|dkrt|t | j|< n*d	|}|d
7 }|d7 }|d7 }t
|qt
d		|qW | S )zqInitialize the dictionaries of exp and log funcs
    By trarnslating inps.exp/log into inps.expDict/logDict.
    r   rB   N   z!NO characteristic time found: {}
z`one or more characteristic time(s) are required for each onset date for the exp function, e.g.:
z--exp 20181026 60 OR
zA--exp 20161231 80.5 200  # append as many char_times as you like!z.input onset time is NOT in YYYYMMDD format: {}z`one or more characteristic time(s) are required for each onset date for the log function, e.g.:
)dictexpDictrL   lennparrayastypefloattolistrG   rN   logDictrM   )rR   Zexp_listZ
onset_timeZ
char_timesmsgZlog_listr=   r=   r>   rP      s6    



rP   c                s  |s
t  }t|}tdtj|   t|jt	
d  fddt| D }x|D ]} |  }|dkr|||< |rb|dkrt|||< qb|dkr|dd	d
d	dd}t| ||< qb|dkrbt|||< qbW d}| kr| rt| |_|S )z.Read input template file into inps.excludeDatez!read options from template file: zmintpy.velocity.c                s    g | ]} |   kr|qS r=   )keys).0i)prefixtemplater=   r>   
<listcomp>   s    z&read_template2inps.<locals>.<listcomp>)r   )r   r   )r   [ ], )r   zmintpy.compute.maxMemory)rS   varsrH   ospathbasenamer   read_templater   utcheck_template_auto_valuelistr_   r	   yyyymmddreplacesplitr5   r[   	maxMemory)r   rR   iDictkeyListrE   valuer=   )rb   rc   r>   rO      s0    

rO   c                s   t d d}tjjs0d}t dj nZt dj tjj}tjj}||krzd}t dj nt dj |dkrt	j t
 fdd	tD rd}t d
t nt dt t d| |S )Nzupdate mode: ONskiprunz1) output file {} NOT found.z!1) output file {} already exists.z02) output file is NOT newer than input file: {}.z,2) output file is newer than input file: {}.c             3   s.   | ]&}t t|  t| d kV  qdS )NoneN)strrj   get
key_prefix)r`   rE   )atrrR   r=   r>   	<genexpr>  s    zrun_or_skip.<locals>.<genexpr>z93) NOT all key configuration parameters are the same: {}.z53) all key configuration parameters are the same: {}.zrun or skip: {}.)rH   rk   rl   isfiler   rG   getmtimer   r   rD   any
configKeys)rR   flagtitor=   )r   rR   r>   run_or_skip   s(    r   c             C   s,  t |}g }|t jt| j|d7 }|r:tdt|  | jrtd| j  t t | j}xFt	t
|D ]6}|| }|| |k rn||krntd|  || qnW | jrtd| j  t t | j}xFt	t
|D ]6}|| }|| |kr||krtd|  || qW ttt|}|S )N)date_list_allzexclude date:zstart date: z  remove date: z
end date: )r	   yyyymmdd2yearsread_date_listrq   r   rH   r|   r   rr   rangerW   appendr   sortedset)rR   dateListAllZyy_list_all
exDateListZyy_minra   r;   Zyy_maxr=   r=   r>   read_exclude_date  s.    
r   c       
         s@   j dkrt j}n* j dkr,t j}n j dkr@t j}|  t |j _t	j
 jdrt j }t j\}}tj|dddk}d|||d < t|dkrtd	t|   jt||  7  _ttt j _ fd
d|jD  _t j _ jd  _ jd  _td td|j|j td t jt|jkr|td ntd j j td tj fdd|jD tjd _  j!s<t	j
"t	j
 jd }d} j dkrt	j
 j#dd }|| }n |dkr.|#dd }	||	 }|d7 }| _! S )zGet inps.excludeDate full list
    Inputs:
        inps          - Namespace,
    Output:
        inps.excludeDate  - list of string for exclude date in YYYYMMDD format
    r   r   r   timeseriesRg)rB   r$   )axisr   REF_DATEz$number of empty dates to exclude: {}c                s   g | ]}| j kr|qS r=   )r   )r`   ra   )rR   r=   r>   rd   L  s    z"read_date_info.<locals>.<listcomp>z2--------------------------------------------------zdates from input file: {}
{}z)using all dates to calculate the velocityz*dates used to estimate the velocity: {}
{}c                s   g | ]}| j kqS r=   )r   )r`   ra   )rR   r=   r>   rd   Z  s    )dtyper   ZPARAMS)r   ZtimeseriesAzz.h5)$rE   r   r   r   r   openr   dateListr   rk   rl   rm   
startswithget_date_listr   readrX   nansumindexsumrH   rG   rY   r\   r   rq   r   rW   numDater   r   bool_dropDater   splitextrt   )
rR   Ztsobj	date_listdatar   r   fbaseoutnamerb   suffixr=   )rR   r>   read_date_info0  sR    




"

r   c             C   s  |s
| j }|d |d  }}t|}t|}| jrxH| jD ]>}t|}||  k r`|k s>n td| d| d| q>W | jrx8| j D ]*}t|}	|	|krtd| d| qW | jrx8| j D ]*}t|}	|	|krtd| d| qW t }
| j	|
d	< | j
|
d
< | j|
d< | j|
d< | j|
d< td x&|
 D ]\}}td|| qLW d	|
 krtd|
d	 d t|
d
 d  t|
d  tdd |
d  D  tdd |
d  D  }|
|fS )zget model info from inpsr   r   zinput step date "z" exceed date list min/max: z, zinput exp onset date "z" >= the last date: zinput log onset date "rI   rJ   rK   rL   rM   zEestimate deformation model with the following assumed time functions:z    {:<10} : {}z7linear/polynomial model is NOT included! Are you sure?!rB   r$   c             S   s   g | ]\}}t |qS r=   )rW   )r`   rE   valr=   r=   r>   rd     s    z#read_inps2model.<locals>.<listcomp>c             S   s   g | ]\}}t |qS r=   )rW   )r`   rE   r   r=   r=   r>   rd     s    )r   r	   r   rK   rN   rV   r_   r]   rU   rI   rJ   rH   itemsrG   rW   r   )rR   r   dmindmaxyminymaxZd_stepy_stepZd_onsetZy_onsetmodelrE   rx   	num_paramr=   r=   r>   read_inps2modell  sB    









\r   c       2   
   C   s  t | j}t|d t|d  }}| j}t| j}|dd}d|	 kr~| j
s~| jd | _
td td| jd  t| \}}t|}	d|	d	< d
|	d< | jd |	d< | jd |	d< d| jd | jd |	d< | jr| jd |	d< | jd |	d< | j
r| j
|	d< tdt x&tD ]}
tt| |
 |	t|
 < q"W t|||fddd  \}}tj| j|	||d | jrt|}x&dD ]}
|
|	 kr||
 qW tj| j|| jd ||d  d | | d }| jr|| j| | | d 7 }tt|d | jd  }tj dd||f|ddd}xht!|D ]Z\}}|d |d  }|d |d  }|| }|dkrtd |d | td!| td"| tj"||ft#d#}tj"||ft#d#}td$| j t j$| j|d%d }| j
rHtd&| j
 | j%| j
}|t&||d d d d f |j'd ddf8 }| jrtd'| jd | jd  | jd | jd | jd d | jd d f}t j$| j|d%d }|t&|(|j'd ddd|j'd |j'd f8 }|| j)d d d d f (| jd}|	d d(kr|d)9 }d }| j*rht j$| j*|d%d }|| j)d d d d f (| jd}d*}||||k < td+ tj+|dd,}t,t-| |d-k}~|d k	rtd. tj.t-|dd,} || dk9 }~ |d d |f }tt.|}!t/|d }"td/|!||!| d0  |!dkr q@td1 | jr4td2| j tj01 }#tj"| j||!ft#d#}$t2j3| jd3}%xvt4| jD ]h}|#j5| j| jdd4}&|&6  t7j8|||& 9 ||& |d5d |$|< |%j:|d d6|d | jd7 qzW |%;  |$j<dd,(|d|d d |f< |$j=dd,(|d|d d |f< ~$nvt7j8|| j||d5\}'|d d |f< }(|d k	rNtd8 t2j3|!d3}%xt4|!D ]}|"| })ydt>d9t?|d d |)f @  }*t>tAB|'jCD|*D|'EtjF}+tG|+|d d |)f< W n* tAjHk
r   tjI|d d |)f< Y nX |%j:|d d:d;|d |!d< qW |%;  n\td= tABtD|'jC|'},|((dd||  }+tGtDt>|,(dd|+|d d |f< |d |d |d |d g}-t||||d>d }.x2|.J D ]&\}/}0tjK| j|0(|||/|-d? qW | jr@d||d |d |d |d g}-tjL||| ftjFd#tjI }1|tD|'|d d |f  |1d d |f< tjK| j|1(|||d@|-d? q@W | jS )ANLENGTHWIDTHCENTER_LINE_UTCr   r   zHWARNING: No REF_DATE found in time-series file or input in command line.z#  Set "--ref-date {}" and continue.r   rA   zm/yearUNIT
START_DATEr   END_DATEz{}_{}DATE12REF_YrB   REF_Xz3add/update the following configuration metadata:
{})ds_shape)metadatads_name_dictds_unit_dict)r   )r   ref_filer$         i   @yT)box	num_split	dimension	print_msgz5
------- processing patch {} out of {} --------------zbox width:  {}zbox length: {})r   zreading data from file {} ...)r   zreferecing to date: {}z%referencing to point (y, x): ({}, {})mmgMbP?gh㈵>z3skip pixels with zero/nan value in all acquisitions)r   g        z1skip pxiels with nan STD value in any acquisitionz2number of pixels to invert: {} out of {} ({:.1f}%)d   z.estimating time functions via linalg.lstsq ...zEestimating time function STD with bootstrap resampling ({} times) ...)maxValue)sizers   )r   r   dis_tssecondsziteration {} / {})r   zDestimating time function STD from time-series STD pixel-by-pixel ...g      ?   z{}/{} pixels)everyr   zBestimating time function STD from time-series fitting residual ...)mask)r   datasetNameblockr   )Mr   rD   r   r5   r   rX   rY   r   r}   r_   r(   rH   rG   r   rU   r!   r   r|   rj   r~   model2hdf5_datasetr   layout_hdf5r   r,   popr-   r   r   ceilru   r   split_box2sub_boxes	enumeratezerosdataTyper   r   tileshapereshaper   r   nanmeanmultiplyisnanr   whererandomdefault_rngr	   progressBarr   choicesortr
   estimate_time_funcr\   updateclosemeanstddiagsquareflattenr   invTdotrZ   float32sqrtLinAlgErrornanr   write_hdf5_blockones)2rR   r   lengthwidthnum_datedatesr   r   r   ZatrVrE   r   r   ZatrRZ	memoryAllnum_boxbox_listra   r   box_widbox_len	num_pixelmm_stdts_dataref_indref_boxref_valts_stdepsilonZts_stackr   Znum_std_nannum_pixel2invidx_pixel2invrngZm_bootprog_barZboot_indGe2idxZC_ts_invZm_varZG_invr   ds_dictds_namer   ts_resr=   r=   r>   run_timeseries2time_func  s$   



.,2""



&  
$&&
,	
 (r  c                s  | d }t | d }t | d }tdd | d  D }i }	i }
i }|dk	r|jdkrb|jd nd}|d	|}|d	|}|dkrtj|tjd
}n|	 }xt
d|d D ]}|dkrd}d}n&|dkrd}d}nd|}d|}|dk	r"||ddf |	|< ||ddf |	|d < t|dg|
|< |||< t|dg|
|d < |||d < qW |d }x`t
|D ]R}| d |  ddg} dkrdd |D }n, dkrdd |D }n fdd|D }|dk	r||d|  ddf }||d|  d ddf }t|d |d  }tj|td
}t|| dks^t|| ||  ||< x$t|||gD ]\}}||	|< qnW t|dg|
|d < d||d < t|dg|
|d < d||d < qnW |d d|  }xt
|D ]}d| d | }|dk	r4||| ddf |	|< ||| ddf |	|d < t|dg|
|< d||< t|dg|
|d < d||d < qW |d d|  | }d}x| d  D ]}x| d | D ]}d||}|dk	r||| ddf |	|< ||| ddf |	|d < t|dg|
|< d||< t|dg|
|d < d||d < |d7 }qW qW |d d|  | | }d}x| d  D ]}x| d | D ]}d||}|dk	r||| ddf |	|< ||| ddf |	|d < t|dg|
|< d||< t|dg|
|d < d||d < |d7 }qxW qfW |	|
|fS ) a#  Prepare the estimated model parameters into a list of dicts for HDF5 dataset writing.
    Parameters: model        - dict,
                m            - 2D np.ndarray in (num_param, num_pixel) where num_pixel = 1 or length * width
                m_std        - 2D np.ndarray in (num_param, num_pixel) where num_pixel = 1 or length * width
                mask         - 1D np.ndarray in (num_pixel), mask of valid pixels
                ds_shape     - tuple of 2 int in (length, width)
    Returns:    ds_dict      - dict, dictionary of dataset values,     input for writefile.write_hdf5_block()
                ds_name_dict - dict, dictionary of dataset initiation, input for writefile.layout_hdf5()
                ds_unit_dict - dict, dictionary of dataset unit,       input for writefile.layout_hdf5()
    Examples:   # read input model parameters into dict
                model = read_inps2model(inps, date_list=inps.date_list)
                # for time series cube
                ds_name_dict, ds_name_dict = model2hdf5_dataset(model, ds_shape=(200,300))[1:]
                ds_dict = model2hdf5_dataset(model, m, m_std, mask=mask)[0]
                # for time series point
                ds_unit_dict = model2hdf5_dataset(model)[2]
                ds_dict = model2hdf5_dataset(model, m, m_std)[0]
    rI   rJ   rK   c             S   s   g | ]\}}t |qS r=   )rW   )r`   rE   r   r=   r=   r>   rd     s    z&model2hdf5_dataset.<locals>.<listcomp>rL   NrB   r   )r   r   zm/yearr$   accelerationzm/year^2zpoly{}z	m/year^{}StdZ	AmplitudePhasec             S   s   g | ]}d | qS )Zannualr=   )r`   xr=   r=   r>   rd     s    g      ?c             S   s   g | ]}d | qS )Z
semiAnnualr=   )r`   r  r=   r=   r>   rd     s    c                s   g | ]}d   | qS )ZperiodYr=   )r`   r  )periodr=   r>   rd     s    r   r  radianzstep{}z
exp{}Tau{}rM   z
log{}Tau{})rW   r   r   ndimr   r   rX   r   r   r   r   rG   r   r   r   allarctanzipr_   )r   r  r  r   r   poly_deg
num_periodnum_stepnum_expr  r   r   r  ra   dsNameunitp0ZdsNameSuffixesdsNamesZcoef_cosZcoef_sinZ
period_ampZ
period_phar   	exp_onsetexp_tau	log_onsetlog_taur=   )r  r>   r     s    








r   c             C   sd   t | }t }t|}|jr0t|dkr0|jS t| tt | d\}}td	|| |jS )Nry   <   z'time used: {:02.0f} mins {:02.1f} secs.)
rS   timer   r   r   r   r  divmodrH   rG   )rQ   rR   
start_timer  sr=   r=   r>   main7  s    r0  __main__rB   )N)N)N)NNNN)N),rk   sysr,  r.   numpyrX   scipyr   mintpy.defaults.templater   mintpy.objectsr   r   r   r   mintpy.utilsr   r	   r
   r   r   r   ro   r   r   r~   r   r1   r2   r3   r7   r?   rS   rP   rO   r   r   r   r   r  r   r0  __name__argvr=   r=   r=   r>   <module>   sD    =
+
 !!<
7 p
 &

