B
    QaJZ                 @   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 d dlmZmZmZmZ dZddd	ZdddZdd ZdddZdd Zdd ZdddZG dd dZdS )     N)Geod)urlretrieve)
coordinate)ptime	time_funcreadfileutils1z7http://geodesy.unr.edu/NGLStationPages/DataHoldings.txtTc             C   s0   t }tj|}| r"td| t|| |S )zdownload DataHoldings.txtz+downloading site list from UNR Geod Lab: {})unr_site_list_fileospathbasenameprintformatr   )	print_msgurlout_file r   H/home/centos/operations/rsmas_insar/sources/MintPy/mintpy/objects/gps.pydload_site_list   s    
r   2   c             C   s  |dkr*t jt}t j|s*t|d tj|tddd	t
}|dddf }|ddddf 	tjj\}}	|	t|	d d 8 }	td	d
 |dddf 	t
D }
tdd
 |dddf 	t
D }|dddf 	tj}|| d k|| d k |	| d k |	| d k }|rDt|gd d }|||k9 }|rjt|gd d }||
|k9 }|dk	r|||k9 }|| || |	| fS )a<  Search available GPS sites within the geo bounding box from UNR website
    Parameters: SNWE       : tuple of 4 float, indicating (South, North, West, East) in degrees
                start_date : string in YYYYMMDD format
                end_date   : string in YYYYMMDD format
                site_list_file : string.
                min_num_solution : int, minimum number of solutions available
    Returns:    site_names : 1D np.array of string for GPS station names
                site_lats  : 1D np.array for lat
                site_lons  : 1D np.array for lon
    N)r      )r   r                        	   
   )dtypeskiprowsusecolsr   r   g     v@c             S   s   g | ]}t j|d qS )z%Y-%m-%d)dtdatetimestrptime).0ir   r   r   
<listcomp><   s    zsearch_gps.<locals>.<listcomp>r   c             S   s   g | ]}t j|d qS )z%Y-%m-%d)r#   r$   r%   )r&   r'   r   r   r   r(   =   s    r   r   r   )r
   r   r   r	   isfiler   nploadtxtbytesastypestrfloat32Troundarrayint16r   date_list2vector)SNWE
start_dateend_datesite_list_fileZmin_num_solutionr   Ztxt_data
site_names	site_lats	site_lonsZt_startZt_endZnum_solutionidxt0t1r   r   r   
search_gps$   s0    
"&&0
r?   c             C   s   t ttt| t|@ }t j|jt jd}	xtt	|D ]}
t 
| ||
 kd d }t 
|||
 kd d }|| ||  d || ||  d  || ||  d  d }||	|
< q>W |	|	d 8 }	t j|	t jd}	||	fS )a  Calculate the baseline change between two GPS displacement time-series
    Parameters: dates1/2     : 1D np.array of datetime.datetime object
                pos_x/y/z1/2 : 1D np.ndarray of displacement in meters in np.float32
    Returns:    dates        : 1D np.array of datetime.datetime object for the common dates
                bases        : 1D np.ndarray of displacement in meters in np.float32 for the common dates
    )r    r   r   g      ?)r*   r2   sortedlistsetzerosshapefloat64rangelenwherer/   )Zdates1Zpos_x1Zpos_y1Zpos_z1Zdates2Zpos_x2Zpos_y2Zpos_z2datesbasesr'   idx1idx2Zbaseir   r   r   get_baseline_changeS   s    @rM   enu2los     VFc          	      sH  |rt ndd }t|}	tj| }
