<jupyter_start><jupyter_text>Problem: Sourcing Fine-Grained Mobility Data from Large Datasets OverviewFinding fine-grained, agent-specific mobility data is very difficult to do. It is extremely important to see how an individual acts, as we can see from various clusters of the COVID-19 pandemic that we have experienced in New Zealand. It would be extremely useful to analyze the mobility of individual people, finding their clusters of movement. A more in-depth overview of the problem, with supporting literature, can be found in the accompanying report. The DataMobility data will be obtained from the GeoLife dataset. This data is available for free [here](https://www.microsoft.com/en-us/download/details.aspx?id=52367). It was collected by Microsoft Asia and consists of GPS trajectories from 182 users over a period of three years. You can download it and load it through the cells that follow. This system depends on data in a PostGIS server - don't worry, this will set that up for you, you will just need a PgAdmin server running with default 'postgres' credentials. The Proposed SolutionA user-friendly analysis system for looking at individual users of the dataset and finding patterns. This system should be approachable, web-based and capable of displaying clusters of mobility. It will also deliver key statistics on the data so that the user can validate the results they are seeing. To this end, it will be open-source. Structure**This Notebook consists of two parts**. The first series of cells go through the project requirement checklist, running each step individually with some maps along the way. They show how the data can be uploaded to PgAdmin etc. The final cell is a 'dashboard' that does all of the analysis of the previous steps, all contained. Code in this final cell is not very well commented as I have done this thoroughly in the cells that preceed it. There are also bugs that exist in the following cells that I have fixed for the final version, so please judge program functionality through the dashboard. To be clear, **sequential cells outline basic functionality and have known bugs that have been fixed for the dashboard, so please run the final cell to assess program functionality** Thanks :)<jupyter_code># import modules                                # used for:
import psycopg2                                 # sql implementation
from sqlalchemy import create_engine            # connection to postgis
import pandas as pd                             # tabular data handling
import geopandas as gpd                         # tabular data handling with geometry (spatial)
import folium                                   # map creation (build on leaflet)
import geoalchemy2                              # support for spatial sql
import matplotlib.pyplot as plt                 # plotting tool
import numpy as np                              # used for getting statistics and complex numerical analyses
from ipyleaflet import *                        # ipython leaflet extension
from ipywidgets import *                        # ipython widgets (buttons, etc.)
import glob                                     # used for getting path names (for finding the geolife folder)
from IPython.display import display             # display function for inline plots etc.
import pandas as pd                             # dataframes
import datetime                                 # getting datetime from delimited data
import os                                       # use local functions and getting system info
import psutil                                   # getting RAM of computer for performance recommendations
import seaborn as sns                           # matplotlib, but fancy
import random                                   # random number generation
from sklearn.cluster import DBSCAN              # clustering algorithm
from sklearn.neighbors import NearestNeighbors  # neigbourhoods used for clustering
from kneed import KneeLocator                   # locating knee point of clustering algorithm
# define constants
SHOW_EXCEPTIONS = True                                                      # bool for debugging exceptions
DATA_EXISTS = True                                                          # bool for not reloading data
TB_RAW = 'geolife_raw_data'                                                 # constant table name
TB_CLEAN = 'geolife_ref'                                                    # ""
TB_BEIJING = 'beijing_raw'                                                  # ""
ENGINE = create_engine('postgresql://postgres:postgres@localhost/postgres') # postgres engine
CON = ENGINE.connect()                                                      # engine connection
CRS = '4326'                                                                # coordinate reference system     
BEIJING = [39.9, 116.41]                                                    # central beijing coords, for map centres
B1 = 115.779419, 39.626846                                                  # bbox limits for beijing extent
B2 = 116.952209, 40.357010                                                  # ""

# method for running certain queries and handling exceptions
def runsql(query):
    try:
        ENGINE.execute(query)
    except Exception as exep:
        if SHOW_EXCEPTIONS:
            print(exep)
# The following code was modified by code developed by HERE Technologies 
# The original can be found: https://heremaps.github.io/pptk/tutorials/viewer/geolife.html
# My edits:
# -made it a class, added getters and setters (to meet assignment guidelines)
# -added customisation by user (eg. custom file locations without constants)
# -general cleanup of code, added nicer progress indicator, commented
class GeoLifeHandler:
    def __init__(self):
        # apply numeric values for transport modes to save space
        self.mode_names = ['walk', 'bike', 'bus', 'car', 'subway','train', 'airplane', 'boat', 
                           'run', 'motorcycle', 'taxi']
        self.mode_ids = {s : i + 1 for i, s in enumerate(self.mode_names)}
        
        # create progress bar
        self.progress_bar = widgets.FloatProgress(
            value=0,
            min=0,
            max=10.0,
            description='Loading:',
            bar_style='info',
            style={'bar_color': '#0000FF'},
            orientation='horizontal')
    
    # read plt file
    def read_plt(self, plt_file):
        points = pd.read_csv(plt_file, skiprows=6, header=None,
                             parse_dates=[[5, 6]], infer_datetime_format=True)
        points.rename(inplace=True, columns={'5_6': 'time', 0: 'lat', 1: 'lon', 3: 'alt'}) # rename cols
        points.drop(inplace=True, columns=[2, 4]) # remove unused columns
        return points

    # read labels (where applicable)
    def read_labels(self, labels_file):
        labels = pd.read_csv(labels_file, skiprows=1, header=None,
                             parse_dates=[[0, 1], [2, 3]],
                             infer_datetime_format=True, delim_whitespace=True)
        labels.columns = ['start_time', 'end_time', 'label']
        labels['label'] = [self.mode_ids[i] for i in labels['label']] # enumerate
        return labels
            
    # apply travel mode labels to points (enumerated to save bits)
    def apply_labels(self, points, labels):
        indices = labels['start_time'].searchsorted(points['time'], side='right') - 1
        no_label = (indices < 0) | (points['time'].values >= labels['end_time'].iloc[indices].values)
        points['label'] = labels['label'].iloc[indices].values
        points['label'][no_label] = 0
            
    # read individual user (folder is preassigned to each user)
    def read_user(self, user_folder):
        labels = None
        plt_files = glob.glob(os.path.join(user_folder, 'Trajectory', '*.plt'))
        df = pd.concat([self.read_plt(f) for f in plt_files])
        labels_file = os.path.join(user_folder, 'labels.txt')
        if os.path.exists(labels_file):
            labels = self.read_labels(labels_file)
            self.apply_labels(df, labels)
        else:
            df['label'] = 0
        return df
    
    # interate through all users (defined by folders in Data folder)
    def read_all_users(self):
        subfolders = os.listdir(self.folder)
        dfs = []
        self.progress_bar.max = len(subfolders)
        display(self.progress_bar)
        for i, sf in enumerate(subfolders):
            #print('processing user {} of {}'.format(i + 1, len(subfolders)))
            self.progress_bar.value = i + 1
            df = self.read_user((os.path.join(self.folder,sf)))
            df['user'] = int(sf)
            dfs.append(df)
        print('Load Complete')
        self.geolife_df = pd.concat(dfs)
    
    # set folder location
    def set_location(self, file):
        self.folder = file
        
    # get dataframe
    def get_df(self):
        return self.geolife_df

    
# function for submitting file location
# starts the reading of all users, saving to geodataframe
def bt_submitevent(b):
    if os.path.exists(file_options.value):
        b.disabled = True
        handler.set_location(file_options.value)
        handler.read_all_users()
    else:
        print('Data file not found. Please try again...')

# submit button for selecting the option, calling the above function when selected
bt_submit = widgets.Button(description='Select') 
bt_submit.on_click(bt_submitevent)

# list of file options
file_options = widgets.Dropdown(
    options=[(os.path.join(os.getcwd(), 'geolife_trajectories\Data')),                   # default zip file location
             (os.path.join(os.path.join(os.environ['USERPROFILE']), 'Desktop\Data')),    # desktop folder
             (os.path.join(os.path.join(os.environ['USERPROFILE']), 'Downloads\Data'))], # downloads folder
    disabled=False,
    layout={'width': '70%'})

# display these widgets, with title
display(VBox([HTML('<h4>GeoLife File Location</h4>'), HBox([file_options, bt_submit])]))

# create instance of handler class
handler = GeoLifeHandler()
# try and get the dataframe information
# good test if the dataframe exists. if it doesn't, assume it is already available on server (see constants)
try:
    geolife_raw = handler.get_df()
    display(geolife_raw.count())
    display(geolife_raw.head(10))
    geolife_raw = geolife_raw.head(10)
except:
    print('Data either not yet extracted, or already available on server')
# class for uploading to postgis and setting data limitations based on memory available
class PostGISUploader:
    def __init__(self):
        # use constants for table names etc.
        self.table_rawdata = TB_RAW
        self.table_cleandata = TB_CLEAN
        self.crs = CRS  
        
        # buttons for setting limitations on heatmap
        self.buttons = [
            widgets.Button(description='100,000'),
            widgets.Button(description='1,000,000'),
            widgets.Button(description='10,000,000'),
            widgets.Button(description='No Limitation')
        ] 
        
        # assign the set_limit function to all buttons (will change output based on value)
        for button in self.buttons:
            button.on_click(self.set_limit)
        self.limit = 'LIMIT 100' # default limit
        
        # get ram of computer
        self.pc_ram = round(psutil.virtual_memory().total / 1073741824, 1)
     
    # upload data to postgis. THIS TAKES A VERY LONG TIME (mine took ~14 hours)
    def upload(self, geolife_df):        
        # create 'raw data' table. this will include all attributes of the dataset
        runsql('DROP TABLE IF EXISTS {}'.format(self.table_rawdata))
        runsql(('CREATE TABLE {}'
               '(gid serial PRIMARY KEY, time timestamp, lat float, lon float, alt float, label int, user_id int)'
                .format(self.table_rawdata)))

        # data must be inserted row by row. using 'to_postgis' will crash as it uses too much memory
        for row in geolife_df.iterrows():
            cur_row = row[1]
            runsql(("INSERT INTO {} (time, lat, lon, alt, label, user_id) VALUES ('{}', {}, {}, {}, {}, {})"
                     .format(self.table_rawdata, cur_row['time'], cur_row['lat'], 
                             cur_row['lon'], cur_row['alt'], cur_row['label'], cur_row['user'])))

        # create 'clean' table. will treat the previous table as a backup
        # this table has a geometry ref instead of lat/lon
        runsql('DROP TABLE IF EXISTS {}'.format(self.table_cleandata))
        runsql(('CREATE TABLE {} AS SELECT gid, time, st_setsrid(st_makepoint(lon, lat), {}) '
            'as geom, label, user_id FROM {}'.format(self.table_cleandata, self.crs, self.table_rawdata)))
    
    # disable buttons, get limitation based on the button value (with some string formatting)
    def set_limit(self, b):
        for button in self.buttons:
            button.disabled = True
        self.limit = b.description.replace(',', '')
        if self.limit == 'No Limitation':
            self.limit = ''
        else:
            self.limit = 'LIMIT ' + self.limit
    
    # get a message based on the amount of RAM your computer has
    def get_limitoptions(self):
        if self.pc_ram < 4:
            message = 'This is quite low (no offense). Recommended: 100,000 row limitation'
        elif 4 < self.pc_ram < 8:
            message = 'This is acceptable. Recommended: 1,000,000 row limitation'
        elif 8 < self.pc_ram < 16:
            message = 'This is good, but still limited. Recommended: 10,000,000 row limitation'
        else:
            message = 'You can handle anything. Try going with the whole dataset!'
        return message
    
    # return buttons
    def get_buttons(self):
        return self.buttons
    
    # return limit suffix
    def get_limit(self):
        return self.limit
    
    # return pc ram
    def get_ram(self):
        return self.pc_ram

    
# create object instance. this is needed even if data is uploaded already to get limitation and ram etc.
uploader = PostGISUploader()

# start the data upload process if the data has not already been uploaded (constant = false)
if not DATA_EXISTS:
    uploader.upload(geolife_raw)
# class for dealing with folium maps
class FoliumCreator:
    # set default values
    def __init__(self, zoom=10, size='600px', center=BEIJING, scale_pos='topleft'):
        self.zoom = zoom
        self.size = size
        self.center = center # using American english spelling as that's what leaflet uses, not bc I like it :(
        self.layout = Layout(height=self.size)
        self.scale_pos = scale_pos
    
    # set map layout
    def set_layout(self, layout):
        self.layout = layout
    
    # set zoom level
    def set_zoom(self, level):
        self.zoom = level
    
    # return map object
    def get_map(self):
        m = Map(center=self.center, zoom=self.zoom, layout=self.layout, basemap=basemaps.Stamen.Terrain)
        self.set_control(m, ScaleControl(position=self.scale_pos))
        return m
    
    # set control (eg. scale bar)
    def set_control(self, instance, control):
        instance.add_control(control)
    
    # set new layer
    def set_layer(self, instance, layer):
        instance.add_layer(layer)<jupyter_output><empty_output><jupyter_text>Showing our Area of Interest A map of Beijing without any data. Basemaps show the topography and cityscape. We have chosen Beijing as most of the GeoLife data has telemetry from this city.<jupyter_code>'''
MAP 1: Area of Interest
-no context, just area
-zoomed into centre based on constants
-uses two basemaps to show urban and environment
'''

# create mapcontroller instance, return a map object
mapcontroller = FoliumCreator()
map1 = mapcontroller.get_map()

# add basemap layers
right_layer = basemap_to_tiles(basemaps.OpenTopoMap)   # OpenTopMap
left_layer = basemap_to_tiles(basemaps.Stamen.Terrain) # Stamen Terrain (my fav)

# add splitmap controller to map
split_control = SplitMapControl(left_layer=left_layer, right_layer=right_layer)
mapcontroller.set_control(map1, split_control)

# show map
display(map1)<jupyter_output><empty_output><jupyter_text>Mapping the Raw DataNow that we have uploaded all of our data to a server we can safely manipulate it. Again, we will focus on Beijing, China as this is where the majority of the data originate from. We will create a table that grabs all points that are within our defined 'bounding box'. Then, we can produce a heatmap of these points, coloured by density of points. This will show where most people moved around in Beijing on the dataset (or at least, where their phones reported location from the most). Note that despite limiting the extent, the data is still over 24 million rows, so this may take some time to process!<jupyter_code># print information and recommendation based on ram 
print(('This dataset is over 24 million rows. It requires high amounts of memory to process in its entirety. '
      'Detected {} gigabytes, not all of which will be dedicated to the interpreter. {}. Of course, this is a free'
      ' world, do what you want - but you have been warned.'
       .format(uploader.get_ram(), uploader.get_limitoptions())))

# add buttons from uploader object, display buttons
buttons = uploader.get_buttons()
display(VBox([HTML('<h4>Please select a data row limitation</h4>'), 
              HBox([buttons[0], buttons[1], buttons[2], buttons[3]])])) 
'''
MAP 2: Movement Heatmap
-focus on Beijing, create table of points within that extent
-saves table in cache, which will use a fair bit of storage on your device (depending on limitation)
'''  
# create beijing table based on constant
if not DATA_EXISTS:
    # drop if exists, just in case
    runsql('DROP TABLE IF EXISTS {}'.format(TB_BEIJING))
    print('Existing table dropped (if it existed). Creating Beijing table. This can take several minutes...')
    
    # selects beijing area from a makeenvelope function (clip to extent essentially)
    runsql('CREATE TABLE {} AS SELECT * FROM {} WHERE geom && st_makeenvelope ({}, {}, {}, {}, {})'
           .format(TB_BEIJING, TB_CLEAN, B1[0], B1[1], B2[0], B2[1], CRS)) 
    print('Table created.') # giving info like this as it can take a while, so good to show something is happening

# get beijing table
print('Querying server. This can take a while...')
query = ('SELECT * FROM {} {}'.format(TB_BEIJING, uploader.get_limit()))

# get geodataframe for heatmap (all beijing)
gdf = gpd.GeoDataFrame.from_postgis(query, CON)
print(gdf.count())

# get map, changing zoom
mapcontroller.set_zoom(10)
map2 = mapcontroller.get_map()

# convert geometry from geodataframe into a list so that it can be used for ipyleaflet heatmap
print('GeoDataFrame created. Converting values to produce heatmap...')
geom_list = [(x, y) for x, y in zip(gdf['geom'].y, gdf['geom'].x)]
heatmap = Heatmap(locations=geom_list, radius=8)

# add heatmap to map object, display map
map2.add_layer(heatmap)
display(map2)<jupyter_output>Querying server. This can take a while...
gid        19567865
time       19567865
geom       19567865
label      19567865
user_id    19567865
dtype: int64
GeoDataFrame created. Converting values to produce heatmap...<jupyter_text>Finding UsersWe can see that, even when restricting to our city of interest, the dataset is far too big to handle. Let's restrict it further and focus in on a single user to have a look at the fine-grained patterns. Initially, I wrote an algorithm to find trips based on the label of transport changing, which worked by:- iterating through every user, grabbing all of their points from the server- checking if they have labels assigned to their travels - calculating their unique 'trips', assigned by the amount of times the label assigned to these trips changes, eg. walking for 10 minutes followed by a taxi ride will be defined as two trips- getting a user who has more than a defined minimum trip count, then using binary search to find the user who has the lowest amount of trips above this minimumHowever, this method has several flaws. As mentioned, not all users have labels assigned for their travels. Another is that two modes of transport may be considered part of a single trip. Instead, we can focus on user-level data, rather than trip-level data. We can focus on clusters of points rather than alternating transport methods, thus removing the reliance on the existence of labels (the algorithm I wrote found that over 70% of users have to labels). Overall, we can see that user-level data will be more effective than trip-level data as a way to handle the data. First, let's select a user to focus on. We will write the analysis in a way that it will work on any user that exists in the dataset.<jupyter_code># class for each individual user
# instead of creating a new class instance each time a new user is chosen, this function can update based on user id
# handles all the displays for interacting with users, as well as giving information to analysis tools
class UserObject:
    
    # initialize with some defaults
    def __init__(self, dashboard_mode=False):
        self.id = -1
        self.id_count = 181
        self.dashboard_mode = dashboard_mode # depends on whether used in final version or inline
        
        # submit user button
        self.bt_submit_user = widgets.Button(description='Select')
        self.bt_submit_user.on_click(self.submit_user)

        # randomize user button
        self.bt_random_user = widgets.Button(description='Randomize')
        self.bt_random_user.on_click(self.randomize_user)

        # get random id for the first user
        self.first_user_id = self.get_random_id()

        # dropdown for all users
        self.user_options = widgets.Dropdown(
            value=self.first_user_id,
            options=self.get_all_ids(),
            disabled=False,
            layout={'width': '30%'})
        # observed function for change of dropdown
        self.user_options.observe(self.user_change, 'value')

        # outputs 
        self.dp_userinfo = widgets.Output(layout={'border': 'solid black 1px',
                                                  'max_width': '360px'})
        self.clusterbox = widgets.Output(layout={'border': 'solid black 1px',
                                                 'max_width': '200px',
                                                 'min_width': '200px'})
        
        # get user info for first user info
        self.get_user_info(self.first_user_id)
        
        # disables submit button for dashboard mode
        if self.dashboard_mode:
            self.bt_submit_user.disabled = True
    
    # display all widgets
    def display(self):
        display(HBox([HTML('<b>User:  '), self.user_options, self.bt_random_user, self.bt_submit_user]))
        display(HBox([self.dp_userinfo, widgets.Output(layout={'min_width': '75px'}), self.clusterbox]))

    # submit user (non-dashboard). disables buttons, sets id based on selected value
    def submit_user(self, b): # 'b' acts as the button event/container
        b.disabled = True
        self.bt_random_user.disabled = True
        self.user_options.disabled = True
        self.set_id(self.user_options.value)
        print('User {} selected!'.format(self.get_id()))
    
    # get random user
    def randomize_user(self, b):
        self.user_options.value = self.get_random_id()

    # change user to value in dropdown (observing dropdown widget)
    def user_change(self, value):
        self.get_user_info(value.new)
        if self.dashboard_mode:
            self.set_id = value.new
            
    # set id
    def set_id(self, new_id):
        self.id = new_id
        
    # return id
    def get_id(self):
        return self.id
    
    # set what to display in 'clusterbox' (print output on the right to show clustering info)
    def set_clusterbox(self, content):
        with self.clusterbox:
            print(content)
    
    # get all possible id values
    def get_all_ids(self):
        return [i for i in range(0, self.id_count)]
    
    # get random id (random number generation)
    def get_random_id(self):
        return random.randint(0, self.id_count)
    
    # get user info (through db queries)
    def get_user_info(self, test_id):
        # show loading process
        self.bt_random_user.description = 'Loading...'
        self.bt_random_user.disabled = True
        
        # get info 
        info = CON.execute('select count(*), count(distinct(label)), min(time), max(time) from {} where user_id = {}'
                           .format(TB_BEIJING, test_id)).first()
        
        # clear any existing info
        self.dp_userinfo.clear_output()
        
        # print formatted summary based on query
        with self.dp_userinfo:
            print('User {} Summary\n--------------------------------------------'.format(test_id))
            print(('Data points:             {}\nUnique Travel Methods:   {}\nData Begin:              {}'
                  '\nData End:                {}'.format(info[0], info[1], info[2], info[3])))
        
        # show loading process
        self.bt_random_user.description = 'Randomize'
        self.bt_random_user.disabled = False
# create user object, display widgets
user = UserObject()
user.display()<jupyter_output><empty_output><jupyter_text>The code below looks at individual user movement and filters these movements to make this fine-grained data easier to handle. We can look at their mobility in a less processing-intensive way by using the function ST_REMOVEREPEATEDPOINTS which does exactly what the name suggests. Repeated points are defined by a tolerance level (essentially a buffer around each point):<jupyter_code># get user information. query is based on collection of geometry (merge)
gdf_raw = gpd.GeoDataFrame.from_postgis((('SELECT st_collect(geom) as all_geom, st_npoints(st_collect(geom)) '
                                     'as point_count from {} where user_id = {}')
             .format(TB_BEIJING, user.get_id())), CON, geom_col='all_geom')

# give up if returns as none, giving some error info
gdf_raw_n = str(gdf_raw['point_count'][0])
if gdf_raw_n == 'None':
    print("No points found. Please ensure you have clicked 'Select' above")
else:
    print('Raw point count: ' + gdf_raw_n)
    
    # get simplified user information. removes repeated points as a simplification of the total user's mobility
    gdf_simple = gpd.GeoDataFrame.from_postgis((('SELECT st_removerepeatedpoints(st_collect(geom)) as simplified_geom, '
                'st_npoints(st_removerepeatedpoints(st_collect(geom), 0.001)) as point_count from {} where user_id = {}')
                .format(TB_BEIJING, user.get_id())), CON, geom_col='simplified_geom')
    gdf_simple_n = str(gdf_simple['point_count'][0])
    print('Simplified point count: ' + gdf_simple_n)

    # plot output
    fig, axes = plt.subplots(1, 2, figsize=(15, 10))
    
    # raw mobility trajectories
    axes[0].set_title('Original Trip Data for User {} (n={})'.format(user.get_id(), gdf_raw_n))
    gdf_raw.plot(ax=axes[0])
    axes[0].set_xlabel('Longitude')
    axes[0].set_ylabel('Latitude')

    # simplified mobility trajectories
    axes[1].set_title('Simplified Trip Data for User {} (n={})'.format(user.get_id(), gdf_simple_n))
    gdf_simple.plot(ax=axes[1])
    axes[1].set_xlabel('Longitude')
    axes[1].set_ylabel('Latitude')

    # show plots
    plt.show()<jupyter_output>Raw point count: 81184
Simplified point count: 2963<jupyter_text>Validation of SimplificationThese two plots look very similar, but we yet can't be sure. We can confirm by testing that the simplified analysis falls entirely within the area of the raw dataset. We will also plot the two datasets together to make sure the movement trends are the same. This takes a bit of time as we need to perform the simplification, split the geometery into latitude and longitude so we can plot it and perform our ST_CONTAINS.Note that the plots will show some distortion as they are not assigned a coordinate reference system. For the purposes of this analysis, this does not matter as we are only concerned about their relative (as opposed to absolute) displacement.<jupyter_code># get raw data, split as x and y for plotting
gdf_raw = gpd.GeoDataFrame.from_postgis((('SELECT st_x((st_dumppoints(st_collect(geom))).geom) as lat, '
        'st_y((st_dumppoints(st_collect(geom))).geom) as lon, (st_dumppoints(st_collect(geom))).geom as geom '
        'FROM {} WHERE user_id = {}'.format(TB_BEIJING, user.get_id()))), CON)
# get simplified data, also split into x and y
gdf_simple = gpd.GeoDataFrame.from_postgis((('SELECT '
        'st_x((st_dumppoints(st_removerepeatedpoints(st_collect(geom), 0.001))).geom) as lat, '
        'st_y((st_dumppoints(st_removerepeatedpoints(st_collect(geom), 0.001))).geom) as lon, '
        '(st_dumppoints(st_removerepeatedpoints(st_collect(geom), 0.001))).geom as geom '
        'FROM {} WHERE user_id = {}'.format(TB_BEIJING, user.get_id()))), CON)

# tests if the simplified data is completely contained by the raw data
contains = CON.execute(('SELECT st_contains(st_collect(geom), st_removerepeatedpoints(st_collect(geom), 0.001)) '
                        'FROM {} where user_id = {}'.format(TB_BEIJING, user.get_id()))).first()[0]

# display some information. gives validation based on what the above query returns
print(('\nNote that size of points is offset for aesthetics and has no meaning.'
      'The ST_CONTAINS was run on the raw data. Also note that the difference in n between this and the previous '
      'plot is due to rounding and is inconsequential to our analysis\n'))
if contains:
    print('Raw data completely contains simplified data.')
else: 
    print('Raw data does not contain simplified data.')

# create new plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)
fig.figsize = (10, 10)

# plot raw and simplified datasets atop one another
ax.scatter(gdf_raw['lat'], gdf_raw['lon'], s=50, c='black', marker="o", label='Raw (n={})'
           .format(gdf_raw.count()[0])) # larger dots to show difference
ax.scatter(gdf_simple['lat'], gdf_simple['lon'], s=2, c='orange', marker="o", label='Simplified (n={})'
           .format(gdf_simple.count()[0]))
ax.set_ylabel('Latitude')
ax.set_xlabel('Longitude')
plt.legend(loc='lower left');

# show plot
plt.show()<jupyter_output>Note that size of points is offset for aesthetics and has no meaning.The ST_CONTAINS was run on the raw data. Also note that the difference in n between this and the previous plot is due to rounding and is inconsequential to our analysis

Raw data completely contains simplified data.<jupyter_text>Pattern of User MovementsWe now have a set of telemetry from an individual user, compressed and validated to be more manageable. Let's do something with this data!   Density-Based Spatial Clustering of Applications with Noise (DBSCAN)DBSCAN will help us find clusters in the data. These clusters will unveil information about the mobility of the user, showing patterns of movement that are clustered together which we can build inferences from. The sklearn package has all the tools we need to run a DBSCAN analysis, starting with calculating the distance of nearest-neighbour points. We can use these distances to find the 'knee point' on the 'elbow diagram' - the point where the distances of neighbours greatly increases (almost) asymptotically. This 'knee point' serves as the epsilon value in the DBSCAN analysis, where we can then classify points into their respective clusters.<jupyter_code># set sample parameter, get list of geometries, create figures for output
samples = 20
geom_list = [(x, y) for x, y in zip(gdf_simple['lat'], gdf_simple['lon'])]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# get nearest neigbours
nn = NearestNeighbors(n_neighbors=samples + 1)
neighbors = nn.fit(geom_list)
distances, indices = neighbors.kneighbors(geom_list)
distances = np.sort(distances[:,samples], axis=0)

# plot nearest neighbours
axes[0].set_title('Nearest Neighbours Plot on User {} Mobility'.format(user.get_id()))
axes[0].plot(distances)
axes[0].set_xlabel('Points')
axes[0].set_ylabel('Distances')

# calculate knee point
i = np.arange(len(distances))
knee = KneeLocator(i, distances, S=1, curve='convex', direction='increasing', interp_method='polynomial')

# plot nearest neighbours with knee point
axes[1].set_title('Knee Point of Nearest Neighbours for User {} Mobility'.format(user.get_id()))
axes[1].plot(distances)
axes[1].axvline(x=knee.knee,color='gray',linestyle='--')
axes[1].text((axes[1].get_xlim()[1] * 0.5), (axes[1].get_ylim()[1] * 0.5), 'knee={}'.format(knee.knee))
axes[1].set_xlabel('Points')
axes[1].set_ylabel('Distances')

# create dbscan instance
db = DBSCAN(eps=distances[knee.knee], min_samples=samples).fit(geom_list)

# get labels (clustered)
labels = db.labels_

# plot figure as scatterplot (seaborn)
fig = plt.figure(figsize=(8, 8))
sns_scatter = sns.scatterplot(x=gdf_simple['lat'], y=gdf_simple['lon'], marker="o", 
            hue=["{cluster}".format(cluster = 'Cluster ' + 
                    str(cur_label) if cur_label != -1 else 'Noise') for cur_label in labels])
sns_scatter.set(xlabel='Latitude', ylabel='Longitude', 
                title='Clusters in User {} Mobility using DBSCAN\n(where -1 is noise)'.format(user.get_id()))

# show plots
plt.show()<jupyter_output><empty_output><jupyter_text>Analyzing Mobility Patterns/ClustersOf course, the simplification will greatly impact the reliability of the clustering. However, it is not viable to load the entire dataset for every user in-situ for clustering as they are far too large. This is a solution that compromises processing time and quality of output. We now have clustered mobility patterns for a specific user. What information can we derive from this? First we will need to consolidate everything about the specific user into one table:<jupyter_code># create dataframe for manipulation, add clustering labels
gdf_user = gdf_simple
gdf_user['label'] = labels
#gdf_user['time'] = 'datetime'

# count with noise
total_count = gdf_user.count()[0]

# get locations of noise
noise_index = gdf_user[gdf_user['label'] == -1].index  

# get counts, drop noise
gdf_user_clean = gdf_user
gdf_user_clean.drop(noise_index, inplace=True)
clean_count = gdf_user_clean.count()[0]
noise_count = total_count - clean_count

# print info of noise, clusters, counts before and after dropping noise
print('Out of {} points, {} were detected as noise ({}%). With noise removed, {} points remain with {} cluster(s)'
     .format(total_count, noise_count, round(noise_count / total_count * 100, 2), clean_count, 
             len(gdf_user['label'].unique())))<jupyter_output>Out of 2963 points, 11 were detected as noise (0.37%). With noise removed, 2952 points remain with 9 cluster(s)<jupyter_text>Visualizing Mobility PatternsWe now have a clean, simplifed dataset. GPS trajectories will be plotted along paths so that we can see exactly where the user moves within clusters.<jupyter_code># Antpath plot
# here we can see the user's mozement along a moving path
colours = ['red', 'green', 'blue', 'orange', 'purple', 'pink', 'brown', 'black', 'cyan', 'gold', 'white', 'crimson', 'olive']
mapcontroller.set_zoom(9)
map2 = mapcontroller.get_map()
                   
by_cluster = gdf_user_clean.groupby('label')
by_cluster = dict(list(by_cluster))

for cluster in by_cluster.keys():
    current_gdf = by_cluster[cluster]
    geom_list = [(x, y) for x, y in zip(current_gdf['lon'], current_gdf['lat'])]
    
    ant_path = AntPath(locations=geom_list,
                      dash_array=[1, 10],
                      delay=3000,
                      color=colours[cluster],
                      pulse_color='black')

    mapcontroller.set_layer(map2, ant_path)

display(map2)<jupyter_output><empty_output><jupyter_text>User InterfaceWe have now outlined the methods required to get fine-grained mobility analysis on this GeoLife dataset. We will finish by combining these methods into a single interface so that those less familiar with GISytems can access this data.This following cell, aside from the collation of data and the requirement of third-party packages, should function independently from the rest of the Notebook. This means it will expect an existing dataset uploaded to PostGIS, which, if you have run the previous cells, should be satisfied.<jupyter_code>'''
GeoLife Mobility Dashboard

-allows the selection of any user, displaying their DBSCAN clusters and trajectories
-based on code above, with a lot of class modification (copied so that it can be run independently))
-for in-depth commenting, see the cells above
-can be run independently (along with imports)
'''
# import modules                                # used for:
import psycopg2                                 # sql implementation
from sqlalchemy import create_engine            # connection to postgis
import pandas as pd                             # tabular data handling
import geopandas as gpd                         # tabular data handling with geometry (spatial)
import folium                                   # map creation (build on leaflet)
import geoalchemy2                              # support for spatial sql
import matplotlib.pyplot as plt                 # plotting tool
import numpy as np                              # used for getting statistics and complex numerical analyses
from ipyleaflet import *                        # ipython leaflet extension
from ipywidgets import *                        # ipython widgets (buttons, etc.)
import glob                                     # used for getting path names (for finding the geolife folder)
from IPython.display import display             # display function for inline plots etc.
import pandas as pd                             # dataframes
import datetime                                 # getting datetime from delimited data
import os                                       # use local functions and getting system info
import psutil                                   # getting RAM of computer for performance recommendations
import seaborn as sns                           # matplotlib, but fancy
import random                                   # random number generation
from sklearn.cluster import DBSCAN              # clustering algorithm
from sklearn.neighbors import NearestNeighbors  # neigbourhoods used for clustering
from kneed import KneeLocator                   # locating knee point of clustering algorithm

DATA_EXISTS = True
SHOW_EXCEPTIONS = True
TB_RAW = 'geolife_raw_data'
TB_CLEAN = 'geolife_ref'
ENGINE = create_engine('postgresql://postgres:postgres@localhost/postgres')
CON = ENGINE.connect()
CRS = '4326'
BEIJING = [39.9, 116.41]
TB_BEIJING = 'beijing_raw'
B1 = 115.779419, 39.626846
B2 = 116.952209, 40.357010
COLOURS = ['green', 'blue', 'orange', 'pink', 'yellow', 'purple', 'red', 'brown', 
           'gray', 'aqua', 'peru', 'gold', 'crimson', 'olive', 'white']

def runsql(query):
    try:
        ENGINE.execute(query)
    except Exception as exep:
        if SHOW_EXCEPTIONS:
            print(exep)

class GeoLifeHandler:
    def __init__(self):
        # apply numeric values for transport modes to save space
        self.mode_names = ['walk', 'bike', 'bus', 'car', 'subway','train', 'airplane', 'boat', 
                           'run', 'motorcycle', 'taxi']
        self.mode_ids = {s : i + 1 for i, s in enumerate(self.mode_names)}
        self.progress_bar = widgets.FloatProgress(
            value=0,
            min=0,
            max=10.0,
            description='Loading:',
            bar_style='info',
            style={'bar_color': '#0000FF'},
            orientation='horizontal')
        self.folder = 'undefined'
        self.geolife_df = []
        self.limit = ''
    
    # read plt file
    def read_plt(self, plt_file):
        points = pd.read_csv(plt_file, skiprows=6, header=None,
                             parse_dates=[[5, 6]], infer_datetime_format=True)
        points.rename(inplace=True, columns={'5_6': 'time', 0: 'lat', 1: 'lon', 3: 'alt'}) # rename cols
        points.drop(inplace=True, columns=[2, 4]) # remove unused columns
        return points

    # read labels (where applicable)
    def read_labels(self, labels_file):
        labels = pd.read_csv(labels_file, skiprows=1, header=None,
                             parse_dates=[[0, 1], [2, 3]],
                             infer_datetime_format=True, delim_whitespace=True)
        labels.columns = ['start_time', 'end_time', 'label']
        labels['label'] = [self.mode_ids[i] for i in labels['label']] # enumerate
        return labels
            
    # apply travel mode labels to points (enumerated to save bits)
    def apply_labels(self, points, labels):
        indices = labels['start_time'].searchsorted(points['time'], side='right') - 1
        no_label = (indices < 0) | (points['time'].values >= labels['end_time'].iloc[indices].values)
        points['label'] = labels['label'].iloc[indices].values
        points['label'][no_label] = 0
            
    # read individual user (folder is preassigned to each user)
    def read_user(self, user_folder):
        labels = None
        plt_files = glob.glob(os.path.join(user_folder, 'Trajectory', '*.plt'))
        df = pd.concat([self.read_plt(f) for f in plt_files])
        labels_file = os.path.join(user_folder, 'labels.txt')
        if os.path.exists(labels_file):
            labels = self.read_labels(labels_file)
            self.apply_labels(df, labels)
        else:
            df['label'] = 0
        return df
    
    # interate through all users (defined by folders in Data folder)
    def read_all_users(self):
        subfolders = os.listdir(self.folder)
        dfs = []
        self.progress_bar.max = len(subfolders)
        display(self.progress_bar)
        for i, sf in enumerate(subfolders):
            #print('processing user {} of {}'.format(i + 1, len(subfolders)))
            self.progress_bar.value = i + 1
            df = self.read_user((os.path.join(self.folder,sf)))
            df['user'] = int(sf)
            dfs.append(df)
        print('Load Complete')
        self.geolife_df = pd.concat(dfs)
    
    def set_location(self, file):
        self.folder = file
        
    def get_df(self):
        return self.geolife_df
    

class FoliumCreator:
    def __init__(self, zoom=10, size='600px', center=BEIJING, scale_pos='topleft'):
        self.zoom = zoom
        self.size = size
        self.center = center # using American english spelling as that's what leaflet uses, not bc I like it :(
        self.layout = Layout(height=self.size)
        self.scale_pos = scale_pos
        
    def set_layout(self, layout):
        self.layout = layout
    
    def set_zoom(self, level):
        self.zoom = level
        
    def get_map(self):
        m = Map(center=self.center, zoom=self.zoom, layout=self.layout, basemap=basemaps.Stamen.Terrain)
        self.set_control(m, ScaleControl(position=self.scale_pos))
        return m
    
    def set_control(self, instance, control):
        instance.add_control(control)
        
    def set_layer(self, instance, layer):
        instance.add_layer(layer)

class UserObject:
    def __init__(self, scanner, panel_plots, folium_map, mapcontrol_object, dashboard_mode=False):
        self.id = -1
        self.trips = {}
        self.id_count = 182
        self.dashboard_mode = dashboard_mode
        self.scan_object = scanner
        self.panel_plots = panel_plots
        self.folium_map = folium_map
        self.mapper = mapcontrol_object
             
        self.bt_submit_user = widgets.Button(description='Select')
        self.bt_submit_user.on_click(self.submit_user)

        self.bt_random_user = widgets.Button(description='Randomize')
        self.bt_random_user.on_click(self.randomize_user)

        self.first_user_id = self.get_random_id()

        self.user_options = widgets.Dropdown(
            value=self.first_user_id,
            options=self.get_all_ids(),
            disabled=False,
            layout={'width': '30%'})
        self.user_options.observe(self.user_change, 'value')

        self.dp_userinfo = widgets.Output(layout={'border': 'solid black 1px',
                                                  'max_width': '360px'})
        self.clusterbox = widgets.Output(layout={'border': 'solid black 1px',
                                                 'max_width': '400px',
                                                 'min_width': '400px'})
        self.get_user_info(self.first_user_id)
        self.scan_object.reset(self.gdf, self, self.panel_plots)
        self.scan_object.show_plot()
        
        if self.dashboard_mode:
            self.bt_submit_user.disabled = True            
        
    def display(self):
        display(HBox([HTML('<b>User:  '), self.user_options, self.bt_random_user, self.bt_submit_user]))
        display(HBox([self.dp_userinfo, widgets.Output(layout={'min_width': '75px'}), self.clusterbox]))

    def submit_user(self, b):
        b.disabled = True
        self.bt_random_user.disabled = True
        self.user_options.disabled = True
        self.set_id(self.user_options.value)
        print('User {} selected!'.format(self.get_id()))
    
    def randomize_user(self, b):
        self.user_options.value = self.get_random_id()

    def user_change(self, value):
        self.get_user_info(value.new)
        
        if self.dashboard_mode:
            try:
                self.set_id(value.new)
                self.scan_object.reset(self.gdf, self, self.panel_plots)
                self.clean_gdf()
                self.get_antline()
                self.scan_object.show_plot()
            except Exception as e:
                print('ERROR during processing. Exception:')
                print(e)
            
    def set_id(self, new_id):
        self.id = new_id
        
    def get_id(self):
        return self.id
    
    # currently unused
    def set_clusterbox(self, content):
        with self.clusterbox:
            print(content)
    
    def get_all_ids(self):
        return [i for i in range(0, self.id_count)]
    
    def get_random_id(self):
        return random.randint(0, self.id_count - 1)
    
    def get_user_info(self, test_id):
        self.bt_random_user.description = 'Loading...'
        self.bt_random_user.disabled = True
        
        info = CON.execute('select count(*), count(distinct(label)), min(time), max(time) from {} where user_id = {}'
                           .format(TB_BEIJING, test_id)).first()
        
        gdf_simplified = gpd.GeoDataFrame.from_postgis((('SELECT '
            'st_x((st_dumppoints(st_removerepeatedpoints(st_collect(geom), 0.001))).geom) as lat, '
            'st_y((st_dumppoints(st_removerepeatedpoints(st_collect(geom), 0.001))).geom) as lon, '
            '(st_dumppoints(st_removerepeatedpoints(st_collect(geom), 0.001))).geom as geom '
            'FROM {} WHERE user_id = {}'.format(TB_BEIJING, test_id))), CON)
        self.gdf = gdf_simplified
        
        
        self.dp_userinfo.clear_output()
        with self.dp_userinfo:
            print('User {} Summary\n--------------------------------------------'.format(test_id))
            print(('Data points:             {}\nUnique Travel Methods:   {}\nData Begin:              {}'
                  '\nData End:                {}'.format(info[0], info[1], info[2], info[3])))

        
        self.bt_random_user.description = 'Randomize'
        self.bt_random_user.disabled = False
        
    def clean_gdf(self):
        self.gdf['label'] = self.scan_object.get_labels()
        self.gdf['time'] = 'datetime'

        total_count = self.gdf.count()[0]

        noise_index = self.gdf[self.gdf['label'] == -1].index  

        self.cleaned_gdf = self.gdf
        self.cleaned_gdf.drop(noise_index, inplace=True)

        clean_count = self.cleaned_gdf.count()[0]
        noise_count = total_count - clean_count
        
        self.clusterbox.clear_output()
        with self.clusterbox:
            print(('Out of {} points, {} were detected as noise ({}%). With noise removed, {} points remain '
                   'with {} cluster(s)\n\nNote that cluster colours between plots are not constant')
                 .format(total_count, noise_count, round(noise_count / total_count * 100, 2), clean_count, 
                         len(self.gdf['label'].unique())))
            
    def get_antline(self):
        self.mapper.set_zoom(9)
        
        self.folium_map.clear_layers()
        self.folium_map.add_layer(basemaps.Stamen.Terrain)

        by_cluster = self.cleaned_gdf.groupby('label')
        by_cluster = dict(list(by_cluster))

        for cluster in by_cluster.keys():
            current_gdf = by_cluster[cluster]
            current_geom_list = [(x, y) for x, y in zip(current_gdf['lon'], current_gdf['lat'])]

            ant_path = AntPath(locations=current_geom_list,
                              dash_array=[1, 10],
                              delay=3000,
                              color=COLOURS[cluster],
                              pulse_color='black')

            self.mapper.set_layer(self.folium_map, ant_path)
        
    
class DBSCANner:

    def __init__(self):
        self.backup_knee = 1
        
    def reset(self, gdf, user, sidebar, samples=20):
        self.samples = samples
        self.gdf = gdf
        self.geom_list = [(x, y) for x, y in zip(self.gdf['lat'], self.gdf['lon'])]
        self.panel = sidebar
        self.user = user
        
        self.panel.clear_output()
        with self.panel:
            self.fig, self.axes = plt.subplots(2, 1, figsize=(4, 8))
        
        self.distances = self.get_distances()
        self.knee = self.get_knee()
        self.labels = self.get_DBSCAN()
        self.gdf['cluster'] = self.labels
        
        self.plot_knee()
        self.plot_clusters()
        
        
    def get_distances(self):
        # get nearest neigbours
        # Americanized spelling as this is the name of the function, not because I like it
        try:
            if (len(self.geom_list) < self.samples + 1):
                self.samples = int(round(len(self.geom_list) / 2, 0))
            nn = NearestNeighbors(n_neighbors=self.samples + 1)
            neighbors = nn.fit(self.geom_list) 
            distances, indices = neighbors.kneighbors(self.geom_list)    
            distances = np.sort(distances[:, self.samples], axis=0)

            return distances
        except Exception as e:
            print('\nERROR during processing. Exception:')
            print(e)
        
    def get_knee(self):
        try:
            # calculate knee point
            i = np.arange(len(self.distances))
            knee = KneeLocator(i, self.distances, S=1, curve='convex', direction='increasing', interp_method='polynomial')

            # make sure knee returns as number
            if (isinstance(knee.knee, int)):
                knee_val = self.backup_knee
            else: 
                knee_val = knee.knee
            # save knee point
            return knee_val
        
        # lots of errors/exceptions occur on 1 or 2 users with very low n. handling a single one and returning
        #    is not feasible. return knee point of 1 instead, better than crashing
        except:
            return 1
    
    def get_DBSCAN(self):
        try:
            # calculate DBSCAN clusters (classify)
            dbscan = DBSCAN(eps=self.distances[self.knee], min_samples=self.samples).fit(self.geom_list)

            # get labels of clusters, return them as list
            labels = dbscan.labels_
            return labels
        except Exception as e:
            print('\nERROR during processing. Exception:')
            print(e)
    
    def plot_knee(self):
        with self.panel:
            try:
                self.axes[0].set_title('Knee Point of Nearest Neighbours for User {} Mobility'.format(self.user.get_id()))
                self.axes[0].plot(self.distances)
                self.axes[0].axvline(x=self.knee, color='gray', linestyle='--')
                #self.axes[0].text((self.knee, (0.5 * self.axes[0].get_ylim()[0])), 
                   #          'knee={}'.format(self.knee))
                self.axes[0].set_xlabel('Points')
                self.axes[0].set_ylabel('Distances')
            except Exception as e:
                print('ERROR during processing. Exception:')
                print(e)
            

    def plot_clusters(self):
        with self.panel:
            
            '''scatter_plot = sns.scatterplot(x=self.gdf['lat'], y=self.gdf['lon'], marker="o", ax=self.axes[1],
                    hue=["{cluster}".format(cluster = 'Cluster ' + 
                    str(cur_label) if cur_label != -1 else 'Noise') for cur_label in self.gdf['labels']])'''
            scatter_plot = sns.scatterplot(x=self.gdf['lat'], y=self.gdf['lon'], marker="o", ax=self.axes[1],
                    palette='Set2', hue=self.gdf['cluster'], legend='full')
            scatter_plot.set(xlabel='Latitude', ylabel='Longitude', 
                             title='Clusters in User {} Mobility using DBSCAN\n(where -1 is noise)'
                             .format(self.user.get_id()))
            
    def get_labels(self):
        return self.labels
    
    def show_plot(self):
        with self.panel:
            plt.tight_layout()
            plt.show()
        
if not DATA_EXISTS:
    bt_submit = widgets.Button(description='Select')

    file_options = widgets.Dropdown(
        options=[(os.path.join(os.getcwd(), 'geolife_trajectories\Data')), 
                 (os.path.join(os.path.join(os.environ['USERPROFILE']), 'Desktop\Data')), 
                 (os.path.join(os.path.join(os.environ['USERPROFILE']), 'Downloads\Data'))],
        disabled=False,
        layout={'width': '70%'})

    display(VBox([HTML('<h4>GeoLife File Location</h4>'), HBox([file_options, bt_submit])]))

    def bt_submitevent(b):
        if os.path.exists(file_options.value):
            b.disabled = True
            handler.set_location(file_options.value)
            handler.read_all_users()
        else:
            print('Data file not found. Please try again...')
    
    handler = GeoLifeHandler()
    bt_submit.on_click(bt_submitevent)
    
    try:
        geolife_raw = handler.get_df()
        display(geolife_raw.count())
        display(geolife_raw.head(10))
        geolife_raw = geolife_raw.head(10)
    except:
        print('Data either not yet extracted, or already available on server')
        
    uploader = PostGISUploader()
    uploader.upload(geolife_raw)
    
mapcontroller = FoliumCreator()
map_main = mapcontroller.get_map()
panel_plots = widgets.Output()

panel_top = HBox([map_main, panel_plots])
display(panel_top)

scanner = DBSCANner()
user = UserObject(scanner, panel_plots, map_main, mapcontroller, dashboard_mode=True)
#scanner.reset(gdf_simple, user, panel_plots)

user.display()
# TO DO
# -get first user working
# -fix issues
# -get time column working
# -sync colours of cluster plot and map<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup