Python package for creating GNSS STD/SWD as well as modelled (GPT/VMF) ZTD
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""    


import numpy as np
import pandas as pd
from GPT3 import gpt3_5_fast
from VMF import read_vmf1_grid,vmf1_ht
from GMF_GradinetMF import gmf,gradient_MF
# File created by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================

def ZHD_Saastamoinen(mjd = None,h_ell = None,h_ortho=None,lat = None,lon = None):
    
#     Input parameters:
#     -------------------    
#     mjd:   modified Julian date (scalar, only one epoch per call is possible)
#     lat:   ellipsoidal latitude in radians [-pi/2:+pi/2] (vector)
#     lon:   longitude in radians [-pi:pi] or [0:2pi] (vector)
#     h_ell: ellipsoidal height in m (vector)
#     it:    case 1: no time variation but static quantities
#     case 0: with time variation (annual and semiannual terms)

    P,T,dT,Tm,e,ah,aw,la,undu,Gn_h,Ge_h,Gn_w,Ge_w= gpt3_5_fast(mjd,lat,lon,h_ell,it =0)
    ZHD=pd.DataFrame(np.zeros((len(P),len(P))))
    
    ZHD = (0.0022768*P) / (1 - (0.00266*np.cos(2*lat)) -(0.00000028*h_ortho))
    
    return ZHD
###############################################################################
# def ZHD_VMF(mjd = None,ell = None):
    
# #     Input parameters:
# #     -------------------     
# #     mjd   ... Modified Julian Date
# #     ell   ... ellipsoidal station coordinates lat [degrees], lon [degrees] ,hell [m] 'Should be in list form e.g 'ell=[42.72,56.13,570]'' 

#     ah,aw,ZHD,ZWD=read_vmf1_grid(mjd,ell)
    
#     return ZHD
    
###############################################################################
#ZHD calculated from GPT model and Saastamoinen Function zhd_GPT
#ZHD derived from VMF1 grid zhd_VMF
def CalcZWD(mjd = None,ZTD = None,h_ell= None,h_ortho=None,lat = None,lon = None,Type_zhd=None,zhd=None):
    
#     Input parameters:
#     -------------------    
#     mjd:   modified Julian date (scalar, only one epoch per call is possible)
#     lat:   ellipsoidal latitude in radians [-pi/2:+pi/2] (vector)
#     lon:   longitude in radians [-pi:pi] or [0:2pi] (vector)
#     h_ell: ellipsoidal height in m (vector)
#     ZTD: Zenith tropospheric delay
#     Type_zhd: 
#     zhd_GPT: ZHD calculated from GPT model and Saastamoinen 
#     zhd_VMF: ZHD derived from VMF1 grid

    if (Type_zhd=='zhd_GPT'):
        ZHD = ZHD_Saastamoinen(mjd,h_ell,h_ortho,lat,lon)
    else:
        ZHD=zhd
        
    
    ZWD = ZTD - ZHD
    return ZWD
###############################################################################
def CalcSWDgrad(mjd = None,ZTD = None,dlat = None,dlon = None,dhgt = None,h_ortho=None,gn = None,ge = None,AZ = None,ELV = None,Type_MF=None,Type_zhd=None,grad=None): 

    #     Input parameters:
#     -------------------    
#     mjd:   modified Julian date (scalar, only one epoch per call is possible)
#     lat:   ellipsoidal latitude in radians [-pi/2:+pi/2] (vector)
#     lon:   longitude in radians [-pi:pi] or [0:2pi] (vector)
#     dhgt: ellipsoidal height in m (vector)
#     ZWD: Zenith wet delay
#     gn & ge: Horizontal gradients
#     az & ELV: Azimuth and elevation angles
#     Type_MF:
#     GMF or VMF mapping function
#     Type_zhd: 
#     zhd_GPT: ZHD calculated from GPT model and Saastamoinen 
#     zhd_VMF: ZHD derived from VMF1 grid
    zhd=0
    if (Type_zhd=='zhd_VMF' or Type_MF=='VMF'):
        ell=[np.degrees(dlat),np.degrees(dlon),dhgt]
        ah,aw,zhd,zwd=read_vmf1_grid(mjd,ell)
        

    ZWD=CalcZWD(mjd ,ZTD,dhgt,h_ortho,dlat,dlon,Type_zhd,zhd)


    zd = np.radians(90 - ELV)    
    
    if (Type_MF=='GMF'):
       
       mf_h,mf_w = gmf(mjd,dlat,dlon,dhgt,zd)
    else:
       mf_h,mf_w = vmf1_ht(mjd,dlat,dlon,dhgt,zd,ah,aw)
    
    
    MF_G=gradient_MF(ELV)
    Grad_part=np.multiply(gn,np.cos(np.radians(AZ)))+np.multiply(ge,np.sin(np.radians(AZ)))
    if grad=='True':
       SWD=ZWD*np.array(mf_w)+np.multiply(np.array(MF_G),np.array(Grad_part))
    else:
       SWD=ZWD*np.array(mf_w)
    return SWD,ZWD,mf_w,MF_G
 ###############################################################################
 def CalcSTD(mjd = None,dlat = None,dlon = None,dhgt = None,h_ortho=None,AZ = None,ELV = None,Type_MF=None,Type_zhd=None): 

    #     Input parameters:
#     -------------------    
#     mjd:   modified Julian date (scalar, only one epoch per call is possible)
#     lat:   ellipsoidal latitude in radians [-pi/2:+pi/2] (vector)
#     lon:   longitude in radians [-pi:pi] or [0:2pi] (vector)
#     dhgt: ellipsoidal height in m (vector)
#     ZWD: Zenith wet delay
#     gn & ge: Horizontal gradients
#     az & ELV: Azimuth and elevation angles
#     Type_MF:
#     GMF or VMF mapping function
#     Type_zhd: 
#     zhd_GPT: ZHD calculated from GPT model and Saastamoinen 
#     zhd_VMF: ZHD derived from VMF1 grid
    zhd=0
    if (Type_zhd=='zhd_VMF' or Type_MF=='VMF'):
        ell=[np.degrees(dlat),np.degrees(dlon),dhgt]
        ah,aw,zhd,zwd=read_vmf1_grid(mjd,ell)
        



    zd = np.radians(90 - ELV)    
    
    if (Type_MF=='GMF'):
       
       mf_h,mf_w = gmf(mjd,dlat,dlon,dhgt,zd)
    else:
       mf_h,mf_w = vmf1_ht(mjd,dlat,dlon,dhgt,zd,ah,aw)
    
    

    STD=zwd*np.array(mf_w)+zhd*np.array(mf_h)
    ZTD=zhd+zwd
    return STD,ZTD
CalcSwd.py
zohrehadavi/STD_SWD_Calc
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""

import numpy as np
import pandas as pd   
def gmf(dmjd = None,dlat = None,dlon = None,dhgt = None,zd = None): 
    # This subroutine determines the Global Mapping Functions GMF
# Reference: Boehm, J., A.E. Niell, P. Tregoning, H. Schuh (2006),
# Global Mapping Functions (GMF): A new empirical mapping function based on numerical weather model data,
# Geoph. Res. Letters, Vol. 33, L07304, doi:10.1029/2005GL025545.
    
    # input data
# ----------
# dmjd: modified julian date
# dlat: ellipsoidal latitude in radians
# dlon: longitude in radians
# dhgt: height in m
# zd:   zenith distance in radians
    
# output data
# -----------
# gmfh: hydrostatic mapping function
# gmfw: wet mapping function
    
# Johannes Boehm, 2005 August 30
    
# ref 2006 Aug. 14: recursions for Legendre polynomials (O. Montenbruck)
# ref 2011 Jul. 21: latitude -> ellipsoidal latitude (J. Boehm)