t| }|d }|dkrBdnd tj|
d| d}d	d
dddg}dgdgt|d   }d}tj|rt	j
||ddd}|j}|stj|r||	kr|d| t	j
||ddd}||   }nRg }|d tjddg|
dd}|r4|}|dtj| n|}|d tj|	|d}xt|D ]\}}|j|d d|d |	|d t|}|j|||||d \}}|jd!kr|d" |d  nt	j}||j|j|j||g qXW |  t	 fd#d$|D }|d%| t|d&$}t|}| | |!| W d'Q R X |S )(aS  Get the GPS LOS observations given the query info.

    Parameters: insar_file - str, InSAR LOS file, e.g. velocity or timeseries
                site_names - list of str, GPS sites, output of search_gps()
                start_date - str, date in YYYYMMDD format
                end_date   - str, date in YYYYMMDD format
                gps_comp   - str, flag of projecting 2/3D GPS into LOS
                             e.g. enu2los, hz2los, up2los
                horz_az_angle - float, azimuth angle of the horizontal motion in degree
                             measured from the north with anti-clockwise as positive
                print_msg  - bool, print verbose info
                redo       - bool, ignore existing CSV file and re-calculate
    Returns:    site_obs   - 1D np.ndarray(), GPS LOS velocity or displacement in m or m/yr
    c              _   s   d S )Nr   )argskwargsr   r   r   <lambda>z       z!get_gps_los_obs.<locals>.<lambda>	FILE_TYPE)velocityr   r   Zgps_z.csvZSiteZLonZLatZDisplacementZVelocityZU10f8r   r   ,T)r    	delimiternamesz#read GPS observations from file: {}zcalculating GPS observation ...incidenceAngleazimuthAnglegeo)work_dircoordz+use incidence / azimuth angle from file: {}z+use incidence / azimuth angle from metadata)maxValuer   z{}/{} {})suffix)r6   r7   gps_comphorz_az_angler   c                s   g | ]}|  qS r   r   )r&   x)obs_indr   r   r(      s    z#get_gps_los_obs.<locals>.<listcomp>z"write GPS observations to file: {}wN)"r   rG   r
   r   dirnamer   read_attributejoinr)   r*   
genfromtxtsizer   utget_geometry_filer   r   progressBar	enumerateupdateGPSget_gps_los_velocitynanappendsitesite_lonsite_latcloser2   opencsvwriterwriterow	writerows)
insar_filer9   r6   r7   ra   rb   r   redovprintZnum_sitefdirmetaZobs_typeZcsv_file	col_namesZ	col_typesnum_rowfcsite_obs	data_list	geom_filegeom_objprog_barr'   	site_nameobjveldis_tsdisZfcwr   )re   r   get_gps_los_obsi   sZ    
 " 

r   c             C   s   dd l }|j| dd}tj|dtdd}|d d df t}|d d df t}|d d df t}d	d
 t|||D }|d d df tj	 }|d d df tj	 }	|d d df tj	 }
|||	|
fS )Nr   cp1252)encoding   )*z-DATA)r!   r    commentsr   r   c             S   s"   g | ]\}}}t j|||d qS ))yearmonthday)r#   r$   )r&   ymdr   r   r   r(      s    z!read_pos_file.<locals>.<listcomp>r   r   r   )
codecsry   r*   r+   r.   r-   intziprE   tolist)fnamer   Zfcpr   ysmsdsrI   XYZr   r   r   read_pos_file   s    r   c             C   s6   t  tj| d|}dd |D }t|}|S )Nz{}.*.posc             S   s"   g | ]}t j|d d qS ).r   )r
   r   r   split)r&   r'   r   r   r   r(      s    z!get_pos_years.<locals>.<listcomp>)globr
   r   ri   r   r   yy2yyyy)gps_dirru   fnamesyearsr   r   r   get_pos_years   s    
r   c          
   C   s6  t |dd }t |dd }|| d }g g g g f\}}}	}
xjt|D ]^}t|| }tj| d||dd  }t|\}}}}||7 }||7 }|	|7 }	|
|7 }
qJW t	|}t	|}t	|	}	t	|
}
t
j|d}t
j|d}tj|jtjd}d|||k < d|||k< || || |	| |
| fS )	Nr   r   r   z	{}.{}.posr   z%Y%m%d)r    F)r   rF   r.   r
   r   ri   r   r   r*   r2   r#   r$   r%   onesrD   bool_)r   ru   r6   r7   Zyear0Zyear1Znum_yearrI   r   r   r   r'   Zyearir   ZdatesiXiZYiZiZdate0date1flagr   r   r   read_GSI_F3   s,    



r   c               @   sx   e Zd ZdZdddZdddZdd	d
Zd ddZd!ddZd"e	e	dddZ
d#ddZd$edddZd%ddZdS )&rq   a  GPS class for GPS time-series of daily solution

    Example:
      import matplotlib.pyplot as plt
      from mintpy.objects.gps import GPS
      from mintpy.utils import utils as ut
      gps_obj = GPS(site='GV05', data_dir='~/insarlab/GPS')
      gps_obj.open()
      dis_los = ut.enu2los(gps_obj.dis_e,
                           gps_obj.dis_n,
                           gps_obj.dis_u)
      dates = gps_obj.dates
      plt.figure()
      plt.scatter(dates, dis_los)
      plt.show()
    ./GPSIGS14c             C   s  || _ tj|| _|| _d| _|dkrDtj|dj||d| _	n0|dkrftj|dj|d| _	nt
d|d	}tj||tj| j	| _tj|d
|| _d}|dkr|d|7 }n|dkr|d|7 }tj|d|| _tjtj|d| _tj| js t  tj| jtdddt}||krTt
d|tx4|tj| jgD ]}tj|sht| qhW d S )NzNevada Geodetic LabZIGS08z{s}.{v}.tenv3)svr   z	{s}.tenv3)r   z"un-recognized GPS data version: {}z+http://geodesy.unr.edu/gps_timeseries/tenv3z
pic/{}.pngzhttp://geodesy.unr.edu/tsplotsz/{0}z/{0}/{0}zTimeSeries/{}.pngzDataHoldings.txtr   r   )r    r!   r"   zSite {} NOT found in file: {})ru   r
   r   abspathdata_dirversionsourceri   r   file
ValueErrorr   file_url	plot_fileplot_file_urlrg   r8   r)   r   r*   r+   r,   r-   r.   r	   isdirmakedirs)selfru   r   r   Z
url_prefixr9   r   r   r   r   __init__  s6    
zGPS.__init__Tc             C   s2   t j| js|   | j|d | j|d d S )N)r   )r
   r   r)   r   
dload_siteget_stat_lat_lonread_displacement)r   r   r   r   r   ry   >  s    zGPS.openc             C   s:   |rt d| j| j t| j| j t| j| j | jS )Nzdownloading {} from {})r   r   ru   r   r   r   r   r   )r   r   r   r   r   r   D  s
    zGPS.dload_sitec             C   s   |rt d tj| js&| j|d tj| jtdd	t
}t|d d }}|ddd	f 	tj\}}}}||7 }||7 }t||tj d
 }	t|d |d  }
