{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"**INTRODUCTION**\n\nRTKLIB is an open source software library used for (among other things) calculating GNSS solutions from raw observation data.  It was originally written by Tomoji Takasu of the Tokyo University of Marine Science and Technology, but there are now multiple forks available, including the demo5 fork which I maintain. When used to generate PPK (post-processing kinematic) solutions it has two advantages over the baseline solutions provided by Google.  First of all, it uses the carrier phase observations (ADR) as well as the pseduorange observations.  The carrier phase observations are more difficult to use but also have smaller errors than the pseudorange observations.  Secondly, the PPK solutions are differential, relative to a nearby known base location, rather than absolute like the Google solutions.  The differential solution allows us to difference raw observations between the rover and base which effectively cancels most of the satellite orbital, clock, and atmospheric errors, resulting in more accurate solutions.  Typically, PPK solutions also use integer ambiguity resolution to further increase accuracy, but in this case, I am using only the float solutions, since the quality of the smartphone observations makes ambiguity resolution extremely challenging. \n\nThis notebook is based on the version of the \"Getting Started with RTKLIB\" notebook which I shared at the end of last year's competition, but it is updated to run with this year's data.  It will generate a score of 1.803 on the public leaderboard.  I have also made some changes to make the code compatible with Linux as well as Windows.  In both cases, you will need to download this code to your computer and run it there.  I have successfully tested it on Windows 11 and WSL2 (Windows Subsystem for Linux 2). In general, the top of each file will have a set of input parameters.  Unless your folder names and paths are identical to mine, you will often need to update these before running.\n\nThe folder structure I use in this solution is:\n\nGSDC_2023\n\n    config\n    \n    data\n    \n      test\n    \n      train\n      \n    python\n   \n      android_rinex\n      \n    rtklib\n\nIt will be easier to follow these instuctions if you use the same folder structure.\n\nMy hope is to provide a platform which will allow competitors to jump right into extending the existing GNSS theory rather than having to build a solution from scratch. In addition to the C version of RTKLIB I describe in this notebook, I have also created rtklib-py, an all python subset of RTKLIB for PPK solutions. The python code runs somewhat slower than the C code but it does make for an easier development platform. This notebook will only cover working with the C code version of RTKLIB but I describe how to use the python version in the \"Getting Started with rtklib-py\" notebook available in the notebook section of last year's competition. ","metadata":{}},{"cell_type":"markdown","source":"**Step 1: Retrieve base observation and satellite navigation files**\n\nSince these are differential solutions, we will need raw observation measurements from a nearby base station for each data set.  Fortunately, these are available from the U.S.National Geodetic Survey (NGS) website.  We will also need satellite orbital data for each data set for the GPS, GLONASS, and Galileo constellations.  These are available from multiple websites.  I chose to retrieve them from the CDDIS site, in part because these files include Galileo navigation data as well as GPS and GLONASS, and because last year they appeared to be more complete than the data from other sources.\n\nTo download the CDDIS files yourself, you will need to setup a free account and then create a .netrc file in your user home directory as described at https://cddis.nasa.gov/Data_and_Derived_Products/CreateNetrcFile.html. \n\nIf you just want the base observation and navigation files necessary for the test and training data sets, I have zipped them up and included them in the input data for this notebook. Unzip them directly into the test or train folder.\n\nThis year there are data sets from the Bay Area and from the LA area, so we will need to select the appropriate base station for each data set and also use the correct location for that base station in the solution.\n\nThe code below simply retrieves the base and navigation data for the full day corresponding to the starting time of each data set.  This works fine for the test data set since all data sets start and end on the same UTC day but but in the training set there is one data set that starts in one (UTC) day and finishes in the next day.  I will be demonstrating this exercise on the test data so will not worry about this issue but if you are trying to retrieve this data for the training data I suggest just removing this one data set.\n\nThe observation files are doubly compressed.  They first need to be decompressed with gzip and then with crx2rnx.  This second step translates from compressed rinex to uncompressed rinex format.  This requires the CRX2RNX executable available for both Windows and Linux at https://terras.gsi.go.jp/ja/crx2rnx.html.  Put this file in the \"rtklib\" folder\n","metadata":{}},{"cell_type":"code","source":"\"\"\" \nget_base_data.py - retrieve base observation and navigation data for the\n    2023 GSDC competition \n\"\"\"\n\nimport os\nfrom os.path import join\nfrom datetime import datetime\nimport numpy as np\nimport requests\nimport gzip\nfrom glob import glob\nimport subprocess\n\n# Input parameters\ndatadir = '../data/test'  # relative to python script\n# List of CORS stations to use\nstas = ['slac', 'vdcy', 'p222']  # Bay Area, LA, backup for Bay Area\n\n# site to retrieve base observation data\nobs_url_base = 'https://geodesy.noaa.gov/corsdata/rinex'  \n\n# site to retrieve satellite navigation data\nnav_url_base = 'https://cddis.nasa.gov/archive/gnss/data/daily' #/2021/342/21p/ \nnav_file_base = 'BRDM00DLR_S_' # 20213420000_01D_MN.rnx.gz\n# Access to CDDIS navigation data requires registering for a free account and \n# setup of a .netrc file as described at \n# https://cddis.nasa.gov/Data_and_Derived_Products/CreateNetrcFile.html.  \n# Make sure this file  is in the users home directory \n\n# Make sure you have downloaded this executable before running this code\ncrx2rnx_bin = '../../rtklib/CRX2RNX' # relative to data directory\n\n# Loop through data sets in the data directory\nos.chdir(datadir)\nfor dataset in np.sort(os.listdir()):\n    if not os.path.isdir(join(dataset)):\n        continue\n    print(dataset)\n    ymd = dataset.split('-')\n    doy = datetime(int(ymd[0]), int(ymd[1]), int(ymd[2])).timetuple().tm_yday # get day of year\n    doy = str(doy).zfill(3)\n    \n    if len(glob(join(dataset,'*.*o'))) == 0:\n        # get obs data\n        i = 1 if '-lax-' in dataset else 0  # use different base for LA\n        fname = stas[i] + doy + '0.' + ymd[0][2:4] + 'd.gz'\n        url = '/'.join([obs_url_base, ymd[0], doy, stas[i], fname])\n        try:\n            obs = gzip.decompress(requests.get(url).content) # get obs and decompress\n            # write obs data\n            open(join(dataset, fname[:-3]), \"wb\").write(obs)\n        except:\n            # try backup CORS station\n            print('Try backup CORS:', dataset)\n            i += 2\n            fname = stas[i] + doy + '0.' + ymd[0][2:4] + 'd.gz'\n            url = '/'.join([obs_url_base, ymd[0], doy, stas[i], fname])\n            try:\n                obs = gzip.decompress(requests.get(url).content) # get obs and decompress\n                # write obs data\n                open(join(dataset, fname[:-3]), \"wb\").write(obs)\n            except:\n                print('Fail obs: %s' % dataset)\n            \n        # convert compact rinex to rinex\n        crx_files = glob(join(dataset,'*.*d'))\n        if len(crx_files) > 0:\n            subprocess.call([crx2rnx_bin, '-f', crx_files[0]])\n    \n    # get nav data\n    if len(glob(join(dataset,'*.rnx'))) > 0:\n           continue  # file already exists\n    fname = nav_file_base + ymd[0] + doy + '0000_01D_MN' + '.rnx.gz'\n    url = '/'.join([nav_url_base, ymd[0], doy, ymd[0][2:4]+'p', fname])\n    try:\n        obs = gzip.decompress(requests.get(url).content) # get obs and decompress    \n        # write nav data\n        open(join(dataset, fname[:-3]), \"wb\").write(obs)\n    except:\n        print('Fail nav: %s' % dataset)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 2: Download android_rinex library and create RTKLIB config file**\n\nYou will need the android_rinex library for converting the raw Android observation files to rinex format.  RTKLIB post processing solutions require that the input files be in the rinex format.  You will need to use my fork of this library which is available at the address shown in the code below. Make sure you have the most recent version which was updated on 8/1/22\n\nPut this in the GSDC_2023/python/android_rinex folder.  There is currently a path issue when running in multi-processing mode.  To avoid this, you will need to copy the files from the android_rinex/src folder into the GSDC_2023/python folder.\n\nYou will also need a configuration file for the RTKLIB PPK solutions. Copy the config file below to the GSDC_2023/config folder with a file name of gsdc_2023_config1.conf","metadata":{}},{"cell_type":"code","source":"\ngit clone https://github.com/rtklibexplorer/android_rinex.git\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# gsdc_2023_config1.conf - config file for RTKLIB PPK solution\n\npos1-posmode       =kinematic  # (0:single,1:dgps,2:kinematic,3:static,4:static-start,5:movingbase,6:fixed,7:ppp-kine,8:ppp-static,9:ppp-fixed)\npos1-frequency     =l1+l2+l5   # (1:l1,2:l1+l2,3:l1+l2+l5,4:l1+l2+l5+l6)\npos1-soltype       =combined-nophasereset # (0:forward,1:backward,2:combined,3:combined-nophasereset)\npos1-elmask        =5          # (deg)\npos1-snrmask_r     =on         # (0:off,1:on)\npos1-snrmask_b     =off        # (0:off,1:on)\npos1-snrmask_L1    =28,28,28,28,28,28,28,28,28\npos1-snrmask_L2    =34,34,34,34,34,34,34,34,34\npos1-snrmask_L5    =20,20,20,20,20,20,20,20,20\npos1-dynamics      =on         # (0:off,1:on)\npos1-tidecorr      =off        # (0:off,1:on,2:otl)\npos1-ionoopt       =brdc       # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc)\npos1-tropopt       =saas       # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad)\npos1-sateph        =brdc       # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom)\npos1-exclsats      =           # (prn ...)\npos1-navsys        =13         # (1:gps+2:sbas+4:glo+8:gal+16:qzs+32:bds+64:navic)\npos2-armode        =off        # (0:off,1:continuous,2:instantaneous,3:fix-and-hold)\npos2-gloarmode     =off        # (0:off,1:on,2:autocal,3:fix-and-hold)\npos2-bdsarmode     =off         # (0:off,1:on)\npos2-arelmask      =15         # (deg)\npos2-arminfix      =10\npos2-armaxiter     =1\npos2-elmaskhold    =15         # (deg)\npos2-aroutcnt      =1\npos2-maxage        =30         # (s)\npos2-syncsol       =off        # (0:off,1:on)\npos2-slipthres     =0.1        # (m)\npos2-dopthres      =10         # (m)\npos2-rejionno      =5          # (m)\npos2-rejcode       =10         # (m)\npos2-niter         =1\npos2-baselen       =0          # (m)\npos2-basesig       =0          # (m)\nout-solformat      =llh        # (0:llh,1:xyz,2:enu,3:nmea)\nout-outhead        =on         # (0:off,1:on)\nout-outopt         =on         # (0:off,1:on)\nout-outvel         =off        # (0:off,1:on)\nout-timesys        =gpst       # (0:gpst,1:utc,2:jst)\nout-timeform       =tow        # (0:tow,1:hms)\nout-timendec       =3\nout-degform        =deg        # (0:deg,1:dms)\nout-fieldsep       =\nout-outsingle      =off        # (0:off,1:on)\nout-maxsolstd      =0          # (m)\nout-height         =ellipsoidal # (0:ellipsoidal,1:geodetic)\nout-geoid          =internal   # (0:internal,1:egm96,2:egm08_2.5,3:egm08_1,4:gsi2000)\nout-solstatic      =all        # (0:all,1:single)\nout-outstat        =residual   # (0:off,1:state,2:residual)\nstats-eratio1      =200\nstats-eratio2      =300\nstats-eratio5      =25\nstats-errphase     =0.005      # (m)\nstats-errphaseel   =0          # (m)\nstats-errphasebl   =0          # (m/10km)\nstats-snrmax       =45         # (dB.Hz)\nstats-errsnr       =0.005      # (m)\nstats-errrcv       =0          # ( )\nstats-stdbias      =30         # (m)\nstats-stdiono      =0.03       # (m)\nstats-stdtrop      =0.3        # (m)\nstats-prnaccelh    =0.5        # (m/s^2)\nstats-prnaccelv    =0.1        # (m/s^2)\nstats-prnbias      =0.001      # (m)\nstats-prniono      =0.01       # (m)\nstats-prntrop      =0.001      # (m)\nstats-prnpos       =0          # (m)\nstats-clkstab      =0          # (s/s)\nant1-postype       =llh        # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw)\nant2-postype       =posfile    # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw)\nant2-maxaveep      =1\nant2-initrst       =on         # (0:off,1:on)\nmisc-timeinterp    =off        # (0:off,1:on)\nfile-satantfile    =\nfile-rcvantfile    =\nfile-staposfile    =../../../../config/bases.sta","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 3: Setup base station locations**\n\nSince we are dealing with multiple base stations, we need a separate file containing the different base locations.  Create a file named bases.sta in the C:\\gps\\GSDC_2023\\config folder and copy the lines below into this file.  RTKLIB will use the first four characters of the base station file to select the correct location from this list.  Note that if you don't use the exact same file name and folder name as I used, you will need to modify the \"file-staposfile\" parameter in the config file above.\n\nThe precise base station locations are continuously changing by small amounts due to tectonic plate movement but I have ignored the relative movement between data sets and the locations below were calculated for roughly the middle of 2022 using the base velocities specified in the coordinate files available for each base station.","metadata":{}},{"cell_type":"code","source":"%  LATITUDE(DEG) LONGITUDE(DEG)    HEIGHT(M)   NAME\n37.41651904  -122.20426828  63.778  SLAC\n34.17856659  -118.22000501  318.230  VDCY\n37.53924080  -122.08326860  53.605  P222","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 4: Download RTKLIB code**\n\nYou can download the RTKLIB executables and source code for the version of RTKLIB that I have optimized for smartphone solutions at https://github.com/rtklibexplorer/RTKLIB/releases/tag/gsdc_2022_v1.0. \n\nYou can download the Windows executables directly from here and put them in the GSDC_2023/rtklib folder.  If you are running in Linux, you will need to build your own executables from the source code. This is decribed in a blog post at https://rtklibexplorer.wordpress.com/2020/12/18/building-rtklib-code-in-linux/.\n\nNote that if you would like to build the Windows executables yourself, the Windows instructions describe using the Embarcadero compiler which is required for the GUI apps.  If you are just compiling the rnx2rtkp app, you can compile it with the VisualStudio compiler using the project file in the \\app\\consapp\\rnx2rtkp\\msc folder.","metadata":{}},{"cell_type":"markdown","source":"**Step 5: Convert the raw observation files and run PPK solutions**\n\nAs configured in the header below, this code will convert the raw Android files to RINEX format and run the RTKLIB PPK solutions for the test set.  Note that these will be float solutions, we are not attempting to resolve the integer ambiguities, since the quality of the smartphone observations is very low.\n\nIn the main code at the bottom of the file, the execution can be set up as either sequential or multiprocessing by commenting or uncommenting the appropriate lines, both for the file conversion, and for running the solutions.  It is easier to debug when run sequentially but is much slower.  I recommend running each step sequentially until you are convinced it's working, then switch it to multiprocessing.  \n\nThe solution files will all be tagged with the \"soltag_rtklib\" parameter defined in the header so you can use this to keep separate the results from multiple runs.  They will be in the \"supplemental\" folders inside each phone folder.\n\nNote that the parameters in the header are configured to overwrite all rinex files and solution files.  If you want to just rerun a subset of the data, set one or both of these parameters to False and then the code will only run when the output file is missing.\n\nThis code is setup to run either the C version of RTKLIB or the python version or both.  In this notebook, I am only addressing the C version, please see my other notebook if you would like to run the python code.\n\nDebugging hint:  If the code runs without error but does not produce any solution files then the error is very likely occurring during the call to the rtklib executable since any errors that occur in that code are not fed back to the python code.  The easiest way to debug this is to place a breakpoint in the \"run_rtklib\" function while in sequential execution mode, open a console window, change the directory to the contents of the \"folder\" variable, then copy and paste the contents of the \"rtkcmd_debug\" variable into the console window and run it.  Most likely you will find that one of the input files is missing or in the wrong location.","metadata":{}},{"cell_type":"code","source":"\"\"\"\nrun_ppk_multi.py - convert raw android files to rinex and run PPK solutions for GDSC_2023\ndata set with RTKLIB and/or rtklib-py.   \n\"\"\"\n\nimport sys\nimport numpy as np\nif 'rtklib-py/src' not in sys.path:\n    sys.path.append('rtklib-py/src')\nif 'android_rinex/src' not in sys.path:\n    sys.path.append('android_rinex/src')\n\nimport os, shutil\nfrom os.path import join, isdir, isfile, abspath\nfrom glob import glob\nfrom multiprocessing import Pool\nimport subprocess\nimport gnsslogger_to_rnx as rnx\nfrom time import time\n\n# set run parameters\nmaxepoch = None # max number of epochs, used for debug, None = no limit\n\n# Set solution choices\nENABLE_PY = False        # Use RTKLIB-PY to generate solutions \nENABLE_RTKLIB = True     # Use RTKLIB to generate solutions\nOVERWRITE_RINEX = True  # overwrite existing rinex filex\nOVERWRITE_SOL = True    # overwrite existing solution files\n\n# specify location of input folder and files\ndatadir = '../data/test'   # relative to python script\n# base and nav file locations are relative to obs files\nbasefiles = '../*0.2*o' # rinex2, use this for rtklib only\n#basefiles = '../base.obs' # rinex3, use this for python only\nnavfiles = '../BRDM*MN.rnx' # navigation files with wild cards\n\n# Setup for RTKLIB,  paths relative to python script\nbinpath_rtklib  = '../rtklib/rnx2rtkp'\ncfgfile_rtklib = '../config/gsdc_2023_config1.conf'\nsoltag_rtklib = '_rtklib' # postfix for solution file names\n\n# Setup for rtklib-py - not supported in this notebook\n#cfgfile = '../config/ppk_phone_0510.py'\nsoltag_py = '_py0510'  # postfix for python solution file names\n\n# convert relative paths to absolute paths\ndatadir = abspath(datadir)\nbinpath_rtklib = abspath(binpath_rtklib)\ncfgfile_rtklib = abspath(cfgfile_rtklib)\n\n# Select phones to process\n# all phones\n# PHONES = ['pixel4', 'pixel4xl', 'pixel5', 'pixel5a', 'pixel6pro', 'pixel7pro',\n#           'mi8', 'xiaomimi8',\n#           'sm-g988b', 'sm-g955f', 'sm-s908b', 'sm-a226b', 'sm-a600t',\n#           'sm-a505g', 'sm-a325f', 'sm-a217m', 'sm-a205u', 'sm-a505u', \n#           'samsungs22ultra', 'samsunga325g', 'samsunga32', 'samsung21ultra']\n# phones in test set\nPHONES = ['pixel4', 'pixel4xl', 'pixel5', 'pixel6pro', 'pixel7pro',\n          'mi8', 'xiaomimi8',\n          'sm-g988b', 'sm-s908b', 'sm-a325f', 'sm-a505u', 'sm-a205u',\n          'samsunga325g', 'samsunga32']\n\n# These are only for rtklib-py, see the bases.sta file described above for RTKLIB base locations\nBASE_POS = {'slac' : [-2703116.3527, -4291766.8501, 3854248.1361],  # WGS84 XYZ coordinates\n            'vdcy' : [-2497836.8748, -4654543.0665, 3563029.0635],\n            'p222' : [-2689640.5799, -4290437.1653, 3865051.0923]}\n\n# input structure for rinex conversion\nclass args:\n    def __init__(self):\n        # Input parameters for conversion to rinex\n        self.slip_mask = 0 # overwritten below\n        self.fix_bias = True\n        self.timeadj = 1e-7\n        self.pseudorange_bias = 0\n        self.filter_mode = 'sync'\n        # Optional hader values for rinex files\n        self.marker_name = ''\n        self.observer = ''\n        self.agency = ''\n        self.receiver_number = ''\n        self.receiver_type = ''\n        self.receiver_version = ''\n        self.antenna_number = ''\n        self.antenna_type = ''\n\n# Copy and read config file\nif ENABLE_PY:\n    shutil.copyfile(cfgfile, '__ppk_config.py')\n    import __ppk_config as cfg\n    import rinex as rn\n    import rtkcmn as gn\n    from rtkpos import rtkinit\n    from postpos import procpos, savesol\n\n# function to convert single rinex file\ndef convert_rnx(folder, rawFile, rovFile, slipMask):\n    os.chdir(folder)\n    argsIn = args()\n    argsIn.input_log = rawFile\n    argsIn.output = os.path.basename(rovFile)\n    argsIn.slip_mask = slipMask\n    rnx.convert2rnx(argsIn)\n\n# function to run single RTKLIB-Py solution\ndef run_ppk(folder, rovfile, basefile, navfile, solfile):\n    # init solution\n    os.chdir(folder)\n    gn.tracelevel(0)\n    nav = rtkinit(cfg)\n    nav.maxepoch = maxepoch\n    print(folder)\n\n    # load rover obs\n    rov = rn.rnx_decode(cfg)\n    print('    Reading rover obs...')\n    if nav.filtertype == 'backward':\n        maxobs = None   # load all obs for backwards\n    else:\n        maxobs = maxepoch\n    rov.decode_obsfile(nav, rovfile, maxobs)\n\n    # load base obs and location\n    base = rn.rnx_decode(cfg)\n    print('   Reading base obs...')\n    base.decode_obsfile(nav, basefile, None)\n    \n    # determine base location from original base obs file name\n    if len(BASE_POS) > 1:\n        baseName = glob('../*.2*o')[0][-12:-8]\n        nav.rb[0:3]  = BASE_POS[baseName]\n    elif nav.rb[0] == 0:\n        nav.rb = base.pos # from obs file\n        \n    # load nav data from rover obs\n    print('   Reading nav data...')\n    rov.decode_nav(navfile, nav)\n\n    # calculate solution\n    print('    Calculating solution...')\n    sol = procpos(nav, rov, base)\n\n    # save solution to file\n    savesol(sol, solfile)\n    return rovfile\n\n# function to run single RTKLIB solution\ndef run_rtklib(binpath_rtklib, cfgfile_rtklib, folder, rovfile, basefile, \n               navfile, solfile):\n    # create command to run solution\n    rtkcmd = ['%s' % binpath_rtklib, '-x', '0', '-y', '2', '-k', cfgfile_rtklib,\n              '-o', solfile, rovfile, basefile, navfile]\n    \n    # use this command line for debug from console, run from path in folder variable\n    rtkcmd_debug = '%s -x 0 -y 2 -k %s -o %s %s %s %s' % (binpath_rtklib, cfgfile_rtklib,\n              solfile, rovfile, basefile, navfile)\n    \n    # run command\n    os.chdir(folder)\n    subprocess.run(rtkcmd, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)   \n\n####### Start of main code ##########################\n\ndef main():\n\n    # get list of data sets in data path\n    datasets = np.sort(os.listdir(datadir))\n\n    # loop through data set folders\n    rinexIn = []\n    ppkIn = []\n    rtklibIn = []\n    for dataset in datasets:\n        for phone in PHONES:\n            # skip if no folder for this phone\n            folder = join(datadir, dataset, phone)\n            if not isdir(folder):  \n                continue\n            os.chdir(folder)\n            rawFile = join('supplemental', 'gnss_log.txt')\n            rovFile = join('supplemental', 'gnss_log.obs')\n\n            rinex = False\n            # check if need rinex conversion\n            if OVERWRITE_RINEX or not isfile(rovFile):\n                # generate list of input parameters for each rinex conversion\n                if phone[:7] == 'samsung': # Use cycle slip flags for Samsung phones\n                    slipMask = 0 # 1 to unmask recevier cycle slips\n                else:\n                    slipMask = 0 \n                rinexIn.append((folder, rawFile, rovFile, slipMask))\n                print(rawFile, '->', rovFile) \n                rinex = True\n            \n            # check if need to create PPK solution\n            try:\n                baseFile = glob(basefiles)[0]\n                navFile = glob(navfiles)[0]\n                solFile = rovFile[:-4] + soltag_py + '.pos'\n                solFile_rtklib = rovFile[:-4] + soltag_rtklib + '.pos'\n            except:\n                print(folder,'  Error: Missing file')\n                continue\n            if ENABLE_PY and (OVERWRITE_SOL == True or len(glob(solFile)) == 0 \n                              or rinex == True):\n                # generate list of input/output files for each python ppk solution\n                print('PY: ', join(dataset, phone))\n                ppkIn.append((folder, rovFile, baseFile, navFile, solFile))\n            if ENABLE_RTKLIB and (OVERWRITE_SOL == True or \n                        len(glob(solFile_rtklib)) == 0 or rinex == True):\n                # generate list of input/output files for each rtklib ppk solution\n                print('RTKLIB: ', join(dataset, phone))\n                rtklibIn.append((binpath_rtklib, cfgfile_rtklib,\n                        folder, rovFile, baseFile, navFile, solFile_rtklib))\n\n    if len(rinexIn) > 0:\n        print('\\nConvert rinex files...')\n        # generate rinx obs files in parallel, does not give error messages\n        #with Pool() as pool: # defaults to using cpu_count for number of procceses\n        #    res = pool.starmap(convert_rnx, rinexIn)\n        # run sequentially, use for debug\n        for input in rinexIn:\n            convert_rnx(input[0],input[1],input[2],input[3])\n\n    if ENABLE_PY and len(ppkIn) > 0:\n        print('Calculate rtklib-py solutions...')\n        # run PPK solutions in parallel, does not give error messages\n        #with Pool() as pool: # defaults to using cpu_count for number of procceses\n        #    res = pool.starmap(run_ppk, ppkIn)\n        # run sequentially, use for debug\n        for input in ppkIn:\n            run_ppk(input[0],input[1],input[2],input[3],input[4])\n\n    if ENABLE_RTKLIB and len(rtklibIn) > 0:\n        print('Calculate RTKLIB solutions...')\n        # run PPK solutions in parallel, does not give error messages\n        #with Pool() as pool: # defaults to using cpu_count for number of procceses\n        #    res = pool.starmap(run_rtklib, rtklibIn)\n        # run sequentially, use for debug\n        for input in rtklibIn:\n            run_rtklib(input[0],input[1],input[2],input[3],input[4],input[5],input[6])\n\nif __name__ == '__main__':\n    t0 = time()\n    main()\n    print('Runtime=%.1f' % (time() - t0))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 6: Combine RTKLIB solutions into a single .csv file**\n\nThe code below will read in all the individual RTKLIB solution files and create a single .csv file in the correct format for submitting to Kaggle.  The time stamps in the RTKLIB solutions may not exactly match the time stamps in the original raw data and may be missing some data points or solution files, so the RTKLIB solution points are interpolated onto the time stamps in a sample submission file. If the solution file is missing, the data from the sample submission file is used instead.  I recommend downloading the current best scoring notebook result from the other notebooks and renaming it to create this file. Make sure this file is named \"best_submission.csv\" and is in the data folder.  The file I used is in the input data and is the result from the baseline notebook from Chirag Chauhan.\n\nThis will work for the test data, but for the training data you will need to generate a reference file for the correct timestamps from the ground truth data.  This is described in the training data section below.\n\nOnly solutions with the same tag will be included so make sure you use the same tag (SOL_TAG) here as you did when creating the solutions in the previous step.\n\nThe output file name will include the test set and date and will be in the datapath folder.","metadata":{}},{"cell_type":"code","source":"\"\"\" create_baseline_csv_from_pos.py -  Create csv file PPK solution files using timestamps in reference file\n\"\"\"\n\nimport os\nfrom os.path import join, isfile\nimport numpy as np\nfrom datetime import date\n\n########### Input parameters ###############################\n\nDATA_SET = 'test'\nSOL_TAG = '_rtklib'\ndatapath = '../data' # relative to python script\nrovfile = 'gnss_log'\nhdrlen = 25    # 25 for RTKLIB, 1 for RTKLIB-py\n\noutThresh = 100   # max horizontal accuracy estimate\n# phones in test set\nPHONES = ['pixel4', 'pixel4xl', 'pixel5', 'pixel6pro', 'pixel7pro',\n          'mi8', 'xiaomimi8',\n          'sm-g988b', 'sm-s908b', 'sm-a325f', 'sm-a505u', 'sm-a205u',\n          'samsunga325g', 'samsunga32']\n# PHONES = []  # use all phones\n\n# Also make sure the appropriate reference file is in the datapath\n#  test: best_submission.csv - best available sample submission\n# train: ground_truths_train.csv - created with create_ground_truths.py\n\n############################################################\n\nGPS_TO_UTC = 315964782  # second\n\ndef create_csv(datapath, DATA_SET, SOL_TAG):\n    # get timestamps from existing baseline file\n    datapath = os.path.abspath(datapath)\n    os.chdir(datapath)\n    if DATA_SET[:5] == 'train':\n        baseline_file = 'ground_truths_' + DATA_SET + '.csv'\n    else: # 'test'\n        baseline_file = 'best_submission.csv'\n    # read data from baseline file\n    base_txt = np.genfromtxt(baseline_file, delimiter=',',invalid_raise=False, \n                             skip_header=1, dtype=str)\n    msecs_base = base_txt[:,1].astype(np.int64)\n    phones_base = base_txt[:,0]\n    pos_base = base_txt[:,2:4].astype(float) # baseline positions\n    \n    # open output file\n    fout =open('locations_' + DATA_SET + '_' + date.today().strftime(\"%m_%d\") + '.csv','w')\n    fout.write('tripId,UnixTimeMillis,LatitudeDegrees,LongitudeDegrees\\n')\n    \n    # get list of data sets in data path\n    os.chdir(join(datapath, DATA_SET))\n    trips = np.sort(os.listdir())\n    \n    # loop through data set folders\n    ix_b, npts = [], 0\n    for trip in trips:\n        if isfile(trip):\n            continue\n        phones = os.listdir(trip)\n        # loop through phone folders\n        for phone in phones:\n            if isinstance(phone, bytearray):\n                phone = phone.decode('utf-8')\n            # check for valid folder and file\n            folder = join(trip, phone)\n            if isfile(folder):\n                continue\n            if PHONES != [] and phone not in PHONES:\n                continue\n            trip_phone = trip + '/' + phone\n            #print(trip_phone)\n    \n            ix_b = np.where(phones_base == trip_phone)[0]\n            sol_path = join(folder, 'supplemental', rovfile + SOL_TAG + '.pos')\n            fields = []\n            if isfile(sol_path):\n                # parse solution file\n                fields = np.genfromtxt(sol_path, invalid_raise=False, skip_header=hdrlen)\n            if len(fields) > 1:\n                if int(fields[0,1]) > int(fields[-1,1]): # invert if backwards solution\n                    fields = fields[::-1]\n                pos = fields[:,2:5]\n                qs = fields[:,5].astype(int)\n                nss = fields[:,6].astype(int)\n                acc = fields[:,7:10]\n                msecs = (1000 * (fields[:,0] * 7 * 24 * 3600 + fields[:,1])).astype(np.int64)\n                msecs += GPS_TO_UTC * 1000\n            # if no data, use baseline data\n            if not isfile(sol_path) or len(fields) == 0:\n                print('Warning: data substitution: ', sol_path)\n                msecs = msecs_base[ix_b].copy()\n                pos = acc = np.zeros((len(msecs), 3))\n                pos[:,:2] = pos_base[ix_b].copy()\n                qs = nss = np.zeros(len(msecs))\n           \n            # interpolate to baseline timestamps to fill in missing samples\n            llhs = []; stds = []\n            for j in range(6):\n                if j < 3:\n                    llhs.append(np.interp(msecs_base[ix_b], msecs, pos[:,j]))\n                    stds.append(np.interp(msecs_base[ix_b], msecs, acc[:,j]))\n            qsi = np.interp(msecs_base[ix_b], msecs, qs)\n            nssi = np.interp(msecs_base[ix_b], msecs, nss)\n    \n            # write results to combined file\n            for i in range(len(ix_b)):\n                fout.write('%s,%d,%.12f,%.12f,%.2f,%.0f,%.0f,%.3f,%.3f,%.3f\\n' % \n                        (trip_phone, msecs_base[ix_b[i]], llhs[0][i], llhs[1][i],\n                         llhs[2][i], qsi[i], nssi[i], stds[0][i], stds[1][i], \n                         stds[2][i]))\n                try:\n                    npts += len(fields)\n                except:\n                    pass\n    \n    fout.close()\n    return npts\n\nif __name__ == '__main__':\n    create_csv(datapath, DATA_SET, SOL_TAG)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 7:  Creating final submission and filtering out problematic RTKLIB solution points**\n\nUnfortunately, there are a couple of data sets in the test data and several more in the training data sets for which the RTKLIB solutions are quite poor. These need further investigation but for now, we will replace the solution points with high estimated errors with values from the \"best_submission.csv\" file mentioned in the previous section.  \n\nNote that the previous step will generate a file with the current date in the name so you will need to rename the input file below.","metadata":{}},{"cell_type":"code","source":"\"\"\"\ncreate_submission.py - convert baseline file into submission file\n\"\"\"\n\nimport numpy as np\nimport os\n\nLOCATIONS_FILE = 'locations_test_09_20.csv'\nOUT_FILE = 'submit_0920.csv'\n# specify data locations\ndatapath = '../data'  # relative to python script\nmax_hstd = 0.5\n\nlowQualityRides = [\n    '2022-06-28-20-56-us-ca-sjc-r/samsunga32',\n    '2022-10-06-20-46-us-ca-sjc-r/sm-a205u']\n\ndatapath = os.path.abspath(datapath)\nos.chdir(datapath)\n\n# load baseline data \nbaseline_file = 'best_submission.csv'\nbase_txt = np.genfromtxt(baseline_file, delimiter=',',invalid_raise=False, \n                         skip_header=1, dtype=str)\nmsecs_base = base_txt[:,1].astype(np.int64)\nphones_base = base_txt[:,0]\npos_base = base_txt[:,2:4].astype(float)\n\n# load test data\nd = np.genfromtxt(LOCATIONS_FILE, delimiter=',',invalid_raise=False, skip_header=1, dtype=str)\nstds = d[:,7:10].astype(float)\nhstds = np.sqrt(stds[:,0]**2 + stds[:,1]**2)\n\n        \n# merge low quality rides with Google baseline\nfor trip_phone in np.unique(d[:,0]):\n    if trip_phone in lowQualityRides:\n        ixt = np.where(d[:,0] == trip_phone)[0]\n        #ix = ixt[np.where(d[ixt,5] != '2')[0]]\n        ix = ixt[np.where(hstds[ixt] >= max_hstd)[0]]\n        d[ix,2:4] = pos_base[ix,0:2]\n\n# save results to file\nfout =open( OUT_FILE,'w')\nfout.write('tripId,UnixTimeMillis,LatitudeDegrees,LongitudeDegrees\\n')\nfor i in range(len(d)):\n    # write results to combined file\n    fout.write('%s, %s, %3.12f, %3.12f\\n' % (d[i,0], d[i,1], float(d[i,2]), float(d[i,3])))\nfout.close()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 8: Submit CSV file to Kaggle**\n\nYou can now submit the csv file created in the previous step to Kaggle.  This should give you a score of 1.803 meters on the public leaderboard.","metadata":{}},{"cell_type":"markdown","source":"**To run this code on the training data**\n\nFor the most part, you can use the same code to generate solutions for the training data simply by changing all references above from the test folder to the train folder.  However you will need to run the code below first to create a reference file for combining the solution files in step 6.","metadata":{}},{"cell_type":"code","source":"\"\"\"\ncreate_groundtruth_csv.py - create csv file from all ground truth files\n\"\"\"\n\nimport os\nfrom os.path import join, isfile\n\n\ndatapath = '../data/train' # relative to python script\n# Select phones to process\n# all phones\n# PHONES = ['pixel4', 'pixel4xl', 'pixel5', 'pixel5a', 'pixel6pro', 'pixel7pro',\n#           'mi8', 'xiaomimi8',\n#           'sm-g988b', 'sm-g955f', 'sm-s908b', 'sm-a226b', 'sm-a600t',\n#           'sm-a505g', 'sm-a325f', 'sm-a217m', 'sm-a205u', 'sm-a505u', \n#           'samsungs22ultra', 'samsunga325g', 'samsunga32', 'samsung21ultra']\n\n# just phones in test set\nPHONES = ['pixel4', 'pixel4xl', 'pixel5', 'pixel6pro', 'pixel7pro',\n          'mi8', 'xiaomimi8',\n          'sm-g988b', 'sm-s908b', 'sm-a325f', 'sm-a505u', \n          'samsunga325g', 'samsunga32']\nGPS_TO_UTC = 315964782  # second\n\n# open output file\ndatapath = os.path.abspath(datapath)\nos.chdir(datapath)\nfout =open('../ground_truths_train.csv','w')\nfout.write('tripId,UnixTimeMillis,LatitudeDegrees,LongitudeDegrees, Height, Heading\\n')\n\n# get list of data sets in data path\ndatasets = sorted(os.listdir(datapath))\n\n# loop through data set folders\nfor dataset in datasets:\n    if isfile(dataset):\n        continue\n    try:\n        phones = sorted(PHONES)\n    except:\n        phones = os.listdir(join(datapath,dataset))\n    for phone in phones:\n        folder = join(datapath, dataset, phone)\n        if isfile(folder):\n            continue\n        \n        csv_file = join(folder, 'ground_truth.csv')\n        if not isfile(csv_file):\n            continue\n\n        # parse ground truth file\n        with open(csv_file) as f:\n            lines = f.readlines()[1:]\n        flag = 0\n        for line in lines:\n            if len(line) <= 1:\n                continue\n            d = line.split(',')\n            t = float(d[8]) # get time stamp\n            if flag == 0:            \n                print('%20s,%16s' % (dataset, phone))\n                flag = 1\n            # write results to combined file\n            fout.write('%s/%s,%.0f,%s,%s,%s, %s\\n' % ((dataset, phone, t, d[2],\n                                                      d[3], d[4][:7], d[7][:5])))\n        \nfout.close()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Final thoughts**\n\nThe intent of this notebook is not to provide a fully optmized solution, but only to get you started with RTKLIB and demonstrate some of its capability.  \n\nFollowing these instructions will provide an improved baseline solution file which can be post-processed with filtering, map-matching, etc to give you a good jump on the competition. This alone, however, will probably not be enough to win the competion.  To do that I believe that you will need to improve the RKTLIB solution itself.  Some of this can be done by modifying the configuration file.  More dramatic changes will require modifying the code.  More information on the configuration file and the code algorithms are available in the [demo5 RTKLIB Users Manual](https://rtkexplorer.com/pdfs/manual_demo5.pdf), particularly section 3.5 and Appendix F for information on configuration, and Appendix E for information on the core algorithms.    \n\nThere are instructions for compiling the code in Windows or Linux on my blog (https://rtklibexplorer.wordpress.com/). Note that the Windows instructions use the Embarcadero compiler which is required for the GUI apps but if you are just compiling the rnx2rtkp app, you can compile it with the VisualStudio compiler using the project file in the \\app\\consapp\\rnx2rtkp\\msc folder.  More involved changes to the code may be done more easily in rtklib-py, the python version of RTKLIB.\n\nI'm happy to answer any questions regarding RTKLIB.  I just ask that, to follow the rules of the competition, you ask your questions in the discussion group here so that the answers are available to all of the competitors.  \n\nIf you find any errors or omission in this code, please let me know and I will update it.\n\nMore details of the optimizations I have made to RTKLIB for smart phone observations are described in these links:\n\n* [RTKLIBexplorer blog post: Google Smartphone Decimeter Challenge](http://https://rtklibexplorer.wordpress.com/2022/01/10/google-smartphone-decimeter-challenge/)\n\n* [Optimizing the Use of RTKLIB for Smartphone-Based GNSS Measurements](http://https://www.mdpi.com/1424-8220/22/10/3825)\n\n* [3rd Place Winner: 2022 Smartphone Decimeter Challenge: An RTKLIB Open-Source Based Solution](http://https://www.ion.org/publications/abstract.cfm?articleID=18376)","metadata":{}}]}