#    File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================  
    
    ah_mean = np.array([+ 125.17,+ 0.8503,+ 0.06936,- 6.76,+ 0.1771,+ 0.0113,+ 0.5963,+ 0.01808,+ 0.002801,- 0.001414,- 1.212,+ 0.093,+ 0.003683,+ 0.001095,+ 4.671e-05,+ 0.3959,- 0.03867,+ 0.005413,- 0.0005289,+ 0.0003229,+ 2.067e-05,+ 0.3,+ 0.02031,+ 0.0059,+ 0.0004573,- 7.619e-05,+ 2.327e-06,+ 3.845e-06,+ 0.1182,+ 0.01158,+ 0.005445,+ 6.219e-05,+ 4.204e-06,- 2.093e-06,+ 1.54e-07,- 4.28e-08,- 0.4751,- 0.0349,+ 0.001758,+ 0.0004019,- 2.799e-06,- 1.287e-06,+ 5.468e-07,+ 7.58e-08,- 6.3e-09,- 0.116,+ 0.008301,+ 0.0008771,+ 9.955e-05,- 1.718e-06,- 2.012e-06,+ 1.17e-08,+ 1.79e-08,- 1.3e-09,+ 1e-10])
    bh_mean = np.array([+ 0.0,+ 0.0,+ 0.03249,+ 0.0,+ 0.03324,+ 0.0185,+ 0.0,- 0.1115,+ 0.02519,+ 0.004923,+ 0.0,+ 0.02737,+ 0.01595,- 0.0007332,+ 0.0001933,+ 0.0,- 0.04796,+ 0.006381,- 0.0001599,- 0.0003685,+ 1.815e-05,+ 0.0,+ 0.07033,+ 0.002426,- 0.001111,- 0.0001357,- 7.828e-06,+ 2.547e-06,+ 0.0,+ 0.005779,+ 0.003133,- 0.0005312,- 2.028e-05,+ 2.323e-07,- 9.1e-08,- 1.65e-08,+ 0.0,+ 0.03688,- 0.0008638,- 8.514e-05,- 2.828e-05,+ 5.403e-07,+ 4.39e-07,+ 1.35e-08,+ 1.8e-09,+ 0.0,- 0.02736,- 0.0002977,+ 8.113e-05,+ 2.329e-07,+ 8.451e-07,+ 4.49e-08,- 8.1e-09,- 1.5e-09,+ 2e-10])
    ah_amp = np.array([- 0.2738,- 2.837,+ 0.01298,- 0.3588,+ 0.02413,+ 0.03427,- 0.7624,+ 0.07272,+ 0.0216,- 0.003385,+ 0.4424,+ 0.03722,+ 0.02195,- 0.001503,+ 0.0002426,+ 0.3013,+ 0.05762,+ 0.01019,- 0.0004476,+ 6.79e-05,+ 3.227e-05,+ 0.3123,- 0.03535,+ 0.00484,+ 3.025e-06,- 4.363e-05,+ 2.854e-07,- 1.286e-06,- 0.6725,- 0.0373,+ 0.0008964,+ 0.0001399,- 3.99e-06,+ 7.431e-06,- 2.796e-07,- 1.601e-07,+ 0.04068,- 0.01352,+ 0.0007282,+ 9.594e-05,+ 2.07e-06,- 9.62e-08,- 2.742e-07,- 6.37e-08,- 6.3e-09,+ 0.08625,- 0.005971,+ 0.0004705,+ 2.335e-05,+ 4.226e-06,+ 2.475e-07,- 8.85e-08,- 3.6e-08,- 2.9e-09,+ 0.0])
    bh_amp = np.array([+ 0.0,+ 0.0,- 0.1136,+ 0.0,- 0.1868,- 0.01399,+ 0.0,- 0.1043,+ 0.01175,- 0.00224,+ 0.0,- 0.03222,+ 0.01333,- 0.002647,- 2.316e-05,+ 0.0,+ 0.05339,+ 0.01107,- 0.003116,- 0.0001079,- 1.299e-05,+ 0.0,+ 0.004861,+ 0.008891,- 0.0006448,- 1.279e-05,+ 6.358e-06,- 1.417e-07,+ 0.0,+ 0.03041,+ 0.00115,- 0.0008743,- 2.781e-05,+ 6.367e-07,- 1.14e-08,- 4.2e-08,+ 0.0,- 0.02982,- 0.003,+ 1.394e-05,- 3.29e-05,- 1.705e-07,+ 7.44e-08,+ 2.72e-08,- 6.6e-09,+ 0.0,+ 0.01236,- 0.0009981,- 3.792e-05,- 1.355e-05,+ 1.162e-06,- 1.789e-07,+ 1.47e-08,- 2.4e-09,- 4e-10])
    aw_mean = np.array([+ 56.4,+ 1.555,- 1.011,- 3.975,+ 0.03171,+ 0.1065,+ 0.6175,+ 0.1376,+ 0.04229,+ 0.003028,+ 1.688,- 0.1692,+ 0.05478,+ 0.02473,+ 0.0006059,+ 2.278,+ 0.006614,- 0.0003505,- 0.006697,+ 0.0008402,+ 0.0007033,- 3.236,+ 0.2184,- 0.04611,- 0.01613,- 0.001604,+ 5.42e-05,+ 7.922e-05,- 0.2711,- 0.4406,- 0.03376,- 0.002801,- 0.000409,- 2.056e-05,+ 6.894e-06,+ 2.317e-06,+ 1.941,- 0.2562,+ 0.01598,+ 0.005449,+ 0.0003544,+ 1.148e-05,+ 7.503e-06,- 5.667e-07,- 3.66e-08,+ 0.8683,- 0.05931,- 0.001864,- 0.0001277,+ 0.0002029,+ 1.269e-05,+ 1.629e-06,+ 9.66e-08,- 1.015e-07,- 5e-10])
    bw_mean = np.array([+ 0.0,+ 0.0,+ 0.2592,+ 0.0,+ 0.02974,- 0.5471,+ 0.0,- 0.5926,- 0.103,- 0.01567,+ 0.0,+ 0.171,+ 0.09025,+ 0.02689,+ 0.002243,+ 0.0,+ 0.3439,+ 0.02402,+ 0.00541,+ 0.001601,+ 9.669e-05,+ 0.0,+ 0.09502,- 0.03063,- 0.001055,- 0.0001067,- 0.000113,+ 2.124e-05,+ 0.0,- 0.3129,+ 0.008463,+ 0.0002253,+ 7.413e-05,- 9.376e-05,- 1.606e-06,+ 2.06e-06,+ 0.0,+ 0.2739,+ 0.001167,- 2.246e-05,- 0.0001287,- 2.438e-05,- 7.561e-07,+ 1.158e-06,+ 4.95e-08,+ 0.0,- 0.1344,+ 0.005342,+ 0.0003775,- 6.756e-05,- 1.686e-06,- 1.184e-06,+ 2.768e-07,+ 2.73e-08,+ 5.7e-09])
    aw_amp = np.array([+ 0.1023,- 2.695,+ 0.3417,- 0.1405,+ 0.3175,+ 0.2116,+ 3.536,- 0.1505,- 0.0166,+ 0.02967,+ 0.3819,- 0.1695,- 0.07444,+ 0.007409,- 0.006262,- 1.836,- 0.01759,- 0.06256,- 0.002371,+ 0.0007947,+ 0.0001501,- 0.8603,- 0.136,- 0.03629,- 0.003706,- 0.0002976,+ 1.857e-05,+ 3.021e-05,+ 2.248,- 0.1178,+ 0.01255,+ 0.001134,- 0.0002161,- 5.817e-06,+ 8.836e-07,- 1.769e-07,+ 0.7313,- 0.1188,+ 0.01145,+ 0.001011,+ 0.0001083,+ 2.57e-06,- 2.14e-06,- 5.71e-08,+ 2e-08,- 1.632,- 0.006948,- 0.003893,+ 0.0008592,+ 7.577e-05,+ 4.539e-06,- 3.852e-07,- 2.213e-07,- 1.37e-08,+ 5.8e-09])
    bw_amp = np.array([+ 0.0,+ 0.0,- 0.08865,+ 0.0,- 0.4309,+ 0.0634,+ 0.0,+ 0.1162,+ 0.06176,- 0.004234,+ 0.0,+ 0.253,+ 0.04017,- 0.006204,+ 0.004977,+ 0.0,- 0.1737,- 0.005638,+ 0.0001488,+ 0.0004857,- 0.0001809,+ 0.0,- 0.1514,- 0.01685,+ 0.005333,- 7.611e-05,+ 2.394e-05,+ 8.195e-06,+ 0.0,+ 0.09326,- 0.01275,- 0.0003071,+ 5.374e-05,- 3.391e-05,- 7.436e-06,+ 6.747e-07,+ 0.0,- 0.08637,- 0.003807,- 0.0006833,- 3.861e-05,- 2.268e-05,+ 1.454e-06,+ 3.86e-07,- 1.068e-07,+ 0.0,- 0.02658,- 0.001947,+ 0.0007131,- 3.506e-05,+ 1.885e-07,+ 5.792e-07,+ 3.99e-08,+ 2e-08,- 5.7e-09])
    np.pi = 3.14159265359
    # reference day is 28 January
# this is taken from Niell (1996) to be consistent
    doy = dmjd - 44239 + 1 - 28
    
    # degree n and order m
    nmax = 9
    mmax = 9
    # unit vector
    x = np.cos(dlat) * np.cos(dlon)
    y = np.cos(dlat) * np.sin(dlon)
    z = np.sin(dlat)
    # Legendre polynomials
    V=pd.DataFrame(np.zeros((nmax+1,nmax+1)))
    W=pd.DataFrame(np.zeros((nmax+1,nmax+1)))
    
    
    V.loc[0,0] = 1
    W.loc[0,0] = 0
    V.loc[1,0] = z * V.loc[0,0]
    W.loc[1,0] = 0
    for n in range(2,nmax+1):
        V.loc[n ,0] = ((2 * n - 1) * z * V.loc[n-1,0] - (n - 1) * V.loc[n-2,0]) / n
        W.loc[n ,0] = 0
    
    for m in range(1,nmax+1):
        V.loc[m ,m ] = (2 * m - 1) * (x * V.loc[m-1,m-1]  - y * W.loc[m-1,m-1] )
        W.loc[m ,m ] = (2 * m - 1) * (x * W.loc[m-1,m-1] +  y * V.loc[m-1,m-1])
        if (m < nmax):
            V.loc[m + 1,m ] = (2 * m + 1) * z * V.loc[m,m]
            W.loc[m + 1,m ] = (2 * m + 1) * z * W.loc[m,m]
        for n in range(m+2,nmax+1):
            V.loc[n ,m ] = ((2 * n - 1) * z * V.loc[n-1,m] - (n + m - 1) * V.loc[n-2,m]) / (n - m)
            W.loc[n ,m ] = ((2 * n - 1) * z * W.loc[n-1,m] - (n + m - 1) * W.loc[n-2,m]) / (n - m)
    
    # (1) hydrostatic mf
    
    bh = 0.0029
    c0h = 0.062
    if (dlat < 0):
        phh = np.pi
        c11h = 0.007
        c10h = 0.002
    else:
        phh = 0
        c11h = 0.005
        c10h = 0.001
    
    ch = c0h + ((np.cos(doy / 365.25 * 2 * np.pi + phh) + 1) * c11h / 2 + c10h) * (1 - np.cos(dlat))
    ahm = 0
    aha = 0
    i = 0
    for n in range(0,nmax+1):
        for m in range(0,n+1):
            
            ahm = ahm + (ah_mean[i] * V.loc[n,m] + bh_mean[i]  * W.loc[n,m])
            aha = aha + (ah_amp[i]  * V.loc[n,m] + bh_amp[i]  * W.loc[n,m])
            i = i + 1
    
    ah = (ahm + aha * np.cos(doy / 365.25 * 2 * np.pi)) * 1e-05
    sine = np.sin(np.pi / 2 - zd)
    cose = np.cos(np.pi / 2 - zd)
    beta = bh / (sine + ch)
    gamma = ah / (sine + beta)
    topcon = (1 + ah / (1 + bh / (1 + ch)))
    gmfh = topcon / (sine + gamma)
    # height correction for hydrostatic mapping function from Niell (1996) in order to reduce the coefficients to sea level
    a_ht = 2.53e-05
    b_ht = 0.00549
    c_ht = 0.00114
    hs_km = dhgt / 1000
    beta = b_ht / (sine + c_ht)
    gamma = a_ht / (sine + beta)
    topcon = (1 + a_ht / (1 + b_ht / (1 + c_ht)))
    ht_corr_coef = 1 / sine - topcon / (sine + gamma)
    ht_corr = ht_corr_coef * hs_km
    gmfh = gmfh + ht_corr
    # (2) wet mf
    
    bw = 0.00146
    cw = 0.04391
    awm = 0
    awa = 0
    i = 0
    for n in range(0,nmax+1):
        for m in range(0,n+1):
            
            awm = awm + (aw_mean[i] * V.loc[n,m] + bw_mean[i] * W.loc[n,m])
            awa = awa + (aw_amp [i] * V.loc[n,m] + bw_amp[i] * W.loc[n,m])
            i = i + 1
    
    aw = (awm + awa * np.cos(doy / 365.25 * 2 * np.pi)) * 1e-05
    beta = bw / (sine + cw)
    gamma = aw / (sine + beta)
    topcon = (1 + aw / (1 + bw / (1 + cw)))
    gmfw = topcon / (sine + gamma)
    return gmfh,gmfw


###############################################################################
def gradient_MF(elev=None):
    
    c=0.0032
    elev=np.radians(elev)
    mfg=1/(np.sin(elev)* np.tan(elev)+c)
    mfg=mfg/1.4
    return mfg
GMF_GradinetMF.py
zohrehadavi/STD_SWD_Calc
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""

import numpy as np
import pandas as pd    
def gpt3_5_fast_readGrid(filename ='gpt3_5.grd'):
# =============================================================================  
#     gpt3_5_fast_readGrid.py
#     This routine reads the grid 'gpt3_5.grd' and overgives the respective
#     parameters in a cell. It must be placed ahead of the for loop, which runs
#    through all observations in order to save a huge time amount. The related
#    function "gpt3_5_fast.py" must then replace the default
#   "gpt3_5.m" at the respective location in the text.
    
    
#    File created by Daniel Landskron, 2017-04-12
    
#    File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================  
    
# read gridfile

    
    file=open(filename)
    C = pd.read_csv(file ,skiprows=[i for i in range(0,1)],sep='\s+',skipfooter =0, skip_blank_lines=True,header = None, index_col= False )
    file.close()   


    
    p_grid = C.loc[:,np.arange(2,7)]
    
    T_grid = C.loc[:,np.arange(7,12)]
    
    Q_grid = C.loc[:,np.arange(12,17)] / 1000
    
    dT_grid= C.loc[:,np.arange(17,22)] / 1000
    
    u_grid = C.loc[:,22]
    
    Hs_grid = C.loc[:,23]
    
    ah_grid = C.loc[:,np.arange(24,29)] / 1000 
    
    aw_grid = C.loc[:,np.arange(29,34)] / 1000 
    
    la_grid = C.loc[:,np.arange(34,39)]
    
    Tm_grid = C.loc[:,np.arange(39,44)]
    
    Gn_h_grid = C.loc[:,np.arange(44,49)] / 100000
    
    Ge_h_grid = C.loc[:,np.arange(49,54)] / 100000
    
    Gn_w_grid = C.loc[:,np.arange(54,59)] / 100000 
    
    Ge_w_grid = C.loc[:,np.arange(59,64)] / 100000
    
    # combine all data to on cell grid
    gpt3_grid={}
    gpt3_grid[0]=p_grid
    gpt3_grid[1]=T_grid
    gpt3_grid[2]=Q_grid
    gpt3_grid[3]=dT_grid
    gpt3_grid[4]=u_grid
    gpt3_grid[5]=Hs_grid
    gpt3_grid[6]=ah_grid
    gpt3_grid[7]=aw_grid
    gpt3_grid[8]=la_grid
    gpt3_grid[9]=Tm_grid
    gpt3_grid[10]=Gn_h_grid
    gpt3_grid[11]=Ge_h_grid
    gpt3_grid[12]=Gn_w_grid
    gpt3_grid[13]=Ge_w_grid
            
    return gpt3_grid

