{"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":"**Update 8/1/22:**\nNow that the competition is over, I thought I would update this notebook to share my best RTKLIB solutions for the test data set.  This includes one of two config files I used in my final submission as well as my most recent version of RTKLIB.  This will provide a score of 1.597 on the private leaderboard and 1.753 on the public leaderboard.  Note that this doesn't include any post-processing of the RTKLIB results other than replacing the data that contains hardware clock discontinuities with other publicly available solutions, so it should not be difficult to improve on these results.  Also note that the private score for this solution is lower than my final submission, since for my submission, I merged these results with a solution from a second config file.  The other solution improved scores for the training data and the public test data but unfortunately was worse for the private test data.\n\nFor those of you who just want the submission file to use as a starting point for post-processing, I have attached it as an input file.\n\nFor those of you who have already gotten an earlier version of this notebook working and just want to update it, then you will need to follow the following steps:\n\n1) Clone the latest competition-specific RTKLIB code from the GDSC_2022 branch of the demo5 repos at https://github.com/rtklibexplorer/RTKLIB/tree/gdsc_2022 and build rnx2rtkp.exe using this code. Replace the existing version of this file in your rtklib folder with the newly built version. \n\n2) Clone the latest version of android_rinex from https://github.com/rtklibexplorer/android_rinex and replace the previous versions of these files with the new files\n\n3) Download the attached input data file (\"brdm_nav_files.zip\") and unzip these into your test folders.  Change the \"nav_files\" variable in step 4 below from \"../*MN.nav to '../BRDM*MN.rnx'.  Or follow the instructions below to download these yourself.  \n\n4) Replace your previous config file for RTKLIB with the one listed below in step 2 and update the config file name in step 4\n\n5) Replace the python code listed in step 6 with the new code below.\n\n6) Run steps 4, 5, and 6 below.\n\nFor those of you running this for the first time, follow the complete set of instructions below.\n\n\n**Update 7/11/22:**\n\nWhen I originally released this notebook a couple of months ago, the resulting score was 3.135 which was good enough to tie for first place.  Since then the scores have improved significantly and this score is no longer low enough to be competitive.  I have made a couple of updates below to improve the resulting score to 2.67 which as of today is good enough to get you into 27th place out of 461 entries.\n\nThe two changes I made were to update the base station locations to account for tetonic plate movement and to replace the rides with hardware clock discontinuities on a point by point basis instead of the entire data set.  I also added a few sort() statements to correct an issue when the code was not run in Windows.\n\n--------------------------------------------------------------------------------------------------------\n\n\n\n**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.  \n\nI describe this process in more detail in [this blog post](https://rtklibexplorer.wordpress.com/2022/01/10/google-smartphone-decimeter-challenge/) in which I share my experience working with last year's GSDC data after the competition was over.  It includes a link to download the code I used to generate RTKLIB solutions for last year's data.  Submitting these solutions to Kaggle (after the competition was over) resulted in a score of 2.15 meters which put it into fifth place on the final Private Leaderboard.\n\nI would like to encourage use of RTKLIB in this year's competition and so I am sharing an updated version of the previous code to duplicate my results with this year's data and the most recent RTKLIB code. This code gives a set of solutions that score 3.135 on this year's Public Leaderboad.  The new code is only slightly modified from the previous version, so anyone interested in using RTKLIB for this challenge can start by becoming familiar with the previous code on the previous data.  I would strongly recommend starting with the older code and data because it comes as a complete package and is easier to get results with than what I am presenting here.\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 described in the above post, I have recently created an all python subset of RTKLIB for PPK solutions. This runs somewhat slower than the C code but does make an easier development platform. In this notebook I will provide code and guidelines for working with the C code version of RTKLIB.  I will also describe working with the Python version in a separate notebook after I complete this one.\n\nIn the interestes of getting this out sooner rather than later, this will not be the \"push the button and you get an answer\" kind of notebook.  It will be more like a set of handwritten notes scribbled down quickly while actually running the experiment.  In it's current form you will need to download the pieces of code to your own system, put them together and run locally since it relies on open source compiled code.  In the future I hope to make this more user friendly, but for now it is assumed that you are fairly familiar with python and can debug simple issues that may crop up when trying to follow these instructions.\n\nI ran this exercise on a Windows PC but you should be able to run it on linux as well.  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_2022\n\n    config\n    \n    data\n    \n      test\n    \n      train\n      \n    python\n   \n      android_rinex\n      \n      rtklib-py\n      \n    rtklib\n\nIt will be easier to follow these instuctions if you use the same folder structure.","metadata":{}},{"cell_type":"markdown","source":"**Step 1: Retrieve base observation and satellite navigation files**\n\n**Update 8/1/22:**\nTo duplicate my latest results, you will need to get the satellite navigation files from CDDIS.  The UNAVCO files are easier to access but seem to be incomplete in some cases.  To download the CDDIS files directly, you will need to setup an account and then setup a .netrc file in the correct location as described at https://cddis.nasa.gov/Data_and_Derived_Products/CreateNetrcFile.html and will need to comment/uncomment the appropriate lines at the beginning of the code below.\nIf you just want the navigation files necessary for this solution, I have zipped them up and included them in the input data. Unzip them directly into the test folder.  Note that you will still need to run this step to download the base observation 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 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 UNAVCO site, in part because these files include Galileo navigation data as well as GPS and GLONASS.  Some of the other sources include only GPS and GLONASS.\n\nLast year it was sufficient to use a single base station for every data set since they were all located in a small geographic area.  This 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 are multiple data sets that start in one (UTC) day and finish 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 you will either need to eliminate the multi-day sets or manually download a file with the correct starting and stopping times.  This is most easily done from the \"User Friendly CORS\" page on the NGS site.\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 and requires an executable file from RTKLIB.\n\nYou can download the RTKLIB executables for Windows at https://github.com/rtklibexplorer/RTKLIB/releases.  Put them in the GSDC_2022/rtklib folder.  Note that these are for the demo5 fork of RTKLIB which I maintain.  You can use any fork of RTKLIB for this step, but later when we are calulating the solutions, you will need to use the demo5 fork.  If you are running in Linux, you will need to build your own executables from the source code.  This is decribed at https://rtklibexplorer.wordpress.com/2020/12/18/building-rtklib-code-in-linux/\n\n","metadata":{}},{"cell_type":"code","source":"\"\"\" \nget_base_data.py - retrieve base observation and navigation data for the\n    2022 GSDC competition \n\"\"\"\n\nimport os\nfrom datetime import datetime\nimport numpy as np\nimport requests\nimport gzip\nfrom glob import glob\n\n\n# Input parameters\n\ndatadir = r'C:\\gps\\GSDC_2022\\data\\test'\nstas = ['slac', 'vdcy', 'p222']  # Bay Area, LA. backup for Bay Area\nobs_url_base = 'https://geodesy.noaa.gov/corsdata/rinex'\n# Use these parameters to load UNAVCO navigation data.  This is not as complete as the BRDM data below but \n# does not require an account to access\nnav_url_base = 'https://data.unavco.org/archive/gnss/rinex3/nav' \nnav_file_base = 'AC0300USA_R_'  # 20210060000_01D_MN.rnx.gz\n# Use these to access the nav files from CDDIS.  This requires and account and setup of a .netrc file as \n# described at https://cddis.nasa.gov/Data_and_Derived_Products/CreateNetrcFile.html.  Make sure this file \n# is in the correct location which will be system dependent\n#nav_url_base = 'https://cddis.nasa.gov/archive/gnss/data/daily' #/2021/342/21p/ \n#nav_file_base = 'BRDM00DLR_S_' # 20213420000_01D_MN.rnx.gz'' \n\n# Make sure you have downloaded this executable before running this code\ncrx2rnx_bin = r'C:\\gps\\GSDC_2022\\rtklib\\crx2rnx.exe'\n\n\nos.chdir(datadir)\n\nfor dataset in np.sort(os.listdir()):\n    if not os.path.isdir(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(os.path.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(os.path.join(dataset, fname[:-3]), \"wb\").write(obs)\n        except:\n            # try backup CORS station\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(os.path.join(dataset, fname[:-3]), \"wb\").write(obs)\n            except:\n                print('Fail obs: %s' % dataset)\n            \n        # convert crx to rnx\n        crx_files = glob(os.path.join(dataset,'*.*d'))\n        if len(crx_files) > 0:\n            os.system(crx2rnx_bin + ' ' + crx_files[0])\n    \n    # get nav data\n    if len(glob(os.path.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, fname])\n    try:\n        obs = gzip.decompress(requests.get(url).content) # get obs and decompress    \n        # write nav data\n        open(os.path.join(dataset, fname[:-3]), \"wb\").write(obs)\n    except:\n        print('Fail nav: %s' % dataset)\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 2: Download and update 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\n\nPut this in the GSDC_2022/python/android_rinex folder.  (Temporary workaround: There is currently a path issue when running in multi-processing mode.  For now, you will need to copy the files from the android_rinex/src folder into the GSDC_2022/python folder)\n\n\n\nYou will also need a configuration file for the solution. Copy the config file below to the GSDC_2022/config folder with a file name of ppk_phone_0625.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":"# ppk_phone_0625.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-posopt1       =off        # (0:off,1:on)\npos1-posopt2       =off        # (0:off,1:on)\npos1-posopt3       =off        # (0:off,1:on,2:precise)\npos1-posopt4       =off        # (0:off,1:on)\npos1-posopt5       =off        # (0:off,1:on)\npos1-posopt6       =off        # (0:off,1:on)\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     =on         # (0:off,1:on)\npos2-arfilter      =on         # (0:off,1:on)\npos2-arthres       =3\npos2-arthresmin    =3\npos2-arthresmax    =3\npos2-arthres1      =1\npos2-arthres2      =0\npos2-arthres3      =1e-09\npos2-arthres4      =1e-05\npos2-varholdamb    =0.1        # (cyc^2)\npos2-gainholdamb   =0.01\npos2-arlockcnt     =5\npos2-minfixsats    =4\npos2-minholdsats   =5\npos2-mindropsats   =10\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-nmeaintv1      =0          # (s)\nout-nmeaintv2      =0          # (s)\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-errdoppler   =1          # (Hz)\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)\nant1-pos1          =0          # (deg|m)\nant1-pos2          =0          # (deg|m)\nant1-pos3          =0          # (m|m)\nant1-anttype       =\nant1-antdele       =0          # (m)\nant1-antdeln       =0          # (m)\nant1-antdelu       =0          # (m)\nant2-postype       =posfile    # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm,6:raw)\nant2-pos1          =0          # (deg|m)\nant2-pos2          =0          # (deg|m)\nant2-pos3          =0          # (m|m)\nant2-anttype       =\nant2-antdele       =0          # (m)\nant2-antdeln       =0          # (m)\nant2-antdelu       =0          # (m)\nant2-maxaveep      =1\nant2-initrst       =on         # (0:off,1:on)\nmisc-timeinterp    =off        # (0:off,1:on)\nmisc-sbasatsel     =0          # (0:all)\nmisc-rnxopt1       =\nmisc-rnxopt2       =\nmisc-pppopt        =\nfile-satantfile    =\nfile-rcvantfile    =\nfile-staposfile    =C:\\gps\\GSDC_2022\\config\\bases.sta\nfile-geoidfile     =\nfile-ionofile      =\nfile-dcbfile       =\nfile-eopfile       =\nfile-blqfile       =\nfile-tempdir       =\nfile-geexefile     =\nfile-solstatfile   =\nfile-tracefile     =\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Step 3: Setup base station locations\n\nSince we are dealing with multiple base stations this year, we need a separate file containing the different base locations.  Create a file named bases.sta in the C:\\gps\\GSDC_2022\\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.","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: Convert raw files and run PPK solutions**\n\nAs configured in the header below, this code will convert the raw Android files to rinex and run the RTKLIB PPK solutions for the test set.  \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, open a console window, change the directory to the contents of the \"folder\" variable, then copy and paste the contents of the \"rtkcmd\" 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_2022\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\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 = r'C:\\gps\\GSDC_2022\\data\\test'\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 \nbinpath_rtklib  = r'C:\\gps\\GSDC_2022\\rtklib\\rnx2rtkp.exe'\ncfgfile_rtklib = r'C:\\gps\\GSDC_2022\\config\\ppk_phone_0625.conf'\nsoltag_rtklib = '_rtklib' # postfix for solution file names\n\n# Setup for rtklib-py - not supported in this notebook\ncfgfile = r'C:\\gps\\GSDC_2022\\config\\ppk_phone_0510.py' # cfgfile must be absolute path\nsoltag_py = '_py0510'  # postfix for solution file names\n\nPHONES = ['GooglePixel4', 'GooglePixel4XL', 'Pixel4Modded', 'GooglePixel5', 'GooglePixel6Pro', 'XiaomiMi8', 'SamsungGalaxyS20Ultra']\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\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(folder, rovfile, basefile, navfile, solfile):\n    # create command to run solution\n    rtkcmd='%s -x 0 -y 2 -k %s -o %s %s %s %s' % \\\n        (binpath_rtklib, cfgfile_rtklib, 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 == 'SamsungS20Ultra': # 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((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 PPK 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])\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 5: 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, so the RTKLIB solution points are interpolated onto the time stamps in the sample submission file that was provided by Google.  Make sure this file (sample_submission.csv) is in the data folder.\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.  The second code block below will do this.\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 = r'C:\\gps\\GSDC_2022\\data'\nhdrlen = 25    # 25 for RTKLIB, 1 for rtklib-py\n\n# Also make sure the appropriate reference file is in the datapath\n#  test: sample_submission.csv - provided in Google data\n# train: ground_truths_train.csv - created with crete_ground_truths.py\n\n############################################################\n\nGPS_TO_UTC = 315964782  # second\n\n\n# get timestamps from existing baseline file\nos.chdir(datapath)\nif DATA_SET == 'train':\n    baseline_file = 'ground_truths_train.csv'\nelse: # 'test'\n    baseline_file = 'sample_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]\n\n# open output file\nfout =open('baseline_locations_' + DATA_SET + '_' + date.today().strftime(\"%m_%d\") + '.csv','w')\nfout.write('tripId,UnixTimeMillis,LatitudeDegrees,LongitudeDegrees\\n')\n\n# get list of data sets in data path\nos.chdir(join(datapath, DATA_SET))\ntrips = np.sort(os.listdir())\n\n# loop through data set folders\nix_b = 0\nfor 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        # check for valid folder and file\n        folder = join(trip, phone)\n        if isfile(folder):\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', 'gnss_log' + SOL_TAG + '.pos')\n        if isfile(sol_path):\n            # parse solution file\n            fields = np.genfromtxt(sol_path, invalid_raise=False, skip_header=hdrlen)\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:13]\n            msecs = (1000 * (fields[:,0] * 7 * 24 * 3600 + fields[:,1])).astype(np.int64)\n            msecs += GPS_TO_UTC * 1000\n        else:\n            print('File not found: ', sol_path)\n            msecs = msecs_base.copy()\n            pos = acc = np.zeros((len(msecs), 3))\n            qs = nss = np.zeros(len(msecs))\n            \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                                     left=1000, right=1000))\n            qsi = np.interp(msecs_base[ix_b], msecs, qs)\n            nssi = np.interp(msecs_base[ix_b], msecs, nss)\n            \n            \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]))\nfout.close()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Step 6:  Filtering out RTKLIB solutions with hardware clock discontinuites**\n\nUnfortunately, there are a couple of data sets in the test data and several more in the training data sets that have corrupted carrier phase measurements.  The RTKLIB solutions are quite poor for these datasets and are typicaly worse than the included Google baseline solutions which do not use the carrier phase measurements.  \n\nThese data sets can be identified by the HardwareClockDiscontinuityCount field in the raw Androiod log files. If the final value in this field is larger than the initial value, then the carrier phase measurements will be corrupted by the clock discontinuities.  \n\nThe code block below will replace all solution points in the bad data sets with values from a file named \"best_submission.csv\"  I recommend downloading the current best scoring notebook result from the other notebooks and renaming it to create this file.\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.\n","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 = 'baseline_locations_test_07_11.csv'\nOUT_FILE = 'submit_0711.csv'\n# specify data locations\ndatapath = r'C:\\gps\\GSDC_2022\\data'\nmax_hstd = 0.5\n\nbadClockRides = [\n    '2021-08-12-US-MTV-1/GooglePixel4',\n    '2022-04-25-US-OAK-2/GooglePixel4']\n\n\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 bad clock rides with Google baseline\nfor trip_phone in np.unique(d[:,0]):\n    if trip_phone in badClockRides:\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 7: 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 close to 2.67 meters.  However, if you are using the latest pre-built RTKLIB demo5 executables, they are missing a couple of robustness improvements and so the score will be a little worse.  If you want to duplicate this score exactly, you will need to use the most recent source code available at https://github.com/rtklibexplorer/RTKLIB and compile the rnx2rktp.exe app yourself.  There are instructions for compiling the code in [Windows](https://rtklibexplorer.wordpress.com/2020/12/04/building-rtklib-code-in-windows/) or [linux](https://rtklibexplorer.wordpress.com/2020/12/18/building-rtklib-code-in-linux/) on my blog.  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.","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 5.\n\nAlso be aware that there are several rides in the training set that cross a UTC day boundary.  This code only pulls a single UTC day of navigation and base station data, so the solutions for these rides will be incomplete.  There are also a few navigation files that contain invalid entries.  The parsing code in this notebook is not robust enough to reject them and will throw an error when it encounters these.  There are a large number of data set in the training data so the simplest solution is to just delete the rides that cause errors. \n","metadata":{}},{"cell_type":"code","source":"\"\"\"\ncreate_groundtruth_csv.py - create csv file from all training set ground truth files\n\"\"\"\n\nimport os\nfrom os.path import join, isfile\n\n\ndatapath = r'C:\\gps\\GSDC_2022\\train'\nGPS_TO_UTC = 315964782  # second\n\n# open output file\nos.chdir(datapath)\nfout =open('../ground_truths_train.csv','w')\nfout.write('tripId,UnixTimeMillis,LatitudeDegrees,LongitudeDegrees\\n')\n\n# get list of data sets in data path\ndatasets = os.listdir(datapath)\n\n# loop through data set folders\nfor dataset in datasets:\n    if isfile(dataset):\n        continue\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            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\\n' % ((dataset, phone, t, d[2], d[3])))\n        \nfout.close()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If you find any errors or omission in this code, please let me know and I will update it.","metadata":{}},{"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.  More involved changes to the code can probably be done more easily in rtklib-py, the python version of RTKLIB.\n\nHere are a few suggestions to get you started.  Most of them do not require any modifications to the RTKLIB code, just the configuration file.\n\n1) **Raw measurement weighting.**  This solution uses only satellite elevation for weighting the input observations.  The code supports weighting with arbitrary combinations of elevation, signal strength, and receiver quality estimates.  Other research has shown that signal strength weighting should outperform elevation weighting for smartphone receivers.  The receiver accuracy estimates are also potentially valuable sources of information that is currently being discarded\n\n2) **Solution ensembles**:  In many cases you will find that adjusting the RTKLIB configuration parameters will improve some solutions and degrade others.  Weighted combinations of multiple solutions will likely out-perform any single solution.  I would recommend the variance based weighting used in my phone merge code example from last year (see link in the introduction).\n\n3) **Hardware clock discontinuites:** In this solution, the two datasets with clock discontinuties were simply replaced with the Google baseline solutions.  It should be possible to use RTKLIB to find an improvement over these.  (Update 7/11:  In the latest update this is done on a point by point basis rather than replacing the entire file)\n\n4) **Cycle slip detection/mitigation**:  This solution entirely ignores flags from the receiver indicating possible cycle slips.  Instead it relies on consistency checks to detect the cycle slips.  While ignoring the slips entirely works better than relying on them completely, there is valuable information being discarded here which could potentially improve the solution.\n\n5) **Tectonic plate movement:**  Base station movement due to tectonic plate movement is not fully accounted for in this solution and there may be some opportunity for improvement with more careful accounting.   (Update 7/11:  This is accounted for in the latest update in the bases.sta code above)\n\n6) **IMU measurements:**  IMU measurements could be incorporated into the RTKLIB solution.  The solutions noticeably degrade when the vehicles are stationary due to increased multipath effects so the simplest solution would be to detect this state and respond accordingly.  More sophisticated approachs could take more advantage of this information.\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\n","metadata":{}}]}