Python, GDAL OGR and RTKLIB

[tracy.love]

Python, GDAL OGR and RTKLIB

hey all, i have written some python code to manage processing shapefiles from the reach receivers and i want to share it. Hopefully if someone else needs something like this maybe my code will help you out. This is Python 3 written on windows 10 and calls GDAL OGR libraries.

WARNING this code is designed to work in my semi automated work arrangement so its not going to work for you out of the box. You will need to rework it to fit your situation. furthermore just because this code works for me doesn’t mean it will work for you. for example anything having to do with date, times and latitudes and longitudes may only be set to work for where i am at (California, USA).

first here is my background local library for functions i use across everything i write. it includes code for leveraging my OSGEO4w install.

http://public.region3dfg.org/apps/python_geo_processing/r3_gis.py

you save it to :

C:<your python install>\Lib\site-packages

then i have a simple script that takes things like the ‘projected’ field and parses it out into separate fields. there’s a bunch in there you may not need.

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_ReachDataProcess/R3GIS_emlid-reach-data-process_1.py

then i wrote a processing front end for RTKLIB PPK. the underlying idea here is that i was never able to get everything i wanted out of one configuration set for rtkpost. so i decided i would rather run configuration files in a serial fashion where configurations run from most conservative to least conservative. the code creates a shapefile from the first PPK process (.pos) and then with subsequent config file runs, if there are fixes they overwrite any remaining float values. then if you have a shapefile of points from your rover the resulting kinematic shapefile will be used to improve the points in your rover file records

here are my config files:

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/rtklib_reach-rs_GNSS_1hz/

numbered 1_.conf where the number represents what i believe to be the most conservative settings and 6_.conf is the least conservative settings.

here is the code

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/R3GIS_rtklib.py

and here is a test set you can use

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/R3GIS_RTKlib.zip

in an effort to develop better techniques to improve productivity, the team i work with has laid out a testing and training site. the site has a number of OPUS based monuments and using a total station we created a monument in a heavy canopy area. The canopy cover is heavy enough that i have been unable to get a real time fix even with the base station within 50 meters. however, i have been able to get PPK float values within 20cm and all of my PPK results that were fixes, where i have a matching RTK value, also show significant improvement (compare to the OPUS control monuments). The main advantage i am looking for is to not have to fiddle at all with the RTKLIB GUI after each project, (though i still do the UBX conversion there).

so there it is, maybe someone will find it useful, if you write Python it should be fairly straight forward reading to see what it does and how.

and i would like to thank Tim Everett for his work at https://rtklibexplorer.wordpress.com

[tracy.love]

well i have done allot of cleaning of the code and added in weighted point averaging (numpy) and calculating the distance each point is adjusted (pyproj).

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_ReachDataProcess/R3GIS_emlid-reach-data-process_1.py

and

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/R3GIS_rtklib.py

[polina.buriak]

Hi @tracy.love,

Thanks for sharing your work and updates on it with us!

It’d be great to know what results our users would be able to achieve with it

May I ask you what you find uncomfortable about RTKLib GUI?

[tracy.love]

Hey Polina

I am really very new to using rtklib so take what i have to say keeping that in mind. RTKLIB is a great program, i love that its portable and that it uses configuration files. but if your using the reach units allot then theres no reason why not to do some automation to make things smoother. with the GUI it is a pain to perform multiple runs on the same file and integrate them together manually. furthermore, i wanted to be able easily PPK improve my rover survey file without much hassel. now if i loose my RTK fix (which is more often than i like) i just sit on my next point for a while till i figure i have enough data to successfully PPK my point. at least for my workflow, the code i wrote really makes it painless.

[RTK_Hunter]

Nice compilation right here. I will give it a spin, thanks for sharing!

[tracy.love]

sure have at it, i just uploaded a new version, chasing down a few more bugs, last set of bugs were due to not testing it on 5hz data. biggest issue left is that it is slow. like run it over night slow on a 30mb raw file. bottle neck is where it goes to integrate each new set of fixes into the composite layer. so the first layer is the product of
http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/rtklib_reach-rs_GNSS_5hz/1_combined-continus_pos2-arthres1_004.conf