# =============================================================================  
    
def gpt3_5_fast(mjd = None,lat = None,lon = None,h_ell = None,it = None): 
    
    grid=gpt3_5_fast_readGrid(filename ='gpt3_5.grd')
    
#    (c) Department of Geodesy and Geoinformation, Vienna University of Technology, 2017
    
# =============================================================================  
#     This subroutine determines pressure, temperature, temperature lapse rate,
#     mean temperature of the water vapor, water vapour pressure, hydrostatic
#     and wet mapping function coefficients ah and aw, water vapour decrease
#     factor, geoid undulation and empirical tropospheric gradients for
#     specific sites near the earth's surface. All output values are valid for
#     the specified ellipsoidal height h_ell.
#     GPT3_5 is based on a 5°x5° external grid file ('gpt3_5.grd') with mean
#     values as well as sine and cosine amplitudes for the annual and
#     semiannual variation of the coefficients.
#     This file goes together with 'gpt3_5_fast_readGrid.m', which reads the
#     grid and saves it in cell structures, what, especially for large numbers
#     of stations, significantly reduces runtime compared to gpt3_5.m.
#     gpt3_5_fast_readGrid' must be placed ahead of the for loop, in which
#     gpt3_5_fast is positioned.
    
    
#     Reference:
#     D. Landskron, J. Böhm (2018), VMF3/GPT3: Refined Discrete and Empirical Troposphere Mapping Functions,
#     J Geod (2018) 92: 349., doi: 10.1007/s00190-017-1066-2.
#     Download at: https://link.springer.com/content/pdf/10.1007#2Fs00190-017-1066-2.pdf
    
    
#     Input parameters:
#     -------------------    
#     mjd:   modified Julian date (scalar, only one epoch per call is possible)
#     lat:   ellipsoidal latitude in radians [-pi/2:+pi/2] (vector)
#     lon:   longitude in radians [-pi:pi] or [0:2pi] (vector)
#     h_ell: ellipsoidal height in m (vector)
#     it:    case 1: no time variation but static quantities
#     case 0: with time variation (annual and semiannual terms)
    
#     Output parameters:
#     -------------------
#     p:    pressure in hPa (vector)
#     T:    temperature in degrees Celsius (vector)
#     dT:   temperature lapse rate in degrees per km (vector)
#     Tm:   mean temperature weighted with the water vapor in degrees Kelvin (vector)
#     e:    water vapour pressure in hPa (vector)
#     ah:   hydrostatic mapping function coefficient (VMF3) (vector)
#     aw:   wet mapping function coefficient (VMF3) (vector)
#     la:   water vapour decrease factor (vector)
#     undu: geoid undulation in m (vector)
#     Gn_h: hydrostatic north gradient in m (vector)
#     Ge_h: hydrostatic east gradient in m (vector)
#     Gn_w: wet north gradient in m (vector)
#     Ge_w: wet east gradient in m (vector)
    
#     File created by Daniel Landskron, 2017-04-12
#     File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, 2023-02-20
# =============================================================================  
    lat=np.array(lat)
    lon=np.array(lon)
    h_ell=np.array(h_ell)
