[GIS] Using PostGIS for a large trajectory dataset

postgispostgis-2.0spatio-temporal-datatrackingtrajectory

I have a spatial trajectory dataset which basically contains about 10E6 (lat,lon,time)-tuples, grouped together by an identifier so they form trajectories.

Currently, SQLite is used to store the raw data.

Lat/Lon are in WGS84 but span only small areas (i.e. European country size).

I would like to store it in a geodatabase for easy descriptive & exploratory analysis, eg. calculating the overall and average distance for trajectories, and basic data preprocessing, eg. dismissing all tuples that are outside a polygon, for example.

Thus, I thought PostGIS would be optimal since I have worked with it before.

Some questions arise though: Should I just store the tuples as POINTS with GEOMETRY data type, or should I use MULTIPOINTS for whole trajectories?

And what about performance? Performance is not critical since the data is used for a research project and not for a production application, but I don't want to wait 5 hours to calculate distances between a million point-pairs.

What about indexing? What type of indices should I use and on what type of data?

I know this question is kind of general but I would like to see if people have experience with storing trajectory data and maybe they can give me some advice.

Best Answer

About different geometry types: From your description it looks like you should absolutely store your trajectories as linestrings. If you store them as points or multipoints you will have to build linestrings in runtime if you don't only want to do the calculations on the points defining the trajectories but also what is between the points. an Example (in meter-based projection)

select ST_Distance('MULTIPOINT(10 10, 10 20)'::geometry,'POINT(0 15)'::geometry)

Will return: 11.1803398874989 which is the distance from the closest point in the multipoint to the point. while:

select ST_Distance('LINESTRING(10 10, 10 20)'::geometry,'POINT(0 15)'::geometry)

will return 10 which is the distance from the point to the closest point to the closest point on the line.

The same things about intersections test. If you only uses the points, a trajectory passing over a polygon without points in the polygon, will not be found if you only calculate on the vertex.points and don't use the points to define linestrings.

About projections you should absolutely transform (with ST_Transform) your data to a local projection if possible. It is faster and you will get a lot more functions to choose from.

About R you should be aware that postgreSQL have have a procedural language handling R directly in the database: PL/R

About indexing: The spatial index used in postgis is GIST which builds an index of the bounding boxes.

Related Solutions

[GIS] From trajectory of GPS points to multiline segment

The best solution up till now implements everything in SQL with postgres. It uses moving averages (windowing function) to smoothen out the trajectory plus some simplification of the computed path. Here's the SQL code:

  select 
  --    st_simplify(
            st_curvetoline(st_linetocurve(st_makeline(the_geom)),1000)
  --    ,0.00001)
  from (
     select orderid, st_makepoint(avg(x) over w, avg(y) over w) as the_geom
     from (
        select orderid, x(the_geom) as x, y(the_geom) as y
        from tracepoints
     ) as moving_averages
  window w as (order by orderid rows between 1 preceding and 1 following)
  ) as make_the_line

uncomment some lines to remove simplification. One may also decide to drop the curvetoline and linetocurve functions

Detecting and Fixing Outliers in GPS Trajectory Data

Algorithm I use.

Calculate Euclidean minimum spanning tree of points:

Find 2 points most far apart from each other on this network

Find shortest route between them:

As one can see it might cut corner on a sharp turn.

I have ArcGIS python implementation of above algorithm, it uses networkx module. Let me know if this is of interest and I'll update my answer with the script

UPDATE:

# Connects points to make polyline. Makes 1 line at a time
# Tool assumes that 1st layer in Table of Conternt is TARGET polyline feature class,
# second layer in TOC is SOURCE point fc.
# If no selection found in SOURCE layer, works on entire dataset

import arcpy, traceback, os, sys
import itertools as itt
from math import sqrt
sys.path.append(r'C:\Users\felix_pertziger\AppData\Roaming\Python\Python27\site-packages')
import networkx as nx
from networkx import dijkstra_path_length

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    def CheckLayerLine(infc):
        d=arcpy.Describe(infc)
        theType=d.shapeType
        if theType!="Polyline":
            arcpy.AddWarning("\nTool designed to work with polylines as TARGET!")
            raise NameError, "Wrong input\n"
        return d
    def CheckLayerPoint(infc):
        d=arcpy.Describe(infc)
        theType=d.shapeType
        if theType!="Point":
            arcpy.AddWarning("\nTool designed to work with points as SOURCE!")
            raise NameError, "Wrong input\n"
        return d
    mxd = arcpy.mapping.MapDocument("CURRENT")
    layers = arcpy.mapping.ListLayers(mxd)
    if len(layers)<=1:
        arcpy.AddWarning("\nNot enough layers in the view!")
        raise NameError, "Wrong input\n"
    destLR, sourceLR=layers[0],layers[1]
    a = CheckLayerPoint(sourceLR);d = CheckLayerLine(destLR)

#  copy all points to manageable list
    g=arcpy.Geometry()
    geometryList=arcpy.CopyFeatures_management(sourceLR,g)
    nPoints=len(geometryList)
    arcpy.AddMessage('Computing minimum spanning tree')
    list2connect=[p.firstPoint for p in geometryList]
#  create network    
    p=list(itt.combinations(range(nPoints), 2))
    arcpy.SetProgressor("step", "", 0, len(p),1)
    G=nx.Graph()
    for f,t in p:
        p1=list2connect[f]
        p2=list2connect[t]
        dX=p2.X-p1.X;dY=p2.Y-p1.Y
        lenV=sqrt(dX*dX+dY*dY)
        G.add_edge(f,t,weight=lenV)
        arcpy.SetProgressorPosition()
    arcpy.AddMessage(len(G.edges()))
    mst=nx.minimum_spanning_tree(G)
    del G

#  find remotest pair
    arcpy.AddMessage(len(mst.edges()))
    length0=nx.all_pairs_dijkstra_path_length(mst)
    lMax=0
    for f,t in p:
        lCur=length0[f][t]
        if lCur>lMax:
            lMax=lCur
            best=(f,t)
    gL=nx.dijkstra_path(mst,best[0],best[1])
    del mst
    nPoints=len(gL)
    ordArray=arcpy.Array()
    for i in gL: ordArray.add(list2connect[i])

#  append line to TARGET
    curT = arcpy.da.InsertCursor(destLR,"SHAPE@")
    curT.insertRow((arcpy.Polyline(ordArray),))
    arcpy.RefreshActiveView()
    del curT

except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

Best Answer

Related Solutions

[GIS] From trajectory of GPS points to multiline segment

Detecting and Fixing Outliers in GPS Trajectory Data

Related Question