added more refraction models

.gitignore

@@ -1,4 +1,5 @@
 env
+gnssrefl/gnssir_old.py
 list_bfg.py
 noaa.txt
 idalia.py

CHANGELOG.md

@@ -3,6 +3,12 @@
 All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
+## 1.8.4
+
+Added refraction model number (refr_model) to gnssir_input and thus station analysis json.
+on the way to implementing more than one refraction model. Default will remain refr_model = 1.
+0 is no model.
+
 ## 1.8.3
 2023 September 26
 

README.md

@@ -1,9 +1,11 @@
 # gnssrefl
 
-**github version: 1.8.3** [![PyPI Version](https://img.shields.io/pypi/v/gnssrefl.svg)](https://pypi.python.org/pypi/gnssrefl) [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.5601495.svg)](http://dx.doi.org/10.5281/zenodo.5601495) [![Documentation Status](https://readthedocs.org/projects/gnssrefl/badge/?version=latest)](https://gnssrefl.readthedocs.io/en/latest/?badge=latest)
+**github version: 1.8.4** [![PyPI Version](https://img.shields.io/pypi/v/gnssrefl.svg)](https://pypi.python.org/pypi/gnssrefl) [![DOI](https://zenodo.org/badge/doi/10.5281/zenodo.5601495.svg)](http://dx.doi.org/10.5281/zenodo.5601495) [![Documentation Status](https://readthedocs.org/projects/gnssrefl/badge/?version=latest)](https://gnssrefl.readthedocs.io/en/latest/?badge=latest)
 
 Our documentation is available [here.](https://gnssrefl.readthedocs.io/en/latest/)
 
+See documentation for gnssir_input for new refraction models.
+
 A long time ago - and in a galaxy far, far away - Bob King told a group of graduate students why he 
 liked GPS observation files. His answer:
 

docs/pages/understand.md

@@ -58,7 +58,7 @@ center above the reflecting surface, or reflector height RH (purple). *The prima
 is to measure RH.* This parameter is directly related to changes in snow height and water levels below
 a GNSS antenna. This is why GNSS-IR can be used as a snow sensor and tide gauge. GNSS-IR can also be 
 used to measure soil moisture, but the code to estimate soil moisture is not as strongly related to RH as
-snow and water. We will be posting the code you need to measure soil moisture later in the year.
+snow and water. 
 
 ![](https://gnss-reflections.org/static/images/overview.png)
 

gnssrefl/gnssir_cl.py

@@ -268,14 +268,16 @@ def gnssir(station: str, year: int, doy: int, snr: int = 66, plt: bool = False,
     # added 2022apr15
     lsp['gzip'] = gzip
 
+    # if refraction model is not assigned, set it to 1
+    if 'refr_model' not in lsp.keys():
+        lsp['refr_model'] = 1
 
     xdir = str(os.environ['REFL_CODE'])
     picklefile = 'gpt_1wA.pickle'
     pname = xdir + '/input/' + picklefile
 
-    if os.path.isfile(pname):
-        print('refraction file exists')
-    else:
+    #    print('refraction file exists')
+    if not os.path.isfile(pname):
         local_copy = 'gnssrefl/' + picklefile
         if os.path.isfile(local_copy):
             print('found local copy of refraction file')

gnssrefl/gnssir_input.py

@@ -36,9 +36,7 @@ def parse_arguments():
     parser.add_argument("-frlist", nargs="*",type=int,  help="User defined frequencies using our nomenclature")
     parser.add_argument("-azlist2", nargs="*",type=float,  default=None,help="list of azimuth regions, default is 0-360") 
     parser.add_argument("-ellist", nargs="*",type=float,  default=None,help="List of elevation angles to allow more complex analysis scenarios-advanced users only!") 
-
-    #print(parser.parse_args())
-
+    parser.add_argument("-refr_model", default=1, type=int, help="refraction model. default is 1, zero turns it off)")
 
 
     args = parser.parse_args().__dict__
@@ -52,10 +50,11 @@ def parse_arguments():
 
 
 def make_gnssir_input(station: str, lat: float=0, lon: float=0, height: float=0, e1: float = 5.0, e2: float = 25.0,
-              h1: float = 0.5, h2: float = 8.0, nr1: float = None, nr2: float = None,
-              peak2noise: float = 2.8, ampl: float = 5.0, allfreq: bool = False,
-              l1: bool = False, l2c: bool = False, xyz: bool = False, refraction: bool = True,
-              extension: str = '', ediff: float=2.0, delTmax: float=75.0, frlist: list=[],azlist2: list=[0,360], ellist : list=[] ):
+       h1: float = 0.5, h2: float = 8.0, nr1: float = None, nr2: float = None, peak2noise: float = 2.8, 
+       ampl: float = 5.0, allfreq: bool = False, l1: bool = False, l2c: bool = False, 
+       xyz: bool = False, refraction: bool = True, extension: str = '', ediff: float=2.0, 
+       delTmax: float=75.0, frlist: list=[], azlist2: list=[0,360], 
+       ellist : list=[], refr_model : int=1 ):
 
     """
     This new script sets the Lomb Scargle analysis strategy you will use in gnssir. It saves your inputs 
@@ -76,6 +75,24 @@ def make_gnssir_input(station: str, lat: float=0, lon: float=0, height: float=0,
     Example of a json file that has both the old and new ways. The arrays within the json file are called azval and azval2.
         https://morefunwithgps.com/public_html/sc02.json
 
+    Refraction  
+    ----------
+        originally we had refraction as a boolean, i.e. on or off. This is stored in the json.  The code however,
+        uses an integer 1 (for a simple non-time-varying Bennett correction) and integer 0 for no correction.
+        From version 1.8.4 we begin to implement more refraction models.  1 (and Bennett) will continue to be 
+        the default.  The "1" is written to the LSP results file so that people can keep track easily of whether
+        they are inadvertently mixing files with different strategies. And that is why it is an integer, because
+        all results in the LSP results files are numbers.  Going forward, we are adding a time-varying capability.
+
+            Model 1: Bennett, static
+            Model 2: Bennett and time-varying
+            Model 3: Ulich, static
+            Model 4: Ulich, time-varying
+
+            gnssir_input will have a new parameter for the json output, refr_model. If it is not set, i.e. you 
+            have an old json, it is assumed to be 1. You can change it be hand if you like. And you can certainly 
+            test out the impact by using -extension option.
+
     Examples
     --------
     gnssir_input p041 
@@ -163,6 +180,11 @@ def make_gnssir_input(station: str, lat: float=0, lon: float=0, height: float=0,
         min and max elevation angles to be used with the azimuth regions you listed, i.e.
         [5 10 6 11 7 12 8 13] would allow overlapping regions - all five degrees long 
         Default is empty list. 
+    refr_model : int
+        refraction model. we are keeping this as integer as it is written to a file withonly
+        numbers in it.  1 is the default simple refraction (just correct elevation angles
+        using standard bending models).  0 is no refraction correction.  As we add more
+        models, they will receiver their own number. 
 
     """
 
@@ -277,6 +299,7 @@ def make_gnssir_input(station: str, lat: float=0, lon: float=0, height: float=0,
     # but the length of the list depends on the length of the list of frequencies
     lsp['reqAmp'] = [reqA]*len(lsp['freqs'])
 
+    # this is true or false.  I think
     lsp['refraction'] = refraction
 
     # write new RH results  each time you run the code
@@ -309,14 +332,16 @@ def make_gnssir_input(station: str, lat: float=0, lon: float=0, height: float=0,
 
     lsp['ellist'] = ellist
 
+    lsp['refr_model'] = refr_model
+
     print('writing out to:', outputfile)
+    print(lsp)
     with open(outputfile, 'w+') as outfile:
         json.dump(lsp, outfile, indent=4)
 
 
 def main():
     args = parse_arguments()
-    print(args)
     make_gnssir_input(**args)
 
 

gnssrefl/gnssir_v2.py

@@ -129,8 +129,14 @@ def gnssir_guts_v2(station,year,doy, snr_type, extension,lsp):
    # this defines the minimum number of points in an arc.  This depends entirely on the sampling
    # rate for the receiver, so you should not assume this value is relevant to your case.
     minNumPts = 20
-    p,T,irefr = set_refraction_params(station, dmjd, lsp)
-    #print('Refraction parameters ',p,T,irefr)
+    p,T,irefr,humidity = set_refraction_params(station, dmjd, lsp)
+    print('Refraction parameters (press, temp, humidity, ModelNum',np.round(p,3),np.round(T,3),np.round(humidity,3),irefr)
+
+    if (irefr == 3) or (irefr == 4):
+        TempK = T + 273.15 # T is in celsius ... I think.
+        # N_ant is the refraction index
+        N_ant=refr.refrc_Rueger(p,humidity,TempK)[0]    #
+
 
 # only doing one day at a time for now - but have started defining the needed inputs for using it
     twoDays = False
@@ -163,10 +169,23 @@ def gnssir_guts_v2(station,year,doy, snr_type, extension,lsp):
                 print('the minimum elevation angle your receiver used. Exiting.')
                 sys.exit()
 
-
-        # apply simple refraction correction
-        ele = snrD[:,1]
-        ele=apply_refraction_corr(lsp,ele,p,T)
+        if (irefr == 3) or (irefr == 4):
+            # elev refraction, lsp, pres, temp, time, sat
+            if irefr == 3:
+                print('Ulich refraction correction')
+            else:
+                print('Ulich refraction correction, time-varying')
+            ele=refr.Ulich_Bending_Angle(snrD[:,1],N_ant,lsp,p,T,snrD[:,3],snrD[:,0])
+        elif irefr == 0:
+            print('No refraction correction applied ')
+            ele = snrD[:,1]
+        else :
+            if irefr == 1:
+                print('Standard Bennett refraction correction')
+            else:
+                print('Standard Bennett refraction correction, time-varying')
+            ele = snrD[:,1]
+            ele=apply_refraction_corr(lsp,ele,p,T)
 
         # apply an elevation mask for all the data for the polynomial fit
         i= (ele >= pele[0]) & (ele < pele[1])
@@ -336,29 +355,58 @@ def set_refraction_params(station, dmjd,lsp):
     dmjd : float
         modified julian date
 
-    lsp : dictionary
+    lsp : dictionary with information about the station
+        lat : float
+            latitude, deg
+        lon : float
+            longitude, deg
+        ht : float
+            height, ellipsoidal
+
+    Returns
+    -------
+    p : float
+        pressure, hPa
+    T : float
+        temperature, Celsius
+    irefr : int
+        refraction model number I believe, which is also sent, so not needed
+    e : float
+        water vapor pressure, hPa
 
     """
+    if 'refr_model' not in lsp.keys():
+        refraction_model = 1
+    else:
+        refraction_model = lsp['refr_model']
+
+
     xdir = os.environ['REFL_CODE']
     p = 0; T = 0; irefr = 0
     if lsp['refraction']:
-        irefr = 1
+        irefr = refraction_model
         refr.readWrite_gpt2_1w(xdir, station, lsp['lat'], lsp['lon'])
-# time varying is set to no for now (it = 1)
-        it = 1
+        # time varying was originally set to no for now (it = 1)
+        # now allows time varying for models 2 and 4
+        if (refraction_model == 2) or (refraction_model == 4):
+            it = 0
+        else:
+            it = 1
         dlat = lsp['lat']*np.pi/180; dlong = lsp['lon']*np.pi/180; ht = lsp['ht']
         p,T,dT,Tm,e,ah,aw,la,undu = refr.gpt2_1w(station, dmjd,dlat,dlong,ht,it)
+        # e is water vapor pressure, so humidity ??
         #print("Pressure {0:8.2f} Temperature {1:6.1f} \n".format(p,T))
 
-    return p,T,irefr
+    #print('refraction model', refraction_model,irefr)
+    return p,T,irefr, e
 
 def apply_refraction_corr(lsp,ele,p,T):
     """
 
     Parameters
     ----------
     lsp : dictionary
-        info from make_json_input
+        info from make_json_input such as station lat and lon
     ele : numpy array of floats
         elevation angles  (deg)
     p : float
@@ -756,7 +804,8 @@ def window_new(snrD, f, satNu,ncols,pele,pfitV,e1,e2,azlist,screenstats):
         ijk = (datatest == 0)
         nzero = len(datatest[ijk])
         if np.sum(datatest) < 1:
-            print('No useful data on frequency ', f , ': all zeros')
+            if screenstats:
+                print('No useful data on frequency ', f , 'and sat ', satNu, ': all zeros')
             good = False
         else:
             if nzero > 0:

gnssrefl/refraction.py

@@ -3,6 +3,7 @@
 from original TU Vienna codes for GMF
 """
 import datetime
+import math
 import os
 import pickle
 import subprocess
@@ -564,3 +565,83 @@ def look_for_pickle_file():
             print('File should be stored in ', inputdir, ' but is not')
 
     return foundit , fullpname
+
+
+def Ulich_Bending_Angle(ele, N0,lsp,p,T,ttime,sat):   #UBA
+    """
+    Ulich, B. L. "Millimeter wave radio telescopes: Gain and pointing characteristics." (1981)
+
+    Author: 20220629, fengpeng
+    modified to use numpy so I can do arrays
+
+    currently writes out corrections to a file for testing
+
+    Parameters
+    ----------
+    ele : numpy array of floats
+        true elevation angle, degrees
+
+    N0 : float
+        antenna refractivity in ppm
+
+    Returns
+    -------
+    De : numpy array of floats
+        corrected elevation angle, deg
+        
+    """
+    e_simple_corr = corr_el_angles(ele, p,T)
+    deg2rad = np.pi/180
+
+    # change to radians
+    ele_rad = deg2rad*(ele)
+    r = N0/1000000.
+    f = np.cos(ele_rad) / (np.sin(ele_rad) + 0.00175 * np.tan(deg2rad*(87.5) - ele_rad))
+    dE = (r * f)/deg2rad
+    fout = open('ulich.txt', 'w+')
+    for i in range(0,len(ele)):
+        fout.write(" {0:8.5f} {1:8.5f} {2:8.5f} {3:10.0f} {4:5.0f} \n ".format( 
+            ele[i], ele[i]+dE[i], e_simple_corr[i], ttime[i], sat[i]))
+
+    fout.close()
+
+    return dE + ele
+
+def refrc_Rueger(drypress,vpress,temp):
+    """
+    Obtains refractivity index suitable for GNSS-IR
+
+    Rüeger, Jean M. "Refractive index formulae for radio waves." Proceedings of the 
+    FIG XXII International Congress, Washington, DC, USA. Vol. 113. 2002.
+
+    Parameters
+    ----------
+    drypress : float
+        dry pressure hPa
+    vpress : float
+        vapor pressure in hPa
+    temp : float
+        temperature in Kelvin
+
+    Returns 
+    -------
+    ref : list of floats
+         [Ntotal, Nhydro, Nwet], which are total, hydrostatic and wet refractivity in ppm
+
+    """
+
+    [K1r,K2r,K3r]=[77.689 ,71.2952 ,375463.]     # Rüeger's "best average", for 375 ppm CO2, 77.690 for 392 ppm CO2
+
+    Nrueger = K1r * drypress / temp + K2r * vpress / temp + K3r * vpress / (temp ** 2)  # Rueger,IAG recommend, in ppm
+
+    [Rgas,Md,Mw] = [8.31446261815324,0.0289644,0.01801528]   # gas constant (SI exact), molar masses of dry air and water, kg/mol
+    drydensity=(drypress*100.)*Md/(Rgas*temp)    #density in kg/m^3
+    vdensity=(vpress*100.)*Mw/(Rgas*temp)   #in PV=nRT, P in Pa, V in m^3, n in g/mol, T in K, R=8.314 Pa*m^3/(k*mol)
+    totaldensity=drydensity+vdensity
+
+    Nhydro=(K1r*Rgas*totaldensity/Md)/100.   #divide 100 because hPa in Rueger formula, result in ppm
+    K2rr=K2r-K1r*(Mw/Md)
+    Nwet=K2rr * vpress / temp + K3r * vpress / (temp ** 2)
+    ref=[round(Nrueger,4),round(Nhydro,4),round(Nwet,4)]    # in ppm
+    return ref
+

gnssrefl/spline_functions.py

@@ -763,8 +763,8 @@ def snr2spline(station,year,doy, azilims, elvlims,rhlims, precision, kdt, snrfit
         if lsp['refraction']:
             imodel = 1 #
         dmjd = 59580 # fake number
-        p,T,irefr = set_refraction_model(station, dmjd,lsp,imodel)
-        print('refraction parameters ', p,T,irefr)
+        p,T,irefr,humidity = set_refraction_model(station, dmjd,lsp,imodel)
+        print('refraction parameters ', p,T,humidity,irefr)
 
     if 'rough_in' in kwargs:
         rough_in  = kwargs.get('rough_in')
@@ -1684,6 +1684,9 @@ def set_refraction_model(station, dmjd,lsp,imodel):
     imodel is 1 for simple refraction model
     eventually will add other refraction models
 
+    Looks like this was copied from other code and should be
+    consolidated ...
+
     Parameters
     ----------
     station : str
@@ -1703,6 +1706,8 @@ def set_refraction_model(station, dmjd,lsp,imodel):
         temperature in deg C
     irefr : int
         number value written to output files to keep track of refraction model
+    e : float
+        water vapor pressure, hPa
 
     """
     xdir = os.environ['REFL_CODE']
@@ -1717,5 +1722,5 @@ def set_refraction_model(station, dmjd,lsp,imodel):
         p,T,dT,Tm,e,ah,aw,la,undu = refr.gpt2_1w(station, dmjd,dlat,dlong,ht,it)
         #print("Pressure {0:8.2f} Temperature {1:6.1f} \n".format(p,T))
 
-    return p,T,irefr
+    return p,T,irefr, e