3D Space-Time Cube.

Technologies Used: Python 2.X (OGR, Matplotlib)

This QGIS Plugin is able to transform a valid GPX file (with a time field) to a space time cube (QGIS Version 2.X). Uses PyQGIS and Matplotlib.

Description

The basic structure of a space time cube includes two of the cube’s axes to represent space (typically, longitude and latitude) and the third axe to visualize time (Kraak, 2003). This was completed in four main steps described in detail below.

Extracting latitude, longitude and time had to be from the GPX file using OGR

The method included first accessing the file’s data source, using OGR.GetDriverByName(“GPX”), and applying it to a user-set GPX file (using OS). Next, a layer file is created from the GPX data source. The default argument in the data_source.GetLayer() is 4, which reflects the track points included in the GPX file. This seems to function with a variety of GPX files from various sources, but can be adjusted if the schema in the GPX is different.

Processing the raw data using Python’s built-in time and date libraries.

Once a layer object has been successfully created, the following arrays are instantiated:

Array Variables	Content
x	Longitude
y	Latitude
z	Time in seconds since epoch
elapsedTime	Cumulative time in seconds since epoch (SSE)
DTArray	Time of each point as a string (example: “18:12:38”)
DateTimeArray	Time & date of each point as a python DateTime object (example: datetime.datetime(2017, 7, 14, 18, 12, 56) )

Storing the data in arrays

After the creation of the layer object, the code loops over each feature (representing a single track point in the GPX) in the layer, and extracts the data for the arrays seen in table 1. Whilst X,Y are extracted ‘as is’ as longitude/latitude respectively, time is handled through several steps. The data stored in the GPX is in ISO-8601 format(“2017/07/14 18:12:56+00”, and this is first converted to a python datetime object using the Strptime function with the following format ("%Y/%m/%d %H:%M:%S+00"). Next, the datetime object is converted to UNIX time (seconds since the epoch, or 00:00:00 Coordinated Universal Time (UTC), Thursday, 1 January 1970). This temporary variable (SSE) is added to an array to store the data – elapsedTime. Finally, every second entry in elapsedTime is converted to total minutes since start time. This is done by first adding 0 in the ‘z’ array, subtracting each entry by the previous entry, and dividing by 60.

The title for the graph is also extracted from the GPX dynamically, with the following format:

 First Time Stamp + “to” + Last Time Stamp 

Visualizing a 3D plot using Matplotlib.

The ultimate step for visualization is plotting the data. After having prepared the figure and the axis via Matplotlib and setting the projection to 3D, the ax.plot method only received the final arrays as argument: x, y and z. Finally, the x, y, and z labels were fixed to longitude, latitude, and time in minutes respectively. The figure was shown to the user using plt.show().

Code

This code is a ready-to-run Python 2 script, with the only required installation of OGR, available via PIP install. See the GitHub version for the QGIS 2.18 plugin version.

            import ogr, os
            import time, datetime
            import matplotlib as mpl
            from mpl_toolkits.mplot3d import Axes3D
            import matplotlib.pyplot as plt

            gpx_path = "C:/GPX/mytrack.gpx"

            #Open GPX
            in_path = os.path.join (gpx_path)
            in_driver= ogr.GetDriverByName("GPX")
            data_source = in_driver.Open(in_path,0)

            #Plot parameters
            mpl.rcParams ['legend.fontsize'] = 10
            fig = plt.figure()
            ax = fig.gca(projection='3d')

            #Create layer from data source
            if data_source is None:
                print "Could not open %s" % (in_path)
            else:
                print "Opened %s" % (in_path)
                layer = data_source.GetLayer(4) # 0 - waypoints, 1 - routes, 2 - tracks, 3 - route points, 4 - track points

            #Get time and spatial data from GPX dynamically and plot
            x=[] #Latitude
            y=[] #Latitude
            z=[] ##Time in minutes since start
            elapsedTime=[] #Time in seconds since epoch
            DTArray=[] ##Array of python strings
            DateTimeArray=[] #Python datetime object of each point
            schema=[]

            ldefn = layer.GetLayerDefn()
            for n in range(ldefn.GetFieldCount()):
                fdefn = ldefn.GetFieldDefn(n)
                schema.append(fdefn.name)

            try:
                timeIndex=schema.index('time')
            except ValueError:
                print "No time field found"

            try:
                for feat in layer:
                    pt = feat.geometry()
                    x.append(pt.GetX())
                    y.append(pt.GetY())
                    dateTime=feat.GetField(timeIndex)
                    try:
                        DT=datetime.datetime.strptime(dateTime, "%Y/%m/%d %H:%M:%S+00" )
                    except ValueError:
                        DT=datetime.datetime.strptime(dateTime, "%Y/%m/%d %H:%M:%S.%f+00" )
                    DateTimeArray.append(DT)
                    SSE= time.mktime(DT.timetuple()) #Seconds since epoch
                    elapsedTime.append(SSE)
                    DTArray.append(DT.strftime("%H:%M:%S"))
            except TypeError:
                print "Time field does contain none or invalid dates"


            #Create sub title for plot
            title=DateTimeArray[0].strftime("%Y/%m/%d")+" " +DTArray[0]+" to " + DateTimeArray[-1].strftime("%Y/%m/%d")+" "+DTArray[-1]

            #Extract elapsed time in minutes
            z.append(0) #First item in elapsedTime array should be 0
            for i in range(1,len(elapsedTime)):
                z.append((elapsedTime[i]-elapsedTime[0])/60)

            #Plot X,Y,Z
            ax.plot(x,y,z)

            #Labels & Plot
            ax.set_xlabel('Longitude')
            ax.set_ylabel('Latitude')
            ax.set_zlabel('Time (Minutes)')
            ax.legend()
            plt.axis("equal")
            fig.suptitle('Space Time Cube', fontsize=12, fontweight='bold')
            plt.title(title,loc='center')
            addImageOverlay(min(x),min(y),min(z))
            plt.show()

            #Optionally save the image
            #plt.savefig("C:/GPX/SpaceTimePlot.jpg", dpi=100, format="jpg")
            print "Done"
          

Installation

This plugin is officially published on the QGIS Plugin Repository, and thus can be installed within the software under Plugins --> Manage and Install Plugins.

Back to my projects →