tdd}||||	|
dd \| _| _| j| jfS )zGet station lat/lonzcalculating station lat/lon)r   r   )r    r!   )r   r   g        r   r      g     f@r   WGS84)ellps)r   r
   r   r)   r   r   r*   r+   r,   r-   r.   floatarctan2pisqrtr   fwdrv   rw   )r   r   dataref_lonref_late0Ze_offn0Zn_offazdistgr   r   r   r   M  s     
 zGPS.get_stat_lat_lonNFc             C   s  t j| js| j|d |r&td tj| jtdd	t
}tdd |dddf D | _|dddf 	tjj\| _| _| _| _| _| _tt| jtj}|rt|gd	 d	 }d	|| j|k < |rt|gd	 d	 }d	|| j|k< | j| | _| j| | _| j| | _| j| | _| j| | _| j| | _| j| | _|rd	dlm}	 |	jd
ddd\}
}|d	 j| j| jddd |d j| j| jddd |d j| j| jddd |	  | j| j| j| j| j| j| jfS )ay   Read GPS displacement time-series (defined by start/end_date)
        Parameters: start/end_date : str in YYYYMMDD format
        Returns:    dates : 1D np.ndarray of datetime.datetime object
                    dis_e/n/u : 1D np.ndarray of displacement in meters in np.float32
                    std_e/n/u : 1D np.ndarray of displacement STD in meters in np.float32
        )r   z>reading time and displacement in east/north/vertical directionr   )r    r!   c             S   s   g | ]}t j|d qS )z%y%b%d)r#   r$   r%   )r&   r'   r   r   r   r(   p  s    z)GPS.read_displacement.<locals>.<listcomp>N)r   r               r   r   T)nrowsncolssharexr   ZEast)r   labelZNorthr   ZUp)r
   r   r)   r   r   r   r*   r+   r,   r-   r.   r2   rI   r/   r0   dis_edis_ndis_ustd_estd_nstd_ur   rG   r   r   r4   matplotlib.pyplotpyplotsubplotsscattershow)r   r6   r7   r   displayr   Zt_flagr=   r>   pltfigaxr   r   r   r   `  s>    "2zGPS.read_displacementrN        V)	inc_angleaz_anglec             C   st  |t jd 9 }|t jd 9 }|t jd 9 }t |t | d t |t | t |g}| }|dkrpn|dkrd|d< n~|dkrd|d< d|d	< nd|d
krt |d |d< t ||d	< d|d< n2|dkrd|d< d|d	< d|d< ntdt| | j|d  | j|d	   | j	|d   | _
| j|d  d | j|d	  d  | j|d  d  d | _| j
| jfS )a  Convert displacement in ENU to LOS direction
        Parameters: inc_angle : float, LOS incidence angle in degree
                    az_angle  : float, LOS aziuth angle in degree
                        from the north, defined as positive in clock-wise direction
                    gps_comp  : string, GPS components used to convert to LOS direction
                    horz_az_angle : float, fault azimuth angle used to convert horizontal to fault-parallel
        Returns:    dis_los  : 1D np.array for displacement in LOS direction
                    std_los  : 1D np.array for displacement standard deviation in LOS direction
        g     f@rc   )rN   )Zen2loshz2losg        r   )Zu2losup2losr   r   )horz)verticalg      ?zUn-known input gps components:g      ?)r*   r   sincoslowerr   r.   r   r   r   dis_losr   r   r   std_los)r   r   r   ra   rb   Zunit_vecr   r   r   displacement_enu2los  s4    




,<zGPS.displacement_enu2losc             C   s  | j |d\}}t|trt|}t||d}|j|||ddd \}}td|}tt	|d d |}td|}tt	|d d |}|||d |d f}	tj
|d|	|d	d d
 }
tj
|d|	|d	d d
 }n>t|trtj|d|d}
tt|d }ntd||
|fS )zNGet the Line-of-Sight geometry info in incidence and azimuth angle in degrees.)r   )lookup_filer   r   LENGTHr   WIDTHrZ   )datasetNameboxr   )r   r   r[   )	dimensionr   HEADINGz)input geom_obj is neight str nor dict: {})r   
isinstancer.   r   rh   r   	geo2radarmaxminr   readdictrl   incidence_angleheading2azimuth_angler   r   r   )r   r   r   latlonatrr^   r   rd   r  r   r   r   r   r   get_los_geometry  s"    


 
 zGPS.get_los_geometry)ra   c             C   st  |  |\}}	| j|||dd }
| j||	||d\}}| j|d}|rbt|| jd}|j|||d | |\}}	|j||	||d |j|d}ttt	t
| jt
|j@ }
t|
jtj}t|
jtj}xtt|
D ]x}t| j|
| kd d }t|j|
| kd d }| j| |j|  ||< | j| d |j| d  d ||< qW nd}|
||||fS )aj  Read GPS displacement in LOS direction
        Parameters: geom_obj : dict / str, metadata of InSAR file, or geometry file path
                    start_date : string in YYYYMMDD format
                    end_date   : string in YYYYMMDD format
                    ref_site   : string, reference GPS site
                    gps_comp   : string, GPS components used to convert to LOS direction
                    az_angle   : float, fault azimuth angle used to convert horizontal to fault-parallel
        Returns:    dates : 1D np.array of datetime.datetime object
                    dis   : 1D np.array of displacement in meters
                    std   : 1D np.array of displacement uncertainty in meters
                    site_lalo : tuple of 2 float, lat/lon of GPS site
                    ref_site_lalo : tuple of 2 float, lat/lon of reference GPS site
        )r   r   )ra   rb   )ru   r   r   g      ?N)r  r   r   r   rq   r   r*   r2   r@   rA   rB   rI   rC   rD   r/   rF   rG   rH   r   r   )r   r   r6   r7   ref_sitera   rb   r   r   r   rI   r   stdZ	site_laloZref_objref_site_lalor'   rK   rL   r   r   r   read_gps_los_displacement  s(    "*zGPS.read_gps_los_displacementc             C   st   | j ||||||dd d \}}dd |D }	t|	dkrbt|	}
ttj|
|d | _ntj	| _| j|fS )N)r6   r7   r  ra   rb   r   c             S   s   g | ]}t j|d qS )z%Y%m%d)r#   r$   strftime)r&   r'   r   r   r   r(     s    z,GPS.get_gps_los_velocity.<locals>.<listcomp>r   )
r  rG   r   get_design_matrix4time_funcr*   dotlinalgpinvrU   rs   )r   r   r6   r7   r  ra   rb   rI   r   	date_listAr   r   r   rr     s    
zGPS.get_gps_los_velocity)r   r   )T)T)T)NNTF)rN   r   )F)NNNrN   r   F)NNNrN   r   )__name__
__module____qualname____doc__r   ry   r   r   r   r   r   r  r.   r  rr   r   r   r   r   rq      s   
,

	

82
 + rq   )T)NNNr   T)rN   rO   TF)NN)r
   rz   r   r$   r#   numpyr*   pyprojr   urllib.requestr   mintpy.objects.coordr   mintpy.utilsr   r   r   r   rl   r	   r   r?   rM   r   r   r   r   rq   r   r   r   r   <module>   s$   


/ 
_
 