the next set of fixes (points not fixed by the conf. above) are then created from
http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/rtklib_reach-rs_GNSS_5hz/2_combined-continus_pos2-arthres1_01.conf

they then have to be integrated into the composite layer

i believe the bottle neck is with
PPKtempFe = outLayer.GetNextFeature()
while PPKtempFe:
daid = PPKtempFe.GetField(“gnss_id”)
combinedPPKlayer.SetAttributeFilter(‘’)
combinedPPKlayer.SetAttributeFilter(“gnss_id = '”+daid+“'”)

which is a bit odd since with each RTKLIB conf file i run i get progressively diminishing returns. and the bit above is pretty standard for GDAL OGR so not sure how to do it otherwise. overall it works fine the way it is, but with that bottle neck it wont make sense to try to move it all to a QGIS plugin which is what i was thinking of doing.

[tracy.love]

i recently had to do some work in redwood canopy canyon conditions with the RS. extremely constrained reception so i used a PPK approach and this code. afterwards i was left with allot of float values that just would not fix no matter what settings i used in RTKLIB. comparing my floats to a LIDAR derived hydrology layer indicated that many were decent but then there would be the occasional crazy point. so i did some redesign and development of this code and i am here to share.

this update follows the same structure as the previous versions, however, now without the bottle neck. so it runs about 3 times faster at least. then it analyzes for velocity between floats to try to find the best floats between RTKLIB processes and then skip the outliers when correcting the rover point file.

caveat: i haven’t retested with processing multiple rover logs against a single base so i may have broken that and i have not tested it at anything but 1hz.

http://public.region3dfg.org/apps/python_geo_processing/R3GIS_RTKlib/

there is a zip there with everything

[tracy.love]

to recap on the function/ strategy

this is a python script to manage processing PPK with RTKLIB and optionally correcting a rover point file with start and stop times.

inputs

1. local RS base raw log file converted to .obs

2. one or more rover raw logs created within the time window of the base file and converted to .obs

3. a folder of RTKLIB .conf files to run each rover raw file through
   (e.g. 2 rover raw logs and 4 .conf files = 8 complete runs to create a combined PPK file)

4. optionally a rover point shapefile with start and stop times which will be corrected against the combined PPK file.

process:

1. the first base raw log + rover raw log is run through RTKLIB and creates a combined PPK file (in shapefile format). the combined PPK file includes two velocity values for each point. the velocity from each previous point and the subsequent point the rover collected. points are given an ID based on the time stamp.

2. each subsequent base raw log + rover raw log run through RTKLIB creates a temporary PPK results file (shapefile) to be compared to the combined PPK file. for a given matching pair between the two layers; if the existing value is float and the new value is fix then it replaces the float. and if the existing is a float and the new float appears better based on velocity constraints then the new one replaces the old one. this develops the combined PPK file, there are attribute fields in the shapefile that track which .conf file produced which points.

3. once all of the rover raw logs have been processed into the combined PPK file you can optionally use the PPK results to correct a rover point file. this rover file could contain RTK corrected points or it could just be autonomous points. the rover file could even come from a completely different app like Collector or QField where you simple mark down start and stop times. points in the combined PPK file that fall between the start and stop times for a given rover point are selected, those outside a maximum velocity are removed, and a weighted point is calculated.

oh well have fun, if anyone has questions ask them here or contact me directly, up to you

[tracy.love]

Being stuck at home means a time for a massive rewrite. new code is basically the same structure but with allot of improvement. bugs fixed, executes faster, and now it calculates velocity, stand deviations, and more. some day it would be good to craft a QGIS plugin for this.

[Pascal_P (Pascal)]

Nice job Tracy, it would be interesting to see where and how you go. Python, Gdal and RTKLib are all three among my preferred tools.

[tracy.love]

Thanks Pascal,

i am also using PyProj and Numpy.

i dont really see a whole lot more to do with it except make a QGIS plugin. which i probably wont get around to doing since it works for me as is. ill continue to tweak it and look for bugs but it feels done to me this time.

working with this really showed me how stable (within say 5cm) results in unconstrained reception areas are when working with the RS+ but also how loose results can get in heavy canopy conifer conditions (like many meters).