<jupyter_start><jupyter_text>Comparison of GPS points within error ellipses**Script prepared by <NAME> - MARUM, University of Bremen**This script compares the distance of two GPS points (Lat/Lon/Height Above Ellipsoid) within their uncertainties.<jupyter_code>import pandas as pd
import geopandas
import matplotlib.patches
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import math
from math import pi<jupyter_output><empty_output><jupyter_text>Input filesInsert the names of the output files. The files must be in *.pos* format, output of RTKlib.<jupyter_code>filename1='28_11_19_BadElf.pos'
filename2='28_11_19_BadElf_CPVG.pos'<jupyter_output><empty_output><jupyter_text>The two files are inserted in a geodataframe and transformed in UTC coordinates. If one uncertainties are null, a minimal uncertainty of 0.0001 is inserted.<jupyter_code>df1 = pd.read_csv(filename1, index_col=None, header=0,skiprows=10,delim_whitespace=True,parse_dates=[['%', 'GPST']])
df2 = pd.read_csv(filename2, index_col=None, header=0,skiprows=10,delim_whitespace=True,parse_dates=[['%', 'GPST']])
points=df1.append(df2)
points.reset_index(drop=True,inplace=True)
gdf = geopandas.GeoDataFrame(points, geometry=geopandas.points_from_xy(points['longitude(deg)'], points['latitude(deg)']),crs={'init':'epsg:4326'})
gdf.to_crs(epsg=3310,inplace=True)
gdf['lon (UTC)'] = gdf.geometry.x
gdf['lat (UTC)'] = gdf.geometry.y
gdf.loc[gdf['sdn(m)'] == 0, 'sdn(m)'] = 0.0001
gdf.loc[gdf['sde(m)'] == 0, 'sde(m)'] = 0.0001
gdf.loc[gdf['sdu(m)'] == 0, 'sdu(m)'] = 0.0001
gdf<jupyter_output><empty_output><jupyter_text>PlotsThe two points are plotted to show how they compare in lat/lon and Ellipsoid Height. Plots include uncertainties, incorporated as 2-sigma from the RTKlib output.<jupyter_code>fig, (ax2, ax1) = plt.subplots(1, 2,figsize=(20,5))
fig.suptitle('Difference between height of points {} (blue) and {} (orange)'.format(filename1,filename2))

llim=min(gdf['height(m)'][0],gdf['height(m)'][1])-4*max(gdf['sdu(m)'][0],gdf['sdu(m)'][1])
ulim=max(gdf['height(m)'][0],gdf['height(m)'][1])+4*max(gdf['sdu(m)'][0],gdf['sdu(m)'][1])

x = np.linspace(llim,ulim, 1000)
ax1.plot(x, stats.norm.pdf(x, gdf['height(m)'][0], gdf['sdu(m)'][0]))
ax1.plot(x, stats.norm.pdf(x, gdf['height(m)'][1], gdf['sdu(m)'][1]))
DeltaH=abs(gdf['height(m)'][0]-gdf['height(m)'][1])
ax1.set(xlabel='Ellipsoid height (m)', title='Delta = {:.2f}m'.format(DeltaH))

f = np.linspace(0, 2*pi, 100)
ax2.plot(gdf['lon (UTC)'][0]+gdf['sde(m)'][0]*np.cos(f) , gdf['lat (UTC)'][0]+gdf['sdn(m)'][0]*np.sin(f))
ax2.plot(gdf['lon (UTC)'][1]+gdf['sde(m)'][1]*np.cos(f) , gdf['lat (UTC)'][1]+gdf['sdn(m)'][1]*np.sin(f))

def distance_cartesian(x1, y1, x2, y2):
    dx = x1 - x2
    dy = y1 - y2
    return np.sqrt(dx * dx + dy * dy)

DeltaPos=(distance_cartesian(gdf['lon (UTC)'][0], gdf['lat (UTC)'][0], gdf['lon (UTC)'][1], gdf['lat (UTC)'][1]))

ax2.set(xlabel='UTM Longitude (m)',ylabel='UTM Latitude (m)',title='Delta = {:.2f}m'.format(DeltaPos))

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2081151832460733, 'fraction_numerical': 0.03926701570680628, 'mean_word_length': 7.991176470588235, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '603236', 'license': 'MIT', 'n_tokens_mistral': 1343, 'n_tokens_neox': 1306, 'n_words': 228, 'path': 'GNSS_data_comparison/.ipynb_checkpoints/Comparing RTKlib points-checkpoint.ipynb', 'repo_name': 'Alerovere/GPS-utilities', 'size': '3059'}