Demo: GNSS Trajectory¶

UW Geospatial Data Analysis
CEE498/CEWA599
David Shean

Pandas Time Series¶

Topics to cover¶

Ingesting and DateTime Index
Resampling
Indexing

import pandas as pd
import geopandas as gpd

fn1='Shean_OASfM_20130909_PosnPnt.csv'
fn2='Shean_OASfM_20130910_PosnPnt.csv'

df = pd.read_csv(fn1)

df

pd.to_datetime(df['GPS Date'] + ' ' + df['GPS Time'])

df = pd.read_csv(fn1, parse_dates=[['GPS Date', 'GPS Time']])

df.columns

df

df.drop(columns=['Latitude.1', 'Longitude.1'], inplace=True)

df.rename(columns={'GPS Date_GPS Time':'GPS DateTime'}, inplace=True)

df.set_index('GPS DateTime', inplace=True)

df.head()

	ID	Latitude	Longitude	HAE	Max PDOP	Max HDOP	Corr Type	GPS Height	Vert Prec	Horz Prec	Std Dev	geometry
GPS DateTime
2013-09-09 20:01:20	1	45.814143	-122.557827	57.791	2.6	1.4	Postprocessed Code	57.791	1.4	0.5	NaN	POINT (-122.55783 45.81414)
2013-09-09 20:01:21	2	45.814157	-122.557834	57.892	2.6	1.4	Postprocessed Code	57.892	1.4	0.6	NaN	POINT (-122.55783 45.81416)
2013-09-09 20:01:22	3	45.814170	-122.557829	53.464	3.3	1.8	Postprocessed Code	53.464	1.6	0.6	NaN	POINT (-122.55783 45.81417)
2013-09-09 20:01:23	4	45.814179	-122.557830	59.237	3.3	1.8	Postprocessed Code	59.237	1.7	0.5	NaN	POINT (-122.55783 45.81418)
2013-09-09 20:01:24	5	45.814199	-122.557822	55.283	10.4	2.5	Postprocessed Code	55.283	2.4	0.6	NaN	POINT (-122.55782 45.81420)

#Define function to do the above wrangling
def parse_trimblexh_export(fn):
    df = pd.read_csv(fn, parse_dates=[['GPS Date', 'GPS Time']])
    df.drop(columns=['Latitude.1', 'Longitude.1'], inplace=True)
    df.rename(columns={'GPS Date_GPS Time':'GPS DateTime'}, inplace=True)
    df.set_index('GPS DateTime', inplace=True)
    return df

df1 = parse_trimblexh_export(fn1)
df2 = parse_trimblexh_export(fn2)

df = pd.concat([df1, df2])

df

import hvplot.pandas

df1.columns

#%matplotlib widget

df1['HAE'].plot();

df1.hvplot()

df.hvplot()

df['2013-09-10 20:26:57':'2013-09-10 20:26:59']

df1['HAE'].hvplot()

gdf1 = gpd.GeoDataFrame(df1, crs='EPSG:4326', geometry=gpd.points_from_xy(df1['Longitude'], df1['Latitude']))

gdf = gpd.GeoDataFrame(df, crs='EPSG:4326', geometry=gpd.points_from_xy(df['Longitude'], df['Latitude']))

gdf1

gdf1.hvplot(aspect='equal')

gdf1_utm = gdf1.to_crs('EPSG:32610')

gdf1_utm.hvplot(aspect='equal', c='HAE')

#https://stackoverflow.com/questions/40452759/pandas-latitude-longitude-to-distance-between-successive-rows

gdf1_utm.head()

	ID	Latitude	Longitude	HAE	Max PDOP	Max HDOP	Corr Type	GPS Height	Vert Prec	Horz Prec	Std Dev	geometry
GPS DateTime
2013-09-09 20:01:20	1	45.814143	-122.557827	57.791	2.6	1.4	Postprocessed Code	57.791	1.4	0.5	NaN	POINT (534352.909 5073492.900)
2013-09-09 20:01:21	2	45.814157	-122.557834	57.892	2.6	1.4	Postprocessed Code	57.892	1.4	0.6	NaN	POINT (534352.393 5073494.477)
2013-09-09 20:01:22	3	45.814170	-122.557829	53.464	3.3	1.8	Postprocessed Code	53.464	1.6	0.6	NaN	POINT (534352.779 5073495.940)
2013-09-09 20:01:23	4	45.814179	-122.557830	59.237	3.3	1.8	Postprocessed Code	59.237	1.7	0.5	NaN	POINT (534352.700 5073496.922)
2013-09-09 20:01:24	5	45.814199	-122.557822	55.283	10.4	2.5	Postprocessed Code	55.283	2.4	0.6	NaN	POINT (534353.303 5073499.172)

Horizontal Velocity (ground speed)¶

Could also look at 3-D vector displacement, vertical velocity, etc.

gdf1_utm.distance(gdf1_utm.shift(1))

gdf1_utm['dm'] = gdf1_utm.distance(gdf1_utm.shift(1))

<ipython-input-115-00434a9b3ef1>:1: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

Left CRS: EPSG:32610
Right CRS: None

  gdf1_utm['dm'] = gdf1_utm.distance(gdf1_utm.shift(1))

gdf1_utm['dm'].max()

gdf1_utm.index

gdf1_utm.index.to_series().diff()

gdf1_utm['dt'] = gdf1_utm.index.to_series().diff().dt.total_seconds()

gdf1_utm['dt']

gdf1_utm['v'] = gdf1_utm['dm']/gdf1_utm['dt']

gdf1_utm['v'].hvplot()

gdf1_utm['a'] = gdf1_utm['v'].diff()

from geoviews import tile_sources as gvts
#map_tiles = gvts.StamenTerrain
map_tiles = gvts.EsriImagery

map_tiles

#kw = {'width':600, 'height':600, 'data_aspect':1, 'alpha':1.0}
kw = {'width':400, 'height':400, 'data_aspect':1, 'alpha':1.0}
kw['colorbar'] = True
kw['cmap'] = 'inferno'
#kw['hover'] = False
#kw['datashade'] = True

Overlay points on map tiles (subset of points)¶

To combine layers in a single plot using holoviews, you use the asterisk (*) to overlay the two objects
- You can use the plus sign (+) to build a layout with two separate subplots
Currently need to use the geo=True option and revert back to the glas_gdf GeoDataFrame with lat/lon geometry, not our reprojected points in glas_gdf_aea.
- In principle, hvplot/geoviews should work with our projected points, but there are some residual issues, and it appears that some underlying code somewhere in the stack is assuming lat/lon

map_tiles * gdf1_utm[::3].to_crs('EPSG:3857').hvplot(c='HAE', size=1, **kw)

gdf1_utm[::3].shape

Using datashader to efficiently render all points¶

map_tiles * gdf1_utm.to_crs('EPSG:3857').hvplot(c='HAE', size=1, datashade=True, **kw)

WARNING:param.Parameterized: Use method 'warning' via param namespace 
WARNING:param.main: There is no reasonable way to use size (or s) with rasterize or datashade. To aggregate along a third dimension, set color (or c) to the desired dimension.

Other useful analysis¶

Distance from starting point
Time spent in different elevation bins
Altitude above ground - sample DEM
Geotagging photos baed on timestamp

University of Washington Geospatial Data Analysis with Python

Demo: GNSS Trajectory¶

Pandas Time Series¶

Topics to cover¶

Horizontal Velocity (ground speed)¶

Overlay points on map tiles (subset of points)¶

Using datashader to efficiently render all points¶

Other useful analysis¶