U
     xbEc                     @   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   s4   t }| s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   )out_file	print_msgurl r   I/home/exouser/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 }|r:t|gd d nd}|rXt|gd d nd}|rn|||k9 }|r||
|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 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_namesZ	site_latsZ	site_lonsZt0sZt1sZnum_solutionidxt0t1r   r   r   
search_gps%   s>    "&&


rB   c                 C   s   t ttt| t|@ }t j|jt jd}	tt	|D ]}
t 
| ||
 kd d }t 
|||
 kd d }|| ||  d || ||  d  || ||  d  d }||	|
< q<|	|	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         ?)r/   r7   sortedlistsetzerosshapefloat64rangelenwherer4   )Zdates1Zpos_x1Zpos_y1Zpos_z1Zdates2Zpos_x2Zpos_y2Zpos_z2datesbasesr,   idx1idx2Zbaseir   r   r   get_baseline_changeU   s     
rR   enu2los     VFc	               	      s  |rt ndd }	t|}
ddg}||kr6td| | dkrFdnd tj| d }tj|d	| }||d
kr|dnd7 }|d7 }dddddg}dgdgt|d   }|	d|  d}tj|rt	j
||ddd}|j}|stj|r||
kr|	d| t	j
||ddd}||   }||
kr||d  }t	j|t	jd}t	j|
t	jdt	j }t|D ](\}}||krx|||k d ||< qxnNg }|	d tjddg|d d!}|r|}|	d"tj| n| }|	d# tj|
|d$}t|D ]\}}|j|d d%|d |
|d& t|}|j|||||d'\}}t	|rdt	jn|d( |d  }||j|j|j||g q
|  |	d)| t |d*$}t!"|}|#| |$| W 5 Q R X t	 fd+d,|D }|S )-a  Get the GPS LOS observations given the query info.

    Parameters: meta       - dict, dictionary of metadata of the InSAR file
                obs_type   - str, GPS observation data type, displacement or velocity.
                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
    Examples:   from mintpy.objects import gps
                from mintpy.utils import readfile, utils as ut
                meta = readfile.read_attribute('geo/geo_velocity.h5')
                SNWE = ut.four_corners(meta)
                site_names = gps.search_gps(SNWE, start_date='20150101', end_date='20190619')
                vel = gps.get_gps_los_obs(meta, 'velocity',     site_names, start_date='20150101', end_date='20190619')
                dis = gps.get_gps_los_obs(meta, 'displacement', site_names, start_date='20150101', end_date='20190619')
    c                  _   s   d S )Nr   )argskwargsr   r   r   <lambda>       z!get_gps_los_obs.<locals>.<lambda>displacementvelocityzun-supported obs_type: r   r   	FILE_PATHZgps_horzz.0f z.csvZSiteZLonZLatZDisplacementZVelocityZU10f8r   z#default GPS observation file name: r   ,T)r"   	delimiternamesz#read GPS observations from file: {}rC   zcalculating GPS observation ...incidenceAngleazimuthAnglegeo)Zwork_dircoordz+use incidence / azimuth angle from file: {}z+use incidence / azimuth angle from metadata)maxValuer   z{}/{} {})suffix)r;   r<   gps_comphorz_az_anglez"write GPS observations to file: {}wc                    s   g | ]}|  qS r   r   r+   xZobs_indr   r   r-      s     z#get_gps_los_obs.<locals>.<listcomp>)%r   rL   
ValueErrorlowerr
   r   dirnamejoinr.   r/   
genfromtxtsizer   r7   r4   rH   nan	enumerateutZget_geometry_filer   r   progressBarupdateGPSget_gps_los_velocityisnanappendsitesite_lonsite_latcloseopencsvwriterZwriterowZ	writerows) metaZobs_typer>   r;   r<   rh   ri   r   ZredovprintZnum_siteZ	obs_typesZfile_dirZcsv_fileZ	col_namesZ	col_typesnum_rowfcZsite_obsZ
temp_namesZtemp_obsr,   Z	site_nameZ	data_list	geom_filegeom_objprog_barobjZveldis_tsdisZfcwr   rn   r   get_gps_los_obsk   sr    

 
"	

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   )
codecsr   r/   r0   r3   r2   intziprJ   tolist)fnamer   Zfcpr   ysmsdsrN   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   splitr*   r   r   r   r-      s     z!get_pos_years.<locals>.<listcomp>)globr
   r   rr   r   r   yy2yyyy)gps_dirr~   fnamesyearsr   r   r   get_pos_years   s    
r   c              
   C   s2  t |dd }t |dd }|| d }g g g g f\}}}	}
t|D ]^}t|| }tj| d||dd  }t|\}}}}||7 }||7 }|	|7 }	|
|7 }
qHt	|}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   %Y%m%drC   F)r   rK   r3   r
   r   rr   r   r   r/   r7   r'   r(   r)   onesrI   bool_)r   r~   r;   r<   Zyear0Zyear1Znum_yearrN   r   r   r   r,   Zyearir   ZdatesiXiZYiZ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 )&rz   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|tj| jfD ]}tj|sft| qfd 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!   zSite {} NOT found in file: {})r~   r
   r   abspathdata_dirversionsourcerr   r   filero   r   file_url	plot_fileplot_file_urlrq   r=   r.   r   r/   r0   r1   r2   r3   r	   isdirmakedirs)selfr~   r   r   Z