# convert mjd to doy
    
    hour = int(np.floor((mjd - int(np.floor(mjd))) * 24))
    
    minu = int(np.floor((((mjd - int(np.floor(mjd))) * 24) - hour) * 60))
    
    sec = (((((mjd - int(np.floor(mjd))) * 24) - hour) * 60) - minu) * 60
    
    # change secs, min hour whose sec==60
    
    
    if sec==60:
        minu=minu+1
        sec=0
        
    if minu==60:
        hour=hour+1
        minu=0
    # calc jd (yet wrong for hour==24)
    jd = mjd + 2400000.5
    # if hr==24, correct jd and set hour==0
    if hour==24:
        jd=jd+1
        hour=0
    
    # integer julian date
    jd_int = int(np.floor(jd + 0.5))
    aa = jd_int + 32044
    bb = int(np.floor((4 * aa + 3) / 146097))
    cc = aa - int(np.floor((bb * 146097) / 4))
    dd = int(np.floor((4 * cc + 3) / 1461))
    ee = cc - int(np.floor((1461 * dd) / 4))
    mm = int(np.floor((5 * ee + 2) / 153))
    day = ee - int(np.floor((153 * mm + 2) / 5)) + 1
    month = mm + 3 - 12 * int(np.floor(mm / 10))
    year = bb * 100 + dd - 4800 + int(np.floor(mm / 10))
    # first check if the specified year is leap year or not (logical output)
    leapYear = (np.logical_or((np.mod(year,4) == np.logical_and(0,np.mod(year,100)) != 0),np.mod(year,400)) == 0)
    days = np.array([31,28,31,30,31,30,31,31,30,31,30,31])
    doy = sum(days[range(0,month-1)]) + day
    if leapYear == 1 and month > 2:
        doy = doy + 1
    
    doy = doy + mjd - int(np.floor(mjd))
    
    # determine the GPT3 coefficients
    
    # mean gravity in m/s**2
    gm = 9.80665
    # molar mass of dry air in kg/mol
    dMtr = 28.965 * 10 ** - 3
    # universal gas constant in J/K/mol
    Rg = 8.3143
    # factors for amplitudes
    if (it == 1):
        cosfy = 0
        coshy = 0
        sinfy = 0
        sinhy = 0
    else:
        cosfy = np.cos(doy / 365.25 * 2 * np.pi)
        coshy = np.cos(doy / 365.25 * 4 * np.pi)
        sinfy = np.sin(doy / 365.25 * 2 * np.pi)
        sinhy = np.sin(doy / 365.25 * 4 * np.pi)
    
    # read the respective data from the grid
    p_grid = grid[0]
    
    T_grid = grid[1]
    
    Q_grid = grid[2]
    
    dT_grid = grid[3]
    
    u_grid = grid[4]
    
    Hs_grid = grid[5]
    
    ah_grid = grid[6]
    
    aw_grid = grid[7]
    
    la_grid = grid[8]
    
    Tm_grid = grid[9]
    
    Gn_h_grid = grid[10]
    
    Ge_h_grid = grid[11]
    
    Gn_w_grid = grid[12]
    
    Ge_w_grid = grid[13]
    
    # determine the number of stations
    
    nstat = lat.size
    # initialization
    p    = np.zeros(np.array([nstat,1]))
    T    = np.zeros(np.array([nstat,1]))
    dT   = np.zeros(np.array([nstat,1]))
    Tm   = np.zeros(np.array([nstat,1]))
    e    = np.zeros(np.array([nstat,1]))
    ah   = np.zeros(np.array([nstat,1]))
    aw   = np.zeros(np.array([nstat,1]))
    la   = np.zeros(np.array([nstat,1]))
    undu = np.zeros(np.array([nstat,1]))
    Gn_h = np.zeros(np.array([nstat,1]))
    Ge_h = np.zeros(np.array([nstat,1]))
    Gn_w = np.zeros(np.array([nstat,1]))
    Ge_w = np.zeros(np.array([nstat,1]))
    # loop over stations
    for k in range(0,nstat):
        # only positive longitude in degrees
        if lon < 0:
            plon = (lon + 2 * np.pi) * 180 / np.pi
        else:
            plon = lon * 180 / np.pi
        # transform to polar distance in degrees
        ppod = (- lat + np.pi / 2) * 180 / np.pi
        # find the index (line in the grid file) of the nearest point
        ipod = int(np.floor((ppod + 5) / 5))
        ilon = int(np.floor((plon + 5) / 5))
        # normalized (to one) differences, can be positive or negative
        diffpod = (ppod - (ipod * 5 - 2.5)) / 5
        difflon = (plon - (ilon * 5 - 2.5)) / 5
        if ipod == 36:
            ipod = 35
        if ilon == 72:
            ilon = 0
        if ilon == -1:
            ilon = 71
        indx=[]
        # get the number of the corresponding line
        indx.append((ipod - 1) * 72 + ilon-1)
        # near the poles: nearest neighbour interpolation, otherwise: bilinear
        bilinear = 0
        if ppod > 2.5 and ppod < 177.5:
            bilinear = 1
        # case of nearest neighbourhood
        if bilinear == 0:
            ix = indx[0]
            # transforming ellipsoidal height to orthometric height
            undu[k] = u_grid[ix]
            hgt = h_ell - undu[k]
            # pressure, temperature at the height of the grid
            T0 = T_grid.iloc[ix,0] + T_grid.iloc[ix,1] * cosfy + T_grid.iloc[ix,2] * sinfy + T_grid.iloc[ix,3] * coshy + T_grid.iloc[ix,4] * sinhy
            p0 = p_grid.iloc[ix,0] + p_grid.iloc[ix,1] * cosfy + p_grid.iloc[ix,2] * sinfy + p_grid.iloc[ix,3] * coshy + p_grid.iloc[ix,4] * sinhy
            # specific humidity
            Q = Q_grid.iloc[ix,0] + Q_grid.iloc[ix,1] * cosfy + Q_grid.iloc[ix,2] * sinfy + Q_grid.iloc[ix,3] * coshy + Q_grid.iloc[ix,4] * sinhy
            # lapse rate of the temperature
            dT[k] = dT_grid.iloc[ix,0] + dT_grid.iloc[ix,1] * cosfy + dT_grid.iloc[ix,2] * sinfy + dT_grid.iloc[ix,3] * coshy + dT_grid.iloc[ix,4] * sinhy
            # station height - grid height
            redh = hgt - Hs_grid[ix]
            # temperature at station height in Celsius
            T[k] = T0 + dT[k] * redh - 273.15
            # temperature lapse rate in degrees / km
            dT[k] = dT[k] * 1000
            # virtual temperature in Kelvin
            Tv = T0 * (1 + 0.6077 * Q)
            c = gm * dMtr / (Rg * Tv)
            # pressure in hPa
            p[k] = (p0 * np.exp(- c * redh)) / 100
            # hydrostatic and wet coefficients ah and aw
            ah[k] = ah_grid.iloc[ix,0] + ah_grid.iloc[ix,1] * cosfy + ah_grid.iloc[ix,2] * sinfy + ah_grid.iloc[ix,3] * coshy + ah_grid.iloc[ix,4] * sinhy
            aw[k] = aw_grid.iloc[ix,0] + aw_grid.iloc[ix,1] * cosfy + aw_grid.iloc[ix,2] * sinfy + aw_grid.iloc[ix,3] * coshy + aw_grid.iloc[ix,4] * sinhy
            # water vapour decrease factor la
            la[k] = la_grid.iloc[ix,0] + la_grid.iloc[ix,1] * cosfy + la_grid.iloc[ix,2] * sinfy + la_grid.iloc[ix,3] * coshy + la_grid.iloc[ix,4] * sinhy
            # mean temperature Tm
            Tm[k] = Tm_grid.iloc[ix,0] + Tm_grid.iloc[ix,1] * cosfy + Tm_grid.iloc[ix,2] * sinfy + Tm_grid.iloc[ix,3] * coshy + Tm_grid.iloc[ix,4] * sinhy
            # north and east gradients (total, hydrostatic and wet)
            Gn_h[k] = Gn_h_grid.iloc[ix,0] + Gn_h_grid.iloc[ix,1] * cosfy + Gn_h_grid.iloc[ix,2] * sinfy + Gn_h_grid.iloc[ix,3] * coshy + Gn_h_grid.iloc[ix,4] * sinhy
            Ge_h[k] = Ge_h_grid.iloc[ix,0] + Ge_h_grid.iloc[ix,1] * cosfy + Ge_h_grid.iloc[ix,2] * sinfy + Ge_h_grid.iloc[ix,3] * coshy + Ge_h_grid.iloc[ix,4] * sinhy
            Gn_w[k] = Gn_w_grid.iloc[ix,0] + Gn_w_grid.iloc[ix,1] * cosfy + Gn_w_grid.iloc[ix,2] * sinfy + Gn_w_grid.iloc[ix,3] * coshy + Gn_w_grid.iloc[ix,4] * sinhy
            Ge_w[k] = Ge_w_grid.iloc[ix,0] + Ge_w_grid.iloc[ix,1] * cosfy + Ge_w_grid.iloc[ix,2] * sinfy + Ge_w_grid.iloc[ix,3] * coshy + Ge_w_grid.iloc[ix,4] * sinhy
            # water vapor pressure in hPa
            e0 = Q * p0 / (0.622 + 0.378 * Q) / 100
            e[k] = e0 * (100 * p[k] / p0) ** (la[k] + 1)
        else:
            ipod1 = int(ipod + np.sign(diffpod))
            ilon1 = int(ilon + np.sign(difflon))
            if ilon1 == 72:
                ilon1 = 0
            if ilon1 == -1:
                ilon1 = 71
            # get the number of the line
            indx.append((ipod1 - 1) * 72 + ilon-1)
            indx.append((ipod - 1) * 72 + ilon1-1)
            indx.append((ipod1 - 1) * 72 + ilon1-1)
            # transforming ellipsoidal height to orthometric height: Hortho = -N + Hell
            undul = u_grid[indx]
            hgt = h_ell - undul
            # pressure, temperature at the height of the grid
            T0 = T_grid.iloc[indx,0] + T_grid.iloc[indx,1] * cosfy + T_grid.iloc[indx,2] * sinfy + T_grid.iloc[indx,3] * coshy + T_grid.iloc[indx,4] * sinhy
            p0 = p_grid.iloc[indx,0] + p_grid.iloc[indx,1] * cosfy + p_grid.iloc[indx,2] * sinfy + p_grid.iloc[indx,3] * coshy + p_grid.iloc[indx,4] * sinhy
            # humidity
            Ql = Q_grid.iloc[indx,0] + Q_grid.iloc[indx,1] * cosfy + Q_grid.iloc[indx,2] * sinfy + Q_grid.iloc[indx,3] * coshy + Q_grid.iloc[indx,4] * sinhy
            # reduction = stationheight - gridheight
            Hs1 = Hs_grid[indx]
            redh = hgt - Hs1
            # lapse rate of the temperature in degree / m
            dTl = dT_grid.iloc[indx,0] + dT_grid.iloc[indx,1] * cosfy + dT_grid.iloc[indx,2] * sinfy + dT_grid.iloc[indx,3] * coshy + dT_grid.iloc[indx,4] * sinhy
            # temperature reduction to station height
            Tl = T0 + np.multiply(dTl,redh) - 273.15
            # virtual temperature
            Tv = np.multiply(T0,(1 + 0.6077 * Ql))
            c = gm * dMtr / (Rg * Tv)
            # pressure in hPa
            pl = (np.multiply(p0,np.exp(np.multiply(- c,redh)))) / 100
            # hydrostatic and wet coefficients ah and aw
            ahl = ah_grid.iloc[indx,0] + ah_grid.iloc[indx,1] * cosfy + ah_grid.iloc[indx,2] * sinfy + ah_grid.iloc[indx,3] * coshy + ah_grid.iloc[indx,4] * sinhy
            awl = aw_grid.iloc[indx,0] + aw_grid.iloc[indx,1] * cosfy + aw_grid.iloc[indx,2] * sinfy + aw_grid.iloc[indx,3] * coshy + aw_grid.iloc[indx,4] * sinhy
            # water vapour decrease factor la
            lal = la_grid.iloc[indx,0] + la_grid.iloc[indx,1] * cosfy + la_grid.iloc[indx,2] * sinfy + la_grid.iloc[indx,3] * coshy + la_grid.iloc[indx,4] * sinhy
            # mean temperature of the water vapor Tm
            Tml = Tm_grid.iloc[indx,0] + Tm_grid.iloc[indx,1] * cosfy + Tm_grid.iloc[indx,2] * sinfy + Tm_grid.iloc[indx,3] * coshy + Tm_grid.iloc[indx,4] * sinhy
            # north and east gradients (total, hydrostatic and wet)
            Gn_hl = Gn_h_grid.iloc[indx,0] + Gn_h_grid.iloc[indx,1] * cosfy + Gn_h_grid.iloc[indx,2] * sinfy + Gn_h_grid.iloc[indx,3] * coshy + Gn_h_grid.iloc[indx,4] * sinhy
            Ge_hl = Ge_h_grid.iloc[indx,0] + Ge_h_grid.iloc[indx,1] * cosfy + Ge_h_grid.iloc[indx,2] * sinfy + Ge_h_grid.iloc[indx,3] * coshy + Ge_h_grid.iloc[indx,4] * sinhy
            Gn_wl = Gn_w_grid.iloc[indx,0] + Gn_w_grid.iloc[indx,1] * cosfy + Gn_w_grid.iloc[indx,2] * sinfy + Gn_w_grid.iloc[indx,3] * coshy + Gn_w_grid.iloc[indx,4] * sinhy
            Ge_wl = Ge_w_grid.iloc[indx,0] + Ge_w_grid.iloc[indx,1] * cosfy + Ge_w_grid.iloc[indx,2] * sinfy + Ge_w_grid.iloc[indx,3] * coshy + Ge_w_grid.iloc[indx,4] * sinhy
            # water vapor pressure in hPa
            e0 = np.multiply(Ql,p0) / (0.622 + 0.378 * Ql) / 100
            el = np.multiply(e0,(100.0 * pl / p0) ** (lal + 1))
            dnpod1 = np.abs(diffpod)
            dnpod2 = 1 - dnpod1
            dnlon1 = np.abs(difflon)
            dnlon2 = 1 - dnlon1
            # pressure
            R1 = dnpod2 * pl.iloc[0] + dnpod1 * pl.iloc[1]
            R2 = dnpod2 * pl.iloc[2] + dnpod1 * pl.iloc[3]
            p[k] = dnlon2 * R1 + dnlon1 * R2
            # temperature
            R1 = dnpod2 * Tl.iloc[0] + dnpod1 * Tl.iloc[1]
            R2 = dnpod2 * Tl.iloc[2] + dnpod1 * Tl.iloc[3]
            T[k] = dnlon2 * R1 + dnlon1 * R2
            # temperature in degree per km
            R1 = dnpod2 * dTl.iloc[0] + dnpod1 * dTl.iloc[1]
            R2 = dnpod2 * dTl.iloc[2] + dnpod1 * dTl.iloc[3]
            dT[k] = (dnlon2 * R1 + dnlon1 * R2) * 1000
            # water vapor pressure in hPa
            R1 = dnpod2 * el.iloc[0] + dnpod1 * el.iloc[1]
            R2 = dnpod2 * el.iloc[2] + dnpod1 * el.iloc[3]
            e[k] = dnlon2 * R1 + dnlon1 * R2
            # ah and aw
            R1 = dnpod2 * ahl.iloc[0] + dnpod1 * ahl.iloc[1]
            R2 = dnpod2 * ahl.iloc[2] + dnpod1 * ahl.iloc[3]
            ah[k] = dnlon2 * R1 + dnlon1 * R2
            R1 = dnpod2 * awl.iloc[0] + dnpod1 * awl.iloc[1]
            R2 = dnpod2 * awl.iloc[2] + dnpod1 * awl.iloc[3]
            aw[k] = dnlon2 * R1 + dnlon1 * R2
            # undulation
            R1 = dnpod2 * undul.iloc[0] + dnpod1 * undul.iloc[1]
            R2 = dnpod2 * undul.iloc[2] + dnpod1 * undul.iloc[3]
            undu[k] = dnlon2 * R1 + dnlon1 * R2
            # water vapor decrease factor la
            R1 = dnpod2 * lal.iloc[0] + dnpod1 * lal.iloc[1]
            R2 = dnpod2 * lal.iloc[2] + dnpod1 * lal.iloc[3]
            la[k] = dnlon2 * R1 + dnlon1 * R2
            # gradients
            R1 = dnpod2 * Gn_hl.iloc[0] + dnpod1 * Gn_hl.iloc[1]
            R2 = dnpod2 * Gn_hl.iloc[2] + dnpod1 * Gn_hl.iloc[3]
            Gn_h[k] = (dnlon2 * R1 + dnlon1 * R2)
            R1 = dnpod2 * Ge_hl.iloc[0] + dnpod1 * Ge_hl.iloc[1]
            R2 = dnpod2 * Ge_hl.iloc[2] + dnpod1 * Ge_hl.iloc[3]
            Ge_h[k] = (dnlon2 * R1 + dnlon1 * R2)
            R1 = dnpod2 * Gn_wl.iloc[0] + dnpod1 * Gn_wl.iloc[1]
            R2 = dnpod2 * Gn_wl.iloc[2] + dnpod1 * Gn_wl.iloc[3]
            Gn_w[k] = (dnlon2 * R1 + dnlon1 * R2)
            R1 = dnpod2 * Ge_wl.iloc[0] + dnpod1 * Ge_wl.iloc[1]
            R2 = dnpod2 * Ge_wl.iloc[2] + dnpod1 * Ge_wl.iloc[3]
            Ge_w[k] = (dnlon2 * R1 + dnlon1 * R2)
            # mean temperature of the water vapor Tm
            R1 = dnpod2 * Tml.iloc[0] + dnpod1 * Tml.iloc[1]
            R2 = dnpod2 * Tml.iloc[2] + dnpod1 * Tml.iloc[3]
            Tm[k] = dnlon2 * R1 + dnlon1 * R2
    
    return p,T,dT,Tm,e,ah,aw,la,undu,Gn_h,Ge_h,Gn_w,Ge_w
GPT3.py
zohrehadavi/STD_SWD_Calc
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""
import pandas as pd
import numpy as np
import time
import pyproj
import linecache
from astropy.time import Time

