updated documentation

Author: kristinemlarson

docs/pages/quick_recall.md

@@ -56,7 +56,7 @@ those discussion links will likely go away in the not very distant future.
 * [gpsweek](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.gpsweek.html)
 * [installexe](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.installexe.html)
 * [llh2xyz](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.llh2xyz.html)
-* [max_resolve_RH](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.max_resolve_cl.html)
+* [max_resolve_RH](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.max_resolve_RH_cl.html)
 * [mjd](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.mjd.html)
 * [query_unr](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.query_unr.html)
 * [quickplt](https://gnssrefl.readthedocs.io/en/latest/api/gnssrefl.quickplt.html)

docs/pages/understand.md

@@ -1,29 +1,21 @@
 # Understanding
 
-
-## Overview
-
 **gnssrefl** is an open source/python version of my GNSS interferometric reflectometry (GNSS-IR) code. 
 
-*If you would like to try out reflectometry without installing the code*
-
-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
-reflectometry parameters using geodetic data. This method is often
-called GNSS-IR, or GNSS Interferometric Reflectometry. There are three main modules:
+reflectometry parameters using standard GNSS data. This method is often
+called GNSS-IR, or GNSS Interferometric Reflectometry. There are three main sections:
 
-* [**rinex2snr**](rinex2snr.md) translates RINEX files into SNR files needed for analysis.
+* Translation: Use either [**rinex2snr**] or [**nmea2snr**] to translate native 
+GNSS formats to what gnssrefl needs. The output is called a *SNR* file.
 
-* [**quickLook**](quickLook.md) gives you a quick (visual) assessment of SNR file without dealing
+* [**quickLook**] gives you a quick (visual) assessment of a SNR file without dealing
 with the details associated with **gnssir**. It is not meant to be used for routine analysis.
 It also helps you pick an appropriate azimuth mask and quality control settings.
 
-* [**gnssir**](gnssir.md) computes reflector heights (RH) from SNR files.
+* [**gnssir**] computes reflector heights (RH) from SNR files.
 
 There are also various [utilities](quick_recall.md) you might find to be useful.
 If you are unsure about why various restrictions are being applied, it is really useful 
@@ -104,8 +96,8 @@ are not violating the Nyquist frequency (**Maximum Resolvable Reflector Height**
 
 ## Quality Control 
 
-This code uses a Lomb Scargle periodogram. This type of periodogram allows the input 
-data to be sampled at uneven periods. The primary inputs are 
+This code uses a Lomb Scargle periodogram (LSP). This type of periodogram allows the input 
+data to be sampled at uneven periods. The primary inputs are :
 
 * how precise (in reflector height units) do you want the periodogram calculated at?
 
@@ -117,11 +109,9 @@ makes no sense (i.e. so small the code takes forever to run). In this code the s
 parameter is the max reflector height (h2). The minimum reflector height is always zero, and then
 the values lower than the minimum reflector height (h1) are thrown out.
 
-
 It is easy to compute a periodogram and pick the maximum value so as to find the reflector height. It is 
 more difficult to determine whether it is one you should trust.
 
-
 * is the peak larger than a user-defined value  (amplitude of the dominant peak in your periodogram)
 
 * is the peak divided by a "noise" metric larger than a user-defined value. This noise metric is defined 
@@ -252,8 +242,8 @@ D. Reflection zones for GPS satellites at
 elevation angles of 5-15 degrees for a reflector height of 4 meters.  
 
 Again using the reflection zone web app, we can plot up the appropriate reflection zones for various options.
-Since <code>ross</code> has been around a long time, [http://gnss-reflections.org](https://gnss-reflections.org) has its coordinates in a 
-database. You can just plug in <code>ross</code> for the station name and leave 
+Since <code>ross</code> has been around a long time, [http://gnss-reflections.org](https://gnss-reflections.org) 
+has its coordinates in a database. You can just plug in <code>ross</code> for the station name and leave 
 latitude/longitude/height blank. You *do* need to plug in a 
 RH of 4 since mean sea level would not be an appropriate 
 reflector height value for this 
@@ -444,3 +434,6 @@ Note that the names of the columns (and units) are provided
 When multiplied by RHdot (meters/hour), you will get a correction in units of meters. For further
 information, see the <code>subdaily</code> code.
 
+Kristine M. Larson
+
+January 11, 2024

gnssrefl/download_noaa.py

@@ -1,10 +1,12 @@
 # -*- coding: utf-8 -*-
 import argparse
 import datetime
+import numpy as np
 import os
 import requests
 import subprocess
 import sys
+import wget
 import gnssrefl.gps as g
 
 #from gnssrefl.utils import validate_input_datatypes, str2bool
@@ -402,6 +404,78 @@ def pickup_from_noaa(station,date1,date2,datum, printmeta):
 
     return data, error
 
+def download_qld(station,year,plt):
+    """
+    Parameters
+    ----------
+    station : str
+        tide gauge station name
+    year : int
+        calendar year
+    plt : bool
+        whether you want a plot to the screen
+    """
+    # unfortunately hardwired for 2022 now
+    #https://www.qld.gov.au/environment/coasts-waterways/beach/storm/storm-sites
+    # last seven days
+    # https://www.data.qld.gov.au/dataset/coastal-data-system-near-real-time-storm-tide-data
+    if year is None:
+        print('Year is required for the Queensland option')
+        sys.exit()
+    else:
+        url = 'https://www.data.qld.gov.au/dataset/179c7cc5-26a7-4f57-9e1e-3ec2ed5dd4de/resource/cacc9e98-be38-44f7-a535-07acd22a3b91/'
+        url2 = 'download/h071004a_' + str(year) + '_' + station + '_10min.csv'
+        tmpfile = station + '_' + str(year) + '_10min.csv'
+        wget.download(url+url2, tmpfile)
+
+    xdir = os.environ['REFL_CODE']
+    outdir = xdir  + '/Files/'
+    if not os.path.exists(outdir) :
+        subprocess.call(['mkdir', outdir])
+    outfile = outdir + station + '_' + str(year) + '.txt'
+
+    fout = open(outfile, 'w+')
+    print('Queensland Tide file written to :', outfile)
+
+    obstimes=[]
+    sl = []
+    fout.write("{0:s} \n".format('%' + ' QLD Station: ' + station ))
+    fout.write("%YYYY MM DD HH MM  SS Water(m) DOY    MJD  \n")
+    fout.write("% 1    2  3  4  5   6   7       8     9\n")
+
+    if os.path.isfile(tmpfile):
+        x=np.loadtxt(tmpfile,usecols=(0,1,3),skiprows=40,dtype='str',delimiter=',')
+        nr,nc=np.shape(x)
+        for i in range(0,nr):
+            sealevel = float(x[i,2])
+            sl.append(sealevel)
+            d = int(x[i,0][0:2]) ; m = int(x[i,0][3:5])
+            y = int(x[i,0][6:10])
+            hr = int(x[i,1][1:3]); mm = int(x[i,1][4:6])
+            bigT = datetime.datetime(year=y, month=m, day=d,hour=hr,minute=mm,second=0)
+            obstimes.append(bigT)
+            modjul, fr = g.mjd(y,m,d,hr,mm,0)
+            mjd = modjul+fr
+            doy,cdoy,cyyyy,cyy = g.ymd2doy(y, m, d )
+            fout.write(" {0:4.0f} {1:2.0f} {2:2.0f} {3:2.0f} {4:2.0f} {5:2.0f} {6:7.3f} {7:3.0f} {8:15.6f} \n".format(y, m, d, hr, mm, 0, sealevel, doy, mjd))
+
+        fout.close()
+    else:
+        print('No file was downloaded')
+        sys.exit()
+
+    if plt:
+        g.quickp(station,obstimes,sl)
+
+    # clean up - remove csv file
+    subprocess.call(['rm','-f',tmpfile])
+
+# this is more direct API call for last seven days ... I think 
+#urlL = 'https://www.data.qld.gov.au/api/3/action/datastore_search?resource_id=7afe7233-fae0-4024-bc98-3a72f05675bd'
+#endL = '&q=' + station + '&limit=' + str(NV)
+
+# https://www.data.qld.gov.au/dataset/coastal-data-system-near-real-time-storm-tide-data
+
 
 if __name__ == "__main__":
     main()

gnssrefl/download_tides.py

@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 import argparse
 import datetime
+import wget
 import os
 import requests
 import subprocess
@@ -17,14 +18,15 @@
 def parse_arguments():
     parser = argparse.ArgumentParser()
     parser.add_argument("station", help="station name, e.g. 8768094", type=str)
-    parser.add_argument("network", help="tidegauge network (noaa,ioc,psmsl,wsv)", type=str)
+    parser.add_argument("network", help="tidegauge network (noaa,ioc,psmsl,wsv,qld)", type=str)
     parser.add_argument("-date1", help="start-date, 20150101", type=str, default=None)
     parser.add_argument("-date2", help="end-date, 20150110", type=str, default=None)
     parser.add_argument("-output", default=None, help="Optional output filename", type=str)
     parser.add_argument("-plt", default=None, help="quick plot to screen", type=str)
     parser.add_argument("-datum", default=None, help="datum for NOAA", type=str)
     parser.add_argument("-sensor", default=None, help="sensor type for IOC", type=str)
     parser.add_argument("-subdir", default=None, help="optional subdirectory name for output", type=str)
+    parser.add_argument("-year", default=None, help="year for archive QLD data", type=str)
     args = parser.parse_args().__dict__
 
     # convert all expected boolean inputs from strings to booleans
@@ -35,7 +37,7 @@ def parse_arguments():
     # only return a dictionary of arguments that were added from the user - all other defaults will be set in code below
     return {key: value for key, value in args.items() if value is not None}
 
-def download_tides(station: str, network : str, date1: str = None, date2: str = None, output: str = None, plt: bool = False, datum: str = 'mllw', subdir: str = None):
+def download_tides(station: str, network : str, date1: str = None, date2: str = None, output: str = None, plt: bool = False, datum: str = 'mllw', subdir: str = None, year: int=None):
     """
     Downloads tide gauge data from four different networks (see below)
 
@@ -114,9 +116,12 @@ def download_tides(station: str, network : str, date1: str = None, date2: str =
         download_wsv.download_wsv(station, plt, output)
     elif network == 'psmsl':
         download_psmsl.download_psmsl(station, output,plt )
+    elif network == 'qld':
+        download_noaa.download_qld(station,year,plt)
     else:
         print('I do not recognize your tide gauge network')
 
+
 def main():
     args = parse_arguments()
     download_tides(**args)

gnssrefl/gps.py

@@ -5303,7 +5303,8 @@ def rinex_jp(station, year, month, day):
 
 def queryUNR_modern(station):
     """
-    Queries the UNR database for station coordinates that has been stored in sql. downloads it if necessary
+    Queries the UNR database for station coordinates that has been stored in sql. downloads 
+    the sql file and stores it locally if necessary
 
     Parameters
     -----------
@@ -5374,6 +5375,8 @@ def queryUNR_modern(station):
 
     # close the database
     conn.close()
+    # some of the Nevada Reno stations have the same names as stations used in GNSS-IR.
+
     if (station == 'moss'):
         lat= -16.434464800 ;lon = 145.403622520 ; ht = 71.418
     elif (station == 'mnis'):

gnssrefl/max_resolve_RH_cl.py

@@ -46,7 +46,7 @@ def max_resolve_RH(station: str, lat: float=None, lon: float=None, el_height: fl
         typical case for most geodetic sites, 30 seconds, elevation angles 5-15 degrees
 
     max_resolve_RH sc02 -samplerate 15 -system galileo
-        receiver sampling rate of 15 seconds and galileo
+        Assume receiver sampling rate of 15 seconds and the Galileo constellation
 
     Parameters
     ----------

gnssrefl/nyquist_libs.py

@@ -108,15 +108,29 @@ def ny_plot(station,allN, info,hires_figs):
 
     """
     # wavelengths in meters, being lazy
+    l1 = 0.19029360
     l2 = 0.24421
     l5 = 0.254828048
-    l1 = 0.19029360
+    # defaults
+    plot_l2 = True
+    plot_l5 = True
+    if 'GALILEO' in info:
+        plot_l2 = False
+    if 'GLONASS' in info:
+        plot_l5 = False
+    if 'BEIDOU' in info:
+        # the beidou signal names are all messed up
+        # but that is how they are stored in the RINEX files...
+        plot_l5 = False
+
     #fig = Figure(figsize=(8,4),dpi=360)
     #fig,ax = fig.subplots()
     fig, ax = plt.subplots(figsize=(10,6))
     ax.plot(allN[:,0],allN[:,1],'bo',label='L1')
-    ax.plot(allN[:,0],l2*allN[:,1]/l1,'ro',label='L2')
-    ax.plot(allN[:,0],l5*allN[:,1]/l1,'co',label='L5')
+    if plot_l2:
+        ax.plot(allN[:,0],l2*allN[:,1]/l1,'ro',label='L2')
+    if plot_l5: 
+        ax.plot(allN[:,0],l5*allN[:,1]/l1,'co',label='L5')
     ax.set_title(station + ': Maximum Resolvable RH /' + info, fontsize=14)
     ax.grid()
     ax.set_xlabel('Azimuth (deg)')

gnssrefl/quicklib.py

@@ -9,7 +9,7 @@
 
 import gnssrefl.gps as g
 
-def trans_time(tvd, ymd, convert_mjd, ydoy ,xcol,ycol,utc_offset):
+def trans_time(tvd, ymd, ymdhm, convert_mjd, ydoy ,xcol,ycol,utc_offset):
     """
     translates time for quickplt
 
@@ -20,6 +20,8 @@ def trans_time(tvd, ymd, convert_mjd, ydoy ,xcol,ycol,utc_offset):
     ymd : bool
         first three columns are year,month,day, hour,
         minute,second
+    ymdhm : bool
+        first five columns are year,month,day, hour, minut,
     convert_mjd : bool
         convert from MJD (column 1 designation)
         time is datetime obj
@@ -54,19 +56,17 @@ def trans_time(tvd, ymd, convert_mjd, ydoy ,xcol,ycol,utc_offset):
         print('You asked to plot column', xcol+1, ' and that column does not exist in the file')
         sys.exit()
 
-    if ymd == True:
-        year = tvd[:,0]; month = tvd[:,1]; day = tvd[:,2];
-        hour = tvd[:,3] ; minute = tvd[:,4]
+    if ymd :
         for i in range(0,len(tvd)):
-            if (tvd[i, 4]) > 0:
-                y = int(year[i]); m = int(month[i]); d = int(day[i])
-                # i am sure there is a better way to do this
-                today=datetime.datetime(y,m,d)
-                doy = (today - datetime.datetime(today.year, 1, 1)).days + 1
-                h = int(hour[i])
-                mi = int(minute[i])
-                tval.append(y + (doy +  h/24 + mi/24/60)/365.25);
-                yval.append( tvd[i,ycol]/1000)
+            bigT = datetime.datetime(year=int(tvd[i,0]), month=int(tvd[i,1]), day=int(tvd[i,2]) )
+            tval.append(bigT)
+            yval.append( tvd[i,ycol])
+    elif ymdhm:
+        for i in range(0,len(tvd)):
+            bigT = datetime.datetime(year=int(tvd[i,0]), month=int(tvd[i,1]), 
+                                     day=int(tvd[i,2]), hour=int(tvd[i,3]), minute=int(tvd[i,4]), second=0)
+            tval.append(bigT)
+            yval.append( tvd[i,ycol])
     else:
         if convert_mjd:
             mm = tvd[:,xcol]
@@ -90,7 +90,6 @@ def trans_time(tvd, ymd, convert_mjd, ydoy ,xcol,ycol,utc_offset):
             tval = tvd[:,xcol] ; yval = tvd[:,ycol]
             x1 = min(tval) ; x2 = max(tval)
 
-
     return tval, yval 
 
 def set_xlimits_ydoy(xlimits):