url_prefixr>   Zfdirr   r   r   __init__0  s6    
zGPS.__init__Tc                 C   s2   t j| js|   | j|d | j|d d S )Nr   )r
   r   r.   r   
dload_siteget_stat_lat_lonread_displacementr   r   r   r   r   r   \  s    zGPS.openc                 C   s:   |rt d| j| j t| j| j t| j| j | jS )Nzdownloading {} from {})r   r   r~   r   r   r   r   r   r   r   r   r   r   b  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/lonr   r   r"   r#   )r   r           r   r           f@r   WGS84)ellps)r   r
   r   r.   r   r   r/   r0   r1   r2   r3   floatarctan2pisqrtr   fwdr   r   )r   r   dataZref_lonZref_latZe0Ze_offZn0Zn_offazdistgr   r   r   r   k  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 | jD | _|ddd	f 	tjj\| _| _| _| _| _| _tt| jtj}|rt|gd
 d
 }d
|| j|k < |rt|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   c                 S   s   g | ]}t j|d qS )z%y%b%dr&   r*   r   r   r   r-     s     z)GPS.read_displacement.<locals>.<listcomp>Nc                 S   s   g | ]}| d qS r   )strftimerl   r   r   r   r-     s     )r   r                r   r   T)nrowsncolssharexr   ZEast)r   labelZNorthr   ZUp) r
   r   r.   r   r   r   r/   r0   r1   r2   r3   r7   rN   	date_listr4   r5   dis_edis_ndis_ustd_estd_nstd_ur   rL   r   r   r9   matplotlib.pyplotpyplotsubplotsscattershow)r   r;   r<   r   displayr   Zt_flagr@   rA   pltfigaxr   r   r   r   ~  sV    "    zGPS.read_displacementrS   rT   )	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
        r   rj   )rS   )Zen2losZhz2losr   r   )Zu2losZup2losr   r   )r\   
horizontal)vertverticalg      ?zUn-known input gps components:rD   )r/   r   sincosrp   ro   r3   r   r   r   dis_losr   r   r   std_los)r   r   r   rh   ri   Zunit_vecr   r   r   displacement_enu2los  sJ    




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   WIDTHrb   )datasetNameboxr   )r   r   rc   )	dimensionr   ZHEADINGz)input geom_obj is neight str nor dict: {})r   
isinstancer3   r   read_attributer   	geo2radarmaxminr   readdictrw   Zincidence_angleZheading2azimuth_angler   ro   r   )r   r   r   latlonatrre   r   rm   r  r   r   r   r   r   get_los_geometry  s"    


 
 zGPS.get_los_geometry)rh   c                 C   sp  |  |\}}	| j|||dd }
| j||	||d\}}| j|d}|r^t|| jd}|j|||d | |\}}	|j||	||d |j|d}ttt	t
| jt
|j@ }
t|
jtj}t|
jtj}tt|
D ]x}t| j|
| kd d }t|j|
| kd d }| j| |j|  ||< | j| d |j| d  d ||< qn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   )rh   ri   )r~   r   r   rD   N)r  r   r   r   rz   r   r/   r7   rE   rF   rG   rN   rH   rI   r4   rK   rL   rM   r   r   )r   r   r;   r<   ref_siterh   ri   r   r   r   rN   r   stdZ	site_laloZref_objref_site_lalor,   rP   rQ   r   r   r   read_gps_los_displacement  s(    "(zGPS.read_gps_los_displacementc                 C   s   | j ||||||dd d \}}d}	t|dkr8d}	nP|r|rt|gd d }
t|gd d }|d |d  ||
 d krd}	|	rdd	 |D }t|}ttj	||d
 | _
ntj| _
| j
|fS )N)r;   r<   r  rh   ri   r   TFr   rj   r   c                 S   s   g | ]}t j|d qS r   )r'   r(   r   r*   r   r   r   r-   J  s     z,GPS.get_gps_los_velocity.<locals>.<listcomp>r   )r  rL   r   r9   r   get_design_matrix4time_funcr/   dotlinalgpinvrZ   ru   )r   r   r;   r<   r  rh   ri   rN   r   Zdis2velr@   rA   r   Ar   r   r   r{   1  s0    

zGPS.get_gps_los_velocity)r   r   )T)T)T)NNTF)rS   rT   )F)NNNrS   rT   F)NNNrS   rT   )__name__
__module____qualname____doc__r   r   r   r   r   r   r   r  r3   r  r{   r   r   r   r   rz     s&   
,

	

92
      ,    rz   )NT)NNNr   T)rS   rT   TF)NN)r
   r   r   r(   r'   numpyr/   pyprojr   urllib.requestr   Zmintpy.objects.coordr   mintpy.utilsr   r   r   r   rw   r	   r   rB   rR   r   r   r   r   rz   r   r   r   r   <module>   s*   

0      
|
 