transformer = pyproj.Transformer.from_crs(
    {"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
    {"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
    )

def Readwrite_TRO_ZTD(file_source=None,file_name=None,file_dest=None,fle_out=None,doy=None,int_need=None,New_rate=None):
# -------------Input Parameters
# file_source: The directory where TRO files are located
# file_name  : The full name of TRO file
# file_dest  : Where the created file should be saved
# fle_out    : The name of output file
# doy        : Day of year
# int_need   : 1 if interpolation is needed and 0 otherwise (default=0)
# New_rate   : New interpolation interval [second] (default=150 s #2.5 minutes)

# File created by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================    

####SAMPLING TROP
 line0='SAMPLING TROP'
 directory_trp=file_source+'/'+file_name
 with open(directory_trp) as f:
          pos = 0
          cur_line = f.readline()
          index_pos=cur_line.find(line0)
          while (index_pos<0):
                pos = pos+1
                cur_line = f.readline()
                index_pos=cur_line.find(line0)
 f.close()
 

 particular_line = linecache.getline(directory_trp, pos+1)
 sample_tro=int(particular_line[50:53])
  
####Reading CRD part
 Columnsmain0=['SITE','lon','lat','h','X','Y','Z']
 Columns_crd=['SITE','PT','SOLN','T','STA_X','STA_Y','STA_Z','SYSTEM','REMRK']
 line0='+TROP/STA_COORDINATES'
 appended_Data0 = []  
 directory_trp=file_source+'/'+file_name
 with open(directory_trp) as f:
          pos = 0
          cur_line = f.readline()
          index_pos=cur_line.find(line0)
          
          while index_pos<0:
                pos = pos+1
                cur_line = f.readline()
                index_pos=cur_line.find(line0)

 f.close()
 
 ##
 line0='-TROP/STA_COORDINATES'
 with open(directory_trp) as f:
          pose = 0
          cur_line = f.readline()
          index_pos=cur_line.find(line0)
          
          while index_pos<0:
                pose = pose+1
                cur_line = f.readline()
                index_pos=cur_line.find(line0)

 f.close()
 
 
 ##
 
 fle=open(directory_trp)
 file_crd = fle.readlines()[pos+2:pose] 
 fle.close()  

 df = pd.DataFrame({'Sta_inf': file_crd})
 df[['SITE','PT','SOLN','T','STA_X','STA_Y','STA_Z','SYSTEM','REMRK']] = df['Sta_inf'].str.split(expand=True)
 df = df.drop('Sta_inf', axis=1)
    
 Mainfle=pd.DataFrame(np.zeros((len(df),len(Columnsmain0))),columns=Columnsmain0) 
 lon1, lat1, alt1 = transformer.transform(df['STA_X'],df['STA_Y'],df['STA_Z'],radians=False) 
 Mainfle.iloc[:,0]=df['SITE']
 Mainfle.iloc[:,1]=lon1
 Mainfle.iloc[:,2]=lat1
 Mainfle.iloc[:,3]=alt1
 Mainfle.iloc[:,4]=df['STA_X']
 Mainfle.iloc[:,5]=df['STA_Y']
 Mainfle.iloc[:,6]=df['STA_Z']
 
 appended_Data0.append(Mainfle)
 del df
 del fle
 del Mainfle
 finaldf0 = pd.concat(appended_Data0)
####Reading TRP part
 Columnsmain=['SITE','Year','DoY','hour','minute','second','ZTD','GN','GE']
 Columns_trp=['SITE','EPOCH','TROTOT','STDDEVt','TGNTOT','STDDEVn','TGETOT','STDDEVe']
 line='+TROP/SOLUTION' 
 directory_trp=file_source+'/'+file_name
 with open(directory_trp) as f:
          pos = 0
          cur_line = f.readline()
         
          while not cur_line.startswith(line):
                pos = pos+1
                cur_line = f.readline()
          f.seek(pos)
          f.close() 
 fle=open(directory_trp)                      
 df = pd.read_csv(fle ,skiprows=[i for i in range(0,pos+1)],sep='\s+',skipfooter = 2, skip_blank_lines=True,header = 0, index_col= False ,names=Columns_trp)
 fle.close()     
 
 
 
 
 
 
 EPOCHall=df['EPOCH'].str[7:12].astype(int)
 DoY=df['EPOCH'].str[3:6].astype(int)
 Year=df['EPOCH'].str[0:2]
 
 finaldf1=pd.DataFrame(np.zeros((0,len(Columnsmain))),columns=Columnsmain)  
 jt=0
 for it in range(0,len(EPOCHall)):
     
     if DoY[it]==doy:
        convert = time.strftime("%H:%M:%S", time.gmtime(EPOCHall[it]))
     

        finaldf1.loc[jt,'SITE']=df.loc[it,'SITE']
        finaldf1.loc[jt,'Year']=int(str(20)+Year[it])
        finaldf1.loc[jt,'DoY']=DoY[it]
        finaldf1.loc[jt,'hour']=convert[0:2]
        finaldf1.loc[jt,'minute']=convert[3:5]
        finaldf1.loc[jt,'second']=convert[6:8]
        finaldf1.loc[jt,'ZTD']=df.loc[it,'TROTOT']
        finaldf1.loc[jt,'GN']=df.loc[it,'TGETOT']
        finaldf1.loc[jt,'GE']=df.loc[it,'TGNTOT']
        jt=jt+1
 
 
 del df
 del fle
########Concating CRD to TRP File
 year=int(finaldf1['Year'].unique())
 Columnsmain_new=['SITE','lon','lat','h','X','Y','Z','hour','minute','second','ZTD','GE','GN']
 
 appended_Data = [] 
 for it in range(0,len(finaldf0)):
    jloc0=finaldf1.loc[(finaldf1['SITE'] ==finaldf0.iloc[it,0])]
    jloc0.reset_index(drop=True, inplace=True)
    Mainfle=pd.DataFrame(np.zeros((len(jloc0),len(Columnsmain_new))),columns=Columnsmain_new)
    Mainfle['SITE']= jloc0['SITE']
    Mainfle['lon']= finaldf0.iloc[it,1]
    Mainfle['lat']= finaldf0.iloc[it,2]  
    Mainfle['h']= finaldf0.iloc[it,3]
    Mainfle['X']= finaldf0.iloc[it,4]
    Mainfle['Y']= finaldf0.iloc[it,5]
    Mainfle['Z']= finaldf0.iloc[it,6]
    
    Mainfle['hour']= jloc0['hour']
    Mainfle['minute']= jloc0['minute']
    Mainfle['second']= jloc0['second']
    Mainfle['ZTD']= jloc0['ZTD']
    Mainfle['GE']= jloc0['GE']
    Mainfle['GN']= jloc0['GN']

    Mainfle.reset_index(drop=True, inplace=True)
    appended_Data.append(Mainfle)
    del jloc0,Mainfle

 finaldf_main = pd.concat(appended_Data)
 
########Interpolation if needed
  #defult: 2.5 minutes
  
 Int_fle=pd.DataFrame(np.zeros((0,len(Columnsmain_new))),columns=Columnsmain_new) 
 nt=0 
 Sta_name=finaldf_main['SITE'].unique()

#######Interpolation of ZTD, GN, and GE
 for it in range(0,len(Sta_name)):
         jloc0=finaldf_main.loc[(finaldf_main['SITE']==Sta_name[it])]
         jloc0.reset_index(drop=True, inplace=True)
         MJDd=[]
         interval_data=[]
         for jt in range(0,len(jloc0)):
             t = Time('{}:{}:{}:{}'.format(year, doy, int(jloc0.loc[jt,'hour']), int(jloc0.loc[jt,'minute']), int(jloc0.loc[jt,'second'])), format='yday', scale='utc')
             MJDd.append(t.mjd)
             interval_data.append((int(jloc0.loc[jt,'hour'])*3600)+(int(jloc0.loc[jt,'minute'])*60)+(int(jloc0.loc[jt,'second'])))
             
         
         intervald=np.array(interval_data)
         dmjd= pd.DataFrame({'values': MJDd})
         del MJDd
         del interval_data
         
         #mjd for 2.5 minutes
         
         convert_par=(New_rate/(3600*24))
         mjd_new=pd.DataFrame(np.arange(float(dmjd.min()),float(dmjd.max()),convert_par,dtype=np.float64))
         interval = range(intervald.min(),intervald.max()+New_rate,New_rate)
         
         ZTD_interp=np.interp(mjd_new.iloc[:,0],dmjd.iloc[:,0], jloc0['ZTD'])
         GE_interp =np.interp(mjd_new.iloc[:,0],dmjd.iloc[:,0], jloc0['GE'])
         GN_interp =np.interp(mjd_new.iloc[:,0],dmjd.iloc[:,0], jloc0['GN'])
         
         
         
         step_fix=int(sample_tro/New_rate)
         lt=0
         def_xs=[]
         for jt in range(0,len(ZTD_interp),step_fix):
             #if len(jloc0)>lt:
              def_xs.append(ZTD_interp[jt]-jloc0.loc[lt,'ZTD'])
              ZTD_interp[jt]=jloc0.loc[lt,'ZTD']
              GE_interp[jt]=jloc0.loc[lt,'GE']
              GN_interp[jt]=jloc0.loc[lt,'GN']
              lt=lt+1
             
         
         for kt in range(0,len(ZTD_interp)):
                convert = time.strftime("%H:%M:%S", time.gmtime(interval[kt]))
                Int_fle.loc[nt,'SITE']=Sta_name[it]
                Int_fle.loc[nt,'lon']=jloc0.iloc[0,1]
                Int_fle.loc[nt,'lat']=jloc0.iloc[0,2]
                Int_fle.loc[nt,'h']=jloc0.iloc[0,3]
                Int_fle.loc[nt,'X']=jloc0.iloc[0,4]
                Int_fle.loc[nt,'Y']=jloc0.iloc[0,5]
                Int_fle.loc[nt,'Z']=jloc0.iloc[0,6]
                
                Int_fle.loc[nt,'hour']=convert[0:2]
                Int_fle.loc[nt,'minute']=convert[3:5]
                Int_fle.loc[nt,'second']=convert[6:8]
                
                Int_fle.loc[nt,'ZTD']=ZTD_interp[kt]
                Int_fle.loc[nt,'GE']=GE_interp[kt]  
                Int_fle.loc[nt,'GN']=GN_interp[kt]  
             
                nt=nt+1
         del ZTD_interp,GE_interp,GN_interp

#########################################################
 if int_need>0:
    finaldf=Int_fle
 else:
    finaldf=finaldf_main
    
 output_path=file_dest+fle_out
 # finaldf.to_csv(output_path, columns = False, index = False)
 with open(output_path, 'a') as f:
    df_string = finaldf.to_string(header=False, index=False)
    f.write(df_string)
    
 return finaldf
Read_Write_TRP.py
zohrehadavi/STD_SWD_Calc
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""
import wget
import os
import numpy as np
import pandas as pd
import os
from scipy.interpolate import griddata
import random
import calendar
#This script download VMF files from VMF server.Here, only operational VMF1 data is downloading. 
# VMF1_OP should be changed to VMF1_FC to download forecating type.
def mjd_to_ymd(mjd):
    jd = mjd + 2400000.5
    jd_int = int(jd)
    fraction = jd - jd_int

    # Calculate year, month, and day
    a = int((jd_int - 1867216.25) / 36524.25)
    b = jd_int + 1 + a - int(a / 4)
    c = b + 1524
    d = int((c - 122.1) / 365.25)
    e = int(365.25 * d)
    f = int((c - e) / 30.6001)

    day = int(c - e - int(30.6001 * f) + fraction)
    month = f - 1 - 12 * int(f / 14)
    year = d - 4715 - int((7 + month) / 10)
    return year,month,day



def download_VMF(year=None,month=None,day=None,directory_vmf=None): 

   LISTt=['00','06','12','18']
   for lstnum in range(0,len(LISTt)):
         try:
            filename='VMFG_'+str(year)+str(month).zfill(2)+str(day).zfill(2)+'.H'+LISTt[lstnum] 
            url="https://vmf.geo.tuwien.ac.at/trop_products/GRID/2.5x2/VMF1/VMF1_OP/"+str(year)+'/'+filename
            if os.path.exists(filename):
               os.remove(filename) # if exist, remove it directly
            wget.download(url,directory_vmf)
         except:
             pass
         
# =============================================================================                   

def read_vmf1_grid(mjd = None,ell = None): 
# =============================================================================
#     Function read_vmf1_grid searches for grid file close in time and computes
#     mapping coefficients ah, aw for selected station

#     Input parameters:
#     -------------------     
#     mjd   ... Modified Julian Date
#     ell   ... ellipsoidal station coordinates lat [degrees], lon [degrees] ,hell [m] 'Should be in list form e.g 'ell=[42.72,56.13,570]'' 
#     File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, 2023-02-20
# =============================================================================
    
# Closest 6 hour period
    mjd=np.array(mjd)
    mjd_6h = np.round(mjd * 4) / 4

# List of relevant 6 hour intervals
    mjd_6h_list = np.unique(mjd_6h)

# Read vmf1 grid of period before and after current date
# (needed for temporal interpolation)
    
    if len(mjd_6h_list) == 1:
        mjd_6h_list = np.array([mjd_6h_list - 0.25,mjd_6h_list,mjd_6h_list + 0.25])
    
# Initialise variable
    ah_ell =[]
    aw_ell =[]
    
    zhd_ell =[]
    zwd_ell = []
    for i in range(0,len(mjd_6h_list)):
        # Convert mjd to date
        yr,mn,dy = mjd_to_ymd(np.rint(np.amin(mjd_6h_list[i])))
        # Hour of the day
        pwd=os.getcwd()
        directory_vmf_down=pwd+"\\Mapping_Fcn\\vmf1\\"
        hh = (mjd_6h_list[i]- np.fix(mjd_6h_list[i])) * 24
        if int(hh)==-6:

         #download_VMF(int(yr),int(mn),int(dy-1),directory_vmf_down)
         hour=18
         day_new=int(dy-1)
        else:
         #download_VMF(int(yr),int(mn),int(dy),directory_vmf_down)
         hour=int(hh)
         day_new=int(dy)
                  
        if day_new==0:
           mn=mn-1
           _,day_new=calendar.monthrange(yr, mn)
        
        # Open grid file with mapping coefficients
        directory_vmf=pwd+"\\Mapping_Fcn\\vmf1\\VMFG_"+str(int(yr))+str(int(mn)).zfill(2)+str(day_new).zfill(2)+'.H'+str(hour).zfill(2)
        
        if not os.path.exists(directory_vmf):
           download_VMF(int(yr),int(mn),int(day_new),directory_vmf_down)
        
        file=open(directory_vmf)
        
        ###############################
        pos = 0
        line0='! Comment:'
        cur_line = file.readline()
        Columns=['lat', 'lon', 'ah', 'aw', 'zhd', 'zwd']
        while not cur_line.startswith(line0):
               pos = pos+1
               cur_line = file.readline()
        
        file.close()  
        file=open(directory_vmf)
        dat = pd.read_csv(file ,skiprows=[i for i in range(0,pos)],sep='\s+',skipfooter =0, skip_blank_lines=True,header = 0, index_col= False ,names=Columns)
        file.close()    

        
        # Restore grid data
        lat = np.array(dat.iloc[:,0])

        lon = np.array(dat.iloc[:,1])
        
        ah = np.array(dat.iloc[:,2])
        
        aw = np.array(dat.iloc[:,3])
        
        zhd = np.array(dat.iloc[:,4])
        
        zwd = np.array(dat.iloc[:,5])
        # Interpolate data to station coordinates
        ah_ell.append(griddata((lat,lon),ah,(ell[0],ell[1])))
        
        aw_ell.append(griddata((lat,lon),aw,(ell[0],ell[1])))
        
        zhd_ell.append(griddata((lat,lon),zhd,(ell[0],ell[1])))
        
        zwd_ell.append(griddata((lat,lon),zwd,(ell[0],ell[1])))
        
        
    mjdd=pd.DataFrame(mjd_6h_list)
    
    ah_ell_new=pd.DataFrame(ah_ell)
    aw_ell_new=pd.DataFrame(aw_ell)
    zhd_ell_new=pd.DataFrame(zhd_ell)
    zwd_ell_new=pd.DataFrame(zwd_ell)
    # Interpolate and store final coefficients to return
    ah_vec=np.interp(np.transpose(mjd),mjdd.iloc[:,0], ah_ell_new.iloc[:,0])
    aw_vec=np.interp(np.transpose(mjd),mjdd.iloc[:,0], aw_ell_new.iloc[:,0])
    
    zhd_vec=np.interp(np.transpose(mjd),mjdd.iloc[:,0], zhd_ell_new.iloc[:,0])
    zwd_vec=np.interp(np.transpose(mjd),mjdd.iloc[:,0], zwd_ell_new.iloc[:,0])
    
    return ah_vec,aw_vec,zhd_vec,zwd_vec
   

# =============================================================================
def vmf1_ht(dmjd = None,dlat = None,dlon=None,ht = None,zd = None,ah= None,aw= None): 
# =============================================================================
#     !!! This is the version with height correction !!!
#     !!! It has to be used with the grid !!!
#
#     This subroutine determines the VMF1 (Vienna Mapping Functions 1)
#     Reference: Boehm, J., B. Werl, H. Schuh (2006),
#     Troposphere mapping functions for GPS and very long baseline interferometry
#     from European Centre for Medium-Range Weather Forecasts operational analysis data,
#     J. Geoph. Res., Vol. 111, B02406, doi:10.1029/2005JB003629.
#
#     Please mind that the coefficients in this paper are wrong. The corrected version of
#     the paper can be found at:
#     http://ggosatm.hg.tuwien.ac.at/DOCS/PAPERS/2006Boehm_etal_VMF1.pdf
#
#     input data
#     ----------
#     ah:   hydrostatic coefficient a (http://ggosatm.hg.tuwien.ac.at/DELAY/GRID/)
#     aw:   wet coefficient a         (http://ggosatm.hg.tuwien.ac.at/DELAY/GRID/)
#     dmjd: modified julian date
#     dlat: ellipsoidal latitude in radians
#     dlon: ellipsoidal longitude in radians
#     ht:   ellipsoidal height in meter
#     zd:   zenith distance in radians
#
#     output data
#     -----------
#     vmf1h: hydrostatic mapping function
#     vmf1w: wet mapping function
#
#     Johannes Boehm, 2005 October 2
#      Rev 2011 July 21: latitude -> ellipsoidal latitude

#     implicit double precision (a-h,o-z)
    
#     pi = 3.14159265359d0
    
#     reference day is 28 January
#     this is taken from Niell (1996) to be consistent
#     File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, 2023-02-20
# =============================================================================
    # ell=[np.degrees(dlat),np.degrees(dlon),ht]
    # ah,aw,zhd,zwd=read_vmf1_grid(dmjd,ell)
    doy = dmjd - 44239.0 + 1 - 28
    bh = 0.0029
    c0h = 0.062
    if (dlat < 0):
        phh = np.pi
        c11h = 0.007
        c10h = 0.002
    else:
        phh = 0.0
        c11h = 0.005
        c10h = 0.001
    
    ch = c0h + ((np.cos(doy / 365.25 * 2.0 * np.pi + phh) + 1.0) * c11h / 2.0 + c10h) * (1.0 - np.cos(dlat))
    sine = np.sin(np.pi / 2.0 - zd)
    beta = bh / (sine + ch)
    gamma = ah / (sine + beta)
    topcon = (1.0 + ah / (1.0 + bh / (1.0 + ch)))
    vmf1h = topcon / (sine + gamma)
    # C  height correction for hydrotatic part [Niell, 1996]
    a_ht = 2.53e-05
    b_ht = 0.00549
    c_ht = 0.00114
    hs_km = ht / 1000.0
    beta = b_ht / (sine + c_ht)
    gamma = a_ht / (sine + beta)
    topcon = (1.0 + a_ht / (1.0 + b_ht / (1.0 + c_ht)))
    ht_corr_coef = 1.0 / sine - topcon / (sine + gamma)
    ht_corr = ht_corr_coef*hs_km
    vmf1h = vmf1h + ht_corr
    bw = 0.00146
    cw = 0.04391
    beta = bw / (sine + cw)
    gamma = aw / (sine + beta)
    topcon = (1.0 + aw / (1.0 + bw / (1.0 + cw)))
    vmf1w = topcon / (sine + gamma)
    return vmf1h,vmf1w

# =============================================================================
def vmf1(dmjd = None,dlat = None,dlon=None,ht = None,zd = None,ah= None,aw= None):
# =============================================================================
#   This subroutine determines the VMF1 (Vienna Mapping Functions 1) for specific sites.
#   Reference: Boehm, J., B. Werl, H. Schuh (2006),
#   Troposphere mapping functions for GPS and very long baseline interferometry
#   from European Centre for Medium-Range Weather Forecasts operational analysis data,
#   J. Geoph. Res., Vol. 111, B02406, doi:10.1029/2005JB003629.
#
#   Please mind that the coefficients in this paper are wrong. The corrected version of
#   the paper can be found at:
#   http://vmf.geo.tuwien.ac.at/documentation/Boehm#20et#20al.,#202006a#20#20#20(VMF1).pdf
#
#
#    input data
#    ----------
#    ah:   hydrostatic coefficient a (http://vmf.geo.tuwien.ac.at/trop_products/)
#    aw:   wet coefficient a         (http://vmf.geo.tuwien.ac.at/trop_products/)
#    dmjd: modified julian date
#    dlat: ellipsoidal latitude in radians
#    dlon: ellipsoidal longitude in radians
#    ht:   ellipsoidal height in meter
#    zd:   zenith distance in radians
#
#    output data
#    -----------
#    vmf1h: hydrostatic mapping function
#    vmf1w: wet mapping function
#
#    Johannes Boehm, 2005 October 2
#    rev 2011 July 21: latitude -> ellipsoidal latitude
#
 
#    implicit double precision (a-h,o-z)    
#    pi = 3.14159265359e0  
#    reference day is 28 January
#    this is taken from Niell (1996) to be consistent
#    File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, 2023-02-20
# =============================================================================
    # ell=[np.degrees(dlat),np.degrees(dlon),ht]
    # ah,aw,zhd,zwd=read_vmf1_grid(dmjd,ell)   
    
    doy = dmjd - 44239.0 + 1 - 28
    bh = 0.0029
    c0h = 0.062
    if (dlat < 0):
        phh = np.pi
        c11h = 0.007
        c10h = 0.002
    else:
        phh = 0.0
        c11h = 0.005
        c10h = 0.001
    
    ch = c0h + ((np.cos(doy / 365.25 * 2.0 * np.pi + phh) + 1.0) * c11h / 2.0 + c10h) * (1.0 - np.cos(dlat))
    sine = np.sin(np.pi / 2.0 - zd)
    beta = bh / (sine + ch)
    gamma = ah / (sine + beta)
    topcon = (1.0 + ah / (1.0 + bh / (1.0 + ch)))
    vmf1h = topcon / (sine + gamma)
    bw = 0.00146
    cw = 0.04391
    beta = bw / (sine + cw)
    gamma = aw / (sine + beta)
    topcon = (1.0 + aw / (1.0 + bw / (1.0 + cw)))
    vmf1w = topcon / (sine + gamma)
    return vmf1h,vmf1w
VMF.py
zohrehadavi/STD_SWD_Calc
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""

#Name_Datasource: e.g: week<2238  igr & week >= 2238 GRG0OPSFIN_
from ftplib import FTP_TLS
import datetime
from datetime import date
import gnsscal
import gzip
import os
import unlzw3
from pathlib import Path
import csv
import time
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings("ignore")

# File created by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================

def merge_files(files_to_merge, output_file):
    merged_data = pd.DataFrame()

# Loop through the files and merge them
    for file in files_to_merge:
    # Read the data from the file
       if file !=files_to_merge[len(files_to_merge)-1]:
           data = pd.read_csv(file,header = None, index_col= False, skip_blank_lines=True,skipfooter = 1)
       else:
           data = pd.read_csv(file,header = None, index_col= False, skip_blank_lines=True,skipfooter = 0)

    # Remove the header from all files except the first one
       if file != files_to_merge[0]:
           data0 = data.loc[22:,0]
       else:
           data0 = data
           
    
    # Merge the data into the merged_data dataframe
           merged_data = pd.concat([merged_data, data0])
           del data,data0
    merged_data.to_csv(output_file, header = None,index=False)
    
###############################################################################
               
def Download_Ultra_Orbit(year=None,month=None,day=None,folder=None,email=None,Name_Datasource=None,Hour=None):

    
    if (Name_Datasource=='igv'):
       try:
          fle_destinations=[]
          if Hour==None:
             Hour=[0,6,12,18,0]
             for it in range(0,len(Hour)):
                if it!=len(Hour)-1:
                   DoY=date(year, month, day).timetuple().tm_yday
                   dayt = datetime.date(year, month, day)
                   D=gnsscal.date2gpswd(dayt)
                   destination = folder +'/'+ "IGR" +str(D[0])+str(D[1])+'.SP3'
                else:
                   DoY=date(year, month, day+1).timetuple().tm_yday
                   dayt = datetime.date(year, month, day+1)
                   D=gnsscal.date2gpswd(dayt)
                   
                directory = "/glonass/products/"+str(D[0])
                file_name =Name_Datasource.lower()+str(D[0])+str(D[1])+'_'+str(Hour[it]).zfill(2)+'.sp3.Z'
                ftps = FTP_TLS(host = 'gdc.cddis.eosdis.nasa.gov')
                ftps.login(user='anonymous', passwd=email)
                ftps.prot_p()
                ftps.nlst()
                ftps.cwd(directory)
                ftps.retrbinary("RETR " + file_name, open(file_name, 'wb').write)
                ftps.quit()
             #####
                uncompressed_data = unlzw3.unlzw(Path(file_name))
                data = uncompressed_data.decode('utf-8').splitlines()
                with open(file_name[0:15], "w", newline="") as csv_file:
                 writer = csv.writer(csv_file)
                 for line in data:
                     # making cells by using , delimiters
                     writer.writerow(line.split(","))
             ######
                destination0 = folder +'/'+ "IGR" +str(D[0])+str(D[1])+'_'+str(Hour[it]).zfill(2)+'.SP3'
                source = folder +'/'+ file_name[0:15]
                os.rename(source, destination0)
                os.remove(file_name)
                fle_destinations.append(destination0)
             
            
             merge_files(fle_destinations, destination)
             #remove downloaded files
             for filename in fle_destinations:
       
                 file_path = filename
                 if os.path.isfile(file_path):
                         os.remove(file_path)
          else:
              
             DoY=date(year, month, day+1).timetuple().tm_yday
             dayt = datetime.date(year, month, day+1)
             D=gnsscal.date2gpswd(dayt)
             
             destination = folder +'/'+ "IGR" +str(D[0])+str(D[1])+'.SP3'
             directory = "/glonass/products/"+str(D[0])
             file_name =Name_Datasource.lower()+str(D[0])+str(D[1])+'_'+str(Hour).zfill(2)+'.sp3.Z'
             ftps = FTP_TLS(host = 'gdc.cddis.eosdis.nasa.gov')
             ftps.login(user='anonymous', passwd=email)
             ftps.prot_p()
             ftps.nlst()
             ftps.cwd(directory)
             ftps.retrbinary("RETR " + file_name, open(file_name, 'wb').write)
             ftps.quit()
          #####
             uncompressed_data = unlzw3.unlzw(Path(file_name))
             data = uncompressed_data.decode('utf-8').splitlines()
             with open(file_name[0:15], "w", newline="") as csv_file:
              writer = csv.writer(csv_file)
              for line in data:
                  # making cells by using , delimiters
                  writer.writerow(line.split(","))
          ######
             source = folder +'/'+ file_name[0:15]
             os.rename(source, destination)
             os.remove(file_name)
       except:
         pass

    else:
        try:
           directory = "/gnss/products/"+str(D[0])
         
           if Name_Datasource=='IGS':
               file_name=Name_Datasource+'0OPSULT_'+str(year)+str(DoY)+'0000_02D_15M_ORB.SP3.gz'
           else:
               file_name=Name_Datasource+'0OPSULT_'+str(year)+str(DoY)+'0000_02D_05M_ORB.SP3.gz'
           
           ftps = FTP_TLS(host = 'gdc.cddis.eosdis.nasa.gov')
           ftps.login(user='anonymous', passwd=email)
           ftps.prot_p()
           ftps.nlst()
           ftps.cwd(directory)
           ftps.retrbinary("RETR " + file_name, open(file_name, 'wb').write)
           ftps.quit()
           #####
           input = gzip.GzipFile(folder+'/'+file_name, 'rb')
           s= input.read()
           input.close()

           output = open(folder+'/'+file_name[0:38], 'wb')
           output.write(s)
           output.close()
           ######
           destination = folder + "IGR" +str(D[0])+str(D[1])+'.SP3'
           source = folder +'/'+ file_name[0:38]
           os.rename(source, destination)
           
        except:
          pass
    return destination
################################################################################
import os
import shutil

def move_and_rename_file(source_file, destination_folder, new_filename):
    # Move the file to the destination folder
    shutil.move(source_file, destination_folder)

    # Rename the file
    old_filename = os.path.basename(source_file)
    old_path = os.path.join(destination_folder, old_filename)
    new_path = os.path.join(destination_folder, new_filename)
    os.rename(old_path, new_path)
################################################################################
#Type_orb='Final' or 'Ultra'
def Int_Orbit(year=None,month=None,day=None,folder=None,email=None,Name_Datasource=None,Hour=None,Type_orb=None,Type_data_source=None,new_int_obs=None,Sat_cons=None):
 
 if Type_data_source=='MGEX':
      destination=Download_MGEX_Orbit(year,month,day,folder,email,Name_Datasource,Type_orb)
 else:
    if (Type_orb=='Final'):
       destination=Download_Final_Orbit(year,month,day,folder,email,Name_Datasource)
    else:
       destination=Download_Ultra_Orbit(year,month,day,folder,email,Name_Datasource,Hour)
                  
 Columnsmain0=['SAT','Year','Month','Day','hour','minute','second','X','Y','Z']
 line0='/* PCV:IGS'
 with open(destination) as f:
        cur_line = f.readline()
        cur_line = f.readline()
        int_obs=int(cur_line[25:29])
        cur_line = f.readline()
        NUMSAT=int(cur_line[4:7])
        f.close()
        
 with open(destination) as f:
        pos = 0
        cur_line = f.readline()
        
        while not cur_line.startswith(line0):
              pos = pos+1
              cur_line = f.readline()
        f.close()
        
 if ((Name_Datasource=='GFZ') or (Name_Datasource=='JAX')):
     pos=pos+3
 elif ((Name_Datasource=='WUM') and (Type_orb!='Final')):
    pos=pos+0
 elif ((Name_Datasource=='WUM') and (Type_orb=='Final')):
     pos=pos+2
 elif (Name_Datasource=='GRG'):
     pos=pos+2
 else:
     pos=pos+0
  
 if ((Name_Datasource=='igv') and (Name_Datasource=='CODE_orbits')):
    const_footer=0
 else:
    const_footer=1
    
 df = pd.read_csv(destination ,skiprows=[i for i in range(0,pos+1)],sep='\s+',skipfooter = const_footer, skip_blank_lines=True,header = None, index_col= False,names=['a','b','c','d','e','f','g','h','i','j','k'] )

 Mainfle0=pd.DataFrame(np.zeros((0,len(Columnsmain0))),columns=Columnsmain0)
 kt=0

 for it in range(0,len(df),NUMSAT+1):
    time_acq=df.iloc[it,0:7]
    for jt in range(it+1,it+NUMSAT+1):
        if jt<len(df):
           Mainfle0.loc[kt,'SAT']=df.iloc[jt,0][1:5]
           Mainfle0.loc[kt,'Year']=time_acq[1]
           Mainfle0.loc[kt,'Month']=time_acq[2]
           Mainfle0.loc[kt,'Day']=time_acq[3]
           Mainfle0.loc[kt,'hour']=time_acq[4]
           Mainfle0.loc[kt,'minute']=time_acq[5]
           Mainfle0.loc[kt,'second']=time_acq[6]
           Mainfle0.loc[kt,'X']=df.iloc[jt,1]
           Mainfle0.loc[kt,'Y']=df.iloc[jt,2]
           Mainfle0.loc[kt,'Z']=df.iloc[jt,3]
           kt=kt+1

 del df
 os.remove(destination) 

################ 
 Mainfle=Mainfle0.where(Mainfle0['Day'] == day)
 Mainfle = Mainfle.dropna()
 Mainfle.reset_index(drop = True, inplace = True)
 
 Mainfle.drop_duplicates(subset=['SAT','Day','hour','minute','second'], keep='last', inplace=True)
 Mainfle.reset_index(drop = True, inplace = True)
##########Interpolation for every 2.5 minutes
 end_time=int(int(Mainfle.iloc[-1,4])*3600+int(Mainfle.iloc[-1,5])*60+int(float(Mainfle.iloc[-1,6])))
 fitTime = range(0,end_time+int_obs,int_obs)



 sorted_Mainfle = Mainfle.sort_values(by=['SAT','hour','minute','second'], ascending=[True,True,True,True])
 sorted_Mainfle.reset_index(drop = True, inplace = True)

  #2.5 minutes
 interval = range(0,end_time+new_int_obs,new_int_obs)


 Int_fle=pd.DataFrame(np.zeros((0,len(Columnsmain0))),columns=Columnsmain0) 
 nt=0 
 mt=0 
 for kt in range(0,NUMSAT):
    
    df_new=sorted_Mainfle.loc[mt:len(fitTime)+mt-1,:]
    df_new.reset_index(drop = True, inplace = True)
    fitx=np.polyfit(fitTime,sorted_Mainfle.loc[mt:mt+len(fitTime)-1,'X'].astype('float'),deg=16)
    fity=np.polyfit(fitTime,sorted_Mainfle.loc[mt:mt+len(fitTime)-1,'Y'].astype('float'),deg=16)
    fitz=np.polyfit(fitTime,sorted_Mainfle.loc[mt:mt+len(fitTime)-1,'Z'].astype('float'),deg=16)
    
    x_interp = np.polyval(fitx, interval)
    y_interp = np.polyval(fity, interval)
    z_interp = np.polyval(fitz, interval)
    
    step_fix=int(int_obs/new_int_obs)
    lt=0
    def_xs=[]
    for jt in range(0,len(x_interp),step_fix):
        x_interp[jt]=df_new.loc[lt,'X']
        y_interp[jt]=df_new.loc[lt,'Y']
        z_interp[jt]=df_new.loc[lt,'Z']
        lt=lt+1
        
    for it in range(0,len(x_interp)):
           convert = time.strftime("%H:%M:%S", time.gmtime(interval[it]))
           Int_fle.loc[nt,'SAT']=df_new.iloc[0,0]
           Int_fle.loc[nt,'Year']=df_new.iloc[0,1]
           Int_fle.loc[nt,'Month']=df_new.iloc[0,2]
           Int_fle.loc[nt,'Day']=df_new.iloc[0,3]
           Int_fle.loc[nt,'hour']=convert[0:2]
           Int_fle.loc[nt,'minute']=convert[3:5]
           Int_fle.loc[nt,'second']=convert[6:8]
           Int_fle.loc[nt,'X']=x_interp[it]
           Int_fle.loc[nt,'Y']=y_interp[it]
           Int_fle.loc[nt,'Z']=z_interp[it]  
        
           nt=nt+1
    
    mt=mt+len(fitTime)
 ############### 
    Int_fle_sat=[]
    for jt in range(0,len(Sat_cons)):
     df_sat=Int_fle.where(Int_fle['SAT'].astype(str).str[0]==Sat_cons[jt])
     Int_fle_sat.append(df_sat)
     del df_sat
  
 Int_fle_final0=pd.concat(Int_fle_sat)
 Int_fle_final = Int_fle_final0.dropna()
 Int_fle_final.reset_index(drop = True, inplace = True) 
 del Mainfle,Mainfle0,sorted_Mainfle   
 return Int_fle_final

##############################################################################
import pyproj
transformer = pyproj.Transformer.from_crs(
    {"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
    {"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
    )

def Orbit_ELV_AZ(year=None,month=None,day=None,folder=None,email=None,Name_Datasource=None,Hour=None,Xo=None,Yo=None,Zo=None,elv=None,Type_orb=None,Type_data_source=None,new_int_obs=None,Sat_cons=None):

    Columnsmain=['SAT','Year','Month','Day','hour','minute','second','X','Y','Z','ELV','AZ']
    df0=Int_Orbit(year,month,day,folder,email,Name_Datasource,Hour,Type_orb,Type_data_source,new_int_obs,Sat_cons)
    lambdaO, phiO, h = transformer.transform(Xo,Yo,Zo,radians=True) 
    Xs=df0.loc[:,'X']*1000
    Ys=df0.loc[:,'Y']*1000
    Zs=df0.loc[:,'Z']*1000
    ## *******************Transformation To ECEF*****************************
    r = np.sqrt((Xs - Xo) ** 2 + (Ys - Yo) ** 2 + (Zs - Zo) ** 2)
    dx = (Xs - Xo) / r
    dy = (Ys - Yo) / r
    dz = (Zs - Zo) / r
    E = (np.multiply((- np.sin(lambdaO)),(dx))) + (np.multiply((np.cos(lambdaO)),(dy)))
    N = (np.multiply((np.multiply(- np.sin(phiO),np.cos(lambdaO))),(dx))) + (np.multiply(np.multiply(- np.sin(phiO),np.sin(lambdaO)),(dy))) + (np.multiply(np.cos(phiO),(dz)))
    U = (np.multiply((np.multiply(np.cos(phiO),np.cos(lambdaO))),(dx))) + (np.multiply((np.multiply(np.cos(phiO),np.sin(lambdaO))),(dy))) + (np.multiply((np.sin(phiO)),(dz)))
    ## ********************Elevation and Azimuth Computation*****************
    elev = np.degrees(np.arctan(U / (np.sqrt((E) ** 2 + (N) ** 2))))
    az = np.degrees(np.arctan2(E,N))
    
    df_new=pd.DataFrame(np.zeros((0,len(Columnsmain))),columns=Columnsmain) 

    
    for it in range(0,len(elev)):
        if (elev[it]>=elv) and (az[it]>=0):
            df_new.loc[it,'SAT']=df0.iloc[it,0]
            df_new.loc[it,'Year']=df0.iloc[it,1]
            df_new.loc[it,'Month']=df0.iloc[it,2]
            df_new.loc[it,'Day']=df0.iloc[it,3]
            df_new.loc[it,'hour']=df0.iloc[it,4]
            df_new.loc[it,'minute']=df0.iloc[it,5]
            df_new.loc[it,'second']=df0.iloc[it,6]
            df_new.loc[it,'X']=df0.iloc[it,7]
            df_new.loc[it,'Y']=df0.iloc[it,8]
            df_new.loc[it,'Z']=df0.iloc[it,9]
            df_new.loc[it,'ELV']=elev[it] 
            df_new.loc[it,'AZ']=az[it] 
            
        elif (elev[it]>=elv) and (az[it]<0):

           df_new.loc[it,'SAT']=df0.iloc[it,0]
           df_new.loc[it,'Year']=df0.iloc[it,1]
           df_new.loc[it,'Month']=df0.iloc[it,2]
           df_new.loc[it,'Day']=df0.iloc[it,3]
           df_new.loc[it,'hour']=df0.iloc[it,4]
           df_new.loc[it,'minute']=df0.iloc[it,5]
           df_new.loc[it,'second']=df0.iloc[it,6]
           df_new.loc[it,'X']=df0.iloc[it,7]
           df_new.loc[it,'Y']=df0.iloc[it,8]
           df_new.loc[it,'Z']=df0.iloc[it,9]
           df_new.loc[it,'ELV']=elev[it] 
           df_new.loc[it,'AZ']=az[it]+360
           
        else:
            df_new.loc[it,'SAT']=df0.iloc[it,0]
            df_new.loc[it,'Year']=df0.iloc[it,1]
            df_new.loc[it,'Month']=df0.iloc[it,2]
            df_new.loc[it,'Day']=df0.iloc[it,3]
            df_new.loc[it,'hour']=df0.iloc[it,4]
            df_new.loc[it,'minute']=df0.iloc[it,5]
            df_new.loc[it,'second']=df0.iloc[it,6]
            df_new.loc[it,'X']=df0.iloc[it,7]
            df_new.loc[it,'Y']=df0.iloc[it,8]
            df_new.loc[it,'Z']=df0.iloc[it,9]
            df_new.loc[it,'ELV']=0 
            df_new.loc[it,'AZ']=0
    
    
    return df_new
satellite_read_int_elv_az.py
zohrehadavi/STD_SWD_Calc
# -*- coding: utf-8 -*-
"""
@author: zadavi
"""

from Read_Write_TRP import Readwrite_TRO_ZTD
from satellite_read_int_elv_az import Orbit_ELV_AZ
from datetime import date
from GPT3 import gpt3_5_fast
from VMF import read_vmf1_grid,vmf1_ht
from GMF_GradinetMF import gmf,gradient_MF
from CalcSwd import CalcSWDgrad

from satellite_read_int_elv_az import *
from Read_Write_TRP import *
from GPT3 import *
from VMF import *
from GMF_GradinetMF  import *
from CalcSwd import *

import pandas as pd
import numpy as np
from astropy.time import Time
import pygeodesy
from pygeodesy.ellipsoidalKarney import LatLon
ginterpolator = pygeodesy.GeoidKarney("./geoids/egm2008-5.pgm")


def write_Sinex_SWD(year=None,month=None,day=None,folder_orb=None,email=None,Name_Datasource_orb=None,Type_orb=None,Hour=0,folder_source_tro=None,file_name_tro=None,folder_dest_snx=None,fle_out_snx_ztd=None,fle_out_snx_swd=None,Type_MF=None,Type_zhd=None,elv=None,Type_data_source=None,New_rate=None,Sat_cons=None,int_need_TRO=None,grad=None):
#     Input parameters:
#     -------------------    
#    folder_orb: The saving path for orbit file
#    email: Registered email for CDDIS
#    Name_Datasource_orb: IGS, GFZ, CODE, GRG, ESA 'For ultra rapid before week 2038 please enter 'igv' if the 'Type_data_source' is not 'MGEX'
#    Type_orb: Final or Ultra
#    Hour: In the case of Ultra, the hour can be selected; otherwise, the orbit file for hour 0,6,12,18 will download and then merge them (it applies for ultra rapid file 'igv')
#    file_source_tro: The path of TRO files
#    file_name_tro: The name of TRO file
#    file_dest_snx: Destination of created ztd or swd SINEX files
#    fle_out_snx_ztd: File name for SINEX ztd
#    fle_out_snx_swd: File name for SINEX swd
#     Type_MF:
#     GMF or VMF1 mapping function
#     Type_zhd: 
#     zhd_GPT: ZHD calculated from GPT model and Saastamoinen 
#     zhd_VMF: ZHD derived from VMF1 grid
#    elv: Cut of angle (defult:5 [deg])
#    New_rate: new rate for interpolating data (ORB and ZTD) and creating finer dataset (defult: 150 [sec] #2.5 minutes)
#    Type_data_source: If you want to download from MGEX please set option as 'MGEX', otherwise it will download from the main directory of CDDIS
#    Sat_cons        : Which GNSS constellation should be used (IT must be consistent with the ORB data source (for instance ['G','R','E'])) 
#    int_need_TRO: 1 if interpolation is needed and 0 otherwise (default=0)
#    File created by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================

  columns_snx_swd=['SITE','lon', 'lat', 'h', 'Xsta', 'Ysta', 'Zsta', 'Xsat', 'Ysat', 'Zsat', 'PRN', 'AZ', 'ELV', 'hour', 'minute', 'second','mjd', 'SWD','ZTD']
  doy=date(year, month, day).timetuple().tm_yday
  
  TRO=Readwrite_TRO_ZTD(folder_source_tro,file_name_tro,folder_dest_snx,fle_out_snx_ztd,doy,int_need_TRO,New_rate)
  Sta_name=TRO['SITE'].unique()
  
  snx_swd=pd.DataFrame(np.zeros((0,len(columns_snx_swd))),columns=columns_snx_swd)
  df_final=[]
  
  for it in range(0,len(Sta_name)):
    jloc0=TRO.loc[(TRO['SITE']==Sta_name[it])]
    jloc0.reset_index(drop=True, inplace=True)
    
    lon=jloc0.iloc[0,1]
    lat=jloc0.iloc[0,2]
    h=jloc0.iloc[0,3]
    Xo=jloc0.iloc[0,4]
    Yo=jloc0.iloc[0,5]
    Zo=jloc0.iloc[0,6]
    

    
    
    single_position=LatLon(lat, lon)
    Ngeo= ginterpolator(single_position)
    h_ortho=h-Ngeo
    
    ORB=Orbit_ELV_AZ(year,month,day,folder_orb,email,Name_Datasource_orb,Hour,float(Xo),float(Yo),float(Zo),elv,Type_orb,Type_data_source,New_rate,Sat_cons)

    for jt in range(0,len(jloc0)):
      ZTD=jloc0.iloc[jt,10]/1000
      ge=jloc0.iloc[jt,11]/1000
      gn=jloc0.iloc[jt,12]/1000
        
      jloc1=ORB.loc[(ORB['hour'] ==jloc0.iloc[jt,7]) & (ORB['minute'] ==jloc0.iloc[jt,8]) & (ORB['second'] == jloc0.iloc[jt,9]) & ((ORB['ELV'])>elv)]
      jloc1.reset_index(drop=True, inplace=True)

      if len(jloc1)>0:
        t = Time('{}-{}-{}T{}:{}:{}'.format(year,month, day, int(jloc0.iloc[jt,7]), int(jloc0.iloc[jt,8]), int(jloc0.iloc[jt,9])), format='isot', scale='utc')
        mjd=t.mjd
         
        ELV=jloc1.loc[:,'ELV']
        AZ =jloc1.loc[:,'AZ']
        SWD,ZWD,mf_w,MF_G=CalcSWDgrad(mjd,float(ZTD),np.radians(float(lat)),np.radians(float(lon)),float(h),h_ortho,float(gn),float(ge),AZ,ELV,Type_MF,Type_zhd,grad)
        df_new=pd.DataFrame(np.zeros((len(ELV),len(columns_snx_swd))),columns=columns_snx_swd)  
         
        df_new.loc[:,'SITE']=jloc0.iloc[0,0]
        df_new.loc[:,'lon']=lon
        df_new.loc[:,'lat']=lat
        df_new.loc[:,'h']=h
         
        df_new.loc[:,'Xsta']=float(Xo)
        df_new.loc[:,'Ysta']=float(Yo)
        df_new.loc[:,'Zsta']=float(Zo)
         
        df_new.loc[:,'Xsat']=jloc1.loc[:,'X']
        df_new.loc[:,'Ysat']=jloc1.loc[:,'Y']
        df_new.loc[:,'Zsat']=jloc1.loc[:,'Z']
            
        df_new.loc[:,'PRN']=jloc1.loc[:,'SAT']
        df_new.loc[:,'AZ']=jloc1.loc[:,'AZ']
        df_new.loc[:,'ELV']=jloc1.loc[:,'ELV']
         
        df_new.loc[:,'hour']=jloc0.iloc[jt,7]
        df_new.loc[:,'minute']=jloc0.iloc[jt,8]
        df_new.loc[:,'second']=jloc0.iloc[jt,9]
        df_new.loc[:,'mjd']=mjd
        
        df_new.loc[:,'SWD']=np.transpose(SWD)
        df_new.loc[:,'ZTD']=float(ZTD)
        df_final.append(df_new) 
        del df_new
    del jloc1
  del jloc0,ORB
  snx_swd=pd.concat(df_final)     
#########################################################
  output_path=folder_dest_snx+'\\'+fle_out_snx_swd
# finaldf.to_csv(output_path, columns = False, index = False)
  with open(output_path, 'a') as f:
        df_string = snx_swd.to_string(header=False, index=False)
        f.write(df_string)
        
  return snx_swd
        
save_Sinex_SWD_ZTD.py
zohrehadavi/STD_SWD_Calc