added documentation and options to daily_avg

CHANGELOG.md

@@ -3,10 +3,19 @@
 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.4.9
+August 10, 2023
+
+New option on daily_avg that plots the median value and the limits from that 
+median value that are used to filter outliers before the daily average is computed.
+
+More documentation about daily_avg has been added online.
+
+
 ## 1.4.8
 
 Another version to get it on pypi
-
 ## 1.4.7
 
 August 7, 2023

README.md

@@ -1,6 +1,6 @@
 # gnssrefl
 
-**github version: 1.4.8** [![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.4.9** [![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)
 
 August 7, 2023, again
 

docs/pages/README_dailyavg.md

@@ -18,6 +18,9 @@ varies a lot depending on the azimuth mask and the number of frequencies availab
 If you are not sure what values to use at your GNSS site, run it once with very minimal constraints.
 The code provides some feedback plots that will let you pick better values.
 
+Note: the computed daily average value should be associated with 12:00 UTC, not midnight.   
+
+
 Here is an example from one of our use cases where there are a few large outliers.  
 I have set the median filter value to 2 meters and the required number of tracks to 12:
 
@@ -29,6 +32,17 @@ You can easily see the outliers.
 <img width=500 src=../_static/mchn-A.png>
 </p>
 
+Is 12 a good choice?  The code also prints out a plot telling you how many
+tracks are available each day:
+
+<p align=center>
+<img width=500 src=../_static/mchn_nvals.png>
+</p>
+
+These can vary quite a bit by year as the station operators change receivers and/or 
+tracking strategies. You should pick the values that are best for your experiment.
+
+
 Next I have rerun the code with a better median filter constraint of 0.25 meters:
 
 <code> daily_avg mchn 0.25 12 </code>
@@ -37,9 +51,34 @@ Next I have rerun the code with a better median filter constraint of 0.25 meters
 <img width=500 src=../_static/mchn-B.png>
 </p>
 
-A daily average plot is also made and a text file of the outputs is created.
+Since this documentation was written, I have added the median value to the plots:
+
+<p align=center>
+<img width=500 src=../_static/mchn_tighter.png>
+</p>
+
+If you still are finding it challenging to see the variations from the median value, you
+can try setting the plot_limits option:
+
+
+<code> daily_avg mchn 0.25 12 -plot_limits T </code>
+
+<p align=center>
+<img width=500 src=../_static/mchn_wlimits.png>
+</p>
+
+
+The daily average plot:
 
 <p align=center>
 <img width=500 src=../_static/mchn-C.png>
 </p>
 
+Two txt files are created. One has all the tracks that fit the QC. The other has the desired daily
+average.  The location of the files is printed to the screen. You are welcome to generate your 
+own quality control codes for the daily average if you find this one does not meet your purposes.
+
+Updated August 10, 2023
+
+Kristine M. Larson
+

docs/pages/README_install.md

@@ -9,7 +9,7 @@ github/pypi package installation.
 
 ## Docker Container
 
-[Install Instructions](docker_cl_instructions.md)
+[Install Instructions](https://gnssrefl.readthedocs.io/en/latest/pages/docker_cl_instructions.html)
 
 ## Local Python Install 
 

docs/pages/signal_issues.md

@@ -80,21 +80,21 @@ shown on the x-axis is NOT the reflector height (RH).  it is
 be obstructed by the new Galileo codes except for E5.
 
 
-This is E5a
 
 <img src="../_static/at01_358_205.png" width="600"/>
+This is E5a (in our software it is called frequency 205; in RINEX it is called L5)
 
-This is E5. Note that instead of nice clean peaks, it is 
+<img src="../_static/at01_358_208.png" width="600"/>
+
+This is E5. In our software it is called frequency 208 and in RINEX it is called L8).
+Note that instead of nice clean peaks, it is 
 spread out. You can also see that the E5 retrievals degrades as elevation angle increases,
 which is exactly what you would expect with the multipath delay 
 increasing with elevation angle. I would just recommend only using
 this signal for RH < 5 meters. And even then, if you are tracking
 L8, you probably also have L5, L6, and L7, so there is not a ton gained 
 by also using L8.
 
-<img src="../_static/at01_358_208.png" width="600"/>
-
-
 ##  What about L1C?
 
 I would be happy to host some results from L1C - please submit a pull request 

docs/pages/understand.md

@@ -7,14 +7,13 @@
 
 *If you would like to try out reflectometry without installing the code*
 
-I recommend you use [this web app](https://gnss-reflections.org). It 
-can show you representative results with minimal constraints. It should provide 
-results in less than 10 seconds. 
+I recommend you use [the GNSS-IR web app](https://gnss-reflections.org/api). It 
+can show you representative results in less than 10 seconds. There are also some 
+[examples](https://gnss-reflections.org/overview).
 
 ## Goals
 
-The goal of the gnssrefl python repository is to help you 
-compute (and evaluate) GNSS-based
+The goal of the gnssrefl python repository is to help you compute (and evaluate) GNSS-based
 reflectometry parameters using geodetic data. This method is often
 called GNSS-IR, or GNSS Interferometric Reflectometry. There are three main modules:
 

gnssrefl/daily_avg.py

@@ -94,10 +94,14 @@ def fbias_daily_avg(station):
 
 
 
-def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,howBig,ReqTracks,azim1,azim2,test,subdir):
+def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,
+        howBig,ReqTracks,azim1,azim2,test,subdir,plot_limits):
     """
     worker code for daily_avg_cl.py
 
+    It reads in RH files created by gnssir. Applies median filter and saves average results 
+    for further analysis
+
     Parameters
     ----------
     station : str
@@ -139,6 +143,9 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
     subdir : str
         subdirectory for output files
 
+    subdir : bool
+        whether plot limits for the median filter are shown
+
     Returns
     -------
     tv : numpy array
@@ -152,7 +159,8 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
 
     """
     xdir = os.environ['REFL_CODE']
-    print('All RH retrievals will be written to: ', alldatafile)
+    print('All RH retrievals that meet your QC criteria will be written to: ' )
+    print(alldatafile)
     allrh = open(alldatafile, 'w+')
     # put in a header
     allrh.write(" {0:s}  \n".format('% year,doy, RH(m), Month, day, azimuth(deg),freq, satNu, LSP amp,pk2noise,UTC(hr)' ))
@@ -161,17 +169,22 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
     fs = 12
     NotEnough = 0
     tv = np.empty(shape=[0, 8])
+    tv_median = np.empty(shape=[0, 9])
     tvall = np.empty(shape=[0, 7])
+    #
     ngps = []; nglo = [] ; ngal = []; nbei = []
     obstimes = []; medRH = []; meanRH = [] ; alltimes = []; meanAmp = []
+    tttimes = []
     fig,ax=plt.subplots()
     year_list = np.arange(year1, year2+1, 1)
     NumFiles = 0
     s1 = time.time()
     for yr in year_list:
         direc = xdir + '/' + str(yr) + '/results/' + station + '/' + extension + '/'
+        # counter for the legends
+        nle = 0
+        # i understand why python people like this - but then the results are not sorted ...
         if os.path.isdir(direc):
-            # i understand why python people like this - but then the results are not sorted ...
             all_files = os.listdir(direc)
             all_files = np.sort(all_files)
             #print('Number of files in ', yr, len(all_files))
@@ -189,6 +202,7 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
                             a = np.loadtxt(fname,skiprows=3,comments='%')
                         nr,nc=a.shape
                         if (nr > 0):
+                            nle = nle + 1
                             # add the new azimuth constraint here ... 2022sep04
                             www = (a[:,5] > azim1 ) & (a[:,5] < azim2 )
                             a = a[www,:]
@@ -203,6 +217,7 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
 
                             medv = np.median(rh)
                             # 0 means use all frequencies.  otherwise, you can specify 
+                            # this is applying the medial filter
                             if fr == 0:
                                 cc = (rh < (medv+howBig))  & (rh > (medv-howBig))
                             else:
@@ -249,7 +264,9 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
                                     ijk = (gsat > 200) * (gsat < 300); # galileo
                                     ngal = np.append(ngal, len(gsat[ijk]))
 
+                                medv_obstime = datetime.datetime(year=yr, month=d.month, day=d.day, hour=12, minute=0, second=0)
                                 obstimes.append(datetime.datetime(year=yr, month=d.month, day=d.day, hour=12, minute=0, second=0))
+
                                 medRH.append(medv)
                                 #medRH =np.append(medRH, medv)
                                 # store the meanRH after the outliers are removed using simple median filter
@@ -264,8 +281,12 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
                                 # updated this to include mean amplitude 2021 november 8
                                 meanAmp.append(np.mean(goodAmp))
                                 newl = [yr, doy, meanRHtoday, len(rh), d.month, d.day, stdRHtoday, np.mean(goodAmp)]
+                                # a new variable is not really needed - but I did not want to oerwrite working code
+                                newl_plus_median = [yr, doy, meanRHtoday, len(rh), d.month, d.day, stdRHtoday, np.mean(goodAmp),medv]
 
                                 tv = np.append(tv, [newl],axis=0)
+                                tv_median = np.append(tv_median, [newl_plus_median],axis=0)
+
                                 k += 1
                             else:
                                 NotEnough = NotEnough + 1
@@ -276,13 +297,31 @@ def readin_plot_daily(station,extension,year1,year2,fr,alldatafile,csvformat,how
     #meanRH = np.asarray(meanRH)
     s2 = time.time()
 
+    # not sure you can sort obstimes
+    dumb_time = tv_median[:,0] + tv_median[:,1]/365.25
+    # sort it ...
+    ii = np.argsort(dumb_time)
+    tv_median = tv_median[ii,:]
+    for www in range(0,len(tv_median)):
+        filler = datetime.datetime(year=int(tv_median[www,0]), 
+                month=int(tv_median[www,4]), day=int(tv_median[www,5]), hour = 12, minute=0, second = 0)
+        tttimes.append(filler)
+
+    # plot the median 
+    ax.plot(tttimes, tv_median[:,8],'ks',markerfacecolor='white',label='median value')
+    # if requested, also show the limits for the median filter
+    if plot_limits:
+        ax.plot(tttimes, tv_median[:,8]+howBig,'-',color='gray',label='median+limit')
+        ax.plot(tttimes, tv_median[:,8]-howBig,'-',color='gray',label='median-limit')
 
-    fig.autofmt_xdate()
     plt.ylabel('meters',fontsize=fs)
     plt.title('All Reflector Heights for Station: ' + station.upper(),fontsize=fs)
     plt.xticks(fontsize=fs)
     plt.yticks(fontsize=fs)
     plt.gca().invert_yaxis()
+    # this command changes the x-axis
+    fig.autofmt_xdate()
+    plt.legend(loc="upper right")
     plt.grid()
 #   new plot
     print('A total of ', NumFiles, ' days were evaluated.')
@@ -460,8 +499,9 @@ def write_out_RH_file(obstimes,tv,outfile,csvformat):
     # 2021 november 8 added another column that has mean amplitude
     fout.write("{0:28s} \n".format( '% calculated on ' + xxx ))
     fout.write("% Daily average reflector heights (RH) calculated using the gnssrefl software \n")
+    fout.write("% Nominally you should assume this average is associated with 12:00 UTC hours \n")
     fout.write("% year doy   RH    numval month day RH-sigma RH-amp\n")
-    fout.write("% year doy   (m)                      (m)    (v/v)\n")
+    fout.write("%            (m)                       (m)    (v/v)\n")
     fout.write("% (1)  (2)   (3)    (4)    (5)  (6)   (7)     (8) \n")
     if csvformat:
         for i in np.arange(0,N,1):